35 NPGParticle(0),DoHadSL(
false),DoEmSL(
false),DeActivatePhysicsProcess(
false),
58 for(
unsigned int i=0;
i<PGParticleIDs.size();
i++) {
59 switch (
int(fabs(PGParticleIDs.at(
i)))) {
73 std::cout <<
"============================================================================"<<std::endl;
74 std::cout <<
"CastorShowerLibraryMaker:: Initialized as observer" << std::endl;
75 std::cout <<
" Event Ntuple will be created" << std::endl;
76 std::cout <<
" Event Ntuple file: " << eventNtFileName << std::endl;
83 std::cout <<
"============================================================================"<<std::endl;
92 int nBinsE,nBinsEta,nBinsPhi,nEvtPerBinPhi;
116 showerholder.
SLCollection.assign(nBinsE,std::vector<std::vector<std::vector<CastorShowerEvent> > >());
118 showerholder.
nEvtInBinEta.assign(nBinsE,std::vector<int>(0));
119 showerholder.
nEvtInBinPhi.assign(nBinsE,std::vector<std::vector<int> >());
120 for(
int i=0;
i<nBinsE;
i++) {
121 showerholder.
SLCollection.at(
i).assign(nBinsEta,std::vector<std::vector<CastorShowerEvent> >());
123 showerholder.
nEvtInBinPhi.at(
i).assign(nBinsEta,std::vector<int>(0));
124 for(
int j=0;
j<nBinsEta;
j++) {
125 showerholder.
SLCollection.at(
i).at(
j).assign(nBinsPhi,std::vector<CastorShowerEvent>());
127 for(
int k=0;
k<nBinsPhi;
k++)
139 std::cout <<
"CastorShowerLibraryMaker: End of process" << std::endl;
146 std::cout <<
" CastorShowerLibraryMaker::Starting new job " << std::endl;
152 std::cout << std::endl <<
"CastorShowerLibraryMaker: Starting Run"<< std::endl;
154 std::cout <<
"CastorShowerLibraryMaker: output event root file created" << std::endl;
157 theFile =
new TFile(eventfilename,
"RECREATE");
158 theTree =
new TTree(
"CastorCherenkovPhotons",
"Cherenkov Photons");
167 theTree->Branch(
"emShowerLibInfo.",
"CastorShowerLibraryInfo", &
emInfo, bsize, split);
168 theTree->Branch(
"emParticles.",
"CastorShowerEvent", &
emShower, bsize, split);
169 theTree->Branch(
"hadShowerLibInfo.",
"CastorShowerLibraryInfo", &
hadInfo, bsize, split);
170 theTree->Branch(
"hadParticles.",
"CastorShowerEvent", &
hadShower, bsize, split);
236 G4EventManager *e_mgr = G4EventManager::GetEventManager();
245 G4PrimaryParticle* thePrim =
thePrims.at(
i);
246 int particleType = thePrim->GetPDGcode();
249 if (particleType==11) {
251 SLType =
"Electromagnetic";
257 double px=0., py=0., pz=0., pInit = 0.,
eta = 0.,
phi = 0.;
263 std::cout <<
"\n========================================================================"
264 <<
"\nBeginOfEvent: E : " << pInit <<
"\t bin : " << ebin
265 <<
"\n Eta : " <<
eta <<
"\t bin : " << etabin
266 <<
"\n Phi : " <<
phi <<
"\t bin : " << phibin
267 <<
"\n========================================================================"
270 if (ebin<0||etabin<0||phibin<0)
continue;
291 if (!accept)
edm::LogInfo(
"CastorShowerLibraryMaker") <<
"Event not accepted for ebin="
292 << ebin<<
",etabin="<<etabin<<
",phibin="<<phibin<<std::endl;
297 if (accept) NAccepted++;
301 const_cast<G4Event*
>((*evt)())->SetEventAborted();
302 const_cast<G4Event*
>((*evt)())->KeepTheEvent((G4bool)
false);
303 e_mgr->AbortCurrentEvent();
307 std::cout <<
"CastorShowerLibraryMaker: Processing Event Number: " <<
eventIndex << std::endl;
312 static int CurrentPrimary = 0;
313 G4Track *trk = aStep->GetTrack();
314 if (trk->GetCurrentStepNumber()==1) {
315 if (trk->GetParentID()==0) {
316 CurrentPrimary = (int)trk->GetDynamicParticle()->GetPDGcode();
317 if (CurrentPrimary==0)
318 SimG4Exception(
"CastorShowerLibraryMaker::update(G4Step) -> Primary particle undefined");
322 G4ProcessManager *p_mgr = trk->GetDefinition()->GetProcessManager();
323 G4ProcessVector *pvec = p_mgr->GetProcessList();
324 for (
int i = 0;
i<pvec->size();
i++) {
325 G4VProcess *
proc = (*pvec)(
i);
326 if (proc->GetProcessName()!=
"Transportation"&&
327 proc->GetProcessName()!=
"Decay") {
328 std::cout <<
"DeActivating process: " << proc->GetProcessName() << std::endl;
329 p_mgr->SetProcessActivation(proc,
false);
336 double t =
abs((pos.z()-trk->GetPosition().z()))/trk->GetVelocity();
337 double r = (pos.z()-trk->GetPosition().z())/trk->GetMomentum().cosTheta();
338 pos.setX(r*
sin(trk->GetMomentum().theta())*
cos(trk->GetMomentum().phi())+trk->GetPosition().x());
339 pos.setY(r*
sin(trk->GetMomentum().theta())*
sin(trk->GetMomentum().phi())+trk->GetPosition().y());
340 trk->SetPosition(pos);
341 trk->SetGlobalTime(trk->GetGlobalTime()+
t);
342 trk->AddTrackLength(r);
345 std::cout<<
"CastorShowerLibraryMaker::update(G4Step) -> Killing spurious track" << std::endl;
346 trk->SetTrackStatus(fKillTrackAndSecondaries);
352 std::string CurVolume = trk->GetVolume()->GetName();
357 CurVolume==
"CAST")) {
362 G4ProcessManager *p_mgr = trk->GetDefinition()->GetProcessManager();
363 G4ProcessVector *pvec = p_mgr->GetProcessList();
364 for (
int i = 0;
i<pvec->size();
i++) {
365 G4VProcess *
proc = (*pvec)(
i);
366 if (proc->GetProcessName()!=
"Transportation"&&
367 proc->GetProcessName()!=
"Decay") {
368 std::cout <<
" Activating process: " << proc->GetProcessName() << std::endl;
369 p_mgr->SetProcessActivation(proc,
true);
375 if (trk->GetMomentum().phi()>
MaxPhi||trk->GetMomentum().eta()>
MaxEta) {
376 trk->SetTrackStatus(fKillTrackAndSecondaries);
386 if (CurrentPrimary!=0&&trk->GetParentID()==0&&!
InsideCastor) {
389 double pre_phi = aStep->GetPreStepPoint()->GetMomentum().phi();
390 double cur_phi = trk->GetMomentum().phi();
391 if (pre_phi!=cur_phi) {
392 std::cout <<
"Primary track phi : " << pre_phi <<
" changed in current step: "
393 << cur_phi <<
" by processes:" << std::endl;
394 const G4VProcess *
proc = aStep->GetPreStepPoint()->GetProcessDefinedStep();
395 std::cout <<
" " << proc->GetProcessName()
396 <<
" In volume " << CurVolume << std::endl;
418 if ((*evt)()->IsAborted()) {
419 std::cout <<
"\n========================================================================"
420 <<
"\nEndOfEvent: EVENT ABORTED"
421 <<
"\n========================================================================"
435 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"No valid primary particle found. Skipping event" << std::endl;
439 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
440 int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID(
"CastorFI");
443 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"Found "<<theCAFI->entries() <<
" hits in G4HitCollection";
444 if (theCAFI->entries()== 0) {
445 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n Empty G4HitCollection";
450 int NEvtAccepted = 0;
453 G4PrimaryParticle* thePrim =
thePrims.at(
i);
455 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"NULL Pointer to the primary" << std::endl;
459 int particleType = thePrim->GetPDGcode();
463 if (particleType==11) {
465 SLType =
"Electromagnetic";
471 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n Primary (thePrim) trackID is " << thePrim->GetTrackID() <<
"\n" ;
474 double px=0., py=0., pz=0., pInit = 0.,
eta = 0.,
phi = 0.;
479 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"Primary did not hit CASTOR" << std::endl;
488 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"Event not accepted for ebin="
489 << ebin<<
",etabin="<<etabin<<
",phibin="<<phibin
490 <<
"("<< pInit <<
","<<
eta<<
","<<
phi<<
")" <<std::endl;
498 <<
"\n CastorShowerLibraryMaker::update(EndOfEvent * evt) - event #"
499 << (*evt)()->GetEventID() ;
528 NHitInEvent+=shower->
getNhit();
531 else {shower->
Clear();}
535 if (NEvtAccepted==
int(
thePrims.size())&&theCAFI->entries()!=NHitInEvent){
536 std::cout <<
"WARNING: Inconsistent Number of Hits -> Hits in collection: "<< theCAFI->entries()
537 <<
" Hits in the showers: " << NHitInEvent << std::endl;
538 double miss_energy=0;
542 std::cout <<
"Total missing energy: " << miss_energy
543 <<
" for an incident energy: " << tot_energy
574 unsigned int ibine,ibineta,ibinphi,ievt;
575 unsigned int jbine,jbineta,jbinphi,jevt;
577 ibine=ibineta=ibinphi=ievt=jbine=jbineta=jbinphi=jevt=0;
582 int maxEvtInTree=
std::max(nEMevt,nHadevt);
587 while(nEvtInTree<maxEvtInTree) {
609 theTree->SetBranchStatus(
"emShowerLibInfo.",0);
610 theTree->SetBranchStatus(
"hadShowerLibInfo.",0);
614 if (run==
NULL)
throw SimG4Exception(
"\n\nNumber of needed trigger events reached in CastorShowerLibraryMaker\n\n");
623 theTree->Write(
"",TObject::kOverwrite);
624 std::cout <<
"CastorShowerLibraryMaker: Ntuple event written" << std::endl;
626 std::cout <<
"CastorShowerLibraryMaker: Event file closed" << std::endl;
640 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nFindEnergyBin can be called only after BeginOfEvent\n\n";
641 throw SimG4Exception(
"\n\nNULL Pointer to the shower library.\n\n");
644 if (energy >= SLenergies.back())
return SLenergies.size()-1;
647 for(;i<SLenergies.size()-1;i++)
648 if (energy >= SLenergies.at(i) && energy < SLenergies.at(i+1))
return (
int)
i;
651 if (energy>=SLenergies.at(i))
return (
int)
i;
661 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nFindEtaBin can be called only after BeginOfEvent\n\n";
662 throw SimG4Exception(
"\n\nNULL Pointer to the shower library.\n\n");
665 if (eta>=SLetas.back())
return SLetas.size()-1;
667 for(;i<SLetas.size()-1;i++)
668 if (eta >= SLetas.at(i) && eta < SLetas.at(i+1))
return (
int)
i;
670 if (eta>=SLetas.at(i))
return (
int)
i;
682 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nFindPhiBin can be called only after BeginOfEvent\n\n";
683 throw SimG4Exception(
"\n\nNULL Pointer to the shower library.\n\n");
686 if (phi>=SLphis.back())
return SLphis.size()-1;
688 for(;i<SLphis.size()-1;i++)
689 if (phi >= SLphis.at(i) && phi < SLphis.at(i+1))
return (
int)
i;
691 if (phi>=SLphis.at(i))
return (
int)
i;
699 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nIsSLReady must be called when a new event starts.\n\n";
700 throw SimG4Exception(
"\n\nNOT NULL Pointer to the shower library.\n\n");
726 if (thePrim==0)
return;
732 if (pInit==0)
return;
733 double costheta = pz/pInit;
735 eta = -
log(
tan(theta/2.0));
736 phi = (px==0 && py==0) ? 0 : atan2(py,px);
741 px=py=pz=phi=eta=0.0;
742 if (thePrim==0)
return;
743 px = thePrim->GetMomentum().x()/GeV;
744 py = thePrim->GetMomentum().y()/GeV;
745 pz = thePrim->GetMomentum().z()/GeV;
746 pInit = thePrim->GetMomentum().mag()/GeV;
748 if (pInit==0)
return;
749 double costheta = pz/pInit;
751 eta = -
log(
tan(theta/2.0));
752 phi = (px==0 && py==0) ? 0 : atan2(py,px);
753 phi = thePrim->GetMomentum().phi();
760 std::vector<G4PrimaryParticle*>
thePrims;
761 G4PrimaryParticle* thePrim = 0;
762 G4int nvertex = evt->GetNumberOfPrimaryVertex();
763 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"Event has " << nvertex <<
" vertex";
765 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"CastorShowerLibraryMaker::GetPrimary ERROR: no vertex";
769 for (
int i = 0 ;
i<nvertex;
i++) {
770 G4PrimaryVertex* avertex = evt->GetPrimaryVertex(
i);
773 <<
"CastorShowerLibraryMaker::GetPrimary ERROR: pointer to vertex = 0";
776 unsigned int npart = avertex->GetNumberOfParticle();
778 for (
unsigned int j=0;
j<
npart;
j++) {
782 thePrim=avertex->GetPrimary(trackID);
785 thePrims.push_back(thePrim);
793 edm::LogInfo(
"CastorShowerLibraryInfo") <<
"NULL shower pointer. Printing both";
794 std::cout <<
"Electromagnetic" << std::endl;
807 for(
int n=0;
n<11+(nBinsEta*nBinsPhi);
n++)
std::cout <<
"=";
809 for(
int i=0;
i<nBinsE;
i++) {
810 std::cout <<
"E bin " << std::setw(6) << SLenergies.at(
i) <<
" : ";
811 for(
int j=0;
j<nBinsEta;
j++) {
812 for(
int k=0;
k<nBinsPhi;
k++) {
819 if (ebin!=
i)
continue;
821 for(
int j=0;
j<nBinsEta;
j++) {
822 for(
int k=0;
k<nBinsPhi;
k++) {
829 for(
int n=0;
n<11+(nBinsEta*nBinsPhi);
n++)
std::cout <<
"=";
835 edm::LogInfo(
"CastorShowerLibraryMaker::SLacceptEvent:") <<
"Error. NULL pointer to CastorShowerEvent";
850 unsigned int volumeID=0;
851 double en_in_fi = 0.;
854 int nentries = theCAFI->entries();
864 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"Error. NULL pointer to CastorShowerEvent";
874 for (
int ihit = 0; ihit < nentries; ihit++) {
879 <<
"Skipping hit from trackID " << hit_particleID;
886 int zside, sector, zmodule;
887 theCastorNumScheme->
unpackIndex(volumeID, zside, sector,zmodule);
891 <<
"\n side , sector , module = " << zside <<
" , "
892 << sector <<
" , " << zmodule
893 <<
"\n nphotons = " << hitEnergy ;
897 << theCastorNumScheme->
packIndex(zside, sector,zmodule);
901 <<
"\n nentries = " << nentries
902 <<
"\n time[" << ihit <<
"] = " << time
903 <<
"\n trackID[" << ihit <<
"] = " << aHit->
getTrackID()
904 <<
"\n volumeID[" << ihit <<
"] = " << volumeID
905 <<
"\n nphotons[" << ihit <<
"] = " << hitEnergy
906 <<
"\n side, sector, module = " << zside <<
", " << sector<<
", " << zmodule
907 <<
"\n packIndex " << theCastorNumScheme->
packIndex(zside,sector,zmodule)
908 <<
"\n X,Y,Z = " << entry.x() <<
","<< entry.y() <<
"," << entry.z();
922 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"No hits found for this track (trackID=" << ipart <<
")." << std::endl;
923 if (theCastorNumScheme)
delete theCastorNumScheme;
929 if (theCastorNumScheme)
delete theCastorNumScheme;
935 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nSLnEvtInBinE can be called only after BeginOfEvent\n\n";
944 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nSLnEvtInBinEta can be called only after BeginOfEvent\n\n";
953 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nSLnEvtInBinPhi can be called only after BeginOfEvent\n\n";
961 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nSLisEBinFilled can be called only after BeginOfEvent\n\n";
970 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nSLisEtaBinFilled can be called only after BeginOfEvent\n\n";
979 edm::LogInfo(
"CastorShowerLibraryMaker") <<
"\n\nSLisPhiBinFilled can be called only after BeginOfEvent\n\n";
986 G4TrackVector *p_sec =
const_cast<G4TrackVector*
>(aStep->GetSecondary());
987 for(
int i=0;
i<int(p_sec->size());
i++) {
988 std::cout <<
"Killing track ID " << p_sec->at(
i)->GetTrackID() <<
" and its secondaries"
989 <<
" Produced by Process " << p_sec->at(
i)->GetCreatorProcess()->GetProcessName()
990 <<
" in the volume " << aStep->GetTrack()->GetVolume()->GetName()
992 p_sec->at(
i)->SetTrackStatus(fKillTrackAndSecondaries);
1000 for (
int ihit = 0; ihit < theCAFI->entries(); ihit++) {
1001 int id = (*theCAFI)[ihit]->getTrackID();
1002 tot_energy+=(*theCAFI)[ihit]->getEnergyDeposit();
1005 int particleType =
thePrims.at(
i)->GetPDGcode();
1007 hit_prim = particleType;
1010 std::cout <<
"Track ID " <<
id <<
" produced a hit which is not associated with a primary."
1012 miss_energy+=(*theCAFI)[ihit]->getEnergyDeposit();
std::string eventNtFileName
std::vector< double > SLEtaBins
T getParameter(std::string const &) const
void setNBins(unsigned int n)
math::XYZPoint getPosition() const
bool DeActivatePhysicsProcess
void setPrimEta(float eta)
int & SLnEvtInBinE(int ebin)
void setDetID(unsigned int id)
CastorShowerEvent * hadShower
std::vector< double > SLPhiBins
int & SLnEvtInBinEta(int ebin, int etabin)
TrainProcessor *const proc
void setPrimPhi(float phi)
std::vector< G4PrimaryParticle * > thePrims
CastorShowerLibraryInfo SLInfo
Sin< T >::type sin(const T &t)
virtual ~CastorShowerLibraryMaker()
Geom::Theta< T > theta() const
bool FillShowerEvent(CaloG4HitCollection *, CastorShowerEvent *, int)
double getIncidentEnergy() const
void KillSecondaries(const G4Step *step)
void InitSLHolder(ShowerLib &)
void printSLstatus(int, int, int)
bool SLacceptEvent(int, int, int)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
bool accept(const edm::Event &event, const edm::TriggerResults &triggerTable, const std::string &triggerPath)
bool SLisEtaBinFilled(int ebin, int etabin)
static void unpackIndex(const uint32_t &idx, int &z, int §or, int &zmodule)
static int position[TOTALCHAMBERS][3]
void setNEvtPerBin(unsigned int n)
std::vector< G4PrimaryParticle * > GetPrimary(const G4Event *)
bool SLisEBinFilled(int ebin)
CastorShowerLibraryInfo * hadInfo
std::vector< int > nEvtInBinE
const T & max(const T &a, const T &b)
int FindEnergyBin(double e)
std::pair< std::string, MonitorElement * > entry
std::map< int, G4ThreeVector > PrimaryMomentum
Cos< T >::type cos(const T &t)
std::map< int, G4ThreeVector > PrimaryPosition
Tan< T >::type tan(const T &t)
std::vector< std::vector< std::vector< int > > > nEvtInBinPhi
unsigned int nEvtPerBinPhi
std::map< int, std::set< int > > MapOfSecondaries
CastorShowerLibraryInfo * emInfo
CastorShowerLibraryMaker(const edm::ParameterSet &p)
void GetMissingEnergy(CaloG4HitCollection *, double &, double &)
void GetKinematics(G4PrimaryParticle *, double &px, double &py, double &pz, double &pInit, double &eta, double &phi)
int & SLnEvtInBinPhi(int ebin, int etabin, int phibin)
XYZPointD XYZPoint
point in space with cartesian internal representation
std::vector< std::vector< int > > nEvtInBinEta
void setHitPosition(Point p)
int FindPhiBin(double phi)
void setNEvts(unsigned int n)
void update(const BeginOfJob *run)
This routine will be called when the appropriate signal arrives.
unsigned int nEvtPerBinEta
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
bool SLisPhiBinFilled(int ebin, int etabin, int phibin)
double getTimeSlice() const
std::vector< double > SLEnergyBins
int FindEtaBin(double eta)
math::XYZPoint getEntry() const
std::vector< int > PGParticleIDs
uint32_t getUnitID() const
static uint32_t packIndex(int z, int sector, int zmodule)
void setNphotons(float np)
CastorShowerEvent * emShower
void setNhit(unsigned int i)
double getEnergyDeposit() const