opdet::SimPhotonCounter Class Reference

#include <SimPhotonCounter.h>

Inheritance diagram for opdet::SimPhotonCounter:

Public Member Functions
	SimPhotonCounter ()
	SimPhotonCounter (size_t s, float t_p1, float t_p2, float t_l1, float t_l2, float min_w=0, float max_w=1e6, float e=1.0)
	SimPhotonCounter (float t_p1, float t_p2, float t_l1, float t_l2, float min_w, float max_w, const std::vector< float > &eV)
void	SetVectorSize (size_t s)
size_t	GetVectorSize () const
void	SetWavelengthRanges (float min_w, float max_w)
float	MinWavelength () const
float	MaxWavelength () const
void	SetTimeRanges (float t_p1, float t_p2, float t_l1, float t_l2)
float	MinPromptTime () const
float	MaxPromptTime () const
float	MinLateTime () const
float	MaxLateTime () const
void	SetQE (size_t i, float e)
float	QE (size_t i) const
void	SetQEVector (const std::vector< float > &eV)
std::vector< float > const &	QEVector () const
void	AddOnePhoton (size_t i_opdet, const sim::OnePhoton &photon)
void	AddSimPhotons (const sim::SimPhotons &photons)
void	ClearVectors ()
const std::vector< float > &	PromptPhotonVector () const
const std::vector< float > &	LatePhotonVector () const
float	PromptPhotonVector (size_t i) const
float	LatePhotonVector (size_t i) const
std::vector< float >	TotalPhotonVector () const
float	TotalPhotonVector (size_t i) const
float	PromptPhotonTotal () const
float	LatePhotonTotal () const
float	PhotonTotal () const
void	Print ()
	SimPhotonCounter (const fhicl::ParameterSet &)
void	analyze (art::Event const &)
void	beginJob ()
void	endJob ()

Private Member Functions
float	Wavelength (const sim::OnePhoton &ph)
void	storeVisibility (int channel, int nDirectPhotons, int nReflectedPhotons, double reflectedT0=0.0) const
bool	isVisible (double wavelength) const
	Returns if we label as "visibile" a photon with specified wavelength [nm]. More...

Private Attributes
std::vector< float >	_photonVector_prompt
std::vector< float >	_photonVector_late
float	_min_prompt_time
float	_max_prompt_time
float	_min_late_time
float	_max_late_time
std::vector< float >	_qeVector
float	_min_wavelength
float	_max_wavelength
TTree *	fThePhotonTreeAll
TTree *	fThePhotonTreeDetected
TTree *	fTheOpDetTree
TTree *	fTheEventTree
std::vector< std::string >	fInputModule
int	fVerbosity
bool	fMakeDetectedPhotonsTree
bool	fMakeAllPhotonsTree
bool	fMakeOpDetsTree
bool	fMakeOpDetEventsTree
TVector3	initialPhotonPosition
TVector3	finalPhotonPosition
Float_t	fWavelength
Float_t	fTime
Int_t	fCountOpDetAll
Int_t	fCountOpDetDetected
Int_t	fCountOpDetReflDetected
Float_t	fT0_vis
Int_t	fCountEventAll
Int_t	fCountEventDetected
Int_t	fEventID
Int_t	fOpChannel
bool	fMakeLightAnalysisTree
std::vector< std::vector < std::vector< double > > >	fSignals_vuv
std::vector< std::vector < std::vector< double > > >	fSignals_vis
TTree *	fLightAnalysisTree = nullptr
int	fRun
int	fTrackID
int	fpdg
int	fmotherTrackID
double	fEnergy
double	fdEdx
std::vector< double >	fPosition0
std::vector< std::vector < double > >	fstepPositions
std::vector< double >	fstepTimes
std::vector< std::vector < double > >	fSignalsvuv
std::vector< std::vector < double > >	fSignalsvis
std::string	fProcess
cheat::ParticleInventoryService const *	pi_serv = nullptr
phot::PhotonVisibilityService const *	fPVS = nullptr
bool const	fUseLitePhotons

Static Private Attributes
static constexpr double	kVisibleThreshold = 200.0
	Threshold used to resolve between visible and ultraviolet light. More...
static constexpr double	kVisibleWavelength = 450.0
	Value used when a typical visible light wavelength is needed. More...
static constexpr double	kVUVWavelength = 128.0
	Value used when a typical ultraviolet light wavelength is needed. More...

Detailed Description

Definition at line 22 of file SimPhotonCounter.h.

Constructor & Destructor Documentation

opdet::SimPhotonCounter::SimPhotonCounter ( )

inline

Definition at line 25 of file SimPhotonCounter.h.

{}

opdet::SimPhotonCounter::SimPhotonCounter	(	size_t	s,
		float	t_p1,
		float	t_p2,
		float	t_l1,
		float	t_l2,
		float	min_w = 0,
		float	max_w = 1e6,
		float	e = 1.0
	)

Definition at line 10 of file SimPhotonCounter.cxx.

{

  SetWavelengthRanges(min_w,max_w);
  SetTimeRanges(t_p1,t_p2,t_l1,t_l2);

  _photonVector_prompt=std::vector<float>(s);
  _photonVector_late=std::vector<float>(s);
  _qeVector = std::vector<float>(s,e);

}

opdet::SimPhotonCounter::SimPhotonCounter	(	float	t_p1,
		float	t_p2,
		float	t_l1,
		float	t_l2,
		float	min_w,
		float	max_w,
		const std::vector< float > &	eV
	)

Definition at line 26 of file SimPhotonCounter.cxx.

{
  SetWavelengthRanges(min_w,max_w);
  SetTimeRanges(t_p1,t_p2,t_l1,t_l2);

  _photonVector_prompt=std::vector<float>(eV.size());
  _photonVector_late=std::vector<float>(eV.size());
  _qeVector = eV;

}

opdet::SimPhotonCounter::SimPhotonCounter

(

const fhicl::ParameterSet &

pset

)

Definition at line 173 of file SimPhotonCounter_module.cc.

    : EDAnalyzer(pset)
    , fPVS(art::ServiceHandle<phot::PhotonVisibilityService const>().get())
    , fUseLitePhotons(art::ServiceHandle<sim::LArG4Parameters const>()->UseLitePhotons())
  {
    fVerbosity=                pset.get<int>("Verbosity");
    try
    {
       fInputModule = pset.get<std::vector<std::string>>("InputModule",{"largeant"});
    }
    catch(...)
    {
       fInputModule.push_back(pset.get<std::string>("InputModule","largeant"));
    }
    fMakeAllPhotonsTree=       pset.get<bool>("MakeAllPhotonsTree");
    fMakeDetectedPhotonsTree=  pset.get<bool>("MakeDetectedPhotonsTree");
    fMakeOpDetsTree=           pset.get<bool>("MakeOpDetsTree");
    fMakeOpDetEventsTree=      pset.get<bool>("MakeOpDetEventsTree");
    fMakeLightAnalysisTree=    pset.get<bool>("MakeLightAnalysisTree", false);
    //fQE=                       pset.get<double>("QuantumEfficiency");
    //fWavelengthCutLow=         pset.get<double>("WavelengthCutLow");
    //fWavelengthCutHigh=        pset.get<double>("WavelengthCutHigh");


    if (fPVS->IsBuildJob() && fPVS->StoreReflected() && fPVS->StoreReflT0() && fUseLitePhotons) {
      throw art::Exception(art::errors::Configuration)
        << "Building a library with reflected light time is not supported when using SimPhotonsLite.\n";
    }

  }

Member Function Documentation

void opdet::SimPhotonCounter::AddOnePhoton	(	size_t	i_opdet,
		const sim::OnePhoton &	photon
	)

Definition at line 66 of file SimPhotonCounter.cxx.

{
  if(i_opdet > GetVectorSize())
    throw std::runtime_error("ERROR in SimPhotonCounter: Opdet requested out of range!");

  if(Wavelength(photon) < _min_wavelength || Wavelength(photon) > _max_wavelength) return;

  if(photon.Time > _min_prompt_time && photon.Time <= _max_prompt_time)
    _photonVector_prompt[i_opdet] += _qeVector[i_opdet];
  else if(photon.Time > _min_late_time && photon.Time < _max_late_time)
    _photonVector_late[i_opdet] += _qeVector[i_opdet];

}

void opdet::SimPhotonCounter::AddSimPhotons

(

const sim::SimPhotons &

photons

)

Definition at line 80 of file SimPhotonCounter.cxx.

{
  for(size_t i_ph=0; i_ph < photons.size(); i_ph++)
    AddOnePhoton(photons.OpChannel(),photons[i_ph]);
}

void opdet::SimPhotonCounter::analyze

(

art::Event const &

evt

)

Definition at line 307 of file SimPhotonCounter_module.cc.

  {

    // Lookup event ID from event
    art::EventNumber_t event = evt.id().event();
    fEventID=Int_t(event);

    // Service for determining opdet responses
    art::ServiceHandle<opdet::OpDetResponseInterface const> odresponse;

    // get the geometry to be able to figure out signal types and chan -> plane mappings
    art::ServiceHandle<geo::Geometry const> geo;

    // GEANT4 info on the particles (only used if making light analysis tree)
    std::vector<simb::MCParticle> const* mcpartVec = nullptr;

    //-------------------------initializing light tree vectors------------------------
    std::vector<double> totalEnergy_track;
    fstepPositions.clear();
    fstepTimes.clear();
    if (fMakeLightAnalysisTree) {
      //mcpartVec = evt.getPointerByLabel<std::vector<simb::MCParticle>>("largeant");
      mcpartVec = evt.getHandle<std::vector<simb::MCParticle>>("largeant").product();

      size_t maxNtracks = 1000U; // mcpartVec->size(); --- { to be fixed soon! ]
      fSignals_vuv.clear();
      fSignals_vuv.resize(maxNtracks);
      fSignals_vis.clear();
      fSignals_vis.resize(maxNtracks);
      for(size_t itrack=0; itrack!=maxNtracks; itrack++) {
        fSignals_vuv[itrack].resize(geo->NOpChannels());
        fSignals_vis[itrack].resize(geo->NOpChannels());
      }
      totalEnergy_track.resize(maxNtracks, 0.);
      //-------------------------stimation of dedx per trackID----------------------
      //get the list of particles from this event
      const sim::ParticleList* plist = pi_serv? &(pi_serv->ParticleList()): nullptr;

      // loop over all sim::SimChannels in the event and make sure there are no
      // sim::IDEs with trackID values that are not in the sim::ParticleList
      std::vector<const sim::SimChannel*> sccol;
      //evt.getView(fG4ModuleLabel, sccol);
     for(auto const& mod : fInputModule){
      evt.getView(mod, sccol);
      double totalCharge=0.0;
      double totalEnergy=0.0;
      //loop over the sim channels collection
      for(size_t sc = 0; sc < sccol.size(); ++sc){
        double numIDEs=0.0;
        double scCharge=0.0;
        double scEnergy=0.0;
        const auto & tdcidemap = sccol[sc]->TDCIDEMap();
        //loop over all of the tdc IDE map objects
        for(auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++){
          const std::vector<sim::IDE> idevec = (*mapitr).second;
          numIDEs += idevec.size();
          //go over all of the IDEs in a given simchannel
          for(size_t iv = 0; iv < idevec.size(); ++iv){
            if (plist) {
              if(plist->find( idevec[iv].trackID ) == plist->end()
                  && idevec[iv].trackID != sim::NoParticleId)
              {
                mf::LogWarning("LArG4Ana") << idevec[iv].trackID << " is not in particle list";
              }
            }
            if(idevec[iv].trackID < 0) continue;
            totalCharge +=idevec[iv].numElectrons;
            scCharge += idevec[iv].numElectrons;
            totalEnergy +=idevec[iv].energy;
            scEnergy += idevec[iv].energy;

            totalEnergy_track[idevec[iv].trackID] += idevec[iv].energy/3.;
          }
        }
      }
     }
    }//End of if(fMakeLightAnalysisTree)


    if(!fUseLitePhotons)
    {

      //Reset counters
      fCountEventAll=0;
      fCountEventDetected=0;

      //Get *ALL* SimPhotonsCollection from Event
      //std::vector< art::Handle< std::vector< sim::SimPhotons > > > photon_handles;
      //evt.getManyByType(photon_handles);
      auto photon_handles = evt.getMany<std::vector<sim::SimPhotons>>();
      if (photon_handles.size() == 0)
        throw art::Exception(art::errors::ProductNotFound)<<"sim SimPhotons retrieved and you requested them.";

      for(auto const& mod : fInputModule){
      // sim::SimPhotonsCollection TheHitCollection = sim::SimListUtils::GetSimPhotonsCollection(evt,mod);
      //switching off to add reading in of labelled collections: Andrzej, 02/26/19

        for (auto const& ph_handle: photon_handles) {
          // Do some checking before we proceed
          if (!ph_handle.isValid()) continue;
          if (ph_handle.provenance()->moduleLabel() != mod) continue;   //not the most efficient way of doing this, but preserves the logic of the module. Andrzej

          bool Reflected = (ph_handle.provenance()->productInstanceName() == "Reflected");

          if((*ph_handle).size()>0)
          {
            if(fMakeLightAnalysisTree) {
            //resetting the signalt to save in the analysis tree per event
              const int maxNtracks = 1000;
              for(size_t itrack=0; itrack!=maxNtracks; itrack++) {
                for(size_t pmt_i=0; pmt_i!=geo->NOpChannels(); pmt_i++) {
                  fSignals_vuv[itrack][pmt_i].clear();
                  fSignals_vis[itrack][pmt_i].clear();
                }
              }
            }
          }


  //      if(fVerbosity > 0) std::cout<<"Found OpDet hit collection of size "<< TheHitCollection.size()<<std::endl;
          if(fVerbosity > 0) std::cout<<"Found OpDet hit collection of size "<< (*ph_handle).size()<<std::endl;

          if((*ph_handle).size()>0)
          {
  //           for(sim::SimPhotonsCollection::const_iterator itOpDet=TheHitCollection.begin(); itOpDet!=TheHitCollection.end(); itOpDet++)
            for(auto const& itOpDet: (*ph_handle) )
            {
              //Reset Counters
              fCountOpDetAll=0;
              fCountOpDetDetected=0;
              fCountOpDetReflDetected=0;
              //Reset t0 for visible light
              fT0_vis = 999.;

              //Get data from HitCollection entry
              fOpChannel=itOpDet.OpChannel();
              const sim::SimPhotons& TheHit=itOpDet;

              //std::cout<<"OpDet " << fOpChannel << " has size " << TheHit.size()<<std::endl;

              // Loop through OpDet phots.
              //   Note we make the screen output decision outside the loop
              //   in order to avoid evaluating large numbers of unnecessary
              //   if conditions.


              for(const sim::OnePhoton& Phot: TheHit)
              {
                // Calculate wavelength in nm
                fWavelength= odresponse->wavelength(Phot.Energy);

                //Get arrival time from phot
                fTime= Phot.Time;

                // special case for LibraryBuildJob: no working "Reflected" handle and all photons stored in single object - must sort using wavelength instead
                if(fPVS->IsBuildJob() && !Reflected) {
                  // all photons contained in object with Reflected = false flag
                  // Increment per OpDet counters and fill per phot trees
                  fCountOpDetAll++;
                  if(fMakeAllPhotonsTree){
                    if (!isVisible(fWavelength) || fPVS->StoreReflected()) {
                        fThePhotonTreeAll->Fill();
                    }
                  }

                  if(odresponse->detected(fOpChannel, Phot))
                  {
                    if(fMakeDetectedPhotonsTree) fThePhotonTreeDetected->Fill();
                    //only store direct direct light
                    if(!isVisible(fWavelength))
                      fCountOpDetDetected++;
                    // reflected and shifted light is in visible range
                    else if(fPVS->StoreReflected()) {
                      fCountOpDetReflDetected++;
                      // find the first visible arrival time
                      if(fPVS->StoreReflT0() && fTime < fT0_vis)
                        fT0_vis = fTime;
                    }
                    if(fVerbosity > 3)
                      std::cout<<"OpDetResponseInterface PerPhoton : Event "<<fEventID<<" OpChannel " <<fOpChannel << " Wavelength " << fWavelength << " Detected 1 "<<std::endl;
                  }
                  else {
                    if(fVerbosity > 3)
                      std::cout<<"OpDetResponseInterface PerPhoton : Event "<<fEventID<<" OpChannel " <<fOpChannel << " Wavelength " << fWavelength << " Detected 0 "<<std::endl;
                  }

                } // if build library and not reflected

                else {
                  // store in appropriate trees using "Reflected" handle and fPVS->StoreReflected() flag
                  // Increment per OpDet counters and fill per phot trees
                  fCountOpDetAll++;
                  if(fMakeAllPhotonsTree){
                    if (!Reflected || (fPVS->StoreReflected() && Reflected)) {
                        fThePhotonTreeAll->Fill();
                    }
                  }

                  if(odresponse->detected(fOpChannel, Phot))
                  {
                    if(fMakeDetectedPhotonsTree) fThePhotonTreeDetected->Fill();
                    //only store direct direct light
                    if(!Reflected)
                      fCountOpDetDetected++;
                    // reflected and shifted light is in visible range
                    else if(fPVS->StoreReflected() && Reflected ) {
                      fCountOpDetReflDetected++;
                      // find the first visible arrival time
                      if(fPVS->StoreReflT0() && fTime < fT0_vis)
                        fT0_vis = fTime;
                    }
                    if(fVerbosity > 3)
                      std::cout<<"OpDetResponseInterface PerPhoton : Event "<<fEventID<<" OpChannel " <<fOpChannel << " Wavelength " << fWavelength << " Detected 1 "<<std::endl;
                  }
                  else {
                    if(fVerbosity > 3)
                      std::cout<<"OpDetResponseInterface PerPhoton : Event "<<fEventID<<" OpChannel " <<fOpChannel << " Wavelength " << fWavelength << " Detected 0 "<<std::endl;
                  }
                }
              } // for each photon in collection



              // If this is a library building job, fill relevant entry
              if(fPVS->IsBuildJob() && !Reflected) // for library build job, both componenents stored in first object with Reflected = false
              {
                storeVisibility(fOpChannel, fCountOpDetDetected, fCountOpDetReflDetected, fT0_vis);
              }

              // Incremenent per event and fill Per OpDet trees
              if(fMakeOpDetsTree) fTheOpDetTree->Fill();
              fCountEventAll+=fCountOpDetAll;
              fCountEventDetected+=fCountOpDetDetected;

              // Give per OpDet output
              if(fVerbosity >2) std::cout<<"OpDetResponseInterface PerOpDet : Event "<<fEventID<<" OpDet " << fOpChannel << " All " << fCountOpDetAll << " Det " <<fCountOpDetDetected<<std::endl;
            }

            // Fill per event tree
            if(fMakeOpDetEventsTree) fTheEventTree->Fill();

            // Give per event output
            if(fVerbosity >1) std::cout<<"OpDetResponseInterface PerEvent : Event "<<fEventID<<" All " << fCountOpDetAll << " Det " <<fCountOpDetDetected<<std::endl;

          }
          else
          {
            // if empty OpDet hit collection,
            // add an empty record to the per event tree
            if(fMakeOpDetEventsTree) fTheEventTree->Fill();
          }
          if(fMakeLightAnalysisTree) {
            assert(mcpartVec);
            assert(fLightAnalysisTree);

            std::cout<<"Filling the analysis tree"<<std::endl;
            //---------------Filling the analysis tree-----------:
            fRun = evt.run();
            std::vector<double> this_xyz;

            //loop over the particles
            for(simb::MCParticle const& pPart: *mcpartVec){

              if(pPart.Process() == "primary")
                fEnergy = pPart.E();

              //resetting the vectors
              fstepPositions.clear();
              fstepTimes.clear();
              fSignalsvuv.clear();
              fSignalsvis.clear();
              fdEdx = -1.;
              //filling the tree fields
              fTrackID = pPart.TrackId();
              fpdg = pPart.PdgCode();
              fmotherTrackID = pPart.Mother();
              fdEdx = totalEnergy_track[fTrackID];
              fSignalsvuv = fSignals_vuv[fTrackID];
              fSignalsvis = fSignals_vis[fTrackID];
              fProcess = pPart.Process();
              //filling the center positions of each step
              for(size_t i_s=1; i_s < pPart.NumberTrajectoryPoints(); i_s++){
                this_xyz.clear();
                this_xyz.resize(3);
                this_xyz[0] = pPart.Position(i_s).X();
                this_xyz[1] = pPart.Position(i_s).Y();
                this_xyz[2] = pPart.Position(i_s).Z();
                fstepPositions.push_back(this_xyz);
                fstepTimes.push_back(pPart.Position(i_s).T());
              }
              //filling the tree per track
              fLightAnalysisTree->Fill();
            }
          } // if fMakeLightAnalysisTree
        }
      }
    }
    if (fUseLitePhotons)
    {

      //Get *ALL* SimPhotonsCollection from Event
      //std::vector< art::Handle< std::vector< sim::SimPhotonsLite > > > photon_handles;
      //evt.getManyByType(photon_handles);
      auto photon_handles = evt.getMany<std::vector<sim::SimPhotonsLite>>();
      if (photon_handles.size() == 0)
        throw art::Exception(art::errors::ProductNotFound)<<"sim SimPhotons retrieved and you requested them.";

     //Get SimPhotonsLite from Event
     for(auto const& mod : fInputModule){
      //art::Handle< std::vector<sim::SimPhotonsLite> > photonHandle;
      //evt.getByLabel(mod, photonHandle);

        // Loop over direct/reflected photons
        for (auto const& ph_handle: photon_handles) {
          // Do some checking before we proceed
          if (!ph_handle.isValid()) continue;
          if (ph_handle.provenance()->moduleLabel() != mod) continue;   //not the most efficient way of doing this, but preserves the logic of the module. Andrzej

          bool Reflected = (ph_handle.provenance()->productInstanceName() == "Reflected");

          //Reset counters
          fCountEventAll=0;
          fCountEventDetected=0;

          if(fVerbosity > 0) std::cout<<"Found OpDet hit collection of size "<< (*ph_handle).size()<<std::endl;


          if((*ph_handle).size()>0)
          {

            for ( auto const& photon : (*ph_handle) )
            {
              //Get data from HitCollection entry
              fOpChannel=photon.OpChannel;
              std::map<int, int> PhotonsMap = photon.DetectedPhotons;

              //Reset Counters
              fCountOpDetAll=0;
              fCountOpDetDetected=0;
              fCountOpDetReflDetected=0;

              for(auto it = PhotonsMap.begin(); it!= PhotonsMap.end(); it++)
              {
                // Calculate wavelength in nm
                if (Reflected) {
                  fWavelength = kVisibleThreshold;
                }
                else {
                  fWavelength= kVUVWavelength;   // original
                }

                //Get arrival time from phot
                fTime= it->first;
                //std::cout<<"Arrival time: " << fTime<<std::endl;

                for(int i = 0; i < it->second ; i++)
                {
                  // Increment per OpDet counters and fill per phot trees
                  fCountOpDetAll++;
                  if(fMakeAllPhotonsTree) fThePhotonTreeAll->Fill();

                  if(odresponse->detectedLite(fOpChannel))
                  {
                    if(fMakeDetectedPhotonsTree) fThePhotonTreeDetected->Fill();
                    // direct light
                    if (!Reflected){
                      fCountOpDetDetected++;
                    }
                    else if (Reflected) {
                      fCountOpDetReflDetected++;
                    }
                    if(fVerbosity > 3)
                    std::cout<<"OpDetResponseInterface PerPhoton : Event "<<fEventID<<" OpChannel " <<fOpChannel << " Wavelength " << fWavelength << " Detected 1 "<<std::endl;
                  }
                  else
                    if(fVerbosity > 3)
                    {
                    std::cout<<"OpDetResponseInterface PerPhoton : Event "<<fEventID<<" OpChannel " <<fOpChannel << " Wavelength " << fWavelength << " Detected 0 "<<std::endl;
                    }
                }
              }



              // Incremenent per event and fill Per OpDet trees
              if(fMakeOpDetsTree) fTheOpDetTree->Fill();
              fCountEventAll+=fCountOpDetAll;
              fCountEventDetected+=fCountOpDetDetected;

              if (fPVS->IsBuildJob())
                storeVisibility(fOpChannel, fCountOpDetDetected, fCountOpDetReflDetected, fT0_vis);


              // Give per OpDet output
              if(fVerbosity >2) std::cout<<"OpDetResponseInterface PerOpDet : Event "<<fEventID<<" OpDet " << fOpChannel << " All " << fCountOpDetAll << " Det " <<fCountOpDetDetected<<std::endl;
            }
            // Fill per event tree
            if(fMakeOpDetEventsTree) fTheEventTree->Fill();

            // Give per event output
            if(fVerbosity >1) std::cout<<"OpDetResponseInterface PerEvent : Event "<<fEventID<<" All " << fCountOpDetAll << " Det " <<fCountOpDetDetected<<std::endl;

          }
          else
          {
            // if empty OpDet hit collection,
            // add an empty record to the per event tree
            if(fMakeOpDetEventsTree) fTheEventTree->Fill();
          }
        }
      }
    }
  } // SimPhotonCounter::analyze()

void opdet::SimPhotonCounter::beginJob

(

)

Definition at line 205 of file SimPhotonCounter_module.cc.

  {
    // Get file service to store trees
    art::ServiceHandle<art::TFileService const> tfs;
    art::ServiceHandle<geo::Geometry const> geo;

    std::cout<<"Optical Channels positions:  "<<geo->Cryostat(0).NOpDet()<<std::endl;
    for(int ch=0; ch!=int(geo->Cryostat(0).NOpDet()); ch++) {
      double OpDetCenter[3];
      geo->OpDetGeoFromOpDet(ch).GetCenter(OpDetCenter);
      std::cout<<ch<<"  "<<OpDetCenter[0]<<"  "<<OpDetCenter[1]<<"  "<<OpDetCenter[2]<<std::endl;
    }

    double CryoBounds[6];
    geo->CryostatBoundaries(CryoBounds);
    std::cout<<"Cryo Boundaries"<<std::endl;
    std::cout<<"Xmin: "<<CryoBounds[0]<<" Xmax: "<<CryoBounds[1]<<" Ymin: "<<CryoBounds[2]
             <<" Ymax: "<<CryoBounds[3]<<" Zmin: "<<CryoBounds[4]<<" Zmax: "<<CryoBounds[5]<<std::endl;

    try {
      pi_serv = &*(art::ServiceHandle<cheat::ParticleInventoryService const>());
    }
    catch (art::Exception const& e) {
      if (e.categoryCode() != art::errors::ServiceNotFound) throw;
      mf::LogError("SimPhotonCounter")
        << "ParticleInventoryService service is not configured!"
        " Please add it in the job configuration."
        " In the meanwhile, some checks to particles will be skipped."
        ;
    }

    // Create and assign branch addresses to required tree
    if(fMakeAllPhotonsTree)
    {
      fThePhotonTreeAll = tfs->make<TTree>("AllPhotons","AllPhotons");
      fThePhotonTreeAll->Branch("EventID",     &fEventID,          "EventID/I");
      fThePhotonTreeAll->Branch("Wavelength",  &fWavelength,       "Wavelength/F");
      fThePhotonTreeAll->Branch("OpChannel",       &fOpChannel,            "OpChannel/I");
      fThePhotonTreeAll->Branch("Time",        &fTime,             "Time/F");
    }

    if(fMakeDetectedPhotonsTree)
    {
      fThePhotonTreeDetected = tfs->make<TTree>("DetectedPhotons","DetectedPhotons");
      fThePhotonTreeDetected->Branch("EventID",     &fEventID,          "EventID/I");
      fThePhotonTreeDetected->Branch("Wavelength",  &fWavelength,       "Wavelength/F");
      fThePhotonTreeDetected->Branch("OpChannel",       &fOpChannel,            "OpChannel/I");
      fThePhotonTreeDetected->Branch("Time",        &fTime,             "Time/F");
    }

    if(fMakeOpDetsTree)
    {
      fTheOpDetTree    = tfs->make<TTree>("OpDets","OpDets");
      fTheOpDetTree->Branch("EventID",        &fEventID,          "EventID/I");
      fTheOpDetTree->Branch("OpChannel",          &fOpChannel,            "OpChannel/I");
      fTheOpDetTree->Branch("CountAll",       &fCountOpDetAll,      "CountAll/I");
      fTheOpDetTree->Branch("CountDetected",  &fCountOpDetDetected, "CountDetected/I");
      if(fPVS->StoreReflected())
        fTheOpDetTree->Branch("CountReflDetected",  &fCountOpDetReflDetected, "CountReflDetected/I");
      fTheOpDetTree->Branch("Time",                   &fTime,                          "Time/F");
    }

    if(fMakeOpDetEventsTree)
    {
      fTheEventTree  = tfs->make<TTree>("OpDetEvents","OpDetEvents");
      fTheEventTree->Branch("EventID",      &fEventID,            "EventID/I");
      fTheEventTree->Branch("CountAll",     &fCountEventAll,     "CountAll/I");
      fTheEventTree->Branch("CountDetected",&fCountEventDetected,"CountDetected/I");
      if(fPVS->StoreReflected())
        fTheOpDetTree->Branch("CountReflDetected",  &fCountOpDetReflDetected, "CountReflDetected/I");

    }

    //generating the tree for the light analysis:
    if(fMakeLightAnalysisTree)
    {
      fLightAnalysisTree = tfs->make<TTree>("LightAnalysis","LightAnalysis");
      fLightAnalysisTree->Branch("RunNumber",&fRun);
      fLightAnalysisTree->Branch("EventID",&fEventID);
      fLightAnalysisTree->Branch("TrackID",&fTrackID);
      fLightAnalysisTree->Branch("PdgCode",&fpdg);
      fLightAnalysisTree->Branch("MotherTrackID",&fmotherTrackID);
      fLightAnalysisTree->Branch("Energy",&fEnergy);
      fLightAnalysisTree->Branch("dEdx",&fdEdx);
      fLightAnalysisTree->Branch("StepPositions",&fstepPositions);
      fLightAnalysisTree->Branch("StepTimes",&fstepTimes);
      fLightAnalysisTree->Branch("SignalsVUV",&fSignalsvuv);
      fLightAnalysisTree->Branch("SignalsVisible",&fSignalsvis);
      fLightAnalysisTree->Branch("Process",&fProcess);
    }

  }

void opdet::SimPhotonCounter::ClearVectors

(

)

Definition at line 86 of file SimPhotonCounter.cxx.

{
  for(size_t i=0; i<GetVectorSize(); i++){
    _photonVector_prompt[i]=0.0;
    _photonVector_late[i]=0.0;
  }
}

void opdet::SimPhotonCounter::endJob

(

)

Definition at line 299 of file SimPhotonCounter_module.cc.

  {
    if(fPVS->IsBuildJob())
    {
      art::ServiceHandle<phot::PhotonVisibilityService>()->StoreLibrary();
    }
  }

size_t opdet::SimPhotonCounter::GetVectorSize ( ) const

inline

Definition at line 32 of file SimPhotonCounter.h.

{ return _photonVector_prompt.size(); }

bool opdet::SimPhotonCounter::isVisible ( double  wavelength ) const

inlineprivate

Returns if we label as "visibile" a photon with specified wavelength [nm].

Definition at line 163 of file SimPhotonCounter_module.cc.

        { return fWavelength < kVisibleThreshold; }

float opdet::SimPhotonCounter::LatePhotonTotal ( ) const

inline

Definition at line 68 of file SimPhotonCounter.h.

    { return std::accumulate(_photonVector_late.begin(),_photonVector_late.end(),0.0); }

const std::vector<float>& opdet::SimPhotonCounter::LatePhotonVector ( ) const

inline

Definition at line 58 of file SimPhotonCounter.h.

{ return _photonVector_late; }

float opdet::SimPhotonCounter::LatePhotonVector ( size_t  i ) const

inline

Definition at line 60 of file SimPhotonCounter.h.

{ return _photonVector_late.at(i); }

float opdet::SimPhotonCounter::MaxLateTime ( ) const

inline

Definition at line 44 of file SimPhotonCounter.h.

{ return _max_late_time; }

float opdet::SimPhotonCounter::MaxPromptTime ( ) const

inline

Definition at line 42 of file SimPhotonCounter.h.

{ return _max_prompt_time; }

float opdet::SimPhotonCounter::MaxWavelength ( ) const

inline

Definition at line 37 of file SimPhotonCounter.h.

{ return _max_wavelength; }

float opdet::SimPhotonCounter::MinLateTime ( ) const

inline

Definition at line 43 of file SimPhotonCounter.h.

{ return _min_late_time; }

float opdet::SimPhotonCounter::MinPromptTime ( ) const

inline

Definition at line 41 of file SimPhotonCounter.h.

{ return _min_prompt_time; }

float opdet::SimPhotonCounter::MinWavelength ( ) const

inline

Definition at line 36 of file SimPhotonCounter.h.

{ return _min_wavelength; }

float opdet::SimPhotonCounter::PhotonTotal ( ) const

inline

Definition at line 70 of file SimPhotonCounter.h.

{ return (PromptPhotonTotal()+LatePhotonTotal()); }

void opdet::SimPhotonCounter::Print

(

)

Definition at line 104 of file SimPhotonCounter.cxx.

{
  std::cout << "Vector size: " << GetVectorSize() << std::endl;
  std::cout << "Time cut ranges: ("
            << MinPromptTime() << "," << MaxPromptTime() << ") , ("
            << MinLateTime() << "," << MaxLateTime() << ")" << std::endl;
  std::cout << "\t" << "i : QE / Prompt / Late / Total" << std::endl;
  for(size_t i=0; i<GetVectorSize(); i++)
    std::cout << "\t" << i << ": " << _qeVector[i] << " / " << _photonVector_prompt[i] << " / " << _photonVector_late[i] << " / " << TotalPhotonVector(i) << std::endl;

}

float opdet::SimPhotonCounter::PromptPhotonTotal ( ) const

inline

Definition at line 66 of file SimPhotonCounter.h.

    { return std::accumulate(_photonVector_prompt.begin(),_photonVector_prompt.end(),0.0); }

const std::vector<float>& opdet::SimPhotonCounter::PromptPhotonVector ( ) const

inline

Definition at line 57 of file SimPhotonCounter.h.

{ return _photonVector_prompt; }

float opdet::SimPhotonCounter::PromptPhotonVector ( size_t  i ) const

inline

Definition at line 59 of file SimPhotonCounter.h.

{ return _photonVector_prompt.at(i); }

float opdet::SimPhotonCounter::QE ( size_t  i ) const

inline

Definition at line 47 of file SimPhotonCounter.h.

{ return _qeVector.at(i); }

std::vector<float> const& opdet::SimPhotonCounter::QEVector ( ) const

inline

Definition at line 51 of file SimPhotonCounter.h.

{ return _qeVector; }

void opdet::SimPhotonCounter::SetQE ( size_t  i, float  e  )

inline

Definition at line 46 of file SimPhotonCounter.h.

{ _qeVector.at(i) = e; }

void opdet::SimPhotonCounter::SetQEVector ( const std::vector< float > &  eV )

inline

Definition at line 49 of file SimPhotonCounter.h.

    { SetVectorSize(eV.size()); _qeVector = eV; }

void opdet::SimPhotonCounter::SetTimeRanges	(	float	t_p1,
		float	t_p2,
		float	t_l1,
		float	t_l2
	)

Definition at line 57 of file SimPhotonCounter.cxx.

{
  if(t_p2<t_p1 || t_l2<t_l1 || t_p2>t_l1 )
    throw std::runtime_error("ERROR in SimPhotonCounter: bad time ranges");

  _min_prompt_time = t_p1; _max_prompt_time = t_p2;
  _min_late_time = t_l1; _max_late_time = t_l2;
}

void opdet::SimPhotonCounter::SetVectorSize ( size_t  s )

inline

Definition at line 30 of file SimPhotonCounter.h.

    { _photonVector_prompt.resize(s); _photonVector_late.resize(s); _qeVector.resize(s); }

void opdet::SimPhotonCounter::SetWavelengthRanges	(	float	min_w,
		float	max_w
	)

Definition at line 40 of file SimPhotonCounter.cxx.

{
  if(min_w >= max_w)
    throw std::runtime_error("ERROR in SimPhotonCounter: bad wavelength range");

  _min_wavelength = min_w;
  _max_wavelength = max_w;
}

void opdet::SimPhotonCounter::storeVisibility ( int  channel, int  nDirectPhotons, int  nReflectedPhotons, double  reflectedT0 = 0.0  ) const

private

Definition at line 723 of file SimPhotonCounter_module.cc.

  {
    phot::PhotonVisibilityService& pvs
      = *(art::ServiceHandle<phot::PhotonVisibilityService>());

    // ask PhotonVisibilityService which voxel was being served,
    // and how many photons where there generated (yikes!!);
    // this value was put there by LightSource (yikes!!)
    int VoxID;
    double NProd;
    fPVS->RetrieveLightProd(VoxID, NProd);

    pvs.SetLibraryEntry(VoxID, channel, double(nDirectPhotons)/NProd);

    //store reflected light
    if(fPVS->StoreReflected()) {
      pvs.SetLibraryEntry(VoxID, channel, double(nReflectedPhotons)/NProd, true);

      //store reflected first arrival time
      if (fPVS->StoreReflT0())
        pvs.SetLibraryReflT0Entry(VoxID, channel, reflectedT0);

    } // if reflected

  } // SimPhotonCounter::storeVisibility()

std::vector< float > opdet::SimPhotonCounter::TotalPhotonVector

(

)

const

Definition at line 94 of file SimPhotonCounter.cxx.

                                                               {

  std::vector<float> totalPhotonVector(GetVectorSize());
  std::transform(PromptPhotonVector().begin(),PromptPhotonVector().end(),
                 LatePhotonVector().begin(),
                 totalPhotonVector.begin(),
                 std::plus<float>());
  return totalPhotonVector;
}

float opdet::SimPhotonCounter::TotalPhotonVector ( size_t  i ) const

inline

Definition at line 63 of file SimPhotonCounter.h.

    { return (PromptPhotonVector(i)+LatePhotonVector(i)); }

float opdet::SimPhotonCounter::Wavelength ( const sim::OnePhoton &  ph )

private

Definition at line 49 of file SimPhotonCounter.cxx.

{
  if(ph.Energy < std::numeric_limits<float>::epsilon())
    throw std::runtime_error("ERROR in SimPhotonCounter: photon energy is zero.");

  return 0.00124 / ph.Energy;
}

Member Data Documentation

float opdet::SimPhotonCounter::_max_late_time

private

Definition at line 82 of file SimPhotonCounter.h.

float opdet::SimPhotonCounter::_max_prompt_time

private

Definition at line 80 of file SimPhotonCounter.h.

float opdet::SimPhotonCounter::_max_wavelength

private

Definition at line 86 of file SimPhotonCounter.h.

float opdet::SimPhotonCounter::_min_late_time

private

Definition at line 81 of file SimPhotonCounter.h.

float opdet::SimPhotonCounter::_min_prompt_time

private

Definition at line 79 of file SimPhotonCounter.h.

float opdet::SimPhotonCounter::_min_wavelength

private

Definition at line 85 of file SimPhotonCounter.h.

std::vector<float> opdet::SimPhotonCounter::_photonVector_late

private

Definition at line 77 of file SimPhotonCounter.h.

std::vector<float> opdet::SimPhotonCounter::_photonVector_prompt

private

Definition at line 76 of file SimPhotonCounter.h.

std::vector<float> opdet::SimPhotonCounter::_qeVector

private

Definition at line 83 of file SimPhotonCounter.h.

Int_t opdet::SimPhotonCounter::fCountEventAll

private

Definition at line 129 of file SimPhotonCounter_module.cc.

Int_t opdet::SimPhotonCounter::fCountEventDetected

private

Definition at line 130 of file SimPhotonCounter_module.cc.

Int_t opdet::SimPhotonCounter::fCountOpDetAll

private

Definition at line 124 of file SimPhotonCounter_module.cc.

Int_t opdet::SimPhotonCounter::fCountOpDetDetected

private

Definition at line 125 of file SimPhotonCounter_module.cc.

Int_t opdet::SimPhotonCounter::fCountOpDetReflDetected

private

Definition at line 126 of file SimPhotonCounter_module.cc.

double opdet::SimPhotonCounter::fdEdx

private

Definition at line 143 of file SimPhotonCounter_module.cc.

double opdet::SimPhotonCounter::fEnergy

private

Definition at line 143 of file SimPhotonCounter_module.cc.

Int_t opdet::SimPhotonCounter::fEventID

private

Definition at line 133 of file SimPhotonCounter_module.cc.

TVector3 opdet::SimPhotonCounter::finalPhotonPosition

private

Definition at line 116 of file SimPhotonCounter_module.cc.

std::vector<std::string> opdet::SimPhotonCounter::fInputModule

private

Definition at line 101 of file SimPhotonCounter_module.cc.

TTree* opdet::SimPhotonCounter::fLightAnalysisTree = nullptr

private

Definition at line 141 of file SimPhotonCounter_module.cc.

bool opdet::SimPhotonCounter::fMakeAllPhotonsTree

private

Definition at line 106 of file SimPhotonCounter_module.cc.

bool opdet::SimPhotonCounter::fMakeDetectedPhotonsTree

private

Definition at line 105 of file SimPhotonCounter_module.cc.

bool opdet::SimPhotonCounter::fMakeLightAnalysisTree

private

Definition at line 137 of file SimPhotonCounter_module.cc.

bool opdet::SimPhotonCounter::fMakeOpDetEventsTree

private

Definition at line 108 of file SimPhotonCounter_module.cc.

bool opdet::SimPhotonCounter::fMakeOpDetsTree

private

Definition at line 107 of file SimPhotonCounter_module.cc.

int opdet::SimPhotonCounter::fmotherTrackID

private

Definition at line 142 of file SimPhotonCounter_module.cc.

Int_t opdet::SimPhotonCounter::fOpChannel

private

Definition at line 134 of file SimPhotonCounter_module.cc.

int opdet::SimPhotonCounter::fpdg

private

Definition at line 142 of file SimPhotonCounter_module.cc.

std::vector<double> opdet::SimPhotonCounter::fPosition0

private

Definition at line 144 of file SimPhotonCounter_module.cc.

std::string opdet::SimPhotonCounter::fProcess

private

Definition at line 149 of file SimPhotonCounter_module.cc.

phot::PhotonVisibilityService const* opdet::SimPhotonCounter::fPVS = nullptr

private

Definition at line 152 of file SimPhotonCounter_module.cc.

int opdet::SimPhotonCounter::fRun

private

Definition at line 142 of file SimPhotonCounter_module.cc.

std::vector<std::vector<std::vector<double> > > opdet::SimPhotonCounter::fSignals_vis

private

Definition at line 139 of file SimPhotonCounter_module.cc.

std::vector<std::vector<std::vector<double> > > opdet::SimPhotonCounter::fSignals_vuv

private

Definition at line 138 of file SimPhotonCounter_module.cc.

std::vector<std::vector<double> > opdet::SimPhotonCounter::fSignalsvis

private

Definition at line 148 of file SimPhotonCounter_module.cc.

std::vector<std::vector<double> > opdet::SimPhotonCounter::fSignalsvuv

private

Definition at line 147 of file SimPhotonCounter_module.cc.

std::vector<std::vector<double> > opdet::SimPhotonCounter::fstepPositions

private

Definition at line 145 of file SimPhotonCounter_module.cc.

std::vector<double> opdet::SimPhotonCounter::fstepTimes

private

Definition at line 146 of file SimPhotonCounter_module.cc.

Float_t opdet::SimPhotonCounter::fT0_vis

private

Definition at line 127 of file SimPhotonCounter_module.cc.

TTree* opdet::SimPhotonCounter::fTheEventTree

private

Definition at line 96 of file SimPhotonCounter_module.cc.

TTree* opdet::SimPhotonCounter::fTheOpDetTree

private

Definition at line 95 of file SimPhotonCounter_module.cc.

TTree* opdet::SimPhotonCounter::fThePhotonTreeAll

private

Definition at line 93 of file SimPhotonCounter_module.cc.

TTree* opdet::SimPhotonCounter::fThePhotonTreeDetected

private

Definition at line 94 of file SimPhotonCounter_module.cc.

Float_t opdet::SimPhotonCounter::fTime

private

Definition at line 122 of file SimPhotonCounter_module.cc.

int opdet::SimPhotonCounter::fTrackID

private

Definition at line 142 of file SimPhotonCounter_module.cc.

bool const opdet::SimPhotonCounter::fUseLitePhotons

private

Definition at line 154 of file SimPhotonCounter_module.cc.

int opdet::SimPhotonCounter::fVerbosity

private

Definition at line 103 of file SimPhotonCounter_module.cc.

Float_t opdet::SimPhotonCounter::fWavelength

private

Definition at line 121 of file SimPhotonCounter_module.cc.

TVector3 opdet::SimPhotonCounter::initialPhotonPosition

private

Definition at line 115 of file SimPhotonCounter_module.cc.

constexpr double opdet::SimPhotonCounter::kVisibleThreshold = 200.0

staticprivate

Threshold used to resolve between visible and ultraviolet light.

Definition at line 83 of file SimPhotonCounter_module.cc.

constexpr double opdet::SimPhotonCounter::kVisibleWavelength = 450.0

staticprivate

Value used when a typical visible light wavelength is needed.

Definition at line 86 of file SimPhotonCounter_module.cc.

constexpr double opdet::SimPhotonCounter::kVUVWavelength = 128.0

staticprivate

Value used when a typical ultraviolet light wavelength is needed.

Definition at line 89 of file SimPhotonCounter_module.cc.

cheat::ParticleInventoryService const* opdet::SimPhotonCounter::pi_serv = nullptr

private

Definition at line 151 of file SimPhotonCounter_module.cc.

---

The documentation for this class was generated from the following files:

• SimPhotonCounter.h
• SimPhotonCounter_module.cc
• SimPhotonCounter.cxx