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  if (fabs(eta0-eta) <= dEta[max-1] && fabs(phi0-phi) <= dPhi[max-1]) {
492  etot += e;
493  if (type == 10 || type == 11 || type == 12) ee += e;
494  if (type == static_cast<int>(HcalBarrel) ||
495  type == static_cast<int>(HcalEndcap) ||
496  type == static_cast<int>(HcalForward)) {
497  he += e;
498  if (type == static_cast<int>(HcalBarrel) && layer > 17) hoe += e;
499 
500  // which concrete i-th square ?
501  for (int i = 0; i < max; i++ ) {
502  if ((eta0-eta) <= dEta[i] && (eta0-eta) >= -dEta[i] &&
503  (phi0-phi) <= dPhi[i] && (phi0-phi) >= -dPhi[i]) {
504 
505  LogDebug("ValidHcal") << "SimG4HcalValidation:: nxNAnalysis eta0,"
506  << " phi0 = " << eta0 << " " << phi0
507  << " type, layer = " << type << ","
508  << layer << " eta, phi = " << eta << " "
509  << phi;
510 
511  product.fillTProfileNxN(e, i, t);
512  break;
513  }
514  }
515  // here comes break ...
516  }
517  }
518  }
519 
520  product.fillEcollectNxN(ee, he, hoe, etot);
521  product.fillHvsE(een, hen, hoen, een+hen);
522 
523  LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
524 
525 
526 }
527 
528 
529 //-----------------------------------------------------------------------------
530 void SimG4HcalValidation::jetAnalysis(PHcalValidInfoJets& product) {
531 
532  std::vector<CaloHit> * hhit = &hitcache;
533 
534  jetf->setInput(hhit);
535 
536  // zeroing cluster list, perfom clustering, fill cluster list & return pntr
537  std::vector<SimG4HcalHitCluster> * result = jetf->getClusters(hcalOnly);
538 
539  std::vector<SimG4HcalHitCluster>::iterator clus_itr;
540 
541  LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size()
542  << " clusters ---------------";
543  for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
544  LogDebug("ValidHcal") << (*clus_itr);
545 
546  std::vector<double> enevec, etavec, phivec;
547 
548  if ((*result).size() > 0) {
549 
550  sort((*result).begin(),(*result).end());
551 
552  clus_itr = result->begin(); // first cluster only
553  double etac = clus_itr->eta();
554  double phic = clus_itr->phi();
555 
556  double ecal_collect = 0.; // collect Ecal energy in the cone
557  if (!hcalOnly) {
558  ecal_collect = clus_itr->collectEcalEnergyR();}
559  else {
560  collectEnergyRdir(etac, phic);
561  ecal_collect = een;
562  }
563  LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect = "
564  << ecal_collect;
565 
566  // eta-phi deviation of the cluster from nominal (eta0, phi0) values
567  double dist = jetf->rDist(eta0, phi0, etac, phic);
568  LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation = " << dist;
569  product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
570 
571  LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
572 
573  std::vector<CaloHit> * hits = clus_itr->getHits() ;
574  std::vector<CaloHit>::iterator hit_itr;
575 
576  double ee = 0., he = 0., hoe = 0., etot = 0.;
577 
578  // cycle over all hits in the FIRST cluster
579  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
580  double e = hit_itr->e();
581  double t = hit_itr->t();
582  double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
583 
584  // energy collection
585  etot += e;
586  if (hit_itr->det() == 10 || hit_itr->det() == 11 ||
587  hit_itr->det() == 12) ee += e;
588  if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
589  hit_itr->det() == static_cast<int>(HcalEndcap) ||
590  hit_itr->det() == static_cast<int>(HcalForward)) {
591  he += e;
592  if (hit_itr->det() == static_cast<int>(HcalBarrel) &&
593  hit_itr->layer() > 17)
594  hoe += e;
595  }
596 
597  if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
598  hit_itr->det() == static_cast<int>(HcalEndcap) ||
599  hit_itr->det() == static_cast<int>(HcalForward)) {
600  product.fillTProfileJet(he, r, t);
601  }
602  }
603 
604  product.fillEcollectJet(ee, he, hoe, etot);
605 
606  LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
607  << "ee/he/hoe/etot " << ee << "/" << he << "/"
608  << hoe << "/" << etot;
609 
610  // Loop over clusters
611  for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
612  if (clus_itr->e() > jetThreshold) {
613  enevec.push_back(clus_itr->e());
614  etavec.push_back(clus_itr->eta());
615  phivec.push_back(clus_itr->phi());
616  }
617  }
618  product.fillJets(enevec, etavec, phivec);
619 
620  LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
621  << " JetAnalysis ===> (*result).size() "
622  << (*result).size();
623 
624  // Di-jet mass
625  if (etavec.size() > 1) {
626  if (etavec[0] > -2.5 && etavec[0] < 2.5 &&
627  etavec[1] > -2.5 && etavec[1] < 2.5) {
628 
629  LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
630  << " JetAnalysis ===> Di-jet vectors ";
631  for (unsigned int i = 0; i < enevec.size(); i++)
632  LogDebug("ValidHcal") << " e, eta, phi = " << enevec[i] << " "
633  << etavec[i] << " " << phivec[i];
634 
635  double et0 = enevec[0] / cosh(etavec[0]);
636  double px0 = et0 * cos(phivec[0]);
637  double py0 = et0 * sin(phivec[0]);
638  double pz0 = et0 * sinh(etavec[0]);
639  double et1 = enevec[1] / cosh(etavec[1]);
640  double px1 = et1 * cos(phivec[1]);
641  double py1 = et1 * sin(phivec[1]);
642  double pz1 = et1 * sinh(etavec[1]);
643 
644  double dijetmass2 = (enevec[0]+enevec[1])*(enevec[0]+enevec[1])
645  - (px1+px0)*(px1+px0) - (py1+py0)*(py1+py0) - (pz1+pz0)*(pz1+pz0);
646 
647  LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ "
648  << dijetmass2;
649 
650  double dijetmass;
651  if(dijetmass2 >= 0.) dijetmass = sqrt(dijetmass2);
652  else dijetmass = -sqrt(-1. * dijetmass2);
653 
654  product.fillDiJets(dijetmass);
655 
656  LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
657  }
658  }
659  }
660 }
661 
662 //---------------------------------------------------
663 void SimG4HcalValidation::fetchHits(PHcalValidInfoLayer& product) {
664 
665  LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with "
666  << hitcache.size() << " hits";
667  int nHit = hitcache.size();
668  int hit = 0;
669  int i;
670  std::vector<CaloHit>::iterator itr;
671  std::vector<CaloHit*> lhits(nHit);
672  for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
673  uint32_t unitID=itr->id();
674  int subdet, zside, group, ieta, iphi, lay;
675  HcalTestNumbering::unpackHcalIndex(unitID,subdet,zside,group,
676  ieta,iphi,lay);
677  subdet = itr->det();
678  lay = itr->layer();
679  group = (subdet&15)<<20;
680  group += ((lay-1)&31)<<15;
681  group += (zside&1)<<14;
682  group += (ieta&127)<<7;
683  group += (iphi&127);
684  itr->setId(group);
685  lhits[i] = &hitcache[i];
686  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i
687  << " " << hitcache[i] << "\n"
688  << "SimG4HcalValidation::fetchHits:Copied " << i
689  << " " << *lhits[i];
690  }
691  sort(lhits.begin(),lhits.end(),CaloHitIdMore());
692  std::vector<CaloHit*>::iterator k1, k2;
693  for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
694  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i
695  << " " << **k1;
696  int nHits = 0;
697  for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
698  double ehit = (**k1).e();
699  double t = (**k1).t();
700  uint32_t unitID= (**k1).id();
701  int jump = 0;
702  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i
703  << " U/T/E" << " 0x" << std::hex << unitID
704  << std::dec << " " << t << " " << ehit;
705  for (k2 = k1+1; k2 != lhits.end() && (t-(**k2).t()) < 1 &&
706  (t-(**k2).t()) > -1 && unitID == (**k2).id(); k2++) {
707  ehit += (**k2).e();
708  LogDebug("ValidHcal") << "\t + " << (**k2).e();
709  jump++;
710  }
711  LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
712 
713  double eta = (*k1)->eta();
714  double phi = (*k1)->phi();
715  int lay = ((unitID>>15)&31) + 1;
716  int subdet = (unitID>>20)&15;
717  int zside = (unitID>>14)&1;
718  int ieta = (unitID>>7)&127;
719  int iphi = (unitID)&127;
720 
721  // All hits in cache
722  product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
723  nHits++;
724 
725  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits
726  << " " << i << " ID 0x" << std::hex << unitID
727  << " det " << std::dec << subdet << " " << lay
728  << " " << zside << " " << ieta << " " << iphi
729  << " Time " << t << " E " << ehit;
730 
731  i += jump;
732  k1 += jump;
733  }
734 
735  LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with "
736  << nHit << " hits" << " and writes out " << nHits
737  << '(' << hit << ") hits";
738 
739 }
740 //---------------------------------------------------
741 void SimG4HcalValidation::clear(){
742  hitcache.erase( hitcache.begin(), hitcache.end());
743 }
744 
745 //---------------------------------------------------
746 void SimG4HcalValidation::collectEnergyRdir (const double eta0,
747  const double phi0) {
748 
749  std::vector<CaloHit> * hits = &hitcache;
750  std::vector<CaloHit>::iterator hit_itr;
751 
752  double sume = 0., sumh = 0., sumho = 0.;
753 
754  for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
755 
756  double e = hit_itr->e();
757  double eta = hit_itr->eta();
758  double phi = hit_itr->phi();
759 
760  int type = hit_itr->det();
761 
762  double r = jetf->rDist(eta0, phi0, eta, phi);
763  if (r < coneSize) {
764  if (type == 10 || type == 11 || type == 12) {
765  sume += e;
766  } else {
767  sumh += e;
768  if (type == static_cast<int>(HcalBarrel) &&
769  hit_itr->layer() > 17) sumho += e;
770  }
771  }
772  }
773 
774  een = sume;
775  hen = sumh;
776  hoen = sumho;
777 }
778 
779 //---------------------------------------------------
780 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
781 
782  double tmp = 0.;
783 
784  if (det == "HB") {
785  tmp = scaleHB[layer];
786  } else if (det == "HES" || det == "HED") {
787  tmp = scaleHE[layer];
788  } else if (det == "HF") {
789  tmp = scaleHF[layer];
790  }
791 
792  return tmp;
793 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
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:780
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:746
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
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1108
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:165
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:172
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:484
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:530
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:663
edm::Event::put
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:85
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
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
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
sd
double sd
Definition: CascadeWrapper.h:137
PHcalValidInfoJets::fillJets
void fillJets(std::vector< double > enj, std::vector< double > etaj, std::vector< double > phij)
Definition: PValidationFormats.cc:1387
SimG4HcalValidation::scaleHE
std::vector< float > scaleHE
Definition: SimG4HcalValidation.h:82
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
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:31
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