CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalTB06Analysis.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HcalTestBeam
4 // Class : HcalTB06Analysis
5 //
6 // Implementation:
7 // Main analysis class for Hcal Test Beam 2006 Analysis
8 //
9 // Original Author:
10 // Created: Tue May 16 10:14:34 CEST 2006
11 // $Id: HcalTB06Analysis.cc,v 1.1 2012/07/02 04:44:40 sunanda Exp $
12 //
13 
14 // system include files
15 #include <cmath>
16 #include <iostream>
17 #include <iomanip>
18 
19 // user include files
21 
26 
27 // to retreive hits
32 
39 
40 #include "G4SDManager.hh"
41 #include "G4VProcess.hh"
42 #include "G4HCofThisEvent.hh"
43 #include "CLHEP/Units/GlobalSystemOfUnits.h"
44 #include "CLHEP/Units/GlobalPhysicalConstants.h"
45 
46 //
47 // constructors and destructor
48 //
49 
51 
52  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB06Analysis");
53  names = m_Anal.getParameter<std::vector<std::string> >("Names");
54  beamOffset =-m_Anal.getParameter<double>("BeamPosition")*cm;
55  double fMinEta = m_Anal.getParameter<double>("MinEta");
56  double fMaxEta = m_Anal.getParameter<double>("MaxEta");
57  double fMinPhi = m_Anal.getParameter<double>("MinPhi");
58  double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
59  double beamEta = (fMaxEta+fMinEta)/2.;
60  double beamPhi = (fMaxPhi+fMinPhi)/2.;
61  double beamThet= 2*atan(exp(-beamEta));
62  if (beamPhi < 0) beamPhi += twopi;
63  iceta = (int)(beamEta/0.087) + 1;
64  icphi = (int)(fabs(beamPhi)/0.087) + 5;
65  if (icphi > 72) icphi -= 73;
66 
67  produces<PHcalTB06Info>();
68 
69  beamline_RM = new G4RotationMatrix;
70  beamline_RM->rotateZ(-beamPhi);
71  beamline_RM->rotateY(-beamThet);
72 
73  edm::LogInfo("HcalTBSim") << "HcalTB06:: Initialised as observer of BeginOf"
74  << "Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
75  << " with Parameter values:\n \tbeamOffset = "
76  << beamOffset << "\ticeta = " << iceta
77  << "\ticphi = " << icphi << "\n\tbeamline_RM = "
78  << *beamline_RM;
79 
80  init();
81 
82  histo = new HcalTB06Histo(m_Anal);
83 }
84 
86 
87  edm::LogInfo("HcalTBSim") << "\n --------> Total number of selected entries"
88  << " : " << count << "\nPointers:: Histo " <<histo;
89  if (histo) {
90  delete histo;
91  histo = 0;
92  }
93 }
94 
95 //
96 // member functions
97 //
98 
100 
101  std::auto_ptr<PHcalTB06Info> product(new PHcalTB06Info);
102  fillEvent(*product);
103  e.put(product);
104 }
105 
107 
108  // counter
109  count = 0;
110  evNum = 0;
111  clear();
112 }
113 
115 
116  int irun = (*run)()->GetRunID();
117  edm::LogInfo("HcalTBSim") <<" =====> Begin of Run = " << irun;
118 
119 }
120 
122 
123  evNum = (*evt) ()->GetEventID ();
124  clear();
125  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis: =====> Begin of event = "
126  << evNum;
127 }
128 
129 void HcalTB06Analysis::update(const G4Step * aStep) {
130 
131  if (aStep != NULL) {
132  //Get Step properties
133  G4ThreeVector thePreStepPoint = aStep->GetPreStepPoint()->GetPosition();
134  G4ThreeVector thePostStepPoint;
135 
136  // Get Tracks properties
137  G4Track* aTrack = aStep->GetTrack();
138  int trackID = aTrack->GetTrackID();
139  int parentID = aTrack->GetParentID();
140  G4ThreeVector position = aTrack->GetPosition();
141  G4ThreeVector momentum = aTrack->GetMomentum();
142  G4String partType = aTrack->GetDefinition()->GetParticleType();
143  G4String partSubType = aTrack->GetDefinition()->GetParticleSubType();
144  int partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
145 #ifdef ddebug
146  bool isPDGStable = aTrack->GetDefinition()->GetPDGStable();
147 #endif
148  double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
149  double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
150 
151  if (!pvFound) { //search for v1
152  double stepDeltaEnergy = aStep->GetDeltaEnergy ();
153  double kinEnergy = aTrack->GetKineticEnergy ();
154 
155  // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
156  if (trackID == 1 && parentID == 0 &&
157  ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
158  pvType = -1;
159  if (kinEnergy == 0.) {
160  pvType = 0;
161  } else {
162  if (fabs (stepDeltaEnergy / kinEnergy) > 0.1) pvType = 1;
163  }
164  pvFound = true;
166  pvMomentum = momentum;
167  // Rotated coord.system:
168  pvUVW = (*beamline_RM)*(pvPosition);
169 
170  //Volume name requires some checks:
171  G4String thePostPVname = "NoName";
172  G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
173  if (thePostPoint) {
174  thePostStepPoint = thePostPoint->GetPosition();
175  G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
176  if (thePostPV) thePostPVname = thePostPV->GetName ();
177  }
178 #ifdef ddebug
179  LogDebug("HcalTBSim") << "HcalTB06Analysis:: V1 found at: "
180  << thePostStepPoint << " G4VPhysicalVolume: "
181  << thePostPVname;
182 #endif
183  LogDebug("HcalTBSim") << "HcalTB06Analysis::fill_v1Pos: Primary Track "
184  << "momentum: " << pvMomentum << " psoition "
185  << pvPosition << " u/v/w " << pvUVW;
186  }
187  } else {
188  // watch for secondaries originating @v1, including the surviving primary
189  if ((trackID != 1 && parentID == 1 &&
190  (aTrack->GetCurrentStepNumber () == 1) &&
191  (thePreStepPoint == pvPosition)) ||
192  (trackID == 1 && thePreStepPoint == pvPosition)) {
193 #ifdef ddebug
194  LogDebug("HcalTBSim") << "HcalTB06Analysis::A secondary... PDG:"
195  << partPDGEncoding << " TrackID:" << trackID
196  << " ParentID:" << parentID << " stable: "
197  << isPDGStable << " Tau: " << pDGlifetime
198  << " cTauGamma="
199  << c_light*pDGlifetime*gammaFactor*1000.
200  << "um" << " GammaFactor: " << gammaFactor;
201 #endif
202  secTrackID.push_back(trackID);
203  secPartID.push_back(partPDGEncoding);
204  secMomentum.push_back(momentum);
205  secEkin.push_back(aTrack->GetKineticEnergy());
206 
207  // Check for short-lived secondaries: cTauGamma<100um
208  double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
209  if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {//short-lived secondary
210  shortLivedSecondaries.push_back(trackID);
211  } else {//normal secondary - enter into the V1-calorimetric tree
212  // histos->fill_v1cSec (aTrack);
213  }
214  }
215  // Also watch for tertiary particles coming from
216  // short-lived secondaries from V1
217  if (aTrack->GetCurrentStepNumber() == 1) {
218  if (shortLivedSecondaries.size() > 0) {
219  int pid = parentID;
220  std::vector<int>::iterator pos1= shortLivedSecondaries.begin();
221  std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
222  std::vector<int>::iterator pos;
223  for (pos = pos1; pos != pos2; pos++) {
224  if (*pos == pid) {//ParentID is on the list of short-lived
225  // secondary
226 #ifdef ddebug
227  LogDebug("HcalTBSim") << "HcalTB06Analysis:: A tertiary... PDG:"
228  << partPDGEncoding << " TrackID:" <<trackID
229  << " ParentID:" << parentID << " stable: "
230  << isPDGStable << " Tau: " << pDGlifetime
231  << " cTauGamma="
232  << c_light*pDGlifetime*gammaFactor*1000.
233  << "um GammaFactor: " << gammaFactor;
234 #endif
235  }
236  }
237  }
238  }
239  }
240  }
241 }
242 
244 
245  count++;
246 
247  //fill the buffer
248  LogDebug("HcalTBSim") << "HcalTB06Analysis::Fill event "
249  << (*evt)()->GetEventID();
250  fillBuffer (evt);
251 
252  //Final Analysis
253  LogDebug("HcalTBSim") << "HcalTB06Analysis::Final analysis";
254  finalAnalysis();
255 
256  int iEvt = (*evt)()->GetEventID();
257  if (iEvt < 10)
258  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
259  else if ((iEvt < 100) && (iEvt%10 == 0))
260  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
261  else if ((iEvt < 1000) && (iEvt%100 == 0))
262  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
263  else if ((iEvt < 10000) && (iEvt%1000 == 0))
264  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
265 }
266 
268 
269  std::vector<CaloHit> hhits;
270  int idHC, j;
271  CaloG4HitCollection* theHC;
272  std::map<int,float,std::less<int> > primaries;
273  double etot1=0, etot2=0;
274 
275  // Look for the Hit Collection of HCal
276  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
277  std::string sdName = names[0];
278  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
279  theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
280  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hit Collection for " << sdName
281  << " of ID " << idHC << " is obtained at " << theHC;
282 
283  if (idHC >= 0 && theHC > 0) {
284  hhits.reserve(theHC->entries());
285  for (j = 0; j < theHC->entries(); j++) {
286  CaloG4Hit* aHit = (*theHC)[j];
287  double e = aHit->getEnergyDeposit()/GeV;
288  double time = aHit->getTimeSlice();
289  math::XYZPoint pos = aHit->getEntry();
290  unsigned int id = aHit->getUnitID();
291  double theta = pos.theta();
292  double eta = -log(tan(theta * 0.5));
293  double phi = pos.phi();
294  CaloHit hit(2,1,e,eta,phi,time,id);
295  hhits.push_back(hit);
296  primaries[aHit->getTrackID()]+= e;
297  etot1 += e;
298 #ifdef ddebug
299  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hcal Hit i/p " << j
300  << " ID 0x" << std::hex << id << std::dec
301  << " time " << std::setw(6) << time << " theta "
302  << std::setw(8) << theta << " eta " << std::setw(8)
303  << eta << " phi " << std::setw(8) << phi << " e "
304  << std::setw(8) << e;
305 #endif
306  }
307  }
308 
309  // Add hits in the same channel within same time slice
310  std::vector<CaloHit>::iterator itr;
311  int nHit = hhits.size();
312  std::vector<CaloHit*> hits(nHit);
313  for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
314  hits[j] = &hhits[j];
315  }
316  sort(hits.begin(),hits.end(),CaloHitIdMore());
317  std::vector<CaloHit*>::iterator k1, k2;
318  int nhit = 0;
319  for (k1 = hits.begin(); k1 != hits.end(); k1++) {
320  int det = (**k1).det();
321  int layer = (**k1).layer();
322  double ehit = (**k1).e();
323  double eta = (**k1).eta();
324  double phi = (**k1).phi();
325  double jitter = (**k1).t();
326  uint32_t unitID = (**k1).id();
327  int jump = 0;
328  for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
329  unitID==(**k2).id(); k2++) {
330  ehit += (**k2).e();
331  jump++;
332  }
333  nhit++;
334  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
335  hcalHitCache.push_back(hit);
336  etot2 += ehit;
337  k1 += jump;
338 #ifdef ddebug
339  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hcal Hit store " << nhit
340  << " ID 0x" << std::hex << unitID << std::dec
341  << " time " << std::setw(6) << jitter << " eta "
342  << std::setw(8) << eta << " phi " << std::setw(8)
343  << phi << " e " << std::setw(8) << ehit;
344 #endif
345  }
346  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Stores " << nhit << " HCal hits"
347  << " from " << nHit << " input hits E(Hcal) " << etot1
348  << " " << etot2;
349 
350  // Look for the Hit Collection of ECal
351  std::vector<CaloHit> ehits;
352  sdName= names[1];
353  idHC = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
354  theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
355  etot1 = etot2 = 0;
356  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hit Collection for " << sdName
357  << " of ID " << idHC << " is obtained at " << theHC;
358  if (idHC >= 0 && theHC > 0) {
359  ehits.reserve(theHC->entries());
360  for (j = 0; j < theHC->entries(); j++) {
361  CaloG4Hit* aHit = (*theHC)[j];
362  double e = aHit->getEnergyDeposit()/GeV;
363  double time = aHit->getTimeSlice();
364  math::XYZPoint pos = aHit->getEntry();
365  unsigned int id = aHit->getUnitID();
366  double theta = pos.theta();
367  double eta = -log(tan(theta * 0.5));
368  double phi = pos.phi();
369  if (e < 0 || e > 100000.) e = 0;
370  CaloHit hit(1,0,e,eta,phi,time,id);
371  ehits.push_back(hit);
372  primaries[aHit->getTrackID()]+= e;
373  etot1 += e;
374 #ifdef ddebug
375  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Ecal Hit i/p " << j
376  << " ID 0x" << std::hex << id << std::dec
377  << " time " << std::setw(6) << time << " theta "
378  << std::setw(8) << theta << " eta " <<std::setw(8)
379  << eta << " phi " << std::setw(8) << phi << " e "
380  << std::setw(8) << e;
381 #endif
382  }
383  }
384 
385  // Add hits in the same channel within same time slice
386  nHit = ehits.size();
387  std::vector<CaloHit*> hite(nHit);
388  for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
389  hite[j] = &ehits[j];
390  }
391  sort(hite.begin(),hite.end(),CaloHitIdMore());
392  nhit = 0;
393  for (k1 = hite.begin(); k1 != hite.end(); k1++) {
394  int det = (**k1).det();
395  int layer = (**k1).layer();
396  double ehit = (**k1).e();
397  double eta = (**k1).eta();
398  double phi = (**k1).phi();
399  double jitter = (**k1).t();
400  uint32_t unitID = (**k1).id();
401  int jump = 0;
402  for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
403  unitID==(**k2).id(); k2++) {
404  ehit += (**k2).e();
405  jump++;
406  }
407  nhit++;
408  CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
409  ecalHitCache.push_back(hit);
410  etot2 += ehit;
411  k1 += jump;
412 #ifdef ddebug
413  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Ecal Hit store " << nhit
414  << " ID 0x" << std::hex << unitID << std::dec
415  << " time " << std::setw(6) << jitter << " eta "
416  << std::setw(8) << eta << " phi " << std::setw(8)
417  << phi << " e " << std::setw(8) << ehit;
418 #endif
419  }
420  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Stores " << nhit << " ECal hits"
421  << " from " << nHit << " input hits E(Ecal) " << etot1
422  << " " << etot2;
423 
424  // Find Primary info:
425  nPrimary = (int)(primaries.size());
426  int trackID = 0;
427  G4PrimaryParticle* thePrim=0;
428  int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
429  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Event has " << nvertex
430  << " verteices";
431  if (nvertex<=0)
432  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERROR: no "
433  << "vertex found for event " << evNum;
434 
435  for (int i = 0 ; i<nvertex; i++) {
436  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
437  if (avertex == 0) {
438  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERR: pointer "
439  << "to vertex = 0 for event " << evNum;
440  } else {
441  LogDebug("HcalTBSim") << "HcalTB06Analysis::Vertex number :" << i << " "
442  << avertex->GetPosition();
443  int npart = avertex->GetNumberOfParticle();
444  if (npart == 0)
445  edm::LogWarning("HcalTBSim") << "HcalTB06Analysis::End Of Event ERR: "
446  << "no primary!";
447  if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
448  }
449  }
450 
451  if (thePrim != 0) {
452  double px = thePrim->GetPx();
453  double py = thePrim->GetPy();
454  double pz = thePrim->GetPz();
455  double p = std::sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
456  pInit = p/GeV;
457  if (p==0)
458  edm::LogWarning("HcalTBSim") << "HcalTB06Analysis:: EndOfEvent ERR: "
459  << "primary has p=0 ";
460  else {
461  double costheta = pz/p;
462  double theta = acos(std::min(std::max(costheta,-1.),1.));
463  etaInit = -log(tan(theta/2));
464  if (px != 0 || py != 0) phiInit = atan2(py,px);
465  }
466  particleType = thePrim->GetPDGcode();
467  } else
468  edm::LogWarning("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERR: could "
469  << "not find primary";
470 
471 }
472 
474 
475  //Beam Information
477 
478  // Total Energy
479  eecals = ehcals = 0.;
480  for (unsigned int i=0; i<hcalHitCache.size(); i++) {
481  ehcals += hcalHitCache[i].e();
482  }
483  for (unsigned int i=0; i<ecalHitCache.size(); i++) {
484  eecals += ecalHitCache[i].e();
485  }
486  etots = eecals + ehcals;
487  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Energy deposit at Sim Level "
488  << "(Total) " << etots << " (ECal) " << eecals
489  << " (HCal) " << ehcals;
490  histo->fillEdep(etots, eecals, ehcals);
491 }
492 
493 
495 
496  //Beam Information
498 
499  // Total Energy
500  product.setEdep(etots, eecals, ehcals);
501 
502  //Energy deposits in the crystals and towers
503  for (unsigned int i=0; i<hcalHitCache.size(); i++)
504  product.saveHit(hcalHitCache[i].id(), hcalHitCache[i].eta(),
505  hcalHitCache[i].phi(), hcalHitCache[i].e(),
506  hcalHitCache[i].t());
507  for (unsigned int i=0; i<ecalHitCache.size(); i++)
508  product.saveHit(ecalHitCache[i].id(), ecalHitCache[i].eta(),
509  ecalHitCache[i].phi(), ecalHitCache[i].e(),
510  ecalHitCache[i].t());
511 
512  //Vertex associated quantities
513  product.setVtxPrim(evNum, pvType, pvPosition.x(), pvPosition.y(),
514  pvPosition.z(), pvUVW.x(), pvUVW.y(), pvUVW.z(),
515  pvMomentum.x(), pvMomentum.y(), pvMomentum.z());
516  for (unsigned int i = 0; i < secTrackID.size(); i++) {
517  product.setVtxSec(secTrackID[i], secPartID[i], secMomentum[i].x(),
518  secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
519  }
520 }
521 
523 
524  pvFound = false;
525  pvType =-2;
526  pvPosition = G4ThreeVector();
527  pvMomentum = G4ThreeVector();
528  pvUVW = G4ThreeVector();
529  secTrackID.clear();
530  secPartID.clear();
531  secMomentum.clear();
532  secEkin.clear();
533  shortLivedSecondaries.clear();
534 
535  ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end());
536  hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end());
537  nPrimary = particleType = 0;
538  pInit = etaInit = phiInit = 0;
539 }
540 
#define LogDebug(id)
T getParameter(std::string const &) const
G4ThreeVector pvMomentum
int i
Definition: DBlmapReader.cc:9
#define DEFINE_SIMWATCHER(type)
void update(const BeginOfRun *run)
This routine will be called when the appropriate signal arrives.
std::vector< std::string > names
std::vector< CaloHit > ecalHitCache
void setEdep(double simtot, double sime, double simh)
void fillEdep(double etots, double eecals, double ehcals)
Geom::Theta< T > theta() const
#define NULL
Definition: scimark2.h:8
double npart
Definition: HydjetWrapper.h:45
#define min(a, b)
Definition: mlp_lapack.h:161
HcalTB06Histo * histo
T eta() const
void setVtxPrim(int evNum, int type, double x, double y, double z, double u, double v, double w, double px, double py, double pz)
float float float z
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
G4ThreeVector pvPosition
virtual void produce(edm::Event &, const edm::EventSetup &)
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:94
T sqrt(T t)
Definition: SSEVec.h:48
void fillEvent(PHcalTB06Info &)
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
int j
Definition: DBlmapReader.cc:9
void setVtxSec(int id, int pdg, double px, double py, double pz, double ek)
void saveHit(unsigned int det, double eta, double phi, double e, double t)
int getTrackID() const
Definition: CaloG4Hit.h:68
std::vector< CaloHit > hcalHitCache
std::vector< G4ThreeVector > secMomentum
std::vector< int > shortLivedSecondaries
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
void setPrimary(int primary, int id, double energy, double eta, double phi)
void fillBuffer(const EndOfEvent *evt)
HcalTB06Analysis(const edm::ParameterSet &p)
std::vector< double > secEkin
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
double getTimeSlice() const
Definition: CaloG4Hit.h:70
math::XYZPoint getEntry() const
Definition: CaloG4Hit.h:50
Definition: DDAxes.h:10
std::vector< int > secPartID
uint32_t getUnitID() const
Definition: CaloG4Hit.h:69
void fillPrimary(double energy, double eta, double phi)
virtual ~HcalTB06Analysis()
std::vector< int > secTrackID
G4RotationMatrix * beamline_RM
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
G4ThreeVector pvUVW
double getEnergyDeposit() const
Definition: CaloG4Hit.h:81
Definition: DDAxes.h:10