#include <CastorShowerLibraryMaker.h>

Inheritance diagram for CastorShowerLibraryMaker:

Classes
struct	ShowerLib

Public Member Functions
	CastorShowerLibraryMaker (const edm::ParameterSet &p)

virtual	~CastorShowerLibraryMaker ()

Public Member Functions inherited from SimWatcher
	SimWatcher ()

virtual	~SimWatcher ()

Public Member Functions inherited from Observer< const BeginOfJob * >
	Observer ()

void	slotForUpdate (const BeginOfJob *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const BeginOfRun * >
	Observer ()

void	slotForUpdate (const BeginOfRun *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfRun * >
	Observer ()

void	slotForUpdate (const EndOfRun *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const BeginOfEvent * >
	Observer ()

void	slotForUpdate (const BeginOfEvent *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const EndOfEvent * >
	Observer ()

void	slotForUpdate (const EndOfEvent *iT)

virtual	~Observer ()

Public Member Functions inherited from Observer< const G4Step * >
	Observer ()

void	slotForUpdate (const G4Step *iT)

virtual	~Observer ()

Private Types
typedef int	ebin

typedef int	etabin

typedef int	phibin

Private Member Functions
bool	FillShowerEvent (CaloG4HitCollection , CastorShowerEvent , int)

int	FindEnergyBin (double e)

int	FindEtaBin (double eta)

int	FindPhiBin (double phi)

void	Finish ()

void	GetKinematics (G4PrimaryParticle *, double &px, double &py, double &pz, double &pInit, double &eta, double &phi)

void	GetKinematics (int, double &px, double &py, double &pz, double &pInit, double &eta, double &phi)

void	GetMissingEnergy (CaloG4HitCollection *, double &, double &)

std::vector< G4PrimaryParticle * >	GetPrimary (const G4Event *)

void	InitSLHolder (ShowerLib &)

bool	IsSLReady ()

void	KillSecondaries (const G4Step *step)

void	printSLstatus (int, int, int)

bool	SLacceptEvent (int, int, int)

bool	SLisEBinFilled (int ebin)

bool	SLisEtaBinFilled (int ebin, int etabin)

bool	SLisPhiBinFilled (int ebin, int etabin, int phibin)

int &	SLnEvtInBinE (int ebin)

int &	SLnEvtInBinEta (int ebin, int etabin)

int &	SLnEvtInBinPhi (int ebin, int etabin, int phibin)

void	update (const BeginOfJob *run)
	This routine will be called when the appropriate signal arrives. More...

void	update (const BeginOfRun *run)
	This routine will be called when the appropriate signal arrives. More...

void	update (const EndOfRun *run)
	This routine will be called when the appropriate signal arrives. More...

void	update (const BeginOfEvent *evt)
	This routine will be called when the appropriate signal arrives. More...

void	update (const EndOfEvent *evt)
	This routine will be called when the appropriate signal arrives. More...

void	update (const G4Step *step)
	This routine will be called when the appropriate signal arrives. More...

Private Attributes
bool	DeActivatePhysicsProcess

bool	DoEmSL

bool	DoHadSL

CastorShowerLibraryInfo *	emInfo

CastorShowerEvent *	emShower

ShowerLib	emSLHolder

int	eventIndex

std::string	eventNtFileName

CastorShowerLibraryInfo *	hadInfo

CastorShowerEvent *	hadShower

ShowerLib	hadSLHolder

bool	InsideCastor

std::map< int, std::set< int > >	MapOfSecondaries

double	MaxEta

double	MaxPhi

unsigned int	NPGParticle

std::vector< int >	PGParticleIDs

std::map< int, G4ThreeVector >	PrimaryMomentum

std::map< int, G4ThreeVector >	PrimaryPosition

ShowerLib *	SLShowerptr

int	stepIndex

TFile *	theFile

std::vector< G4PrimaryParticle * >	thePrims

TTree *	theTree

int	verbosity

Additional Inherited Members

Detailed Description

Definition at line 83 of file CastorShowerLibraryMaker.h.

Member Typedef Documentation

typedef int CastorShowerLibraryMaker::ebin

private

Definition at line 97 of file CastorShowerLibraryMaker.h.

typedef int CastorShowerLibraryMaker::etabin

private

Definition at line 98 of file CastorShowerLibraryMaker.h.

typedef int CastorShowerLibraryMaker::phibin

private

Definition at line 99 of file CastorShowerLibraryMaker.h.

Constructor & Destructor Documentation

CastorShowerLibraryMaker::CastorShowerLibraryMaker ( const edm::ParameterSet & p )

Definition at line 37 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout, DeActivatePhysicsProcess, DoEmSL, DoHadSL, emInfo, emSLHolder, eventNtFileName, edm::ParameterSet::getParameter(), hadInfo, hadSLHolder, i, InitSLHolder(), MapOfSecondaries, MaxEta, MaxPhi, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinE, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinEta, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinPhi, NPGParticle, NULL, PGParticleIDs, CastorShowerLibraryMaker::ShowerLib::SLEnergyBins, CastorShowerLibraryMaker::ShowerLib::SLEtaBins, CastorShowerLibraryMaker::ShowerLib::SLPhiBins, AlCaHLTBitMon_QueryRunRegistry::string, and verbosity.

                                                                            : 
                              NPGParticle(0),DoHadSL(false),DoEmSL(false),DeActivatePhysicsProcess(false),
                              emShower(NULL) , hadShower(NULL) {
 
   MapOfSecondaries.clear();
   hadInfo = NULL;
   emInfo  = NULL;
   edm::ParameterSet p_SLM   = p.getParameter<edm::ParameterSet>("CastorShowerLibraryMaker");
   verbosity                 = p_SLM.getParameter<int>("Verbosity");
   eventNtFileName           = p_SLM.getParameter<std::string>("EventNtupleFileName");
   hadSLHolder.nEvtPerBinPhi = p_SLM.getParameter<int>("nhadEvents");
   emSLHolder.nEvtPerBinPhi  = p_SLM.getParameter<int>("nemEvents");
   hadSLHolder.SLEnergyBins  = p_SLM.getParameter<std::vector<double> >("SLhadEnergyBins");
   hadSLHolder.SLEtaBins     = p_SLM.getParameter<std::vector<double> >("SLhadEtaBins");
   hadSLHolder.SLPhiBins     = p_SLM.getParameter<std::vector<double> >("SLhadPhiBins");
   emSLHolder.SLEnergyBins   = p_SLM.getParameter<std::vector<double> >("SLemEnergyBins");
   emSLHolder.SLEtaBins      = p_SLM.getParameter<std::vector<double> >("SLemEtaBins");
   emSLHolder.SLPhiBins      = p_SLM.getParameter<std::vector<double> >("SLemPhiBins");
   PGParticleIDs             = p_SLM.getParameter<std::vector<int> >("PartID");
   DeActivatePhysicsProcess  = p_SLM.getParameter<bool>("DeActivatePhysicsProcess");
   MaxPhi                    = p_SLM.getParameter<double>("SLMaxPhi");
   MaxEta                    = p_SLM.getParameter<double>("SLMaxEta");
 //
   NPGParticle               = PGParticleIDs.size(); 
   for(unsigned int i=0;i<PGParticleIDs.size();i++) {
      switch (int(fabs(PGParticleIDs.at(i)))) {
         case 11:
         case 22:
                DoEmSL = true;
                break;
         default:
                DoHadSL = true;
      }
   }
   hadSLHolder.nEvtPerBinEta = (hadSLHolder.nEvtPerBinPhi)*(hadSLHolder.SLPhiBins.size());
   hadSLHolder.nEvtPerBinE   = (hadSLHolder.nEvtPerBinEta)*(hadSLHolder.SLEtaBins.size());
   emSLHolder.nEvtPerBinEta = (emSLHolder.nEvtPerBinPhi)*(emSLHolder.SLPhiBins.size());
   emSLHolder.nEvtPerBinE   = (emSLHolder.nEvtPerBinEta)*(emSLHolder.SLEtaBins.size());
 
   std::cout << "============================================================================"<<std::endl;
   std::cout << "CastorShowerLibraryMaker:: Initialized as observer" << std::endl;
   std::cout << " Event Ntuple will be created" << std::endl;
   std::cout << " Event Ntuple file: " << eventNtFileName << std::endl;
   std::cout << " Number of Hadronic events in E   bins: " << hadSLHolder.nEvtPerBinE << std::endl;
   std::cout << " Number of Hadronic events in Eta bins: " << hadSLHolder.nEvtPerBinEta << std::endl;
   std::cout << " Number of Hadronic events in Phi bins: " << hadSLHolder.nEvtPerBinPhi << std::endl;
   std::cout << " Number of Electromag. events in E   bins: " << emSLHolder.nEvtPerBinE << std::endl;
   std::cout << " Number of Electromag. events in Eta bins: " << emSLHolder.nEvtPerBinEta << std::endl;
   std::cout << " Number of Electromag. events in Phi bins: " << emSLHolder.nEvtPerBinPhi << std::endl;
   std::cout << "============================================================================"<<std::endl;
   std::cout << std::endl;
 
   // Initializing the SL collections
   InitSLHolder(hadSLHolder);
   InitSLHolder(emSLHolder);
 }

CastorShowerLibraryMaker::~CastorShowerLibraryMaker ( )

virtual

Definition at line 138 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout, and Finish().

                                                     {
 
   Finish();
 
   std::cout << "CastorShowerLibraryMaker: End of process" << std::endl;
 
 }

Member Function Documentation

bool CastorShowerLibraryMaker::FillShowerEvent	(	CaloG4HitCollection *	theCAFI,
		CastorShowerEvent *	shower,
		int	ipart
	)

private

Definition at line 851 of file CastorShowerLibraryMaker.cc.

References end, spr::find(), CaloG4Hit::getEnergyDeposit(), CaloG4Hit::getEntry(), CaloG4Hit::getIncidentEnergy(), CaloG4Hit::getPosition(), CaloG4Hit::getTimeSlice(), CaloG4Hit::getTrackID(), CaloG4Hit::getUnitID(), MapOfSecondaries, CastorNumberingScheme::packIndex(), position, CastorShowerEvent::setDetID(), CastorShowerEvent::setHitPosition(), CastorShowerEvent::setNhit(), CastorShowerEvent::setNphotons(), CastorShowerEvent::setTime(), cond::rpcobgas::time, CastorNumberingScheme::unpackIndex(), verbosity, and ecaldqm::zside().

Referenced by update().

 {
      unsigned int volumeID=0;
      double en_in_fi = 0.;
      //double totalEnergy = 0;
 
      int nentries = theCAFI->entries();
 
 // Compute Total Energy in CastorFI volume
 /*
      for(int ihit = 0; ihit < nentries; ihit++) {
        CaloG4Hit* aHit = (*theCAFI)[ihit];
        totalEnergy += aHit->getEnergyDeposit();
      }
 */
      if (!shower) {
         edm::LogInfo("CastorShowerLibraryMaker") << "Error. NULL pointer to CastorShowerEvent";
         return false;
      }
 
      CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
 // Hit position
      math::XYZPoint entry;
      math::XYZPoint position;
      int nHits;
      nHits=0;
      for (int ihit = 0; ihit < nentries; ihit++) {
        CaloG4Hit* aHit  = (*theCAFI)[ihit];
        int hit_particleID = aHit->getTrackID();
        if (MapOfSecondaries[ipart].find(hit_particleID)==MapOfSecondaries[ipart].end()) {
           if (verbosity) edm::LogInfo("CastorShowerLibraryMaker")
               << "Skipping hit from trackID " << hit_particleID;
           continue;
        }
        volumeID         = aHit->getUnitID();
        double hitEnergy = aHit->getEnergyDeposit();
        en_in_fi        += aHit->getEnergyDeposit();
        float time       = aHit->getTimeSlice();
        int zside, sector, zmodule;
        theCastorNumScheme->unpackIndex(volumeID, zside, sector,zmodule);
        entry    = aHit->getEntry();
        position = aHit->getPosition();
        if (verbosity) edm::LogInfo("CastorShowerLibraryMaker")
           << "\n side , sector , module = " << zside << " , " 
           << sector << " , " << zmodule 
           << "\n nphotons = " << hitEnergy ; 
 
        if (verbosity) edm::LogInfo("CastorShowerLibraryMaker")
           << "\n packIndex = " 
           << theCastorNumScheme->packIndex(zside, sector,zmodule);
 
        if(verbosity&&time>100.) {
          edm::LogInfo("CastorShowerLibraryMaker")
             << "\n nentries = " << nentries 
             << "\n     time[" << ihit << "] = " << time
             << "\n  trackID[" << ihit << "] = " << aHit->getTrackID()
             << "\n volumeID[" << ihit << "] = " << volumeID 
             << "\n nphotons[" << ihit << "] = " << hitEnergy 
             << "\n side, sector, module  = " << zside <<", " << sector<<", " << zmodule
             << "\n packIndex " << theCastorNumScheme->packIndex(zside,sector,zmodule)
             << "\n X,Y,Z = " << entry.x() << ","<< entry.y() << "," << entry.z();
        }
        if (verbosity) edm::LogInfo("CastorShowerLibraryMaker") << "\n    Incident Energy = "  
                                                 << aHit->getIncidentEnergy() << " \n" ;
 
 //  CaloG4Hit information 
        shower->setDetID(volumeID);
        shower->setHitPosition(position);
        shower->setNphotons(hitEnergy);
        shower->setTime(time);
        nHits++;
      }
 // Write number of hits to CastorShowerEvent instance
      if (nHits==0) {
         edm::LogInfo("CastorShowerLibraryMaker") << "No hits found for this track (trackID=" << ipart << ")." << std::endl;
         if (theCastorNumScheme) delete theCastorNumScheme;
         return false;
      }
      shower->setNhit(nHits);
 
 // update the event counters
      if (theCastorNumScheme) delete theCastorNumScheme;
      return true;
 }

int CastorShowerLibraryMaker::FindEnergyBin ( double e )

private

Definition at line 637 of file CastorShowerLibraryMaker.cc.

References i, CastorShowerLibraryMaker::ShowerLib::SLEnergyBins, and SLShowerptr.

Referenced by update().

                                                          {
   //
   // returns the integer index of the energy bin, taken from SLenergies vector
   // returns -1 if ouside valid range
   //
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nFindEnergyBin can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.\n\n");
   }
   const std::vector<double>& SLenergies = SLShowerptr->SLEnergyBins;
   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 CastorShowerLibraryMaker::FindEtaBin ( double eta )

private

Definition at line 658 of file CastorShowerLibraryMaker.cc.

References i, CastorShowerLibraryMaker::ShowerLib::SLEtaBins, and SLShowerptr.

Referenced by update().

                                                    {
   //
   // returns the integer index of the eta bin, taken from SLetas vector
   // returns -1 if ouside valid range
   //
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nFindEtaBin can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.\n\n");
   }
   const std::vector<double>& SLetas = SLShowerptr->SLEtaBins;
   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 CastorShowerLibraryMaker::FindPhiBin ( double phi )

private

Definition at line 677 of file CastorShowerLibraryMaker.cc.

References i, CastorShowerLibraryMaker::ShowerLib::SLPhiBins, and SLShowerptr.

Referenced by update().

                                                    {
   //
   // 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 (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nFindPhiBin can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.\n\n");
   }
   const std::vector<double>& SLphis = SLShowerptr->SLPhiBins;
   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 CastorShowerLibraryMaker::Finish ( )

private

Definition at line 621 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout, theFile, and theTree.

Referenced by ~CastorShowerLibraryMaker().

                                       {
 
   // if (doNTcastorevent) {
 
   theFile->cd();
   theTree->Write("",TObject::kOverwrite);
   std::cout << "CastorShowerLibraryMaker: Ntuple event written" << std::endl;   
   theFile->Close();
   std::cout << "CastorShowerLibraryMaker: Event file closed" << std::endl;
 
   // Delete pointers to objects, now that TTree has been written and TFile closed
 //  delete      info; 
 //  delete  emShower;
 //  delete hadShower;
   // }
 } 

void CastorShowerLibraryMaker::GetKinematics	(	G4PrimaryParticle *	thePrim,
		double &	px,
		double &	py,
		double &	pz,
		double &	pInit,
		double &	eta,
		double &	phi
	)

private

Definition at line 742 of file CastorShowerLibraryMaker.cc.

References GeV, fff_deleter::log, bookConverter::max, min(), funct::tan(), and theta().

Referenced by update().

 {
     px=py=pz=phi=eta=0.0;
     if (thePrim==0) return;
     px = thePrim->GetMomentum().x()/GeV;
     py = thePrim->GetMomentum().y()/GeV;
     pz = thePrim->GetMomentum().z()/GeV;
     pInit = thePrim->GetMomentum().mag()/GeV;
     //pInit = sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
     if (pInit==0) return;
     double costheta = pz/pInit;
     double theta = acos(std::min(std::max(costheta,double(-1.)),double(1.)));
     eta = -log(tan(theta/2.0));
     phi = (px==0 && py==0) ? 0 : atan2(py,px); // the recommended way of calculating phi
     phi = thePrim->GetMomentum().phi();
     //if (px!=0) phi=atan(py/px);
 }

void CastorShowerLibraryMaker::GetKinematics	(	int	thePrim,
		double &	px,
		double &	py,
		double &	pz,
		double &	pInit,
		double &	eta,
		double &	phi
	)

private

Definition at line 727 of file CastorShowerLibraryMaker.cc.

References GeV, fff_deleter::log, bookConverter::max, min(), PrimaryMomentum, funct::tan(), and theta().

 {
     if (thePrim==0) return;
     if (PrimaryMomentum.find(thePrim)==PrimaryMomentum.end()) return;
     px = PrimaryMomentum[thePrim].x()/GeV;
     py = PrimaryMomentum[thePrim].y()/GeV;
     pz = PrimaryMomentum[thePrim].z()/GeV;
     pInit=PrimaryMomentum[thePrim].mag()/GeV;
     if (pInit==0) return;
     double costheta = pz/pInit;
     double theta = acos(std::min(std::max(costheta,double(-1.)),double(1.)));
     eta = -log(tan(theta/2.0));
     phi = (px==0 && py==0) ? 0 : atan2(py,px); // the recommended way of calculating phi
     phi = PrimaryMomentum[thePrim].phi();
 }

void CastorShowerLibraryMaker::GetMissingEnergy	(	CaloG4HitCollection *	theCAFI,
		double &	miss_energy,
		double &	tot_energy
	)

private

Definition at line 999 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout, end, spr::find(), i, MapOfSecondaries, and thePrims.

Referenced by update().

                                                                                                                   {
 // Get the total deposited energy and the one from hit not associated to a primary
      miss_energy = 0;
      tot_energy = 0;
      for (int ihit = 0; ihit < theCAFI->entries(); ihit++) {
          int id = (*theCAFI)[ihit]->getTrackID();
          tot_energy+=(*theCAFI)[ihit]->getEnergyDeposit();
          int hit_prim = 0;
          for(unsigned int i=0;i<thePrims.size();i++) {
             int particleType = thePrims.at(i)->GetPDGcode();
             if (MapOfSecondaries[particleType].find(id)!=MapOfSecondaries[particleType].end())
                hit_prim = particleType;
          }
          if (hit_prim==0) {
             std::cout << "Track ID " << id << " produced a hit which is not associated with a primary."
             << std::endl;
             miss_energy+=(*theCAFI)[ihit]->getEnergyDeposit();
          }
      }
 }

std::vector< G4PrimaryParticle * > CastorShowerLibraryMaker::GetPrimary ( const G4Event * evt )

private

Definition at line 759 of file CastorShowerLibraryMaker.cc.

References i, j, gen::k, npart, NPGParticle, PGParticleIDs, and thePrims.

Referenced by update().

 {
   // Find Primary info:
   int trackID = 0;
   std::vector<G4PrimaryParticle*> thePrims;
   G4PrimaryParticle* thePrim = 0;
   G4int nvertex = evt->GetNumberOfPrimaryVertex();
   edm::LogInfo("CastorShowerLibraryMaker")  << "Event has " << nvertex << " vertex";   
   if (nvertex!=1) {
      edm::LogInfo("CastorShowerLibraryMaker") << "CastorShowerLibraryMaker::GetPrimary ERROR: no vertex";
      return thePrims;
   }
 
   for (int i = 0 ; i<nvertex; i++) {
     G4PrimaryVertex* avertex = evt->GetPrimaryVertex(i);
     if (avertex == 0) {
        edm::LogInfo("CastorShowerLibraryMaker")
           << "CastorShowerLibraryMaker::GetPrimary ERROR: pointer to vertex = 0";
        continue;
     }
     unsigned int npart = avertex->GetNumberOfParticle();
     if (npart!=NPGParticle) continue;
     for (unsigned int j=0;j<npart;j++) {
         unsigned int k = 0;
         //int test_pID = 0;
         trackID = j;
         thePrim=avertex->GetPrimary(trackID);
         while(k<NPGParticle&&PGParticleIDs.at(k++)!=thePrim->GetPDGcode()){;};
         if (k>NPGParticle) continue; // ID not found in the requested particles
         thePrims.push_back(thePrim);
     }
   }
   return thePrims;
 }

void CastorShowerLibraryMaker::InitSLHolder ( ShowerLib & showerholder )

private

Definition at line 93 of file CastorShowerLibraryMaker.cc.

Referenced by CastorShowerLibraryMaker().

 {
   int nBinsE,nBinsEta,nBinsPhi,nEvtPerBinPhi;
   nBinsE      = showerholder.SLEnergyBins.size();
   nBinsEta    = showerholder.SLEtaBins.size();
   nBinsPhi    = showerholder.SLPhiBins.size();
   nEvtPerBinPhi=showerholder.nEvtPerBinPhi;
 //
 // Info
 //
   showerholder.SLInfo.Energy.setNEvts(nEvtPerBinPhi*nBinsPhi*nBinsEta*nBinsE);
   showerholder.SLInfo.Energy.setNEvtPerBin(nEvtPerBinPhi*nBinsPhi*nBinsEta);
   showerholder.SLInfo.Energy.setNBins(nBinsE);
   showerholder.SLInfo.Energy.setBin(showerholder.SLEnergyBins);
 //
   showerholder.SLInfo.Eta.setNEvts(nEvtPerBinPhi*nBinsPhi*nBinsEta);
   showerholder.SLInfo.Eta.setNEvtPerBin(nEvtPerBinPhi*nBinsPhi);
   showerholder.SLInfo.Eta.setNBins(nBinsEta);
   showerholder.SLInfo.Eta.setBin(showerholder.SLEtaBins);
 //
   showerholder.SLInfo.Phi.setNEvts(nEvtPerBinPhi*nBinsPhi);
   showerholder.SLInfo.Phi.setNEvtPerBin(nEvtPerBinPhi);
   showerholder.SLInfo.Phi.setNBins(nBinsPhi);
   showerholder.SLInfo.Phi.setBin(showerholder.SLPhiBins);
 //
 // Shower
   showerholder.SLCollection.assign(nBinsE,std::vector<std::vector<std::vector<CastorShowerEvent> > >());
   showerholder.nEvtInBinE.assign(nBinsE,0);
   showerholder.nEvtInBinEta.assign(nBinsE,std::vector<int>(0));
   showerholder.nEvtInBinPhi.assign(nBinsE,std::vector<std::vector<int> >());
   for(int i=0;i<nBinsE;i++) {
      showerholder.SLCollection.at(i).assign(nBinsEta,std::vector<std::vector<CastorShowerEvent> >());
      showerholder.nEvtInBinEta.at(i).assign(nBinsEta,0);
      showerholder.nEvtInBinPhi.at(i).assign(nBinsEta,std::vector<int>(0));
      for(int j=0;j<nBinsEta;j++) {
         showerholder.SLCollection.at(i).at(j).assign(nBinsPhi,std::vector<CastorShowerEvent>());
         showerholder.nEvtInBinPhi.at(i).at(j).assign(nBinsPhi,0);
         for(int k=0;k<nBinsPhi;k++) 
            showerholder.SLCollection.at(i).at(j).at(k).assign(nEvtPerBinPhi,CastorShowerEvent());
      }
   }
 }

bool CastorShowerLibraryMaker::IsSLReady ( )

private

Definition at line 698 of file CastorShowerLibraryMaker.cc.

References DoEmSL, DoHadSL, emSLHolder, hadSLHolder, i, NULL, CastorShowerLibraryMaker::ShowerLib::SLEnergyBins, SLisEBinFilled(), and SLShowerptr.

Referenced by update().

 {
 // at this point, the pointer to the shower library should be NULL
   if (SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nIsSLReady must be called when a new event starts.\n\n";
      throw SimG4Exception("\n\nNOT NULL Pointer to the shower library.\n\n");
   }
 // it is enough to check if all the energy bin is filled
   if (DoEmSL) {
      SLShowerptr = &emSLHolder;
      for(unsigned int i=0;i<SLShowerptr->SLEnergyBins.size();i++) {
         if (!SLisEBinFilled(i)) {
            SLShowerptr=NULL;
            return false;
         }
      }
   }
   if (DoHadSL) {
      SLShowerptr = &hadSLHolder;
      for(unsigned int i=0;i<SLShowerptr->SLEnergyBins.size();i++) {
         if (!SLisEBinFilled(i)) {
            SLShowerptr=NULL;
            return false;
         }
      }
   }
   SLShowerptr=NULL;
   return true;
 }

void CastorShowerLibraryMaker::KillSecondaries ( const G4Step * step )

private

Definition at line 988 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout, and i.

Referenced by update().

                                                                    {
    G4TrackVector *p_sec = const_cast<G4TrackVector*>(aStep->GetSecondary());
    for(int i=0;i<int(p_sec->size());i++) {
       /*if (verbosity)*/ std::cout << "Killing track ID " << p_sec->at(i)->GetTrackID() << " and its secondaries"
           << " Produced by Process " << p_sec->at(i)->GetCreatorProcess()->GetProcessName()
           << " in the volume " << aStep->GetTrack()->GetVolume()->GetName()
           << std::endl;
       p_sec->at(i)->SetTrackStatus(fKillTrackAndSecondaries);
    }
 }

void CastorShowerLibraryMaker::printSLstatus	(	int	ebin,
		int	etabin,
		int	phibin
	)

private

Definition at line 793 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout, emSLHolder, hadSLHolder, i, j, gen::k, n, NULL, CastorShowerLibraryMaker::ShowerLib::SLEnergyBins, CastorShowerLibraryMaker::ShowerLib::SLEtaBins, SLisPhiBinFilled(), SLnEvtInBinE(), CastorShowerLibraryMaker::ShowerLib::SLPhiBins, and SLShowerptr.

Referenced by update().

 {
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryInfo") << "NULL shower pointer. Printing both";
      std::cout << "Electromagnetic" << std::endl;
      SLShowerptr = &emSLHolder;
      this->printSLstatus(ebin,etabin,phibin);
      std::cout << "Hadronic" << std::endl;
      SLShowerptr = &hadSLHolder;
      this->printSLstatus(ebin,etabin,phibin);
      SLShowerptr = NULL;
      return;
   }
   int nBinsE  =SLShowerptr->SLEnergyBins.size();
   int nBinsEta=SLShowerptr->SLEtaBins.size();
   int nBinsPhi=SLShowerptr->SLPhiBins.size();
   std::vector<double> SLenergies = SLShowerptr->SLEnergyBins;
   for(int n=0;n<11+(nBinsEta*nBinsPhi);n++) std::cout << "=";
   std::cout << std::endl;
   for(int i=0;i<nBinsE;i++) {
      std::cout << "E bin " << std::setw(6) << SLenergies.at(i) << " : ";
      for(int j=0;j<nBinsEta;j++) {
         for(int k=0;k<nBinsPhi;k++) {
            (SLisPhiBinFilled(i,j,k))?std::cout << "1":std::cout << "-";
         }
         if (j<nBinsEta-1) std::cout << "|";
      }
      std::cout << " (" << SLnEvtInBinE(i) << " events)";
      std::cout << std::endl;
      if (ebin!=i) continue;
      std::cout << "               ";
      for(int j=0;j<nBinsEta;j++) {
         for(int k=0;k<nBinsPhi;k++) {
            (ebin==i&&etabin==j&&phibin==k)?std::cout <<  "^":std::cout << " ";
         }
         if (j<nBinsEta-1) std::cout << " ";
      }
      std::cout << std::endl;
   }
   for(int n=0;n<11+(nBinsEta*nBinsPhi);n++) std::cout << "=";
   std::cout << std::endl;
 }

bool CastorShowerLibraryMaker::SLacceptEvent	(	int	ebin,
		int	etabin,
		int	phibin
	)

private

Definition at line 835 of file CastorShowerLibraryMaker.cc.

References NULL, CastorShowerLibraryMaker::ShowerLib::SLEnergyBins, CastorShowerLibraryMaker::ShowerLib::SLEtaBins, SLisEBinFilled(), SLisEtaBinFilled(), SLisPhiBinFilled(), CastorShowerLibraryMaker::ShowerLib::SLPhiBins, and SLShowerptr.

Referenced by update().

 {
      if (SLShowerptr==NULL) {
         edm::LogInfo("CastorShowerLibraryMaker::SLacceptEvent:") << "Error. NULL pointer to CastorShowerEvent";
         return false;
      }
      if (ebin<0||ebin>=int(SLShowerptr->SLEnergyBins.size())) return false;
      if (SLisEBinFilled(ebin)) return false;
 
      if (etabin<0||etabin>=int(SLShowerptr->SLEtaBins.size())) return false;
      if (SLisEtaBinFilled(ebin,etabin)) return false;
 
      if (phibin<0||phibin>=int(SLShowerptr->SLPhiBins.size())) return false;
      if (SLisPhiBinFilled(ebin,etabin,phibin)) return false;
      return true;
 }

bool CastorShowerLibraryMaker::SLisEBinFilled ( int ebin )

private

Definition at line 961 of file CastorShowerLibraryMaker.cc.

References CastorShowerLibraryMaker::ShowerLib::nEvtInBinE, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinE, and SLShowerptr.

Referenced by IsSLReady(), and SLacceptEvent().

 {
    if (!SLShowerptr) {
       edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLisEBinFilled can be called only after BeginOfEvent\n\n";
       throw SimG4Exception("\n\nNULL Pointer to the shower library.");
    }
    if (SLShowerptr->nEvtInBinE.at(ebin)<(int)SLShowerptr->nEvtPerBinE) return false;
    return true;
 }

bool CastorShowerLibraryMaker::SLisEtaBinFilled	(	int	ebin,
		int	etabin
	)

private

Definition at line 970 of file CastorShowerLibraryMaker.cc.

References CastorShowerLibraryMaker::ShowerLib::nEvtInBinEta, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinEta, and SLShowerptr.

Referenced by SLacceptEvent().

 {
    if (!SLShowerptr) {
       edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLisEtaBinFilled can be called only after BeginOfEvent\n\n";
       throw SimG4Exception("\n\nNULL Pointer to the shower library.");
    }
    if (SLShowerptr->nEvtInBinEta.at(ebin).at(etabin)<(int)SLShowerptr->nEvtPerBinEta) return false;
    return true;
 }

bool CastorShowerLibraryMaker::SLisPhiBinFilled	(	int	ebin,
		int	etabin,
		int	phibin
	)

private

Definition at line 979 of file CastorShowerLibraryMaker.cc.

References CastorShowerLibraryMaker::ShowerLib::nEvtInBinPhi, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinPhi, and SLShowerptr.

Referenced by printSLstatus(), and SLacceptEvent().

 {
    if (!SLShowerptr) {
       edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLisPhiBinFilled can be called only after BeginOfEvent\n\n";
       throw SimG4Exception("\n\nNULL Pointer to the shower library.");
    }
    if (SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin)<(int)SLShowerptr->nEvtPerBinPhi) return false;
    return true;
 }

int & CastorShowerLibraryMaker::SLnEvtInBinE ( int ebin )

private

Definition at line 935 of file CastorShowerLibraryMaker.cc.

References CastorShowerLibraryMaker::ShowerLib::nEvtInBinE, and SLShowerptr.

Referenced by printSLstatus(), and update().

 {
    if (!SLShowerptr) {
       edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinE can be called only after BeginOfEvent\n\n";
       throw SimG4Exception("\n\nNULL Pointer to the shower library.");
    }
    return SLShowerptr->nEvtInBinE.at(ebin);
 }

int & CastorShowerLibraryMaker::SLnEvtInBinEta	(	int	ebin,
		int	etabin
	)

private

Definition at line 944 of file CastorShowerLibraryMaker.cc.

References CastorShowerLibraryMaker::ShowerLib::nEvtInBinEta, and SLShowerptr.

Referenced by update().

 {
    if (!SLShowerptr) {
       edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinEta can be called only after BeginOfEvent\n\n";
       throw SimG4Exception("\n\nNULL Pointer to the shower library.");
    }
    return SLShowerptr->nEvtInBinEta.at(ebin).at(etabin);
 }

int & CastorShowerLibraryMaker::SLnEvtInBinPhi	(	int	ebin,
		int	etabin,
		int	phibin
	)

private

Definition at line 953 of file CastorShowerLibraryMaker.cc.

References CastorShowerLibraryMaker::ShowerLib::nEvtInBinPhi, and SLShowerptr.

Referenced by update().

 {
    if (!SLShowerptr) {
       edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinPhi can be called only after BeginOfEvent\n\n";
       throw SimG4Exception("\n\nNULL Pointer to the shower library.");
    }
    return SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin);
 }

void CastorShowerLibraryMaker::update ( const BeginOfJob * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfJob * >.

Definition at line 147 of file CastorShowerLibraryMaker.cc.

References gather_cfg::cout.

                                                             {
 
   std::cout << " CastorShowerLibraryMaker::Starting new job " << std::endl;
 }

void CastorShowerLibraryMaker::update ( const BeginOfRun * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfRun * >.

Definition at line 153 of file CastorShowerLibraryMaker.cc.

                                                             {
 
   std::cout << std::endl << "CastorShowerLibraryMaker: Starting Run"<< std::endl; 
 
   std::cout << "CastorShowerLibraryMaker: output event root file created" << std::endl;
 
   TString eventfilename = eventNtFileName;
   theFile = new TFile(eventfilename,"RECREATE");
   theTree = new TTree("CastorCherenkovPhotons", "Cherenkov Photons");
 
   Int_t split = 1;
   Int_t bsize = 64000;
   emInfo      = new CastorShowerLibraryInfo();
   emShower    = new CastorShowerEvent();
   hadInfo     = new CastorShowerLibraryInfo();
   hadShower   = new CastorShowerEvent();
   // Create Branchs
   theTree->Branch("emShowerLibInfo.", "CastorShowerLibraryInfo", &emInfo, bsize, split);
   theTree->Branch("emParticles.", "CastorShowerEvent", &emShower, bsize, split);
   theTree->Branch("hadShowerLibInfo.", "CastorShowerLibraryInfo", &hadInfo, bsize, split);
   theTree->Branch("hadParticles.", "CastorShowerEvent", &hadShower, bsize, split);
 
 // set the Info for electromagnetic shower
 // set the energy bins info
   emInfo->Energy.setNEvts(emSLHolder.nEvtPerBinE*emSLHolder.SLEnergyBins.size());
   emInfo->Energy.setNBins(emSLHolder.SLEnergyBins.size());
   emInfo->Energy.setNEvtPerBin(emSLHolder.nEvtPerBinE);
   emInfo->Energy.setBin(emSLHolder.SLEnergyBins);
 // set the eta bins info
   emInfo->Eta.setNEvts(emSLHolder.nEvtPerBinEta*emSLHolder.SLEtaBins.size());
   emInfo->Eta.setNBins(emSLHolder.SLEtaBins.size());
   emInfo->Eta.setNEvtPerBin(emSLHolder.nEvtPerBinEta);
   emInfo->Eta.setBin(emSLHolder.SLEtaBins);
 // set the eta bins info
   emInfo->Phi.setNEvts(emSLHolder.nEvtPerBinPhi*emSLHolder.SLPhiBins.size());
   emInfo->Phi.setNBins(emSLHolder.SLPhiBins.size());
   emInfo->Phi.setNEvtPerBin(emSLHolder.nEvtPerBinPhi);
   emInfo->Phi.setBin(emSLHolder.SLPhiBins);
 // The same for the hadronic shower
 // set the energy bins info
   hadInfo->Energy.setNEvts(hadSLHolder.nEvtPerBinE*hadSLHolder.SLEnergyBins.size());
   hadInfo->Energy.setNBins(hadSLHolder.SLEnergyBins.size());
   hadInfo->Energy.setNEvtPerBin(hadSLHolder.nEvtPerBinE);
   hadInfo->Energy.setBin(hadSLHolder.SLEnergyBins);
 // set the eta bins info
   hadInfo->Eta.setNEvts(hadSLHolder.nEvtPerBinEta*hadSLHolder.SLEtaBins.size());
   hadInfo->Eta.setNBins(hadSLHolder.SLEtaBins.size());
   hadInfo->Eta.setNEvtPerBin(hadSLHolder.nEvtPerBinEta);
   hadInfo->Eta.setBin(hadSLHolder.SLEtaBins);
 // set the eta bins info
   hadInfo->Phi.setNEvts(hadSLHolder.nEvtPerBinPhi*hadSLHolder.SLPhiBins.size());
   hadInfo->Phi.setNBins(hadSLHolder.SLPhiBins.size());
   hadInfo->Phi.setNEvtPerBin(hadSLHolder.nEvtPerBinPhi);
   hadInfo->Phi.setBin(hadSLHolder.SLPhiBins);
   // int flag = theTree->GetBranch("CastorShowerLibInfo")->Fill();
   // Loop on all leaves of this branch to fill Basket buffer.
   // The function returns the number of bytes committed to the memory basket.
   // If a write error occurs, the number of bytes returned is -1.
   // If no data are written, because e.g. the branch is disabled,
   // the number of bytes returned is 0.
   // if(flag==-1) {
   //    edm::LogInfo("CastorAnalyzer") << " WARNING: Error writing to Branch \"CastorShowerLibInfo\" \n" ; 
   // } else 
   // if(flag==0) {
   //    edm::LogInfo("CastorAnalyzer") << " WARNING: No data written to Branch \"CastorShowerLibInfo\" \n" ; 
   // }
 
   // Initialize "accounting" variables
 
   eventIndex = 0;
 
 }

void CastorShowerLibraryMaker::update ( const EndOfRun * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfRun * >.

Definition at line 572 of file CastorShowerLibraryMaker.cc.

References CastorShowerEvent::Clear(), emInfo, emShower, emSLHolder, hadInfo, hadShower, hadSLHolder, IsSLReady(), bookConverter::max, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinE, CastorShowerLibraryMaker::ShowerLib::nEvtPerBinPhi, NULL, CastorShowerLibraryMaker::ShowerLib::SLCollection, CastorShowerLibraryMaker::ShowerLib::SLEnergyBins, CastorShowerLibraryMaker::ShowerLib::SLEtaBins, CastorShowerLibraryMaker::ShowerLib::SLInfo, CastorShowerLibraryMaker::ShowerLib::SLPhiBins, and theTree.

 {
 // Fill the tree with the collected objects
   if (!IsSLReady()) SimG4Exception("\n\nShower Library  NOT READY.\n\n");
 
   unsigned int  ibine,ibineta,ibinphi,ievt; // indexes for em shower
   unsigned int  jbine,jbineta,jbinphi,jevt;// indexes for had shower
 
   ibine=ibineta=ibinphi=ievt=jbine=jbineta=jbinphi=jevt=0;
 
   int  nEvtInTree = 0;
   int  nEMevt  = emSLHolder.nEvtPerBinE*emSLHolder.SLEnergyBins.size();
   int  nHadevt = hadSLHolder.nEvtPerBinE*hadSLHolder.SLEnergyBins.size();
   int  maxEvtInTree=std::max(nEMevt,nHadevt);
 
   emInfo = &emSLHolder.SLInfo;
   hadInfo= &hadSLHolder.SLInfo;
 
   while(nEvtInTree<maxEvtInTree) {
     if (emShower) emShower->Clear();
     if (hadShower) hadShower->Clear();
     while(ibine<emSLHolder.SLEnergyBins.size()&&nEMevt>0){
       emShower = &(emSLHolder.SLCollection.at(ibine).at(ibineta).at(ibinphi).at(ievt));
       ievt++;
       if (ievt==emSLHolder.nEvtPerBinPhi) {ievt=0;ibinphi++;}
       if (ibinphi==emSLHolder.SLPhiBins.size()) {ibinphi=0;ibineta++;}
       if (ibineta==emSLHolder.SLEtaBins.size()) {ibineta=0;ibine++;}
       break;
     }
     while(jbine<hadSLHolder.SLEnergyBins.size()&&nHadevt>0){
       hadShower = &(hadSLHolder.SLCollection.at(jbine).at(jbineta).at(jbinphi).at(jevt));
       jevt++;
       if (jevt==hadSLHolder.nEvtPerBinPhi) {jevt=0;jbinphi++;}
       if (jbinphi==hadSLHolder.SLPhiBins.size()) {jbinphi=0;jbineta++;}
       if (jbineta==hadSLHolder.SLEtaBins.size()) {jbineta=0;jbine++;}
       break;
     }
     theTree->Fill();
     nEvtInTree++;
     if (nEvtInTree==1) {
        theTree->SetBranchStatus("emShowerLibInfo.",0);
        theTree->SetBranchStatus("hadShowerLibInfo.",0);
     }
   }
 // check if run is NULL and exit
   if (run==NULL) throw SimG4Exception("\n\nNumber of needed trigger events reached in CastorShowerLibraryMaker\n\n");
 }

void CastorShowerLibraryMaker::update ( const BeginOfEvent * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 227 of file CastorShowerLibraryMaker.cc.

References accept(), gather_cfg::cout, emSLHolder, eta(), eventIndex, FindEnergyBin(), FindEtaBin(), FindPhiBin(), GetKinematics(), GetPrimary(), hadSLHolder, i, InsideCastor, IsSLReady(), MapOfSecondaries, NULL, phi, PrimaryMomentum, PrimaryPosition, printSLstatus(), SLacceptEvent(), SLShowerptr, stepIndex, AlCaHLTBitMon_QueryRunRegistry::string, thePrims, update(), and verbosity.

                                                               {
 
   eventIndex++;
   stepIndex = 0;
   InsideCastor = false;
   PrimaryMomentum.clear();
   PrimaryPosition.clear();
   int NAccepted =0;
 // reset the pointers to the shower objects
   SLShowerptr = NULL;
   MapOfSecondaries.clear();
   thePrims.clear();
   G4EventManager *e_mgr = G4EventManager::GetEventManager();
   if (IsSLReady()) {
       printSLstatus(-1,-1,-1);
       update((EndOfRun*)NULL);
       return;
   }
 
   thePrims = GetPrimary((*evt)());
   for(unsigned int i=0;i<thePrims.size();i++) {
      G4PrimaryParticle* thePrim = thePrims.at(i);
      int particleType = thePrim->GetPDGcode();
 
      std::string SLType("");
      if (particleType==11) {
         SLShowerptr = &emSLHolder;
         SLType = "Electromagnetic";
      }
      else {
         SLShowerptr = &hadSLHolder;
         SLType = "Hadronic";
      }
      double px=0., py=0., pz=0., pInit = 0., eta = 0., phi = 0.;
      GetKinematics(thePrim,px,py,pz,pInit,eta,phi);
      int ebin = FindEnergyBin(pInit);
      int etabin= FindEtaBin(eta);
      int phibin = FindPhiBin(phi);
      if (verbosity)
         std::cout << "\n========================================================================"
                   << "\nBeginOfEvent: E   : " << pInit <<"\t  bin : " << ebin
                   << "\n              Eta : " << eta   <<"\t  bin : " << etabin
                   << "\n              Phi : " << phi   <<"\t  bin : " << phibin
                   << "\n========================================================================"
                   << std::endl;
 
      if (ebin<0||etabin<0||phibin<0) continue;
      bool accept=false;
      if (!(SLacceptEvent(ebin,etabin,phibin))) {
 /*
 // To increase the chance of a particle arriving at CASTOR inside a not full bin,
 // check if there is available phase space in the neighboring bins
         unsigned int ebin_min = std::max(0,ebin-3);
         unsigned int eta_bin_min = std::max(0,etabin-2);
         unsigned int eta_bin_max = std::min(etabin,etabin+2);
         unsigned int phi_bin_min = std::max(0,phibin-2);
         unsigned int phi_bin_max = std::min(phibin,phibin+2);
         for(unsigned int i_ebin=ebin_min;i_ebin<=(unsigned int)ebin;i_ebin++) {
            for (unsigned int i_etabin=eta_bin_min;i_etabin<=eta_bin_max;i_etabin++) {
                for (unsigned int i_phibin=phi_bin_min;i_phibin<=phi_bin_max;i_phibin++) {
                    if (SLacceptEvent((int)i_ebin,(int)i_etabin,(int)i_phibin)) {accept=true;break;}
                }
                if (accept) break;
            }
            if (accept) break;
         }
 */
         if (!accept) edm::LogInfo("CastorShowerLibraryMaker") << "Event not accepted for ebin="
              << ebin<<",etabin="<<etabin<<",phibin="<<phibin<<std::endl;
      }
      else {
         accept=true;
      }
      if (accept) NAccepted++;
   }
   
   if (NAccepted==0) {
      const_cast<G4Event*>((*evt)())->SetEventAborted();
      const_cast<G4Event*>((*evt)())->KeepTheEvent((G4bool)false);
      e_mgr->AbortCurrentEvent();
   }
   SLShowerptr=NULL;
 //
   std::cout << "CastorShowerLibraryMaker: Processing Event Number: " << eventIndex << std::endl;
 }

void CastorShowerLibraryMaker::update ( const EndOfEvent * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 418 of file CastorShowerLibraryMaker.cc.

                                                             {
 
 // check if the job is done!
   if ((*evt)()->IsAborted()) {
      std::cout << "\n========================================================================"
                << "\nEndOfEvent: EVENT ABORTED"
                << "\n========================================================================"
                << std::endl;
      return;
   }
   //DynamicRangeFlatRandomEGunProducer* pgun = edm::DynamicRangeFlatRandomEGunKernel::get_instance();
   //std::cout << pgun->EGunMaxE() << std::endl;
 /*
   std::cout << "Minimum energy in Particle Gun : " << pgun->EGunMinE() << "\n"
             << "Maximum energy in Particle Gun : " << pgun->EGunMaxE() << std::endl;
 */
   if (verbosity) 
      std::cout << "CastorShowerLibraryMaker: End of Event: " << eventIndex << std::endl;
 // Get the pointer to the primary particle
   if (thePrims.size() == 0) {
      edm::LogInfo("CastorShowerLibraryMaker") << "No valid primary particle found. Skipping event" << std::endl;
      return;
   }
   // access to the G4 hit collections 
   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
   int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
   CaloG4HitCollection* theCAFI = (CaloG4HitCollection*) allHC->GetHC(CAFIid);
   if (verbosity) edm::LogInfo("CastorShowerLibraryMaker") << " update(*evt) --> accessed all HC ";
   edm::LogInfo("CastorShowerLibraryMaker") << "Found "<<theCAFI->entries() << " hits in G4HitCollection";
   if (theCAFI->entries()== 0) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n Empty G4HitCollection";
      return;
   }
 
 // Loop over primaries
   int NEvtAccepted = 0;
   int NHitInEvent  = 0;
   for(unsigned int i=0;i<thePrims.size();i++) {
      G4PrimaryParticle* thePrim = thePrims.at(i);
      if (!thePrim) {
         edm::LogInfo("CastorShowerLibraryMaker") << "NULL Pointer to the primary" << std::endl;
         continue;
      }
 // Check primary particle type
      int particleType = thePrim->GetPDGcode();
 
 // set the pointer to the shower collection
      std::string SLType("");
      if (particleType==11) {
         SLShowerptr = &emSLHolder;
         SLType = "Electromagnetic";
      }
      else {
         SLShowerptr = &hadSLHolder;
         SLType = "Hadronic";
      }
      edm::LogInfo("CastorShowerLibraryMaker") << "\n Primary (thePrim) trackID is " << thePrim->GetTrackID() << "\n" ;
 
 // Obtain primary particle's initial momentum (pInit)
      double px=0., py=0., pz=0., pInit = 0., eta = 0., phi = 0.;
      GetKinematics(particleType,px,py,pz,pInit,eta,phi);
 // Check if current event falls into any bin
 // first: energy
      if (pInit==0) {
         edm::LogInfo("CastorShowerLibraryMaker") << "Primary did not hit CASTOR" << std::endl;
         continue;
      }
      int ebin = FindEnergyBin(pInit);
      int etabin= FindEtaBin(eta);
      int phibin = FindPhiBin(phi);
      std::cout << SLType << std::endl;
      printSLstatus(ebin,etabin,phibin);
      if (!SLacceptEvent(ebin,etabin,phibin)) {
         edm::LogInfo("CastorShowerLibraryMaker") << "Event not accepted for ebin="
              << ebin<<",etabin="<<etabin<<",phibin="<<phibin
              <<"("<< pInit <<","<<eta<<","<<phi<<")" <<std::endl;
         continue;
      }
 //
 // event passed. Fill the vector accordingly
 //  
 // Look for the Hit Collection 
      edm::LogInfo("CastorShowerLibraryMaker")
         << "\n CastorShowerLibraryMaker::update(EndOfEvent * evt) - event #" 
         << (*evt)()->GetEventID() ;
 
 /*
      std::cout << "Number of collections : " << allHC->GetNumberOfCollections() << std::endl;
      for(int ii = 0;ii<allHC->GetNumberOfCollections();ii++) 
         std::cout << "Name of collection " << ii << " : " << allHC->GetHC(ii)->GetName() << std::endl;
 */
 
      CastorShowerEvent* shower=NULL;
      int cur_evt_idx = SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin);
      shower = &(SLShowerptr->SLCollection.at(ebin).at(etabin).at(phibin).at(cur_evt_idx));
 
 // Get Hit information
      if (FillShowerEvent(theCAFI,shower,particleType)) { 
 //  Primary particle information
 /*
         edm::LogInfo("CastorShowerLibraryMaker") << "New SL event: Primary = " << particleType
              << "; Energy = " << pInit << "; Eta = " << eta << "; Phi = " << phi
              << "; Nhits = " << shower->getNhit() << std::endl;
 */
         shower->setPrimE(pInit);
         shower->setPrimEta(eta);
         shower->setPrimPhi(phi);
         shower->setPrimX(PrimaryPosition[particleType].x());
         shower->setPrimY(PrimaryPosition[particleType].y());
         shower->setPrimZ(PrimaryPosition[particleType].z());
         SLnEvtInBinE(ebin)++;
         SLnEvtInBinEta(ebin,etabin)++;
         SLnEvtInBinPhi(ebin,etabin,phibin)++;
         NHitInEvent+=shower->getNhit();
         NEvtAccepted++;
       }
       else {shower->Clear();}
   }
 // Check for unassociated energy
      
   if (NEvtAccepted==int(thePrims.size())&&theCAFI->entries()!=NHitInEvent){
      std::cout << "WARNING: Inconsistent Number of Hits -> Hits in collection: "<< theCAFI->entries()
           << "   Hits in the showers: " << NHitInEvent << std::endl;
      double miss_energy=0;
      double tot_energy=0;
      GetMissingEnergy(theCAFI,miss_energy,tot_energy);
      if (miss_energy>0) {
         std::cout << "Total missing energy: " << miss_energy
                << " for an incident energy: " << tot_energy
                << std::endl;
      }
   }
          
 /*
   for (int i=emSLHolder.SLEnergyBins.size()-1;i>0;i--) {
       if (emSLHolder.nEvtInBinE.at(i)==(int)emSLHolder.nEvtPerBinE) {
          std::ostringstream out;
          out << emSLHolder.SLEnergyBins.at(i);
          cout << "Bin Limit: " << out.str() << std::endl;
          setenv("CASTOR_SL_PG_MAXE",out.str().c_str(),1);
        }
        break;
    }
 */
   //int iEvt = (*evt)()->GetEventID();
   //double xint;
 /*
   if (modf(log10(iEvt),&xint)==0) 
     std::cout << " CastorShowerLibraryMaker Event " << iEvt << std::endl;
 */
   // std::cout << std::endl << "===>>> Done writing user histograms " << std::endl;
 }

void CastorShowerLibraryMaker::update ( const G4Step * )

privatevirtual

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 314 of file CastorShowerLibraryMaker.cc.

References funct::abs(), funct::cos(), gather_cfg::cout, DeActivatePhysicsProcess, i, InsideCastor, KillSecondaries(), MapOfSecondaries, MaxEta, MaxPhi, PrimaryMomentum, PrimaryPosition, proc, alignCSCRings::r, funct::sin(), AlCaHLTBitMon_QueryRunRegistry::string, edmStreamStallGrapher::t, and verbosity.

                                                           {
    static int CurrentPrimary = 0;
    G4Track *trk = aStep->GetTrack();
    if (trk->GetCurrentStepNumber()==1) {
       if (trk->GetParentID()==0) {
          CurrentPrimary = (int)trk->GetDynamicParticle()->GetPDGcode();
          if (CurrentPrimary==0) 
              SimG4Exception("CastorShowerLibraryMaker::update(G4Step) -> Primary particle undefined");
          InsideCastor=false;
 // Deactivate the physics process
          if (DeActivatePhysicsProcess) {
             G4ProcessManager *p_mgr = trk->GetDefinition()->GetProcessManager();
             G4ProcessVector *pvec = p_mgr->GetProcessList();
             for (int i = 0;i<pvec->size();i++) {
                 G4VProcess *proc = (*pvec)(i);
                 if (proc->GetProcessName()!="Transportation"&&
                     proc->GetProcessName()!="Decay") {
                     std::cout << "DeActivating process: " << proc->GetProcessName() << std::endl;
                     p_mgr->SetProcessActivation(proc,false);
                 }
             } 
          }
 // move track to z of CASTOR
          G4ThreeVector pos;
          pos.setZ(-14390);
          double t = abs((pos.z()-trk->GetPosition().z()))/trk->GetVelocity();
          double r = (pos.z()-trk->GetPosition().z())/trk->GetMomentum().cosTheta();
          pos.setX(r*sin(trk->GetMomentum().theta())*cos(trk->GetMomentum().phi())+trk->GetPosition().x());
          pos.setY(r*sin(trk->GetMomentum().theta())*sin(trk->GetMomentum().phi())+trk->GetPosition().y());
          trk->SetPosition(pos);
          trk->SetGlobalTime(trk->GetGlobalTime()+t);
          trk->AddTrackLength(r);
       }
       else if (!InsideCastor) {
          std::cout<<"CastorShowerLibraryMaker::update(G4Step) -> Killing spurious track" << std::endl;
          trk->SetTrackStatus(fKillTrackAndSecondaries);
          return;
       }
       MapOfSecondaries[CurrentPrimary].insert((int)trk->GetTrackID());
    }
 // Checks if primary already inside CASTOR
    std::string CurVolume = trk->GetVolume()->GetName();
    if (!InsideCastor&&(
        //CurVolume=="C3EF"||CurVolume=="C4EF"||CurVolume=="CAEL"||
        //CurVolume=="CAHL"||CurVolume=="C3HF"||CurVolume=="C4HF")) {
        //CurVolume=="CastorB"||
        CurVolume=="CAST")) {
        //CurVolume=="CAIR")) {
        InsideCastor=true;
 // Activate the physics process
        if (trk->GetParentID()==0&&DeActivatePhysicsProcess) {
           G4ProcessManager *p_mgr = trk->GetDefinition()->GetProcessManager();
           G4ProcessVector *pvec = p_mgr->GetProcessList();
           for (int i = 0;i<pvec->size();i++) {
               G4VProcess *proc = (*pvec)(i);
               if (proc->GetProcessName()!="Transportation"&&
                  proc->GetProcessName()!="Decay") {
                  std::cout << "  Activating process: " << proc->GetProcessName() << std::endl;
                  p_mgr->SetProcessActivation(proc,true);
               }
           }
        }
        //PrimaryMomentum[CurrentPrimary]=aStep->GetPreStepPoint()->GetMomentum();
 // check fiducial eta and phi
        if (trk->GetMomentum().phi()>MaxPhi||trk->GetMomentum().eta()>MaxEta) {
           trk->SetTrackStatus(fKillTrackAndSecondaries);
           InsideCastor=false;
           return;
        }
        PrimaryMomentum[CurrentPrimary]=trk->GetMomentum();
        PrimaryPosition[CurrentPrimary]=trk->GetPosition();
        KillSecondaries(aStep);
        return;
    } 
 // Kill the secondaries if they have been produced before entering castor
    if (CurrentPrimary!=0&&trk->GetParentID()==0&&!InsideCastor) {
       KillSecondaries(aStep);
       if (verbosity) {
          double pre_phi = aStep->GetPreStepPoint()->GetMomentum().phi();
          double cur_phi = trk->GetMomentum().phi();
          if (pre_phi!=cur_phi) {
             std::cout << "Primary track phi :  " << pre_phi << " changed in current step: "
                       << cur_phi << " by processes:" << std::endl;
             const G4VProcess *proc = aStep->GetPreStepPoint()->GetProcessDefinedStep();
             std::cout << "           " << proc->GetProcessName()
                       << "  In volume " << CurVolume << std::endl;
          }
       }
    }
 
 //==============================================
 /*
 */
 /*
   if(aStep->IsFirstStepInVolume()) { 
     edm::LogInfo("CastorShowerLibraryMaker") << "CastorShowerLibraryMaker::update(const G4Step * aStep):"
                                              << "\n IsFirstStepInVolume , " 
                                              << "time = " << aStep->GetTrack()->GetGlobalTime() ; 
   }
   stepIndex++;
 */
 }