LightSource_module.cc

Go to the documentation of this file.

// C++ includes.
#include <cassert>
#include <cmath>
#include <fstream>
#include <memory>
#include <set>
#include <string>
#include <vector>

// ART includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "canvas/Utilities/Exception.h"
#include "cetlib/pow.h" // cet::square()
#include "cetlib_except/exception.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// art extensions
#include "nurandom/RandomUtils/NuRandomService.h"

// nusimdata includes
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// lar includes
#include "larcore/CoreUtils/ServiceUtil.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/GeometryCore.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
#include "larcoreobj/SummaryData/RunData.h"
#include "larsim/PhotonPropagation/PhotonVisibilityService.h"
#include "larsim/Simulation/PhotonVoxels.h"

#include "TGeoManager.h"
#include "TGeoMaterial.h"
#include "TGeoNavigator.h"
#include "TList.h"
#include "TLorentzVector.h"
#include "TTree.h"
#include "TVector3.h"

#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGaussQ.h"

namespace evgen {

  class LightSource : public art::EDProducer {
  public:
    explicit LightSource(fhicl::ParameterSet const& pset);

    void produce(art::Event& evt);
    void beginRun(art::Run& run);

  private:
    class MaterialPointFilter {
    public:
      MaterialPointFilter(geo::GeometryCore const& geom,
                          std::set<std::string> const& materialNames);

      ~MaterialPointFilter();

      // @{
      bool accept(geo::Point_t const& point);
      bool operator()(geo::Point_t const& point) { return accept(point); }
      // @}

    private:
      TGeoManager* fManager = nullptr;     
      TGeoNavigator* fNavigator = nullptr; 

      std::set<std::string> const& fMaterials;

      TGeoMaterial const* materialAt(geo::Point_t const& point);

      TGeoMaterial const* findMaterial(std::string const& name) const;

    }; // MaterialPointFilter

    simb::MCTruth Sample();

    void checkMaterials() const;

    bool readParametersFromInputFile();

    // for c2: fSeed is unused
    //int               fSeed;              //random number seed
    std::string fVersion; //version of the configuration

    // Flags to mark module modes
    static const int kUNIF = 0;
    static const int kGAUS = 1;
    static const int kFILE = 0;
    static const int kSCAN = 1;

    // File stream, filename and empty string for file processing
    std::ifstream fInputFile;
    std::string fFileName;
    char fDummyString[256];

    // A ttree to keep track of where particles have been shot - ends up in histos.root
    TTree* fPhotonsGenerated;
    TLorentzVector fShotPos;
    TLorentzVector fShotMom;
    Int_t fEvID;

    // Parameters loaded from config - both modes
    int fSourceMode; // Mode to run in - scan or file
    bool fFillTree;  // Do we want to create a TTree of shot particles?
    int fPosDist;    //
    int fTDist;      // Random distributions to use : 1= gauss, 0= uniform
    int fPDist;      //

    std::set<std::string> const fSelectedMaterials;

    //Scan mode specific parameters
    int fXSteps; //
    int fYSteps; //  Number of steps to take in each dimension
    int fZSteps; //

    sim::PhotonVoxelDef fThePhotonVoxelDef; // The photon voxel definition object for scan mode

    int fVoxelCount;   // Total number of voxels
    int fCurrentVoxel; // Counter to keep track of vox ID

    //  TPC Measurements
    geo::Vector_t fTPCCenter;
    std::vector<double> fRegionMin;
    std::vector<double> fRegionMax;
    bool fUseCustomRegion;

    // Parameters used to shoot in distributions
    geo::Point_t fCenter; 
    bool fPointSource;    // Point-like light source in fCenter
    double fT;            // central t position of source
    double fSigmaX;       // x width
    double fSigmaY;       // y width
    double fSigmaZ;       // z width
    double fSigmaT;       // t width
    double fP;            // central momentm of photon
    double fSigmaP;       // mom width;

    // Number of photons per event
    int fN; // number of photons per event

    double const fNMaxF;

    int fFirstVoxel;
    int fLastVoxel;

    CLHEP::HepRandomEngine& fEngine;
    geo::GeometryCore const& fGeom; 
  };
}

namespace {

  template <typename Coll>
  std::set<typename Coll::value_type> makeSet(Coll const& coll)
  {
    return {begin(coll), end(coll)};
  }

} // local namespace

namespace evgen {

  //----------------------------------------------------------------
  LightSource::LightSource(fhicl::ParameterSet const& pset)
    : art::EDProducer{pset}
    , fSourceMode{pset.get<int>("SourceMode")}
    , fFillTree{pset.get<bool>("FillTree")}
    , fPosDist{pset.get<int>("PosDist")}
    , fTDist{pset.get<int>("TDist")}
    , fPDist{pset.get<int>("PDist")}
    , fSelectedMaterials{makeSet(pset.get<std::vector<std::string>>("SelectMaterials", {}))}
    , fNMaxF{pset.get<double>("NMaxFactor", 100.0)}
    // create a default random engine; obtain the random seed from NuRandomService,
    // unless overridden in configuration with key "Seed"
    , fEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(createEngine(0),
                                                                                 pset,
                                                                                 "Seed"))
    , fGeom(*lar::providerFrom<geo::Geometry const>())
  {

    checkMaterials();

    // load optional parameters in function
    produces<sumdata::RunData, art::InRun>();
    produces<std::vector<simb::MCTruth>>();

    if (fSourceMode == kFILE) {
      fFileName = pset.get<std::string>("SteeringFile");
      fInputFile.open(fFileName.c_str());
      fInputFile.getline(fDummyString, 256);
    }
    else if (fSourceMode == kSCAN) {
      fT = pset.get<double>("T0");
      fSigmaT = pset.get<double>("SigmaT");
      fN = pset.get<int>("N");

      fFirstVoxel = pset.get<int>("FirstVoxel");
      fLastVoxel = pset.get<int>("LastVoxel");

      fP = pset.get<double>("P");
      fSigmaP = pset.get<double>("SigmaP");

      fUseCustomRegion = pset.get<bool>("UseCustomRegion");
      fPointSource = pset.get<bool>("PointSource", false);

      if (fUseCustomRegion) {
        fRegionMin = pset.get<std::vector<double>>("RegionMin");
        fRegionMax = pset.get<std::vector<double>>("RegionMax");
        fXSteps = pset.get<int>("XSteps");
        fYSteps = pset.get<int>("YSteps");
        fZSteps = pset.get<int>("ZSteps");
      }

      // get TPC dimensions removed. -TA

      fCurrentVoxel = 0;

      // define voxelization based on parameters read from config.
      // There are two modes - either read the dimensions of the TPC from
      // the geometry, or use values specified by the user.
      if (!fUseCustomRegion) {
        art::ServiceHandle<phot::PhotonVisibilityService const> vis;
        fThePhotonVoxelDef = vis->GetVoxelDef();
      }
      else {
        fThePhotonVoxelDef = sim::PhotonVoxelDef(fRegionMin[0],
                                                 fRegionMax[0],
                                                 fXSteps,
                                                 fRegionMin[1],
                                                 fRegionMax[1],
                                                 fYSteps,
                                                 fRegionMin[2],
                                                 fRegionMax[2],
                                                 fZSteps);
      }

      // Set distribution widths to voxel size

      fSigmaX = fThePhotonVoxelDef.GetVoxelSize().X() / 2.0;
      fSigmaY = fThePhotonVoxelDef.GetVoxelSize().Y() / 2.0;
      fSigmaZ = fThePhotonVoxelDef.GetVoxelSize().Z() / 2.0;

      // Get number of voxels we will step through

      fVoxelCount = fThePhotonVoxelDef.GetNVoxels();

      if (fLastVoxel < 0) fLastVoxel = fVoxelCount;

      mf::LogVerbatim("LightSource") << "Light Source : Determining voxel params : " << fVoxelCount
                                     << " " << fSigmaX << " " << fSigmaY << " " << fSigmaZ;
    }
    else {
      throw cet::exception("LightSource")
        << "EVGEN Light Source : Unrecognised light source mode\n";
    }

    if (fFillTree) {
      art::ServiceHandle<art::TFileService const> tfs;
      fPhotonsGenerated = tfs->make<TTree>("PhotonsGenerated", "PhotonsGenerated");
      fPhotonsGenerated->Branch("X", &(fShotPos[0]), "X/D");
      fPhotonsGenerated->Branch("Y", &(fShotPos[1]), "Y/D");
      fPhotonsGenerated->Branch("Z", &(fShotPos[2]), "Z/D");
      fPhotonsGenerated->Branch("T", &(fShotPos[3]), "T/D");
      fPhotonsGenerated->Branch("PX", &(fShotMom[0]), "PX/D");
      fPhotonsGenerated->Branch("PY", &(fShotMom[1]), "PY/D");
      fPhotonsGenerated->Branch("PZ", &(fShotMom[2]), "PZ/D");
      fPhotonsGenerated->Branch("PT", &(fShotMom[3]), "PT/D");
      fPhotonsGenerated->Branch("EventID", &fEvID, "EventID/I");
    }
  }

  //____________________________________________________________________________
  void LightSource::beginRun(art::Run& run)
  {
    run.put(std::make_unique<sumdata::RunData>(fGeom.DetectorName()), art::fullRun());

    fCurrentVoxel = fFirstVoxel;
  }

  //----------------------------------------------------------------
  void LightSource::produce(art::Event& evt)
  {
    if (fSourceMode == kFILE) {
      //  Each event, read coordinates of gun and number of photons to shoot from file

      // read in one line
      if (!readParametersFromInputFile()) {
        // Loop file if required
        mf::LogWarning("LightSource") << "EVGEN Light Source : Warning, reached end of file,"
                                      << " looping back to beginning";
        fInputFile.clear();
        fInputFile.seekg(0, std::ios::beg);
        fInputFile.getline(fDummyString, 256);

        if (!readParametersFromInputFile()) {
          throw cet::exception("LightSource")
            << "EVGEN Light Source : File error in " << fFileName << "\n";
        }
      }

      fThePhotonVoxelDef = sim::PhotonVoxelDef(fCenter.X() - fSigmaX,
                                               fCenter.X() + fSigmaX,
                                               1,
                                               fCenter.Y() - fSigmaY,
                                               fCenter.Y() + fSigmaY,
                                               1,
                                               fCenter.Z() - fSigmaZ,
                                               fCenter.Z() + fSigmaZ,
                                               1);

      fCurrentVoxel = 0;
    }
    else if (fSourceMode == kSCAN) {
      //  Step through detector using a number of steps provided in the config file
      //  firing a constant number of photons from each point
      fCenter = fThePhotonVoxelDef.GetPhotonVoxel(fCurrentVoxel).GetCenter();
    }
    else {
      //  Neither file or scan mode, probably a config file error
      throw cet::exception("LightSource") << "EVGEN : Light Source, unrecognised source mode\n";
    }

    auto truthcol = std::make_unique<std::vector<simb::MCTruth>>();

    truthcol->push_back(Sample());
    int const nPhotons = truthcol->back().NParticles();

    evt.put(std::move(truthcol));

    phot::PhotonVisibilityService* vis = nullptr;
    try {
      vis = art::ServiceHandle<phot::PhotonVisibilityService>().get();
    }
    catch (art::Exception const& e) {
      // if the service is not configured, then this is not a build job
      // (it is a simple generation job instead)
      if (e.categoryCode() != art::errors::ServiceNotFound) throw;
    }

    if (vis && vis->IsBuildJob()) {
      mf::LogVerbatim("LightSource") << "Light source : Stowing voxel params ";
      vis->StoreLightProd(fCurrentVoxel, nPhotons);
    }

    if (fCurrentVoxel != fLastVoxel) { ++fCurrentVoxel; }
    else {
      mf::LogVerbatim("LightSource")
        << "EVGEN Light Source fully scanned detector.  Starting over.";
      fCurrentVoxel = fFirstVoxel;
    }
  }

  simb::MCTruth LightSource::Sample()
  {
    mf::LogVerbatim("LightSource") << "Light source debug : Shooting at " << fCenter;

    CLHEP::RandFlat flat(fEngine, -1.0, 1.0);
    CLHEP::RandGaussQ gauss(fEngine);

    MaterialPointFilter filter(fGeom, fSelectedMaterials);

    simb::MCTruth mct;
    mct.SetOrigin(simb::kSingleParticle);

    unsigned long long int const nMax = static_cast<unsigned long long int>(double(fN) * fNMaxF);
    unsigned long long int fired = 0ULL;
    while (mct.NParticles() < fN) {
      if (fired >= nMax) break;

      // Choose momentum (supplied in eV, convert to GeV)
      double const p =
        1e-9 * ((fPDist == kGAUS) ? gauss.fire(fP, fSigmaP) : fP + fSigmaP * flat.fire());

      // Choose position
      ++fired;
      geo::Point_t x;
      if (fPointSource) { x = fCenter; }
      else {
        if (fPosDist == kGAUS) {
          x = {gauss.fire(fCenter.X(), fSigmaX),
               gauss.fire(fCenter.Y(), fSigmaY),
               gauss.fire(fCenter.Z(), fSigmaZ)};
        }
        else {
          x = {fCenter.X() + fSigmaX * flat.fire(),
               fCenter.Y() + fSigmaY * flat.fire(),
               fCenter.Z() + fSigmaZ * flat.fire()};
        }

        if (!filter.accept(x)) continue;
      }

      // Choose time
      double t;
      if (fTDist == kGAUS) { t = gauss.fire(fT, fSigmaT); }
      else {
        t = fT + fSigmaT * flat.fire();
      }

      //assume the position is relative to the center of the TPC
      //x += fTPCCenter;

      fShotPos = TLorentzVector(x.X(), x.Y(), x.Z(), t);

      // Choose angles
      double costh = flat.fire();
      double sinth = std::sqrt(1.0 - cet::square(costh));
      double phi = 2 * M_PI * flat.fire();

      // Generate momentum 4-vector

      fShotMom = TLorentzVector(p * sinth * cos(phi), p * sinth * sin(phi), p * costh, p);

      int trackid = -(mct.NParticles() +
                      1); // set track id to -i as these are all primary particles and have id <= 0
      std::string primary("primary");
      int PDG = 0; //optical photons have PDG 0

      simb::MCParticle part(trackid, PDG, primary);
      part.AddTrajectoryPoint(fShotPos, fShotMom);

      if (fFillTree) fPhotonsGenerated->Fill();

      mct.Add(part);
    }

    mf::LogInfo("LightSource") << "Generated " << mct.NParticles() << " photons after " << fired
                               << " tries.";
    if (mct.NParticles() < fN) {
      // this may mean `NMaxFactor` is too small, or the volume is wrong;
      // or it may be just expected
      mf::LogWarning("LightSource")
        << "Warning: " << mct.NParticles() << " photons generated after " << fired << " tries, but "
        << fN << " were requested.";
    }

    return mct;
  } // LightSource::Sample()

  // ---------------------------------------------------------------------------
  void LightSource::checkMaterials() const
  {

    TGeoManager const& manager = *(fGeom.ROOTGeoManager());

    { // start scope
      mf::LogDebug log("LightSource");
      auto const& matList = *(manager.GetListOfMaterials());
      log << matList.GetSize() << " elements/materials in the geometry:";
      for (auto const* obj : matList) {
        auto const mat = dynamic_cast<TGeoMaterial const*>(obj);
        log << "\n  '" << mat->GetName() << "' (Z=" << mat->GetZ() << " A=" << mat->GetA() << ")";
      } // for
    }   // end scope

    std::set<std::string> missingMaterials;
    for (auto const& matName : fSelectedMaterials) {
      if (!manager.GetMaterial(matName.c_str())) missingMaterials.insert(matName);
    }
    if (missingMaterials.empty()) return;

    art::Exception e(art::errors::Configuration);
    e << "Requested filtering on " << missingMaterials.size()
      << " materials which are not present in the geometry:";
    for (auto const& matName : missingMaterials)
      e << "\n  '" << matName << "'";
    throw e << "\n";

  } // LightSource::checkMaterials()

  // ---------------------------------------------------------------------------
  // --- LightSource::MaterialPointFilter
  // ---------------------------------------------------------------------------
  LightSource::MaterialPointFilter::MaterialPointFilter(geo::GeometryCore const& geom,
                                                        std::set<std::string> const& materialNames)
    : fManager(geom.ROOTGeoManager())
    , fNavigator(fManager->AddNavigator())
    , fMaterials(materialNames)
  {
    assert(fManager);
    assert(fNavigator);
  }

  // ---------------------------------------------------------------------------
  LightSource::MaterialPointFilter::~MaterialPointFilter()
  {
    fManager->RemoveNavigator(fNavigator); // this deletes the navigator
    fNavigator = nullptr;
  } // LightSource::MaterialPointFilter::~MaterialPointFilter()

  // ---------------------------------------------------------------------------
  TGeoMaterial const* LightSource::MaterialPointFilter::materialAt(geo::Point_t const& point)
  {
    TGeoNode const* node = fNavigator->FindNode(point.X(), point.Y(), point.Z());
    return node ? node->GetVolume()->GetMaterial() : nullptr;
  } // LightSource::MaterialPointFilter::materialAt()

  // ---------------------------------------------------------------------------
  bool LightSource::MaterialPointFilter::accept(geo::Point_t const& point)
  {
    if (fMaterials.empty()) return true;
    TGeoMaterial const* material = materialAt(point);
    MF_LOG_TRACE("LightSource") << "Material at " << point << ": "
                                << (material ? material->GetName() : "not found");
    return material ? (fMaterials.count(material->GetName()) > 0) : false;
  } // LightSource::MaterialPointFilter::accept()

  //----------------------------------------------------------------------------
  bool LightSource::readParametersFromInputFile()
  {
    double x, y, z;
    fInputFile >> x >> y >> z >> fT >> fSigmaX >> fSigmaY >> fSigmaZ >> fSigmaT >> fP >> fSigmaP >>
      fN;
    fCenter = {x, y, z};
    if (!fInputFile.good()) return false;

    std::string dummy;
    std::getline(fInputFile, dummy); // this can fail for what I care
    return true;
  } // LightSource::readParametersFromInputFile()

  // ---------------------------------------------------------------------------

} // namespace evgen

DEFINE_ART_MODULE(evgen::LightSource)