CMS 3D CMS Logo

ZdcTestAnalysis.cc
Go to the documentation of this file.
1 // File: ZdcTestAnalysis.cc
3 // Date: 03.06 Edmundo Garcia
4 // Description: simulation analysis steering code
5 //
10 
14 
22 
23 #include "G4SDManager.hh"
24 #include "G4Step.hh"
25 #include "G4Track.hh"
26 #include "G4Event.hh"
27 #include "G4PrimaryVertex.hh"
28 #include "G4VProcess.hh"
29 #include "G4HCofThisEvent.hh"
30 #include "G4UserEventAction.hh"
31 
32 #include <CLHEP/Units/GlobalSystemOfUnits.h>
33 #include <CLHEP/Units/GlobalPhysicalConstants.h>
34 #include <CLHEP/Random/Randomize.h>
35 
36 #include "TROOT.h"
37 #include "TH1.h"
38 #include "TH2.h"
39 #include "TProfile.h"
40 #include "TNtuple.h"
41 #include "TRandom.h"
42 #include "TLorentzVector.h"
43 #include "TUnixSystem.h"
44 #include "TSystem.h"
45 #include "TMath.h"
46 #include "TF1.h"
47 #include "TFile.h"
48 
49 #include <cassert>
50 #include <cmath>
51 #include <iostream>
52 #include <iomanip>
53 #include <map>
54 #include <string>
55 #include <vector>
56 
57 class ZdcTestAnalysis : public SimWatcher,
58  public Observer<const BeginOfJob*>,
59  public Observer<const BeginOfRun*>,
60  public Observer<const EndOfRun*>,
61  public Observer<const BeginOfEvent*>,
62  public Observer<const EndOfEvent*>,
63  public Observer<const G4Step*> {
64 public:
66  ~ZdcTestAnalysis() override;
67 
68 private:
69  // observer classes
70  void update(const BeginOfJob* run) override;
71  void update(const BeginOfRun* run) override;
72  void update(const EndOfRun* run) override;
73  void update(const BeginOfEvent* evt) override;
74  void update(const EndOfEvent* evt) override;
75  void update(const G4Step* step) override;
76 
77  void finish();
78 
79  int verbosity;
84 
87 
88  TNtuple* zdcstepntuple;
89  TNtuple* zdceventntuple;
90 
92  int stepIndex;
93 
94  Float_t zdcsteparray[18];
95  Float_t zdceventarray[16];
96 
98 };
99 
119 };
120 
138 };
139 
141  //constructor
142  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("ZdcTestAnalysis");
143  verbosity = m_Anal.getParameter<int>("Verbosity");
144  doNTzdcstep = m_Anal.getParameter<int>("StepNtupleFlag");
145  doNTzdcevent = m_Anal.getParameter<int>("EventNtupleFlag");
146  stepNtFileName = m_Anal.getParameter<std::string>("StepNtupleFileName");
147  eventNtFileName = m_Anal.getParameter<std::string>("EventNtupleFileName");
148 
149  if (verbosity > 0)
150  edm::LogVerbatim("ZdcAnalysis") << "\n============================================================================";
151  edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis:: Initialized as observer";
152  if (doNTzdcstep > 0) {
153  edm::LogVerbatim("ZdcAnalysis") << " Step Ntuple will be created";
154  edm::LogVerbatim("ZdcAnalysis") << " Step Ntuple file: " << stepNtFileName;
155  }
156  if (doNTzdcevent > 0) {
157  edm::LogVerbatim("ZdcAnalysis") << " Event Ntuple will be created";
158  edm::LogVerbatim("ZdcAnalysis") << " Step Ntuple file: " << stepNtFileName;
159  }
160  edm::LogVerbatim("ZdcAnalysis") << "============================================================================"
161  << std::endl;
162 
163  if (doNTzdcstep > 0)
164  zdcstepntuple =
165  new TNtuple("NTzdcstep",
166  "NTzdcstep",
167  "evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz:idx:idl:pvtype:ncherphot");
168 
169  if (doNTzdcevent > 0)
171  new TNtuple("NTzdcevent",
172  "NTzdcevent",
173  "evt:ihit:fiberid:zside:subdet:layer:fiber:channel:enem:enhad:hitenergy:x:y:z:time:etot");
174 
175  theZdcNumScheme = nullptr;
176  //theZdcSD = new ZdcSD("ZDCHITSB", new ZdcNumberingScheme());
177 }
178 
180  // destructor
181  finish();
182  delete theZdcNumScheme;
183 }
184 
186  //job
187  edm::LogVerbatim("ZdcAnalysis") << "beggining of job";
188 }
189 
190 //==================================================================== per RUN
192  //run
193 
194  edm::LogVerbatim("ZdcAnalysis") << "\nZdcTestAnalysis: Begining of Run";
195  if (doNTzdcstep) {
196  edm::LogVerbatim("ZdcAnalysis") << "ZDCTestAnalysis: output step file created";
197  TString stepfilename = stepNtFileName;
198  zdcOutputStepFile = new TFile(stepfilename, "RECREATE");
199  }
200 
201  if (doNTzdcevent) {
202  edm::LogVerbatim("ZdcAnalysis") << "ZDCTestAnalysis: output event file created";
203  TString stepfilename = eventNtFileName;
204  zdcOutputEventFile = new TFile(stepfilename, "RECREATE");
205  }
206 
207  eventIndex = 0;
208 }
209 
211  //event
212  edm::LogVerbatim("ZdcAnalysis") << "ZdcTest: Processing Event Number: " << eventIndex;
213  eventIndex++;
214  stepIndex = 0;
215 }
216 
217 void ZdcTestAnalysis::update(const G4Step* aStep) {
218  //step;
219  stepIndex++;
220 
221  if (doNTzdcstep) {
222  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
223  // G4StepPoint * postStepPoint= aStep->GetPostStepPoint();
224  G4double stepL = aStep->GetStepLength();
225  G4double stepE = aStep->GetTotalEnergyDeposit();
226 
227  if (verbosity >= 2)
228  edm::LogVerbatim("ZdcAnalysis") << "Step " << stepL << "," << stepE;
229 
230  G4Track* theTrack = aStep->GetTrack();
231  G4int theTrackID = theTrack->GetTrackID();
232  G4double theCharge = theTrack->GetDynamicParticle()->GetCharge();
233  G4String particleType = theTrack->GetDefinition()->GetParticleName();
234  G4int pdgcode = theTrack->GetDefinition()->GetPDGEncoding();
235 
236  const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
237  G4double vpx = vert_mom.x();
238  G4double vpy = vert_mom.y();
239  G4double vpz = vert_mom.z();
240  double eta = 0.5 * log((1. + vpz) / (1. - vpz));
241  double phi = atan2(vpy, vpx);
242 
243  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
244  G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
245 
246  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
247  int idx = touch->GetReplicaNumber(0);
248  int idLayer = -1;
249  int thePVtype = -1;
250 
251  int historyDepth = touch->GetHistoryDepth();
252 
253  if (historyDepth > 0) {
254  std::vector<int> theReplicaNumbers(historyDepth);
255  std::vector<G4VPhysicalVolume*> thePhysicalVolumes(historyDepth);
256  std::vector<G4String> thePVnames(historyDepth);
257  std::vector<G4LogicalVolume*> theLogicalVolumes(historyDepth);
258  std::vector<G4String> theLVnames(historyDepth);
259  std::vector<G4Material*> theMaterials(historyDepth);
260  std::vector<G4String> theMaterialNames(historyDepth);
261 
262  for (int jj = 0; jj < historyDepth; jj++) {
263  theReplicaNumbers[jj] = touch->GetReplicaNumber(jj);
264  thePhysicalVolumes[jj] = touch->GetVolume(jj);
265  thePVnames[jj] = thePhysicalVolumes[jj]->GetName();
266  theLogicalVolumes[jj] = thePhysicalVolumes[jj]->GetLogicalVolume();
267  theLVnames[jj] = theLogicalVolumes[jj]->GetName();
268  theMaterials[jj] = theLogicalVolumes[jj]->GetMaterial();
269  theMaterialNames[jj] = theMaterials[jj]->GetName();
270  if (verbosity >= 2)
271  edm::LogVerbatim("ZdcAnalysis") << " GHD " << jj << ": " << theReplicaNumbers[jj] << "," << thePVnames[jj]
272  << "," << theLVnames[jj] << "," << theMaterialNames[jj];
273  }
274 
275  idLayer = theReplicaNumbers[1];
276  if (thePVnames[0] == "ZDC_EMLayer")
277  thePVtype = 1;
278  else if (thePVnames[0] == "ZDC_EMAbsorber")
279  thePVtype = 2;
280  else if (thePVnames[0] == "ZDC_EMFiber")
281  thePVtype = 3;
282  else if (thePVnames[0] == "ZDC_HadLayer")
283  thePVtype = 7;
284  else if (thePVnames[0] == "ZDC_HadAbsorber")
285  thePVtype = 8;
286  else if (thePVnames[0] == "ZDC_HadFiber")
287  thePVtype = 9;
288  else if (thePVnames[0] == "ZDC_PhobosLayer")
289  thePVtype = 11;
290  else if (thePVnames[0] == "ZDC_PhobosAbsorber")
291  thePVtype = 12;
292  else if (thePVnames[0] == "ZDC_PhobosFiber")
293  thePVtype = 13;
294  else {
295  thePVtype = 0;
296  if (verbosity >= 2)
297  edm::LogVerbatim("ZdcAnalysis") << " pvtype=0 hd=" << historyDepth << " " << theReplicaNumbers[0] << ","
298  << thePVnames[0] << "," << theLVnames[0] << "," << theMaterialNames[0];
299  }
300  } else if (historyDepth == 0) {
301  int theReplicaNumber = touch->GetReplicaNumber(0);
302  G4VPhysicalVolume* thePhysicalVolume = touch->GetVolume(0);
303  const G4String& thePVname = thePhysicalVolume->GetName();
304  G4LogicalVolume* theLogicalVolume = thePhysicalVolume->GetLogicalVolume();
305  const G4String& theLVname = theLogicalVolume->GetName();
306  G4Material* theMaterial = theLogicalVolume->GetMaterial();
307  const G4String& theMaterialName = theMaterial->GetName();
308  if (verbosity >= 2)
309  edm::LogVerbatim("ZdcAnalysis") << " hd=0 " << theReplicaNumber << "," << thePVname << "," << theLVname << ","
310  << theMaterialName;
311  } else {
312  edm::LogVerbatim("ZdcAnalysis") << " hd<0: hd=" << historyDepth;
313  }
314 
315  double NCherPhot = -1;
317  zdcsteparray[ntzdcs_trackid] = (float)theTrackID;
319  zdcsteparray[ntzdcs_pdgcode] = (float)pdgcode;
320  zdcsteparray[ntzdcs_x] = hitPoint.x();
321  zdcsteparray[ntzdcs_y] = hitPoint.y();
322  zdcsteparray[ntzdcs_z] = hitPoint.z();
323  zdcsteparray[ntzdcs_stepl] = stepL;
324  zdcsteparray[ntzdcs_stepe] = stepE;
327  zdcsteparray[ntzdcs_vpx] = vpx;
328  zdcsteparray[ntzdcs_vpy] = vpy;
329  zdcsteparray[ntzdcs_vpz] = vpz;
331  zdcsteparray[ntzdcs_idl] = (float)idLayer;
332  zdcsteparray[ntzdcs_pvtype] = thePVtype;
333  zdcsteparray[ntzdcs_ncherphot] = NCherPhot;
335  }
336 }
337 
339  //end of event
340 
341  // Look for the Hit Collection
342  edm::LogVerbatim("ZdcAnalysis") << "\nZdcTest::upDate(const EndOfEvent * evt) - event #" << (*evt)()->GetEventID()
343  << "\n # of aSteps followed in event = " << stepIndex;
344 
345  // access to the G4 hit collections
346  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
347  edm::LogVerbatim("ZdcAnalysis") << " accessed all HC";
348 
349  int theZDCHCid = G4SDManager::GetSDMpointer()->GetCollectionID("ZDCHITS");
350  edm::LogVerbatim("ZdcAnalysis") << " - theZDCHCid = " << theZDCHCid;
351 
352  CaloG4HitCollection* theZDCHC = (CaloG4HitCollection*)allHC->GetHC(theZDCHCid);
353  edm::LogVerbatim("ZdcAnalysis") << " - theZDCHC = " << theZDCHC;
354 
355  if (!theZdcNumScheme) {
357  }
358 
359  float ETot = 0., SEnergy = 0.;
360  int maxTime = 0;
361  int fiberID = 0;
362  unsigned int unsignedfiberID = 0;
363  std::map<int, float, std::less<int> > energyInFibers;
364  std::map<int, float, std::less<int> > primaries;
365  float totalEnergy = 0;
366  int nentries = theZDCHC->entries();
367  edm::LogVerbatim("ZdcAnalysis") << " theZDCHC has " << nentries << " entries";
368 
369  if (doNTzdcevent) {
370  if (nentries > 0) {
371  for (int ihit = 0; ihit < nentries; ihit++) {
372  CaloG4Hit* caloHit = (*theZDCHC)[ihit];
373  totalEnergy += caloHit->getEnergyDeposit();
374  }
375 
376  for (int ihit = 0; ihit < nentries; ihit++) {
377  CaloG4Hit* aHit = (*theZDCHC)[ihit];
378  fiberID = aHit->getUnitID();
379  unsignedfiberID = aHit->getUnitID();
380  double enEm = aHit->getEM();
381  double enHad = aHit->getHadr();
382  math::XYZPoint hitPoint = aHit->getPosition();
383  double hitEnergy = aHit->getEnergyDeposit();
384  if (verbosity >= 1)
385  edm::LogVerbatim("ZdcAnalysis")
386  << " entry #" << ihit << ": fiberID=0x" << std::hex << fiberID << std::dec << "; enEm=" << enEm
387  << "; enHad=" << enHad << "; hitEnergy=" << hitEnergy << "z=" << hitPoint.z();
388  energyInFibers[fiberID] += enEm + enHad;
389  primaries[aHit->getTrackID()] += enEm + enHad;
390  float time = aHit->getTimeSliceID();
391  if (time > maxTime)
392  maxTime = (int)time;
393 
394  int thesubdet, thelayer, thefiber, thechannel, thez;
395  theZdcNumScheme->unpackZdcIndex(fiberID, thesubdet, thelayer, thefiber, thechannel, thez);
396  int unsignedsubdet, unsignedlayer, unsignedfiber, unsignedchannel, unsignedz;
398  unsignedfiberID, unsignedsubdet, unsignedlayer, unsignedfiber, unsignedchannel, unsignedz);
399 
400  // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
401  // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
402  // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
403  // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
404 
407  zdceventarray[ntzdce_fiberid] = (float)fiberID;
409  zdceventarray[ntzdce_subdet] = (float)thesubdet;
410  zdceventarray[ntzdce_layer] = (float)thelayer;
411  zdceventarray[ntzdce_fiber] = (float)thefiber;
412  zdceventarray[ntzdce_channel] = (float)thechannel;
413  zdceventarray[ntzdce_enem] = enEm;
414  zdceventarray[ntzdce_enhad] = enHad;
415  zdceventarray[ntzdce_hitenergy] = hitEnergy;
416  zdceventarray[ntzdce_x] = hitPoint.x();
417  zdceventarray[ntzdce_y] = hitPoint.y();
418  zdceventarray[ntzdce_z] = hitPoint.z();
420  zdceventarray[ntzdce_etot] = totalEnergy;
422  }
423 
424  for (std::map<int, float, std::less<int> >::iterator is = energyInFibers.begin(); is != energyInFibers.end();
425  is++) {
426  ETot = (*is).second;
427  SEnergy += ETot;
428  }
429 
430  // Find Primary info:
431  int trackID = 0;
432  G4PrimaryParticle* thePrim = nullptr;
433  G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
434  edm::LogVerbatim("ZdcAnalysis") << "Event has " << nvertex << " vertex";
435  if (nvertex == 0)
436  edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERROR: no vertex";
437 
438  for (int i = 0; i < nvertex; i++) {
439  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
440  if (avertex == nullptr) {
441  edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: pointer to vertex = 0";
442  } else {
443  edm::LogVerbatim("ZdcAnalysis") << "Vertex number :" << i;
444  int npart = avertex->GetNumberOfParticle();
445  if (npart == 0)
446  edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: no primary!";
447  if (thePrim == nullptr)
448  thePrim = avertex->GetPrimary(trackID);
449  }
450  }
451 
452  double px = 0., py = 0., pz = 0.;
453  double pInit = 0.;
454 
455  if (thePrim != nullptr) {
456  px = thePrim->GetPx();
457  py = thePrim->GetPy();
458  pz = thePrim->GetPz();
459  pInit = sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
460  if (pInit == 0) {
461  edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: primary has p=0 ";
462  }
463  } else {
464  edm::LogVerbatim("ZdcAnalysis") << "ZdcTest End Of Event ERR: could not find primary ";
465  }
466 
467  } // nentries > 0
468 
469  } // createNTzdcevent
470 
471  int iEvt = (*evt)()->GetEventID();
472  if (iEvt < 10)
473  edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
474  else if ((iEvt < 100) && (iEvt % 10 == 0))
475  edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
476  else if ((iEvt < 1000) && (iEvt % 100 == 0))
477  edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
478  else if ((iEvt < 10000) && (iEvt % 1000 == 0))
479  edm::LogVerbatim("ZdcAnalysis") << " ZdcTest Event " << iEvt;
480 }
481 
483 
485  if (doNTzdcstep) {
486  zdcOutputStepFile->cd();
487  zdcstepntuple->Write();
488  edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Ntuple step written for event: " << eventIndex;
489  zdcOutputStepFile->Close();
490  edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Step file closed";
491  }
492 
493  if (doNTzdcevent) {
494  zdcOutputEventFile->cd();
495  zdceventntuple->Write("", TObject::kOverwrite);
496  edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Ntuple event written for event: " << eventIndex;
497  zdcOutputEventFile->Close();
498  edm::LogVerbatim("ZdcAnalysis") << "ZdcTestAnalysis: Event file closed";
499  }
500 }
501 
504 
Log< level::Info, true > LogVerbatim
void update(const BeginOfJob *run) override
This routine will be called when the appropriate signal arrives.
int getTrackID() const
Definition: CaloG4Hit.h:64
#define DEFINE_SIMWATCHER(type)
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
Float_t zdceventarray[16]
TFile * zdcOutputEventFile
double getEM() const
Definition: CaloG4Hit.h:55
double npart
Definition: HydjetWrapper.h:46
ntzdce_elements
std::string stepNtFileName
double getHadr() const
Definition: CaloG4Hit.h:58
TFile * zdcOutputStepFile
~ZdcTestAnalysis() override
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:52
T sqrt(T t)
Definition: SSEVec.h:19
TNtuple * zdceventntuple
ZdcNumberingScheme * theZdcNumScheme
TNtuple * zdcstepntuple
ntzdcs_elements
int getTimeSliceID() const
Definition: CaloG4Hit.h:68
ZdcTestAnalysis(const edm::ParameterSet &p)
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
static void unpackZdcIndex(const unsigned int &idx, int &subDet, int &layer, int &fiber, int &channel, int &z)
double getEnergyDeposit() const
Definition: CaloG4Hit.h:79
uint32_t getUnitID() const
Definition: CaloG4Hit.h:66
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
std::string eventNtFileName
Float_t zdcsteparray[18]
step
Definition: StallMonitor.cc:98
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29