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/ESTransientHandle.h"
22 #include "FWCore/MessageLogger/interface/MessageLogger.h"
23 
24 #include "Geometry/Records/interface/IdealGeometryRecord.h"
25 #include "DetectorDescription/Core/interface/DDCompactView.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::ESTransientHandle<DDCompactView> pDD;
157  (*job)()->get<IdealGeometryRecord>().get(pDD);
158  edm::LogInfo("ValidHcal") << "HcalTestAnalysis:: Initialise "
159  << "HcalNumberingFromDDD for " << names[0];
160  numberingFromDDD = new HcalNumberingFromDDD(names[0], (*pDD));
161 
162  // Numbering scheme
163  org = new HcalTestNumberingScheme(false);
164 
165 }
166 
167 void SimG4HcalValidation::update(const BeginOfRun * run) {
168 
169  int irun = (*run)()->GetRunID();
170 
171  edm::LogInfo("ValidHcal") <<" =====> Begin of Run = " << irun;
172 
173  std::string sdname = names[0];
174  G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
175  if (sd != 0) {
176  G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
177  if (aSD==0) {
178  edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: No SD"
179  << " with name " << sdname << " in this "
180  << "Setup";
181  } else {
182  HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
183  edm::LogInfo("ValidHcal") << "SimG4HcalValidation::beginOfRun: Finds SD"
184  << "with name " << theCaloSD->GetName()
185  << " in this Setup";
186  if (org) {
187  theCaloSD->setNumberingScheme(org);
188  edm::LogInfo("ValidHcal") << "SimG4HcalValidation::beginOfRun: set a "
189  << "new numbering scheme";
190  }
191  }
192  } else {
193  edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: Could "
194  << "not get SD Manager!";
195  }
196 
197 }
198 
199 //=================================================================== per EVENT
200 void SimG4HcalValidation::update(const BeginOfEvent * evt) {
201 
202  int i = 0;
203  edepEB = edepEE = edepHB = edepHE = edepHO = 0.;
204  for (i = 0; i < 5; i++) edepd[i] = 0.;
205  for (i = 0; i < 20; i++) edepl[i] = 0.;
206  vhitec = vhithc = enEcal = enHcal = 0;
207  // Cache reset
208  clear();
209 
210  int iev = (*evt)()->GetEventID();
211  LogDebug("ValidHcal") << "SimG4HcalValidation: =====> Begin of event = "
212  << iev;
213 }
214 
215 //=================================================================== each STEP
216 void SimG4HcalValidation::update(const G4Step * aStep) {
217 
218  if (aStep != NULL) {
219  G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
220  G4String name = curPV->GetName();
221  name.assign(name,0,3);
222  double edeposit = aStep->GetTotalEnergyDeposit();
223  int layer=-1, depth=-1;
224  if (name == "EBR") {
225  depth = 0;
226  edepEB += edeposit;
227  } else if (name == "EFR") {
228  depth = 0;
229  edepEE += edeposit;
230  } else if (name == "HBS") {
231  layer = (curPV->GetCopyNo()/10)%100;
232  depth = (curPV->GetCopyNo())%10 + 1;
233  if (depth > 0 && depth < 4 && layer >= 0 && layer < 17) {
234  edepHB += edeposit;
235  } else {
236  edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error "
237  << curPV->GetName() << curPV->GetCopyNo();
238  depth = -1; layer = -1;
239  }
240  } else if (name == "HES") {
241  layer = (curPV->GetCopyNo()/10)%100;
242  depth = (curPV->GetCopyNo())%10 + 1;
243  if (depth > 0 && depth < 3 && layer >= 0 && layer < 19) {
244  edepHE += edeposit;
245  } else {
246  edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error "
247  << curPV->GetName() << curPV->GetCopyNo();
248  depth = -1; layer = -1;
249  }
250  } else if (name == "HTS") {
251  layer = (curPV->GetCopyNo()/10)%100;
252  depth = (curPV->GetCopyNo())%10 + 1;
253  if (depth > 3 && depth < 5 && layer >= 17 && layer < 20) {
254  edepHO += edeposit;
255  } else {
256  edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error "
257  << curPV->GetName() <<curPV->GetCopyNo();
258  depth = -1; layer = -1;
259  }
260  }
261  if (depth >= 0 && depth < 5) edepd[depth] += edeposit;
262  if (layer >= 0 && layer < 20) edepl[layer] += edeposit;
263 
264  if (layer >= 0 && layer < 20) {
265  LogDebug("ValidHcal") << "SimG4HcalValidation:: G4Step: " << name
266  << " Layer " << std::setw(3) << layer << " Depth "
267  << std::setw(2) << depth << " Edep " <<std::setw(6)
268  << edeposit/MeV << " MeV";
269  }
270  }
271 }
272 
273 //================================================================ End of EVENT
274 void SimG4HcalValidation::update(const EndOfEvent * evt) {
275 
276  count++;
277 
278  // Fill hits cache for jetfinding etc.
279  fill(evt);
280  LogDebug("ValidHcal") << "SimG4HcalValidation:: --- after Fill ";
281 }
282 
283 //---------------------------------------------------
284 void SimG4HcalValidation::fill(const EndOfEvent * evt) {
285 
286  LogDebug("ValidHcal") << "SimG4HcalValidation:Fill event "
287  << (*evt)()->GetEventID();
288 
289  // access to the G4 hit collections
290  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
291 
292  int nhc = 0, j = 0;
293 
294  // Hcal
295  int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[0]);
296  CaloG4HitCollection* theHCHC = (CaloG4HitCollection*) allHC->GetHC(HCHCid);
297  LogDebug("ValidHcal") << "SimG4HcalValidation :: Hit Collection for "
298  << names[0] << " of ID " << HCHCid <<" is obtained at "
299  << theHCHC;
300  if (HCHCid >= 0 && theHCHC > 0) {
301  for (j = 0; j < theHCHC->entries(); j++) {
302 
303  CaloG4Hit* aHit = (*theHCHC)[j];
304 
305  double e = aHit->getEnergyDeposit()/GeV;
306  double time = aHit->getTimeSlice();
307 
308  math::XYZPoint pos = aHit->getPosition();
309  double theta = pos.theta();
310  double eta = -log(tan(theta * 0.5));
311  double phi = pos.phi();
312 
313  uint32_t unitID = aHit->getUnitID();
314  int subdet, zside, layer, etaIndex, phiIndex, lay;
315  org->unpackHcalIndex(unitID,subdet,zside,layer,etaIndex,phiIndex,lay);
316 
317  // some logic to separate HO ...
318  layer--;
319  std::string det = "HB";
320  if (subdet == static_cast<int>(HcalForward)) {
321  det = "HF";
322  uint16_t depth = aHit->getDepth();
323  if (depth != 0) { layer += 2; lay += 2; }
324  if (layer == 1) vhithc += e;
325  else if (layer == 0) vhitec += e;
326  else
327  edm::LogInfo("ValidHcal") << "SimG4HcalValidation::HitPMT "
328  << subdet << " " << (2*zside-1)*etaIndex
329  << " " << phiIndex << " " << layer+1
330  << " R " << pos.rho() << " Phi " << phi/deg
331  << " Edep " << e << " Time " << time;
332  } else if (subdet == static_cast<int>(HcalEndcap)) {
333  if (etaIndex <= 20) {
334  det = "HES";
335  } else {
336  det = "HED";
337  }
338  }
339  LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Hcal "
340  << det << " layer " << std::setw(2) << layer
341  << " lay " << std::setw(2) << lay << " time "
342  << std::setw(6) << time << " theta "
343  << std::setw(8) << theta << " eta " << std::setw(8)
344  << eta << " phi " << std::setw(8) << phi << " e "
345  << std::setw(8) << e << " ID 0x" << std::hex
346  << unitID << " ID dec " << std::dec << (int)unitID;
347 
348  // if desired, apply sampling factors in HCAL !!!
349  if (applySampling) e *= getHcalScale(det,layer);
350 
351  // filter on time & energy
352  if (time >= timeLowlim && time <= timeUplim && e > hhitThreshold ) {
353  enHcal += e;
354  CaloHit ahit(subdet,lay,e,eta,phi,time,unitID);
355  hitcache.push_back(ahit); // fill cache
356  ++nhc;
357  }
358  }
359  }
360  LogDebug("ValidHcal") << "SimG4HcalValidation:: HCAL hits : " << nhc;
361 
362  if (!hcalOnly) { //-------------------------- ECAL hits --------------------
363  int ndets = names.size();
364  if (ndets > 3) ndets = 3;
365  for (int idty = 1; idty < ndets; idty++) {
366  std::string det = "EB";
367  if (idty == 2) det = "EE";
368  else if (idty > 2) det = "ES";
369 
370  int nec = 0;
371  int ECHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[idty]);
372  CaloG4HitCollection* theECHC =(CaloG4HitCollection*)allHC->GetHC(ECHCid);
373  LogDebug("ValidHcal") << "SimG4HcalValidation:: Hit Collection for "
374  << names[idty] << " of ID " << ECHCid
375  << " is obtained at " << theECHC;
376  if (ECHCid >= 0 && theECHC > 0) {
377  for (j = 0; j < theECHC->entries(); j++) {
378 
379  CaloG4Hit* aHit = (*theECHC)[j];
380 
381  double e = aHit->getEnergyDeposit()/GeV;
382  double time = aHit->getTimeSlice();
383 
384  math::XYZPoint pos = aHit->getPosition();
385  double theta = pos.theta();
386  double eta = -log(tan(theta/2.));
387  double phi = pos.phi();
388  HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta,phi,1,1);
389  uint32_t unitID = org->getUnitID(id);
390  int subdet, zside, layer, ieta, iphi, lay;
391  org->unpackHcalIndex(unitID,subdet,zside,layer,ieta,iphi,lay);
392  subdet = idty+9;
393  layer = 0;
394  unitID = org->packHcalIndex(subdet,zside,layer,ieta,iphi,lay);
395 
396  // filter on time & energy
397  if (time >= timeLowlim && time <= timeUplim && e > ehitThreshold ) {
398  enEcal += e;
399  CaloHit ahit(subdet,lay,e,eta,phi,time,unitID);
400  hitcache.push_back(ahit); // fill cache
401  ++nec;
402  }
403 
404  LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Ecal "
405  << det << " layer " << std::setw(2) << layer
406  << " lay " << std::setw(2) << lay << " time "
407  << std::setw(6) << time << " theta "
408  << std::setw(8) << theta << " eta "
409  << std::setw(8) << eta << " phi "
410  << std::setw(8) << phi << " e " << std::setw(8)
411  << e << " ID 0x" << std::hex << unitID
412  << " ID dec " << std::dec << (int)unitID;
413 
414  }
415  }
416 
417  LogDebug("ValidHcal") << "SimG4HcalValidation:: " << det << " hits : "
418  << nec;
419  }
420  } // end of if(!hcalOnly)
421 
422 }
423 
424 void SimG4HcalValidation::layerAnalysis(PHcalValidInfoLayer& product) {
425 
426  int i = 0;
427  LogDebug("ValidHcal") << "\n ===>>> SimG4HcalValidation: Energy deposit "
428  << "in MeV\n at EB : " << std::setw(6) << edepEB/MeV
429  << "\n at EE : " << std::setw(6) << edepEE/MeV
430  << "\n at HB : " << std::setw(6) << edepHB/MeV
431  << "\n at HE : " << std::setw(6) << edepHE/MeV
432  << "\n at HO : " << std::setw(6) << edepHO/MeV
433  << "\n ---- SimG4HcalValidation: Energy deposit in";
434  for (i = 0; i < 5; i++)
435  LogDebug("ValidHcal") << " Depth " << std::setw(2) << i << " E "
436  << std::setw(8) << edepd[i]/MeV << " MeV";
437  LogDebug("ValidHcal") << " ---- SimG4HcalValidation: Energy deposit in"
438  << "layers";
439  for (i = 0; i < 20; i++)
440  LogDebug("ValidHcal") << " Layer " << std::setw(2) << i << " E "
441  << std::setw(8) << edepl[i]/MeV << " MeV";
442 
443  product.fillLayers(edepl, edepd, edepHO, edepHB+edepHE, edepEB+edepEE);
444 
445  // Hits in HF
446  product.fillHF(vhitec, vhithc, enEcal, enHcal);
447  LogDebug("ValidHcal") << "SimG4HcalValidation::HF hits " << vhitec
448  << " in EC and " << vhithc << " in HC\n"
449  << " HB/HE " << enEcal
450  << " in EC and " << enHcal << " in HC";
451 
452  // Another HCAL hist to porcess and store separately (a bit more complicated)
453  fetchHits(product);
454 
455  LogDebug("ValidHcal") << " layerAnalysis ===> after fetchHits";
456 
457 }
458 
459 //-----------------------------------------------------------------------------
460 void SimG4HcalValidation::nxNAnalysis(PHcalValidInfoNxN& product) {
461 
462  std::vector<CaloHit> * hits = &hitcache;
463  std::vector<CaloHit>::iterator hit_itr;
464 
465  LogDebug("ValidHcal") << "SimG4HcalValidation::NxNAnalysis : entrance ";
466 
467  collectEnergyRdir(eta0,phi0); // HCAL and ECAL energy in SimHitCache
468  // around (eta0,phi0)
469 
470  LogDebug("ValidHcal") << " NxNAnalysis : coolectEnergyRdir - Ecal " << een
471  << " Hcal " << hen;
472 
473  double etot = 0.; // total e deposited in "cluster"
474  double ee = 0.; // ECAL e deposited in "cluster"
475  double he = 0.; // HCAL e deposited in "cluster"
476  double hoe = 0.; // HO e deposited in "cluster"
477 
478  int max = dEta.size(); // 4
479 
480  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
481 
482  double e = hit_itr->e();
483  double t = hit_itr->t();
484  double eta = hit_itr->eta();
485  double phi = hit_itr->phi();
486  int type = hit_itr->det();
487  int layer = hit_itr->layer();
488 
489  // NxN calulation
490 
491  int save_i = 0;
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  save_i = i;
514  break;
515  }
516  }
517  // here comes break ...
518  }
519  }
520  }
521 
522  product.fillEcollectNxN(ee, he, hoe, etot);
523  product.fillHvsE(een, hen, hoen, een+hen);
524 
525  LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
526 
527 
528 }
529 
530 
531 //-----------------------------------------------------------------------------
532 void SimG4HcalValidation::jetAnalysis(PHcalValidInfoJets& product) {
533 
534  std::vector<CaloHit> * hhit = &hitcache;
535 
536  jetf->setInput(hhit);
537 
538  // zeroing cluster list, perfom clustering, fill cluster list & return pntr
539  std::vector<SimG4HcalHitCluster> * result = jetf->getClusters(hcalOnly);
540 
541  std::vector<SimG4HcalHitCluster>::iterator clus_itr;
542 
543  LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size()
544  << " clusters ---------------";
545  for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
546  LogDebug("ValidHcal") << (*clus_itr);
547 
548  std::vector<double> enevec, etavec, phivec;
549 
550  if ((*result).size() > 0) {
551 
552  sort((*result).begin(),(*result).end());
553 
554  clus_itr = result->begin(); // first cluster only
555  double etac = clus_itr->eta();
556  double phic = clus_itr->phi();
557 
558  double ecal_collect = 0.; // collect Ecal energy in the cone
559  if (!hcalOnly) {
560  ecal_collect = clus_itr->collectEcalEnergyR();}
561  else {
562  collectEnergyRdir(etac, phic);
563  ecal_collect = een;
564  }
565  LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect = "
566  << ecal_collect;
567 
568  // eta-phi deviation of the cluster from nominal (eta0, phi0) values
569  double dist = jetf->rDist(eta0, phi0, etac, phic);
570  LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation = " << dist;
571  product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
572 
573  LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
574 
575  std::vector<CaloHit> * hits = clus_itr->getHits() ;
576  std::vector<CaloHit>::iterator hit_itr;
577 
578  double ee = 0., he = 0., hoe = 0., etot = 0.;
579 
580  // cycle over all hits in the FIRST cluster
581  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
582  double e = hit_itr->e();
583  double t = hit_itr->t();
584  double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
585 
586  // energy collection
587  etot += e;
588  if (hit_itr->det() == 10 || hit_itr->det() == 11 ||
589  hit_itr->det() == 12) ee += e;
590  if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
591  hit_itr->det() == static_cast<int>(HcalEndcap) ||
592  hit_itr->det() == static_cast<int>(HcalForward)) {
593  he += e;
594  if (hit_itr->det() == static_cast<int>(HcalBarrel) &&
595  hit_itr->layer() > 17)
596  hoe += e;
597  }
598 
599  if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
600  hit_itr->det() == static_cast<int>(HcalEndcap) ||
601  hit_itr->det() == static_cast<int>(HcalForward)) {
602  product.fillTProfileJet(he, r, t);
603  }
604  }
605 
606  product.fillEcollectJet(ee, he, hoe, etot);
607 
608  LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
609  << "ee/he/hoe/etot " << ee << "/" << he << "/"
610  << hoe << "/" << etot;
611 
612  // Loop over clusters
613  for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
614  if (clus_itr->e() > jetThreshold) {
615  enevec.push_back(clus_itr->e());
616  etavec.push_back(clus_itr->eta());
617  phivec.push_back(clus_itr->phi());
618  }
619  }
620  product.fillJets(enevec, etavec, phivec);
621 
622  LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
623  << " JetAnalysis ===> (*result).size() "
624  << (*result).size();
625 
626  // Di-jet mass
627  if (etavec.size() > 1) {
628  if (etavec[0] > -2.5 && etavec[0] < 2.5 &&
629  etavec[1] > -2.5 && etavec[1] < 2.5) {
630 
631  LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
632  << " JetAnalysis ===> Di-jet vectors ";
633  for (unsigned int i = 0; i < enevec.size(); i++)
634  LogDebug("ValidHcal") << " e, eta, phi = " << enevec[i] << " "
635  << etavec[i] << " " << phivec[i];
636 
637  double et0 = enevec[0] / cosh(etavec[0]);
638  double px0 = et0 * cos(phivec[0]);
639  double py0 = et0 * sin(phivec[0]);
640  double pz0 = et0 * sinh(etavec[0]);
641  double et1 = enevec[1] / cosh(etavec[1]);
642  double px1 = et1 * cos(phivec[1]);
643  double py1 = et1 * sin(phivec[1]);
644  double pz1 = et1 * sinh(etavec[1]);
645 
646  double dijetmass2 = (enevec[0]+enevec[1])*(enevec[0]+enevec[1])
647  - (px1+px0)*(px1+px0) - (py1+py0)*(py1+py0) - (pz1+pz0)*(pz1+pz0);
648 
649  LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ "
650  << dijetmass2;
651 
652  double dijetmass;
653  if(dijetmass2 >= 0.) dijetmass = sqrt(dijetmass2);
654  else dijetmass = -sqrt(-1. * dijetmass2);
655 
656  product.fillDiJets(dijetmass);
657 
658  LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
659  }
660  }
661  }
662 }
663 
664 //---------------------------------------------------
665 void SimG4HcalValidation::fetchHits(PHcalValidInfoLayer& product) {
666 
667  LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with "
668  << hitcache.size() << " hits";
669  int nHit = hitcache.size();
670  int hit = 0;
671  int i;
672  std::vector<CaloHit>::iterator itr;
673  std::vector<CaloHit*> lhits(nHit);
674  for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
675  uint32_t unitID=itr->id();
676  int subdet, zside, group, ieta, iphi, lay;
677  HcalTestNumbering::unpackHcalIndex(unitID,subdet,zside,group,
678  ieta,iphi,lay);
679  subdet = itr->det();
680  lay = itr->layer();
681  group = (subdet&15)<<20;
682  group += ((lay-1)&31)<<15;
683  group += (zside&1)<<14;
684  group += (ieta&127)<<7;
685  group += (iphi&127);
686  itr->setId(group);
687  lhits[i] = &hitcache[i];
688  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i
689  << " " << hitcache[i] << "\n"
690  << "SimG4HcalValidation::fetchHits:Copied " << i
691  << " " << *lhits[i];
692  }
693  sort(lhits.begin(),lhits.end(),CaloHitIdMore());
694  std::vector<CaloHit*>::iterator k1, k2;
695  for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
696  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i
697  << " " << **k1;
698  int nHits = 0;
699  for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
700  double ehit = (**k1).e();
701  double t = (**k1).t();
702  uint32_t unitID= (**k1).id();
703  int jump = 0;
704  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i
705  << " U/T/E" << " 0x" << std::hex << unitID
706  << std::dec << " " << t << " " << ehit;
707  for (k2 = k1+1; k2 != lhits.end() && (t-(**k2).t()) < 1 &&
708  (t-(**k2).t()) > -1 && unitID == (**k2).id(); k2++) {
709  ehit += (**k2).e();
710  LogDebug("ValidHcal") << "\t + " << (**k2).e();
711  jump++;
712  }
713  LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
714 
715  double eta = (*k1)->eta();
716  double phi = (*k1)->phi();
717  int lay = ((unitID>>15)&31) + 1;
718  int subdet = (unitID>>20)&15;
719  int zside = (unitID>>14)&1;
720  int ieta = (unitID>>7)&127;
721  int iphi = (unitID)&127;
722 
723  // All hits in cache
724  product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
725  nHits++;
726 
727  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits
728  << " " << i << " ID 0x" << std::hex << unitID
729  << " det " << std::dec << subdet << " " << lay
730  << " " << zside << " " << ieta << " " << iphi
731  << " Time " << t << " E " << ehit;
732 
733  i += jump;
734  k1 += jump;
735  }
736 
737  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with "
738  << nHit << " hits" << " and writes out " << nHits
739  << '(' << hit << ") hits";
740 
741 }
742 //---------------------------------------------------
743 void SimG4HcalValidation::clear(){
744  hitcache.erase( hitcache.begin(), hitcache.end());
745 }
746 
747 //---------------------------------------------------
748 void SimG4HcalValidation::collectEnergyRdir (const double eta0,
749  const double phi0) {
750 
751  std::vector<CaloHit> * hits = &hitcache;
752  std::vector<CaloHit>::iterator hit_itr;
753 
754  double sume = 0., sumh = 0., sumho = 0.;
755 
756  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
757 
758  double e = hit_itr->e();
759  double eta = hit_itr->eta();
760  double phi = hit_itr->phi();
761 
762  int type = hit_itr->det();
763 
764  double r = jetf->rDist(eta0, phi0, eta, phi);
765  if (r < coneSize) {
766  if (type == 10 || type == 11 || type == 12) {
767  sume += e;
768  } else {
769  sumh += e;
770  if (type == static_cast<int>(HcalBarrel) &&
771  hit_itr->layer() > 17) sumho += e;
772  }
773  }
774  }
775 
776  een = sume;
777  hen = sumh;
778  hoen = sumho;
779 }
780 
781 //---------------------------------------------------
782 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
783 
784  double tmp = 0.;
785 
786  if (det == "HB") {
787  tmp = scaleHB[layer];
788  } else if (det == "HES" || det == "HED") {
789  tmp = scaleHE[layer];
790  } else if (det == "HF") {
791  tmp = scaleHF[layer];
792  }
793 
794  return tmp;
795 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:516
type
type
Definition: HCALResponse.h:22
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
SimG4HcalValidation::numberingFromDDD
HcalNumberingFromDDD * numberingFromDDD
Definition: SimG4HcalValidation.h:72
i
int i
Definition: DBlmapReader.cc:9
CaloG4Hit::getPosition
math::XYZPoint getPosition() const
Definition: CaloG4Hit.h:56
PHcalValidInfoNxN::fillHvsE
void fillHvsE(double ee, double he, double hoe, double etot)
Definition: PValidationFormats.cc:1452
PHcalValidInfoJets::fillTProfileJet
void fillTProfileJet(double e, double r, double t)
Definition: PValidationFormats.cc:1359
SimG4HcalValidation::getHcalScale
double getHcalScale(std::string, int) const
Definition: SimG4HcalValidation.cc:782
SimG4HcalValidation::fill
void fill(const EndOfEvent *ev)
Definition: SimG4HcalValidation.cc:284
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:1409
SimG4HcalValidation::scaleHF
std::vector< float > scaleHF
Definition: SimG4HcalValidation.h:83
PHcalValidInfoJets::fillEcollectJet
void fillEcollectJet(double ee, double he, double hoe, double etot)
Definition: PValidationFormats.cc:1369
SimG4HcalValidation::collectEnergyRdir
void collectEnergyRdir(const double, const double)
Definition: SimG4HcalValidation.cc:748
SimG4HcalValidation::names
std::vector< std::string > names
Definition: SimG4HcalValidation.h:86
PHcalValidInfoJets::fillEtaPhiProfileJet
void fillEtaPhiProfileJet(double eta0, double phi0, double eta, double phi, double dist)
Definition: PValidationFormats.cc:1379
SimG4HcalValidation::hitcache
std::vector< CaloHit > hitcache
Definition: SimG4HcalValidation.h:78
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PHcalValidInfoJets::fillDiJets
void fillDiJets(double mass)
Definition: PValidationFormats.cc:1401
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:161
CaloG4Hit::getDepth
uint16_t getDepth() const
Definition: CaloG4Hit.h:72
SimG4HcalValidation::dEta
std::vector< double > dEta
Definition: SimG4HcalValidation.h:94
SimG4HcalHitJetFinder::setInput
void setInput(std::vector< CaloHit > *)
Definition: SimG4HcalHitJetFinder.cc:23
NULL
#define NULL
Definition: scimark2.h:8
eta
T eta() const
Definition: Basic3DVectorLD.h:168
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:488
PHcalValidInfoLayer::fillHits
void fillHits(int Nhits, int lay, int unitID, double eta, double phi, double ehit, double t)
Definition: PValidationFormats.cc:1433
HCalSD
Definition: HCalSD.h:31
SimG4HcalValidation::jetAnalysis
void jetAnalysis(PHcalValidInfoJets &)
Definition: SimG4HcalValidation.cc:532
SimG4HcalValidation::nxNAnalysis
void nxNAnalysis(PHcalValidInfoNxN &)
Definition: SimG4HcalValidation.cc:460
SimG4HcalValidation::jetf
SimG4HcalHitJetFinder * jetf
Definition: SimG4HcalValidation.h:69
SimG4HcalValidation::dPhi
std::vector< double > dPhi
Definition: SimG4HcalValidation.h:95
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
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:665
edm::Event::put
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
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:81
PHcalValidInfoNxN::fillEcollectNxN
void fillEcollectNxN(double een, double hen, double hoen, double etotn)
Definition: PValidationFormats.cc:1460
SimG4HcalHitJetFinder::getClusters
std::vector< SimG4HcalHitCluster > * getClusters(bool)
Definition: SimG4HcalHitJetFinder.cc:27
funct::log
Log< T >::type log(const T &t)
Definition: Log.h:22
PHcalValidInfoNxN::fillTProfileNxN
void fillTProfileNxN(double e, int i, double t)
Definition: PValidationFormats.cc:1468
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
PHcalValidInfoJets::fillJets
void fillJets(std::vector< double > enj, std::vector< double > etaj, std::vector< double > phij)
Definition: PValidationFormats.cc:1387
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
SimG4HcalValidation::scaleHE
std::vector< float > scaleHE
Definition: SimG4HcalValidation.h:82
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
CaloG4HitCollection
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
Definition: CaloG4HitCollection.h:11
launcher.job
tuple job
Definition: launcher.py:57
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
HcalNumberingFromDDD::unitID
HcalID unitID(int det, CLHEP::Hep3Vector pos, int depth, int lay=-1) const
Definition: HcalNumberingFromDDD.cc:32
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:69
SimG4HcalValidation::org
HcalTestNumberingScheme * org
Definition: SimG4HcalValidation.h:75
SimG4HcalValidation::layerAnalysis
void layerAnalysis(PHcalValidInfoLayer &)
Definition: SimG4HcalValidation.cc:424
HcalTestNumberingScheme::getUnitID
virtual uint32_t getUnitID(const HcalNumberingFromDDD::HcalID id)
Definition: HcalTestNumberingScheme.cc:21
PHcalValidInfoLayer::fillHF
void fillHF(double fibl, double fibs, double enec, double enhc)
Definition: PValidationFormats.cc:1425
SimG4HcalValidation::produce
void produce(edm::Event &, const edm::EventSetup &)
Definition: SimG4HcalValidation.cc:99
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5