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];
54  theVolumeName = new std::string[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];
69  theMaterialName = new std::string[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  // rr
103  std::cout << "MaterialBudget Analysis of Event #" << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << std::endl;
104  // rr
105 
106  const G4ThreeVector& dir = aTrack->GetMomentum() ;
107 
109 
110  theStepN=0;
111  theTotalMB=0;
112  theTotalIL=0;
113  theEta=0;
114  thePhi=0;
115 
116  // rr
117  theID=0;
118  thePt=0;
119  theEnergy=0;
120  theMass=0;
121 
122  theSupportMB = 0.;
123  theSensitiveMB = 0.;
124  theCablesMB = 0.;
125  theCoolingMB = 0.;
126  theElectronicsMB = 0.;
127  theOtherMB = 0.;
128  theAirMB = 0.;
129  theSupportIL = 0.;
130  theSensitiveIL = 0.;
131  theCablesIL = 0.;
132  theCoolingIL = 0.;
133  theElectronicsIL = 0.;
134  theOtherIL = 0.;
135  theAirIL = 0.;
138  theCablesFractionMB = 0.;
141  theOtherFractionMB = 0.;
142  theAirFractionMB = 0.;
145  theCablesFractionIL = 0.;
148  theOtherFractionIL = 0.;
149  theAirFractionIL = 0.;
150  // rr
151 
152  theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
153  thePt = dir.perp();
154  if( dir.theta() != 0 ) {
155  theEta = dir.eta();
156  } else {
157  theEta = -99;
158  }
159  // thePhi = dir.phi()/deg; // better not to store in deg
160  thePhi = dir.phi();
161  theEnergy = aTrack->GetTotalEnergy();
162  theMass = aTrack->GetDefinition()->GetPDGMass();
163 
164 }
165 
166 
167 void MaterialBudgetData::dataEndTrack( const G4Track* aTrack )
168 {
169  //- std::cout << "[OVAL] MaterialBudget " << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " " << theEta << " " << thePhi << " " << theTotalMB << std::endl;
170  // rr
171  std::cout << "Recorded steps " << theStepN << std::endl;
172  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;
173  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;
174  // rr
175 }
176 
177 
178 void MaterialBudgetData::dataPerStep( const G4Step* aStep )
179 {
180  G4Material * theMaterialPre = aStep->GetPreStepPoint()->GetMaterial();
181  // G4Material * theMaterialPost = aStep->GetPostStepPoint()->GetMaterial();
182 
183  G4StepPoint* prePoint = aStep->GetPreStepPoint();
184  G4StepPoint* postPoint = aStep->GetPostStepPoint();
185 
186  CLHEP::Hep3Vector prePos = prePoint->GetPosition();
187  CLHEP::Hep3Vector postPos = postPoint->GetPosition();
188 
189  G4double steplen = aStep->GetStepLength();
190 
191  G4double radlen;
192  G4double intlen;
193  G4double density;
194 
195  radlen = theMaterialPre->GetRadlen();
196  intlen = theMaterialPre->GetNuclearInterLength();
197  density = theMaterialPre->GetDensity() / densityConvertionFactor; // always g/cm3
198 
199  G4String name = theMaterialPre->GetName();
200  // std::cout << " steplen " << steplen << " radlen " << radlen << " mb " << steplen/radlen << " mate " << theMaterialPre->GetName() << std::endl;
201 
202  G4LogicalVolume* lv = aStep->GetTrack()->GetVolume()->GetLogicalVolume();
203 
204  // instantiate the categorizer
205  int volumeID = myMaterialBudgetCategorizer->volume( lv->GetName() );
206  int materialID = myMaterialBudgetCategorizer->material( lv->GetMaterial()->GetName() );
207  // rr
208  std::string volumeName = lv->GetName();
209  std::string materialName = lv->GetMaterial()->GetName();
210  // rr
211 
212  // rr
213  /*
214  std::cout << " Volume " << lv->GetName() << std::endl;
215  std::cout << " Material " << lv->GetMaterial()->GetName() << std::endl;
216  */
217  theSupportFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[0];
218  theSensitiveFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[1];
219  theCablesFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[2];
220  theCoolingFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[3];
221  theElectronicsFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[4];
222  theOtherFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[5];
223  theAirFractionMB = myMaterialBudgetCategorizer->x0fraction(lv->GetMaterial()->GetName())[6];
224  if(theOtherFractionMB!=0) std::cout << " material found with no category " << lv->GetMaterial()->GetName()
225  << " in volume " << lv->GetName() << std::endl;
226  theSupportFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[0];
227  theSensitiveFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[1];
228  theCablesFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[2];
229  theCoolingFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[3];
230  theElectronicsFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[4];
231  theOtherFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[5];
232  theAirFractionIL = myMaterialBudgetCategorizer->l0fraction(lv->GetMaterial()->GetName())[6];
233  if(theOtherFractionIL!=0) std::cout << " material found with no category " << lv->GetMaterial()->GetName()
234  << " in volume " << lv->GetName() << std::endl;
235  // if(theOtherFractionMB!=0) LogDebug("MaterialBudgetData") << " material found with no category " << lv->GetMaterial()->GetName()
236  // << " in volume " << lv->GetName();
237  // rr
238 
239  float dmb = steplen/radlen;
240  float dil = steplen/intlen;
241 
242  G4VPhysicalVolume* pv = aStep->GetPreStepPoint()->GetPhysicalVolume();
243  const G4VTouchable* t = aStep->GetPreStepPoint()->GetTouchable();
244  G4ThreeVector objectTranslation = t->GetTranslation();
245  const G4RotationMatrix* objectRotation = t->GetRotation();
246  const G4VProcess* interactionPre = prePoint->GetProcessDefinedStep();
247  const G4VProcess* interactionPost = postPoint->GetProcessDefinedStep();
248 
249  G4Track* track = aStep->GetTrack();
250  if(theStepN==0) std::cout << " Simulated Particle " << theID << "\tMass " << theMass << " MeV/c2"
251  << "\tPt = " << thePt << " MeV/c" << "\tEta = " << theEta << "\tPhi = " << thePhi
252  << "\tEnergy = " << theEnergy << " MeV"
253  // << std::endl
254  // << "\tMagnetic Field at (0,0,0): (" << B000[0] << "," < B000[1] << "," << B000[2] << ")"
255  << std::endl;
256 
257  //fill data per step
258  if( allStepsToTree ){
260  theDmb[theStepN] = dmb;
261  theDil[theStepN] = dil;
276  theInitialX[theStepN] = prePos.x();
277  theInitialY[theStepN] = prePos.y();
278  theInitialZ[theStepN] = prePos.z();
279  theFinalX[theStepN] = postPos.x();
280  theFinalY[theStepN] = postPos.y();
281  theFinalZ[theStepN] = postPos.z();
282  theVolumeID[theStepN] = volumeID;
283  theVolumeName[theStepN] = volumeName;
284  theVolumeCopy[theStepN] = pv->GetCopyNo();
285  theVolumeX[theStepN] = objectTranslation.x();
286  theVolumeY[theStepN] = objectTranslation.y();
287  theVolumeZ[theStepN] = objectTranslation.z();
288  theVolumeXaxis1[theStepN] = objectRotation->xx();
289  theVolumeXaxis2[theStepN] = objectRotation->xy();
290  theVolumeXaxis3[theStepN] = objectRotation->xz();
291  theVolumeYaxis1[theStepN] = objectRotation->yx();
292  theVolumeYaxis2[theStepN] = objectRotation->yy();
293  theVolumeYaxis3[theStepN] = objectRotation->yz();
294  theVolumeZaxis1[theStepN] = objectRotation->zx();
295  theVolumeZaxis2[theStepN] = objectRotation->zy();
296  theVolumeZaxis3[theStepN] = objectRotation->zz();
297  theMaterialID[theStepN] = materialID;
298  theMaterialName[theStepN] = materialName;
299  theMaterialX0[theStepN] = radlen;
300  theMaterialLambda0[theStepN] = intlen;
301  theMaterialDensity[theStepN] = density;
302  theStepID[theStepN] = track->GetDefinition()->GetPDGEncoding();
303  theStepInitialPt[theStepN] = prePoint->GetMomentum().perp();
304  theStepInitialEta[theStepN] = prePoint->GetMomentum().eta();
305  theStepInitialPhi[theStepN] = prePoint->GetMomentum().phi();
306  theStepInitialEnergy[theStepN] = prePoint->GetTotalEnergy();
307  theStepInitialPx[theStepN] = prePoint->GetMomentum().x();
308  theStepInitialPy[theStepN] = prePoint->GetMomentum().y();
309  theStepInitialPz[theStepN] = prePoint->GetMomentum().z();
310  theStepInitialBeta[theStepN] = prePoint->GetBeta();
311  theStepInitialGamma[theStepN] = prePoint->GetGamma();
312  theStepInitialMass[theStepN] = prePoint->GetMass();
313  theStepFinalPt[theStepN] = postPoint->GetMomentum().perp();
314  theStepFinalEta[theStepN] = postPoint->GetMomentum().eta();
315  theStepFinalPhi[theStepN] = postPoint->GetMomentum().phi();
316  theStepFinalEnergy[theStepN] = postPoint->GetTotalEnergy();
317  theStepFinalPx[theStepN] = postPoint->GetMomentum().x();
318  theStepFinalPy[theStepN] = postPoint->GetMomentum().y();
319  theStepFinalPz[theStepN] = postPoint->GetMomentum().z();
320  theStepFinalBeta[theStepN] = postPoint->GetBeta();
321  theStepFinalGamma[theStepN] = postPoint->GetGamma();
322  theStepFinalMass[theStepN] = postPoint->GetMass();
323  int preProcType = -99;
324  int postProcType = -99;
325  if (interactionPre) preProcType = interactionPre->GetProcessType();
326  theStepPreProcess[theStepN] = preProcType;
327  if (interactionPost) postProcType = interactionPost->GetProcessType();
328  theStepPostProcess[theStepN] = postProcType;
329 #ifdef TREE_DEBUG
330  std::cout << " step " << theStepN
331  << "\tDelta MB = " << theDmb[theStepN]
332  << std::endl
333  << "\t\tSupport " << theSupportDmb[theStepN]
334  << " Sensitive " << theSensitiveDmb[theStepN]
335  << " Cables " << theCablesDmb[theStepN]
336  << " Cooling " << theCoolingDmb[theStepN]
337  << " Electronics " << theElectronicsDmb[theStepN]
338  << " Other " << theOtherDmb[theStepN]
339  << " Air " << theAirDmb[theStepN]
340  << std::endl
341  << "\tDelta IL = " << theDil[theStepN]
342  << std::endl
343  << "\t\tSupport " << theSupportDil[theStepN]
344  << " Sensitive " << theSensitiveDil[theStepN]
345  << " Cables " << theCablesDil[theStepN]
346  << " Cooling " << theCoolingDil[theStepN]
347  << " Electronics " << theElectronicsDil[theStepN]
348  << " Other " << theOtherDil[theStepN]
349  << " Air " << theAirDil[theStepN]
350  << std::endl;
351  if (interactionPre)
352  std::cout << "\tProcess Pre " << interactionPre->GetProcessName()
353  << " type " << theStepPreProcess[theStepN] << " " << interactionPre->GetProcessTypeName(G4ProcessType(theStepPreProcess[theStepN]))
354  << std::endl;
355  if (interactionPost)
356  std::cout << "\tProcess Post " << interactionPost->GetProcessName()
357  << " type " << theStepPostProcess[theStepN] << " "
358  << interactionPost->GetProcessTypeName(G4ProcessType(theStepPostProcess[theStepN]))
359  << std::endl;
360  std::cout << "\tPre x = " << theInitialX[theStepN]
361  << "\ty = " << theInitialY[theStepN]
362  << "\tz = " << theInitialZ[theStepN]
363  << "\tPolar Radius = " << sqrt(prePos.x()*prePos.x()+prePos.y()*prePos.y())
364  << "\tPt = " << theStepInitialPt[theStepN]
365  << "\tEnergy = " << theStepInitialEnergy[theStepN]
366  // << std::endl
367  // << "B-field(T) at Pre (cm): " << field->fieldInTesla(GlobalPoint(pos.x()/10.,pos.y()/10.,pos.z()/10.))
368  << std::endl;
369  std::cout << "\tPost x = " << theFinalX[theStepN]
370  << "\ty = " << theFinalY[theStepN]
371  << "\tz = " << theFinalZ[theStepN]
372  << "\tPolar Radius = " << sqrt(postPos.x()*postPos.x()+postPos.y()*postPos.y())
373  << "\tPt = " << theStepFinalPt[theStepN]
374  << "\tEnergy = " << theStepFinalEnergy[theStepN]
375  << std::endl;
376  std::cout << "\tvolume " << volumeID << " " << theVolumeName[theStepN]
377  << " copy number " << theVolumeCopy[theStepN]
378  << "\tmaterial " << theMaterialID[theStepN] << " " << theMaterialName[theStepN]
379  << "\tDensity = " << theMaterialDensity[theStepN] << " g/cm3"
380  << "\tX0 = " << theMaterialX0[theStepN] << " mm"
381  << "\tLambda0 = " << theMaterialLambda0[theStepN] << " mm"
382  << std::endl;
383  std::cout << "\t\tParticle " << theStepID[theStepN]
384  << " Initial Pt = " << theStepInitialPt[theStepN] << " MeV/c"
385  << " eta = " << theStepInitialEta[theStepN]
386  << " phi = " << theStepInitialPhi[theStepN]
387  << " Energy = " << theStepInitialEnergy[theStepN] << " MeV"
388  << " Mass = " << theStepInitialMass[theStepN] << " MeV/c2"
389  << " Beta = " << theStepInitialBeta[theStepN]
390  << " Gamma = " << theStepInitialGamma[theStepN]
391  << std::endl
392  << "\t\tParticle " << theStepID[theStepN]
393  << " Final Pt = " << theStepFinalPt[theStepN] << " MeV/c"
394  << " eta = " << theStepFinalEta[theStepN]
395  << " phi = " << theStepFinalPhi[theStepN]
396  << " Energy = " << theStepFinalEnergy[theStepN] << " MeV"
397  << " Mass = " << theStepFinalMass[theStepN] << " MeV/c2"
398  << " Beta = " << theStepFinalBeta[theStepN]
399  << " Gamma = " << theStepFinalGamma[theStepN]
400  << std::endl;
401  std:: cout << "\tVolume Centre x = " << theVolumeX[theStepN]
402  << "\ty = " << theVolumeY[theStepN]
403  << "\tz = " << theVolumeZ[theStepN]
404  << "\tPolar Radius = " << sqrt( theVolumeX[theStepN]*theVolumeX[theStepN] +
405  theVolumeY[theStepN]*theVolumeY[theStepN] )
406  << std::endl;
407  std::cout << "\tx axis = ("
408  << theVolumeXaxis1[theStepN] << ","
409  << theVolumeXaxis2[theStepN] << ","
410  << theVolumeXaxis3[theStepN] << ")"
411  << std::endl;
412  std::cout << "\ty axis = ("
413  << theVolumeYaxis1[theStepN] << ","
414  << theVolumeYaxis2[theStepN] << ","
415  << theVolumeYaxis3[theStepN] << ")"
416  << std::endl;
417  std::cout << "\tz axis = ("
418  << theVolumeZaxis1[theStepN] << ","
419  << theVolumeZaxis2[theStepN] << ","
420  << theVolumeZaxis3[theStepN] << ")"
421  << std::endl;
422  std::cout << "\tSecondaries"
423  << std::endl;
424  for(G4TrackVector::iterator iSec = aStep->GetSecondary()->begin(); iSec!=aStep->GetSecondary()->end(); iSec++) {
425  std::cout << "\t\tid " << (*iSec)->GetDefinition()->GetPDGEncoding()
426  << " created through process "
427  << " type " << (*iSec)->GetCreatorProcess()->GetProcessType()
428  << " " << (*iSec)->GetCreatorProcess()->GetProcessTypeName(G4ProcessType((*iSec)->GetCreatorProcess()->GetProcessType()))
429  << std::endl;
430  }
431 #endif
432  }
433 
434  theTrkLen = aStep->GetTrack()->GetTrackLength();
435  //- std::cout << " theTrkLen " << theTrkLen << " theTrkLen2 " << theTrkLen2 << " postPos " << postPos.mag() << postPos << std::endl;
436  thePVname = pv->GetName();
437  thePVcopyNo = pv->GetCopyNo();
438  theRadLen = radlen;
439  theIntLen = intlen;
440  theTotalMB += dmb;
441  theTotalIL += dil;
442 
443  // rr
446  theCablesMB += (dmb * theCablesFractionMB);
449  theOtherMB += (dmb * theOtherFractionMB);
450  theAirMB += (dmb * theAirFractionMB);
453  theCablesIL += (dil * theCablesFractionIL);
456  theOtherIL += (dil * theOtherFractionIL);
457  theAirIL += (dil * theAirFractionIL);
458  // rr
459 
460  theStepN++;
461 
462 }
463 
464 
std::vector< float > x0fraction(std::string s)
std::string * theMaterialName
T sqrt(T t)
Definition: SSEVec.h:46
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