CMS 3D CMS Logo

SimG4HcalValidation.cc
Go to the documentation of this file.
1 // File: SimG4Validation.cc
3 // Description: Main analysis class for Hcal Validation of G4 Hits
6 
11 
12 // to retreive hits
19 
24 
27 
28 #include "G4HCofThisEvent.hh"
29 #include "G4SDManager.hh"
30 #include "G4Step.hh"
31 #include "G4Track.hh"
32 #include "G4VProcess.hh"
33 
34 using namespace geant_units::operators;
35 
37  : jetf(nullptr), numberingFromDDD(nullptr), org(nullptr) {
38  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("SimG4HcalValidation");
39  infolevel = m_Anal.getParameter<int>("InfoLevel");
40  hcalOnly = m_Anal.getParameter<bool>("HcalClusterOnly");
41  applySampling = m_Anal.getParameter<bool>("HcalSampling");
42  coneSize = m_Anal.getParameter<double>("ConeSize");
43  ehitThreshold = m_Anal.getParameter<double>("EcalHitThreshold");
44  hhitThreshold = m_Anal.getParameter<double>("HcalHitThreshold");
45  timeLowlim = m_Anal.getParameter<double>("TimeLowLimit");
46  timeUplim = m_Anal.getParameter<double>("TimeUpLimit");
47  jetThreshold = m_Anal.getParameter<double>("JetThreshold");
48  eta0 = m_Anal.getParameter<double>("Eta0");
49  phi0 = m_Anal.getParameter<double>("Phi0");
50  names = m_Anal.getParameter<std::vector<std::string>>("Names");
51  labelLayer = m_Anal.getUntrackedParameter<std::string>("LabelLayerInfo", "HcalInfoLayer");
52  labelNxN = m_Anal.getUntrackedParameter<std::string>("LabelNxNInfo", "HcalInfoNxN");
53  labelJets = m_Anal.getUntrackedParameter<std::string>("LabelJetsInfo", "HcalInfoJets");
54 
55  produces<PHcalValidInfoLayer>(labelLayer);
56  if (infolevel > 0)
57  produces<PHcalValidInfoNxN>(labelNxN);
58  if (infolevel > 1)
59  produces<PHcalValidInfoJets>(labelJets);
60 
61  edm::LogVerbatim("ValidHcal") << "HcalTestAnalysis:: Initialised as observer of begin/end events and "
62  << "of G4step with Parameter values: \n\tInfoLevel = " << infolevel
63  << "\n\thcalOnly = " << hcalOnly << "\n\tapplySampling = " << applySampling
64  << "\n\tconeSize = " << coneSize << "\n\tehitThreshold = " << ehitThreshold
65  << "\n\thhitThreshold = " << hhitThreshold << "\n\tttimeLowlim = " << timeLowlim
66  << "\n\tttimeUplim = " << timeUplim << "\n\teta0 = " << eta0
67  << "\n\tphi0 = " << phi0 << "\nLabels (Layer): " << labelLayer
68  << " (NxN): " << labelNxN << " (Jets): " << labelJets;
69 
70  init();
71 }
72 
74  edm::LogVerbatim("ValidHcal") << "\n --------> Total number of selected entries"
75  << " : " << count << "\nPointers:: JettFinder " << jetf << ", Numbering Scheme " << org
76  << " and FromDDD " << numberingFromDDD;
77  if (jetf) {
78  edm::LogVerbatim("ValidHcal") << "Delete Jetfinder";
79  delete jetf;
80  jetf = nullptr;
81  }
82  if (numberingFromDDD) {
83  edm::LogVerbatim("ValidHcal") << "Delete HcalNumberingFromDDD";
84  delete numberingFromDDD;
85  numberingFromDDD = nullptr;
86  }
87 }
88 
90  std::unique_ptr<PHcalValidInfoLayer> productLayer(new PHcalValidInfoLayer);
91  layerAnalysis(*productLayer);
92  e.put(std::move(productLayer), labelLayer);
93 
94  if (infolevel > 0) {
95  std::unique_ptr<PHcalValidInfoNxN> productNxN(new PHcalValidInfoNxN);
96  nxNAnalysis(*productNxN);
97  e.put(std::move(productNxN), labelNxN);
98  }
99 
100  if (infolevel > 1) {
101  std::unique_ptr<PHcalValidInfoJets> productJets(new PHcalValidInfoJets);
102  jetAnalysis(*productJets);
103  e.put(std::move(productJets), labelJets);
104  }
105 }
106 
107 //==================================================================== per RUN
109  float sHB[4] = {117., 117., 178., 217.};
110  float sHE[3] = {178., 178., 178.};
111  float sHF[3] = {2.84, 2.09, 0.}; //
112 
113  float deta[4] = {0.0435, 0.1305, 0.2175, 0.3045};
114  float dphi[4] = {0.0436, 0.1309, 0.2182, 0.3054};
115 
116  int i = 0;
117  for (i = 0; i < 4; i++) {
118  scaleHB.push_back(sHB[i]);
119  }
120  for (i = 0; i < 3; i++) {
121  scaleHE.push_back(sHE[i]);
122  }
123  for (i = 0; i < 3; i++) {
124  scaleHF.push_back(sHF[i]);
125  }
126 
127  // window steps;
128  for (i = 0; i < 4; i++) {
129  dEta.push_back(deta[i]);
130  dPhi.push_back(dphi[i]);
131  }
132 
133  // jetfinder conse size setting
135 
136  // counter
137  count = 0;
138 }
139 
141  // Numbering From DDD
143  (*job)()->get<HcalSimNumberingRecord>().get(hdc);
144  const HcalDDDSimConstants *hcons = hdc.product();
145  edm::LogVerbatim("ValidHcal") << "HcalTestAnalysis:: Initialise "
146  << "HcalNumberingFromDDD";
148 
149  // Numbering scheme
150  org = new HcalTestNumberingScheme(false);
151 }
152 
154  int irun = (*run)()->GetRunID();
155 
156  edm::LogVerbatim("ValidHcal") << " =====> Begin of Run = " << irun;
157 
158  std::string sdname = names[0];
159  G4SDManager *sd = G4SDManager::GetSDMpointerIfExist();
160  if (sd != nullptr) {
161  G4VSensitiveDetector *aSD = sd->FindSensitiveDetector(sdname);
162  if (aSD == nullptr) {
163  edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: No SD"
164  << " with name " << sdname << " in this "
165  << "Setup";
166  } else {
167  HCalSD *theCaloSD = dynamic_cast<HCalSD *>(aSD);
168  edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::beginOfRun: Finds SD with name " << theCaloSD->GetName()
169  << " in this Setup";
170  if (org) {
171  theCaloSD->setNumberingScheme(org);
172  edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::beginOfRun: set a new numbering scheme";
173  }
174  }
175  } else {
176  edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: Could "
177  << "not get SD Manager!";
178  }
179 }
180 
181 //=================================================================== per EVENT
183  int i = 0;
184  edepEB = edepEE = edepHB = edepHE = edepHO = 0.;
185  for (i = 0; i < 5; i++)
186  edepd[i] = 0.;
187  for (i = 0; i < 20; i++)
188  edepl[i] = 0.;
189  vhitec = vhithc = enEcal = enHcal = 0;
190  // Cache reset
191  clear();
192 
193  int iev = (*evt)()->GetEventID();
194  LogDebug("ValidHcal") << "SimG4HcalValidation: =====> Begin of event = " << iev;
195 }
196 
197 //=================================================================== each STEP
198 void SimG4HcalValidation::update(const G4Step *aStep) {
199  if (aStep != nullptr) {
200  G4VPhysicalVolume *curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
201  G4String name = curPV->GetName();
202  name.assign(name, 0, 3);
203  double edeposit = aStep->GetTotalEnergyDeposit();
204  int layer = -1, depth = -1;
205  if (name == "EBR") {
206  depth = 0;
207  edepEB += edeposit;
208  } else if (name == "EFR") {
209  depth = 0;
210  edepEE += edeposit;
211  } else if (name == "HBS") {
212  layer = (curPV->GetCopyNo() / 10) % 100;
213  depth = (curPV->GetCopyNo()) % 10 + 1;
214  if (depth > 0 && depth < 4 && layer >= 0 && layer < 17) {
215  edepHB += edeposit;
216  } else {
217  edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
218  depth = -1;
219  layer = -1;
220  }
221  } else if (name == "HES") {
222  layer = (curPV->GetCopyNo() / 10) % 100;
223  depth = (curPV->GetCopyNo()) % 10 + 1;
224  if (depth > 0 && depth < 3 && layer >= 0 && layer < 19) {
225  edepHE += edeposit;
226  } else {
227  edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
228  depth = -1;
229  layer = -1;
230  }
231  } else if (name == "HTS") {
232  layer = (curPV->GetCopyNo() / 10) % 100;
233  depth = (curPV->GetCopyNo()) % 10 + 1;
234  if (depth > 3 && depth < 5 && layer >= 17 && layer < 20) {
235  edepHO += edeposit;
236  } else {
237  edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error " << curPV->GetName() << curPV->GetCopyNo();
238  depth = -1;
239  layer = -1;
240  }
241  }
242  if (depth >= 0 && depth < 5)
243  edepd[depth] += edeposit;
244  if (layer >= 0 && layer < 20)
245  edepl[layer] += edeposit;
246 
247  if (layer >= 0 && layer < 20) {
248  LogDebug("ValidHcal") << "SimG4HcalValidation:: G4Step: " << name << " Layer " << std::setw(3) << layer
249  << " Depth " << std::setw(2) << depth << " Edep " << std::setw(6) << edeposit / MeV
250  << " MeV";
251  }
252  }
253 }
254 
255 //================================================================ End of EVENT
257  count++;
258 
259  // Fill hits cache for jetfinding etc.
260  fill(evt);
261  LogDebug("ValidHcal") << "SimG4HcalValidation:: --- after Fill ";
262 }
263 
264 //---------------------------------------------------
266  LogDebug("ValidHcal") << "SimG4HcalValidation:Fill event " << (*evt)()->GetEventID();
267 
268  // access to the G4 hit collections
269  G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
270 
271  int nhc = 0, j = 0;
272 
273  // Hcal
274  int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[0]);
275  CaloG4HitCollection *theHCHC = (CaloG4HitCollection *)allHC->GetHC(HCHCid);
276  LogDebug("ValidHcal") << "SimG4HcalValidation :: Hit Collection for " << names[0] << " of ID " << HCHCid
277  << " is obtained at " << theHCHC;
278  if (HCHCid >= 0 && theHCHC != nullptr) {
279  int nhits = theHCHC->entries();
280  for (j = 0; j < nhits; j++) {
281  CaloG4Hit *aHit = (*theHCHC)[j];
282 
283  double e = aHit->getEnergyDeposit() / GeV;
284  double time = aHit->getTimeSlice();
285 
286  math::XYZPoint pos = aHit->getPosition();
287  double theta = pos.theta();
288  double eta = -log(tan(theta * 0.5));
289  double phi = pos.phi();
290 
291  uint32_t unitID = aHit->getUnitID();
292  int subdet, zside, layer, etaIndex, phiIndex, lay;
293  org->unpackHcalIndex(unitID, subdet, zside, layer, etaIndex, phiIndex, lay);
294 
295  // some logic to separate HO ...
296  layer--;
297  std::string det = "HB";
298  if (subdet == static_cast<int>(HcalForward)) {
299  det = "HF";
300  uint16_t depth = aHit->getDepth();
301  if (depth != 0) {
302  layer += 2;
303  lay += 2;
304  }
305  if (layer == 1)
306  vhithc += e;
307  else if (layer == 0)
308  vhitec += e;
309  else
310  edm::LogVerbatim("ValidHcal") << "SimG4HcalValidation::HitPMT " << subdet << " " << (2 * zside - 1) * etaIndex
311  << " " << phiIndex << " " << layer + 1 << " R " << pos.rho() << " Phi "
312  << convertRadToDeg(phi) << " Edep " << e << " Time " << time;
313  } else if (subdet == static_cast<int>(HcalEndcap)) {
314  if (etaIndex <= 20) {
315  det = "HES";
316  } else {
317  det = "HED";
318  }
319  }
320  LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Hcal " << det << " layer " << std::setw(2) << layer
321  << " lay " << std::setw(2) << lay << " time " << std::setw(6) << time << " theta "
322  << std::setw(8) << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8) << phi
323  << " e " << std::setw(8) << e << " ID 0x" << std::hex << unitID << " ID dec " << std::dec
324  << (int)unitID;
325 
326  // if desired, apply sampling factors in HCAL !!!
327  if (applySampling)
328  e *= getHcalScale(det, layer);
329 
330  // filter on time & energy
331  if (time >= timeLowlim && time <= timeUplim && e > hhitThreshold) {
332  enHcal += e;
333  CaloHit ahit(subdet, lay, e, eta, phi, time, unitID);
334  hitcache.push_back(ahit); // fill cache
335  ++nhc;
336  }
337  }
338  }
339  LogDebug("ValidHcal") << "SimG4HcalValidation:: HCAL hits : " << nhc;
340 
341  if (!hcalOnly) { //-------------------------- ECAL hits --------------------
342  int ndets = names.size();
343  if (ndets > 3)
344  ndets = 3;
345  for (int idty = 1; idty < ndets; idty++) {
346  std::string det = "EB";
347  if (idty == 2)
348  det = "EE";
349  else if (idty > 2)
350  det = "ES";
351 
352  int nec = 0;
353  int ECHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[idty]);
354  CaloG4HitCollection *theECHC = (CaloG4HitCollection *)allHC->GetHC(ECHCid);
355  LogDebug("ValidHcal") << "SimG4HcalValidation:: Hit Collection for " << names[idty] << " of ID " << ECHCid
356  << " is obtained at " << theECHC;
357  int theechc_entries = theECHC->entries();
358  if (ECHCid >= 0 && theECHC != nullptr) {
359  for (j = 0; j < theechc_entries; j++) {
360  CaloG4Hit *aHit = (*theECHC)[j];
361 
362  double e = aHit->getEnergyDeposit() / GeV;
363  double time = aHit->getTimeSlice();
364 
365  math::XYZPoint pos = aHit->getPosition();
366  double theta = pos.theta();
367  double eta = -log(tan(theta / 2.));
368  double phi = pos.phi();
370  uint32_t unitID = org->getUnitID(id);
371  int subdet, zside, layer, ieta, iphi, lay;
372  org->unpackHcalIndex(unitID, subdet, zside, layer, ieta, iphi, lay);
373  subdet = idty + 9;
374  layer = 0;
375  unitID = org->packHcalIndex(subdet, zside, layer, ieta, iphi, lay);
376 
377  // filter on time & energy
378  if (time >= timeLowlim && time <= timeUplim && e > ehitThreshold) {
379  enEcal += e;
380  CaloHit ahit(subdet, lay, e, eta, phi, time, unitID);
381  hitcache.push_back(ahit); // fill cache
382  ++nec;
383  }
384 
385  LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Ecal " << det << " layer " << std::setw(2) << layer
386  << " lay " << std::setw(2) << lay << " time " << std::setw(6) << time << " theta "
387  << std::setw(8) << theta << " eta " << std::setw(8) << eta << " phi " << std::setw(8)
388  << phi << " e " << std::setw(8) << e << " ID 0x" << std::hex << unitID << " ID dec "
389  << std::dec << (int)unitID;
390  }
391  }
392 
393  LogDebug("ValidHcal") << "SimG4HcalValidation:: " << det << " hits : " << nec;
394  }
395  } // end of if(!hcalOnly)
396 }
397 
399  int i = 0;
400  LogDebug("ValidHcal") << "\n ===>>> SimG4HcalValidation: Energy deposit "
401  << "in MeV\n at EB : " << std::setw(6) << edepEB / MeV << "\n at EE : " << std::setw(6)
402  << edepEE / MeV << "\n at HB : " << std::setw(6) << edepHB / MeV
403  << "\n at HE : " << std::setw(6) << edepHE / MeV << "\n at HO : " << std::setw(6)
404  << edepHO / MeV << "\n ---- SimG4HcalValidation: Energy deposit in";
405  for (i = 0; i < 5; i++)
406  LogDebug("ValidHcal") << " Depth " << std::setw(2) << i << " E " << std::setw(8) << edepd[i] / MeV << " MeV";
407  LogDebug("ValidHcal") << " ---- SimG4HcalValidation: Energy deposit in"
408  << "layers";
409  for (i = 0; i < 20; i++)
410  LogDebug("ValidHcal") << " Layer " << std::setw(2) << i << " E " << std::setw(8) << edepl[i] / MeV << " MeV";
411 
413 
414  // Hits in HF
415  product.fillHF(vhitec, vhithc, enEcal, enHcal);
416  LogDebug("ValidHcal") << "SimG4HcalValidation::HF hits " << vhitec << " in EC and " << vhithc << " in HC\n"
417  << " HB/HE " << enEcal << " in EC and " << enHcal << " in HC";
418 
419  // Another HCAL hist to porcess and store separately (a bit more complicated)
420  fetchHits(product);
421 
422  LogDebug("ValidHcal") << " layerAnalysis ===> after fetchHits";
423 }
424 
425 //-----------------------------------------------------------------------------
427  std::vector<CaloHit> *hits = &hitcache;
428  std::vector<CaloHit>::iterator hit_itr;
429 
430  LogDebug("ValidHcal") << "SimG4HcalValidation::NxNAnalysis : entrance ";
431 
432  collectEnergyRdir(eta0, phi0); // HCAL and ECAL energy in SimHitCache
433  // around (eta0,phi0)
434 
435  LogDebug("ValidHcal") << " NxNAnalysis : coolectEnergyRdir - Ecal " << een << " Hcal " << hen;
436 
437  double etot = 0.; // total e deposited in "cluster"
438  double ee = 0.; // ECAL e deposited in "cluster"
439  double he = 0.; // HCAL e deposited in "cluster"
440  double hoe = 0.; // HO e deposited in "cluster"
441 
442  int max = dEta.size(); // 4
443 
444  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
445  double e = hit_itr->e();
446  double t = hit_itr->t();
447  double eta = hit_itr->eta();
448  double phi = hit_itr->phi();
449  int type = hit_itr->det();
450  int layer = hit_itr->layer();
451 
452  // NxN calulation
453 
454  if (fabs(eta0 - eta) <= dEta[max - 1] && fabs(phi0 - phi) <= dPhi[max - 1]) {
455  etot += e;
456  if (type == 10 || type == 11 || type == 12)
457  ee += e;
458  if (type == static_cast<int>(HcalBarrel) || type == static_cast<int>(HcalEndcap) ||
459  type == static_cast<int>(HcalForward)) {
460  he += e;
461  if (type == static_cast<int>(HcalBarrel) && layer > 17)
462  hoe += e;
463 
464  // which concrete i-th square ?
465  for (int i = 0; i < max; i++) {
466  if ((eta0 - eta) <= dEta[i] && (eta0 - eta) >= -dEta[i] && (phi0 - phi) <= dPhi[i] &&
467  (phi0 - phi) >= -dPhi[i]) {
468  LogDebug("ValidHcal") << "SimG4HcalValidation:: nxNAnalysis eta0,"
469  << " phi0 = " << eta0 << " " << phi0 << " type, layer = " << type << "," << layer
470  << " eta, phi = " << eta << " " << phi;
471 
472  product.fillTProfileNxN(e, i, t);
473  break;
474  }
475  }
476  // here comes break ...
477  }
478  }
479  }
480 
481  product.fillEcollectNxN(ee, he, hoe, etot);
482  product.fillHvsE(een, hen, hoen, een + hen);
483 
484  LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
485 }
486 
487 //-----------------------------------------------------------------------------
489  std::vector<CaloHit> *hhit = &hitcache;
490 
491  jetf->setInput(hhit);
492 
493  // zeroing cluster list, perfom clustering, fill cluster list & return pntr
494  std::vector<SimG4HcalHitCluster> *result = jetf->getClusters(hcalOnly);
495 
496  std::vector<SimG4HcalHitCluster>::iterator clus_itr;
497 
498  LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size() << " clusters ---------------";
499  for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
500  LogDebug("ValidHcal") << (*clus_itr);
501 
502  std::vector<double> enevec, etavec, phivec;
503 
504  if (!(*result).empty()) {
505  sort((*result).begin(), (*result).end());
506 
507  clus_itr = result->begin(); // first cluster only
508  double etac = clus_itr->eta();
509  double phic = clus_itr->phi();
510 
511  double ecal_collect = 0.; // collect Ecal energy in the cone
512  if (!hcalOnly) {
513  ecal_collect = clus_itr->collectEcalEnergyR();
514  } else {
515  collectEnergyRdir(etac, phic);
516  ecal_collect = een;
517  }
518  LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect = " << ecal_collect;
519 
520  // eta-phi deviation of the cluster from nominal (eta0, phi0) values
521  double dist = jetf->rDist(eta0, phi0, etac, phic);
522  LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation = " << dist;
523  product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
524 
525  LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
526 
527  std::vector<CaloHit> *hits = clus_itr->getHits();
528  std::vector<CaloHit>::iterator hit_itr;
529 
530  double ee = 0., he = 0., hoe = 0., etot = 0.;
531 
532  // cycle over all hits in the FIRST cluster
533  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
534  double e = hit_itr->e();
535  double t = hit_itr->t();
536  double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
537 
538  // energy collection
539  etot += e;
540  if (hit_itr->det() == 10 || hit_itr->det() == 11 || hit_itr->det() == 12)
541  ee += e;
542  if (hit_itr->det() == static_cast<int>(HcalBarrel) || hit_itr->det() == static_cast<int>(HcalEndcap) ||
543  hit_itr->det() == static_cast<int>(HcalForward)) {
544  he += e;
545  if (hit_itr->det() == static_cast<int>(HcalBarrel) && hit_itr->layer() > 17)
546  hoe += e;
547  }
548 
549  if (hit_itr->det() == static_cast<int>(HcalBarrel) || hit_itr->det() == static_cast<int>(HcalEndcap) ||
550  hit_itr->det() == static_cast<int>(HcalForward)) {
551  product.fillTProfileJet(he, r, t);
552  }
553  }
554 
555  product.fillEcollectJet(ee, he, hoe, etot);
556 
557  LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
558  << "ee/he/hoe/etot " << ee << "/" << he << "/" << hoe << "/" << etot;
559 
560  // Loop over clusters
561  for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
562  if (clus_itr->e() > jetThreshold) {
563  enevec.push_back(clus_itr->e());
564  etavec.push_back(clus_itr->eta());
565  phivec.push_back(clus_itr->phi());
566  }
567  }
568  product.fillJets(enevec, etavec, phivec);
569 
570  LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
571  << " JetAnalysis ===> (*result).size() " << (*result).size();
572 
573  // Di-jet mass
574  if (etavec.size() > 1) {
575  if (etavec[0] > -2.5 && etavec[0] < 2.5 && etavec[1] > -2.5 && etavec[1] < 2.5) {
576  LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
577  << " JetAnalysis ===> Di-jet vectors ";
578  for (unsigned int i = 0; i < enevec.size(); i++)
579  LogDebug("ValidHcal") << " e, eta, phi = " << enevec[i] << " " << etavec[i] << " " << phivec[i];
580 
581  double et0 = enevec[0] / cosh(etavec[0]);
582  double px0 = et0 * cos(phivec[0]);
583  double py0 = et0 * sin(phivec[0]);
584  double pz0 = et0 * sinh(etavec[0]);
585  double et1 = enevec[1] / cosh(etavec[1]);
586  double px1 = et1 * cos(phivec[1]);
587  double py1 = et1 * sin(phivec[1]);
588  double pz1 = et1 * sinh(etavec[1]);
589 
590  double dijetmass2 = (enevec[0] + enevec[1]) * (enevec[0] + enevec[1]) - (px1 + px0) * (px1 + px0) -
591  (py1 + py0) * (py1 + py0) - (pz1 + pz0) * (pz1 + pz0);
592 
593  LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ " << dijetmass2;
594 
595  double dijetmass;
596  if (dijetmass2 >= 0.)
597  dijetmass = sqrt(dijetmass2);
598  else
599  dijetmass = -sqrt(-1. * dijetmass2);
600 
601  product.fillDiJets(dijetmass);
602 
603  LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
604  }
605  }
606  }
607 }
608 
609 //---------------------------------------------------
611  LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with " << hitcache.size() << " hits";
612  int nHit = hitcache.size();
613  int hit = 0;
614  int i;
615  std::vector<CaloHit>::iterator itr;
616  std::vector<CaloHit *> lhits(nHit);
617  for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
618  uint32_t unitID = itr->id();
619  int subdet, zside, group, ieta, iphi, lay;
620  HcalTestNumbering::unpackHcalIndex(unitID, subdet, zside, group, ieta, iphi, lay);
621  subdet = itr->det();
622  lay = itr->layer();
623  group = (subdet & 15) << 20;
624  group += ((lay - 1) & 31) << 15;
625  group += (zside & 1) << 14;
626  group += (ieta & 127) << 7;
627  group += (iphi & 127);
628  itr->setId(group);
629  lhits[i] = &hitcache[i];
630  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i << " " << hitcache[i] << "\n"
631  << "SimG4HcalValidation::fetchHits:Copied " << i << " " << *lhits[i];
632  }
633  sort(lhits.begin(), lhits.end(), CaloHitIdMore());
634  std::vector<CaloHit *>::iterator k1, k2;
635  for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
636  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i << " " << **k1;
637  int nHits = 0;
638  for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
639  double ehit = (**k1).e();
640  double t = (**k1).t();
641  uint32_t unitID = (**k1).id();
642  int jump = 0;
643  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i << " U/T/E"
644  << " 0x" << std::hex << unitID << std::dec << " " << t << " " << ehit;
645  for (k2 = k1 + 1; k2 != lhits.end() && (t - (**k2).t()) < 1 && (t - (**k2).t()) > -1 && unitID == (**k2).id();
646  k2++) {
647  ehit += (**k2).e();
648  LogDebug("ValidHcal") << "\t + " << (**k2).e();
649  jump++;
650  }
651  LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
652 
653  double eta = (*k1)->eta();
654  double phi = (*k1)->phi();
655  int lay = ((unitID >> 15) & 31) + 1;
656  int subdet = (unitID >> 20) & 15;
657  int zside = (unitID >> 14) & 1;
658  int ieta = (unitID >> 7) & 127;
659  int iphi = (unitID)&127;
660 
661  // All hits in cache
662  product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
663  nHits++;
664 
665  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits << " " << i << " ID 0x" << std::hex
666  << unitID << " det " << std::dec << subdet << " " << lay << " " << zside << " " << ieta
667  << " " << iphi << " Time " << t << " E " << ehit;
668 
669  i += jump;
670  k1 += jump;
671  }
672 
673  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with " << nHit << " hits"
674  << " and writes out " << nHits << '(' << hit << ") hits";
675 }
676 //---------------------------------------------------
677 void SimG4HcalValidation::clear() { hitcache.erase(hitcache.begin(), hitcache.end()); }
678 
679 //---------------------------------------------------
680 void SimG4HcalValidation::collectEnergyRdir(const double eta0, const double phi0) {
681  std::vector<CaloHit> *hits = &hitcache;
682  std::vector<CaloHit>::iterator hit_itr;
683 
684  double sume = 0., sumh = 0., sumho = 0.;
685 
686  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
687  double e = hit_itr->e();
688  double eta = hit_itr->eta();
689  double phi = hit_itr->phi();
690 
691  int type = hit_itr->det();
692 
693  double r = jetf->rDist(eta0, phi0, eta, phi);
694  if (r < coneSize) {
695  if (type == 10 || type == 11 || type == 12) {
696  sume += e;
697  } else {
698  sumh += e;
699  if (type == static_cast<int>(HcalBarrel) && hit_itr->layer() > 17)
700  sumho += e;
701  }
702  }
703  }
704 
705  een = sume;
706  hen = sumh;
707  hoen = sumho;
708 }
709 
710 //---------------------------------------------------
711 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
712  double tmp = 0.;
713 
714  if (det == "HB") {
715  tmp = scaleHB[layer];
716  } else if (det == "HES" || det == "HED") {
717  tmp = scaleHE[layer];
718  } else if (det == "HF") {
719  tmp = scaleHF[layer];
720  }
721 
722  return tmp;
723 }
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
SimG4HcalValidation::scaleHF
std::vector< float > scaleHF
Definition: SimG4HcalValidation.h:85
SimG4HcalValidation::ehitThreshold
double ehitThreshold
Definition: SimG4HcalValidation.h:89
PHcalValidInfoNxN::fillHvsE
void fillHvsE(double ee, double he, double hoe, double etot)
Definition: PValidationFormats.cc:1323
SimG4HcalValidation::edepl
double edepl[20]
Definition: SimG4HcalValidation.h:102
HcalTestNumberingScheme::getUnitID
uint32_t getUnitID(const HcalNumberingFromDDD::HcalID &id) override
Definition: HcalTestNumberingScheme.cc:21
mps_fire.i
i
Definition: mps_fire.py:355
MessageLogger.h
HcalNumberingFromDDD::HcalID
Definition: HcalNumberingFromDDD.h:21
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
ESHandle.h
HcalSimNumberingRecord.h
BeginOfJob.h
PHcalValidInfoLayer::fillHF
void fillHF(double fibl, double fibs, double enec, double enhc)
Definition: PValidationFormats.cc:1300
SimG4HcalValidation::timeLowlim
float timeLowlim
Definition: SimG4HcalValidation.h:90
SimG4HcalValidation::edepHB
double edepHB
Definition: SimG4HcalValidation.h:101
SimG4HcalHitJetFinder
Definition: SimG4HcalHitJetFinder.h:13
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
CaloG4HitCollection.h
HCalSD
Definition: HCalSD.h:38
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:65
PHcalValidInfoLayer
Definition: PValidationFormats.h:1328
SimG4HcalValidation::scaleHB
std::vector< float > scaleHB
Definition: SimG4HcalValidation.h:83
angle_units::operators::convertRadToDeg
constexpr NumType convertRadToDeg(NumType radians)
Definition: angle_units.h:21
SimG4HcalValidation::dEta
std::vector< double > dEta
Definition: SimG4HcalValidation.h:96
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HcalNumberingFromDDD
Definition: HcalNumberingFromDDD.h:16
SimG4HcalHitJetFinder::getClusters
std::vector< SimG4HcalHitCluster > * getClusters(bool)
Definition: SimG4HcalHitJetFinder.cc:21
pos
Definition: PixelAliasList.h:18
SimG4HcalValidation::een
double een
Definition: SimG4HcalValidation.h:103
SimG4HcalValidation::edepEB
double edepEB
Definition: SimG4HcalValidation.h:101
PHcalValidInfoLayer::fillHits
void fillHits(int Nhits, int lay, int unitID, double eta, double phi, double ehit, double t)
Definition: PValidationFormats.cc:1307
MeV
const double MeV
SimG4HcalValidation::labelJets
std::string labelJets
Definition: SimG4HcalValidation.h:93
EndOfEvent.h
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
geant_units::operators
Definition: GeantUnits.h:18
SimG4HcalValidation::labelNxN
std::string labelNxN
Definition: SimG4HcalValidation.h:93
HcalBarrel
Definition: HcalAssistant.h:33
SimG4HcalValidation::hitcache
std::vector< CaloHit > hitcache
Definition: SimG4HcalValidation.h:80
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
PHcalValidInfoNxN
Definition: PValidationFormats.h:1393
CaloG4Hit::getTimeSlice
double getTimeSlice() const
Definition: CaloG4Hit.h:66
SimG4HcalValidation::hcalOnly
bool hcalOnly
Definition: SimG4HcalValidation.h:91
PHcalValidInfoJets::fillTProfileJet
void fillTProfileJet(double e, double r, double t)
Definition: PValidationFormats.cc:1243
SimG4HcalValidation::SimG4HcalValidation
SimG4HcalValidation(const edm::ParameterSet &p)
Definition: SimG4HcalValidation.cc:36
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
LEDCalibrationChannels.iphi
iphi
Definition: LEDCalibrationChannels.py:64
SimG4HcalValidation::scaleHE
std::vector< float > scaleHE
Definition: SimG4HcalValidation.h:84
HcalDDDSimConstants
Definition: HcalDDDSimConstants.h:24
SimG4HcalValidation::enHcal
double enHcal
Definition: SimG4HcalValidation.h:104
BeginOfRun.h
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
PHcalValidInfoJets::fillEtaPhiProfileJet
void fillEtaPhiProfileJet(double eta0, double phi0, double eta, double phi, double dist)
Definition: PValidationFormats.cc:1261
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
SimG4HcalValidation::org
HcalTestNumberingScheme * org
Definition: SimG4HcalValidation.h:77
SimG4HcalValidation::jetThreshold
float jetThreshold
Definition: SimG4HcalValidation.h:90
PVValHelper::eta
Definition: PVValidationHelpers.h:69
SimG4HcalValidation::eta0
float eta0
Definition: SimG4HcalValidation.h:90
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
SimG4HcalValidation::collectEnergyRdir
void collectEnergyRdir(const double, const double)
Definition: SimG4HcalValidation.cc:680
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
BeginOfJob
Definition: BeginOfJob.h:8
SimG4HcalValidation::jetAnalysis
void jetAnalysis(PHcalValidInfoJets &)
Definition: SimG4HcalValidation.cc:488
PHcalValidInfoNxN::fillEcollectNxN
void fillEcollectNxN(double een, double hen, double hoen, double etotn)
Definition: PValidationFormats.cc:1330
OrderedSet.t
t
Definition: OrderedSet.py:90
nhits
Definition: HIMultiTrackSelector.h:42
EndOfEvent
Definition: EndOfEvent.h:6
CaloHitIdMore
Definition: CaloHit.h:41
SimG4HcalValidation::layerAnalysis
void layerAnalysis(PHcalValidInfoLayer &)
Definition: SimG4HcalValidation.cc:398
LEDCalibrationChannels.depth
depth
Definition: LEDCalibrationChannels.py:65
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::LogWarning
Definition: MessageLogger.h:141
SimG4HcalValidation::clear
void clear()
Definition: SimG4HcalValidation.cc:677
SimG4HcalValidation::edepHO
double edepHO
Definition: SimG4HcalValidation.h:101
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
SimG4HcalValidation::hen
double hen
Definition: SimG4HcalValidation.h:103
SimG4HcalValidation::dPhi
std::vector< double > dPhi
Definition: SimG4HcalValidation.h:97
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:670
SimG4HcalHitJetFinder::setInput
void setInput(std::vector< CaloHit > *)
Definition: SimG4HcalHitJetFinder.cc:19
edm::ParameterSet
Definition: ParameterSet.h:36
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
SimG4HcalValidation::edepEE
double edepEE
Definition: SimG4HcalValidation.h:101
SimG4HcalValidation::init
void init()
Definition: SimG4HcalValidation.cc:108
SimG4HcalValidation::nxNAnalysis
void nxNAnalysis(PHcalValidInfoNxN &)
Definition: SimG4HcalValidation.cc:426
Event.h
HcalTestNumberingScheme
Definition: HcalTestNumberingScheme.h:11
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
GeantUnits.h
SimG4HcalValidation::hoen
double hoen
Definition: SimG4HcalValidation.h:103
GeV
const double GeV
Definition: MathUtil.h:16
PHcalValidInfoJets::fillJets
void fillJets(const std::vector< double > &enj, const std::vector< double > &etaj, const std::vector< double > &phij)
Definition: PValidationFormats.cc:1267
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
createfilelist.int
int
Definition: createfilelist.py:10
SimG4HcalValidation::fill
void fill(const EndOfEvent *ev)
Definition: SimG4HcalValidation.cc:265
PHcalValidInfoJets
Definition: PValidationFormats.h:1255
PHcalValidInfoLayer::fillLayers
void fillLayers(double el[], double ed[], double ho, double hbhe, double ebee)
Definition: PValidationFormats.cc:1286
edm::LogVerbatim
Definition: MessageLogger.h:297
BeginOfEvent
Definition: BeginOfEvent.h:6
BeginOfRun
Definition: BeginOfRun.h:6
CaloHit
Definition: CaloHit.h:12
edm::EventSetup
Definition: EventSetup.h:57
SimG4HcalValidation::update
void update(const BeginOfJob *job) override
This routine will be called when the appropriate signal arrives.
Definition: SimG4HcalValidation.cc:140
HcalSubdetector.h
CaloG4Hit
Definition: CaloG4Hit.h:32
CaloG4Hit::getPosition
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:52
SimG4HcalValidation::count
unsigned int count
Definition: SimG4HcalValidation.h:100
itr
std::vector< std::pair< float, float > >::iterator itr
Definition: HGCDigitizer.cc:28
PHcalValidInfoJets::fillDiJets
void fillDiJets(double mass)
Definition: PValidationFormats.cc:1280
SimG4HcalValidation::coneSize
double coneSize
Definition: SimG4HcalValidation.h:89
SimG4HcalValidation::names
std::vector< std::string > names
Definition: SimG4HcalValidation.h:88
HCalSD.h
PHcalValidInfoJets::fillEcollectJet
void fillEcollectJet(double ee, double he, double hoe, double etot)
Definition: PValidationFormats.cc:1252
SimG4HcalValidation::labelLayer
std::string labelLayer
Definition: SimG4HcalValidation.h:93
electrons_cff.hoe
hoe
Definition: electrons_cff.py:404
alignCSCRings.r
r
Definition: alignCSCRings.py:93
HcalForward
Definition: HcalAssistant.h:36
DDAxes::phi
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:550
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
SimG4HcalValidation::vhithc
double vhithc
Definition: SimG4HcalValidation.h:104
SimG4HcalValidation::applySampling
bool applySampling
Definition: SimG4HcalValidation.h:91
SimG4HcalHitJetFinder::rDist
double rDist(const SimG4HcalHitCluster *, const CaloHit *) const
Definition: SimG4HcalHitJetFinder.cc:111
type
type
Definition: HCALResponse.h:21
eostools.move
def move(src, dest)
Definition: eostools.py:511
writedatasetfile.run
run
Definition: writedatasetfile.py:27
PHcalValidInfoNxN::fillTProfileNxN
void fillTProfileNxN(double e, int i, double t)
Definition: PValidationFormats.cc:1337
SimG4HcalValidation::edepHE
double edepHE
Definition: SimG4HcalValidation.h:101
hcalSimParameters_cfi.he
he
Definition: hcalSimParameters_cfi.py:75
SimG4HcalValidation::timeUplim
float timeUplim
Definition: SimG4HcalValidation.h:90
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
HcalNumberingFromDDD::unitID
HcalID unitID(int det, const math::XYZVectorD &pos, int depth, int lay=-1) const
Definition: HcalNumberingFromDDD.cc:30
Point3D.h
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
SimG4HcalValidation::~SimG4HcalValidation
~SimG4HcalValidation() override
Definition: SimG4HcalValidation.cc:73
SimG4HcalValidation::phi0
float phi0
Definition: SimG4HcalValidation.h:90
CaloG4HitCollection
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
Definition: CaloG4HitCollection.h:11
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
sd
double sd
Definition: CascadeWrapper.h:113
SimG4HcalValidation::getHcalScale
double getHcalScale(std::string, int) const
Definition: SimG4HcalValidation.cc:711
SimG4HcalValidation::vhitec
double vhitec
Definition: SimG4HcalValidation.h:104
SimG4HcalValidation::hhitThreshold
double hhitThreshold
Definition: SimG4HcalValidation.h:89
mps_fire.result
result
Definition: mps_fire.py:303
SimG4HcalValidation::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: SimG4HcalValidation.cc:89
SimG4HcalValidation::enEcal
double enEcal
Definition: SimG4HcalValidation.h:104
SimG4HcalValidation::edepd
double edepd[5]
Definition: SimG4HcalValidation.h:102
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
ntuplemaker.time
time
Definition: ntuplemaker.py:310
SimG4HcalValidation::numberingFromDDD
HcalNumberingFromDDD * numberingFromDDD
Definition: SimG4HcalValidation.h:74
edm::Event
Definition: Event.h:73
CaloG4Hit.h
TauDecayModes.dec
dec
Definition: TauDecayModes.py:143
SimG4HcalValidation::jetf
SimG4HcalHitJetFinder * jetf
Definition: SimG4HcalValidation.h:71
hit
Definition: SiStripHitEffFromCalibTree.cc:88
CaloG4Hit::getDepth
uint16_t getDepth() const
Definition: CaloG4Hit.h:68
SimG4HcalValidation.h
SimG4HcalValidation::infolevel
int infolevel
Definition: SimG4HcalValidation.h:92
HcalDDDSimConstants.h
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
SimG4HcalValidation::fetchHits
void fetchHits(PHcalValidInfoLayer &)
Definition: SimG4HcalValidation.cc:610
watchdog.group
group
Definition: watchdog.py:82