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