CMS 3D CMS Logo

CastorTestAnalysis.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Forward
4 // Class : CastorTestAnalysis
5 //
6 // Implementation:
7 // <Notes on implementation>
8 //
9 // Original Author: P. Katsas
10 // Created: 02/2007
11 //
12 
15 
18 
19 #include "TFile.h"
20 #include <cmath>
21 #include <iostream>
22 #include <iomanip>
23 
24 //#define EDM_ML_DEBUG
25 
41 };
42 
55 };
56 
58  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("CastorTestAnalysis");
59  verbosity = m_Anal.getParameter<int>("Verbosity");
60  doNTcastorstep = m_Anal.getParameter<int>("StepNtupleFlag");
61  doNTcastorevent = m_Anal.getParameter<int>("EventNtupleFlag");
62  stepNtFileName = m_Anal.getParameter<std::string>("StepNtupleFileName");
63  eventNtFileName = m_Anal.getParameter<std::string>("EventNtupleFileName");
64 
65  if (verbosity > 0) {
66  edm::LogVerbatim("ForwardSim") << std::endl;
67  edm::LogVerbatim("ForwardSim") << "============================================================================";
68  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis:: Initialized as observer";
69  if (doNTcastorstep > 0) {
70  edm::LogVerbatim("ForwardSim") << " Step Ntuple will be created";
71  edm::LogVerbatim("ForwardSim") << " Step Ntuple file: " << stepNtFileName;
72  }
73  if (doNTcastorevent > 0) {
74  edm::LogVerbatim("ForwardSim") << " Event Ntuple will be created";
75  edm::LogVerbatim("ForwardSim") << " Step Ntuple file: " << stepNtFileName;
76  }
77  edm::LogVerbatim("ForwardSim") << "============================================================================";
78  edm::LogVerbatim("ForwardSim") << std::endl;
79  }
80  if (doNTcastorstep > 0)
82  new TNtuple("NTcastorstep", "NTcastorstep", "evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz");
83 
84  if (doNTcastorevent > 0)
85  castoreventntuple = new TNtuple(
86  "NTcastorevent", "NTcastorevent", "evt:ihit:detector:sector:module:enem:totalenergy:hitenergy:x:y:z");
87 }
88 
90  //destructor of CastorTestAnalysis
91 
92  Finish();
93  if (verbosity > 0) {
94  edm::LogVerbatim("ForwardSim") << std::endl << "End of CastorTestAnalysis";
95  }
96 
97  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: End of process";
98 }
99 
100 //=================================================================== per EVENT
101 void CastorTestAnalysis::update(const BeginOfJob* job) { edm::LogVerbatim("ForwardSim") << " Starting new job "; }
102 
103 //==================================================================== per RUN
105  edm::LogVerbatim("ForwardSim") << std::endl << "CastorTestAnalysis: Starting Run";
106  if (doNTcastorstep) {
107  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: output step root file created";
108  TString stepfilename = stepNtFileName;
109  castorOutputStepFile = new TFile(stepfilename, "RECREATE");
110  }
111 
112  if (doNTcastorevent) {
113  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: output event root file created";
114  TString stepfilename = eventNtFileName;
115  castorOutputEventFile = new TFile(stepfilename, "RECREATE");
116  }
117 
118  eventIndex = 0;
119 }
120 
122  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Processing Event Number: " << eventIndex;
123  eventIndex++;
124  stepIndex = 0;
125 }
126 
127 void CastorTestAnalysis::update(const G4Step* aStep) {
128  stepIndex++;
129 
130  if (doNTcastorstep) {
131  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
132  // G4StepPoint * postStepPoint= aStep->GetPostStepPoint();
133  G4double stepL = aStep->GetStepLength();
134  G4double stepE = aStep->GetTotalEnergyDeposit();
135 
136  if (verbosity >= 2)
137  edm::LogVerbatim("ForwardSim") << "Step " << stepL << ", " << stepE;
138 
139  G4Track* theTrack = aStep->GetTrack();
140  G4int theTrackID = theTrack->GetTrackID();
141  G4double theCharge = theTrack->GetDynamicParticle()->GetCharge();
142  // G4String particleType = theTrack->GetDefinition()->GetParticleName();
143  G4int pdgcode = theTrack->GetDefinition()->GetPDGEncoding();
144 
145  const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
146  G4double vpx = vert_mom.x();
147  G4double vpy = vert_mom.y();
148  G4double vpz = vert_mom.z();
149  double eta = 0.5 * log((1. + vpz) / (1. - vpz));
150  double phi = atan2(vpy, vpx);
151 
152  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
153 
154  // Fill-in ntuple
155  // castorsteparray[ntcastors_evt] = (*evt)()->GetEventID();
157  castorsteparray[ntcastors_trackid] = (float)theTrackID;
160  castorsteparray[ntcastors_x] = hitPoint.x();
161  castorsteparray[ntcastors_y] = hitPoint.y();
162  castorsteparray[ntcastors_z] = hitPoint.z();
170 
171  /*
172  edm::LogVerbatim("ForwardSim") << "TrackID: " << theTrackID;
173  edm::LogVerbatim("ForwardSim") << " StepN: "<< theTrack->GetCurrentStepNumber();
174  edm::LogVerbatim("ForwardSim") << " ParentID: " << aStep->GetTrack()->GetParentID();
175  edm::LogVerbatim("ForwardSim") << " PDG: " << pdgcode;
176  edm::LogVerbatim("ForwardSim") << " X,Y,Z (mm): " << theTrack->GetPosition().x() << "," << theTrack->GetPosition().y() << "," << theTrack->GetPosition().z();
177  edm::LogVerbatim("ForwardSim") << " KE (MeV): " << theTrack->GetKineticEnergy();
178  edm::LogVerbatim("ForwardSim") << " Total EDep (MeV): " << aStep->GetTotalEnergyDeposit();
179  edm::LogVerbatim("ForwardSim") << " StepLength (mm): " << aStep->GetStepLength();
180  edm::LogVerbatim("ForwardSim") << " TrackLength (mm): " << theTrack->GetTrackLength();
181 
182  if ( theTrack->GetNextVolume() != 0 )
183  edm::LogVerbatim("ForwardSim") <<" NextVolume: " << theTrack->GetNextVolume()->GetName();
184  else
185  edm::LogVerbatim("ForwardSim") <<" NextVolume: OutOfWorld";
186 
187  if(aStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL)
188  edm::LogVerbatim("ForwardSim") << " ProcessName: "<< aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
189  else
190  edm::LogVerbatim("ForwardSim") <<" ProcessName: UserLimit";
191 
192 
193  edm::LogVerbatim("ForwardSim") << std::endl;
194  */
195 
196 #ifdef EDM_ML_DEBUG
197  if (theTrack->GetNextVolume() != 0)
198  edm::LogVerbatim("ForwardSim") << " NextVolume: " << theTrack->GetNextVolume()->GetName();
199  else
200  edm::LogVerbatim("ForwardSim") << " NextVolume: OutOfWorld";
201 #endif
202 
203  //fill ntuple with step level information
205  }
206 }
207 
208 //================= End of EVENT ===============
210  // Look for the Hit Collection
211  edm::LogVerbatim("ForwardSim") << std::endl
212  << "CastorTest::update(EndOfEvent * evt) - event #" << (*evt)()->GetEventID();
213 
214  // access to the G4 hit collections
215  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
216  edm::LogVerbatim("ForwardSim") << "update(*evt) --> accessed all HC";
217 
218  int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
219 
220  CaloG4HitCollection* theCAFI = (CaloG4HitCollection*)allHC->GetHC(CAFIid);
221 
223  // CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
224 
225  /*
226  unsigned int volumeID=0;
227  int det, zside, sector, zmodule;
228  std::map<int,float,std::less<int> > themap;
229  double totalEnergy = 0;
230  double hitEnergy = 0;
231  double en_in_fi = 0.;
232  double en_in_pl = 0.;
233 */
234  // double en_in_bu = 0.;
235  // double en_in_tu = 0.;
236 
237  if (doNTcastorevent) {
238  eventGlobalHit = 0;
239  // int eventGlobalHit = 0 ;
240 
241  // Check FI TBranch for Hits
242  if (theCAFI->entries() > 0)
243  getCastorBranchData(theCAFI);
244 
245  // Find Primary info:
246  int trackID = 0;
247 #ifdef EDM_ML_DEBUG
248  int particleType = 0;
249 #endif
250  G4PrimaryParticle* thePrim = nullptr;
251  G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
252  edm::LogVerbatim("ForwardSim") << "Event has " << nvertex << " vertex";
253  if (nvertex == 0)
254  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERROR: no vertex";
255 
256  for (int i = 0; i < nvertex; i++) {
257  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
258  if (avertex == nullptr)
259  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: pointer to vertex = 0";
260  edm::LogVerbatim("ForwardSim") << "Vertex number :" << i;
261  int npart = avertex->GetNumberOfParticle();
262  if (npart == 0)
263  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: no primary!";
264  if (thePrim == nullptr)
265  thePrim = avertex->GetPrimary(trackID);
266  }
267 
268  double px = 0., py = 0., pz = 0., pInit = 0;
269 #ifdef EDM_ML_DEBUG
270  double eta = 0., phi = 0.;
271 #endif
272  if (thePrim != nullptr) {
273  px = thePrim->GetPx();
274  py = thePrim->GetPy();
275  pz = thePrim->GetPz();
276  pInit = sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
277  if (pInit == 0) {
278  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: primary has p=0 ";
279 #ifdef EDM_ML_DEBUG
280  } else {
281  float costheta = pz / pInit;
282  float theta = acos(std::min(std::max(costheta, float(-1.)), float(1.)));
283  eta = -log(tan(theta / 2));
284 
285  if (px != 0)
286  phi = atan(py / px);
287 #endif
288  }
289 #ifdef EDM_ML_DEBUG
290  particleType = thePrim->GetPDGcode();
291 #endif
292  } else {
293  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: could not find primary ";
294  }
295 #ifdef EDM_ML_DEBUG
296  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Particle Type " << particleType << " p/eta/phi " << pInit
297  << ", " << eta << ", " << phi;
298 #endif
299  }
300 
301  int iEvt = (*evt)()->GetEventID();
302  if (iEvt < 10)
303  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
304  else if ((iEvt < 100) && (iEvt % 10 == 0))
305  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
306  else if ((iEvt < 1000) && (iEvt % 100 == 0))
307  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
308  else if ((iEvt < 10000) && (iEvt % 1000 == 0))
309  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
310 
311  edm::LogVerbatim("ForwardSim") << std::endl << "===>>> Done writing user histograms ";
312 }
313 
315 
316 //===================================================================
318  int nentries = hc->entries();
319 
320  if (nentries > 0) {
321  unsigned int volumeID = 0;
322  int det = 0, zside, sector, zmodule;
323  std::map<int, float, std::less<int> > themap;
324  double totalEnergy = 0;
325  double hitEnergy = 0;
326  double en_in_sd = 0.;
327 
328  for (int ihit = 0; ihit < nentries; ihit++) {
329  CaloG4Hit* aHit = (*hc)[ihit];
330  totalEnergy += aHit->getEnergyDeposit();
331  }
332 
333  for (int ihit = 0; ihit < nentries; ihit++) {
334  CaloG4Hit* aHit = (*hc)[ihit];
335  volumeID = aHit->getUnitID();
336  hitEnergy = aHit->getEnergyDeposit();
337  en_in_sd += aHit->getEnergyDeposit();
338  // double enEm = aHit->getEM();
339  // double enHad = aHit->getHadr();
340 
341  themap[volumeID] += aHit->getEnergyDeposit();
342  // int det, zside, sector, zmodule;
343  theCastorNumScheme->unpackIndex(volumeID, zside, sector, zmodule);
344 
345  // det = 2 ; // det=2/3 for CAFI/CAPL
346 
348  // castoreventarray[ntcastore_ihit] = (float)ihit;
353  castoreventarray[ntcastore_enem] = en_in_sd;
354  castoreventarray[ntcastore_enhad] = totalEnergy;
359  // castoreventarray[ntcastore_x] = aHit->getEntry().x();
360  // castoreventarray[ntcastore_y] = aHit->getEntry().y();
361  // castoreventarray[ntcastore_z] = aHit->getEntry().z();
362 
364 
365  eventGlobalHit++;
366  }
367  } // nentries > 0
368 }
369 
370 //===================================================================
371 
373  if (doNTcastorstep) {
374  castorOutputStepFile->cd();
375  castorstepntuple->Write();
376  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Ntuple step written";
377  castorOutputStepFile->Close();
378  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Step file closed";
379  }
380 
381  if (doNTcastorevent) {
382  castorOutputEventFile->cd();
383  castoreventntuple->Write("", TObject::kOverwrite);
384  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Ntuple event written";
385  castorOutputEventFile->Close();
386  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Event file closed";
387  }
388 }
mps_fire.i
i
Definition: mps_fire.py:428
MessageLogger.h
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
ntcastore_elements
ntcastore_elements
Definition: CastorTestAnalysis.cc:43
CastorTestAnalysis::eventGlobalHit
int eventGlobalHit
Definition: CastorTestAnalysis.h:105
min
T min(T a, T b)
Definition: MathUtil.h:58
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:66
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
ntcastore_enem
Definition: CastorTestAnalysis.cc:49
CastorTestAnalysis::CastorTestAnalysis
CastorTestAnalysis(const edm::ParameterSet &p)
Definition: CastorTestAnalysis.cc:57
ntcastors_evt
Definition: CastorTestAnalysis.cc:27
ntcastors_elements
ntcastors_elements
Definition: CastorTestAnalysis.cc:26
CastorTestAnalysis::eventIndex
int eventIndex
Definition: CastorTestAnalysis.h:103
ntcastore_enhad
Definition: CastorTestAnalysis.cc:50
ntcastore_y
Definition: CastorTestAnalysis.cc:53
ntcastore_x
Definition: CastorTestAnalysis.cc:52
ntcastors_trackid
Definition: CastorTestAnalysis.cc:28
npart
double npart
Definition: HydjetWrapper.h:46
ntcastors_x
Definition: CastorTestAnalysis.cc:31
ntcastore_hitenergy
Definition: CastorTestAnalysis.cc:51
CastorTestAnalysis::theCastorNumScheme
CastorNumberingScheme * theCastorNumScheme
Definition: CastorTestAnalysis.h:101
ntcastore_module
Definition: CastorTestAnalysis.cc:48
CastorNumberingScheme.h
PVValHelper::eta
Definition: PVValidationHelpers.h:70
CastorTestAnalysis::doNTcastorstep
int doNTcastorstep
Definition: CastorTestAnalysis.h:90
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:78
ntcastors_vpx
Definition: CastorTestAnalysis.cc:38
CastorTestAnalysis::getCastorBranchData
void getCastorBranchData(const CaloG4HitCollection *hc)
Definition: CastorTestAnalysis.cc:317
CastorTestAnalysis::eventNtFileName
std::string eventNtFileName
Definition: CastorTestAnalysis.h:93
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
BeginOfJob
Definition: BeginOfJob.h:8
CastorTestAnalysis::castorOutputStepFile
TFile * castorOutputStepFile
Definition: CastorTestAnalysis.h:96
EndOfEvent
Definition: EndOfEvent.h:6
CastorTestAnalysis::castoreventarray
Float_t castoreventarray[11]
Definition: CastorTestAnalysis.h:108
CastorTestAnalysis::stepNtFileName
std::string stepNtFileName
Definition: CastorTestAnalysis.h:92
CastorTestAnalysis::stepIndex
int stepIndex
Definition: CastorTestAnalysis.h:104
CastorTestAnalysis::castorstepntuple
TNtuple * castorstepntuple
Definition: CastorTestAnalysis.h:98
CastorTestAnalysis::castorsteparray
Float_t castorsteparray[14]
Definition: CastorTestAnalysis.h:107
ntcastors_stepe
Definition: CastorTestAnalysis.cc:35
ntcastors_phi
Definition: CastorTestAnalysis.cc:37
edm::ParameterSet
Definition: ParameterSet.h:47
CastorTestAnalysis::Finish
void Finish()
Definition: CastorTestAnalysis.cc:372
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
ntcastors_y
Definition: CastorTestAnalysis.cc:32
CastorNumberingScheme
Definition: CastorNumberingScheme.h:30
ntcastore_ihit
Definition: CastorTestAnalysis.cc:45
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
CastorNumberingScheme::unpackIndex
static void unpackIndex(const uint32_t &idx, int &z, int &sector, int &zmodule)
Definition: CastorNumberingScheme.cc:193
CastorTestAnalysis::~CastorTestAnalysis
~CastorTestAnalysis() override
Definition: CastorTestAnalysis.cc:89
CastorTestAnalysis::castoreventntuple
TNtuple * castoreventntuple
Definition: CastorTestAnalysis.h:99
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
ntcastors_vpz
Definition: CastorTestAnalysis.cc:40
BeginOfEvent
Definition: BeginOfEvent.h:6
BeginOfRun
Definition: BeginOfRun.h:6
EndOfRun
Definition: EndOfRun.h:6
ntcastors_stepl
Definition: CastorTestAnalysis.cc:34
CaloG4Hit
Definition: CaloG4Hit.h:32
CaloG4Hit::getPosition
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:52
ntcastors_z
Definition: CastorTestAnalysis.cc:33
ntcastors_vpy
Definition: CastorTestAnalysis.cc:39
ntcastore_evt
Definition: CastorTestAnalysis.cc:44
ntcastore_z
Definition: CastorTestAnalysis.cc:54
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
DDAxes::phi
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
CastorTestAnalysis::verbosity
int verbosity
Definition: CastorTestAnalysis.h:89
AlignmentTrackSelector_cfi.theCharge
theCharge
Definition: AlignmentTrackSelector_cfi.py:20
CastorTestAnalysis::doNTcastorevent
int doNTcastorevent
Definition: CastorTestAnalysis.h:91
writedatasetfile.run
run
Definition: writedatasetfile.py:27
ntcastors_pdgcode
Definition: CastorTestAnalysis.cc:30
CastorTestAnalysis::castorOutputEventFile
TFile * castorOutputEventFile
Definition: CastorTestAnalysis.h:95
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
ntcastors_charge
Definition: CastorTestAnalysis.cc:29
ntcastore_detector
Definition: CastorTestAnalysis.cc:46
Point3D.h
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
CaloG4HitCollection
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
Definition: CaloG4HitCollection.h:11
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
PbPb_ZMuSkimMuonDPG_cff.particleType
particleType
Definition: PbPb_ZMuSkimMuonDPG_cff.py:27
ntcastore_sector
Definition: CastorTestAnalysis.cc:47
ntcastors_eta
Definition: CastorTestAnalysis.cc:36
edm::Log
Definition: MessageLogger.h:70
CastorTestAnalysis::update
void update(const BeginOfJob *run) override
This routine will be called when the appropriate signal arrives.
Definition: CastorTestAnalysis.cc:101
CastorTestAnalysis.h