Fast simulation of propagating the photons created from SimEnergyDeposits. More...

Inheritance diagram for phot::PhotonLibraryPropagation:

Public Member Functions
	PhotonLibraryPropagation (fhicl::ParameterSet const &)

	PhotonLibraryPropagation (PhotonLibraryPropagation const &)=delete

	PhotonLibraryPropagation (PhotonLibraryPropagation &&)=delete

PhotonLibraryPropagation &	operator= (PhotonLibraryPropagation const &)=delete

PhotonLibraryPropagation &	operator= (PhotonLibraryPropagation &&)=delete

Private Member Functions
void	produce (art::Event &) override

Private Attributes
double	fRiseTimeFast

double	fRiseTimeSlow

bool	fDoSlowComponent

vector< art::InputTag >	fEDepTags

larg4::ISCalcSeparate	fISAlg

CLHEP::HepRandomEngine &	fPhotonEngine

CLHEP::HepRandomEngine &	fScintTimeEngine

Detailed Description

Fast simulation of propagating the photons created from SimEnergyDeposits.

This module does a fast simulation of propagating the photons created from SimEnergyDeposits, which is the Geant4 output after each step, to each of the optical detectors. This simulation is done using the PhotonLibrary, which stores the visibilities of each optical channel with respect to each optical voxel in the TPC volume, to avoid propagating single photons using Geant4. At the end of this module a collection of the propagated photons either as sim::SimPhotonsLite or as sim::SimPhotons is placed into the art event.

Keep in mind that at this stage the LArG4 main module is not capable of running the full optical simulation, because the necessary code has not yet been written.

In the future when the PhotonLibrary has the propagation time included, it could be possible to enhance sim::SimPhotons and sim::SimPhotonsLite to contain the propagation time. At this point the time recorded for the photon is the creation time of the photon.

The steps this module takes are:

to take sim::SimEnergyDeposits produced by larg4Main,
use Ionisation and Scintillation to calculate the amount of scintillated photons,
use the PhotonLibrary (visibilities) to determine the amount of visible photons at each optical channel,
- visible photons = the amount of scintillated photons calculated times the visibility at the middle of the Geant4 step for a given optical channel.
and if sim::SimPhotonsLite produced:
- since a sim::SimPhotonsLite only keeps a set of times with the number of photons produced at each time for a given OpChannel number:
  - for each time (as an integer in [ns]) photons are produced along the Geant4 step (taking into account the rise time and decay time of the fast and slow components of the scintillation process),
  - count the amount of photons emitted at that time.

the total amount of visible photons produced during the current Geant4 step equals the sum of counts for each time.

the total amount of visible photons produced during the current Geant4 step is determined by throwing a random number from a Poisson distribution with a mean of the amount of visible photons calculated above.

and if sim::SimPhotons produced:

since a sim::SimPhotons keeps a collection of photons (sim::OnePhoton) for a given OpChannel number:
- each single photon produced by this algorithm is just a copy containing the same information about:
  - energy (set to a constant value = 9.7e-6, which is equivalent to a wavelength of 128 nm, it should actually be 126.6 nm!!),
  - initial position,
  - time (as an integer in [ns]) the photon is produced along the Geant4 Step (taking into account the rise time and decay time of the fast and slow components of the scintillation process),
- the total amount of photon copies produced during the current Geant4 step is determined by throwing a random number from a Poisson distribution with a mean of the amount of visible photons calculated above.

This module should only be run for the fast optical simulation even though it can create sim::SimPhotonsLite and sim::SimPhotons as data products. If there is need to create sim::SimPhotons, there are some considerations you must be aware of. Since the amount of sim::SimPhotons produced even at low energies and in small geometries quickly exceeds the memory capacity of the job, right now it is actually impossible to produce sim::SimPhotons for any realistic geometry. A possible way around the problem is to implement a scaling of the produced sim::SimPhotons, to only produce a fraction of them.

Definition at line 145 of file PhotonLibraryPropagation_module.cc.

Constructor & Destructor Documentation

phot::PhotonLibraryPropagation::PhotonLibraryPropagation ( fhicl::ParameterSet const & p )

explicit

Definition at line 165 of file PhotonLibraryPropagation_module.cc.

     : art::EDProducer{p}
     , fRiseTimeFast{p.get<double>("RiseTimeFast", 0.0)}
     , fRiseTimeSlow{p.get<double>("RiseTimeSlow", 0.0)}
     , fDoSlowComponent{p.get<bool>("DoSlowComponent")}
     , fEDepTags{p.get<vector<art::InputTag>>("EDepModuleLabels")}
     , fPhotonEngine(art::ServiceHandle<rndm::NuRandomService> {}

phot::PhotonLibraryPropagation::PhotonLibraryPropagation ( PhotonLibraryPropagation const & )

delete

phot::PhotonLibraryPropagation::PhotonLibraryPropagation ( PhotonLibraryPropagation && )

delete

Member Function Documentation

PhotonLibraryPropagation& phot::PhotonLibraryPropagation::operator= ( PhotonLibraryPropagation const & )

delete

PhotonLibraryPropagation& phot::PhotonLibraryPropagation::operator= ( PhotonLibraryPropagation && )

delete

void phot::PhotonLibraryPropagation::produce ( art::Event & e )

overrideprivate

Definition at line 185 of file PhotonLibraryPropagation_module.cc.

   {
     art::ServiceHandle<PhotonVisibilityService const> pvs;
     art::ServiceHandle<sim::LArG4Parameters const> lgp;
     auto const* larp = lar::providerFrom<detinfo::LArPropertiesService>();
     auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e);
     CLHEP::RandPoissonQ randpoisphot{fPhotonEngine};
     CLHEP::RandFlat randflatscinttime{fScintTimeEngine};
     auto const nOpChannels = pvs->NOpChannels();
     unique_ptr<vector<sim::SimPhotons>> photCol{new vector<sim::SimPhotons>{}};
     auto& photonCollection{*photCol};
     photonCollection.resize(nOpChannels);
     unique_ptr<vector<sim::SimPhotonsLite>> photLiteCol{new vector<sim::SimPhotonsLite>{}};
     auto& photonLiteCollection{*photLiteCol};
     photonLiteCollection.resize(nOpChannels);
     for (unsigned int i = 0; i < nOpChannels; ++i) {
       photonLiteCollection[i].OpChannel = i;
       photonCollection[i].SetChannel(i);
     }
     vector<vector<sim::SimEnergyDeposit> const*> edep_vecs;
     for (auto label : fEDepTags) {
       auto const& edep_handle = e.getValidHandle<vector<sim::SimEnergyDeposit>>(label);
       edep_vecs.push_back(edep_handle);
     }
     for (auto const& edeps : edep_vecs) { //loop over modules
       for (auto const& edep : *edeps) {   //loop over energy deposits: one per step
         //int count_onePhot =0; // unused
         auto const& p = edep.MidPoint();
         auto const& Visibilities = pvs->GetAllVisibilities(p);
         if (!Visibilities) {
           throw cet::exception("PhotonLibraryPropagation")
             << "There is no entry in the PhotonLibrary for this position in space. "
                "Position: "
             << edep.MidPoint();
         }
         auto const isCalcData = fISAlg.CalcIonAndScint(detProp, edep);
         //total amount of scintillation photons
         double nphot = static_cast<int>(isCalcData.numPhotons);
         //amount of scintillated photons created via the fast scintillation process
         double nphot_fast = static_cast<int>(GetScintYield(edep, *larp) * nphot);
         //amount of scintillated photons created via the slow scintillation process
         double nphot_slow = nphot - nphot_fast;
         for (unsigned int channel = 0; channel < nOpChannels; ++channel) {
           auto visibleFraction = Visibilities[channel];
           if (visibleFraction == 0.0) {
             // Voxel is not visible at this optical channel, skip doing anything for this channel.
             continue;
           }
           if (lgp->UseLitePhotons()) {
             if (nphot_fast > 0) {
               //throwing a random number from a poisson distribution with a mean of the amount of photons visible at this channel
               auto n = static_cast<int>(randpoisphot.fire(nphot_fast * visibleFraction));
               for (long i = 0; i < n; ++i) {
                 //calculates the time at which the photon was produced
                 auto time = static_cast<int>(edep.T0() + GetScintTime(fRiseTimeFast,
                                                                       larp->ScintFastTimeConst(),
                                                                       randflatscinttime));
                 ++photonLiteCollection[channel].DetectedPhotons[time];
               }
             }
             if ((nphot_slow > 0) && fDoSlowComponent) {
               //throwing a random number from a poisson distribution with a mean of the amount of photons visible at this channel
               auto n = randpoisphot.fire(nphot_slow * visibleFraction);
               for (long i = 0; i < n; ++i) {
                 //calculates the time at which the photon was produced
                 auto time = static_cast<int>(edep.T0() + GetScintTime(fRiseTimeSlow,
                                                                       larp->ScintSlowTimeConst(),
                                                                       randflatscinttime));
                 ++photonLiteCollection[channel].DetectedPhotons[time];
               }
             }
           }
           else {
             sim::OnePhoton photon;
             photon.SetInSD = false;
             photon.InitialPosition = edep.End();
             photon.Energy = 9.7e-6;
             if (nphot_fast > 0) {
               //throwing a random number from a poisson distribution with a mean of the amount of photons visible at this channel
               auto n = randpoisphot.fire(nphot_fast * visibleFraction);
               if (n > 0) {
                 //calculates the time at which the photon was produced
                 photon.Time =
                   edep.T0() +
                   GetScintTime(fRiseTimeFast, larp->ScintFastTimeConst(), randflatscinttime);
                 // add n copies of sim::OnePhoton photon to the photon collection for a given OpChannel
                 photonCollection[channel].insert(photonCollection[channel].end(), n, photon);
               }
             }
             if ((nphot_slow > 0) && fDoSlowComponent) {
               //throwing a random number from a poisson distribution with a mean of the amount of photons visible at this channel
               auto n = randpoisphot.fire(nphot_slow * visibleFraction);
               if (n > 0) {
                 //calculates the time at which the photon was produced
                 photon.Time =
                   edep.T0() +
                   GetScintTime(fRiseTimeSlow, larp->ScintSlowTimeConst(), randflatscinttime);
                 // add n copies of sim::OnePhoton photon to the photon collection for a given OpChannel
                 photonCollection[channel].insert(photonCollection[channel].end(), n, photon);
               }
             }
           }
         }
       }
     }
     if (lgp->UseLitePhotons()) {
       // put the photon collection of LitePhotons into the art event
       e.put(move(photLiteCol));
     }
     else {
       //put the photon collection of SimPhotons into the art event
       e.put(move(photCol));
     }
   }