CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MaterialBudgetData.cc
Go to the documentation of this file.
2 
4 
5 #include "G4Step.hh"
6 #include "G4Material.hh"
7 #include "G4EventManager.hh"
8 #include "G4Event.hh"
9 
10 //#define TREE_DEBUG
11 
13 {
14  //instantiate categorizer to assign an ID to volumes and materials
16  allStepsToTree = false;
17  densityConvertionFactor = 6.24E18;
18 }
19 
21 }
22 
24 {
25  allStepsToTree = true;
26  MAXNUMBERSTEPS = 0;
27  MAXNUMBERSTEPS = 10000;
28  theDmb = new float[MAXNUMBERSTEPS];
29  theDil = new float[MAXNUMBERSTEPS];
30  // rr
31  theSupportDmb = new float[MAXNUMBERSTEPS];
32  theSensitiveDmb = new float[MAXNUMBERSTEPS];
33  theCablesDmb = new float[MAXNUMBERSTEPS];
34  theCoolingDmb = new float[MAXNUMBERSTEPS];
35  theElectronicsDmb = new float[MAXNUMBERSTEPS];
36  theOtherDmb = new float[MAXNUMBERSTEPS];
37  theAirDmb = new float[MAXNUMBERSTEPS];
38  theSupportDil = new float[MAXNUMBERSTEPS];
39  theSensitiveDil = new float[MAXNUMBERSTEPS];
40  theCablesDil = new float[MAXNUMBERSTEPS];
41  theCoolingDil = new float[MAXNUMBERSTEPS];
42  theElectronicsDil = new float[MAXNUMBERSTEPS];
43  theOtherDil = new float[MAXNUMBERSTEPS];
44  theAirDil = new float[MAXNUMBERSTEPS];
45  // rr
46  theInitialX = new double[MAXNUMBERSTEPS];
47  theInitialY = new double[MAXNUMBERSTEPS];
48  theInitialZ = new double[MAXNUMBERSTEPS];
49  theFinalX = new double[MAXNUMBERSTEPS];
50  theFinalY = new double[MAXNUMBERSTEPS];
51  theFinalZ = new double[MAXNUMBERSTEPS];
52  // rr
53  theVolumeID = new int[MAXNUMBERSTEPS];
55  theVolumeCopy = new int[MAXNUMBERSTEPS];
56  theVolumeX = new float[MAXNUMBERSTEPS];
57  theVolumeY = new float[MAXNUMBERSTEPS];
58  theVolumeZ = new float[MAXNUMBERSTEPS];
59  theVolumeXaxis1 = new float[MAXNUMBERSTEPS];
60  theVolumeXaxis2 = new float[MAXNUMBERSTEPS];
61  theVolumeXaxis3 = new float[MAXNUMBERSTEPS];
62  theVolumeYaxis1 = new float[MAXNUMBERSTEPS];
63  theVolumeYaxis2 = new float[MAXNUMBERSTEPS];
64  theVolumeYaxis3 = new float[MAXNUMBERSTEPS];
65  theVolumeZaxis1 = new float[MAXNUMBERSTEPS];
66  theVolumeZaxis2 = new float[MAXNUMBERSTEPS];
67  theVolumeZaxis3 = new float[MAXNUMBERSTEPS];
68  theMaterialID = new int[MAXNUMBERSTEPS];
70  theMaterialX0 = new float[MAXNUMBERSTEPS];
73  theStepID = new int[MAXNUMBERSTEPS];
74  theStepInitialPt = new float[MAXNUMBERSTEPS];
75  theStepInitialEta = new float[MAXNUMBERSTEPS];
76  theStepInitialPhi = new float[MAXNUMBERSTEPS];
78  theStepInitialPx = new float[MAXNUMBERSTEPS];
79  theStepInitialPy = new float[MAXNUMBERSTEPS];
80  theStepInitialPz = new float[MAXNUMBERSTEPS];
84  theStepFinalPt = new float[MAXNUMBERSTEPS];
85  theStepFinalEta = new float[MAXNUMBERSTEPS];
86  theStepFinalPhi = new float[MAXNUMBERSTEPS];
88  theStepFinalPx = new float[MAXNUMBERSTEPS];
89  theStepFinalPy = new float[MAXNUMBERSTEPS];
90  theStepFinalPz = new float[MAXNUMBERSTEPS];
91  theStepFinalBeta = new float[MAXNUMBERSTEPS];
92  theStepFinalGamma = new float[MAXNUMBERSTEPS];
93  theStepFinalMass = new float[MAXNUMBERSTEPS];
96  // rr
97 }
98 
99 
100 void MaterialBudgetData::dataStartTrack( const G4Track* aTrack )
101 {
102  const G4ThreeVector& dir = aTrack->GetMomentum() ;
103 
105 
106  theStepN=0;
107  theTotalMB=0;
108  theTotalIL=0;
109  theEta=0;
110  thePhi=0;
111 
112  // rr
113  theID=0;
114  thePt=0;
115  theEnergy=0;
116  theMass=0;
117 
118  theSupportMB = 0.;
119  theSensitiveMB = 0.;
120  theCablesMB = 0.;
121  theCoolingMB = 0.;
122  theElectronicsMB = 0.;
123  theOtherMB = 0.;
124  theAirMB = 0.;
125  theSupportIL = 0.;
126  theSensitiveIL = 0.;
127  theCablesIL = 0.;
128  theCoolingIL = 0.;
129  theElectronicsIL = 0.;
130  theOtherIL = 0.;
131  theAirIL = 0.;
134  theCablesFractionMB = 0.;
137  theOtherFractionMB = 0.;
138  theAirFractionMB = 0.;
141  theCablesFractionIL = 0.;
144  theOtherFractionIL = 0.;
145  theAirFractionIL = 0.;
146  // rr
147 
148  theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
149  thePt = dir.perp();
150  if( dir.theta() != 0 ) {
151  theEta = dir.eta();
152  } else {
153  theEta = -99;
154  }
155  // thePhi = dir.phi()/deg; // better not to store in deg
156  thePhi = dir.phi();
157  theEnergy = aTrack->GetTotalEnergy();
158  theMass = aTrack->GetDefinition()->GetPDGMass();
159 
160 }
161 
162 
163 void MaterialBudgetData::dataEndTrack( const G4Track* aTrack )
164 {
165  //- std::cout << "[OVAL] MaterialBudget " << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " " << theEta << " " << thePhi << " " << theTotalMB << std::endl;
166  // rr
167  std::cout << "Recorded steps " << theStepN << std::endl;
168  std::cout << " Material Budget: Radiation Length " << "G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()" << " eta " << theEta << " phi " << thePhi << " total X " << theTotalMB << " SUP " << theSupportMB << " SEN " << theSensitiveMB << " CAB " << theCablesMB << " COL " << theCoolingMB << " ELE " << theElectronicsMB << " other " << theOtherMB << " Air " << theAirMB << std::endl;
169  std::cout << " Material Budget: Interaction Length " << "G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()" << " eta " << theEta << " phi " << thePhi << " total L " << theTotalIL << " SUP " << theSupportIL << " SEN " << theSensitiveIL << " CAB " << theCablesIL << " COL " << theCoolingIL << " ELE " << theElectronicsIL << " other " << theOtherIL << " Air " << theAirIL << std::endl;
170  // rr
171 }
172 
173 void MaterialBudgetData::dataPerStep( const G4Step* aStep )
174 {
175  assert(aStep);
176  G4StepPoint* prePoint = aStep->GetPreStepPoint();
177  G4StepPoint* postPoint = aStep->GetPostStepPoint();
178  assert(prePoint);
179  assert(postPoint);
180  G4Material * theMaterialPre = prePoint->GetMaterial();
181  assert(theMaterialPre);
182 
183  CLHEP::Hep3Vector prePos = prePoint->GetPosition();
184  CLHEP::Hep3Vector postPos = postPoint->GetPosition();
185 
186  G4double steplen = aStep->GetStepLength();
187 
188  G4double radlen = theMaterialPre->GetRadlen();
189  G4double intlen = theMaterialPre->GetNuclearInterLength();
190  G4double density = theMaterialPre->GetDensity() / densityConvertionFactor; // always g/cm3
191 
192  G4String materialName = theMaterialPre->GetName();
193  std::cout << " steplen " << steplen << " radlen " << radlen << " mb " << steplen/radlen << " mate " << theMaterialPre->GetName() << std::endl;
194 
195  G4String volumeName = aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0)->GetLogicalVolume()->GetName();
196  std::cout << " Volume " << volumeName << "\n";
197  std::cout << " Material " << materialName << "\n";
198 
199  // instantiate the categorizer
201  int volumeID = myMaterialBudgetCategorizer->volume( volumeName );
202  int materialID = myMaterialBudgetCategorizer->material( materialName );
203  std::cout << "Volume ID " << volumeID << " and material ID " << materialID << "\n";
204 
205  // FIXME: Both volume ID and material ID are zeros, so this part is not executed leaving all
206  // values as zeros.
207  if( volumeID != 0 )
208  {
216 
217  if(theOtherFractionMB!=0) std::cout << " material found with no category " << materialName
218  << " in volume " << volumeName << std::endl;
226  if(theOtherFractionIL!=0) std::cout << " material found with no category " << materialName
227  << " in volume " << volumeName << std::endl;
228  }
229  // if(theOtherFractionMB!=0) LogDebug("MaterialBudgetData") << " material found with no category " << name
230  // << " in volume " << lv->GetName();
231  // rr
232 
233  float dmb = steplen/radlen;
234  float dil = steplen/intlen;
235 
236  G4VPhysicalVolume* pv = aStep->GetPreStepPoint()->GetPhysicalVolume();
237  const G4VTouchable* t = aStep->GetPreStepPoint()->GetTouchable();
238  G4ThreeVector objectTranslation = t->GetTranslation();
239  const G4RotationMatrix* objectRotation = t->GetRotation();
240  const G4VProcess* interactionPre = prePoint->GetProcessDefinedStep();
241  const G4VProcess* interactionPost = postPoint->GetProcessDefinedStep();
242 
243  G4Track* track = aStep->GetTrack();
244  if(theStepN==0) std::cout << " Simulated Particle " << theID << "\tMass " << theMass << " MeV/c2"
245  << "\tPt = " << thePt << " MeV/c" << "\tEta = " << theEta << "\tPhi = " << thePhi
246  << "\tEnergy = " << theEnergy << " MeV"
247  // << std::endl
248  // << "\tMagnetic Field at (0,0,0): (" << B000[0] << "," < B000[1] << "," << B000[2] << ")"
249  << std::endl;
250 
251  //fill data per step
252  if( allStepsToTree ){
254  theDmb[theStepN] = dmb;
255  theDil[theStepN] = dil;
270  theInitialX[theStepN] = prePos.x();
271  theInitialY[theStepN] = prePos.y();
272  theInitialZ[theStepN] = prePos.z();
273  theFinalX[theStepN] = postPos.x();
274  theFinalY[theStepN] = postPos.y();
275  theFinalZ[theStepN] = postPos.z();
276  theVolumeID[theStepN] = volumeID;
277  theVolumeName[theStepN] = volumeName;
278  theVolumeCopy[theStepN] = pv->GetCopyNo();
279  theVolumeX[theStepN] = objectTranslation.x();
280  theVolumeY[theStepN] = objectTranslation.y();
281  theVolumeZ[theStepN] = objectTranslation.z();
282  theVolumeXaxis1[theStepN] = objectRotation->xx();
283  theVolumeXaxis2[theStepN] = objectRotation->xy();
284  theVolumeXaxis3[theStepN] = objectRotation->xz();
285  theVolumeYaxis1[theStepN] = objectRotation->yx();
286  theVolumeYaxis2[theStepN] = objectRotation->yy();
287  theVolumeYaxis3[theStepN] = objectRotation->yz();
288  theVolumeZaxis1[theStepN] = objectRotation->zx();
289  theVolumeZaxis2[theStepN] = objectRotation->zy();
290  theVolumeZaxis3[theStepN] = objectRotation->zz();
291  theMaterialID[theStepN] = materialID;
292  theMaterialName[theStepN] = materialName;
293  theMaterialX0[theStepN] = radlen;
294  theMaterialLambda0[theStepN] = intlen;
295  theMaterialDensity[theStepN] = density;
296  theStepID[theStepN] = track->GetDefinition()->GetPDGEncoding();
297  theStepInitialPt[theStepN] = prePoint->GetMomentum().perp();
298  theStepInitialEta[theStepN] = prePoint->GetMomentum().eta();
299  theStepInitialPhi[theStepN] = prePoint->GetMomentum().phi();
300  theStepInitialEnergy[theStepN] = prePoint->GetTotalEnergy();
301  theStepInitialPx[theStepN] = prePoint->GetMomentum().x();
302  theStepInitialPy[theStepN] = prePoint->GetMomentum().y();
303  theStepInitialPz[theStepN] = prePoint->GetMomentum().z();
304  theStepInitialBeta[theStepN] = prePoint->GetBeta();
305  theStepInitialGamma[theStepN] = prePoint->GetGamma();
306  theStepInitialMass[theStepN] = prePoint->GetMass();
307  theStepFinalPt[theStepN] = postPoint->GetMomentum().perp();
308  theStepFinalEta[theStepN] = postPoint->GetMomentum().eta();
309  theStepFinalPhi[theStepN] = postPoint->GetMomentum().phi();
310  theStepFinalEnergy[theStepN] = postPoint->GetTotalEnergy();
311  theStepFinalPx[theStepN] = postPoint->GetMomentum().x();
312  theStepFinalPy[theStepN] = postPoint->GetMomentum().y();
313  theStepFinalPz[theStepN] = postPoint->GetMomentum().z();
314  theStepFinalBeta[theStepN] = postPoint->GetBeta();
315  theStepFinalGamma[theStepN] = postPoint->GetGamma();
316  theStepFinalMass[theStepN] = postPoint->GetMass();
317  int preProcType = -99;
318  int postProcType = -99;
319  if (interactionPre) preProcType = interactionPre->GetProcessType();
320  theStepPreProcess[theStepN] = preProcType;
321  if (interactionPost) postProcType = interactionPost->GetProcessType();
322  theStepPostProcess[theStepN] = postProcType;
323 #ifdef TREE_DEBUG
324  std::cout << " step " << theStepN
325  << "\tDelta MB = " << theDmb[theStepN]
326  << std::endl
327  << "\t\tSupport " << theSupportDmb[theStepN]
328  << " Sensitive " << theSensitiveDmb[theStepN]
329  << " Cables " << theCablesDmb[theStepN]
330  << " Cooling " << theCoolingDmb[theStepN]
331  << " Electronics " << theElectronicsDmb[theStepN]
332  << " Other " << theOtherDmb[theStepN]
333  << " Air " << theAirDmb[theStepN]
334  << std::endl
335  << "\tDelta IL = " << theDil[theStepN]
336  << std::endl
337  << "\t\tSupport " << theSupportDil[theStepN]
338  << " Sensitive " << theSensitiveDil[theStepN]
339  << " Cables " << theCablesDil[theStepN]
340  << " Cooling " << theCoolingDil[theStepN]
341  << " Electronics " << theElectronicsDil[theStepN]
342  << " Other " << theOtherDil[theStepN]
343  << " Air " << theAirDil[theStepN]
344  << std::endl;
345  if (interactionPre)
346  std::cout << "\tProcess Pre " << interactionPre->GetProcessName()
347  << " type " << theStepPreProcess[theStepN] << " " << interactionPre->GetProcessTypeName(G4ProcessType(theStepPreProcess[theStepN]))
348  << std::endl;
349  if (interactionPost)
350  std::cout << "\tProcess Post " << interactionPost->GetProcessName()
351  << " type " << theStepPostProcess[theStepN] << " "
352  << interactionPost->GetProcessTypeName(G4ProcessType(theStepPostProcess[theStepN]))
353  << std::endl;
354  std::cout << "\tPre x = " << theInitialX[theStepN]
355  << "\ty = " << theInitialY[theStepN]
356  << "\tz = " << theInitialZ[theStepN]
357  << "\tPolar Radius = " << sqrt(prePos.x()*prePos.x()+prePos.y()*prePos.y())
358  << "\tPt = " << theStepInitialPt[theStepN]
359  << "\tEnergy = " << theStepInitialEnergy[theStepN]
360  // << std::endl
361  // << "B-field(T) at Pre (cm): " << field->fieldInTesla(GlobalPoint(pos.x()/10.,pos.y()/10.,pos.z()/10.))
362  << std::endl;
363  std::cout << "\tPost x = " << theFinalX[theStepN]
364  << "\ty = " << theFinalY[theStepN]
365  << "\tz = " << theFinalZ[theStepN]
366  << "\tPolar Radius = " << sqrt(postPos.x()*postPos.x()+postPos.y()*postPos.y())
367  << "\tPt = " << theStepFinalPt[theStepN]
368  << "\tEnergy = " << theStepFinalEnergy[theStepN]
369  << std::endl;
370  std::cout << "\tvolume " << volumeID << " " << theVolumeName[theStepN]
371  << " copy number " << theVolumeCopy[theStepN]
372  << "\tmaterial " << theMaterialID[theStepN] << " " << theMaterialName[theStepN]
373  << "\tDensity = " << theMaterialDensity[theStepN] << " g/cm3"
374  << "\tX0 = " << theMaterialX0[theStepN] << " mm"
375  << "\tLambda0 = " << theMaterialLambda0[theStepN] << " mm"
376  << std::endl;
377  std::cout << "\t\tParticle " << theStepID[theStepN]
378  << " Initial Pt = " << theStepInitialPt[theStepN] << " MeV/c"
379  << " eta = " << theStepInitialEta[theStepN]
380  << " phi = " << theStepInitialPhi[theStepN]
381  << " Energy = " << theStepInitialEnergy[theStepN] << " MeV"
382  << " Mass = " << theStepInitialMass[theStepN] << " MeV/c2"
383  << " Beta = " << theStepInitialBeta[theStepN]
384  << " Gamma = " << theStepInitialGamma[theStepN]
385  << std::endl
386  << "\t\tParticle " << theStepID[theStepN]
387  << " Final Pt = " << theStepFinalPt[theStepN] << " MeV/c"
388  << " eta = " << theStepFinalEta[theStepN]
389  << " phi = " << theStepFinalPhi[theStepN]
390  << " Energy = " << theStepFinalEnergy[theStepN] << " MeV"
391  << " Mass = " << theStepFinalMass[theStepN] << " MeV/c2"
392  << " Beta = " << theStepFinalBeta[theStepN]
393  << " Gamma = " << theStepFinalGamma[theStepN]
394  << std::endl;
395  std:: cout << "\tVolume Centre x = " << theVolumeX[theStepN]
396  << "\ty = " << theVolumeY[theStepN]
397  << "\tz = " << theVolumeZ[theStepN]
398  << "\tPolar Radius = " << sqrt( theVolumeX[theStepN]*theVolumeX[theStepN] +
399  theVolumeY[theStepN]*theVolumeY[theStepN] )
400  << std::endl;
401  std::cout << "\tx axis = ("
402  << theVolumeXaxis1[theStepN] << ","
403  << theVolumeXaxis2[theStepN] << ","
404  << theVolumeXaxis3[theStepN] << ")"
405  << std::endl;
406  std::cout << "\ty axis = ("
407  << theVolumeYaxis1[theStepN] << ","
408  << theVolumeYaxis2[theStepN] << ","
409  << theVolumeYaxis3[theStepN] << ")"
410  << std::endl;
411  std::cout << "\tz axis = ("
412  << theVolumeZaxis1[theStepN] << ","
413  << theVolumeZaxis2[theStepN] << ","
414  << theVolumeZaxis3[theStepN] << ")"
415  << std::endl;
416  std::cout << "\tSecondaries"
417  << std::endl;
418  for(G4TrackVector::iterator iSec = aStep->GetSecondary()->begin(); iSec!=aStep->GetSecondary()->end(); iSec++) {
419  std::cout << "\t\tid " << (*iSec)->GetDefinition()->GetPDGEncoding()
420  << " created through process "
421  << " type " << (*iSec)->GetCreatorProcess()->GetProcessType()
422  << " " << (*iSec)->GetCreatorProcess()->GetProcessTypeName(G4ProcessType((*iSec)->GetCreatorProcess()->GetProcessType()))
423  << std::endl;
424  }
425 #endif
426  }
427 
428  theTrkLen = aStep->GetTrack()->GetTrackLength();
429  //- std::cout << " theTrkLen " << theTrkLen << " theTrkLen2 " << theTrkLen2 << " postPos " << postPos.mag() << postPos << std::endl;
430  thePVname = pv->GetName();
431  thePVcopyNo = pv->GetCopyNo();
432  theRadLen = radlen;
433  theIntLen = intlen;
434  theTotalMB += dmb;
435  theTotalIL += dil;
436 
437  // rr
440  theCablesMB += (dmb * theCablesFractionMB);
443  theOtherMB += (dmb * theOtherFractionMB);
444  theAirMB += (dmb * theAirFractionMB);
447  theCablesIL += (dil * theCablesFractionIL);
450  theOtherIL += (dil * theOtherFractionIL);
451  theAirIL += (dil * theAirFractionIL);
452  // rr
453 
454  theStepN++;
455 
456 }
457 
458 
tuple t
Definition: tree.py:139
std::vector< float > x0fraction(std::string s)
assert(m_qm.get())
std::string * theMaterialName
T sqrt(T t)
Definition: SSEVec.h:48
void dataEndTrack(const G4Track *aTrack)
void dataPerStep(const G4Step *aStep)
MaterialBudgetCategorizer * myMaterialBudgetCategorizer
void dataStartTrack(const G4Track *aTrack)
std::string * theVolumeName
std::vector< float > l0fraction(std::string s)
tuple cout
Definition: gather_cfg.py:121
dbl *** dir
Definition: mlp_gen.cc:35