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