CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalTB04Analysis.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HcalTestBeam
4 // Class : HcalTB04Analysis
5 //
6 // Implementation:
7 // Main analysis class for Hcal Test Beam 2004 Analysis
8 //
9 // Original Author:
10 // Created: Tue May 16 10:14:34 CEST 2006
11 //
12 
13 // system include files
14 #include <cmath>
15 #include <iostream>
16 #include <iomanip>
17 
18 // user include files
19 #include "SimG4CMS/HcalTestBeam/interface/HcalTB04Analysis.h"
20 
21 #include "SimG4Core/Notification/interface/BeginOfRun.h"
22 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
23 #include "SimG4Core/Notification/interface/EndOfEvent.h"
24 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
25 
26 // to retreive hits
27 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
28 #include "SimG4CMS/Calo/interface/ECalSD.h"
29 #include "SimG4CMS/Calo/interface/HCalSD.h"
30 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
31 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
32 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
33 #include "SimG4CMS/HcalTestBeam/interface/HcalTBNumberingScheme.h"
34 #include "SimG4CMS/HcalTestBeam/interface/HcalTB04XtalNumberingScheme.h"
35 #include "DataFormats/Math/interface/Point3D.h"
36 
37 #include "FWCore/Framework/interface/Event.h"
38 #include "FWCore/Framework/interface/EventSetup.h"
39 #include "FWCore/Framework/interface/ESHandle.h"
40 #include "FWCore/Framework/interface/MakerMacros.h"
41 #include "FWCore/PluginManager/interface/ModuleDef.h"
42 #include "FWCore/MessageLogger/interface/MessageLogger.h"
43 
44 #include "CLHEP/Random/RandGaussQ.h"
45 
46 #include "G4SDManager.hh"
47 #include "G4VProcess.hh"
48 #include "G4HCofThisEvent.hh"
49 #include "CLHEP/Units/GlobalSystemOfUnits.h"
50 #include "CLHEP/Units/GlobalPhysicalConstants.h"
51 #include "CLHEP/Random/Random.h"
52 
53 //
54 // constructors and destructor
55 //
56 
57 HcalTB04Analysis::HcalTB04Analysis(const edm::ParameterSet &p): myQie(0),
58  histo(0) {
59 
60  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB04Analysis");
61  hcalOnly = m_Anal.getParameter<bool>("HcalOnly");
62  mode = m_Anal.getParameter<int>("Mode");
63  type = m_Anal.getParameter<int>("Type");
64  ecalNoise = m_Anal.getParameter<double>("EcalNoise");
65  scaleHB0 = m_Anal.getParameter<double>("ScaleHB0");
66  scaleHB16 = m_Anal.getParameter<double>("ScaleHB16");
67  scaleHO = m_Anal.getParameter<double>("ScaleHO");
68  scaleHE0 = m_Anal.getParameter<double>("ScaleHE0");
69  names = m_Anal.getParameter<std::vector<std::string> >("Names");
70  beamOffset =-m_Anal.getParameter<double>("BeamPosition")*cm;
71  double fMinEta = m_Anal.getParameter<double>("MinEta");
72  double fMaxEta = m_Anal.getParameter<double>("MaxEta");
73  double fMinPhi = m_Anal.getParameter<double>("MinPhi");
74  double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
75  double beamEta = (fMaxEta+fMinEta)/2.;
76  double beamPhi = (fMaxPhi+fMinPhi)/2.;
77  double beamThet= 2*atan(exp(-beamEta));
78  if (beamPhi < 0) beamPhi += twopi;
79  iceta = (int)(beamEta/0.087) + 1;
80  icphi = (int)(fabs(beamPhi)/0.087) + 5;
81  if (icphi > 72) icphi -= 73;
82 
83  produces<PHcalTB04Info>();
84 
85  beamline_RM = new G4RotationMatrix;
86  beamline_RM->rotateZ(-beamPhi);
87  beamline_RM->rotateY(-beamThet);
88 
89  edm::LogInfo("HcalTBSim") << "HcalTB04:: Initialised as observer of BeginOf"
90  << "Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
91  << " with Parameter values:\n \thcalOnly = "
92  << hcalOnly << "\tecalNoise = " << ecalNoise
93  << "\n\tMode = " << mode << " (0: HB2 Standard; "
94  << "1:HB2 Segmented)" << "\tType = " << type
95  << " (0: HB; 1 HE; 2 HB+HE)\n\tbeamOffset = "
96  << beamOffset << "\ticeta = " << iceta
97  << "\ticphi = " << icphi << "\n\tbeamline_RM = "
98  << *beamline_RM;
99 
100  init();
101 
102  myQie = new HcalQie(p);
103  histo = new HcalTB04Histo(m_Anal);
104 }
105 
106 HcalTB04Analysis::~HcalTB04Analysis() {
107 
108  edm::LogInfo("HcalTBSim") << "\n --------> Total number of selected entries"
109  << " : " << count << "\nPointers:: QIE " << myQie
110  << " Histo " << histo;
111  if (myQie) {
112  delete myQie;
113  myQie = 0;
114  }
115  if (histo) {
116  delete histo;
117  histo = 0;
118  }
119 }
120 
121 //
122 // member functions
123 //
124 
125 void HcalTB04Analysis::produce(edm::Event& e, const edm::EventSetup&) {
126 
127  std::auto_ptr<PHcalTB04Info> product(new PHcalTB04Info);
128  fillEvent(*product);
129  e.put(product);
130 }
131 
132 void HcalTB04Analysis::init() {
133 
134  idTower = HcalTBNumberingScheme::getUnitIDs(type, mode);
135  nTower = idTower.size();
136  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Save information from "
137  << nTower << " HCal towers";
138  idHcal.reserve(nTower);
139  for (int i=0; i<nTower; i++) {
140  int id = unitID(idTower[i]);
141  idHcal.push_back(id);
142  LogDebug("HcalTBSim") << "\tTower[" << i << "] Original " << std::hex
143  << idTower[i] << " Stored " << idHcal[i] << std::dec;
144  }
145 
146  if (!hcalOnly) {
147  int det = 10;
148  uint32_t id1;
149  nCrystal = 0;
150  for (int lay=1; lay<8; lay++) {
151  for (int icr=1; icr<8; icr++) {
152  id1 = HcalTestNumbering::packHcalIndex(det,0,1,icr,lay,1);
153  int id = unitID(id1);
154  idEcal.push_back(id1);
155  idXtal.push_back(id);
156  nCrystal++;
157  }
158  }
159  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Save information from "
160  << nCrystal << " ECal Crystals";
161  for (int i=0; i<nCrystal; i++) {
162  LogDebug("HcalTBSim") << "\tCrystal[" << i << "] Original " << std::hex
163  << idEcal[i] << " Stored " << idXtal[i] <<std::dec;
164  }
165  }
166  // Profile vectors
167  eseta.reserve(5);
168  eqeta.reserve(5);
169  esphi.reserve(3);
170  eqphi.reserve(3);
171  eslay.reserve(20);
172  eqlay.reserve(20);
173  for (int i=0; i<5; i++) {
174  eseta.push_back(0.);
175  eqeta.push_back(0.);
176  }
177  for (int i=0; i<3; i++) {
178  esphi.push_back(0.);
179  eqphi.push_back(0.);
180  }
181  for (int i=0; i<20; i++) {
182  eslay.push_back(0.);
183  eqlay.push_back(0.);
184  }
185 
186  // counter
187  count = 0;
188  evNum = 0;
189  clear();
190 }
191 
192 void HcalTB04Analysis::update(const BeginOfRun * run) {
193 
194  int irun = (*run)()->GetRunID();
195  edm::LogInfo("HcalTBSim") <<" =====> Begin of Run = " << irun;
196 
197  G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
198  if (sd != 0) {
199  std::string sdname = names[0];
200  G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
201  if (aSD==0) {
202  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
203  << " with name " << sdname << " in this "
204  << "Setup";
205  } else {
206  HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
207  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD "
208  << "with name " << theCaloSD->GetName()
209  << " in this Setup";
210  HcalNumberingScheme* org = new HcalTestNumberingScheme(false);
211  theCaloSD->setNumberingScheme(org);
212  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a "
213  << "new numbering scheme";
214  }
215  if (!hcalOnly) {
216  sdname = names[1];
217  aSD = sd->FindSensitiveDetector(sdname);
218  if (aSD==0) {
219  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
220  << " with name " << sdname << " in this "
221  << "Setup";
222  } else {
223  ECalSD* theCaloSD = dynamic_cast<ECalSD*>(aSD);
224  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD "
225  << "with name " << theCaloSD->GetName()
226  << " in this Setup";
227  EcalNumberingScheme* org = new HcalTB04XtalNumberingScheme();
228  theCaloSD->setNumberingScheme(org);
229  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a "
230  << "new numbering scheme";
231  }
232  }
233  } else {
234  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Could "
235  << "not get SD Manager!";
236  }
237 
238 }
239 
240 void HcalTB04Analysis::update(const BeginOfEvent * evt) {
241 
242  evNum = (*evt) ()->GetEventID ();
243  clear();
244  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis: =====> Begin of event = "
245  << evNum;
246 }
247 
248 void HcalTB04Analysis::update(const G4Step * aStep) {
249 
250  if (aStep != NULL) {
251  //Get Step properties
252  G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
253  G4ThreeVector thePostStepPoint;
254 
255  // Get Tracks properties
256  G4Track* aTrack = aStep->GetTrack();
257  int trackID = aTrack->GetTrackID();
258  int parentID = aTrack->GetParentID();
259  G4ThreeVector position = aTrack->GetPosition();
260  G4ThreeVector momentum = aTrack->GetMomentum();
261  G4String partType = aTrack->GetDefinition()->GetParticleType();
262  G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
263  int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
264 #ifdef ddebug
265  bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
266 #endif
267  double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
268  double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
269 
270  if (!pvFound) { //search for v1
271  double stepDeltaEnergy = aStep->GetDeltaEnergy ();
272  double kinEnergy = aTrack->GetKineticEnergy ();
273 
274  // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
275  if (trackID == 1 && parentID == 0 &&
276  ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
277  pvType = -1;
278  if (kinEnergy == 0.) {
279  pvType = 0;
280  } else {
281  if (fabs (stepDeltaEnergy / kinEnergy) > 0.1) pvType = 1;
282  }
283  pvFound = true;
284  pvPosition = position;
285  pvMomentum = momentum;
286  // Rotated coord.system:
287  pvUVW = (*beamline_RM)*(pvPosition);
288 
289  //Volume name requires some checks:
290  G4String thePostPVname = "NoName";
291  G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
292  if (thePostPoint) {
293  thePostStepPoint = thePostPoint->GetPosition();
294  G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
295  if (thePostPV) thePostPVname = thePostPV->GetName ();
296  }
297 #ifdef ddebug
298  LogDebug("HcalTBSim") << "HcalTB04Analysis:: V1 found at: "
299  << thePostStepPoint << " G4VPhysicalVolume: "
300  << thePostPVname;
301 #endif
302  LogDebug("HcalTBSim") << "HcalTB04Analysis::fill_v1Pos: Primary Track "
303  << "momentum: " << pvMomentum << " psoition "
304  << pvPosition << " u/v/w " << pvUVW;
305  }
306  } else {
307  // watch for secondaries originating @v1, including the surviving primary
308  if ((trackID != 1 && parentID == 1 &&
309  (aTrack->GetCurrentStepNumber () == 1) &&
310  (thePreStepPoint == pvPosition)) ||
311  (trackID == 1 && thePreStepPoint == pvPosition)) {
312 #ifdef ddebug
313  LogDebug("HcalTBSim") << "HcalTB04Analysis::A secondary... PDG:"
314  << partPDGEncoding << " TrackID:" << trackID
315  << " ParentID:" << parentID << " stable: "
316  << isPDGStable << " Tau: " << pDGlifetime
317  << " cTauGamma="
318  << c_light*pDGlifetime*gammaFactor*1000.
319  << "um" << " GammaFactor: " << gammaFactor;
320 #endif
321  secTrackID.push_back(trackID);
322  secPartID.push_back(partPDGEncoding);
323  secMomentum.push_back(momentum);
324  secEkin.push_back(aTrack->GetKineticEnergy());
325 
326  // Check for short-lived secondaries: cTauGamma<100um
327  double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
328  if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {//short-lived secondary
329  shortLivedSecondaries.push_back(trackID);
330  } else {//normal secondary - enter into the V1-calorimetric tree
331  // histos->fill_v1cSec (aTrack);
332  }
333  }
334  // Also watch for tertiary particles coming from
335  // short-lived secondaries from V1
336  if (aTrack->GetCurrentStepNumber() == 1) {
337  if (shortLivedSecondaries.size() > 0) {
338  int pid = parentID;
339  std::vector<int>::iterator pos1= shortLivedSecondaries.begin();
340  std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
341  std::vector<int>::iterator pos;
342  for (pos = pos1; pos != pos2; pos++) {
343  if (*pos == pid) {//ParentID is on the list of short-lived
344  // secondary
345 #ifdef ddebug
346  LogDebug("HcalTBSim") << "HcalTB04Analysis:: A tertiary... PDG:"
347  << partPDGEncoding << " TrackID:" <<trackID
348  << " ParentID:" << parentID << " stable: "
349  << isPDGStable << " Tau: " << pDGlifetime
350  << " cTauGamma="
351  << c_light*pDGlifetime*gammaFactor*1000.
352  << "um GammaFactor: " << gammaFactor;
353 #endif
354  }
355  }
356  }
357  }
358  }
359  }
360 }
361 
362 void HcalTB04Analysis::update(const EndOfEvent * evt) {
363 
364  count++;
365 
366  //fill the buffer
367  LogDebug("HcalTBSim") << "HcalTB04Analysis::Fill event "
368  << (*evt)()->GetEventID();
369  fillBuffer (evt);
370 
371  //QIE analysis
372  LogDebug("HcalTBSim") << "HcalTB04Analysis::Do QIE analysis with "
373  << hcalHitCache.size() << " hits";
374  CLHEP::HepRandomEngine* engine = CLHEP::HepRandom::getTheEngine();
375  qieAnalysis(engine);
376 
377  //Energy in Crystal Matrix
378  if (!hcalOnly) {
379  LogDebug("HcalTBSim") << "HcalTB04Analysis::Do Xtal analysis with "
380  << ecalHitCache.size() << " hits";
381  xtalAnalysis(engine);
382  }
383 
384  //Final Analysis
385  LogDebug("HcalTBSim") << "HcalTB04Analysis::Final analysis";
386  finalAnalysis();
387 
388  int iEvt = (*evt)()->GetEventID();
389  if (iEvt < 10)
390  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
391  else if ((iEvt < 100) && (iEvt%10 == 0))
392  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
393  else if ((iEvt < 1000) && (iEvt%100 == 0))
394  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
395  else if ((iEvt < 10000) && (iEvt%1000 == 0))
396  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
397 }
398 
399 void HcalTB04Analysis::fillBuffer(const EndOfEvent * evt) {
400 
401  std::vector<CaloHit> hhits, hhitl;
402  int idHC, j;
403  CaloG4HitCollection* theHC;
404  std::map<int,float,std::less<int> > primaries;
405  double etot1=0, etot2=0;
406 
407  // Look for the Hit Collection of HCal
408  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
409  std::string sdName = names[0];
410  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
411  theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
412  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName
413  << " of ID " << idHC << " is obtained at " << theHC;
414 
415  if (idHC >= 0 && theHC > 0) {
416  hhits.reserve(theHC->entries());
417  hhitl.reserve(theHC->entries());
418  for (j = 0; j < theHC->entries(); j++) {
419  CaloG4Hit* aHit = (*theHC)[j];
420  double e = aHit->getEnergyDeposit()/GeV;
421  double time = aHit->getTimeSlice();
422  math::XYZPoint pos = aHit->getEntry();
423  unsigned int id = aHit->getUnitID();
424  double theta = pos.theta();
425  double eta = -log(tan(theta * 0.5));
426  double phi = pos.phi();
427  int det, z, group, ieta, iphi, layer;
428  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
429  double jitter = time-timeOfFlight(det,layer,eta);
430  if (jitter<0) jitter = 0;
431  if (e < 0 || e > 1.) e = 0;
432  double escl = e * scale(det,layer);
433  unsigned int idx= HcalTBNumberingScheme::getUnitID(id,mode);
434  CaloHit hit(det,layer,escl,eta,phi,jitter,idx);
435  hhits.push_back(hit);
436  CaloHit hitl(det,layer,escl,eta,phi,jitter,id);
437  hhitl.push_back(hitl);
438  primaries[aHit->getTrackID()]+= e;
439  etot1 += escl;
440 #ifdef ddebug
441  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j
442  << " ID 0x" << std::hex << id << " 0x" << idx
443  << std::dec << " time " << std::setw(6) << time
444  << " " << std::setw(6) << jitter << " theta "
445  << std::setw(8) << theta << " eta " << std::setw(8)
446  << eta << " phi " << std::setw(8) << phi << " e "
447  << std::setw(8) << e << " " << std::setw(8) <<escl;
448 #endif
449  }
450  }
451 
452  // Add hits in the same channel within same time slice
453  std::vector<CaloHit>::iterator itr;
454  int nHit = hhits.size();
455  std::vector<CaloHit*> hits(nHit);
456  for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
457  hits[j] = &hhits[j];
458  }
459  sort(hits.begin(),hits.end(),CaloHitIdMore());
460  std::vector<CaloHit*>::iterator k1, k2;
461  int nhit = 0;
462  for (k1 = hits.begin(); k1 != hits.end(); k1++) {
463  int det = (**k1).det();
464  int layer = (**k1).layer();
465  double ehit = (**k1).e();
466  double eta = (**k1).eta();
467  double phi = (**k1).phi();
468  double jitter = (**k1).t();
469  uint32_t unitID = (**k1).id();
470  int jump = 0;
471  for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
472  unitID==(**k2).id(); k2++) {
473  ehit += (**k2).e();
474  jump++;
475  }
476  nhit++;
477  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
478  hcalHitCache.push_back(hit);
479  etot2 += ehit;
480  k1 += jump;
481 #ifdef ddebug
482  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit
483  << " ID 0x" << std::hex << unitID << std::dec
484  << " time " << std::setw(6) << jitter << " eta "
485  << std::setw(8) << eta << " phi " << std::setw(8)
486  << phi << " e " << std::setw(8) << ehit;
487 #endif
488  }
489  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits"
490  << " from " << nHit << " input hits E(Hcal) " << etot1
491  << " " << etot2;
492 
493  //Repeat for Hit in each layer (hhits and hhitl sizes are the same)
494  for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
495  hits[j] = &hhitl[j];
496  }
497  sort(hits.begin(),hits.end(),CaloHitIdMore());
498  int nhitl = 0;
499  double etotl = 0;
500  for (k1 = hits.begin(); k1 != hits.end(); k1++) {
501  int det = (**k1).det();
502  int layer = (**k1).layer();
503  double ehit = (**k1).e();
504  double eta = (**k1).eta();
505  double phi = (**k1).phi();
506  double jitter = (**k1).t();
507  uint32_t unitID = (**k1).id();
508  int jump = 0;
509  for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
510  unitID==(**k2).id(); k2++) {
511  ehit += (**k2).e();
512  jump++;
513  }
514  nhitl++;
515  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
516  hcalHitLayer.push_back(hit);
517  etotl += ehit;
518  k1 += jump;
519 #ifdef ddebug
520  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl
521  << " ID 0x" << std::hex << unitID << std::dec
522  << " time " << std::setw(6) << jitter << " eta "
523  << std::setw(8) << eta << " phi " << std::setw(8)
524  << phi << " e " << std::setw(8) << ehit;
525 #endif
526  }
527  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal "
528  << "hits from " << nHit << " input hits E(Hcal) "
529  << etot1 << " " << etotl;
530 
531  // Look for the Hit Collection of ECal
532  std::vector<CaloHit> ehits;
533  sdName= names[1];
534  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
535  theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
536  etot1 = etot2 = 0;
537  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName
538  << " of ID " << idHC << " is obtained at " << theHC;
539  if (idHC >= 0 && theHC > 0) {
540  ehits.reserve(theHC->entries());
541  for (j = 0; j < theHC->entries(); j++) {
542  CaloG4Hit* aHit = (*theHC)[j];
543  double e = aHit->getEnergyDeposit()/GeV;
544  double time = aHit->getTimeSlice();
545  math::XYZPoint pos = aHit->getEntry();
546  unsigned int id = aHit->getUnitID();
547  double theta = pos.theta();
548  double eta = -log(tan(theta * 0.5));
549  double phi = pos.phi();
550  if (e < 0 || e > 100000.) e = 0;
551  int det, z, group, ieta, iphi, layer;
552  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
553  CaloHit hit(det,0,e,eta,phi,time,id);
554  ehits.push_back(hit);
555  primaries[aHit->getTrackID()]+= e;
556  etot1 += e;
557 #ifdef ddebug
558  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j
559  << " ID 0x" << std::hex << id << std::dec
560  << " time " << std::setw(6) << time << " theta "
561  << std::setw(8) << theta << " eta " <<std::setw(8)
562  << eta << " phi " << std::setw(8) << phi << " e "
563  << std::setw(8) << e;
564 #endif
565  }
566  }
567 
568  // Add hits in the same channel within same time slice
569  nHit = ehits.size();
570  std::vector<CaloHit*> hite(nHit);
571  for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
572  hite[j] = &ehits[j];
573  }
574  sort(hite.begin(),hite.end(),CaloHitIdMore());
575  nhit = 0;
576  for (k1 = hite.begin(); k1 != hite.end(); k1++) {
577  int det = (**k1).det();
578  int layer = (**k1).layer();
579  double ehit = (**k1).e();
580  double eta = (**k1).eta();
581  double phi = (**k1).phi();
582  double jitter = (**k1).t();
583  uint32_t unitID = (**k1).id();
584  int jump = 0;
585  for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
586  unitID==(**k2).id(); k2++) {
587  ehit += (**k2).e();
588  jump++;
589  }
590  nhit++;
591  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
592  ecalHitCache.push_back(hit);
593  etot2 += ehit;
594  k1 += jump;
595 #ifdef ddebug
596  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit
597  << " ID 0x" << std::hex << unitID << std::dec
598  << " time " << std::setw(6) << jitter << " eta "
599  << std::setw(8) << eta << " phi " << std::setw(8)
600  << phi << " e " << std::setw(8) << ehit;
601 #endif
602  }
603  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits"
604  << " from " << nHit << " input hits E(Ecal) " << etot1
605  << " " << etot2;
606 
607  // Find Primary info:
608  nPrimary = (int)(primaries.size());
609  int trackID = 0;
610  G4PrimaryParticle* thePrim=0;
611  int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
612  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex
613  << " verteices";
614  if (nvertex<=0)
615  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no "
616  << "vertex found for event " << evNum;
617 
618  for (int i = 0 ; i<nvertex; i++) {
619  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
620  if (avertex == 0) {
621  edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer "
622  << "to vertex = 0 for event " << evNum;
623  } else {
624  LogDebug("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " "
625  << avertex->GetPosition();
626  int npart = avertex->GetNumberOfParticle();
627  if (npart == 0)
628  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: "
629  << "no primary!";
630  if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
631  }
632  }
633 
634  if (thePrim != 0) {
635  double px = thePrim->GetPx();
636  double py = thePrim->GetPy();
637  double pz = thePrim->GetPz();
638  double p = std::sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
639  pInit = p/GeV;
640  if (p==0)
641  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: "
642  << "primary has p=0 ";
643  else {
644  double costheta = pz/p;
645  double theta = acos(std::min(std::max(costheta,-1.),1.));
646  etaInit = -log(tan(theta/2));
647  if (px != 0 || py != 0) phiInit = atan2(py,px);
648  }
649  particleType = thePrim->GetPDGcode();
650  } else
651  edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could "
652  << "not find primary";
653 
654 }
655 
656 void HcalTB04Analysis::qieAnalysis(CLHEP::HepRandomEngine* engine) {
657 
658  int hittot = hcalHitCache.size();
659  if (hittot<=0) hittot = 1;
660  std::vector<CaloHit> hits(hittot);
661  std::vector<int> todo(nTower,0);
662 
663  LogDebug("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size "
664  << hits.size() << " " << todo.size() << " "
665  << idTower.size() << " " << esimh.size() << " "
666  << eqie.size();
667  // Loop over all HCal hits
668  for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
669  CaloHit hit = hcalHitCache[k1];
670  uint32_t id = hit.id();
671  int nhit = 0;
672  double esim = hit.e();
673  hits[nhit] = hit;
674  for (unsigned int k2 = k1+1; k2 < hcalHitCache.size(); k2++) {
675  hit = hcalHitCache[k2];
676  if (hit.id() == id) {
677  nhit++;
678  hits[nhit] = hit;
679  esim += hit.e();
680  }
681  }
682  k1 += nhit;
683  nhit++;
684  std::vector<int> cd = myQie->getCode(nhit, hits, engine);
685  double eq = myQie->getEnergy(cd);
686  LogDebug("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id
687  << std::dec << " registers " << esim << " energy "
688  << "from " << nhit << " hits starting with hit # "
689  << k1 << " energy with noise " << eq;
690  for (int k2 = 0; k2 < nTower; k2++) {
691  if (id == idTower[k2]) {
692  todo[k2] = 1;
693  esimh[k2] = esim;
694  eqie[k2] = eq;
695  }
696  }
697  }
698 
699  // Towers with no hit
700  for (int k2 = 0; k2 < nTower; k2++) {
701  if (todo[k2] == 0) {
702  std::vector<int> cd = myQie->getCode(0, hits, engine);
703  double eq = myQie->getEnergy(cd);
704  esimh[k2] = 0;
705  eqie[k2] = eq;
706 #ifdef ddebug
707  LogDebug("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex
708  << idTower[k2] << std::dec << " registers "
709  << esimh[k2] << " energy from hits and energy "
710  << "after QIE analysis " << eqie[k2];
711 #endif
712  }
713  }
714 }
715 
716 void HcalTB04Analysis::xtalAnalysis(CLHEP::HepRandomEngine* engine) {
717 
718  CLHEP::RandGaussQ randGauss(*engine);
719 
720  // Crystal Data
721  std::vector<int> iok(nCrystal,0);
722  LogDebug("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " <<iok.size()
723  << " " << idEcal.size() << " " << esime.size() << " "
724  << enois.size();
725  for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
726  uint32_t id = ecalHitCache[k1].id();
727  int nhit = 0;
728  double esim = ecalHitCache[k1].e();
729  for (unsigned int k2 = k1+1; k2 < ecalHitCache.size(); k2++) {
730  if (ecalHitCache[k2].id() == id) {
731  nhit++;
732  esim += ecalHitCache[k2].e();
733  }
734  }
735  k1 += nhit;
736  nhit++;
737  double eq = esim + randGauss.fire(0., ecalNoise);
738 #ifdef ddebug
739  LogDebug("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex << id
740  << std::dec << " registers " << esim << " energy "
741  << "from " << nhit << " hits starting with hit # "
742  << k1 << " energy with noise " << eq;
743 #endif
744  for (int k2 = 0; k2 < nCrystal; k2++) {
745  if (id == idEcal[k2]) {
746  iok[k2] = 1;
747  esime[k2] = esim;
748  enois[k2] = eq;
749  }
750  }
751  }
752 
753  // Crystals with no hit
754  for (int k2 = 0; k2 < nCrystal; k2++) {
755  if (iok[k2] == 0) {
756  esime[k2] = 0;
757  enois[k2] = randGauss.fire(0., ecalNoise);
758 #ifdef ddebug
759  LogDebug("HcalTBSim") << "HcalTB04Analysis:: ID 0x" << std::hex
760  << idEcal[k2] << std::dec << " registers "
761  << esime[k2] << " energy from hits and energy from"
762  << " noise " << enois[k2];
763 #endif
764  }
765  }
766 }
767 
768 void HcalTB04Analysis::finalAnalysis() {
769 
770  //Beam Information
771  histo->fillPrimary(pInit, etaInit, phiInit);
772 
773  // Total Energy
774  eecals = ehcals = eecalq = ehcalq = 0.;
775  for (int i=0; i<nTower; i++) {
776  ehcals += esimh[i];
777  ehcalq += eqie[i];
778  }
779  for (int i=0; i<nCrystal; i++) {
780  eecals += esime[i];
781  eecalq += enois[i];
782  }
783  etots = eecals + ehcals;
784  etotq = eecalq + ehcalq;
785  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level "
786  << "(Total) " << etots << " (ECal) " << eecals
787  << " (HCal) " << ehcals << "\nHcalTB04Analysis:: "
788  << "Energy deposit at Qie Level (Total) " << etotq
789  << " (ECal) " << eecalq << " (HCal) " << ehcalq;
790  histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
791 
792  // Lateral Profile
793  for (int i=0; i<5; i++) {
794  eseta[i] = 0.;
795  eqeta[i] = 0.;
796  }
797  for (int i=0; i<3; i++) {
798  esphi[i] = 0.;
799  eqphi[i] = 0.;
800  }
801  double e1=0, e2=0;
802  unsigned int id;
803  for (int i=0; i<nTower; i++) {
804  int det, z, group, ieta, iphi, layer;
805  id = idTower[i];
806  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
807  iphi -= (icphi - 1);
808  if (icphi > 4) {
809  if (ieta == 0) ieta = 2;
810  else ieta =-1;
811  } else {
812  ieta = ieta - iceta + 2;
813  }
814  if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
815  eseta[ieta] += esimh[i];
816  esphi[iphi] += esimh[i];
817  e1 += esimh[i];
818  eqeta[ieta] += eqie[i];
819  eqphi[iphi] += eqie[i];
820  e2 += eqie[i];
821  }
822  }
823  for (int i=0; i<3; i++) {
824  if (e1>0) esphi[i] /= e1;
825  if (e2>0) eqphi[i] /= e2;
826  }
827  for (int i=0; i<5; i++) {
828  if (e1>0) eseta[i] /= e1;
829  if (e2>0) eqeta[i] /= e2;
830  }
831  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and"
832  << " Phi (Sim/Qie)";
833  for (int i=0; i<5; i++)
834  LogDebug("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = "
835  << eseta[i] << " Qie = " << eqeta[i] << " Phi Sim = "
836  << esphi[i] << " Qie = " << eqphi[i];
837  histo->fillTrnsProf(eseta,eqeta,esphi,eqphi);
838 
839  // Longitudianl profile
840  for (int i=0; i<20; i++) {
841  eslay[i] = 0.;
842  eqlay[i] = 0.;
843  }
844  e1=0; e2=0;
845  for (int i=0; i<nTower; i++) {
846  int det, z, group, ieta, iphi, layer;
847  id = idTower[i];
848  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
849  iphi -= (icphi - 1);
850  layer -= 1;
851  if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
852  eslay[layer] += esimh[i];
853  e1 += esimh[i];
854  eqlay[layer] += eqie[i];
855  e2 += eqie[i];
856  }
857  }
858  for (int i=0; i<20; i++) {
859  if (e1>0) eslay[i] /= e1;
860  if (e2>0) eqlay[i] /= e2;
861  }
862  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
863  for (int i=0; i<20; i++)
864  LogDebug("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = "
865  << eslay[i] << " Qie = " << eqlay[i];
866  histo->fillLongProf(eslay, eqlay);
867 }
868 
869 
870 void HcalTB04Analysis::fillEvent (PHcalTB04Info& product) {
871 
872  //Setup the ID's
873  product.setIDs(idHcal, idXtal);
874 
875  //Beam Information
876  product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
877 
878  //Energy deposits in the crystals and towers
879  product.setEdepHcal(esimh, eqie);
880  product.setEdepHcal(esime, enois);
881 
882  // Total Energy
883  product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
884 
885  // Lateral Profile
886  product.setTrnsProf(eseta,eqeta,esphi,eqphi);
887 
888  // Longitudianl profile
889  product.setLongProf(eslay, eqlay);
890 
891  //Save Hits
892  int i, nhit=0;
893  std::vector<CaloHit>::iterator itr;
894  for (i=0, itr=ecalHitCache.begin(); itr!=ecalHitCache.end(); i++,itr++) {
895  uint32_t id = itr->id();
896  int det, z, group, ieta, iphi, lay;
897  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,lay);
898  product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
899  nhit++;
900 #ifdef debug
901  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3)
902  << i+1 << " ID 0x" << std::hex << group << std::dec
903  << " " << std::setw(2) << det << " " << std::setw(2)
904  << lay << " " << std::setw(1) << z << " "
905  << std::setw(3) << ieta << " " << std::setw(3) <<iphi
906  << " T " << std::setw(6) << itr->t() << " E "
907  << std::setw(6) << itr->e();
908 #endif
909  }
910  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit
911  << " hits from Crystals";
912  int hit = nhit;
913  nhit = 0;
914 
915  for (i=hit, itr=hcalHitCache.begin(); itr!=hcalHitCache.end(); i++,itr++) {
916  uint32_t id = itr->id();
917  int det, z, group, ieta, iphi, lay;
918  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,lay);
919  product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
920  nhit++;
921 #ifdef debug
922  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3)
923  << i+1 << " ID 0x" << std::hex << group << std::dec
924  << " " << std::setw(2) << det << " " << std::setw(2)
925  << lay << " " << std::setw(1) << z << " "
926  << std::setw(3) << ieta << " " << std::setw(3) <<iphi
927  << " T " << std::setw(6) << itr->t() << " E "
928  << std::setw(6) << itr->e();
929 #endif
930  }
931  LogDebug("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit
932  << " hits from HCal";
933 
934  //Vertex associated quantities
935  product.setVtxPrim(evNum, pvType, pvPosition.x(), pvPosition.y(),
936  pvPosition.z(), pvUVW.x(), pvUVW.y(), pvUVW.z(),
937  pvMomentum.x(), pvMomentum.y(), pvMomentum.z());
938  for (unsigned int i = 0; i < secTrackID.size(); i++) {
939  product.setVtxSec(secTrackID[i], secPartID[i], secMomentum[i].x(),
940  secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
941  }
942 }
943 
944 void HcalTB04Analysis::clear(){
945  pvFound = false;
946  pvType =-2;
947  pvPosition = G4ThreeVector();
948  pvMomentum = G4ThreeVector();
949  pvUVW = G4ThreeVector();
950  secTrackID.clear();
951  secPartID.clear();
952  secMomentum.clear();
953  secEkin.clear();
954  shortLivedSecondaries.clear();
955 
956  ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
957  hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
958  hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
959  nPrimary = particleType = 0;
960  pInit = etaInit = phiInit = 0;
961 
962  esimh.clear();
963  eqie.clear();
964  esimh.reserve(nTower);
965  eqie.reserve(nTower);
966  for (int i=0; i<nTower; i++) {
967  esimh.push_back(0.);
968  eqie.push_back(0.);
969  }
970  esime.clear();
971  enois.clear();
972  esime.reserve(nCrystal);
973  enois.reserve(nCrystal);
974  for (int i=0; i<nCrystal; i++) {
975  esime.push_back(0.);
976  enois.push_back(0.);
977  }
978 }
979 
980 int HcalTB04Analysis::unitID(uint32_t id) {
981 
982  int det, z, group, ieta, iphi, lay;
983  HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,lay);
984  group = (det&15)<<20;
985  group += ((lay-1)&31)<<15;
986  group += (z&1)<<14;
987  group += (ieta&127)<<7;
988  group += (iphi&127);
989  return group;
990 }
991 
992 double HcalTB04Analysis::scale(int det, int layer) {
993 
994  double tmp = 1.;
995  if (det == static_cast<int>(HcalBarrel)) {
996  if (layer == 1) tmp = scaleHB0;
997  else if (layer == 17) tmp = scaleHB16;
998  else if (layer > 17) tmp = scaleHO;
999  } else {
1000  if (layer <= 2) tmp = scaleHE0;
1001  }
1002  return tmp;
1003 }
1004 
1005 double HcalTB04Analysis::timeOfFlight(int det, int layer, double eta) {
1006 
1007  double theta = 2.0*atan(exp(-eta));
1008  double dist = beamOffset;
1009  if (det == static_cast<int>(HcalBarrel)) {
1010  const double rLay[19] = {
1011  1836.0, 1902.0, 1962.0, 2022.0, 2082.0, 2142.0, 2202.0, 2262.0, 2322.0,
1012  2382.0, 2448.0, 2514.0, 2580.0, 2646.0, 2712.0, 2776.0, 2862.5, 3847.0,
1013  4052.0};
1014  if (layer>0 && layer<=19) dist += rLay[layer-1]*mm/sin(theta);
1015  } else {
1016  const double zLay[19] = {
1017  4034.0, 4032.0, 4123.0, 4210.0, 4297.0, 4384.0, 4471.0, 4558.0, 4645.0,
1018  4732.0, 4819.0, 4906.0, 4993.0, 5080.0, 5167.0, 5254.0, 5341.0, 5428.0,
1019  5515.0};
1020  if (layer>0 && layer<=19) dist += zLay[layer-1]*mm/cos(theta);
1021  }
1022 
1023  double tmp = dist/c_light/ns;
1024 #ifdef ddebug
1025  LogDebug("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp
1026  << " for det/lay " << det << " " << layer
1027  << " eta/theta " << eta << " " << theta/deg
1028  << " dist " << dist;
1029 #endif
1030  return tmp;
1031 }
1032 
1033 DEFINE_SIMWATCHER (HcalTB04Analysis);
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HcalTB04Analysis::update
void update(const BeginOfRun *run)
This routine will be called when the appropriate signal arrives.
Definition: HcalTB04Analysis.cc:192
i
int i
Definition: DBlmapReader.cc:9
DEFINE_SIMWATCHER
#define DEFINE_SIMWATCHER(type)
Definition: SimWatcherFactory.h:13
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1087
HcalTB04Analysis::secEkin
std::vector< double > secEkin
Definition: HcalTB04Analysis.h:123
PHcalTB04Info::setLongProf
void setLongProf(const std::vector< double > &es, const std::vector< double > &eq)
Definition: PHcalTB04Info.cc:175
HcalTB04Analysis::idHcal
std::vector< int > idHcal
Definition: HcalTB04Analysis.h:106
HcalTB04Histo::fillPrimary
void fillPrimary(double energy, double eta, double phi)
Definition: HcalTB04Histo.cc:65
HcalTB04Analysis::hcalHitLayer
std::vector< CaloHit > hcalHitLayer
Definition: HcalTB04Analysis.h:113
HcalTB04Analysis::eqeta
std::vector< double > eqeta
Definition: HcalTB04Analysis.h:115
HcalTB04Analysis::eqie
std::vector< double > eqie
Definition: HcalTB04Analysis.h:114
HcalTB04Histo::fillTrnsProf
void fillTrnsProf(const std::vector< double > &es1, const std::vector< double > &eq1, const std::vector< double > &es2, const std::vector< double > &eq2)
Definition: HcalTB04Histo.cc:91
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HcalTB04Analysis::produce
virtual void produce(edm::Event &, const edm::EventSetup &)
Definition: HcalTB04Analysis.cc:125
HcalTB04Analysis::fillBuffer
void fillBuffer(const EndOfEvent *evt)
Definition: HcalTB04Analysis.cc:399
HcalTB04Analysis::shortLivedSecondaries
std::vector< int > shortLivedSecondaries
Definition: HcalTB04Analysis.h:124
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:170
HcalTB04Analysis::secTrackID
std::vector< int > secTrackID
Definition: HcalTB04Analysis.h:121
fff_deleter.log
tuple log
Definition: fff_deleter.py:13
NULL
#define NULL
Definition: scimark2.h:8
npart
double npart
Definition: HydjetWrapper.h:44
HcalTB04Analysis::eqphi
std::vector< double > eqphi
Definition: HcalTB04Analysis.h:115
PHcalTB04Info::setPrimary
void setPrimary(int primary, int id, double energy, double eta, double phi)
Definition: PHcalTB04Info.cc:102
eta
T eta() const
Definition: Basic3DVectorLD.h:177
HCalSD::setNumberingScheme
void setNumberingScheme(HcalNumberingScheme *)
Definition: HCalSD.cc:602
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:15
position
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
PHcalTB04Info::setEdep
void setEdep(double simtot, double sime, double simh, double digtot, double dige, double digh)
Definition: PHcalTB04Info.cc:113
PHcalTB04Info::setTrnsProf
void setTrnsProf(const std::vector< double > &es1, const std::vector< double > &eq1, const std::vector< double > &es2, const std::vector< double > &eq2)
Definition: PHcalTB04Info.cc:154
PHcalTB04Info::saveHit
void saveHit(int det, int lay, int eta, int phi, double e, double t)
Definition: PHcalTB04Info.cc:190
HcalTB04Analysis::secPartID
std::vector< int > secPartID
Definition: HcalTB04Analysis.h:121
HCalSD
Definition: HCalSD.h:33
HcalTestNumbering::packHcalIndex
static uint32_t packHcalIndex(int det, int z, int depth, int eta, int phi, int lay)
Definition: HcalTestNumbering.cc:7
HcalTB04Analysis::idEcal
std::vector< uint32_t > idEcal
Definition: HcalTB04Analysis.h:107
HcalTB04Analysis::eseta
std::vector< double > eseta
Definition: HcalTB04Analysis.h:115
PHcalTB04Info::setIDs
void setIDs(const std::vector< int > &, const std::vector< int > &)
Definition: PHcalTB04Info.cc:74
CaloHit::id
uint32_t id() const
Definition: CaloHit.h:28
PHcalTB04Info::setVtxSec
void setVtxSec(int id, int pdg, double px, double py, double pz, double ek)
Definition: PHcalTB04Info.cc:224
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
HcalTB04Analysis::pvUVW
G4ThreeVector pvUVW
Definition: HcalTB04Analysis.h:120
edm::Event::put
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:116
PHcalTB04Info::setEdepHcal
void setEdepHcal(const std::vector< double > &esim, const std::vector< double > &edig)
Definition: PHcalTB04Info.cc:140
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
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
HcalTB04Analysis::beamline_RM
G4RotationMatrix * beamline_RM
Definition: HcalTB04Analysis.h:101
HcalTB04Analysis::esime
std::vector< double > esime
Definition: HcalTB04Analysis.h:114
HcalTBNumberingScheme::getUnitIDs
static std::vector< uint32_t > getUnitIDs(const int type, const int mode)
Definition: HcalTBNumberingScheme.cc:75
HcalTB04Analysis::scale
double scale(int det, int layer)
Definition: HcalTB04Analysis.cc:992
HcalTB04Analysis::esphi
std::vector< double > esphi
Definition: HcalTB04Analysis.h:115
HcalTB04Analysis::timeOfFlight
double timeOfFlight(int det, int layer, double eta)
Definition: HcalTB04Analysis.cc:1005
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:89
ECalSD
Definition: ECalSD.h:26
reco::e1
Float e1
Definition: deltaR.h:22
watchdog.group
tuple group
Definition: watchdog.py:82
HcalTB04Analysis::hcalHitCache
std::vector< CaloHit > hcalHitCache
Definition: HcalTB04Analysis.h:113
CaloG4Hit::getTrackID
int getTrackID() const
Definition: CaloG4Hit.h:68
HcalTB04Analysis::xtalAnalysis
void xtalAnalysis(CLHEP::HepRandomEngine *)
Definition: HcalTB04Analysis.cc:716
HcalTBNumberingScheme::getUnitID
static uint32_t getUnitID(const uint32_t id, const int mode)
Definition: HcalTBNumberingScheme.cc:26
HcalTB04Analysis::HcalTB04Analysis
HcalTB04Analysis(const edm::ParameterSet &p)
Definition: HcalTB04Analysis.cc:57
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
HcalTB04Analysis::idXtal
std::vector< int > idXtal
Definition: HcalTB04Analysis.h:106
HcalTB04Analysis::ecalHitCache
std::vector< CaloHit > ecalHitCache
Definition: HcalTB04Analysis.h:112
sysUtil.pid
tuple pid
Definition: sysUtil.py:22
CaloHit::e
double e() const
Definition: CaloHit.h:24
HcalTB04Analysis::enois
std::vector< double > enois
Definition: HcalTB04Analysis.h:114
reco::e2
Float e2
Definition: deltaR.h:23
PHcalTB04Info::setVtxPrim
void setVtxPrim(int evNum, int type, double x, double y, double z, double u, double v, double w, double px, double py, double pz)
Definition: PHcalTB04Info.cc:204
HcalTB04Analysis::secMomentum
std::vector< G4ThreeVector > secMomentum
Definition: HcalTB04Analysis.h:122
HcalTB04Analysis::unitID
int unitID(uint32_t id)
Definition: HcalTB04Analysis.cc:980
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
HcalTB04Analysis::histo
HcalTB04Histo * histo
Definition: HcalTB04Analysis.h:92
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
HcalTB04Analysis::pvMomentum
G4ThreeVector pvMomentum
Definition: HcalTB04Analysis.h:120
CaloG4Hit::getTimeSlice
double getTimeSlice() const
Definition: CaloG4Hit.h:70
HcalTB04Analysis::~HcalTB04Analysis
virtual ~HcalTB04Analysis()
Definition: HcalTB04Analysis.cc:106
HcalTB04Analysis::idTower
std::vector< uint32_t > idTower
Definition: HcalTB04Analysis.h:107
HcalTB04Analysis::qieAnalysis
void qieAnalysis(CLHEP::HepRandomEngine *)
Definition: HcalTB04Analysis.cc:656
HcalTB04Analysis::eslay
std::vector< double > eslay
Definition: HcalTB04Analysis.h:115
CaloG4Hit::getEntry
math::XYZPoint getEntry() const
Definition: CaloG4Hit.h:50
HcalTB04Analysis::eqlay
std::vector< double > eqlay
Definition: HcalTB04Analysis.h:115
x
Definition: DDAxes.h:10
CaloG4Hit::getUnitID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:69
HcalTB04Analysis::pvPosition
G4ThreeVector pvPosition
Definition: HcalTB04Analysis.h:120
ECalSD::setNumberingScheme
void setNumberingScheme(EcalNumberingScheme *)
Definition: ECalSD.cc:266
HcalTB04Analysis::names
std::vector< std::string > names
Definition: HcalTB04Analysis.h:100
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
HcalTB04Histo::fillEdep
void fillEdep(double etots, double eecals, double ehcals, double etotq, double eecalq, double ehcalq)
Definition: HcalTB04Histo.cc:74
HcalTB04Analysis::fillEvent
void fillEvent(PHcalTB04Info &)
Definition: HcalTB04Analysis.cc:870
CaloG4Hit::getEnergyDeposit
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
HcalTB04Analysis::esimh
std::vector< double > esimh
Definition: HcalTB04Analysis.h:114
HcalTB04Histo::fillLongProf
void fillLongProf(const std::vector< double > &es, const std::vector< double > &eq)
Definition: HcalTB04Histo.cc:121
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5