CMS 3D CMS Logo
HcalTestAnalysis.cc
Go to the documentation of this file.
1 #include "SimG4Core/Notification/interface/BeginOfJob.h"
2 #include "SimG4Core/Notification/interface/BeginOfRun.h"
3 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
4 #include "SimG4Core/Notification/interface/EndOfEvent.h"
5 
6 // to retreive hits
7 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
8 #include "SimG4CMS/Calo/interface/HcalTestAnalysis.h"
9 #include "SimG4CMS/Calo/interface/HCalSD.h"
10 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
11 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
12 #include "DataFormats/Math/interface/Point3D.h"
13 
14 #include "FWCore/Framework/interface/EventSetup.h"
15 #include "FWCore/Framework/interface/ESHandle.h"
16 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
17 
18 #include "G4SDManager.hh"
19 #include "G4Step.hh"
20 #include "G4Track.hh"
21 #include "G4VProcess.hh"
22 #include "G4HCofThisEvent.hh"
23 #include "CLHEP/Units/GlobalSystemOfUnits.h"
24 #include "CLHEP/Units/GlobalPhysicalConstants.h"
25 #include "Randomize.hh"
26 
27 #include <cmath>
28 #include <iostream>
29 #include <iomanip>
30 
31 HcalTestAnalysis::HcalTestAnalysis(const edm::ParameterSet& p)
32  : addTower_(3), tuples_(nullptr), hcons_(nullptr), org_(nullptr) {
33  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTestAnalysis");
34  eta0_ = m_Anal.getParameter<double>("Eta0");
35  phi0_ = m_Anal.getParameter<double>("Phi0");
36  int laygroup = m_Anal.getParameter<int>("LayerGrouping");
37  centralTower_ = m_Anal.getParameter<int>("CentralTower");
38  names_ = m_Anal.getParameter<std::vector<std::string> >("Names");
39  fileName_ = m_Anal.getParameter<std::string>("FileName");
40 
41  tuplesManager_.reset(nullptr);
42  numberingFromDDD_.reset(nullptr);
43  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Initialised as observer of begin/end events"
44  << " and of G4step";
45 
46  count_ = 0;
47  group_ = layerGrouping(laygroup);
48  nGroup_ = 0;
49  for (unsigned int i = 0; i < group_.size(); i++)
50  if (group_[i] > nGroup_)
51  nGroup_ = group_[i];
52  tower_ = towersToAdd(centralTower_, addTower_);
53  nTower_ = tower_.size() / 2;
54 
55  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: initialised for " << nGroup_ << " Longitudinal groups and "
56  << nTower_ << " towers";
57 
58  // qie
59  myqie_.reset(new HcalQie(p));
60 }
61 
62 HcalTestAnalysis::~HcalTestAnalysis() {
63  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: --------> Total number of selected entries : " << count_;
64  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Pointers:: HistoClass " << tuples_ << ", Numbering Scheme " << org_;
65 }
66 
67 std::vector<int> HcalTestAnalysis::layerGrouping(int group) {
68  std::vector<int> temp(19);
69  if (group <= 1) {
70  int grp[19] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2};
71  for (int i = 0; i < 19; i++)
72  temp[i] = grp[i];
73  } else if (group == 2) {
74  int grp[19] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
75  for (int i = 0; i < 19; i++)
76  temp[i] = grp[i];
77  } else if (group == 3) {
78  int grp[19] = {1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7};
79  for (int i = 0; i < 19; i++)
80  temp[i] = grp[i];
81  } else if (group == 4) {
82  int grp[19] = {1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7};
83  for (int i = 0; i < 19; i++)
84  temp[i] = grp[i];
85  } else {
86  int grp[19] = {1, 1, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7};
87  for (int i = 0; i < 19; i++)
88  temp[i] = grp[i];
89  }
90 
91  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Layer Grouping ";
92  for (int i = 0; i < 19; i++)
93  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Group[" << i << "] = " << temp[i];
94  return temp;
95 }
96 
97 std::vector<int> HcalTestAnalysis::towersToAdd(int centre, int nadd) {
98  int etac = (centre / 100) % 100;
99  int phic = (centre % 100);
100  int etamin, etamax, phimin, phimax;
101  if (etac > 0) {
102  etamin = etac - nadd;
103  etamax = etac + nadd;
104  } else {
105  etamin = etac;
106  etamax = etac;
107  }
108  if (phic > 0) {
109  phimin = phic - nadd;
110  phimax = phic + nadd;
111  } else {
112  phimin = phic;
113  phimax = phic;
114  }
115 
116  int nbuf, kount = 0;
117  nbuf = (etamax - etamin + 1) * (phimax - phimin + 1);
118  std::vector<int> temp(2 * nbuf);
119  for (int eta = etamin; eta <= etamax; eta++) {
120  for (int phi = phimin; phi <= phimax; phi++) {
121  temp[kount] = (eta * 100 + phi);
122  temp[kount + nbuf] = std::max(abs(eta - etac), abs(phi - phic));
123  kount++;
124  }
125  }
126 
127  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Towers to be considered for Central " << centre << " and " << nadd
128  << " on either side";
129  for (int i = 0; i < nbuf; i++)
130  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Tower[" << std::setw(3) << i << "] " << temp[i] << " "
131  << temp[nbuf + i];
132  return temp;
133 }
134 
135 //==================================================================== per JOB
136 
137 void HcalTestAnalysis::update(const BeginOfJob* job) {
138  // Numbering From DDD
139  edm::ESHandle<HcalDDDSimConstants> hdc;
140  (*job)()->get<HcalSimNumberingRecord>().get(hdc);
141  hcons_ = hdc.product();
142  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Initialise HcalNumberingFromDDD for " << names_[0];
143  numberingFromDDD_.reset(new HcalNumberingFromDDD(hcons_));
144 
145  // Ntuples
146  tuplesManager_.reset(new HcalTestHistoManager(fileName_));
147 
148  // Numbering scheme
149  org_ = new HcalTestNumberingScheme(false);
150 }
151 
152 //==================================================================== per RUN
153 void HcalTestAnalysis::update(const BeginOfRun* run) {
154  int irun = (*run)()->GetRunID();
155  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Begin of Run = " << irun;
156 
157  bool loop = true, eta = true, phi = true;
158  int etac = (centralTower_ / 100) % 100;
159  if (etac == 0) {
160  etac = 1;
161  eta = false;
162  }
163  int phic = (centralTower_ % 100);
164  if (phic == 0) {
165  phic = 1;
166  phi = false;
167  }
168  int idet = static_cast<int>(HcalBarrel);
169  while (loop) {
170  HcalCellType::HcalCell tmp = hcons_->cell(idet, 1, 1, etac, phic);
171  if (tmp.ok) {
172  if (eta)
173  eta0_ = tmp.eta;
174  if (phi)
175  phi0_ = tmp.phi;
176  loop = false;
177  } else if (idet == static_cast<int>(HcalBarrel)) {
178  idet = static_cast<int>(HcalEndcap);
179  } else if (idet == static_cast<int>(HcalEndcap)) {
180  idet = static_cast<int>(HcalForward);
181  } else {
182  loop = false;
183  }
184  }
185 
186  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: Central Tower " << centralTower_
187  << " corresponds to eta0 = " << eta0_ << " phi0 = " << phi0_;
188 
189  std::string sdname = names_[0];
190  G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
191  if (sd != nullptr) {
192  G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
193  if (aSD == nullptr) {
194  edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: No SD with "
195  << "name " << sdname << " in this Setup";
196  } else {
197  HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
198  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::beginOfRun: Finds SD with name " << theCaloSD->GetName()
199  << " in this Setup";
200  if (org_) {
201  theCaloSD->setNumberingScheme(org_);
202  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::beginOfRun: set a new numbering scheme";
203  }
204  }
205  } else {
206  edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: Could not get"
207  << " SD Manager!";
208  }
209 }
210 
211 //=================================================================== per EVENT
212 void HcalTestAnalysis::update(const BeginOfEvent* evt) {
213  // create tuple object
214  tuples_ = new HcalTestHistoClass();
215  // Reset counters
216  tuples_->setCounters();
217 
218  int i = 0;
219  edepEB_ = edepEE_ = edepHB_ = edepHE_ = edepHO_ = 0.;
220  for (i = 0; i < 20; i++)
221  edepl_[i] = 0.;
222  for (i = 0; i < 20; i++)
223  mudist_[i] = -1.;
224 
225  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Begin of event = " << (*evt)()->GetEventID();
226 }
227 
228 //=================================================================== each STEP
229 void HcalTestAnalysis::update(const G4Step* aStep) {
230  if (aStep != nullptr) {
231  G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
232  G4String name = curPV->GetName();
233  name.assign(name, 0, 3);
234  double edeposit = aStep->GetTotalEnergyDeposit();
235  int layer = -1;
236  if (name == "EBR") {
237  edepEB_ += edeposit;
238  } else if (name == "EFR") {
239  edepEE_ += edeposit;
240  } else if (name == "HBS") {
241  layer = (curPV->GetCopyNo() / 10) % 100;
242  if (layer >= 0 && layer < 17) {
243  edepHB_ += edeposit;
244  } else {
245  edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HB " << curPV->GetName() << curPV->GetCopyNo();
246  layer = -1;
247  }
248  } else if (name == "HES") {
249  layer = (curPV->GetCopyNo() / 10) % 100;
250  if (layer >= 0 && layer < 19) {
251  edepHE_ += edeposit;
252  } else {
253  edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HE " << curPV->GetName() << curPV->GetCopyNo();
254  layer = -1;
255  }
256  } else if (name == "HTS") {
257  layer = (curPV->GetCopyNo() / 10) % 100;
258  if (layer >= 17 && layer < 20) {
259  edepHO_ += edeposit;
260  } else {
261  edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HO " << curPV->GetName() << curPV->GetCopyNo();
262  layer = -1;
263  }
264  }
265  if (layer >= 0 && layer < 20) {
266  edepl_[layer] += edeposit;
267 
268  // Calculate the distance if it is a muon
269  G4String part = aStep->GetTrack()->GetDefinition()->GetParticleName();
270  if ((part == "mu-" || part == "mu+") && mudist_[layer] < 0) {
271  math::XYZPoint pos(aStep->GetPreStepPoint()->GetPosition().x(),
272  aStep->GetPreStepPoint()->GetPosition().y(),
273  aStep->GetPreStepPoint()->GetPosition().z());
274  double theta = pos.theta();
275  double eta = -log(tan(theta * 0.5));
276  double phi = pos.phi();
277  double dist = sqrt((eta - eta0_) * (eta - eta0_) + (phi - phi0_) * (phi - phi0_));
278  mudist_[layer] = dist * std::sqrt(pos.perp2());
279  }
280  }
281 
282  if (layer >= 0 && layer < 20) {
283  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: G4Step: " << name << " Layer " << std::setw(3) << layer
284  << " Edep " << std::setw(6) << edeposit / MeV << " MeV";
285  }
286  } else {
287  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: G4Step: Null Step";
288  }
289 }
290 
291 //================================================================ End of EVENT
292 void HcalTestAnalysis::update(const EndOfEvent* evt) {
293  ++count_;
294  // Fill event input
295  fill(evt);
296  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after Fill";
297 
298  // Qie analysis
299  CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
300  qieAnalysis(engine);
301  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after QieAnalysis";
302 
303  // Layers tuples filling
304  layerAnalysis();
305  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after LayerAnalysis";
306 
307  // Writing the data to the Tree
308  tuplesManager_.get()->fillTree(tuples_); // (no need to delete it...)
309  tuples_ = nullptr; // but avoid to reuse it...
310  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis:: --- after fillTree";
311 }
312 
313 //---------------------------------------------------
314 void HcalTestAnalysis::fill(const EndOfEvent* evt) {
315  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis: Fill event " << (*evt)()->GetEventID();
316 
317  // access to the G4 hit collections
318  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
319 
320  int nhc = 0, neb = 0, nef = 0, j = 0;
321  caloHitCache_.erase(caloHitCache_.begin(), caloHitCache_.end());
322 
323  // Hcal
324  int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[0]);
325  CaloG4HitCollection* theHCHC = (CaloG4HitCollection*)allHC->GetHC(HCHCid);
326  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[0] << " of ID " << HCHCid
327  << " is obtained at " << theHCHC;
328  int hchc_entries = theHCHC->entries();
329  if (HCHCid >= 0 && theHCHC != nullptr) {
330  for (j = 0; j < hchc_entries; j++) {
331  CaloG4Hit* aHit = (*theHCHC)[j];
332 
333  double e = aHit->getEnergyDeposit() / GeV;
334  double time = aHit->getTimeSlice();
335 
336  math::XYZPoint pos = aHit->getPosition();
337  double theta = pos.theta();
338  double eta = -log(tan(theta * 0.5));
339  double phi = pos.phi();
340 
341  uint32_t unitID = aHit->getUnitID();
342  int subdet, zside, layer, etaIndex, phiIndex, lay;
343  org_->unpackHcalIndex(unitID, subdet, zside, layer, etaIndex, phiIndex, lay);
344  double jitter = time - timeOfFlight(subdet, lay, eta);
345  if (jitter < 0)
346  jitter = 0;
347  CaloHit hit(subdet, lay, e, eta, phi, jitter, unitID);
348  caloHitCache_.push_back(hit);
349  nhc++;
350 
351  std::string det = "HB";
352  if (subdet == static_cast<int>(HcalForward)) {
353  det = "HF";
354  } else if (subdet == static_cast<int>(HcalEndcap)) {
355  if (etaIndex <= 20) {
356  det = "HES";
357  } else {
358  det = "HED";
359  }
360  }
361 
362  edm::LogVerbatim("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2) << layer << " time "
363  << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
364  << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
365  << e;
366  }
367  }
368  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::HCAL hits : " << nhc;
369 
370  // EB
371  int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[1]);
372  CaloG4HitCollection* theEBHC = (CaloG4HitCollection*)allHC->GetHC(EBHCid);
373  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[1] << " of ID " << EBHCid
374  << " is obtained at " << theEBHC;
375  int ebhc_entries = theEBHC->entries();
376  if (EBHCid >= 0 && theEBHC != nullptr) {
377  for (j = 0; j < ebhc_entries; j++) {
378  CaloG4Hit* aHit = (*theEBHC)[j];
379 
380  double e = aHit->getEnergyDeposit() / GeV;
381  double time = aHit->getTimeSlice();
382  std::string det = "EB";
383 
384  math::XYZPoint pos = aHit->getPosition();
385  double theta = pos.theta();
386  double eta = -log(tan(theta / 2.));
387  double phi = pos.phi();
388 
389  HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta, phi, 1, 1);
390  uint32_t unitID = org_->getUnitID(id);
391  int subdet, zside, layer, ieta, iphi, lay;
392  org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
393  subdet = 10;
394  layer = 0;
395  unitID = org_->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
396  CaloHit hit(subdet, lay, e, eta, phi, time, unitID);
397  caloHitCache_.push_back(hit);
398  neb++;
399  edm::LogVerbatim("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2) << layer << " time "
400  << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
401  << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
402  << e;
403  }
404  }
405  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::EB hits : " << neb;
406 
407  // EE
408  int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names_[2]);
409  CaloG4HitCollection* theEEHC = (CaloG4HitCollection*)allHC->GetHC(EEHCid);
410  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names_[2] << " of ID " << EEHCid
411  << " is obtained at " << theEEHC;
412  int eehc_entries = theEEHC->entries();
413  if (EEHCid >= 0 && theEEHC != nullptr) {
414  for (j = 0; j < eehc_entries; j++) {
415  CaloG4Hit* aHit = (*theEEHC)[j];
416 
417  double e = aHit->getEnergyDeposit() / GeV;
418  double time = aHit->getTimeSlice();
419  std::string det = "EE";
420 
421  math::XYZPoint pos = aHit->getPosition();
422  double theta = pos.theta();
423  double eta = -log(tan(theta / 2.));
424  double phi = pos.phi();
425 
426  HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta, phi, 1, 1);
427  uint32_t unitID = org_->getUnitID(id);
428  int subdet, zside, layer, ieta, iphi, lay;
429  org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
430  subdet = 11;
431  layer = 0;
432  unitID = org_->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
433  CaloHit hit(subdet, lay, e, eta, phi, time, unitID);
434  caloHitCache_.push_back(hit);
435  nef++;
436  edm::LogVerbatim("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2) << layer << " time "
437  << std::setw(6) << time << " theta " << std::setw(8) << theta << " eta "
438  << std::setw(8) << eta << " phi " << std::setw(8) << phi << " e " << std::setw(8)
439  << e;
440  }
441  }
442  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::EE hits : " << nef;
443 }
444 
445 //-----------------------------------------------------------------------------
446 void HcalTestAnalysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
447  //Fill tuple with hit information
448  int hittot = caloHitCache_.size();
449  tuples_->fillHits(caloHitCache_);
450 
451  //Get the index of the central tower
452  HcalNumberingFromDDD::HcalID id = numberingFromDDD_->unitID(eta0_, phi0_, 1, 1);
453  uint32_t unitID = org_->getUnitID(id);
454  int subdet, zside, layer, ieta, iphi, lay;
455  org_->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
456  int laymax = 0;
457  std::string det = "Unknown";
458  if (subdet == static_cast<int>(HcalBarrel)) {
459  laymax = 4;
460  det = "HB";
461  } else if (subdet == static_cast<int>(HcalEndcap)) {
462  laymax = 2;
463  det = "HES";
464  }
465  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: " << det << " Eta " << ieta << " Phi " << iphi << " Laymax "
466  << laymax << " Hits " << hittot;
467 
468  if (laymax > 0 && hittot > 0) {
469  std::vector<CaloHit> hits(hittot);
470  std::vector<double> eqielay(80, 0.0), esimlay(80, 0.0), esimtot(4, 0.0);
471  std::vector<double> eqietow(200, 0.0), esimtow(200, 0.0), eqietot(4, 0.0);
472  int etac = (centralTower_ / 100) % 100;
473  int phic = (centralTower_ % 100);
474 
475  for (int layr = 0; layr < nGroup_; layr++) {
476  /*
477  int layx, layy=20;
478  for (int i=0; i<20; i++)
479  if (group_[i] == layr+1 && i < layy) layy = i+1;
480  if (subdet == static_cast<int>(HcalBarrel)) {
481  if (layy < 2) layx = 0;
482  else if (layy < 17) layx = 1;
483  else if (layy == 17) layx = 2;
484  else layx = 3;
485  } else {
486  if (layy < 2) layx = 0;
487  else layx = 1;
488  }
489  */
490  for (int it = 0; it < nTower_; it++) {
491  int nhit = 0;
492  double esim = 0;
493  for (int k1 = 0; k1 < hittot; k1++) {
494  CaloHit hit = caloHitCache_[k1];
495  int subdetc = hit.det();
496  int layer = hit.layer();
497  int group = 0;
498  if (layer > 0 && layer < 20)
499  group = group_[layer];
500  if (subdetc == subdet && group == layr + 1) {
501  int zsidec, ietac, iphic, idx;
502  unitID = hit.id();
503  org_->unpackHcalIndex(unitID, subdetc, zsidec, layer, ietac, iphic, lay);
504  if (etac > 0 && phic > 0) {
505  idx = ietac * 100 + iphic;
506  } else if (etac > 0) {
507  idx = ietac * 100;
508  } else if (phic > 0) {
509  idx = iphic;
510  } else {
511  idx = 0;
512  }
513  if (zsidec == zside && idx == tower_[it]) {
514  hits[nhit] = hit;
515  edm::LogVerbatim("HcalSim") << "HcalTest: Hit " << nhit << " " << hit;
516  nhit++;
517  esim += hit.e();
518  }
519  }
520  }
521 
522  std::vector<int> cd = myqie_->getCode(nhit, hits, engine);
523  double eqie = myqie_->getEnergy(cd);
524 
525  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: Energy in layer " << layr << " Sim " << esim
526  << " After QIE " << eqie;
527  for (int i = 0; i < 4; i++) {
528  if (tower_[nTower_ + it] <= i) {
529  esimtot[i] += esim;
530  eqietot[i] += eqie;
531  esimlay[20 * i + layr] += esim;
532  eqielay[20 * i + layr] += eqie;
533  esimtow[50 * i + it] += esim;
534  eqietow[50 * i + it] += eqie;
535  }
536  }
537  }
538  }
539  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::Qie: Total energy " << esimtot[3] << " (SimHit) " << eqietot[3]
540  << " (After QIE)";
541 
542  std::vector<double> latphi(10);
543  int nt = 2 * addTower_ + 1;
544  for (int it = 0; it < nt; it++)
545  latphi[it] = it - addTower_;
546  for (int i = 0; i < 4; i++) {
547  double scals = 1, scalq = 1;
548  std::vector<double> latfs(10, 0.), latfq(10, 0.), longs(20), longq(20);
549  if (esimtot[i] > 0)
550  scals = 1. / esimtot[i];
551  if (eqietot[i] > 0)
552  scalq = 1. / eqietot[i];
553  for (int it = 0; it < nTower_; it++) {
554  int phib = it % nt;
555  latfs[phib] += scals * esimtow[50 * i + it];
556  latfq[phib] += scalq * eqietow[50 * i + it];
557  }
558  for (int layr = 0; layr <= nGroup_; layr++) {
559  longs[layr] = scals * esimlay[20 * i + layr];
560  longq[layr] = scalq * eqielay[20 * i + layr];
561  }
562  tuples_->fillQie(i, esimtot[i], eqietot[i], nGroup_, longs, longq, nt, latphi, latfs, latfq);
563  }
564  }
565 }
566 
567 //---------------------------------------------------
568 void HcalTestAnalysis::layerAnalysis() {
569  int i = 0;
570  edm::LogVerbatim("HcalSim") << "\n ===>>> HcalTestAnalysis: Energy deposit "
571  << "\n at EB : " << std::setw(6) << edepEB_ / MeV << "\n at EE : " << std::setw(6)
572  << edepEE_ / MeV << "\n at HB : " << std::setw(6) << edepHB_ / MeV
573  << "\n at HE : " << std::setw(6) << edepHE_ / MeV << "\n at HO : " << std::setw(6)
574  << edepHO_ / MeV << "\n ---- HcalTestAnalysis: Energy deposit in Layers";
575  for (i = 0; i < 20; i++)
576  edm::LogVerbatim("HcalSim") << " Layer " << std::setw(2) << i << " E " << std::setw(8) << edepl_[i] / MeV << " MeV";
577 
578  tuples_->fillLayers(edepl_, edepHO_, edepHB_ + edepHE_, mudist_);
579 }
580 
581 //---------------------------------------------------
582 double HcalTestAnalysis::timeOfFlight(int det, int layer, double eta) {
583  double theta = 2.0 * atan(exp(-eta));
584  double dist = 0.;
585  if (det == static_cast<int>(HcalBarrel)) {
586  const double rLay[19] = {1836.0,
587  1902.0,
588  1962.0,
589  2022.0,
590  2082.0,
591  2142.0,
592  2202.0,
593  2262.0,
594  2322.0,
595  2382.0,
596  2448.0,
597  2514.0,
598  2580.0,
599  2646.0,
600  2712.0,
601  2776.0,
602  2862.5,
603  3847.0,
604  4052.0};
605  if (layer > 0 && layer < 20)
606  dist += rLay[layer - 1] * mm / sin(theta);
607  } else {
608  const double zLay[19] = {4034.0,
609  4032.0,
610  4123.0,
611  4210.0,
612  4297.0,
613  4384.0,
614  4471.0,
615  4558.0,
616  4645.0,
617  4732.0,
618  4819.0,
619  4906.0,
620  4993.0,
621  5080.0,
622  5167.0,
623  5254.0,
624  5341.0,
625  5428.0,
626  5515.0};
627  if (layer > 0 && layer < 20)
628  dist += zLay[layer - 1] * mm / cos(theta);
629  }
630  double tmp = dist / c_light / ns;
631  edm::LogVerbatim("HcalSim") << "HcalTestAnalysis::timeOfFlight " << tmp << " for det/lay " << det << " " << layer
632  << " eta/theta " << eta << " " << theta / deg << " dist " << dist;
633  return tmp;
634 }
HcalTestAnalysis::HcalTestAnalysis
HcalTestAnalysis(const edm::ParameterSet &p)
Definition: HcalTestAnalysis.cc:31
HcalTestAnalysis::edepl_
double edepl_[20]
Definition: HcalTestAnalysis.h:90
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
HcalTestAnalysis::timeOfFlight
double timeOfFlight(int det, int layer, double eta)
Definition: HcalTestAnalysis.cc:582
HcalTestAnalysis::count_
unsigned int count_
Definition: HcalTestAnalysis.h:88
HcalTestNumberingScheme::getUnitID
uint32_t getUnitID(const HcalNumberingFromDDD::HcalID &id) override
Definition: HcalTestNumberingScheme.cc:21
mps_fire.i
i
Definition: mps_fire.py:355
HcalTestAnalysis::myqie_
std::unique_ptr< HcalQie > myqie_
Definition: HcalTestAnalysis.h:65
HcalTestHistoClass::fillLayers
void fillLayers(double el[], double ho, double hbhe, double muxy[])
Definition: HcalTestHistoClass.cc:21
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:92
HcalNumberingFromDDD::HcalID
Definition: HcalNumberingFromDDD.h:21
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
ESHandle.h
HcalSimNumberingRecord.h
BeginOfJob.h
HcalTestAnalysis::tower_
std::vector< int > tower_
Definition: HcalTestAnalysis.h:79
nt
int nt
Definition: AMPTWrapper.h:42
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
CaloG4HitCollection.h
HCalSD
Definition: HCalSD.h:38
phimin
float phimin
Definition: ReggeGribovPartonMCHadronizer.h:107
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:65
HcalTestAnalysis::~HcalTestAnalysis
~HcalTestAnalysis() override
Definition: HcalTestAnalysis.cc:62
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HcalNumberingFromDDD
Definition: HcalNumberingFromDDD.h:16
pos
Definition: PixelAliasList.h:18
HcalTestHistoManager
Definition: HcalTestHistoManager.h:22
MeV
const double MeV
EndOfEvent.h
HcalBarrel
Definition: HcalAssistant.h:33
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
CaloG4Hit::getTimeSlice
double getTimeSlice() const
Definition: CaloG4Hit.h:66
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
training_settings.idx
idx
Definition: training_settings.py:16
HcalTestAnalysis::edepEE_
double edepEE_
Definition: HcalTestAnalysis.h:89
HcalTestAnalysis::layerAnalysis
void layerAnalysis()
Definition: HcalTestAnalysis.cc:568
HcalDDDSimConstants::cell
HcalCellType::HcalCell cell(const int &det, const int &zside, const int &depth, const int &etaR, const int &iphi) const
Definition: HcalDDDSimConstants.cc:28
muonTiming_cfi.etamin
etamin
Definition: muonTiming_cfi.py:30
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
LEDCalibrationChannels.iphi
iphi
Definition: LEDCalibrationChannels.py:64
BeginOfRun.h
part
part
Definition: HCALResponse.h:20
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HcalTestAnalysis::edepHO_
double edepHO_
Definition: HcalTestAnalysis.h:90
PVValHelper::eta
Definition: PVValidationHelpers.h:69
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:77
edm::ESHandle
Definition: DTSurvey.h:22
HcalTestAnalysis::caloHitCache_
std::vector< CaloHit > caloHitCache_
Definition: HcalTestAnalysis.h:78
HcalTestAnalysis::edepHB_
double edepHB_
Definition: HcalTestAnalysis.h:89
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
BeginOfJob
Definition: BeginOfJob.h:8
HcalTestAnalysis::fileName_
std::string fileName_
Definition: HcalTestAnalysis.h:62
HcalTestAnalysis::centralTower_
int centralTower_
Definition: HcalTestAnalysis.h:85
EndOfEvent
Definition: EndOfEvent.h:6
HcalTestAnalysis::group_
std::vector< int > group_
Definition: HcalTestAnalysis.h:79
HcalTestHistoClass
Definition: HcalTestHistoClass.h:15
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::LogWarning
Definition: MessageLogger.h:141
HcalTestAnalysis::tuplesManager_
std::unique_ptr< HcalTestHistoManager > tuplesManager_
Definition: HcalTestAnalysis.h:69
HcalTestAnalysis::qieAnalysis
void qieAnalysis(CLHEP::HepRandomEngine *)
Definition: HcalTestAnalysis.cc:446
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
edm::ParameterSet
Definition: ParameterSet.h:36
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
HcalTestNumberingScheme
Definition: HcalTestNumberingScheme.h:11
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
GeV
const double GeV
Definition: MathUtil.h:16
heppy_loop.loop
loop
Definition: heppy_loop.py:28
HcalTestAnalysis::org_
HcalTestNumberingScheme * org_
Definition: HcalTestAnalysis.h:75
BeginOfEvent.h
HcalTestNumberingScheme::packHcalIndex
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
Definition: HcalTestNumberingScheme.cc:47
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
HcalTestAnalysis::layerGrouping
std::vector< int > layerGrouping(int)
Definition: HcalTestAnalysis.cc:67
HcalTestAnalysis::eta0_
double eta0_
Definition: HcalTestAnalysis.h:84
HcalTestAnalysis::numberingFromDDD_
std::unique_ptr< HcalNumberingFromDDD > numberingFromDDD_
Definition: HcalTestAnalysis.h:73
edm::LogVerbatim
Definition: MessageLogger.h:297
BeginOfEvent
Definition: BeginOfEvent.h:6
BeginOfRun
Definition: BeginOfRun.h:6
CaloHit
Definition: CaloHit.h:12
HcalSubdetector.h
HcalCellType::HcalCell
Definition: HcalCellType.h:15
CaloG4Hit
Definition: CaloG4Hit.h:32
HcalTestAnalysis::addTower_
int addTower_
Definition: HcalTestAnalysis.h:66
CaloG4Hit::getPosition
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:52
HcalTestAnalysis::mudist_
double mudist_[20]
Definition: HcalTestAnalysis.h:91
HCalSD.h
HcalForward
Definition: HcalAssistant.h:36
phimax
float phimax
Definition: ReggeGribovPartonMCHadronizer.h:106
DDAxes::phi
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:550
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
writedatasetfile.run
run
Definition: writedatasetfile.py:27
HcalTestAnalysis::edepEB_
double edepEB_
Definition: HcalTestAnalysis.h:89
HcalTestAnalysis::names_
std::vector< std::string > names_
Definition: HcalTestAnalysis.h:83
HcalEndcap
Definition: HcalAssistant.h:34
HcalTestNumberingScheme::unpackHcalIndex
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
Definition: HcalTestNumberingScheme.cc:51
HcalTestAnalysis::update
void update(const BeginOfJob *run) override
This routine will be called when the appropriate signal arrives.
Definition: HcalTestAnalysis.cc:137
HcalTestHistoClass::fillHits
void fillHits(std::vector< CaloHit > &)
Definition: HcalTestHistoClass.cc:38
muonTiming_cfi.etamax
etamax
Definition: muonTiming_cfi.py:23
HcalTestAnalysis::phi0_
double phi0_
Definition: HcalTestAnalysis.h:84
HcalTestAnalysis.h
Point3D.h
HcalTestAnalysis::hcons_
const HcalDDDSimConstants * hcons_
Definition: HcalTestAnalysis.h:74
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
CaloG4HitCollection
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
Definition: CaloG4HitCollection.h:11
HcalTestAnalysis::nTower_
int nTower_
Definition: HcalTestAnalysis.h:80
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
sd
double sd
Definition: CascadeWrapper.h:113
hippyaddtobaddatafiles.cd
def cd(newdir)
Definition: hippyaddtobaddatafiles.py:40
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HcalTestHistoClass::setCounters
void setCounters()
Definition: HcalTestHistoClass.cc:13
HcalTestHistoClass::fillQie
void fillQie(int id, double esimtot, double eqietot, int nGroup, const std::vector< double > &longs, const std::vector< double > &longq, int nTower, const std::vector< double > &latphi, const std::vector< double > &latfs, const std::vector< double > &latfq)
Definition: HcalTestHistoClass.cc:113
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
CaloG4Hit.h
HcalTestAnalysis::towersToAdd
std::vector< int > towersToAdd(int centre, int nadd)
Definition: HcalTestAnalysis.cc:97
HcalTestAnalysis::tuples_
HcalTestHistoClass * tuples_
Definition: HcalTestAnalysis.h:70
HcalTestAnalysis::nGroup_
int nGroup_
Definition: HcalTestAnalysis.h:80
hit
Definition: SiStripHitEffFromCalibTree.cc:88
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HcalTestAnalysis::fill
void fill(const EndOfEvent *ev)
Definition: HcalTestAnalysis.cc:314
watchdog.group
group
Definition: watchdog.py:82
HcalTestAnalysis::edepHE_
double edepHE_
Definition: HcalTestAnalysis.h:89
Generated for CMSSW Reference Manual by doxygen 1.8.16