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