CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 
13 #include "SimG4Core/Notification/interface/Observer.h"
14 #include "SimG4Core/Notification/interface/BeginOfJob.h"
15 #include "SimG4Core/Notification/interface/BeginOfRun.h"
16 #include "SimG4Core/Notification/interface/EndOfRun.h"
17 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
18 #include "SimG4Core/Notification/interface/EndOfEvent.h"
19 #include "SimG4Core/Watcher/interface/SimWatcher.h"
20 
21 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
22 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
23 #include "SimG4CMS/Forward/interface/CastorNumberingScheme.h"
24 
25 #include "DataFormats/Math/interface/Point3D.h"
26 #include "FWCore/MessageLogger/interface/MessageLogger.h"
27 #include "FWCore/ParameterSet/interface/ParameterSet.h"
28 
29 #include "G4SDManager.hh"
30 #include "G4Step.hh"
31 #include "G4Track.hh"
32 #include "G4Event.hh"
33 #include "G4PrimaryVertex.hh"
34 #include "G4VProcess.hh"
35 #include "G4HCofThisEvent.hh"
36 #include "G4UserEventAction.hh"
37 
38 #include <CLHEP/Units/GlobalSystemOfUnits.h>
39 #include <CLHEP/Units/GlobalPhysicalConstants.h>
40 #include <CLHEP/Random/Randomize.h>
41 
42 #include "TROOT.h"
43 #include "TFile.h"
44 #include "TH1.h"
45 #include "TH2.h"
46 #include "TProfile.h"
47 #include "TNtuple.h"
48 #include "TRandom.h"
49 #include "TLorentzVector.h"
50 #include "TUnixSystem.h"
51 #include "TSystem.h"
52 #include "TMath.h"
53 #include "TF1.h"
54 
55 #include <cassert>
56 #include <cmath>
57 #include <iostream>
58 #include <iomanip>
59 #include <map>
60 #include <string>
61 #include <vector>
62 
63 //#define EDM_ML_DEBUG
64 
65 class CastorTestAnalysis : public SimWatcher,
66  public Observer<const BeginOfJob *>,
67  public Observer<const BeginOfRun *>,
68  public Observer<const EndOfRun *>,
69  public Observer<const BeginOfEvent *>,
70  public Observer<const EndOfEvent *>,
71  public Observer<const G4Step *> {
72 public:
73  CastorTestAnalysis(const edm::ParameterSet &p);
74  ~CastorTestAnalysis() override;
75 
76 private:
77  // observer classes
78  void update(const BeginOfJob *run) override;
79  void update(const BeginOfRun *run) override;
80  void update(const EndOfRun *run) override;
81  void update(const BeginOfEvent *evt) override;
82  void update(const EndOfEvent *evt) override;
83  void update(const G4Step *step) override;
84 
85 private:
86  void getCastorBranchData(const CaloG4HitCollection *hc);
87  void Finish();
88 
89  int verbosity;
90  int doNTcastorstep;
91  int doNTcastorevent;
92  std::string stepNtFileName;
93  std::string eventNtFileName;
94 
95  TFile *castorOutputEventFile;
96  TFile *castorOutputStepFile;
97 
98  TNtuple *castorstepntuple;
99  TNtuple *castoreventntuple;
100 
101  CastorNumberingScheme *theCastorNumScheme;
102 
103  int eventIndex;
104  int stepIndex;
105  int eventGlobalHit;
106 
107  Float_t castorsteparray[14];
108  Float_t castoreventarray[11];
109 };
110 
111 enum ntcastors_elements {
112  ntcastors_evt,
113  ntcastors_trackid,
114  ntcastors_charge,
115  ntcastors_pdgcode,
116  ntcastors_x,
117  ntcastors_y,
118  ntcastors_z,
119  ntcastors_stepl,
120  ntcastors_stepe,
121  ntcastors_eta,
122  ntcastors_phi,
123  ntcastors_vpx,
124  ntcastors_vpy,
125  ntcastors_vpz
126 };
127 
128 enum ntcastore_elements {
129  ntcastore_evt,
130  ntcastore_ihit,
131  ntcastore_detector,
132  ntcastore_sector,
133  ntcastore_module,
134  ntcastore_enem,
135  ntcastore_enhad,
136  ntcastore_hitenergy,
137  ntcastore_x,
138  ntcastore_y,
139  ntcastore_z
140 };
141 
142 CastorTestAnalysis::CastorTestAnalysis(const edm::ParameterSet &p) {
143  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("CastorTestAnalysis");
144  verbosity = m_Anal.getParameter<int>("Verbosity");
145  doNTcastorstep = m_Anal.getParameter<int>("StepNtupleFlag");
146  doNTcastorevent = m_Anal.getParameter<int>("EventNtupleFlag");
147  stepNtFileName = m_Anal.getParameter<std::string>("StepNtupleFileName");
148  eventNtFileName = m_Anal.getParameter<std::string>("EventNtupleFileName");
149 
150  if (verbosity > 0) {
151  edm::LogVerbatim("ForwardSim") << std::endl;
152  edm::LogVerbatim("ForwardSim") << "============================================================================";
153  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis:: Initialized as observer";
154  if (doNTcastorstep > 0) {
155  edm::LogVerbatim("ForwardSim") << " Step Ntuple will be created";
156  edm::LogVerbatim("ForwardSim") << " Step Ntuple file: " << stepNtFileName;
157  }
158  if (doNTcastorevent > 0) {
159  edm::LogVerbatim("ForwardSim") << " Event Ntuple will be created";
160  edm::LogVerbatim("ForwardSim") << " Step Ntuple file: " << stepNtFileName;
161  }
162  edm::LogVerbatim("ForwardSim") << "============================================================================";
163  edm::LogVerbatim("ForwardSim") << std::endl;
164  }
165  if (doNTcastorstep > 0)
166  castorstepntuple =
167  new TNtuple("NTcastorstep", "NTcastorstep", "evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz");
168 
169  if (doNTcastorevent > 0)
170  castoreventntuple = new TNtuple(
171  "NTcastorevent", "NTcastorevent", "evt:ihit:detector:sector:module:enem:totalenergy:hitenergy:x:y:z");
172 }
173 
174 CastorTestAnalysis::~CastorTestAnalysis() {
175  //destructor of CastorTestAnalysis
176 
177  Finish();
178  if (verbosity > 0) {
179  edm::LogVerbatim("ForwardSim") << std::endl << "End of CastorTestAnalysis";
180  }
181 
182  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: End of process";
183 }
184 
185 //=================================================================== per EVENT
186 void CastorTestAnalysis::update(const BeginOfJob *job) { edm::LogVerbatim("ForwardSim") << " Starting new job "; }
187 
188 //==================================================================== per RUN
189 void CastorTestAnalysis::update(const BeginOfRun *run) {
190  edm::LogVerbatim("ForwardSim") << std::endl << "CastorTestAnalysis: Starting Run";
191  if (doNTcastorstep) {
192  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: output step root file created";
193  TString stepfilename = stepNtFileName;
194  castorOutputStepFile = new TFile(stepfilename, "RECREATE");
195  }
196 
197  if (doNTcastorevent) {
198  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: output event root file created";
199  TString stepfilename = eventNtFileName;
200  castorOutputEventFile = new TFile(stepfilename, "RECREATE");
201  }
202 
203  eventIndex = 0;
204 }
205 
206 void CastorTestAnalysis::update(const BeginOfEvent *evt) {
207  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Processing Event Number: " << eventIndex;
208  eventIndex++;
209  stepIndex = 0;
210 }
211 
212 void CastorTestAnalysis::update(const G4Step *aStep) {
213  stepIndex++;
214 
215  if (doNTcastorstep) {
216  G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
217  // G4StepPoint * postStepPoint= aStep->GetPostStepPoint();
218  G4double stepL = aStep->GetStepLength();
219  G4double stepE = aStep->GetTotalEnergyDeposit();
220 
221  if (verbosity >= 2)
222  edm::LogVerbatim("ForwardSim") << "Step " << stepL << ", " << stepE;
223 
224  G4Track *theTrack = aStep->GetTrack();
225  G4int theTrackID = theTrack->GetTrackID();
226  G4double theCharge = theTrack->GetDynamicParticle()->GetCharge();
227  // G4String particleType = theTrack->GetDefinition()->GetParticleName();
228  G4int pdgcode = theTrack->GetDefinition()->GetPDGEncoding();
229 
230  const G4ThreeVector &vert_mom = theTrack->GetVertexMomentumDirection();
231  G4double vpx = vert_mom.x();
232  G4double vpy = vert_mom.y();
233  G4double vpz = vert_mom.z();
234  double eta = 0.5 * log((1. + vpz) / (1. - vpz));
235  double phi = atan2(vpy, vpx);
236 
237  const G4ThreeVector &hitPoint = preStepPoint->GetPosition();
238 
239  // Fill-in ntuple
240  // castorsteparray[ntcastors_evt] = (*evt)()->GetEventID();
241  castorsteparray[ntcastors_evt] = (float)eventIndex;
242  castorsteparray[ntcastors_trackid] = (float)theTrackID;
243  castorsteparray[ntcastors_charge] = theCharge;
244  castorsteparray[ntcastors_pdgcode] = pdgcode;
245  castorsteparray[ntcastors_x] = hitPoint.x();
246  castorsteparray[ntcastors_y] = hitPoint.y();
247  castorsteparray[ntcastors_z] = hitPoint.z();
248  castorsteparray[ntcastors_stepl] = stepL;
249  castorsteparray[ntcastors_stepe] = stepE;
250  castorsteparray[ntcastors_eta] = eta;
251  castorsteparray[ntcastors_phi] = phi;
252  castorsteparray[ntcastors_vpx] = vpx;
253  castorsteparray[ntcastors_vpy] = vpy;
254  castorsteparray[ntcastors_vpz] = vpz;
255 
256  /*
257  edm::LogVerbatim("ForwardSim") << "TrackID: " << theTrackID;
258  edm::LogVerbatim("ForwardSim") << " StepN: "<< theTrack->GetCurrentStepNumber();
259  edm::LogVerbatim("ForwardSim") << " ParentID: " << aStep->GetTrack()->GetParentID();
260  edm::LogVerbatim("ForwardSim") << " PDG: " << pdgcode;
261  edm::LogVerbatim("ForwardSim") << " X,Y,Z (mm): " << theTrack->GetPosition().x() << "," << theTrack->GetPosition().y() << "," << theTrack->GetPosition().z();
262  edm::LogVerbatim("ForwardSim") << " KE (MeV): " << theTrack->GetKineticEnergy();
263  edm::LogVerbatim("ForwardSim") << " Total EDep (MeV): " << aStep->GetTotalEnergyDeposit();
264  edm::LogVerbatim("ForwardSim") << " StepLength (mm): " << aStep->GetStepLength();
265  edm::LogVerbatim("ForwardSim") << " TrackLength (mm): " << theTrack->GetTrackLength();
266 
267  if ( theTrack->GetNextVolume() != 0 )
268  edm::LogVerbatim("ForwardSim") <<" NextVolume: " << theTrack->GetNextVolume()->GetName();
269  else
270  edm::LogVerbatim("ForwardSim") <<" NextVolume: OutOfWorld";
271 
272  if(aStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL)
273  edm::LogVerbatim("ForwardSim") << " ProcessName: "<< aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
274  else
275  edm::LogVerbatim("ForwardSim") <<" ProcessName: UserLimit";
276 
277 
278  edm::LogVerbatim("ForwardSim") << std::endl;
279  */
280 
281 #ifdef EDM_ML_DEBUG
282  if (theTrack->GetNextVolume() != 0)
283  edm::LogVerbatim("ForwardSim") << " NextVolume: " << theTrack->GetNextVolume()->GetName();
284  else
285  edm::LogVerbatim("ForwardSim") << " NextVolume: OutOfWorld";
286 #endif
287 
288  //fill ntuple with step level information
289  castorstepntuple->Fill(castorsteparray);
290  }
291 }
292 
293 //================= End of EVENT ===============
294 void CastorTestAnalysis::update(const EndOfEvent *evt) {
295  // Look for the Hit Collection
296  edm::LogVerbatim("ForwardSim") << std::endl
297  << "CastorTest::update(EndOfEvent * evt) - event #" << (*evt)()->GetEventID();
298 
299  // access to the G4 hit collections
300  G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
301  edm::LogVerbatim("ForwardSim") << "update(*evt) --> accessed all HC";
302 
303  int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
304 
305  CaloG4HitCollection *theCAFI = (CaloG4HitCollection *)allHC->GetHC(CAFIid);
306 
307  theCastorNumScheme = new CastorNumberingScheme();
308  // CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
309 
310  /*
311  unsigned int volumeID=0;
312  int det, zside, sector, zmodule;
313  std::map<int,float,std::less<int> > themap;
314  double totalEnergy = 0;
315  double hitEnergy = 0;
316  double en_in_fi = 0.;
317  double en_in_pl = 0.;
318 */
319  // double en_in_bu = 0.;
320  // double en_in_tu = 0.;
321 
322  if (doNTcastorevent) {
323  eventGlobalHit = 0;
324  // int eventGlobalHit = 0 ;
325 
326  // Check FI TBranch for Hits
327  if (theCAFI->entries() > 0)
328  getCastorBranchData(theCAFI);
329 
330  // Find Primary info:
331  int trackID = 0;
332 #ifdef EDM_ML_DEBUG
333  int particleType = 0;
334 #endif
335  G4PrimaryParticle *thePrim = nullptr;
336  G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
337  edm::LogVerbatim("ForwardSim") << "Event has " << nvertex << " vertex";
338  if (nvertex == 0)
339  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERROR: no vertex";
340 
341  for (int i = 0; i < nvertex; i++) {
342  G4PrimaryVertex *avertex = (*evt)()->GetPrimaryVertex(i);
343  if (avertex == nullptr)
344  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: pointer to vertex = 0";
345  edm::LogVerbatim("ForwardSim") << "Vertex number :" << i;
346  int npart = avertex->GetNumberOfParticle();
347  if (npart == 0)
348  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: no primary!";
349  if (thePrim == nullptr)
350  thePrim = avertex->GetPrimary(trackID);
351  }
352 
353  double px = 0., py = 0., pz = 0., pInit = 0;
354 #ifdef EDM_ML_DEBUG
355  double eta = 0., phi = 0.;
356 #endif
357  if (thePrim != nullptr) {
358  px = thePrim->GetPx();
359  py = thePrim->GetPy();
360  pz = thePrim->GetPz();
361  pInit = sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
362  if (pInit == 0) {
363  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: primary has p=0 ";
364 #ifdef EDM_ML_DEBUG
365  } else {
366  float costheta = pz / pInit;
367  float theta = acos(std::min(std::max(costheta, float(-1.)), float(1.)));
368  eta = -log(tan(theta / 2));
369 
370  if (px != 0)
371  phi = atan(py / px);
372 #endif
373  }
374 #ifdef EDM_ML_DEBUG
375  particleType = thePrim->GetPDGcode();
376 #endif
377  } else {
378  edm::LogVerbatim("ForwardSim") << "CASTORTest End Of Event ERR: could not find primary ";
379  }
380 #ifdef EDM_ML_DEBUG
381  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Particle Type " << particleType << " p/eta/phi " << pInit
382  << ", " << eta << ", " << phi;
383 #endif
384  }
385 
386  int iEvt = (*evt)()->GetEventID();
387  if (iEvt < 10)
388  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
389  else if ((iEvt < 100) && (iEvt % 10 == 0))
390  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
391  else if ((iEvt < 1000) && (iEvt % 100 == 0))
392  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
393  else if ((iEvt < 10000) && (iEvt % 1000 == 0))
394  edm::LogVerbatim("ForwardSim") << " CastorTest Event " << iEvt;
395 
396  edm::LogVerbatim("ForwardSim") << std::endl << "===>>> Done writing user histograms ";
397 }
398 
399 void CastorTestAnalysis::update(const EndOfRun *run) { ; }
400 
401 //===================================================================
402 void CastorTestAnalysis::getCastorBranchData(const CaloG4HitCollection *hc) {
403  int nentries = hc->entries();
404 
405  if (nentries > 0) {
406  unsigned int volumeID = 0;
407  int det = 0, zside, sector, zmodule;
408  std::map<int, float, std::less<int> > themap;
409  double totalEnergy = 0;
410  double hitEnergy = 0;
411  double en_in_sd = 0.;
412 
413  for (int ihit = 0; ihit < nentries; ihit++) {
414  CaloG4Hit *aHit = (*hc)[ihit];
415  totalEnergy += aHit->getEnergyDeposit();
416  }
417 
418  for (int ihit = 0; ihit < nentries; ihit++) {
419  CaloG4Hit *aHit = (*hc)[ihit];
420  volumeID = aHit->getUnitID();
421  hitEnergy = aHit->getEnergyDeposit();
422  en_in_sd += aHit->getEnergyDeposit();
423  // double enEm = aHit->getEM();
424  // double enHad = aHit->getHadr();
425 
426  themap[volumeID] += aHit->getEnergyDeposit();
427  // int det, zside, sector, zmodule;
428  theCastorNumScheme->unpackIndex(volumeID, zside, sector, zmodule);
429 
430  // det = 2 ; // det=2/3 for CAFI/CAPL
431 
432  castoreventarray[ntcastore_evt] = (float)eventIndex;
433  // castoreventarray[ntcastore_ihit] = (float)ihit;
434  castoreventarray[ntcastore_ihit] = (float)eventGlobalHit;
435  castoreventarray[ntcastore_detector] = (float)det;
436  castoreventarray[ntcastore_sector] = (float)sector;
437  castoreventarray[ntcastore_module] = (float)zmodule;
438  castoreventarray[ntcastore_enem] = en_in_sd;
439  castoreventarray[ntcastore_enhad] = totalEnergy;
440  castoreventarray[ntcastore_hitenergy] = hitEnergy;
441  castoreventarray[ntcastore_x] = aHit->getPosition().x();
442  castoreventarray[ntcastore_y] = aHit->getPosition().y();
443  castoreventarray[ntcastore_z] = aHit->getPosition().z();
444  // castoreventarray[ntcastore_x] = aHit->getEntry().x();
445  // castoreventarray[ntcastore_y] = aHit->getEntry().y();
446  // castoreventarray[ntcastore_z] = aHit->getEntry().z();
447 
448  castoreventntuple->Fill(castoreventarray);
449 
450  eventGlobalHit++;
451  }
452  } // nentries > 0
453 }
454 
455 //===================================================================
456 
457 void CastorTestAnalysis::Finish() {
458  if (doNTcastorstep) {
459  castorOutputStepFile->cd();
460  castorstepntuple->Write();
461  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Ntuple step written";
462  castorOutputStepFile->Close();
463  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Step file closed";
464  }
465 
466  if (doNTcastorevent) {
467  castorOutputEventFile->cd();
468  castoreventntuple->Write("", TObject::kOverwrite);
469  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Ntuple event written";
470  castorOutputEventFile->Close();
471  edm::LogVerbatim("ForwardSim") << "CastorTestAnalysis: Event file closed";
472  }
473 }
474 
475 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
476 #include "FWCore/PluginManager/interface/ModuleDef.h"
477 
478 DEFINE_SIMWATCHER(CastorTestAnalysis);
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
ntcastors_elements
ntcastors_elements
Definition: CastorTestAnalysis.cc:111
DEFINE_SIMWATCHER
#define DEFINE_SIMWATCHER(type)
Definition: SimWatcherFactory.h:13
CaloG4Hit::getPosition
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:52
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
mps_fire.i
i
Definition: mps_fire.py:428
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
npart
double npart
Definition: HydjetWrapper.h:46
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
ntcastore_elements
ntcastore_elements
Definition: CastorTestAnalysis.cc:128
CastorNumberingScheme::unpackIndex
static void unpackIndex(const uint32_t &idx, int &z, int &sector, int &zmodule)
Definition: CastorNumberingScheme.cc:193
CastorTestAnalysis::theCastorNumScheme
CastorNumberingScheme * theCastorNumScheme
Definition: CastorTestAnalysis.cc:101
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
CastorTestAnalysis::update
void update(const BeginOfJob *run) override
This routine will be called when the appropriate signal arrives.
Definition: CastorTestAnalysis.cc:186
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
CastorTestAnalysis::CastorTestAnalysis
CastorTestAnalysis(const edm::ParameterSet &p)
Definition: CastorTestAnalysis.cc:142
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
CaloG4HitCollection
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
Definition: CaloG4HitCollection.h:11
step
step
Definition: StallMonitor.cc:98
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:66
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
CastorTestAnalysis::getCastorBranchData
void getCastorBranchData(const CaloG4HitCollection *hc)
Definition: CastorTestAnalysis.cc:402
HLT_FULL_cff.particleType
tuple particleType
Definition: HLT_FULL_cff.py:18671
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:79
Generated for CMSSW Reference Manual by doxygen 1.8.5