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