CMS 3D CMS Logo
HcalTB04Analysis.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HcalTestBeam
4 // Class : HcalTB04Analysis
5 //
6 // Implementation:
7 // Main analysis class for Hcal Test Beam 2004 Analysis
8 //
9 // Usage: A Simwatcher class and can be activated from Oscarproducer module
10 //
11 // Original Author:
12 // Created: Tue May 16 10:14:34 CEST 2006
13 //
14 
15 // system include files
16 #include <cmath>
17 #include <iomanip>
18 #include <iostream>
19 #include <memory>
20 #include <vector>
21 #include <string>
22 
23 // user include files
24 #include "SimG4Core/Notification/interface/Observer.h"
25 #include "SimG4Core/Notification/interface/BeginOfRun.h"
26 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
27 #include "SimG4Core/Notification/interface/EndOfEvent.h"
28 #include "SimG4Core/Watcher/interface/SimProducer.h"
29 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
30 
31 // to retreive hits
32 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
33 #include "SimDataFormats/CaloHit/interface/CaloHit.h"
34 #include "SimDataFormats/HcalTestBeam/interface/PHcalTB04Info.h"
35 #include "SimG4CMS/Calo/interface/ECalSD.h"
36 #include "SimG4CMS/Calo/interface/HCalSD.h"
37 #include "SimG4CMS/Calo/interface/HcalQie.h"
38 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
39 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
40 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
41 #include "DataFormats/Math/interface/Point3D.h"
42 
43 #include "FWCore/Framework/interface/Event.h"
44 #include "FWCore/Framework/interface/EventSetup.h"
45 #include "FWCore/Framework/interface/ESHandle.h"
46 #include "FWCore/Framework/interface/MakerMacros.h"
47 #include "FWCore/PluginManager/interface/ModuleDef.h"
48 #include "FWCore/MessageLogger/interface/MessageLogger.h"
49 
50 #include "HcalTBNumberingScheme.h"
51 #include "HcalTB04Histo.h"
52 #include "HcalTB04XtalNumberingScheme.h"
53 
54 #include "G4SDManager.hh"
55 #include "G4Step.hh"
56 #include "G4Track.hh"
57 #include "G4ThreeVector.hh"
58 #include "G4VProcess.hh"
59 #include "G4HCofThisEvent.hh"
60 #include "CLHEP/Random/RandGaussQ.h"
61 #include "CLHEP/Units/GlobalSystemOfUnits.h"
62 #include "CLHEP/Units/GlobalPhysicalConstants.h"
63 #include "Randomize.hh"
64 #include <cstdint>
65 
66 //#define EDM_ML_DEBUG
67 
68 namespace CLHEP {
69  class HepRandomEngine;
70 }
71 
72 class HcalTB04Analysis : public SimProducer,
73  public Observer<const BeginOfRun*>,
74  public Observer<const BeginOfEvent*>,
75  public Observer<const EndOfEvent*>,
76  public Observer<const G4Step*> {
77 public:
78  HcalTB04Analysis(const edm::ParameterSet& p);
79  ~HcalTB04Analysis() override;
80 
81  void produce(edm::Event&, const edm::EventSetup&) override;
82 
83 private:
84  HcalTB04Analysis(const HcalTB04Analysis&) = delete; // stop default
85  const HcalTB04Analysis& operator=(const HcalTB04Analysis&) = delete;
86 
87  void init();
88 
89  // observer methods
90  void update(const BeginOfRun* run) override;
91  void update(const BeginOfEvent* evt) override;
92  void update(const G4Step* step) override;
93  void update(const EndOfEvent* evt) override;
94 
95  //User methods
96  void fillBuffer(const EndOfEvent* evt);
97  void qieAnalysis(CLHEP::HepRandomEngine*);
98  void xtalAnalysis(CLHEP::HepRandomEngine*);
99  void finalAnalysis();
100  void fillEvent(PHcalTB04Info&);
101 
102  void clear();
103  int unitID(uint32_t id);
104  double scale(int det, int layer);
105  double timeOfFlight(int det, int layer, double eta);
106 
107 private:
108  HcalQie* myQie;
109  HcalTB04Histo* histo;
110 
111  // to read from parameter set
112  bool hcalOnly;
113  int mode, type;
114  double ecalNoise, beamOffset;
115  int iceta, icphi;
116  double scaleHB0, scaleHB16, scaleHO, scaleHE0;
117  std::vector<std::string> names;
118  G4RotationMatrix* beamline_RM;
119 
120  // Constants for the run
121  int count;
122  int nTower, nCrystal;
123  std::vector<int> idHcal, idXtal;
124  std::vector<uint32_t> idTower, idEcal;
125 
126  // Constants for the event
127  int nPrimary, particleType;
128  double pInit, etaInit, phiInit;
129  std::vector<CaloHit> ecalHitCache;
130  std::vector<CaloHit> hcalHitCache, hcalHitLayer;
131  std::vector<double> esimh, eqie, esime, enois;
132  std::vector<double> eseta, eqeta, esphi, eqphi, eslay, eqlay;
133  double etots, eecals, ehcals, etotq, eecalq, ehcalq;
134 
135  bool pvFound;
136  int evNum, pvType;
137  G4ThreeVector pvPosition, pvMomentum, pvUVW;
138  std::vector<int> secTrackID, secPartID;
139  std::vector<G4ThreeVector> secMomentum;
140  std::vector<double> secEkin;
141  std::vector<int> shortLivedSecondaries;
142 };
143 
144 //
145 // constructors and destructor
146 //
147 
148 HcalTB04Analysis::HcalTB04Analysis(const edm::ParameterSet& p) : myQie(nullptr), histo(nullptr) {
149  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB04Analysis");
150  hcalOnly = m_Anal.getParameter<bool>("HcalOnly");
151  mode = m_Anal.getParameter<int>("Mode");
152  type = m_Anal.getParameter<int>("Type");
153  ecalNoise = m_Anal.getParameter<double>("EcalNoise");
154  scaleHB0 = m_Anal.getParameter<double>("ScaleHB0");
155  scaleHB16 = m_Anal.getParameter<double>("ScaleHB16");
156  scaleHO = m_Anal.getParameter<double>("ScaleHO");
157  scaleHE0 = m_Anal.getParameter<double>("ScaleHE0");
158  names = m_Anal.getParameter<std::vector<std::string> >("Names");
159  beamOffset = -m_Anal.getParameter<double>("BeamPosition") * cm;
160  double fMinEta = m_Anal.getParameter<double>("MinEta");
161  double fMaxEta = m_Anal.getParameter<double>("MaxEta");
162  double fMinPhi = m_Anal.getParameter<double>("MinPhi");
163  double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
164  double beamEta = (fMaxEta + fMinEta) / 2.;
165  double beamPhi = (fMaxPhi + fMinPhi) / 2.;
166  double beamThet = 2 * atan(exp(-beamEta));
167  if (beamPhi < 0)
168  beamPhi += twopi;
169  iceta = (int)(beamEta / 0.087) + 1;
170  icphi = (int)(fabs(beamPhi) / 0.087) + 5;
171  if (icphi > 72)
172  icphi -= 73;
173 
174  produces<PHcalTB04Info>();
175 
176  beamline_RM = new G4RotationMatrix;
177  beamline_RM->rotateZ(-beamPhi);
178  beamline_RM->rotateY(-beamThet);
179 
180  edm::LogVerbatim("HcalTBSim")
181  << "HcalTB04:: Initialised as observer of BeginOf Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent with Parameter "
182  "values:\n \thcalOnly = "
183  << hcalOnly << "\tecalNoise = " << ecalNoise << "\n\tMode = " << mode << " (0: HB2 Standard; 1:HB2 Segmented)"
184  << "\tType = " << type << " (0: HB; 1 HE; 2 HB+HE)\n\tbeamOffset = " << beamOffset << "\ticeta = " << iceta
185  << "\ticphi = " << icphi << "\n\tbeamline_RM = " << *beamline_RM;
186 
187  init();
188 
189  myQie = new HcalQie(p);
190  histo = new HcalTB04Histo(m_Anal);
191 }
192 
193 HcalTB04Analysis::~HcalTB04Analysis() {
194  edm::LogVerbatim("HcalTBSim") << "\n --------> Total number of selected entries : " << count << "\nPointers:: QIE "
195  << myQie << " Histo " << histo;
196  if (myQie) {
197  delete myQie;
198  myQie = nullptr;
199  }
200  if (histo) {
201  delete histo;
202  histo = nullptr;
203  }
204 }
205 
206 //
207 // member functions
208 //
209 
210 void HcalTB04Analysis::produce(edm::Event& e, const edm::EventSetup&) {
211  std::unique_ptr<PHcalTB04Info> product(new PHcalTB04Info);
212  fillEvent(*product);
213  e.put(std::move(product));
214 }
215 
216 void HcalTB04Analysis::init() {
217  idTower = HcalTBNumberingScheme::getUnitIDs(type, mode);
218  nTower = idTower.size();
219  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save information from " << nTower << " HCal towers";
220  idHcal.reserve(nTower);
221  for (int i = 0; i < nTower; i++) {
222  int id = unitID(idTower[i]);
223  idHcal.push_back(id);
224 #ifdef EDM_ML_DEBUG
225  edm::LogVerbatim("HcalTBSim") << "\tTower[" << i << "] Original " << std::hex << idTower[i] << " Stored "
226  << idHcal[i] << std::dec;
227 #endif
228  }
229 
230  if (!hcalOnly) {
231  int det = 10;
232  uint32_t id1;
233  nCrystal = 0;
234  for (int lay = 1; lay < 8; lay++) {
235  for (int icr = 1; icr < 8; icr++) {
236  id1 = HcalTestNumbering::packHcalIndex(det, 0, 1, icr, lay, 1);
237  int id = unitID(id1);
238  idEcal.push_back(id1);
239  idXtal.push_back(id);
240  nCrystal++;
241  }
242  }
243  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save information from " << nCrystal << " ECal Crystals";
244 #ifdef EDM_ML_DEBUG
245  for (int i = 0; i < nCrystal; i++) {
246  edm::LogVerbatim("HcalTBSim") << "\tCrystal[" << i << "] Original " << std::hex << idEcal[i] << " Stored "
247  << idXtal[i] << std::dec;
248  }
249 #endif
250  }
251  // Profile vectors
252  eseta.reserve(5);
253  eqeta.reserve(5);
254  esphi.reserve(3);
255  eqphi.reserve(3);
256  eslay.reserve(20);
257  eqlay.reserve(20);
258  for (int i = 0; i < 5; i++) {
259  eseta.push_back(0.);
260  eqeta.push_back(0.);
261  }
262  for (int i = 0; i < 3; i++) {
263  esphi.push_back(0.);
264  eqphi.push_back(0.);
265  }
266  for (int i = 0; i < 20; i++) {
267  eslay.push_back(0.);
268  eqlay.push_back(0.);
269  }
270 
271  // counter
272  count = 0;
273  evNum = 0;
274  clear();
275 }
276 
277 void HcalTB04Analysis::update(const BeginOfRun* run) {
278  int irun = (*run)()->GetRunID();
279  edm::LogVerbatim("HcalTBSim") << " =====> Begin of Run = " << irun;
280 
281  G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
282  if (sd != nullptr) {
283  std::string sdname = names[0];
284  G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
285  if (aSD == nullptr) {
286  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
287  << " with name " << sdname << " in this "
288  << "Setup";
289  } else {
290  HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
291  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
292  << " in this Setup";
293  HcalNumberingScheme* org = new HcalTestNumberingScheme(false);
294  theCaloSD->setNumberingScheme(org);
295  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a new numbering scheme";
296  }
297  if (!hcalOnly) {
298  sdname = names[1];
299  aSD = sd->FindSensitiveDetector(sdname);
300  if (aSD == nullptr) {
301  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
302  << " with name " << sdname << " in this "
303  << "Setup";
304  } else {
305  ECalSD* theCaloSD = dynamic_cast<ECalSD*>(aSD);
306  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
307  << " in this Setup";
308  EcalNumberingScheme* org = new HcalTB04XtalNumberingScheme();
309  theCaloSD->setNumberingScheme(org);
310  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a new numbering scheme";
311  }
312  }
313  } else {
314  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Could "
315  << "not get SD Manager!";
316  }
317 }
318 
319 void HcalTB04Analysis::update(const BeginOfEvent* evt) {
320  evNum = (*evt)()->GetEventID();
321  clear();
322  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis: =====> Begin of event = " << evNum;
323 }
324 
325 void HcalTB04Analysis::update(const G4Step* aStep) {
326  if (aStep != nullptr) {
327  //Get Step properties
328  G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
329  G4ThreeVector thePostStepPoint;
330 
331  // Get Tracks properties
332  G4Track* aTrack = aStep->GetTrack();
333  int trackID = aTrack->GetTrackID();
334  int parentID = aTrack->GetParentID();
335  const G4ThreeVector& position = aTrack->GetPosition();
336  G4ThreeVector momentum = aTrack->GetMomentum();
337  G4String partType = aTrack->GetDefinition()->GetParticleType();
338  G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
339  int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
340 #ifdef EDM_ML_DEBUG
341  bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
342 #endif
343  double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
344  double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
345 
346  if (!pvFound) { //search for v1
347  double stepDeltaEnergy = aStep->GetDeltaEnergy();
348  double kinEnergy = aTrack->GetKineticEnergy();
349 
350  // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
351  if (trackID == 1 && parentID == 0 && ((kinEnergy == 0.) || (fabs(stepDeltaEnergy / kinEnergy) > 0.1))) {
352  pvType = -1;
353  if (kinEnergy == 0.) {
354  pvType = 0;
355  } else {
356  if (fabs(stepDeltaEnergy / kinEnergy) > 0.1)
357  pvType = 1;
358  }
359  pvFound = true;
360  pvPosition = position;
361  pvMomentum = momentum;
362  // Rotated coord.system:
363  pvUVW = (*beamline_RM) * (pvPosition);
364 
365  //Volume name requires some checks:
366  G4String thePostPVname = "NoName";
367  G4StepPoint* thePostPoint = aStep->GetPostStepPoint();
368  if (thePostPoint) {
369  thePostStepPoint = thePostPoint->GetPosition();
370  G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
371  if (thePostPV)
372  thePostPVname = thePostPV->GetName();
373  }
374 #ifdef EDM_ML_DEBUG
375  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: V1 found at: " << thePostStepPoint
376  << " G4VPhysicalVolume: " << thePostPVname;
377  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::fill_v1Pos: Primary Track momentum: " << pvMomentum
378  << " psoition " << pvPosition << " u/v/w " << pvUVW;
379 #endif
380  }
381  } else {
382  // watch for secondaries originating @v1, including the surviving primary
383  if ((trackID != 1 && parentID == 1 && (aTrack->GetCurrentStepNumber() == 1) && (thePreStepPoint == pvPosition)) ||
384  (trackID == 1 && thePreStepPoint == pvPosition)) {
385 #ifdef EDM_ML_DEBUG
386  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::A secondary... PDG:" << partPDGEncoding
387  << " TrackID:" << trackID << " ParentID:" << parentID
388  << " stable: " << isPDGStable << " Tau: " << pDGlifetime
389  << " cTauGamma=" << c_light * pDGlifetime * gammaFactor * 1000.
390  << "um GammaFactor: " << gammaFactor;
391 #endif
392  secTrackID.push_back(trackID);
393  secPartID.push_back(partPDGEncoding);
394  secMomentum.push_back(momentum);
395  secEkin.push_back(aTrack->GetKineticEnergy());
396 
397  // Check for short-lived secondaries: cTauGamma<100um
398  double ctaugamma_um = CLHEP::c_light * pDGlifetime * gammaFactor * 1000.;
399  if ((ctaugamma_um > 0.) && (ctaugamma_um < 100.)) { //short-lived secondary
400  shortLivedSecondaries.push_back(trackID);
401  } else { //normal secondary - enter into the V1-calorimetric tree
402  // histos->fill_v1cSec (aTrack);
403  }
404  }
405  // Also watch for tertiary particles coming from
406  // short-lived secondaries from V1
407  if (aTrack->GetCurrentStepNumber() == 1) {
408  if (!shortLivedSecondaries.empty()) {
409  int pid = parentID;
410  std::vector<int>::iterator pos1 = shortLivedSecondaries.begin();
411  std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
412  std::vector<int>::iterator pos;
413  for (pos = pos1; pos != pos2; pos++) {
414  if (*pos == pid) { //ParentID is on the list of short-lived
415  // secondary
416 #ifdef EDM_ML_DEBUG
417  edm::LogVerbatim("HcalTBSim")
418  << "HcalTB04Analysis:: A tertiary... PDG:" << partPDGEncoding << " TrackID:" << trackID
419  << " ParentID:" << parentID << " stable: " << isPDGStable << " Tau: " << pDGlifetime
420  << " cTauGamma=" << c_light * pDGlifetime * gammaFactor * 1000. << "um GammaFactor: " << gammaFactor;
421 #endif
422  }
423  }
424  }
425  }
426  }
427  }
428 }
429 
430 void HcalTB04Analysis::update(const EndOfEvent* evt) {
431  count++;
432 
433  //fill the buffer
434 #ifdef EDM_ML_DEBUG
435  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Fill event " << (*evt)()->GetEventID();
436 #endif
437  fillBuffer(evt);
438 
439  //QIE analysis
440 #ifdef EDM_ML_DEBUG
441  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Do QIE analysis with " << hcalHitCache.size() << " hits";
442 #endif
443  CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
444  qieAnalysis(engine);
445 
446  //Energy in Crystal Matrix
447  if (!hcalOnly) {
448 #ifdef EDM_ML_DEBUG
449  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Do Xtal analysis with " << ecalHitCache.size() << " hits";
450 #endif
451  xtalAnalysis(engine);
452  }
453 
454  //Final Analysis
455 #ifdef EDM_ML_DEBUG
456  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Final analysis";
457 #endif
458  finalAnalysis();
459 
460  int iEvt = (*evt)()->GetEventID();
461  if (iEvt < 10)
462  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
463  else if ((iEvt < 100) && (iEvt % 10 == 0))
464  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
465  else if ((iEvt < 1000) && (iEvt % 100 == 0))
466  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
467  else if ((iEvt < 10000) && (iEvt % 1000 == 0))
468  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
469 }
470 
471 void HcalTB04Analysis::fillBuffer(const EndOfEvent* evt) {
472  std::vector<CaloHit> hhits, hhitl;
473  int idHC, j;
474  CaloG4HitCollection* theHC;
475  std::map<int, float, std::less<int> > primaries;
476  double etot1 = 0, etot2 = 0;
477 
478  // Look for the Hit Collection of HCal
479  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
480  std::string sdName = names[0];
481  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
482  theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
483 #ifdef EDM_ML_DEBUG
484  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
485  << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
486 #endif
487  int thehc_entries = theHC->entries();
488  if (idHC >= 0 && theHC != nullptr) {
489  hhits.reserve(theHC->entries());
490  hhitl.reserve(theHC->entries());
491  for (j = 0; j < thehc_entries; j++) {
492  CaloG4Hit* aHit = (*theHC)[j];
493  double e = aHit->getEnergyDeposit() / CLHEP::GeV;
494  double time = aHit->getTimeSlice();
495  math::XYZPoint pos = aHit->getEntry();
496  unsigned int id = aHit->getUnitID();
497  double theta = pos.theta();
498  double eta = -log(tan(theta * 0.5));
499  double phi = pos.phi();
500  int det, z, group, ieta, iphi, layer;
501  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
502  double jitter = time - timeOfFlight(det, layer, eta);
503  if (jitter < 0)
504  jitter = 0;
505  if (e < 0 || e > 1.)
506  e = 0;
507  double escl = e * scale(det, layer);
508  unsigned int idx = HcalTBNumberingScheme::getUnitID(id, mode);
509  CaloHit hit(det, layer, escl, eta, phi, jitter, idx);
510  hhits.push_back(hit);
511  CaloHit hitl(det, layer, escl, eta, phi, jitter, id);
512  hhitl.push_back(hitl);
513  primaries[aHit->getTrackID()] += e;
514  etot1 += escl;
515 #ifdef EDM_ML_DEBUG
516  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j << " ID 0x" << std::hex << id << " 0x"
517  << idx << std::dec << " time " << std::setw(6) << time << " " << std::setw(6)
518  << jitter << " theta " << std::setw(8) << theta << " eta " << std::setw(8) << eta
519  << " phi " << std::setw(8) << phi << " e " << std::setw(8) << e << " "
520  << std::setw(8) << escl;
521 #endif
522  }
523  }
524 
525  // Add hits in the same channel within same time slice
526  std::vector<CaloHit>::iterator itr;
527  int nHit = hhits.size();
528  std::vector<CaloHit*> hits(nHit);
529  for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
530  hits[j] = &hhits[j];
531  }
532  sort(hits.begin(), hits.end(), CaloHitIdMore());
533  std::vector<CaloHit*>::iterator k1, k2;
534  int nhit = 0;
535  for (k1 = hits.begin(); k1 != hits.end(); k1++) {
536  int det = (**k1).det();
537  int layer = (**k1).layer();
538  double ehit = (**k1).e();
539  double eta = (**k1).eta();
540  double phi = (**k1).phi();
541  double jitter = (**k1).t();
542  uint32_t unitID = (**k1).id();
543  int jump = 0;
544  for (k2 = k1 + 1; k2 != hits.end() && fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
545  ehit += (**k2).e();
546  jump++;
547  }
548  nhit++;
549  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
550  hcalHitCache.push_back(hit);
551  etot2 += ehit;
552  k1 += jump;
553 #ifdef EDM_ML_DEBUG
554  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit << " ID 0x" << std::hex << unitID
555  << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
556  << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
557 #endif
558  }
559 #ifdef EDM_ML_DEBUG
560  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits from " << nHit
561  << " input hits E(Hcal) " << etot1 << " " << etot2;
562 #endif
563  //Repeat for Hit in each layer (hhits and hhitl sizes are the same)
564  for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
565  hits[j] = &hhitl[j];
566  }
567  sort(hits.begin(), hits.end(), CaloHitIdMore());
568  int nhitl = 0;
569  double etotl = 0;
570  for (k1 = hits.begin(); k1 != hits.end(); k1++) {
571  int det = (**k1).det();
572  int layer = (**k1).layer();
573  double ehit = (**k1).e();
574  double eta = (**k1).eta();
575  double phi = (**k1).phi();
576  double jitter = (**k1).t();
577  uint32_t unitID = (**k1).id();
578  int jump = 0;
579  for (k2 = k1 + 1; k2 != hits.end() && fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
580  ehit += (**k2).e();
581  jump++;
582  }
583  nhitl++;
584  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
585  hcalHitLayer.push_back(hit);
586  etotl += ehit;
587  k1 += jump;
588 #ifdef EDM_ML_DEBUG
589  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl << " ID 0x" << std::hex << unitID
590  << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
591  << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
592 #endif
593  }
594 #ifdef EDM_ML_DEBUG
595  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal hits from " << nHit
596  << " input hits E(Hcal) " << etot1 << " " << etotl;
597 #endif
598  // Look for the Hit Collection of ECal
599  std::vector<CaloHit> ehits;
600  sdName = names[1];
601  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
602  theHC = (CaloG4HitCollection*)allHC->GetHC(idHC);
603  etot1 = etot2 = 0;
604 #ifdef EDM_ML_DEBUG
605  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName << " of ID " << idHC
606  << " is obtained at " << theHC << " with " << theHC->entries() << " entries";
607 #endif
608  if (idHC >= 0 && theHC != nullptr) {
609  thehc_entries = theHC->entries();
610  ehits.reserve(theHC->entries());
611  for (j = 0; j < thehc_entries; j++) {
612  CaloG4Hit* aHit = (*theHC)[j];
613  double e = aHit->getEnergyDeposit() / CLHEP::GeV;
614  double time = aHit->getTimeSlice();
615  if (e < 0 || e > 100000.)
616  e = 0;
617  if (e > 0) {
618  math::XYZPoint pos = aHit->getEntry();
619  unsigned int id = aHit->getUnitID();
620  double theta = pos.theta();
621  double eta = -log(tan(theta * 0.5));
622  double phi = pos.phi();
623  int det, z, group, ieta, iphi, layer;
624  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
625  CaloHit hit(det, 0, e, eta, phi, time, id);
626  ehits.push_back(hit);
627  primaries[aHit->getTrackID()] += e;
628  etot1 += e;
629 #ifdef EDM_ML_DEBUG
630  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j << " ID 0x" << std::hex << id
631  << std::dec << " time " << std::setw(6) << time << " theta " << std::setw(8)
632  << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
633  << " e " << std::setw(8) << e;
634 #endif
635  }
636  }
637  }
638 
639  // Add hits in the same channel within same time slice
640  nHit = ehits.size();
641  std::vector<CaloHit*> hite(nHit);
642  for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
643  hite[j] = &ehits[j];
644  }
645  sort(hite.begin(), hite.end(), CaloHitIdMore());
646  nhit = 0;
647  for (k1 = hite.begin(); k1 != hite.end(); k1++) {
648  int det = (**k1).det();
649  int layer = (**k1).layer();
650  double ehit = (**k1).e();
651  double eta = (**k1).eta();
652  double phi = (**k1).phi();
653  double jitter = (**k1).t();
654  uint32_t unitID = (**k1).id();
655  int jump = 0;
656  for (k2 = k1 + 1; k2 != hite.end() && fabs(jitter - (**k2).t()) < 1 && unitID == (**k2).id(); k2++) {
657  ehit += (**k2).e();
658  jump++;
659  }
660  nhit++;
661  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
662  ecalHitCache.push_back(hit);
663  etot2 += ehit;
664  k1 += jump;
665 #ifdef EDM_ML_DEBUG
666  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit << " ID 0x" << std::hex << unitID
667  << std::dec << " time " << std::setw(6) << jitter << " eta " << std::setw(8) << eta
668  << " phi " << std::setw(8) << phi << " e " << std::setw(8) << ehit;
669 #endif
670  }
671 #ifdef EDM_ML_DEBUG
672  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits from " << nHit
673  << " input hits E(Ecal) " << etot1 << " " << etot2;
674 #endif
675  // Find Primary info:
676  nPrimary = (int)(primaries.size());
677  int trackID = 0;
678  G4PrimaryParticle* thePrim = nullptr;
679  int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
680 #ifdef EDM_ML_DEBUG
681  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex << " verteices";
682 #endif
683  if (nvertex <= 0)
684  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no vertex found for event " << evNum;
685  for (int i = 0; i < nvertex; i++) {
686  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
687  if (avertex == nullptr) {
688  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer to vertex = 0 for event " << evNum;
689  } else {
690  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " " << avertex->GetPosition();
691  int npart = avertex->GetNumberOfParticle();
692  if (npart == 0)
693  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: no primary!";
694  if (thePrim == nullptr)
695  thePrim = avertex->GetPrimary(trackID);
696  }
697  }
698 
699  if (thePrim != nullptr) {
700  double px = thePrim->GetPx();
701  double py = thePrim->GetPy();
702  double pz = thePrim->GetPz();
703  double p = std::sqrt(pow(px, 2.) + pow(py, 2.) + pow(pz, 2.));
704  pInit = p / CLHEP::GeV;
705  if (p == 0)
706  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: primary has p=0 ";
707  else {
708  double costheta = pz / p;
709  double theta = acos(std::min(std::max(costheta, -1.), 1.));
710  etaInit = -log(tan(theta / 2));
711  if (px != 0 || py != 0)
712  phiInit = atan2(py, px);
713  }
714  particleType = thePrim->GetPDGcode();
715  } else
716  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could not find primary";
717 }
718 
719 void HcalTB04Analysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
720  int hittot = hcalHitCache.size();
721  if (hittot <= 0)
722  hittot = 1;
723  std::vector<CaloHit> hits(hittot);
724  std::vector<int> todo(nTower, 0);
725 
726 #ifdef EDM_ML_DEBUG
727  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size " << hits.size() << " " << todo.size() << " "
728  << idTower.size() << " " << esimh.size() << " " << eqie.size();
729 #endif
730  // Loop over all HCal hits
731  for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
732  CaloHit hit = hcalHitCache[k1];
733  uint32_t id = hit.id();
734  int nhit = 0;
735  double esim = hit.e();
736  hits[nhit] = hit;
737  for (unsigned int k2 = k1 + 1; k2 < hcalHitCache.size(); k2++) {
738  hit = hcalHitCache[k2];
739  if (hit.id() == id) {
740  nhit++;
741  hits[nhit] = hit;
742  esim += hit.e();
743  }
744  }
745  k1 += nhit;
746  nhit++;
747  std::vector<int> cd = myQie->getCode(nhit, hits, engine);
748  double eq = myQie->getEnergy(cd);
749 #ifdef EDM_ML_DEBUG
750  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id << std::dec << " registers " << esim
751  << " energy from " << nhit << " hits starting with hit # " << k1
752  << " energy with noise " << eq;
753 #endif
754  for (int k2 = 0; k2 < nTower; k2++) {
755  if (id == idTower[k2]) {
756  todo[k2] = 1;
757  esimh[k2] = esim;
758  eqie[k2] = eq;
759  }
760  }
761  }
762 
763  // Towers with no hit
764  for (int k2 = 0; k2 < nTower; k2++) {
765  if (todo[k2] == 0) {
766  std::vector<int> cd = myQie->getCode(0, hits, engine);
767  double eq = myQie->getEnergy(cd);
768  esimh[k2] = 0;
769  eqie[k2] = eq;
770 #ifdef EDM_ML_DEBUG
771  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << idTower[k2] << std::dec
772  << " registers " << esimh[k2] << " energy from hits and energy after QIE analysis "
773  << eqie[k2];
774 #endif
775  }
776  }
777 }
778 
779 void HcalTB04Analysis::xtalAnalysis(CLHEP::HepRandomEngine* engine) {
780  CLHEP::RandGaussQ randGauss(*engine);
781 
782  // Crystal Data
783  std::vector<int> iok(nCrystal, 0);
784 #ifdef EDM_ML_DEBUG
785  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " << iok.size() << " " << idEcal.size() << " "
786  << esime.size() << " " << enois.size();
787 #endif
788  for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
789  uint32_t id = ecalHitCache[k1].id();
790  int nhit = 0;
791  double esim = ecalHitCache[k1].e();
792  for (unsigned int k2 = k1 + 1; k2 < ecalHitCache.size(); k2++) {
793  if (ecalHitCache[k2].id() == id) {
794  nhit++;
795  esim += ecalHitCache[k2].e();
796  }
797  }
798  k1 += nhit;
799  nhit++;
800  double eq = esim + randGauss.fire(0., ecalNoise);
801 #ifdef EDM_ML_DEBUG
802  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id << std::dec << " registers " << esim
803  << " energy from " << nhit << " hits starting with hit # " << k1
804  << " energy with noise " << eq;
805 #endif
806  for (int k2 = 0; k2 < nCrystal; k2++) {
807  if (id == idEcal[k2]) {
808  iok[k2] = 1;
809  esime[k2] = esim;
810  enois[k2] = eq;
811  }
812  }
813  }
814 
815  // Crystals with no hit
816  for (int k2 = 0; k2 < nCrystal; k2++) {
817  if (iok[k2] == 0) {
818  esime[k2] = 0;
819  enois[k2] = randGauss.fire(0., ecalNoise);
820 #ifdef EDM_ML_DEBUG
821  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << idEcal[k2] << std::dec
822  << " registers " << esime[k2] << " energy from hits and energy from noise "
823  << enois[k2];
824 #endif
825  }
826  }
827 }
828 
829 void HcalTB04Analysis::finalAnalysis() {
830  //Beam Information
831  histo->fillPrimary(pInit, etaInit, phiInit);
832 
833  // Total Energy
834  eecals = ehcals = eecalq = ehcalq = 0.;
835  for (int i = 0; i < nTower; i++) {
836  ehcals += esimh[i];
837  ehcalq += eqie[i];
838  }
839  for (int i = 0; i < nCrystal; i++) {
840  eecals += esime[i];
841  eecalq += enois[i];
842  }
843  etots = eecals + ehcals;
844  etotq = eecalq + ehcalq;
845 #ifdef EDM_ML_DEBUG
846  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level (Total) " << etots << " (ECal) "
847  << eecals << " (HCal) " << ehcals
848  << "\nHcalTB04Analysis:: Energy deposit at Qie Level (Total) " << etotq << " (ECal) "
849  << eecalq << " (HCal) " << ehcalq;
850 #endif
851  histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
852 
853  // Lateral Profile
854  for (int i = 0; i < 5; i++) {
855  eseta[i] = 0.;
856  eqeta[i] = 0.;
857  }
858  for (int i = 0; i < 3; i++) {
859  esphi[i] = 0.;
860  eqphi[i] = 0.;
861  }
862  double e1 = 0, e2 = 0;
863  unsigned int id;
864  for (int i = 0; i < nTower; i++) {
865  int det, z, group, ieta, iphi, layer;
866  id = idTower[i];
867  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
868  iphi -= (icphi - 1);
869  if (icphi > 4) {
870  if (ieta == 0)
871  ieta = 2;
872  else
873  ieta = -1;
874  } else {
875  ieta = ieta - iceta + 2;
876  }
877  if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
878  eseta[ieta] += esimh[i];
879  esphi[iphi] += esimh[i];
880  e1 += esimh[i];
881  eqeta[ieta] += eqie[i];
882  eqphi[iphi] += eqie[i];
883  e2 += eqie[i];
884  }
885  }
886  for (int i = 0; i < 3; i++) {
887  if (e1 > 0)
888  esphi[i] /= e1;
889  if (e2 > 0)
890  eqphi[i] /= e2;
891  }
892  for (int i = 0; i < 5; i++) {
893  if (e1 > 0)
894  eseta[i] /= e1;
895  if (e2 > 0)
896  eqeta[i] /= e2;
897  }
898 #ifdef EDM_ML_DEBUG
899  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and Phi (Sim/Qie)";
900  for (int i = 0; i < 5; i++)
901  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = " << eseta[i] << " Qie = " << eqeta[i]
902  << " Phi Sim = " << esphi[i] << " Qie = " << eqphi[i];
903 #endif
904  histo->fillTrnsProf(eseta, eqeta, esphi, eqphi);
905 
906  // Longitudianl profile
907  for (int i = 0; i < 20; i++) {
908  eslay[i] = 0.;
909  eqlay[i] = 0.;
910  }
911  e1 = 0;
912  e2 = 0;
913  for (int i = 0; i < nTower; i++) {
914  int det, z, group, ieta, iphi, layer;
915  id = idTower[i];
916  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, layer);
917  iphi -= (icphi - 1);
918  layer -= 1;
919  if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
920  eslay[layer] += esimh[i];
921  e1 += esimh[i];
922  eqlay[layer] += eqie[i];
923  e2 += eqie[i];
924  }
925  }
926  for (int i = 0; i < 20; i++) {
927  if (e1 > 0)
928  eslay[i] /= e1;
929  if (e2 > 0)
930  eqlay[i] /= e2;
931  }
932 #ifdef EDM_ML_DEBUG
933  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
934  for (int i = 0; i < 20; i++)
935  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = " << eslay[i] << " Qie = " << eqlay[i];
936 #endif
937  histo->fillLongProf(eslay, eqlay);
938 }
939 
940 void HcalTB04Analysis::fillEvent(PHcalTB04Info& product) {
941  //Setup the ID's
942  product.setIDs(idHcal, idXtal);
943 
944  //Beam Information
945  product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
946 
947  //Energy deposits in the crystals and towers
948  product.setEdepHcal(esimh, eqie);
949  product.setEdepHcal(esime, enois);
950 
951  // Total Energy
952  product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
953 
954  // Lateral Profile
955  product.setTrnsProf(eseta, eqeta, esphi, eqphi);
956 
957  // Longitudianl profile
958  product.setLongProf(eslay, eqlay);
959 
960  //Save Hits
961  int i, nhit = 0;
962  std::vector<CaloHit>::iterator itr;
963  for (i = 0, itr = ecalHitCache.begin(); itr != ecalHitCache.end(); i++, itr++) {
964  uint32_t id = itr->id();
965  int det, z, group, ieta, iphi, lay;
966  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
967  product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
968  nhit++;
969 #ifdef EDM_ML_DEBUG
970  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
971  << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
972  << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
973  << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
974 #endif
975  }
976  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from Crystals";
977  int hit = nhit;
978  nhit = 0;
979 
980  for (i = hit, itr = hcalHitCache.begin(); itr != hcalHitCache.end(); i++, itr++) {
981  uint32_t id = itr->id();
982  int det, z, group, ieta, iphi, lay;
983  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
984  product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
985  nhit++;
986 #ifdef EDM_ML_DEBUG
987  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) << i + 1 << " ID 0x" << std::hex
988  << group << std::dec << " " << std::setw(2) << det << " " << std::setw(2) << lay
989  << " " << std::setw(1) << z << " " << std::setw(3) << ieta << " " << std::setw(3)
990  << iphi << " T " << std::setw(6) << itr->t() << " E " << std::setw(6) << itr->e();
991 #endif
992  }
993  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit << " hits from HCal";
994 
995  //Vertex associated quantities
996  product.setVtxPrim(evNum,
997  pvType,
998  pvPosition.x(),
999  pvPosition.y(),
1000  pvPosition.z(),
1001  pvUVW.x(),
1002  pvUVW.y(),
1003  pvUVW.z(),
1004  pvMomentum.x(),
1005  pvMomentum.y(),
1006  pvMomentum.z());
1007  for (unsigned int i = 0; i < secTrackID.size(); i++) {
1008  product.setVtxSec(
1009  secTrackID[i], secPartID[i], secMomentum[i].x(), secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
1010  }
1011 }
1012 
1013 void HcalTB04Analysis::clear() {
1014  pvFound = false;
1015  pvType = -2;
1016  pvPosition = G4ThreeVector();
1017  pvMomentum = G4ThreeVector();
1018  pvUVW = G4ThreeVector();
1019  secTrackID.clear();
1020  secPartID.clear();
1021  secMomentum.clear();
1022  secEkin.clear();
1023  shortLivedSecondaries.clear();
1024 
1025  ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
1026  hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
1027  hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
1028  nPrimary = particleType = 0;
1029  pInit = etaInit = phiInit = 0;
1030 
1031  esimh.clear();
1032  eqie.clear();
1033  esimh.reserve(nTower);
1034  eqie.reserve(nTower);
1035  for (int i = 0; i < nTower; i++) {
1036  esimh.push_back(0.);
1037  eqie.push_back(0.);
1038  }
1039  esime.clear();
1040  enois.clear();
1041  esime.reserve(nCrystal);
1042  enois.reserve(nCrystal);
1043  for (int i = 0; i < nCrystal; i++) {
1044  esime.push_back(0.);
1045  enois.push_back(0.);
1046  }
1047 }
1048 
1049 int HcalTB04Analysis::unitID(uint32_t id) {
1050  int det, z, group, ieta, iphi, lay;
1051  HcalTestNumbering::unpackHcalIndex(id, det, z, group, ieta, iphi, lay);
1052  group = (det & 15) << 20;
1053  group += ((lay - 1) & 31) << 15;
1054  group += (z & 1) << 14;
1055  group += (ieta & 127) << 7;
1056  group += (iphi & 127);
1057  return group;
1058 }
1059 
1060 double HcalTB04Analysis::scale(int det, int layer) {
1061  double tmp = 1.;
1062  if (det == static_cast<int>(HcalBarrel)) {
1063  if (layer == 1)
1064  tmp = scaleHB0;
1065  else if (layer == 17)
1066  tmp = scaleHB16;
1067  else if (layer > 17)
1068  tmp = scaleHO;
1069  } else {
1070  if (layer <= 2)
1071  tmp = scaleHE0;
1072  }
1073  return tmp;
1074 }
1075 
1076 double HcalTB04Analysis::timeOfFlight(int det, int layer, double eta) {
1077  double theta = 2.0 * atan(exp(-eta));
1078  double dist = beamOffset;
1079  if (det == static_cast<int>(HcalBarrel)) {
1080  const double rLay[19] = {1836.0,
1081  1902.0,
1082  1962.0,
1083  2022.0,
1084  2082.0,
1085  2142.0,
1086  2202.0,
1087  2262.0,
1088  2322.0,
1089  2382.0,
1090  2448.0,
1091  2514.0,
1092  2580.0,
1093  2646.0,
1094  2712.0,
1095  2776.0,
1096  2862.5,
1097  3847.0,
1098  4052.0};
1099  if (layer > 0 && layer <= 19)
1100  dist += rLay[layer - 1] * mm / sin(theta);
1101  } else {
1102  const double zLay[19] = {4034.0,
1103  4032.0,
1104  4123.0,
1105  4210.0,
1106  4297.0,
1107  4384.0,
1108  4471.0,
1109  4558.0,
1110  4645.0,
1111  4732.0,
1112  4819.0,
1113  4906.0,
1114  4993.0,
1115  5080.0,
1116  5167.0,
1117  5254.0,
1118  5341.0,
1119  5428.0,
1120  5515.0};
1121  if (layer > 0 && layer <= 19)
1122  dist += zLay[layer - 1] * mm / cos(theta);
1123  }
1124 
1125  double tmp = dist / c_light / ns;
1126 #ifdef EDM_ML_DEBUG
1127  edm::LogVerbatim("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
1128  << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
1129 #endif
1130  return tmp;
1131 }
1132 
1133 DEFINE_SIMWATCHER(HcalTB04Analysis);
HcalTB04Analysis::esime
std::vector< double > esime
Definition: HcalTB04Analysis.cc:131
HcalTB04Analysis::particleType
int particleType
Definition: HcalTB04Analysis.cc:127
DDAxes::y
HcalTB04Analysis::xtalAnalysis
void xtalAnalysis(CLHEP::HepRandomEngine *)
Definition: HcalTB04Analysis.cc:779
CaloG4Hit::getTrackID
int getTrackID() const
Definition: CaloG4Hit.h:64
mps_fire.i
i
Definition: mps_fire.py:428
HcalTB04Analysis::idXtal
std::vector< int > idXtal
Definition: HcalTB04Analysis.cc:123
Observer
Definition: Observer.h:23
HcalTB04Analysis::iceta
int iceta
Definition: HcalTB04Analysis.cc:115
HcalTB04Analysis::timeOfFlight
double timeOfFlight(int det, int layer, double eta)
Definition: HcalTB04Analysis.cc:1076
MessageLogger.h
DQMOfflineHeavyIons_cff.todo
todo
Definition: DQMOfflineHeavyIons_cff.py:112
HcalTB04Analysis::beamline_RM
G4RotationMatrix * beamline_RM
Definition: HcalTB04Analysis.cc:118
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:92
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
ESHandle.h
step
step
Definition: StallMonitor.cc:94
HcalTB04Analysis::pvUVW
G4ThreeVector pvUVW
Definition: HcalTB04Analysis.cc:137
HcalTB04Analysis::hcalHitCache
std::vector< CaloHit > hcalHitCache
Definition: HcalTB04Analysis.cc:130
min
T min(T a, T b)
Definition: MathUtil.h:58
CaloG4HitCollection.h
HCalSD
Definition: HCalSD.h:38
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:65
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
SimProducer.h
HcalTB04Analysis::ecalHitCache
std::vector< CaloHit > ecalHitCache
Definition: HcalTB04Analysis.cc:129
HcalTB04Histo::fillTrnsProf
void fillTrnsProf(const std::vector< double > &es1, const std::vector< double > &eq1, const std::vector< double > &es2, const std::vector< double > &eq2)
Definition: HcalTB04Histo.cc:98
HcalTB04Analysis::histo
HcalTB04Histo * histo
Definition: HcalTB04Analysis.cc:109
HcalTB04Analysis::beamOffset
double beamOffset
Definition: HcalTB04Analysis.cc:114
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HcalQie.h
globals_cff.id1
id1
Definition: globals_cff.py:33
pos
Definition: PixelAliasList.h:18
HcalTB04Analysis::ehcals
double ehcals
Definition: HcalTB04Analysis.cc:133
HcalTB04Analysis::eecalq
double eecalq
Definition: HcalTB04Analysis.cc:133
PHcalTB04Info::setIDs
void setIDs(const std::vector< int > &, const std::vector< int > &)
Definition: PHcalTB04Info.cc:69
HcalTB04Analysis::ehcalq
double ehcalq
Definition: HcalTB04Analysis.cc:133
timingPdfMaker.histo
histo
Definition: timingPdfMaker.py:279
HcalTB04Analysis::idTower
std::vector< uint32_t > idTower
Definition: HcalTB04Analysis.cc:124
PHcalTB04Info::setPrimary
void setPrimary(int primary, int id, double energy, double eta, double phi)
Definition: PHcalTB04Info.cc:95
HcalTB04Analysis::pvType
int pvType
Definition: HcalTB04Analysis.cc:136
HcalTB04Analysis::eqlay
std::vector< double > eqlay
Definition: HcalTB04Analysis.cc:132
HcalTB04Analysis::scale
double scale(int det, int layer)
Definition: HcalTB04Analysis.cc:1060
EndOfEvent.h
HcalTB04Analysis::unitID
int unitID(uint32_t id)
Definition: HcalTB04Analysis.cc:1049
DDAxes::x
PHcalTB04Info::saveHit
void saveHit(int det, int lay, int eta, int phi, double e, double t)
Definition: PHcalTB04Info.cc:185
HcalTB04Analysis::HcalTB04Analysis
HcalTB04Analysis(const edm::ParameterSet &p)
Definition: HcalTB04Analysis.cc:148
HcalTB04Histo
Definition: HcalTB04Histo.h:31
HcalBarrel
Definition: HcalAssistant.h:33
HcalTB04Analysis::esphi
std::vector< double > esphi
Definition: HcalTB04Analysis.cc:132
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
CaloG4Hit::getTimeSlice
double getTimeSlice() const
Definition: CaloG4Hit.h:66
PHcalTB04Info::setEdep
void setEdep(double simtot, double sime, double simh, double digtot, double dige, double digh)
Definition: PHcalTB04Info.cc:105
PHcalTB04Info::setTrnsProf
void setTrnsProf(const std::vector< double > &es1, const std::vector< double > &eq1, const std::vector< double > &es2, const std::vector< double > &eq2)
Definition: PHcalTB04Info.cc:149
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
HcalTB04Analysis::pvMomentum
G4ThreeVector pvMomentum
Definition: HcalTB04Analysis.cc:137
npart
double npart
Definition: HydjetWrapper.h:46
Observer.h
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
LEDCalibrationChannels.iphi
iphi
Definition: LEDCalibrationChannels.py:64
HcalTB04Analysis::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: HcalTB04Analysis.cc:210
DEFINE_SIMWATCHER
#define DEFINE_SIMWATCHER(type)
Definition: SimWatcherFactory.h:13
HcalTB04Analysis::operator=
const HcalTB04Analysis & operator=(const HcalTB04Analysis &)=delete
HcalTB04Analysis::nTower
int nTower
Definition: HcalTB04Analysis.cc:122
HcalTB04Histo.h
HcalTB04Analysis
Definition: HcalTB04Analysis.cc:72
EcalNumberingScheme
Definition: EcalNumberingScheme.h:14
MakerMacros.h
HcalTB04Analysis::~HcalTB04Analysis
~HcalTB04Analysis() override
Definition: HcalTB04Analysis.cc:193
BeginOfRun.h
HcalTB04Analysis::eecals
double eecals
Definition: HcalTB04Analysis.cc:133
HcalTB04Analysis::init
void init()
Definition: HcalTB04Analysis.cc:216
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HcalTestNumbering::unpackHcalIndex
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
Definition: HcalTestNumbering.cc:18
HcalTB04Histo::fillPrimary
void fillPrimary(double energy, double eta, double phi)
Definition: HcalTB04Histo.cc:78
HcalTB04Analysis::ecalNoise
double ecalNoise
Definition: HcalTB04Analysis.cc:114
HcalTB04Analysis::enois
std::vector< double > enois
Definition: HcalTB04Analysis.cc:131
HcalTB04Analysis::secMomentum
std::vector< G4ThreeVector > secMomentum
Definition: HcalTB04Analysis.cc:139
PVValHelper::eta
Definition: PVValidationHelpers.h:69
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HcalTBNumberingScheme.h
ECalSD::setNumberingScheme
void setNumberingScheme(EcalNumberingScheme *)
Definition: ECalSD.cc:333
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:77
DDAxes::z
HcalTB04Analysis::pvPosition
G4ThreeVector pvPosition
Definition: HcalTB04Analysis.cc:137
HcalTB04Analysis::myQie
HcalQie * myQie
Definition: HcalTB04Analysis.cc:108
PHcalTB04Info::setLongProf
void setLongProf(const std::vector< double > &es, const std::vector< double > &eq)
Definition: PHcalTB04Info.cc:171
HcalTB04Analysis::phiInit
double phiInit
Definition: HcalTB04Analysis.cc:128
HcalTB04Analysis::icphi
int icphi
Definition: HcalTB04Analysis.cc:115
HcalTB04Analysis::esimh
std::vector< double > esimh
Definition: HcalTB04Analysis.cc:131
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
HcalTB04Analysis::etots
double etots
Definition: HcalTB04Analysis.cc:133
HcalTB04Analysis::type
int type
Definition: HcalTB04Analysis.cc:113
EndOfEvent
Definition: EndOfEvent.h:6
CaloHitIdMore
Definition: CaloHit.h:41
ECalSD
Definition: ECalSD.h:29
PHcalTB04Info.h
HcalTB04Analysis::qieAnalysis
void qieAnalysis(CLHEP::HepRandomEngine *)
Definition: HcalTB04Analysis.cc:719
CLHEP
Definition: CocoaGlobals.h:27
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HcalTB04Analysis::hcalHitLayer
std::vector< CaloHit > hcalHitLayer
Definition: HcalTB04Analysis.cc:130
HcalTB04Analysis::eslay
std::vector< double > eslay
Definition: HcalTB04Analysis.cc:132
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
edm::ParameterSet
Definition: ParameterSet.h:47
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
HcalTestNumbering::packHcalIndex
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
Definition: HcalTestNumbering.cc:7
HcalTB04Analysis::secEkin
std::vector< double > secEkin
Definition: HcalTB04Analysis.cc:140
HcalTB04Analysis::names
std::vector< std::string > names
Definition: HcalTB04Analysis.cc:117
HcalTB04Analysis::finalAnalysis
void finalAnalysis()
Definition: HcalTB04Analysis.cc:829
Event.h
HcalTestNumberingScheme
Definition: HcalTestNumberingScheme.h:11
HcalTB04Analysis::secTrackID
std::vector< int > secTrackID
Definition: HcalTB04Analysis.cc:138
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
SimProducer
Definition: SimProducer.h:64
PHcalTB04Info::setVtxSec
void setVtxSec(int id, int pdg, double px, double py, double pz, double ek)
Definition: PHcalTB04Info.cc:218
type
type
Definition: SiPixelVCal_PayloadInspector.cc:37
SimWatcherFactory.h
GeV
const double GeV
Definition: MathUtil.h:16
BeginOfEvent.h
HcalTB04Analysis::fillEvent
void fillEvent(PHcalTB04Info &)
Definition: HcalTB04Analysis.cc:940
PHcalTB04Info::setVtxPrim
void setVtxPrim(int evNum, int type, double x, double y, double z, double u, double v, double w, double px, double py, double pz)
Definition: PHcalTB04Info.cc:198
StorageManager_cfg.e1
e1
Definition: StorageManager_cfg.py:16
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
ModuleDef.h
createfilelist.int
int
Definition: createfilelist.py:10
HcalTB04Analysis::eqphi
std::vector< double > eqphi
Definition: HcalTB04Analysis.cc:132
HcalTB04XtalNumberingScheme.h
HcalTB04Analysis::fillBuffer
void fillBuffer(const EndOfEvent *evt)
Definition: HcalTB04Analysis.cc:471
BeginOfEvent
Definition: BeginOfEvent.h:6
BeginOfRun
Definition: BeginOfRun.h:6
HcalTB04Analysis::update
void update(const BeginOfRun *run) override
This routine will be called when the appropriate signal arrives.
Definition: HcalTB04Analysis.cc:277
HcalQie::getCode
std::vector< int > getCode(int, const std::vector< CaloHit > &, CLHEP::HepRandomEngine *)
Definition: HcalQie.cc:268
CaloHit
Definition: CaloHit.h:12
edm::EventSetup
Definition: EventSetup.h:57
HcalTB04Histo::fillEdep
void fillEdep(double etots, double eecals, double ehcals, double etotq, double eecalq, double ehcalq)
Definition: HcalTB04Histo.cc:85
HcalSubdetector.h
CaloG4Hit
Definition: CaloG4Hit.h:32
HcalTB04Analysis::nPrimary
int nPrimary
Definition: HcalTB04Analysis.cc:127
itr
std::vector< std::pair< float, float > >::iterator itr
Definition: HGCDigitizer.cc:29
HcalTB04Analysis::shortLivedSecondaries
std::vector< int > shortLivedSecondaries
Definition: HcalTB04Analysis.cc:141
ECalSD.h
HcalNumberingScheme
Definition: HcalNumberingScheme.h:13
HCalSD.h
HcalTBNumberingScheme::getUnitIDs
static std::vector< uint32_t > getUnitIDs(const int type, const int mode)
Definition: HcalTBNumberingScheme.cc:69
HcalTB04Analysis::etotq
double etotq
Definition: HcalTB04Analysis.cc:133
HcalTB04Analysis::count
int count
Definition: HcalTB04Analysis.cc:121
PHcalTB04Info::setEdepHcal
void setEdepHcal(const std::vector< double > &esim, const std::vector< double > &edig)
Definition: PHcalTB04Info.cc:133
HcalTB04Analysis::secPartID
std::vector< int > secPartID
Definition: HcalTB04Analysis.cc:138
DDAxes::phi
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:564
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
HcalTB04Analysis::eqie
std::vector< double > eqie
Definition: HcalTB04Analysis.cc:131
eostools.move
def move(src, dest)
Definition: eostools.py:511
writedatasetfile.run
run
Definition: writedatasetfile.py:27
HcalTB04Analysis::scaleHB0
double scaleHB0
Definition: HcalTB04Analysis.cc:116
HcalTB04Analysis::idHcal
std::vector< int > idHcal
Definition: HcalTB04Analysis.cc:123
HcalTB04Analysis::evNum
int evNum
Definition: HcalTB04Analysis.cc:136
HcalTB04Analysis::eqeta
std::vector< double > eqeta
Definition: HcalTB04Analysis.cc:132
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:31
Point3D.h
EventSetup.h
HcalTB04Analysis::hcalOnly
bool hcalOnly
Definition: HcalTB04Analysis.cc:112
HcalTB04Analysis::idEcal
std::vector< uint32_t > idEcal
Definition: HcalTB04Analysis.cc:124
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
CaloG4HitCollection
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
Definition: CaloG4HitCollection.h:11
HcalTB04Analysis::clear
void clear()
Definition: HcalTB04Analysis.cc:1013
HcalTB04Analysis::mode
int mode
Definition: HcalTB04Analysis.cc:113
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
HcalTB04Analysis::scaleHB16
double scaleHB16
Definition: HcalTB04Analysis.cc:116
HcalQie::getEnergy
double getEnergy(const std::vector< int > &)
Definition: HcalQie.cc:354
sd
double sd
Definition: CascadeWrapper.h:113
hippyaddtobaddatafiles.cd
def cd(newdir)
Definition: hippyaddtobaddatafiles.py:40
HcalTB04Analysis::nCrystal
int nCrystal
Definition: HcalTB04Analysis.cc:122
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
HcalTB04Analysis::scaleHE0
double scaleHE0
Definition: HcalTB04Analysis.cc:116
PHcalTB04Info
Definition: PHcalTB04Info.h:33
HcalTB04Analysis::pvFound
bool pvFound
Definition: HcalTB04Analysis.cc:135
HcalTBNumberingScheme::getUnitID
static uint32_t getUnitID(const uint32_t id, const int mode)
Definition: HcalTBNumberingScheme.cc:26
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
ntuplemaker.time
time
Definition: ntuplemaker.py:310
HcalQie
Definition: HcalQie.h:18
edm::Event
Definition: Event.h:73
HcalTestNumberingScheme.h
edm::Log
Definition: MessageLogger.h:70
CaloG4Hit.h
HcalTB04Analysis::eseta
std::vector< double > eseta
Definition: HcalTB04Analysis.cc:132
HcalTB04Analysis::pInit
double pInit
Definition: HcalTB04Analysis.cc:128
CaloHit.h
TauDecayModes.dec
dec
Definition: TauDecayModes.py:143
HcalTB04Histo::fillLongProf
void fillLongProf(const std::vector< double > &es, const std::vector< double > &eq)
Definition: HcalTB04Histo.cc:126
HcalTB04XtalNumberingScheme
Definition: HcalTB04XtalNumberingScheme.h:24
hit
Definition: SiStripHitEffFromCalibTree.cc:88
CaloG4Hit::getEntry
math::XYZPoint getEntry() const
Definition: CaloG4Hit.h:46
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
watchdog.group
group
Definition: watchdog.py:82
HcalTB04Analysis::scaleHO
double scaleHO
Definition: HcalTB04Analysis.cc:116
HcalTB04Analysis::etaInit
double etaInit
Definition: HcalTB04Analysis.cc:128
Generated for CMSSW Reference Manual by doxygen 1.8.16