CMS 3D CMS Logo
BscTest.cc
Go to the documentation of this file.
1 // system include files
2 #include <cmath>
3 #include <iostream>
4 #include <iomanip>
5 #include <map>
6 #include <string>
7 #include <vector>
8 
9 //
10 #include "FWCore/Framework/interface/Event.h"
11 #include "FWCore/Framework/interface/EventSetup.h"
12 #include "FWCore/Framework/interface/Frameworkfwd.h"
13 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
14 #include "FWCore/MessageLogger/interface/MessageLogger.h"
15 #include "FWCore/ParameterSet/interface/ParameterSet.h"
16 
17 #include "SimG4Core/Notification/interface/TrackInformation.h"
18 #include "SimG4Core/Notification/interface/Observer.h"
19 #include "SimG4Core/Notification/interface/BeginOfJob.h"
20 #include "SimG4Core/Notification/interface/BeginOfRun.h"
21 #include "SimG4Core/Notification/interface/EndOfRun.h"
22 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
23 #include "SimG4Core/Notification/interface/EndOfEvent.h"
24 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
25 #include "SimG4Core/Notification/interface/EndOfTrack.h"
26 #include "SimG4Core/Watcher/interface/SimWatcher.h"
27 
28 // to retreive hits
29 #include "SimG4CMS/Forward/interface/BscG4Hit.h"
30 #include "SimG4CMS/Forward/interface/BscNumberingScheme.h"
31 #include "SimG4CMS/Forward/interface/BscG4HitCollection.h"
32 
33 // G4 stuff
34 #include "G4SDManager.hh"
35 #include "G4Step.hh"
36 #include "G4Track.hh"
37 #include "G4VProcess.hh"
38 #include "G4HCofThisEvent.hh"
39 #include "G4UserEventAction.hh"
40 #include "G4TransportationManager.hh"
41 #include "G4ProcessManager.hh"
42 #include "G4VTouchable.hh"
43 
44 #include <CLHEP/Vector/ThreeVector.h>
45 #include <CLHEP/Vector/LorentzVector.h>
46 #include <CLHEP/Random/Randomize.h>
47 #include <CLHEP/Units/SystemOfUnits.h>
48 #include <CLHEP/Units/GlobalPhysicalConstants.h>
49 
50 // ----------------------------------------------------------------
51 // Includes needed for Root ntupling
52 //
53 #include "TROOT.h"
54 #include "TFile.h"
55 #include "TH1.h"
56 #include "TH2.h"
57 #include "TProfile.h"
58 #include "TNtuple.h"
59 #include "TRandom.h"
60 #include "TLorentzVector.h"
61 #include "TUnixSystem.h"
62 #include "TSystem.h"
63 #include "TMath.h"
64 #include "TF1.h"
65 #include "TObjArray.h"
66 #include "TObjString.h"
67 #include "TNamed.h"
68 
69 //#define EDM_ML_DEBUG
70 
71 class BscAnalysisHistManager : public TNamed {
72 public:
73  BscAnalysisHistManager(const TString& managername);
74  ~BscAnalysisHistManager() override;
75 
76  TH1F* getHisto(Int_t Number);
77  TH1F* getHisto(const TObjString& histname);
78 
79  TH2F* getHisto2(Int_t Number);
80  TH2F* getHisto2(const TObjString& histname);
81 
82  void writeToFile(const TString& fOutputFile, const TString& fRecreateFile);
83 
84 private:
85  void bookHistos();
86  void storeWeights();
87  void histInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup);
88  void histInit(
89  const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup);
90 
91  const char* fTypeTitle;
92  TObjArray* fHistArray;
93  TObjArray* fHistNamesArray;
94 };
95 
96 class BscTest : public SimWatcher,
97  public Observer<const BeginOfJob*>,
98  public Observer<const BeginOfRun*>,
99  public Observer<const EndOfRun*>,
100  public Observer<const BeginOfEvent*>,
101  public Observer<const BeginOfTrack*>,
102  public Observer<const G4Step*>,
103  public Observer<const EndOfTrack*>,
104  public Observer<const EndOfEvent*> {
105 public:
106  BscTest(const edm::ParameterSet& p);
107  ~BscTest() override;
108 
109 private:
110  // observer classes
111  void update(const BeginOfJob* run) override;
112  void update(const BeginOfRun* run) override;
113  void update(const EndOfRun* run) override;
114  void update(const BeginOfEvent* evt) override;
115  void update(const BeginOfTrack* trk) override;
116  void update(const G4Step* step) override;
117  void update(const EndOfTrack* trk) override;
118  void update(const EndOfEvent* evt) override;
119 
120  int iev;
121  int itrk;
122  G4double entot0, tracklength0;
123 
124  // Utilities to get detector levels during a step
125 
126  int detLevels(const G4VTouchable*) const;
127  G4String detName(const G4VTouchable*, int, int) const;
128  void detectorLevel(const G4VTouchable*, int&, int*, G4String*) const;
129 
130  double rinCalo, zinCalo;
131  int lastTrackID;
132  int verbosity;
133 
134  // sumEnerDeposit - all deposited energy on all steps ; sumStepl - length in steel !!!
135  G4double sumEnerDeposit, sumStepl, sumStepc;
136  // numofpart - # particles produced along primary track
137  int numofpart;
138  // last point of the primary track
139  G4ThreeVector lastpo;
140 
141  // z:
142  double z1, z2, z3, z4;
143 
144 private:
145  Float_t bsceventarray[1];
146  TNtuple* bsceventntuple;
147  TFile bscOutputFile;
148  int whichevent;
149 
150  BscAnalysisHistManager* TheHistManager; //Histogram Manager of the analysis
151  std::string fDataLabel; // Data type label
152  std::string fOutputFile; // The output file name
153  std::string fRecreateFile; // Recreate the file flag, default="RECREATE"
154 };
155 
156 enum ntbsc_elements { ntbsc_evt };
157 
158 //================================================================
159 BscTest::BscTest(const edm::ParameterSet& p) {
160  //constructor
161  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("BscTest");
162  verbosity = m_Anal.getParameter<int>("Verbosity");
163  //verbosity = 1;
164 
165  fDataLabel = m_Anal.getParameter<std::string>("FDataLabel");
166  fOutputFile = m_Anal.getParameter<std::string>("FOutputFile");
167  fRecreateFile = m_Anal.getParameter<std::string>("FRecreateFile");
168 
169  if (verbosity > 0) {
170  edm::LogVerbatim("BscTest") << "============================================================================";
171  edm::LogVerbatim("BscTest") << "BscTestconstructor :: Initialized as observer";
172  }
173  // Initialization:
174 
175  bsceventntuple = new TNtuple("NTbscevent", "NTbscevent", "evt");
176  whichevent = 0;
177  TheHistManager = new BscAnalysisHistManager(fDataLabel);
178 
179  if (verbosity > 0) {
180  edm::LogVerbatim("BscTest") << "BscTest constructor :: Initialized BscAnalysisHistManager";
181  }
182 }
183 
184 BscTest::~BscTest() {
185  // delete UserNtuples;
186 
187  TFile bscOutputFile("newntbsc.root", "RECREATE");
188  edm::LogVerbatim("BscTest") << "Bsc output root file has been created";
189  bsceventntuple->Write();
190  edm::LogVerbatim("BscTest") << ", written";
191  bscOutputFile.Close();
192  edm::LogVerbatim("BscTest") << ", closed";
193  delete bsceventntuple;
194  edm::LogVerbatim("BscTest") << ", and deleted";
195 
196  //------->while end
197 
198  // Write histograms to file
199  TheHistManager->writeToFile(fOutputFile, fRecreateFile);
200  if (verbosity > 0) {
201  edm::LogVerbatim("BscTest") << std::endl << "BscTest Destructor --------> End of BscTest : ";
202  }
203 
204  edm::LogVerbatim("BscTest") << "BscTest: End of process";
205 }
206 
207 //================================================================
208 // Histoes:
209 //-----------------------------------------------------------------------------
210 
211 BscAnalysisHistManager::BscAnalysisHistManager(const TString& managername) {
212  // The Constructor
213 
214  fTypeTitle = managername;
215  fHistArray = new TObjArray(); // Array to store histos
216  fHistNamesArray = new TObjArray(); // Array to store histos's names
217 
218  bookHistos();
219 
220  fHistArray->Compress(); // Removes empty space
221  fHistNamesArray->Compress();
222 }
223 //-----------------------------------------------------------------------------
224 
225 BscAnalysisHistManager::~BscAnalysisHistManager() {
226  // The Destructor
227 
228  if (fHistArray) {
229  fHistArray->Delete();
230  delete fHistArray;
231  }
232 
233  if (fHistNamesArray) {
234  fHistNamesArray->Delete();
235  delete fHistNamesArray;
236  }
237 }
238 //-----------------------------------------------------------------------------
239 
240 void BscAnalysisHistManager::bookHistos() {
241  // at Start: (mm)
242  histInit("TrackPhi", "Primary Phi", 100, 0., 360.);
243  histInit("TrackTheta", "Primary Theta", 100, 0., 180.);
244  histInit("TrackP", "Track XY position Z+ ", 80, -80., 80., 80, -80., 80.);
245  histInit("TrackM", "Track XY position Z-", 80, -80., 80., 80, -80., 80.);
246  histInit("DetIDs", "Track DetId - vs +", 16, -0.5, 15.5, 16, 15.5, 31.5);
247 }
248 
249 //-----------------------------------------------------------------------------
250 
251 void BscAnalysisHistManager::writeToFile(const TString& fOutputFile, const TString& fRecreateFile) {
252  //Write to file = fOutputFile
253 
254  edm::LogVerbatim("BscTest") << "================================================================";
255  edm::LogVerbatim("BscTest") << " Write this Analysis to File " << fOutputFile;
256  edm::LogVerbatim("BscTest") << "================================================================";
257 
258  TFile* file = new TFile(fOutputFile, fRecreateFile);
259 
260  fHistArray->Write();
261  file->Close();
262 }
263 //-----------------------------------------------------------------------------
264 
265 void BscAnalysisHistManager::histInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup) {
266  // Add histograms and histograms names to the array
267 
268  char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
269  strcpy(newtitle, title);
270  strcat(newtitle, " (");
271  strcat(newtitle, fTypeTitle);
272  strcat(newtitle, ") ");
273  fHistArray->AddLast((new TH1F(name, newtitle, nbinsx, xlow, xup)));
274  fHistNamesArray->AddLast(new TObjString(name));
275 }
276 //-----------------------------------------------------------------------------
277 
278 void BscAnalysisHistManager::histInit(
279  const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup) {
280  // Add histograms and histograms names to the array
281 
282  char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
283  strcpy(newtitle, title);
284  strcat(newtitle, " (");
285  strcat(newtitle, fTypeTitle);
286  strcat(newtitle, ") ");
287  fHistArray->AddLast((new TH2F(name, newtitle, nbinsx, xlow, xup, nbinsy, ylow, yup)));
288  fHistNamesArray->AddLast(new TObjString(name));
289 }
290 //-----------------------------------------------------------------------------
291 
292 TH1F* BscAnalysisHistManager::getHisto(Int_t Number) {
293  // Get a histogram from the array with index = Number
294 
295  if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
296  return (TH1F*)(fHistArray->At(Number));
297 
298  } else {
299  edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!!!getHisto!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
300  edm::LogVerbatim("BscTest") << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number;
301  edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
302 
303  return (TH1F*)(fHistArray->At(0));
304  }
305 }
306 //-----------------------------------------------------------------------------
307 
308 TH2F* BscAnalysisHistManager::getHisto2(Int_t Number) {
309  // Get a histogram from the array with index = Number
310 
311  if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
312  return (TH2F*)(fHistArray->At(Number));
313 
314  } else {
315  edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!getHisto2!!!!!!!!!!!!!!!!!!!!!!!!!!!";
316  edm::LogVerbatim("BscTest") << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number;
317  edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
318 
319  return (TH2F*)(fHistArray->At(0));
320  }
321 }
322 //-----------------------------------------------------------------------------
323 
324 TH1F* BscAnalysisHistManager::getHisto(const TObjString& histname) {
325  // Get a histogram from the array with name = histname
326 
327  Int_t index = fHistNamesArray->IndexOf(&histname);
328  return getHisto(index);
329 }
330 //-----------------------------------------------------------------------------
331 
332 TH2F* BscAnalysisHistManager::getHisto2(const TObjString& histname) {
333  // Get a histogram from the array with name = histname
334 
335  Int_t index = fHistNamesArray->IndexOf(&histname);
336  return getHisto2(index);
337 }
338 //-----------------------------------------------------------------------------
339 
340 void BscAnalysisHistManager::storeWeights() {
341  // Add structure to each histogram to store the weights
342 
343  for (int i = 0; i < fHistArray->GetEntries(); i++) {
344  ((TH1F*)(fHistArray->At(i)))->Sumw2();
345  }
346 }
347 
348 //==================================================================== per JOB
349 void BscTest::update(const BeginOfJob* job) {
350  //job
351  edm::LogVerbatim("BscTest") << "BscTest:beggining of job";
352  ;
353 }
354 
355 //==================================================================== per RUN
356 void BscTest::update(const BeginOfRun* run) {
357  //run
358 
359  edm::LogVerbatim("BscTest") << std::endl << "BscTest:: Begining of Run";
360 }
361 
362 void BscTest::update(const EndOfRun* run) { ; }
363 
364 //=================================================================== per EVENT
365 void BscTest::update(const BeginOfEvent* evt) {
366  iev = (*evt)()->GetEventID();
367  if (verbosity > 0) {
368  edm::LogVerbatim("BscTest") << "BscTest:update Event number = " << iev;
369  }
370  whichevent++;
371 }
372 
373 //=================================================================== per Track
374 void BscTest::update(const BeginOfTrack* trk) {
375  itrk = (*trk)()->GetTrackID();
376  if (verbosity > 1) {
377  edm::LogVerbatim("BscTest") << "BscTest:update BeginOfTrack number = " << itrk;
378  }
379  if (itrk == 1) {
380  sumEnerDeposit = 0.;
381  numofpart = 0;
382  sumStepl = 0.;
383  sumStepc = 0.;
384  tracklength0 = 0.;
385  }
386 }
387 
388 //=================================================================== per EndOfTrack
389 void BscTest::update(const EndOfTrack* trk) {
390  itrk = (*trk)()->GetTrackID();
391  if (verbosity > 1) {
392  edm::LogVerbatim("BscTest") << "BscTest:update EndOfTrack number = " << itrk;
393  }
394  if (itrk == 1) {
395  G4double tracklength = (*trk)()->GetTrackLength(); // Accumulated track length
396 
397  TheHistManager->getHisto("SumEDep")->Fill(sumEnerDeposit);
398  TheHistManager->getHisto("TrackL")->Fill(tracklength);
399 
400  // direction !!!
401  G4ThreeVector vert_mom = (*trk)()->GetVertexMomentumDirection();
402  G4ThreeVector vert_pos = (*trk)()->GetVertexPosition(); // vertex ,where this track was created
403 
404  // last step information
405  const G4Step* aStep = (*trk)()->GetStep();
406  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
407  lastpo = preStepPoint->GetPosition();
408  }
409 }
410 
411 // ====================================================
412 
413 //=================================================================== each STEP
414 void BscTest::update(const G4Step* aStep) {
415  // ==========================================================================
416 
417  if (verbosity > 2) {
418  G4int stepnumber = aStep->GetTrack()->GetCurrentStepNumber();
419  edm::LogVerbatim("BscTest") << "BscTest:update Step number = " << stepnumber;
420  }
421  // track on aStep: !
422  G4Track* theTrack = aStep->GetTrack();
423  TrackInformation* trkInfo = dynamic_cast<TrackInformation*>(theTrack->GetUserInformation());
424  if (trkInfo == nullptr) {
425  edm::LogVerbatim("BscTest") << "BscTest on aStep: No trk info !!!! abort ";
426  }
427  G4int id = theTrack->GetTrackID();
428  G4String particleType = theTrack->GetDefinition()->GetParticleName(); // !!!
429  G4int parentID = theTrack->GetParentID(); // !!!
430  G4TrackStatus trackstatus = theTrack->GetTrackStatus(); // !!!
431  G4double tracklength = theTrack->GetTrackLength(); // Accumulated track length
432  G4ThreeVector trackmom = theTrack->GetMomentum();
433  G4double entot = theTrack->GetTotalEnergy(); // !!! deposited on step
434  G4int curstepnumber = theTrack->GetCurrentStepNumber();
435  G4double stepl = aStep->GetStepLength();
436  G4double EnerDeposit = aStep->GetTotalEnergyDeposit();
437  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
438  const G4ThreeVector& preposition = preStepPoint->GetPosition();
439  G4ThreeVector prelocalpoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(preposition);
440  G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
441  const G4String& prename = currentPV->GetName();
442 
443  const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
444  int pre_levels = detLevels(pre_touch);
445  G4String name1[20];
446  int copyno1[20];
447  for (int i = 0; i < 20; ++i) {
448  name1[i] = "";
449  copyno1[i] = 0;
450  }
451  if (pre_levels > 0) {
452  detectorLevel(pre_touch, pre_levels, copyno1, name1);
453  }
454 
455  if (id == 1) {
456  // on 1-st step:
457  if (curstepnumber == 1) {
458  entot0 = entot;
459  //UserNtuples->fillg519(entot0,1.);
460  }
461 
462  // on every step:
463 
464  // for Copper:
465  if (prename == "SBST") {
466  sumStepc += stepl;
467  // =========
468  }
469  // for ststeel:
470  // if(prename == "SBSTs") {
471  if (prename == "SBSTs") {
472  sumStepl += stepl;
473  // =========
474  }
475  // =========
476  // =========
477 
478  // exclude last track point if it is in SD (MI was started their)
479  if (trackstatus != 2) {
480  // for SD: Si Det.: SISTATION:SIPLANE:(SIDETL+BOUNDDET +SIDETR + CERAMDET)
481  if (prename == "SIDETL" || prename == "SIDETR") {
482  if (prename == "SIDETL") {
483  //UserNtuples->fillg569(EnerDeposit,1.);
484  }
485  if (prename == "SIDETR") {
486  //UserNtuples->fillg570(EnerDeposit,1.);
487  }
488 
489  G4String posname = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
490  if ((prename == "SIDETL" && posname != "SIDETL") || (prename == "SIDETR" && posname != "SIDETR")) {
491  if (name1[2] == "SISTATION") {
492  //UserNtuples->fillg539(copyno1[2],1.);
493  }
494  if (name1[3] == "SIPLANE") {
495  //UserNtuples->fillg540(copyno1[3],1.);
496  }
497 
498  if (prename == "SIDETL") {
499  //UserNtuples->fillg541(EnerDeposit,1.);
500  //UserNtuples->fillg561(numbcont,1.);
501  if (copyno1[2] < 2) {
502  //UserNtuples->fillg571(dx,1.);
503  } else if (copyno1[2] < 3) {
504  //UserNtuples->fillg563(dx,1.);
505  if (copyno1[3] < 2) {
506  } else if (copyno1[3] < 3) {
507  //UserNtuples->fillg572(dx,1.);
508  } else if (copyno1[3] < 4) {
509  //UserNtuples->fillg573(dx,1.);
510  } else if (copyno1[3] < 5) {
511  //UserNtuples->fillg574(dx,1.);
512  } else if (copyno1[3] < 6) {
513  //UserNtuples->fillg575(dx,1.);
514  } else if (copyno1[3] < 7) {
515  //UserNtuples->fillg576(dx,1.);
516  } else if (copyno1[3] < 8) {
517  //UserNtuples->fillg577(dx,1.);
518  } else if (copyno1[3] < 9) {
519  //UserNtuples->fillg578(dx,1.);
520  } else if (copyno1[3] < 10) {
521  //UserNtuples->fillg579(dx,1.);
522  }
523  } else if (copyno1[2] < 4) {
524  //UserNtuples->fillg565(dx,1.);
525  } else if (copyno1[2] < 5) {
526  //UserNtuples->fillg567(dx,1.);
527  }
528  }
529  if (prename == "SIDETR") {
530  //UserNtuples->fillg542(EnerDeposit,1.);
531  //UserNtuples->fillg562(numbcont,1.);
532  if (copyno1[2] < 2) {
533  //UserNtuples->fillg581(dy,1.);
534  } else if (copyno1[2] < 3) {
535  //UserNtuples->fillg564(dy,1.);
536  if (copyno1[3] < 2) {
537  } else if (copyno1[3] < 3) {
538  //UserNtuples->fillg582(dy,1.);
539  } else if (copyno1[3] < 4) {
540  //UserNtuples->fillg583(dy,1.);
541  } else if (copyno1[3] < 5) {
542  //UserNtuples->fillg584(dy,1.);
543  } else if (copyno1[3] < 6) {
544  //UserNtuples->fillg585(dy,1.);
545  } else if (copyno1[3] < 7) {
546  //UserNtuples->fillg586(dy,1.);
547  } else if (copyno1[3] < 8) {
548  //UserNtuples->fillg587(dy,1.);
549  } else if (copyno1[3] < 9) {
550  //UserNtuples->fillg588(dy,1.);
551  } else if (copyno1[3] < 10) {
552  //UserNtuples->fillg589(dy,1.);
553  }
554  } else if (copyno1[2] < 4) {
555  //UserNtuples->fillg566(dy,1.);
556  } else if (copyno1[2] < 5) {
557  //UserNtuples->fillg568(dy,1.);
558  }
559  }
560  }
561  }
562  // end of prenames SIDETL // SIDETR
563  }
564  // end of trackstatus != 2
565 
566  sumEnerDeposit += EnerDeposit;
567  if (trackstatus == 2) {
568  // primary track length
569  // //UserNtuples->fillg508(tracklength,1.);
570  tracklength0 = tracklength;
571  }
572  }
573  // end of primary track !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574 
575  if (parentID == 1 && curstepnumber == 1) {
576  // particles deposit their energy along primary track
577  numofpart += 1;
578  if (prename == "SBST") {
579  //UserNtuples->filld225(vert_pos.y(),vert_pos.z(),1.);
580  }
581  if (prename == "SBSTs") {
582  //UserNtuples->filld226(vert_pos.y(),vert_pos.z(),1.);
583  }
584  }
585 }
586 // ==========================================================================
587 // ==========================================================================
588 int BscTest::detLevels(const G4VTouchable* touch) const {
589  //Return number of levels
590  if (touch)
591  return ((touch->GetHistoryDepth()) + 1);
592  else
593  return 0;
594 }
595 // ==========================================================================
596 
597 G4String BscTest::detName(const G4VTouchable* touch, int level, int currentlevel) const {
598  //Go down to current level
599  if (level > 0 && level >= currentlevel) {
600  int ii = level - currentlevel;
601  return touch->GetVolume(ii)->GetName();
602  } else {
603  return "NotFound";
604  }
605 }
606 
607 void BscTest::detectorLevel(const G4VTouchable* touch, int& level, int* copyno, G4String* name) const {
608  //Get name and copy numbers
609  if (level > 0) {
610  for (int ii = 0; ii < level; ii++) {
611  int i = level - ii - 1;
612  G4VPhysicalVolume* pv = touch->GetVolume(i);
613  if (pv != nullptr)
614  name[ii] = pv->GetName();
615  else
616  name[ii] = "Unknown";
617  copyno[ii] = touch->GetReplicaNumber(i);
618  }
619  }
620 }
621 // ==========================================================================
622 
623 //=================================================================== End Of Event
624 void BscTest::update(const EndOfEvent* evt) {
625  // ==========================================================================
626 
627  if (verbosity > 1) {
628  iev = (*evt)()->GetEventID();
629  edm::LogVerbatim("BscTest") << "BscTest:update EndOfEvent = " << iev;
630  }
631  // Fill-in ntuple
632  bsceventarray[ntbsc_evt] = (float)whichevent;
633 
634  //
635  int trackID = 0;
636  G4PrimaryParticle* thePrim = nullptr;
637 
638  // prim.vertex:
639  G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
640  if (nvertex != 1)
641  edm::LogVerbatim("BscTest") << "BscTest: My warning: NumberOfPrimaryVertex != 1 --> = " << nvertex;
642 
643  for (int i = 0; i < nvertex; i++) {
644  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
645  if (avertex == nullptr) {
646  edm::LogVerbatim("BscTest") << "BscTest End Of Event ERR: pointer to vertex = 0";
647  continue;
648  }
649  G4int npart = avertex->GetNumberOfParticle();
650  if (npart != 1)
651  edm::LogVerbatim("BscTest") << "BscTest: My warning: NumberOfPrimaryPart != 1 --> = " << npart;
652  if (npart == 0)
653  edm::LogVerbatim("BscTest") << "BscTest End Of Event ERR: no NumberOfParticle";
654 
655  if (thePrim == nullptr)
656  thePrim = avertex->GetPrimary(trackID);
657 
658  if (thePrim != nullptr) {
659  // primary vertex:
660  G4double vx = 0., vy = 0., vz = 0.;
661  vx = avertex->GetX0();
662  vy = avertex->GetY0();
663  vz = avertex->GetZ0();
664  //UserNtuples->fillh01(vx);
665  //UserNtuples->fillh02(vy);
666  //UserNtuples->fillh03(vz);
667  TheHistManager->getHisto("VtxX")->Fill(vx);
668  TheHistManager->getHisto("VtxY")->Fill(vy);
669  TheHistManager->getHisto("VtxZ")->Fill(vz);
670  }
671  }
672  // prim.vertex loop end
673 
674  //=========================== thePrim != 0 ================================================================================
675  if (thePrim != nullptr) {
676  //
677  // number of secondary particles deposited their energy along primary track
678  //UserNtuples->fillg518(numofpart,1.);
679  if (lastpo.z() < z4 && lastpo.perp() < 100.) {
680  //UserNtuples->fillg536(numofpart,1.);
681  }
682  //
683 
684  // direction !!!
685  G4ThreeVector mom = thePrim->GetMomentum();
686 
687  double phi = atan2(mom.y(), mom.x());
688  if (phi < 0.)
689  phi += twopi;
690  double phigrad = phi * 180. / pi;
691 
692  double th = mom.theta();
693  double eta = -log(tan(th / 2));
694  TheHistManager->getHisto("PrimaryEta")->Fill(eta);
695  TheHistManager->getHisto("PrimaryPhigrad")->Fill(phigrad);
696  TheHistManager->getHisto("PrimaryTh")->Fill(th * 180. / pi);
697 
698  TheHistManager->getHisto("PrimaryLastpoZ")->Fill(lastpo.z());
699  if (lastpo.z() < z4) {
700  TheHistManager->getHisto("PrimaryLastpoX")->Fill(lastpo.x());
701  TheHistManager->getHisto("PrimaryLastpoY")->Fill(lastpo.y());
702  }
703  if (numofpart > 4) {
704  TheHistManager->getHisto("XLastpoNumofpart")->Fill(lastpo.x());
705  TheHistManager->getHisto("YLastpoNumofpart")->Fill(lastpo.y());
706  }
707 
708  // ==========================================================================
709 
710  // hit map for Bsc
711  // ==================================
712 
713  std::map<int, float, std::less<int> > themap;
714  std::map<int, float, std::less<int> > themap1;
715 
716  std::map<int, float, std::less<int> > themapxy;
717  std::map<int, float, std::less<int> > themapz;
718  // access to the G4 hit collections: -----> this work OK:
719 #ifdef EDM_ML_DEBUG
720  edm::LogVerbatim("BscTest") << "1";
721 #endif
722  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
723 #ifdef EDM_ML_DEBUG
724  edm::LogVerbatim("BscTest") << "2";
725 #endif
726  if (verbosity > 0) {
727  edm::LogVerbatim("BscTest") << "BscTest: accessed all HC";
728  ;
729  }
730  int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("BSCHits");
731 
732  BscG4HitCollection* theCAFI = (BscG4HitCollection*)allHC->GetHC(CAFIid);
733  if (verbosity > 0) {
734  edm::LogVerbatim("BscTest") << "BscTest: theCAFI->entries = " << theCAFI->entries();
735  }
736  int varia; // = 0 -all; =1 - MI; =2 - noMI
737  //varia = 0;
738  if (lastpo.z() < z4) {
739  varia = 1;
740  } else {
741  varia = 2;
742  } // no MI end:
743  int nhits = theCAFI->entries();
744  for (int j = 0; j < nhits; j++) {
745  BscG4Hit* aHit = (*theCAFI)[j];
746  const CLHEP::Hep3Vector& hitPoint = aHit->getEntry();
747  double zz = hitPoint.z();
748  TheHistManager->getHisto("zHits")->Fill(zz);
749  if (tracklength0 > 8300.)
750  TheHistManager->getHisto("zHitsTrLoLe")->Fill(zz);
751  }
752 
753  if (varia == 2) {
754  //int nhit11 = 0, nhit12 = 0, nhit13 = 0;
755  double totallosenergy = 0.;
756  for (int j = 0; j < nhits; j++) {
757  BscG4Hit* aHit = (*theCAFI)[j];
758 
759  const CLHEP::Hep3Vector& hitEntryLocalPoint = aHit->getEntryLocalP();
760  const CLHEP::Hep3Vector& hitExitLocalPoint = aHit->getExitLocalP();
761  const CLHEP::Hep3Vector& hitPoint = aHit->getEntry();
762  //int trackIDhit = aHit->getTrackID();
763  unsigned int unitID = aHit->getUnitID();
764  double losenergy = aHit->getEnergyLoss();
765 
766  double zz = hitPoint.z();
767 
768  TheHistManager->getHisto("zHitsnoMI")->Fill(zz);
769 
770  if (verbosity > 2) {
771  edm::LogVerbatim("BscTest") << "BscTest:zHits = " << zz;
772  }
773 
774  themap[unitID] += losenergy;
775  totallosenergy += losenergy;
776 
777  int zside;
778  //int sector;
779  BscNumberingScheme::unpackBscIndex(unitID);
780  zside = (unitID & 32) >> 5;
781  //sector = (unitID & 7);
782 
783  //
784  //=======================================
785  G4ThreeVector middle = (hitExitLocalPoint + hitEntryLocalPoint) / 2.;
786  themapz[unitID] = hitPoint.z() + middle.z();
787  //=======================================
788  // Y
789  if (zside == 1) {
790  //UserNtuples->fillg24(losenergy,1.);
791  if (losenergy > 0.00003) {
792  themap1[unitID] += 1.;
793  }
794  }
795  //X
796  else if (zside == 2) {
797  //UserNtuples->fillg25(losenergy,1.);
798  if (losenergy > 0.00005) {
799  themap1[unitID] += 1.;
800  }
801  }
802  // }
803  //
804  /*
805  if (sector == 1) {
806  nhit11 += 1;
807  //UserNtuples->fillg33(rr,1.);
808  //UserNtuples->fillg11(yy,1.);
809  }
810  if (sector == 2) {
811  nhit12 += 1;
812  //UserNtuples->fillg34(rr,1.);
813  //UserNtuples->fillg86(yy,1.);
814  }
815  if (sector == 3) {
816  nhit13 += 1;
817  //UserNtuples->fillg35(rr,1.);
818  //UserNtuples->fillg87(yy,1.);
819  }
820 
821  if (lastpo.z() < z4 && lastpo.perp() < 120.) {
822  // MIonly:
823  //UserNtuples->fillg16(lastpo.z(),1.);
824  //UserNtuples->fillg18(zz,1.);
825  // Station I
826  if (zz < z2) {
827  //UserNtuples->fillg54(dx,1.);
828  //UserNtuples->fillg55(dy,1.);
829  }
830  // Station II
831  if (zz < z3 && zz > z2) {
832  //UserNtuples->fillg50(dx,1.);
833  //UserNtuples->fillg51(dy,1.);
834  }
835  // Station III
836  if (zz < z4 && zz > z3) {
837  //UserNtuples->fillg64(dx,1.);
838  //UserNtuples->fillg65(dy,1.);
839  //UserNtuples->filld209(xx,yy,1.);
840  }
841  } else {
842  // no MIonly:
843  //UserNtuples->fillg17(lastpo.z(),1.);
844  //UserNtuples->fillg19(zz,1.);
845  //UserNtuples->fillg74(incidentEnergyHit,1.);
846  //UserNtuples->fillg75(float(trackIDhit),1.);
847  // Station I
848  if (zz < z2) {
849  //UserNtuples->fillg56(dx,1.);
850  //UserNtuples->fillg57(dy,1.);
851  //UserNtuples->fillg20(numofpart,1.);
852  //UserNtuples->fillg21(sumEnerDeposit,1.);
853  if (zside == 1) {
854  //UserNtuples->fillg26(losenergy,1.);
855  }
856  if (zside == 2) {
857  //UserNtuples->fillg76(losenergy,1.);
858  }
859  if (trackIDhit == 1) {
860  //UserNtuples->fillg70(dx,1.);
861  //UserNtuples->fillg71(incidentEnergyHit,1.);
862  //UserNtuples->fillg79(losenergy,1.);
863  } else {
864  //UserNtuples->fillg82(dx,1.);
865  }
866  }
867  // Station II
868  if (zz < z3 && zz > z2) {
869  //UserNtuples->fillg52(dx,1.);
870  //UserNtuples->fillg53(dy,1.);
871  //UserNtuples->fillg22(numofpart,1.);
872  //UserNtuples->fillg23(sumEnerDeposit,1.);
873  //UserNtuples->fillg80(incidentEnergyHit,1.);
874  //UserNtuples->fillg81(float(trackIDhit),1.);
875  if (zside == 1) {
876  //UserNtuples->fillg27(losenergy,1.);
877  }
878  if (zside == 2) {
879  //UserNtuples->fillg77(losenergy,1.);
880  }
881  if (trackIDhit == 1) {
882  //UserNtuples->fillg72(dx,1.);
883  //UserNtuples->fillg73(incidentEnergyHit,1.);
884  } else {
885  //UserNtuples->fillg83(dx,1.);
886  }
887  }
888  // Station III
889  if (zz < z4 && zz > z3) {
890  if (zside == 1) {
891  //UserNtuples->fillg28(losenergy,1.);
892  }
893  if (zside == 2) {
894  //UserNtuples->fillg78(losenergy,1.);
895  }
896  }
897  }
898  */
899  } // MIonly or noMIonly ENDED
900  if (totallosenergy == 0.0) {
901  edm::LogVerbatim("BscTest") << "BscTest: number of hits = " << theCAFI->entries();
902  for (int j = 0; j < nhits; j++) {
903  BscG4Hit* aHit = (*theCAFI)[j];
904  double losenergy = aHit->getEnergyLoss();
905  edm::LogVerbatim("BscTest") << " j hits = " << j << "losenergy = " << losenergy;
906  }
907  }
908  // FIBRE Hit collected analysis
909  /*
910  double totalEnergy = 0.;
911  int nhitsX = 0, nhitsY = 0, nsumhit = 0;
912  for (int sector = 1; sector < 4; sector++) {
913  int nhitsecX = 0, nhitsecY = 0;
914  for (int zmodule = 1; zmodule < 11; zmodule++) {
915  for (int zside = 1; zside < 3; zside++) {
916  int det = 1; // nhit = 0;
917  // int sScale = 20;
918  int index = BscNumberingScheme::packBscIndex(det, zside, sector);
919  double theTotalEnergy = themap[index];
920  // X planes
921  if (zside < 2) {
922  //UserNtuples->fillg47(theTotalEnergy,1.);
923  if (theTotalEnergy > 0.00003) {
924  nhitsX += 1;
925  // nhitsecX += themap1[index];
926  // nhit=1;
927  }
928  }
929  // Y planes
930  else {
931  //UserNtuples->fillg49(theTotalEnergy,1.);
932  if (theTotalEnergy > 0.00005) {
933  nhitsY += 1;
934  // nhitsecY += themap1[index];
935  // nhit=1;
936  }
937  }
938 
939  totalEnergy += themap[index];
940  } // for
941  } // for
942  //UserNtuples->fillg39(nhitsecY,1.);
943  if (nhitsecX > 10 || nhitsecY > 10) {
944  nsumhit += 1;
945  //UserNtuples->fillp213(float(sector),float(1.),1.);
946  } else { //UserNtuples->fillp213(float(sector),float(0.),1.);
947  }
948  } // for
949 
950  if (nsumhit >= 2) { //UserNtuples->fillp212(vy,float(1.),1.);
951  } else { //UserNtuples->fillp212(vy,float(0.),1.);
952  }
953  */
954  } // MI or no MI or all - end
955  } // primary end
956 
957  if (verbosity > 0) {
958  edm::LogVerbatim("BscTest") << "BscTest: END OF Event " << (*evt)()->GetEventID();
959  }
960 }
961 
962 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
963 #include "FWCore/PluginManager/interface/ModuleDef.h"
964 
965 DEFINE_SIMWATCHER(BscTest);
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:134
BscTest::numofpart
int numofpart
Definition: BscTest.cc:137
BscTest::z2
double z2
Definition: BscTest.cc:142
DEFINE_SIMWATCHER
#define DEFINE_SIMWATCHER(type)
Definition: SimWatcherFactory.h:13
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
totemT2DQMSource_cfi.nbinsy
nbinsy
Definition: totemT2DQMSource_cfi.py:7
BscTest::z1
double z1
Definition: BscTest.cc:142
BscTest::detLevels
int detLevels(const G4VTouchable *) const
Definition: BscTest.cc:588
mps_fire.i
i
Definition: mps_fire.py:429
BscG4Hit::getUnitID
uint32_t getUnitID() const
Definition: BscG4Hit.h:56
BscTest::fOutputFile
std::string fOutputFile
Definition: BscTest.cc:152
BscAnalysisHistManager::getHisto
TH1F * getHisto(Int_t Number)
Definition: BscTest.cc:292
npart
double npart
Definition: HydjetWrapper.h:48
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
BscTest::lastpo
G4ThreeVector lastpo
Definition: BscTest.cc:139
BscAnalysisHistManager::fHistNamesArray
TObjArray * fHistNamesArray
Definition: BscTest.cc:93
HcalTBWriter_cfi.Number
Number
Definition: HcalTBWriter_cfi.py:19
BscTest::iev
int iev
Definition: BscTest.cc:120
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
BscAnalysisHistManager::fHistArray
TObjArray * fHistArray
Definition: BscTest.cc:92
BscTest::sumStepc
G4double sumStepc
Definition: BscTest.cc:135
BscTest::fRecreateFile
std::string fRecreateFile
Definition: BscTest.cc:153
pi
const Double_t pi
Definition: trackSplitPlot.h:36
BscTest::detName
G4String detName(const G4VTouchable *, int, int) const
Definition: BscTest.cc:597
BscTest::rinCalo
double rinCalo
Definition: BscTest.cc:130
BscNumberingScheme::unpackBscIndex
void unpackBscIndex(const unsigned int &idx)
Definition: BscNumberingScheme.cc:96
BscTest::verbosity
int verbosity
Definition: BscTest.cc:132
BscTest::z4
double z4
Definition: BscTest.cc:142
nano_mu_digi_cff.float
float
Definition: nano_mu_digi_cff.py:14
BscAnalysisHistManager::fTypeTitle
const char * fTypeTitle
Definition: BscTest.cc:91
BscTest::whichevent
int whichevent
Definition: BscTest.cc:148
BscAnalysisHistManager::histInit
void histInit(const char *name, const char *title, Int_t nbinsx, Axis_t xlow, Axis_t xup)
Definition: BscTest.cc:265
BscTest::detectorLevel
void detectorLevel(const G4VTouchable *, int &, int *, G4String *) const
Definition: BscTest.cc:607
BscTest::fDataLabel
std::string fDataLabel
Definition: BscTest.cc:151
BscTest::entot0
G4double entot0
Definition: BscTest.cc:122
BscTest::TheHistManager
BscAnalysisHistManager * TheHistManager
Definition: BscTest.cc:150
BscTest::z3
double z3
Definition: BscTest.cc:142
totemT2DQMSource_cfi.nbinsx
nbinsx
Definition: totemT2DQMSource_cfi.py:6
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
BscTest::update
void update(const BeginOfJob *run) override
This routine will be called when the appropriate signal arrives.
Definition: BscTest.cc:349
BscTest::BscTest
BscTest(const edm::ParameterSet &p)
Definition: BscTest.cc:159
btvMC_cff.vy
vy
Definition: btvMC_cff.py:26
BscTest::zinCalo
double zinCalo
Definition: BscTest.cc:130
cuy.ii
ii
Definition: cuy.py:589
BscAnalysisHistManager::getHisto2
TH2F * getHisto2(Int_t Number)
Definition: BscTest.cc:308
PbPb_ZMuSkimMuonDPG_cff.particleType
particleType
Definition: PbPb_ZMuSkimMuonDPG_cff.py:27
BscTest::sumEnerDeposit
G4double sumEnerDeposit
Definition: BscTest.cc:135
BscTest::sumStepl
G4double sumStepl
Definition: BscTest.cc:135
BscTest::bsceventntuple
TNtuple * bsceventntuple
Definition: BscTest.cc:146
BscTest::~BscTest
~BscTest() override
Definition: BscTest.cc:184
BscTest::lastTrackID
int lastTrackID
Definition: BscTest.cc:131
BscG4Hit::getEnergyLoss
float getEnergyLoss() const
Definition: BscG4Hit.h:69
BscG4Hit::getEntry
const G4ThreeVector & getEntry() const
Definition: BscG4Hit.h:29
btvMC_cff.vz
vz
Definition: btvMC_cff.py:27
BscTest::itrk
int itrk
Definition: BscTest.cc:121
BscG4Hit::getExitLocalP
const G4ThreeVector & getExitLocalP() const
Definition: BscG4Hit.h:35
TrackingDataMCValidation_Standalone_cff.nhits
nhits
Definition: TrackingDataMCValidation_Standalone_cff.py:40
btvMC_cff.vx
vx
Definition: btvMC_cff.py:25
BscTest::bsceventarray
Float_t bsceventarray[1]
Definition: BscTest.cc:145
step
step
Definition: StallMonitor.cc:83
BscTest::bscOutputFile
TFile bscOutputFile
Definition: BscTest.cc:147
BscG4HitCollection
G4THitsCollection< BscG4Hit > BscG4HitCollection
Definition: BscG4HitCollection.h:13
ntbsc_elements
ntbsc_elements
Definition: BscTest.cc:156
BscAnalysisHistManager::BscAnalysisHistManager
BscAnalysisHistManager(const TString &managername)
Definition: BscTest.cc:211
BscG4Hit::getEntryLocalP
const G4ThreeVector & getEntryLocalP() const
Definition: BscG4Hit.h:32
BscAnalysisHistManager::writeToFile
void writeToFile(const TString &fOutputFile, const TString &fRecreateFile)
Definition: BscTest.cc:251
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
BscTest::tracklength0
G4double tracklength0
Definition: BscTest.cc:122
BscAnalysisHistManager::~BscAnalysisHistManager
~BscAnalysisHistManager() override
Definition: BscTest.cc:225
Generated for CMSSW Reference Manual by doxygen 1.8.14