#include <CastorShowerLibrary.h>

Public Member Functions
	CastorShowerLibrary (std::string &name, edm::ParameterSet const &p)

int	FindEnergyBin (double)

int	FindEtaBin (double)

int	FindPhiBin (double)

CastorShowerEvent	getShowerHits (G4Step *, bool &)

void	initParticleTable (G4ParticleTable *)

	~CastorShowerLibrary ()

Protected Member Functions
void	getRecord (int, int)

void	initFile (edm::ParameterSet const &)

void	loadEventInfo (TBranchObject *)

void	select (int, double, double=0, double=9)

Private Attributes
int	anuePDG

int	anumuPDG

int	anutauPDG

TBranchObject *	emBranch

int	emPDG

int	epPDG

int	etaPDG

CastorShowerLibraryInfo *	eventInfo

TBranchObject *	evtInfo

unsigned int	evtPerBin

int	gammaPDG

int	geantinoPDG

TBranchObject *	hadBranch

TFile *	hf

int	mumPDG

int	mupPDG

unsigned int	nBinsE

unsigned int	nBinsEta

unsigned int	nBinsPhi

unsigned int	nEvtPerBinE

unsigned int	nEvtPerBinEta

unsigned int	nEvtPerBinPhi

unsigned int	nMomBin

int	nuePDG

int	numuPDG

int	nutauPDG

int	pi0PDG

std::vector< double >	pmom

CastorShowerEvent *	showerEvent

std::vector< double >	SLenergies

std::vector< double >	SLetas

std::vector< double >	SLphis

unsigned int	totEvents

bool	verbose

Detailed Description

Definition at line 32 of file CastorShowerLibrary.h.

Constructor & Destructor Documentation

CastorShowerLibrary::CastorShowerLibrary	(	std::string &	name,
		edm::ParameterSet const &	p
	)

Definition at line 22 of file CastorShowerLibrary.cc.

References initFile().

                                           : hf(0), evtInfo(0), emBranch(0), hadBranch(0),
                                             nMomBin(0), totEvents(0), evtPerBin(0),
                                             nBinsE(0),nBinsEta(0),nBinsPhi(0),
                                             nEvtPerBinE(0),nEvtPerBinEta(0),nEvtPerBinPhi(0),
                                             SLenergies(),SLetas(),SLphis() {
   
   initFile(p);
   
 }

CastorShowerLibrary::~CastorShowerLibrary ( )

Definition at line 36 of file CastorShowerLibrary.cc.

References hf.

                                           {
   if (hf)     hf->Close();
 }

Member Function Documentation

int CastorShowerLibrary::FindEnergyBin ( double energy )

Definition at line 399 of file CastorShowerLibrary.cc.

References i, and SLenergies.

Referenced by select().

                                                     {
   //
   // returns the integer index of the energy bin, taken from SLenergies vector
   // returns -1 if ouside valid range
   //
   if (energy >= SLenergies.back()) return SLenergies.size()-1;
 
   unsigned int i = 0;
   for(;i<SLenergies.size()-1;i++)
     if (energy >= SLenergies.at(i) && energy < SLenergies.at(i+1)) return (int)i;
 
   // now i points to the last but 1 bin
   if (energy>=SLenergies.at(i)) return (int)i;
   // energy outside bin range
   return -1;
 }

int CastorShowerLibrary::FindEtaBin ( double eta )

Definition at line 415 of file CastorShowerLibrary.cc.

References i, and SLetas.

Referenced by select().

                                               {
   //
   // returns the integer index of the eta bin, taken from SLetas vector
   // returns -1 if ouside valid range
   //
   if (eta>=SLetas.back()) return SLetas.size()-1;
   unsigned int i = 0;
   for(;i<SLetas.size()-1;i++)
      if (eta >= SLetas.at(i) && eta < SLetas.at(i+1)) return (int)i;
   // now i points to the last but 1 bin
   if (eta>=SLetas.at(i)) return (int)i;
   // eta outside bin range
   return -1;
 }

int CastorShowerLibrary::FindPhiBin ( double phi )

Definition at line 429 of file CastorShowerLibrary.cc.

References i, and SLphis.

Referenced by select().

                                               {
   //
   // returns the integer index of the phi bin, taken from SLphis vector
   // returns -1 if ouside valid range
   //
   // needs protection in case phi is outside range -pi,pi
   //
   if (phi>=SLphis.back()) return SLphis.size()-1;
   unsigned int i = 0;
   for(;i<SLphis.size()-1;i++)
      if (phi >= SLphis.at(i) && phi < SLphis.at(i+1)) return (int)i;
   // now i points to the last but 1 bin
   if (phi>=SLphis.at(i)) return (int)i;
   // phi outside bin range
   return -1;
 }

void CastorShowerLibrary::getRecord	(	int	type,
		int	record
	)

protected

Definition at line 272 of file CastorShowerLibrary.cc.

References emBranch, CastorShowerEvent::getNhit(), hadBranch, LogDebug, record, and showerEvent.

Referenced by select().

                                                         {
 //
 //  Retrieve event # "record" from the library and stores it   
 //  into  vector<CastorShowerHit> showerHit;
 //
 //  Based on HFShowerLibrary::getRecord
 //
 //  Modified 02/02/09 by W. Carvalho
 //
 
 #ifdef DebugLog
   LogDebug("CastorShower") << "CastorShowerLibrary::getRecord: ";
 #endif  
   int nrc  = record;
   int nHit = 0;
   showerEvent = new CastorShowerEvent();
   if (type > 0) {
     hadBranch->SetAddress(&showerEvent);
     hadBranch->GetEntry(nrc);
   } else {
     emBranch->SetAddress(&showerEvent);
     emBranch->GetEntry(nrc);
   }
   nHit = showerEvent->getNhit();
 
 #ifdef DebugLog
   LogDebug("CastorShower") << "CastorShowerLibrary::getRecord: Record " << record
                    << " of type " << type << " with " << nHit 
                    << " CastorShowerHits";
 
 #endif
 
 }

CastorShowerEvent CastorShowerLibrary::getShowerHits	(	G4Step *	aStep,
		bool &	ok
	)

Definition at line 211 of file CastorShowerLibrary.cc.

References abs, anuePDG, anumuPDG, anutauPDG, CastorShowerEvent::Clear(), emPDG, epPDG, etaPDG, gammaPDG, geantinoPDG, funct::log(), M_PI, nuePDG, numuPDG, nutauPDG, pi0PDG, dttmaxenums::R, select(), mathSSE::sqrt(), funct::tan(), and theta().

Referenced by CastorSD::getFromLibrary().

                                                                               {
 
   G4StepPoint * preStepPoint  = aStep->GetPreStepPoint(); 
   // G4StepPoint * postStepPoint = aStep->GetPostStepPoint(); 
   G4Track *     track         = aStep->GetTrack();
   // Get Z-direction 
   const G4DynamicParticle *aParticle = track->GetDynamicParticle();
   G4ThreeVector               momDir = aParticle->GetMomentumDirection();
   //  double mom = aParticle->GetTotalMomentum();
 
   G4ThreeVector hitPoint = preStepPoint->GetPosition();   
   G4String      partType = track->GetDefinition()->GetParticleName();
   int           parCode  = track->GetDefinition()->GetPDGEncoding();
 
   CastorShowerEvent hit;
   hit.Clear();
   
   ok = false;
   if (parCode == pi0PDG   || parCode == etaPDG    || parCode == nuePDG  ||
       parCode == numuPDG  || parCode == nutauPDG  || parCode == anuePDG ||
       parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG) 
     return hit;
   ok = true;
 
   double pin    = preStepPoint->GetTotalEnergy();
   double etain  = momDir.getEta();
   double phiin  = momDir.getPhi();
   
   double zint = hitPoint.z();
   double R=sqrt(hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
   double theta = atan2(R,std::abs(zint));
   phiin = atan2(hitPoint.y(),hitPoint.x());
   etain = -1*(std::log(std::tan((M_PI-theta)/2)));
 
   // Replace "interpolation/extrapolation" by new method "select" that just randomly 
   // selects a record from the appropriate energy bin and fills its content to  
   // "showerEvent" instance of class CastorShowerEvent
   
   if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
     select(0, pin, etain, phiin);
     // if (pin<SLenergies[nBinsE-1]) {
     //   interpolate(0, pin);
     // } else {
     //   extrapolate(0, pin);
     // }
   } else {
     select(1, pin, etain, phiin);
     // if (pin<SLenergies[nBinsE-1]) {
     //   interpolate(1, pin);
     // } else {
     //   extrapolate(1, pin);
     // }
   }
     
   hit = (*showerEvent);
   return hit;
 
 }

void CastorShowerLibrary::initFile ( edm::ParameterSet const & p )

protected

Definition at line 43 of file CastorShowerLibrary.cc.

References emBranch, event(), evtInfo, edm::hlt::Exception, edm::FileInPath::fullPath(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hadBranch, hf, and loadEventInfo().

Referenced by CastorShowerLibrary().

                                                             {
 //
 //  Init TFile and associated TBranch's of CASTOR Root file 
 //  holding library events 
 //
 
   //
   //  Read PSet for Castor shower library
   //
 
   edm::ParameterSet m_CS   = p.getParameter<edm::ParameterSet>("CastorShowerLibrary");
   edm::FileInPath fp       = m_CS.getParameter<edm::FileInPath>("FileName");
   std::string pTreeName    = fp.fullPath();
   std::string branchEvInfo = m_CS.getUntrackedParameter<std::string>("BranchEvt");
   std::string branchEM     = m_CS.getUntrackedParameter<std::string>("BranchEM");
   std::string branchHAD    = m_CS.getUntrackedParameter<std::string>("BranchHAD");
   verbose                  = m_CS.getUntrackedParameter<bool>("Verbosity",false);
 
   // Open TFile 
   if (pTreeName.find(".") == 0) pTreeName.erase(0,2);
   const char* nTree = pTreeName.c_str();
   hf                = TFile::Open(nTree);
 
   // Check that TFile has been successfully opened
   if (!hf->IsOpen()) { 
      edm::LogError("CastorShower") << "CastorShowerLibrary: opening " << nTree << " failed";
      throw cms::Exception("Unknown", "CastorShowerLibrary") 
                                    << "Opening of " << pTreeName << " fails\n";
   } else {
      edm::LogInfo("CastorShower")  << "CastorShowerLibrary: opening " << nTree << " successfully"; 
   }
 
   // Check for the TBranch holding EventInfo in "Events" TTree
   TTree* event = (TTree *) hf ->Get("CastorCherenkovPhotons");
   if (event) {
      evtInfo = (TBranchObject *) event->GetBranch(branchEvInfo.c_str());
      if (evtInfo) {
         loadEventInfo(evtInfo);
      } else {
         edm::LogError("CastorShower") << "CastorShowerLibrary: " << branchEvInfo.c_str()
                                       << " Branch does not exit in Event";
         throw cms::Exception("Unknown", "CastorShowerLibrary") << "Event information absent\n";
      }
   } else {
      edm::LogError("CastorShower") << "CastorShowerLibrary: Events Tree does not exist";
      throw cms::Exception("Unknown", "CastorShowerLibrary") << "Events tree absent\n";
   }
   
   // Get EM and HAD Branchs
   emBranch         = (TBranchObject *) event->GetBranch(branchEM.c_str());
   if (verbose) emBranch->Print();
   hadBranch        = (TBranchObject *) event->GetBranch(branchHAD.c_str());
   if (verbose) hadBranch->Print();
   edm::LogInfo("CastorShower") << "CastorShowerLibrary: Branch " << branchEM 
                    << " has " << emBranch->GetEntries() 
                    << " entries and Branch " << branchHAD 
                    << " has " << hadBranch->GetEntries() 
                    << " entries";
 
 }

void CastorShowerLibrary::initParticleTable ( G4ParticleTable * theParticleTable )

Definition at line 162 of file CastorShowerLibrary.cc.

References anuePDG, anumuPDG, anutauPDG, emPDG, epPDG, etaPDG, gammaPDG, geantinoPDG, LogDebug, mumPDG, mupPDG, nuePDG, numuPDG, nutauPDG, and pi0PDG.

Referenced by CastorSD::initRun().

                                                                               {
 //
 //  Set particle codes according to PDG encoding 
 //
 //  Based on HFShowerLibrary::initRun
 //
 
   G4String particleName;
 
   edm::LogInfo("CastorShower") << "CastorShowerLibrary::initParticleTable"
                                << " ***  Accessing PDGEncoding  ***" ;
 
   emPDG       = theParticleTable->FindParticle(particleName="e-")->GetPDGEncoding();
   epPDG       = theParticleTable->FindParticle(particleName="e+")->GetPDGEncoding();
   gammaPDG    = theParticleTable->FindParticle(particleName="gamma")->GetPDGEncoding();
   pi0PDG      = theParticleTable->FindParticle(particleName="pi0")->GetPDGEncoding();
   etaPDG      = theParticleTable->FindParticle(particleName="eta")->GetPDGEncoding();
   nuePDG      = theParticleTable->FindParticle(particleName="nu_e")->GetPDGEncoding();
   numuPDG     = theParticleTable->FindParticle(particleName="nu_mu")->GetPDGEncoding();
   nutauPDG    = theParticleTable->FindParticle(particleName="nu_tau")->GetPDGEncoding();
   anuePDG     = theParticleTable->FindParticle(particleName="anti_nu_e")->GetPDGEncoding();
   anumuPDG    = theParticleTable->FindParticle(particleName="anti_nu_mu")->GetPDGEncoding();
   anutauPDG   = theParticleTable->FindParticle(particleName="anti_nu_tau")->GetPDGEncoding();
   geantinoPDG = theParticleTable->FindParticle(particleName="geantino")->GetPDGEncoding();
   mumPDG      = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
   mupPDG      = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
   
 //#ifdef DebugLog
   LogDebug("CastorShower") << "CastorShowerLibrary: Particle codes for e- = " << emPDG
                    << ", e+ = " << epPDG << ", gamma = " << gammaPDG 
                    << ", pi0 = " << pi0PDG << ", eta = " << etaPDG
                    << ", geantino = " << geantinoPDG << "\n        nu_e = "
                    << nuePDG << ", nu_mu = " << numuPDG << ", nu_tau = "
                    << nutauPDG << ", anti_nu_e = " << anuePDG
                    << ", anti_nu_mu = " << anumuPDG << ", anti_nu_tau = "
                    << anutauPDG << ", mu- = " << mumPDG << ", mu+ = "
                    << mupPDG;
 //#endif
 
     edm::LogInfo("CastorShower") << "  *****   Successfully called:  " 
                                  << "CastorShowerLibrary::initParticleTable()   ***** " ;
 
 }

void CastorShowerLibrary::loadEventInfo ( TBranchObject * branch )

protected

Definition at line 108 of file CastorShowerLibrary.cc.

References CastorShowerLibraryInfo::Energy, CastorShowerLibraryInfo::Eta, eventInfo, SLBin::getBin(), SLBin::getNBins(), SLBin::getNEvtPerBin(), SLBin::getNEvts(), i, nBinsE, nBinsEta, nBinsPhi, nEvtPerBinE, nEvtPerBinEta, nEvtPerBinPhi, CastorShowerLibraryInfo::Phi, SLenergies, SLetas, SLphis, and totEvents.

Referenced by initFile().

                                                              {
 //
 //  Get EventInfo from the "TBranch* branch" of Root file 
 //  holding library events 
 //
 //  Based on HFShowerLibrary::loadEventInfo
 //
 
   eventInfo = new CastorShowerLibraryInfo();
   branch->SetAddress(&eventInfo);
   branch->GetEntry(0);
   // Initialize shower library general parameters
 
   totEvents   = eventInfo->Energy.getNEvts();
 //  nMomBin     = eventInfo->Energy.getNBins();
 //  evtPerBin   = eventInfo->Energy.getNEvtPerBin();
 //  pmom        = eventInfo->Energy.getBin();
   nBinsE        = eventInfo->Energy.getNBins();
   nEvtPerBinE   = eventInfo->Energy.getNEvtPerBin();
   SLenergies    = eventInfo->Energy.getBin();
   nBinsEta      = eventInfo->Eta.getNBins();
   nEvtPerBinEta = eventInfo->Eta.getNEvtPerBin();
   SLetas        = eventInfo->Eta.getBin();
   nBinsPhi      = eventInfo->Phi.getNBins();
   nEvtPerBinPhi = eventInfo->Phi.getNEvtPerBin();
   SLphis        = eventInfo->Phi.getBin();
   
   // Convert from GeV to MeV
   for (unsigned int i=0; i<SLenergies.size(); i++) SLenergies[i] *= GeV;
   
   edm::LogInfo("CastorShower") << " CastorShowerLibrary::loadEventInfo : " 
                    << "\n \n Total number of events     :  " << totEvents 
                    <<    "\n   Number of bins  (E)       :  " << nBinsE
                    <<    "\n   Number of events/bin (E)  :  " << nEvtPerBinE
                    <<    "\n   Number of bins  (Eta)       :  " << nBinsEta
                    <<    "\n   Number of events/bin (Eta)  :  " << nEvtPerBinEta 
                    <<    "\n   Number of bins  (Phi)       :  " << nBinsPhi
                    <<    "\n   Number of events/bin (Phi)  :  " << nEvtPerBinPhi << "\n";
   
   for (unsigned int i=0; i<nBinsE; i++)
      edm::LogInfo("CastorShower") << "CastorShowerLibrary: SLenergies[" << i << "] = "
                       << SLenergies[i]/GeV << " GeV";
   for (unsigned int i=0; i<nBinsEta; i++)
      edm::LogInfo("CastorShower") << "CastorShowerLibrary: SLetas[" << i << "] = "
                       << SLetas[i];
   for (unsigned int i=0; i<nBinsPhi; i++)
      edm::LogInfo("CastorShower") << "CastorShowerLibrary: SLphis[" << i << "] = "
                       << SLphis[i] << " rad";
 }

void CastorShowerLibrary::select	(	int	type,
		double	pin,
		double	etain = `0`,
		double	phiin = `9`
	)

protected

Definition at line 310 of file CastorShowerLibrary.cc.

References FindEnergyBin(), FindEtaBin(), FindPhiBin(), getRecord(), M_PI, nEvtPerBinE, nEvtPerBinEta, nEvtPerBinPhi, jptDQMConfig_cff::phiMax, jptDQMConfig_cff::phiMin, and csvReporter::r.

Referenced by getShowerHits().

                                                                                  {
 //
 //  Selects an event from the library based on 
 //
 //    type:   0 --> em
 //           >0 --> had
 //    pin  :  momentum
 //    etain:  eta (if not given, disregard the eta binning
 //    phiin:  phi (if not given, disregard the phi binning
 //
 //  Created 30/01/09 by W. Carvalho
 //
 
   int irec;                         // to hold record number
   double r = G4UniformRand();       // to randomly select within an energy bin (r=[0,1])
 
   // Randomly select a record from the right energy bin in the library, 
   // based on track momentum (pin) 
    
   //   pin < SLenergies[MIN]
 /*
   if (pin<SLenergies[0]) {
      edm::LogWarning("CastorShower") << "CastorShowerLibrary: Warning, pin = " << pin 
                                      << " less than minimum SLenergies " << SLenergies[0] << " in library." 
                                      << " For the moment, selecting hit from the first bin" ;
      irec = int(nEvtPerBinE*r) + 1;
   //   pin > SLenergies[MAX]
   } else if (pin>SLenergies[nBinsE]) {
 
 // This part needs rethinking because the last element of SLenergies is no longer the upper limit of the bins
     edm::LogWarning("CastorShower") << "CastorShowerLibrary: Warning, pin = " << pin 
                                     << " greater than maximum SLenergies " << SLenergies[nBinsE] << " in library." 
                                     << " For the moment, selecting hit from the last bin";
     irec = (nBinsE-1)*nEvtPerBinE + int(evtPerBin*r) + 1;
   //   SLenergies[MIN] < pin < SLenergies[MAX]
   } else {
      for (unsigned int j=0; j<nBinsE; j++) {
         if (pin >= SLenergies[j] && pin < SLenergies[j+1]) {
        irec = j*evtPerBin + int(evtPerBin*r) + 1 ;
        if (irec<0) {
           edm::LogWarning("CastorShower") << "CastorShowerLibrary:: Illegal irec = "
                               << irec << " now set to 0";
           irec = 0;
        } else if (irec > totEvents) {
           edm::LogWarning("CastorShower") << "CastorShowerLibrary:: Illegal irec = "
                               << irec << " now set to "<< totEvents;
           irec = totEvents;
        }
         }
      }
   }
 */
   int ienergy = FindEnergyBin(pin);
   int ieta    = FindEtaBin(etain);
 #ifdef DebugLog
   if (verbose) edm::LogInfo("CastorShower") << " ienergy = " << ienergy ;
   if (verbose) edm::LogInfo("CastorShower") << " ieta = " << ieta;
 #endif
 
   int iphi;
   double phiMin = 0. ;
   double phiMax = M_PI/4 ;     // 45 * (pi/180)  rad 
   if(phiin < phiMin) phiin = phiin + M_PI ;
   if(phiin >= phiMin && phiin < phiMax) {
      iphi = FindPhiBin(phiin) ;
   } else {
      double remainder = fmod(phiin , M_PI/4) ;
      phiin = phiMin + remainder ;
      iphi = FindPhiBin(phiin) ;
   }
 #ifdef DebugLog
   if (verbose) edm::LogInfo("CastorShower") << " iphi = " << iphi;
 #endif
   if (ienergy==-1) ienergy=0;    // if pin < first bin, choose an event in the first one
   if (ieta==-1)    ieta=0; // idem for eta
   if (iphi!=-1) irec = int(nEvtPerBinE*ienergy+nEvtPerBinEta*ieta+nEvtPerBinPhi*(iphi+r));
   else irec = int(nEvtPerBinE*(ienergy+r));
 
 #ifdef DebugLog
   edm::LogInfo("CastorShower") << "CastorShowerLibrary:: Select record " << irec 
                                << " of type " << type ; 
 #endif
 
   //  Retrieve record number "irec" from the library
   getRecord (type, irec);
   
 }