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/ESTransientHandle.h"
16 #include "Geometry/Records/interface/IdealGeometryRecord.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 "CLHEP/Random/Random.h"
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::ESTransientHandle<DDCompactView> pDD;
155  (*job)()->get<IdealGeometryRecord>().get(pDD);
156  edm::LogInfo("HcalSim") << "HcalTestAnalysis:: Initialise "
157  << "HcalNumberingFromDDD for " << names[0];
158  numberingFromDDD = new HcalNumberingFromDDD(names[0], (*pDD));
159 
160  // Ntuples
161  tuplesManager.reset(new HcalTestHistoManager(fileName));
162 
163  // Numbering scheme
164  org = new HcalTestNumberingScheme(false);
165 
166 }
167 
168 //==================================================================== per RUN
169 void HcalTestAnalysis::update(const BeginOfRun * run) {
170 
171  int irun = (*run)()->GetRunID();
172  edm::LogInfo("HcalSim") << "HcalTestAnalysis:: Begin of Run = " << irun;
173 
174  bool loop = true, eta = true, phi = true;
175  int etac = (centralTower/100)%100;
176  if (etac == 0) {
177  etac = 1;
178  eta = false;
179  }
180  int phic = (centralTower%100);
181  if (phic == 0) {
182  phic = 1;
183  phi = false;
184  }
185  int idet = static_cast<int>(HcalBarrel);
186  while (loop) {
187  HcalCellType::HcalCell tmp = numberingFromDDD->cell(idet,1,1,etac,phic);
188  if (tmp.ok) {
189  if (eta) eta0 = tmp.eta;
190  if (phi) phi0 = tmp.phi;
191  loop = false;
192  } else if (idet == static_cast<int>(HcalBarrel)) {
193  idet = static_cast<int>(HcalEndcap);
194  } else if (idet == static_cast<int>(HcalEndcap)) {
195  idet = static_cast<int>(HcalForward);
196  } else {
197  loop = false;
198  }
199  }
200 
201  edm::LogInfo("HcalSim") << "HcalTestAnalysis:: Central Tower "
202  << centralTower << " corresponds to eta0 = " << eta0
203  << " phi0 = " << phi0;
204 
205  std::string sdname = names[0];
206  G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
207  if (sd != 0) {
208  G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
209  if (aSD==0) {
210  edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: No SD with "
211  << "name " << sdname << " in this Setup";
212  } else {
213  HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
214  edm::LogInfo("HcalSim") << "HcalTestAnalysis::beginOfRun: Finds SD with "
215  << "name " << theCaloSD->GetName()
216  << " in this Setup";
217  if (org) {
218  theCaloSD->setNumberingScheme(org);
219  edm::LogInfo("HcalSim") << "HcalTestAnalysis::beginOfRun: set a new "
220  << "numbering scheme";
221  }
222  }
223  } else {
224  edm::LogWarning("HcalSim") << "HcalTestAnalysis::beginOfRun: Could not get"
225  << " SD Manager!";
226  }
227 
228 }
229 
230 //=================================================================== per EVENT
231 void HcalTestAnalysis::update(const BeginOfEvent * evt) {
232 
233  // create tuple object
234  tuples = new HcalTestHistoClass();
235  // Reset counters
236  tuples->setCounters();
237 
238  int i = 0;
239  edepEB = edepEE = edepHB = edepHE = edepHO = 0.;
240  for (i = 0; i < 20; i++) edepl[i] = 0.;
241  for (i = 0; i < 20; i++) mudist[i] = -1.;
242 
243  int iev = (*evt)()->GetEventID();
244  LogDebug("HcalSim") <<"HcalTestAnalysis: Begin of event = " << iev;
245 }
246 
247 //=================================================================== each STEP
248 void HcalTestAnalysis::update(const G4Step * aStep) {
249 
250  if (aStep != NULL) {
251  G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
252  G4String name = curPV->GetName();
253  name.assign(name,0,3);
254  double edeposit = aStep->GetTotalEnergyDeposit();
255  int layer=-1;
256  if (name == "EBR") {
257  edepEB += edeposit;
258  } else if (name == "EFR") {
259  edepEE += edeposit;
260  } else if (name == "HBS") {
261  layer = (curPV->GetCopyNo()/10)%100;
262  if (layer >= 0 && layer < 17) {
263  edepHB += edeposit;
264  } else {
265  edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HB "
266  << curPV->GetName() << curPV->GetCopyNo();
267  layer = -1;
268  }
269  } else if (name == "HES") {
270  layer = (curPV->GetCopyNo()/10)%100;
271  if (layer >= 0 && layer < 19) {
272  edepHE += edeposit;
273  } else {
274  edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HE "
275  << curPV->GetName() << curPV->GetCopyNo();
276  layer = -1;
277  }
278  } else if (name == "HTS") {
279  layer = (curPV->GetCopyNo()/10)%100;
280  if (layer >= 17 && layer < 20) {
281  edepHO += edeposit;
282  } else {
283  edm::LogWarning("HcalSim") << "HcalTestAnalysis::Error in HO "
284  << curPV->GetName() << curPV->GetCopyNo();
285  layer = -1;
286  }
287  }
288  if (layer >= 0 && layer < 20) {
289  edepl[layer] += edeposit;
290 
291  // Calculate the distance if it is a muon
292  G4String part = aStep->GetTrack()->GetDefinition()->GetParticleName();
293  if ((part == "mu-" || part == "mu+") && mudist[layer] < 0) {
294  math::XYZPoint pos(aStep->GetPreStepPoint()->GetPosition().x(),
295  aStep->GetPreStepPoint()->GetPosition().y(),
296  aStep->GetPreStepPoint()->GetPosition().z());
297  double theta = pos.theta();
298  double eta = -log(tan(theta * 0.5));
299  double phi = pos.phi();
300  double dist = sqrt ((eta-eta0)*(eta-eta0) + (phi-phi0)*(phi-phi0));
301  mudist[layer] = dist*std::sqrt(pos.perp2());
302  }
303  }
304 
305  if (layer >= 0 && layer < 20) {
306  LogDebug("HcalSim") << "HcalTestAnalysis:: G4Step: " << name << " Layer "
307  << std::setw(3) << layer << " Edep " << std::setw(6)
308  << edeposit/MeV << " MeV";
309  }
310  } else {
311  edm::LogInfo("HcalSim") << "HcalTestAnalysis:: G4Step: Null Step";
312  }
313 }
314 
315 //================================================================ End of EVENT
316 void HcalTestAnalysis::update(const EndOfEvent * evt) {
317 
318  count++;
319  // Fill event input
320  fill(evt);
321  LogDebug("HcalSim") << "HcalTestAnalysis:: --- after Fill";
322 
323  // Qie analysis
324  CLHEP::HepRandomEngine* engine = CLHEP::HepRandom::getTheEngine();
325  qieAnalysis(engine);
326  LogDebug("HcalSim") << "HcalTestAnalysis:: --- after QieAnalysis";
327 
328  // Layers tuples filling
329  layerAnalysis();
330  LogDebug("HcalSim") << "HcalTestAnalysis:: --- after LayerAnalysis";
331 
332  // Writing the data to the Tree
333  tuplesManager->fillTree(tuples); // (no need to delete it...)
334  tuples = 0; // but avoid to reuse it...
335  LogDebug("HcalSim") << "HcalTestAnalysis:: --- after fillTree";
336 
337 }
338 
339 //---------------------------------------------------
340 void HcalTestAnalysis::fill(const EndOfEvent * evt) {
341 
342  LogDebug("HcalSim") << "HcalTestAnalysis: Fill event "
343  << (*evt)()->GetEventID();
344 
345  // access to the G4 hit collections
346  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
347 
348  int nhc = 0, neb = 0, nef = 0, j = 0;
349  caloHitCache.erase (caloHitCache.begin(), caloHitCache.end());
350 
351  // Hcal
352  int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[0]);
353  CaloG4HitCollection* theHCHC = (CaloG4HitCollection*) allHC->GetHC(HCHCid);
354  LogDebug("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names[0]
355  << " of ID " << HCHCid << " is obtained at " << theHCHC;
356  if (HCHCid >= 0 && theHCHC > 0) {
357  for (j = 0; j < theHCHC->entries(); j++) {
358 
359  CaloG4Hit* aHit = (*theHCHC)[j];
360 
361  double e = aHit->getEnergyDeposit()/GeV;
362  double time = aHit->getTimeSlice();
363 
364  math::XYZPoint pos = aHit->getPosition();
365  double theta = pos.theta();
366  double eta = -log(tan(theta * 0.5));
367  double phi = pos.phi();
368 
369  uint32_t unitID = aHit->getUnitID();
370  int subdet, zside, layer, etaIndex, phiIndex, lay;
371  org->unpackHcalIndex(unitID,subdet,zside,layer,etaIndex,phiIndex,lay);
372  double jitter = time-timeOfFlight(subdet,lay,eta);
373  if (jitter<0) jitter = 0;
374  CaloHit hit(subdet,lay,e,eta,phi,jitter,unitID);
375  caloHitCache.push_back(hit);
376  nhc++;
377 
378  std::string det = "HB";
379  if (subdet == static_cast<int>(HcalForward)) {
380  det = "HF";
381  } else if (subdet == static_cast<int>(HcalEndcap)) {
382  if (etaIndex <= 20) {
383  det = "HES";
384  } else {
385  det = "HED";
386  }
387  }
388 
389  LogDebug("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2)
390  << layer << " time " << std::setw(6) << time
391  << " theta " << std::setw(8) << theta << " eta "
392  << std::setw(8) << eta << " phi " << std::setw(8)
393  << phi << " e " << std::setw(8) << e;
394  }
395  }
396  LogDebug("HcalSim") << "HcalTestAnalysis::HCAL hits : " << nhc << "\n";
397 
398  // EB
399  int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[1]);
400  CaloG4HitCollection* theEBHC = (CaloG4HitCollection*) allHC->GetHC(EBHCid);
401  LogDebug("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names[1]
402  << " of ID " << EBHCid << " is obtained at " << theEBHC;
403  if (EBHCid >= 0 && theEBHC > 0) {
404  for (j = 0; j < theEBHC->entries(); j++) {
405 
406  CaloG4Hit* aHit = (*theEBHC)[j];
407 
408  double e = aHit->getEnergyDeposit()/GeV;
409  double time = aHit->getTimeSlice();
410  std::string det = "EB";
411 
412  math::XYZPoint pos = aHit->getPosition();
413  double theta = pos.theta();
414  double eta = -log(tan(theta/2.));
415  double phi = pos.phi();
416 
417  HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta,phi,1,1);
418  uint32_t unitID = org->getUnitID(id);
419  int subdet, zside, layer, ieta, iphi, lay;
420  org->unpackHcalIndex(unitID,subdet,zside,layer,ieta,iphi,lay);
421  subdet = 10;
422  layer = 0;
423  unitID = org->packHcalIndex(subdet,zside,layer,ieta,iphi,lay);
424  CaloHit hit(subdet,lay,e,eta,phi,time,unitID);
425  caloHitCache.push_back(hit);
426  neb++;
427  LogDebug("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2)
428  << layer << " time " << std::setw(6) << time
429  << " theta " << std::setw(8) << theta << " eta "
430  << std::setw(8) << eta << " phi " << std::setw(8)
431  << phi << " e " << std::setw(8) << e;
432  }
433  }
434  LogDebug("HcalSim") << "HcalTestAnalysis::EB hits : " << neb << "\n";
435 
436  // EE
437  int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[2]);
438  CaloG4HitCollection* theEEHC = (CaloG4HitCollection*) allHC->GetHC(EEHCid);
439  LogDebug("HcalSim") << "HcalTestAnalysis :: Hit Collection for " << names[2]
440  << " of ID " << EEHCid << " is obtained at " << theEEHC;
441  if (EEHCid >= 0 && theEEHC > 0) {
442  for (j = 0; j < theEEHC->entries(); j++) {
443 
444  CaloG4Hit* aHit = (*theEEHC)[j];
445 
446  double e = aHit->getEnergyDeposit()/GeV;
447  double time = aHit->getTimeSlice();
448  std::string det = "EE";
449 
450  math::XYZPoint pos = aHit->getPosition();
451  double theta = pos.theta();
452  double eta = -log(tan(theta/2.));
453  double phi = pos.phi();
454 
455  HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta,phi,1,1);
456  uint32_t unitID = org->getUnitID(id);
457  int subdet, zside, layer, ieta, iphi, lay;
458  org->unpackHcalIndex(unitID,subdet,zside,layer,ieta,iphi,lay);
459  subdet = 11;
460  layer = 0;
461  unitID = org->packHcalIndex(subdet,zside,layer,ieta,iphi,lay);
462  CaloHit hit(subdet,lay,e,eta,phi,time,unitID);
463  caloHitCache.push_back(hit);
464  nef++;
465  LogDebug("HcalSim") << "HcalTest: " << det << " layer " << std::setw(2)
466  << layer << " time " << std::setw(6) << time
467  << " theta " << std::setw(8) << theta << " eta "
468  << std::setw(8) << eta << " phi " << std::setw(8)
469  << phi << " e " << std::setw(8) << e;
470  }
471  }
472  LogDebug("HcalSim") << "HcalTestAnalysis::EE hits : " << nef << "\n";
473 }
474 
475 //-----------------------------------------------------------------------------
476 void HcalTestAnalysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
477 
478  //Fill tuple with hit information
479  int hittot = caloHitCache.size();
480  tuples->fillHits(caloHitCache);
481 
482  //Get the index of the central tower
483  HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta0,phi0,1,1);
484  uint32_t unitID = org->getUnitID(id);
485  int subdet, zside, layer, ieta, iphi, lay;
486  org->unpackHcalIndex(unitID,subdet,zside,layer,ieta,iphi,lay);
487  int laymax = 0;
488  std::string det = "Unknown";
489  if (subdet == static_cast<int>(HcalBarrel)) {
490  laymax = 4; det = "HB";
491  } else if (subdet == static_cast<int>(HcalEndcap)) {
492  laymax = 2; det = "HES";
493  }
494  LogDebug("HcalSim") << "HcalTestAnalysis::Qie: " << det << " Eta " << ieta
495  << " Phi " << iphi << " Laymax " << laymax << " Hits "
496  << hittot;
497 
498  if (laymax>0 && hittot>0) {
499  std::vector<CaloHit> hits(hittot);
500  std::vector<double> eqielay(80,0.0), esimlay(80,0.0), esimtot(4,0.0);
501  std::vector<double> eqietow(200,0.0), esimtow(200,0.0), eqietot(4,0.0);
502  int etac = (centralTower/100)%100;
503  int phic = (centralTower%100);
504 
505  for (int layr=0; layr<nGroup; layr++) {
506  /*
507  int layx, layy=20;
508  for (int i=0; i<20; i++)
509  if (group_[i] == layr+1 && i < layy) layy = i+1;
510  if (subdet == static_cast<int>(HcalBarrel)) {
511  if (layy < 2) layx = 0;
512  else if (layy < 17) layx = 1;
513  else if (layy == 17) layx = 2;
514  else layx = 3;
515  } else {
516  if (layy < 2) layx = 0;
517  else layx = 1;
518  }
519  */
520  for (int it=0; it<nTower; it++) {
521  int nhit = 0;
522  double esim = 0;
523  for (int k1 = 0; k1 < hittot; k1++) {
524  CaloHit hit = caloHitCache[k1];
525  int subdetc = hit.det();
526  int layer = hit.layer();
527  int group = 0;
528  if (layer > 0 && layer < 20) group = group_[layer];
529  if (subdetc == subdet && group == layr+1) {
530  int zsidec, ietac, iphic, idx;
531  unitID = hit.id();
532  org->unpackHcalIndex(unitID,subdetc,zsidec,layer,ietac,iphic,lay);
533  if (etac > 0 && phic > 0) {
534  idx = ietac*100 + iphic;
535  } else if (etac > 0) {
536  idx = ietac*100;
537  } else if (phic > 0) {
538  idx = iphic;
539  } else {
540  idx = 0;
541  }
542  if (zsidec==zside && idx==tower_[it]) {
543  hits[nhit] = hit;
544  LogDebug("HcalSim") << "HcalTest: Hit " << nhit << " " << hit;
545  nhit++;
546  esim += hit.e();
547  }
548  }
549  }
550 
551  std::vector<int> cd = myqie->getCode(nhit, hits, engine);
552  double eqie = myqie->getEnergy(cd);
553 
554  LogDebug("HcalSim") << "HcalTestAnalysis::Qie: Energy in layer "
555  << layr << " Sim " << esim << " After QIE " <<eqie;
556  for (int i=0; i<4; i++) {
557  if (tower_[nTower+it] <= i) {
558  esimtot[i] += esim;
559  eqietot[i] += eqie;
560  esimlay[20*i+layr] += esim;
561  eqielay[20*i+layr] += eqie;
562  esimtow[50*i+it] += esim;
563  eqietow[50*i+it] += eqie;
564  }
565  }
566  }
567  }
568  LogDebug("HcalSim") << "HcalTestAnalysis::Qie: Total energy " << esimtot[3]
569  << " (SimHit) " << eqietot[3] << " (After QIE)";
570 
571  std::vector<double> latphi(10);
572  int nt = 2*addTower + 1;
573  for (int it=0; it<nt; it++)
574  latphi[it] = it-addTower;
575  for (int i=0; i<4; i++) {
576  double scals=1, scalq=1;
577  std::vector<double> latfs(10,0.), latfq(10,0.), longs(20), longq(20);
578  if (esimtot[i]>0) scals = 1./esimtot[i];
579  if (eqietot[i]>0) scalq = 1./eqietot[i];
580  for (int it=0; it<nTower; it++) {
581  int phib = it%nt;
582  latfs[phib] += scals*esimtow[50*i+it];
583  latfq[phib] += scalq*eqietow[50*i+it];
584  }
585  for (int layr=0; layr<=nGroup; layr++) {
586  longs[layr] = scals*esimlay[20*i+layr];
587  longq[layr] = scalq*eqielay[20*i+layr];
588  }
589  tuples->fillQie(i,esimtot[i],eqietot[i],nGroup,longs,longq,
590  nt,latphi,latfs,latfq);
591  }
592  }
593 }
594 
595 //---------------------------------------------------
596 void HcalTestAnalysis::layerAnalysis(){
597 
598  int i = 0;
599  LogDebug("HcalSim") << "\n ===>>> HcalTestAnalysis: Energy deposit in MeV "
600  << "\n at EB : " << std::setw(6) << edepEB/MeV
601  << "\n at EE : " << std::setw(6) << edepEE/MeV
602  << "\n at HB : " << std::setw(6) << edepHB/MeV
603  << "\n at HE : " << std::setw(6) << edepHE/MeV
604  << "\n at HO : " << std::setw(6) << edepHO/MeV
605  << "\n ---- HcalTestAnalysis: Energy deposit in Layers";
606  for (i = 0; i < 20; i++)
607  LogDebug("HcalSim") << " Layer " << std::setw(2) << i << " E "
608  << std::setw(8) << edepl[i]/MeV << " MeV";
609 
610  tuples->fillLayers(edepl, edepHO, edepHB+edepHE, mudist);
611 }
612 
613 
614 //---------------------------------------------------
615 double HcalTestAnalysis::timeOfFlight(int det, int layer, double eta) {
616 
617  double theta = 2.0*atan(exp(-eta));
618  double dist = 0.;
619  if (det == static_cast<int>(HcalBarrel)) {
620  const double rLay[19] = {
621  1836.0, 1902.0, 1962.0, 2022.0, 2082.0, 2142.0, 2202.0, 2262.0, 2322.0,
622  2382.0, 2448.0, 2514.0, 2580.0, 2646.0, 2712.0, 2776.0, 2862.5, 3847.0,
623  4052.0};
624  if (layer>0 && layer<20) dist += rLay[layer-1]*mm/sin(theta);
625  } else {
626  const double zLay[19] = {
627  4034.0, 4032.0, 4123.0, 4210.0, 4297.0, 4384.0, 4471.0, 4558.0, 4645.0,
628  4732.0, 4819.0, 4906.0, 4993.0, 5080.0, 5167.0, 5254.0, 5341.0, 5428.0,
629  5515.0};
630  if (layer>0 && layer<20) dist += zLay[layer-1]*mm/cos(theta);
631  }
632  double tmp = dist/c_light/ns;
633  LogDebug("HcalSim") << "HcalTestAnalysis::timeOfFlight " << tmp
634  << " for det/lay " << det << " " << layer
635  << " eta/theta " << eta << " " << theta/deg << " dist "
636  << dist;
637  return tmp;
638 }
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
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1087
HcalTestAnalysis::numberingFromDDD
HcalNumberingFromDDD * numberingFromDDD
Definition: HcalTestAnalysis.h:75
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:21
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:89
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:170
HcalTestAnalysis::tuplesManager
std::auto_ptr< HcalTestHistoManager > tuplesManager
Definition: HcalTestAnalysis.h:71
HcalTestHistoClass::fillHits
void fillHits(std::vector< CaloHit > &)
Definition: HcalTestHistoClass.cc:41
fff_deleter.log
tuple log
Definition: fff_deleter.py:13
NULL
#define NULL
Definition: scimark2.h:8
cmsHarvester.loop
int loop
CMSSW
Definition: cmsHarvester.py:4391
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:210
HcalTestAnalysis::names
std::vector< std::string > names
Definition: HcalTestAnalysis.h:84
eta
T eta() const
Definition: Basic3DVectorLD.h:177
HcalTestAnalysis::org
HcalTestNumberingScheme * org
Definition: HcalTestAnalysis.h:76
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:602
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
HcalTestAnalysis::timeOfFlight
double timeOfFlight(int det, int layer, double eta)
Definition: HcalTestAnalysis.cc:615
HCalSD
Definition: HCalSD.h:33
CaloHit::id
uint32_t id() const
Definition: CaloHit.h:28
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
HcalTestAnalysis::fill
void fill(const EndOfEvent *ev)
Definition: HcalTestAnalysis.cc:340
HcalTestNumberingScheme::packHcalIndex
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
Definition: HcalTestNumberingScheme.cc:53
HcalTestAnalysis::tower_
std::vector< int > tower_
Definition: HcalTestAnalysis.h:80
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:48
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:72
nt
int nt
Definition: AMPTWrapper.h:32
HcalTestAnalysis::fileName
std::string fileName
Definition: HcalTestAnalysis.h:64
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
HcalTestAnalysis::caloHitCache
std::vector< CaloHit > caloHitCache
Definition: HcalTestAnalysis.h:79
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
HcalNumberingFromDDD::cell
HcalCellType::HcalCell cell(int det, int zside, int depth, int etaR, int iphi, bool corr=true) const
Definition: HcalNumberingFromDDD.cc:210
HcalTestAnalysis::qieAnalysis
void qieAnalysis(CLHEP::HepRandomEngine *)
Definition: HcalTestAnalysis.cc:476
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:31
HcalTestAnalysis::group_
std::vector< int > group_
Definition: HcalTestAnalysis.h:80
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5