Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:32:15

0001 #include "Validation/Geometry/interface/MaterialBudgetData.h"
0002 
0003 #include "G4Step.hh"
0004 #include "G4Material.hh"
0005 #include "G4EventManager.hh"
0006 #include "G4Event.hh"
0007 #include "DD4hep/Filter.h"
0008 
0009 MaterialBudgetData::MaterialBudgetData() {
0010   //instantiate categorizer to assign an ID to volumes and materials
0011   myMaterialBudgetCategorizer = nullptr;
0012   allStepsToTree = false;
0013   isHGCal = false;
0014   isMtd = false;
0015   densityConvertionFactor = 6.24E18;
0016 }
0017 
0018 MaterialBudgetData::~MaterialBudgetData() = default;
0019 
0020 void MaterialBudgetData::SetAllStepsToTree() { allStepsToTree = true; }
0021 
0022 void MaterialBudgetData::dataStartTrack(const G4Track* aTrack) {
0023   const G4ThreeVector& dir = aTrack->GetMomentum();
0024 
0025   if (myMaterialBudgetCategorizer == nullptr) {
0026     if (isHGCal) {
0027       myMaterialBudgetCategorizer = std::make_unique<MaterialBudgetCategorizer>("HGCal");
0028     } else if (isMtd) {
0029       myMaterialBudgetCategorizer = std::make_unique<MaterialBudgetCategorizer>("Mtd");
0030     } else {
0031       myMaterialBudgetCategorizer = std::make_unique<MaterialBudgetCategorizer>("Tracker");
0032     }
0033   }
0034 
0035   theStepN = 0;
0036   theTotalMB = 0;
0037   theTotalIL = 0;
0038   theEta = 0;
0039   thePhi = 0;
0040   theID = 0;
0041   thePt = 0;
0042   theEnergy = 0;
0043   theMass = 0;
0044 
0045   theSupportMB = 0.f;
0046   theSensitiveMB = 0.f;
0047   theCoolingMB = 0.f;
0048   theElectronicsMB = 0.f;
0049   theOtherMB = 0.f;
0050 
0051   //HGCal
0052   theAirMB = 0.f;
0053   theCablesMB = 0.f;
0054   theCopperMB = 0.f;
0055   theH_ScintillatorMB = 0.f;
0056   theLeadMB = 0.f;
0057   theEpoxyMB = 0.f;
0058   theKaptonMB = 0.f;
0059   theAluminiumMB = 0.f;
0060   theHGC_G10_FR4MB = 0.f;
0061   theSiliconMB = 0.f;
0062   theStainlessSteelMB = 0.f;
0063   theWCuMB = 0.f;
0064   thePolystyreneMB = 0.f;
0065   theHGC_EEConnectorMB = 0.f;
0066   theHGC_HEConnectorMB = 0.f;
0067 
0068   theSupportIL = 0.f;
0069   theSensitiveIL = 0.f;
0070   theCoolingIL = 0.f;
0071   theElectronicsIL = 0.f;
0072   theOtherIL = 0.f;
0073 
0074   //HGCal
0075   theAirIL = 0.f;
0076   theCablesIL = 0.f;
0077   theCopperIL = 0.f;
0078   theH_ScintillatorIL = 0.f;
0079   theLeadIL = 0.f;
0080   theEpoxyIL = 0.f;
0081   theKaptonIL = 0.f;
0082   theAluminiumIL = 0.f;
0083   theHGC_G10_FR4IL = 0.f;
0084   theSiliconIL = 0.f;
0085   theStainlessSteelIL = 0.f;
0086   theWCuIL = 0.f;
0087   thePolystyreneIL = 0.f;
0088   theHGC_EEConnectorIL = 0.f;
0089   theHGC_HEConnectorIL = 0.f;
0090 
0091   theSupportFractionMB = 0.f;
0092   theSensitiveFractionMB = 0.f;
0093   theCoolingFractionMB = 0.f;
0094   theElectronicsFractionMB = 0.f;
0095   theOtherFractionMB = 0.f;
0096   //HGCal
0097   theAirFractionMB = 0.f;
0098   theCablesFractionMB = 0.f;
0099   theCopperFractionMB = 0.f;
0100   theH_ScintillatorFractionMB = 0.f;
0101   theLeadFractionMB = 0.f;
0102   theEpoxyFractionMB = 0.f;
0103   theKaptonFractionMB = 0.f;
0104   theAluminiumFractionMB = 0.f;
0105   theHGC_G10_FR4FractionMB = 0.f;
0106   theSiliconFractionMB = 0.f;
0107   theStainlessSteelFractionMB = 0.f;
0108   theWCuFractionMB = 0.f;
0109   thePolystyreneFractionMB = 0.f;
0110   theHGC_EEConnectorFractionMB = 0.f;
0111   theHGC_HEConnectorFractionMB = 0.f;
0112 
0113   theSupportFractionIL = 0.f;
0114   theSensitiveFractionIL = 0.f;
0115   theCoolingFractionIL = 0.f;
0116   theElectronicsFractionIL = 0.f;
0117   theOtherFractionIL = 0.f;
0118   //HGCal
0119   theAirFractionIL = 0.f;
0120   theCablesFractionIL = 0.f;
0121   theCopperFractionIL = 0.f;
0122   theH_ScintillatorFractionIL = 0.f;
0123   theLeadFractionIL = 0.f;
0124   theEpoxyFractionIL = 0.f;
0125   theKaptonFractionIL = 0.f;
0126   theAluminiumFractionIL = 0.f;
0127   theHGC_G10_FR4FractionIL = 0.f;
0128   theSiliconFractionIL = 0.f;
0129   theStainlessSteelFractionIL = 0.f;
0130   theWCuFractionIL = 0.f;
0131   thePolystyreneFractionIL = 0.f;
0132   theHGC_EEConnectorFractionIL = 0.f;
0133   theHGC_HEConnectorFractionIL = 0.f;
0134 
0135   theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
0136   thePt = dir.perp();
0137   if (dir.theta() != 0) {
0138     theEta = dir.eta();
0139   } else {
0140     theEta = -99;
0141   }
0142   thePhi = dir.phi();
0143   theEnergy = aTrack->GetTotalEnergy();
0144   theMass = aTrack->GetDefinition()->GetPDGMass();
0145 }
0146 
0147 void MaterialBudgetData::dataEndTrack(const G4Track* aTrack) {
0148   LogDebug("MaterialBudget") << "MaterialBudgetData: [OVAL] MaterialBudget "
0149                              << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()
0150                              << " Eta:" << theEta << " Phi:" << thePhi << " TotalMB" << theTotalMB;
0151 
0152   LogDebug("MaterialBudget") << "MaterialBudgetData:" << theStepN << "Recorded steps ";
0153 
0154   if (!isHGCal) {
0155     LogDebug("Material Budget") << "MaterialBudgetData: Radiation Length "
0156                                 << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " Eta"
0157                                 << theEta << " Phi" << thePhi << " TotalMB" << theTotalMB << " SUP " << theSupportMB
0158                                 << " SEN " << theSensitiveMB << " CAB " << theCablesMB << " COL " << theCoolingMB
0159                                 << " ELE " << theElectronicsMB << " other " << theOtherMB << " Air " << theAirMB;
0160 
0161     LogDebug("Material Budget") << "MaterialBudgetData: Interaction Length "
0162                                 << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " Eta "
0163                                 << theEta << " Phi " << thePhi << " TotalIL " << theTotalIL << " SUP " << theSupportIL
0164                                 << " SEN " << theSensitiveIL << " CAB " << theCablesIL << " COL " << theCoolingIL
0165                                 << " ELE " << theElectronicsIL << " other " << theOtherIL << " Air " << theAirIL
0166                                 << std::endl;
0167 
0168   } else {
0169     LogDebug("MaterialBudget") << "MaterialBudgetData: HGCal Material Budget: Radiation Length "
0170                                << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " Eta "
0171                                << theEta << " Phi " << thePhi << " TotaLMB" << theTotalMB << " theCopperMB "
0172                                << theCopperMB << " theH_ScintillatorMB " << theH_ScintillatorMB << " CAB "
0173                                << theCablesMB << " theLeadMB " << theLeadMB << " theEpoxyMB " << theEpoxyMB
0174                                << " theKaptonMB " << theKaptonMB << " theAluminiumMB " << theAluminiumMB
0175                                << " theHGC_G10_FR4MB " << theHGC_G10_FR4MB << " theSiliconMB " << theSiliconMB
0176                                << " theAirMB " << theAirMB << " theStainlessSteelMB " << theStainlessSteelMB
0177                                << " theWCuMB " << theWCuMB << " thePolystyreneMB " << thePolystyreneMB
0178                                << " theHGC_EEConnectorMB " << theHGC_EEConnectorMB << " theHGC_HEConnectorMB "
0179                                << theHGC_HEConnectorMB;
0180 
0181     LogDebug("MaterialBudget") << "MaterialBudgetData: HGCal Material Budget: Interaction Length "
0182                                << G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID() << " Eta "
0183                                << theEta << " Phi " << thePhi << " TotalIL " << theTotalIL << " theCopperIL "
0184                                << theCopperIL << " theH_ScintillatorIL " << theH_ScintillatorIL << " CAB "
0185                                << theCablesIL << " theLeadIL " << theLeadIL << " theEpoxyIL " << theEpoxyIL
0186                                << " theKaptonIL " << theKaptonIL << " theAluminiumIL " << theAluminiumIL
0187                                << " theHGC_G10_FR4IL " << theHGC_G10_FR4IL << " theSiliconIL " << theSiliconIL
0188                                << " Air " << theAirIL << " theStainlessSteelIL " << theStainlessSteelIL << " theWCuIL "
0189                                << theWCuIL << " thePolystyreneIL " << thePolystyreneIL << " theHGC_EEConnectorIL "
0190                                << theHGC_EEConnectorIL << " theHGC_HEConnectorIL " << theHGC_HEConnectorIL << std::endl;
0191   }
0192 }
0193 
0194 void MaterialBudgetData::dataPerStep(const G4Step* aStep) {
0195   assert(aStep);
0196   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0197   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0198   assert(prePoint);
0199   assert(postPoint);
0200   G4Material* theMaterialPre = prePoint->GetMaterial();
0201   assert(theMaterialPre);
0202 
0203   CLHEP::Hep3Vector prePos = prePoint->GetPosition();
0204   CLHEP::Hep3Vector postPos = postPoint->GetPosition();
0205 
0206   G4double steplen = aStep->GetStepLength();
0207 
0208   G4double radlen = theMaterialPre->GetRadlen();
0209   G4double intlen = theMaterialPre->GetNuclearInterLength();
0210   G4double density = theMaterialPre->GetDensity() / densityConvertionFactor;  // always g/cm3
0211 
0212   G4String materialName = static_cast<std::string>(dd4hep::dd::noNamespace(theMaterialPre->GetName()));
0213 
0214   LogDebug("MaterialBudget") << "MaterialBudgetData: Material " << materialName << " steplen " << steplen << " radlen "
0215                              << radlen << " mb " << steplen / radlen;
0216 
0217   G4String volumeName = static_cast<std::string>(
0218       dd4hep::dd::noNamespace(aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0)->GetLogicalVolume()->GetName()));
0219 
0220   LogDebug("MaterialBudget") << "MaterialBudgetData: Volume " << volumeName << " Material " << materialName;
0221 
0222   // instantiate the categorizer
0223   assert(myMaterialBudgetCategorizer);
0224   int volumeID = myMaterialBudgetCategorizer->volume(volumeName);
0225   int materialID = myMaterialBudgetCategorizer->material(materialName);
0226 
0227   LogDebug("MaterialBudget") << "MaterialBudgetData: Volume ID " << volumeID << " Material ID " << materialID;
0228 
0229   // FIXME: Both volume ID and material ID are zeros, so this part is not executed leaving all
0230   // values as zeros.
0231 
0232   if (!isHGCal) {
0233     bool isCtgOk = !myMaterialBudgetCategorizer->x0fraction(materialName).empty() &&
0234                    !myMaterialBudgetCategorizer->l0fraction(materialName).empty() &&
0235                    (myMaterialBudgetCategorizer->x0fraction(materialName).size() == 7) /*7 Categories*/
0236                    && (myMaterialBudgetCategorizer->l0fraction(materialName).size() == 7);
0237 
0238     if (!isCtgOk) {
0239       if (materialName.compare("Air") == 0) {
0240         theAirFractionMB = 1.f;
0241         theAirFractionIL = 1.f;
0242       } else {
0243         theOtherFractionMB = 1.f;
0244         theOtherFractionIL = 1.f;
0245         edm::LogVerbatim("MaterialBudget") << "MaterialBudgetData: Material forced to 'Other': " << materialName
0246                                            << " in volume " << volumeName << ". Check Categorization.";
0247       }
0248     } else {
0249       theSupportFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[0];
0250       theSensitiveFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[1];
0251       theCablesFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[2];
0252       theCoolingFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[3];
0253       theElectronicsFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[4];
0254       theOtherFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[5];
0255       theAirFractionMB = myMaterialBudgetCategorizer->x0fraction(materialName)[6];
0256 
0257       if (theOtherFractionMB > 0.f)
0258         edm::LogVerbatim("MaterialBudget")
0259             << "MaterialBudgetData: Material found with no category: " << materialName << " in volume " << volumeName;
0260 
0261       theSupportFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[0];
0262       theSensitiveFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[1];
0263       theCablesFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[2];
0264       theCoolingFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[3];
0265       theElectronicsFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[4];
0266       theOtherFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[5];
0267       theAirFractionIL = myMaterialBudgetCategorizer->l0fraction(materialName)[6];
0268 
0269       if (theOtherFractionIL > 0.f)
0270         edm::LogVerbatim("MaterialBudget")
0271             << "MaterialBudgetData: Material found with no category: " << materialName << " in volume " << volumeName;
0272     }
0273   } else {  // isHGCal
0274 
0275     bool isHGCalx0fractionEmpty = myMaterialBudgetCategorizer->HGCalx0fraction(materialName).empty();
0276     bool isHGCall0fractionEmpty = myMaterialBudgetCategorizer->HGCall0fraction(materialName).empty();
0277 
0278     if (isHGCalx0fractionEmpty && isHGCall0fractionEmpty) {
0279       theOtherFractionMB = 1.f;
0280       theOtherFractionIL = 1.f;
0281 
0282       edm::LogVerbatim("MaterialBudget") << "MaterialBudgetData: Material forced to 'Other': " << materialName
0283                                          << " in volume " << volumeName;
0284     } else {
0285       theAirFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[0];
0286       theCablesFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[1];
0287       theCopperFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[2];
0288       theH_ScintillatorFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[3];
0289       theLeadFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[4];
0290       theHGC_G10_FR4FractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[5];
0291       theSiliconFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[6];
0292       theStainlessSteelFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[7];
0293       theWCuFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[8];
0294       theOtherFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[9];
0295       theEpoxyFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[10];
0296       theKaptonFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[11];
0297       theAluminiumFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[12];
0298       thePolystyreneFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[13];
0299       theHGC_EEConnectorFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[14];
0300       theHGC_HEConnectorFractionMB = myMaterialBudgetCategorizer->HGCalx0fraction(materialName)[15];
0301 
0302       if (theOtherFractionMB != 0)
0303         // edm::LogVerbatim("MaterialBudget") << "MaterialBudgetData: Material found with no category: " << materialName
0304         //                << " in volume " << volumeName << std::endl;
0305         std::cout << "MaterialBudgetData: Material found with no category: " << materialName << " in volume "
0306                   << volumeName << std::endl;
0307 
0308       theAirFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[0];
0309       theCablesFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[1];
0310       theCopperFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[2];
0311       theH_ScintillatorFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[3];
0312       theLeadFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[4];
0313       theHGC_G10_FR4FractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[5];
0314       theSiliconFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[6];
0315       theStainlessSteelFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[7];
0316       theWCuFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[8];
0317       theOtherFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[9];
0318       theEpoxyFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[10];
0319       theKaptonFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[11];
0320       theAluminiumFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[12];
0321       thePolystyreneFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[13];
0322       theHGC_EEConnectorFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[14];
0323       theHGC_HEConnectorFractionIL = myMaterialBudgetCategorizer->HGCall0fraction(materialName)[15];
0324 
0325       if (theOtherFractionIL != 0)
0326         edm::LogVerbatim("MaterialBudget") << "MaterialBudgetData: Material found with no category " << materialName
0327                                            << " in volume " << volumeName << std::endl;
0328     }
0329   }
0330 
0331   float dmb = steplen / radlen;
0332   float dil = steplen / intlen;
0333 
0334   G4VPhysicalVolume* pv = aStep->GetPreStepPoint()->GetPhysicalVolume();
0335   const G4VTouchable* t = aStep->GetPreStepPoint()->GetTouchable();
0336   const G4ThreeVector& objectTranslation = t->GetTranslation();
0337   const G4RotationMatrix* objectRotation = t->GetRotation();
0338   const G4VProcess* interactionPre = prePoint->GetProcessDefinedStep();
0339   const G4VProcess* interactionPost = postPoint->GetProcessDefinedStep();
0340 
0341   G4Track* track = aStep->GetTrack();
0342   if (theStepN == 0)
0343     LogDebug("MaterialBudget") << "MaterialBudgetData: Simulated Particle " << theID << "\tMass " << theMass
0344                                << " MeV/c2"
0345                                << "\tPt = " << thePt << " MeV/c"
0346                                << "\tEta = " << theEta << "\tPhi = " << thePhi << "\tEnergy = " << theEnergy << " MeV";
0347 
0348   //fill data per step
0349   if (allStepsToTree) {
0350     assert(theStepN < MAXNUMBERSTEPS);
0351     if (theStepN > MAXNUMBERSTEPS)
0352       theStepN = MAXNUMBERSTEPS - 1;
0353     theDmb[theStepN] = dmb;
0354     theDil[theStepN] = dil;
0355     theSupportDmb[theStepN] = (dmb * theSupportFractionMB);
0356     theSensitiveDmb[theStepN] = (dmb * theSensitiveFractionMB);
0357     theCoolingDmb[theStepN] = (dmb * theCoolingFractionMB);
0358     theElectronicsDmb[theStepN] = (dmb * theElectronicsFractionMB);
0359     theOtherDmb[theStepN] = (dmb * theOtherFractionMB);
0360     //HGCal
0361     theAirDmb[theStepN] = (dmb * theAirFractionMB);
0362     theCablesDmb[theStepN] = (dmb * theCablesFractionMB);
0363     theCopperDmb[theStepN] = (dmb * theCopperFractionMB);
0364     theH_ScintillatorDmb[theStepN] = (dmb * theH_ScintillatorFractionMB);
0365     theLeadDmb[theStepN] = (dmb * theLeadFractionMB);
0366     theEpoxyDmb[theStepN] = (dmb * theEpoxyFractionMB);
0367     theKaptonDmb[theStepN] = (dmb * theKaptonFractionMB);
0368     theAluminiumDmb[theStepN] = (dmb * theAluminiumFractionMB);
0369     theHGC_G10_FR4Dmb[theStepN] = (dmb * theHGC_G10_FR4FractionMB);
0370     theSiliconDmb[theStepN] = (dmb * theSiliconFractionMB);
0371     theStainlessSteelDmb[theStepN] = (dmb * theStainlessSteelFractionMB);
0372     theWCuDmb[theStepN] = (dmb * theWCuFractionMB);
0373     thePolystyreneDmb[theStepN] = (dmb * thePolystyreneFractionMB);
0374     theHGC_EEConnectorDmb[theStepN] = (dmb * theHGC_EEConnectorFractionMB);
0375     theHGC_HEConnectorDmb[theStepN] = (dmb * theHGC_HEConnectorFractionMB);
0376 
0377     theSupportDil[theStepN] = (dil * theSupportFractionIL);
0378     theSensitiveDil[theStepN] = (dil * theSensitiveFractionIL);
0379     theCoolingDil[theStepN] = (dil * theCoolingFractionIL);
0380     theElectronicsDil[theStepN] = (dil * theElectronicsFractionIL);
0381     theOtherDil[theStepN] = (dil * theOtherFractionIL);
0382     //HGCal
0383     theAirDil[theStepN] = (dil * theAirFractionIL);
0384     theCablesDil[theStepN] = (dil * theCablesFractionIL);
0385     theCopperDil[theStepN] = (dil * theCopperFractionIL);
0386     theH_ScintillatorDil[theStepN] = (dil * theH_ScintillatorFractionIL);
0387     theLeadDil[theStepN] = (dil * theLeadFractionIL);
0388     theEpoxyDil[theStepN] = (dil * theEpoxyFractionIL);
0389     theKaptonDil[theStepN] = (dil * theKaptonFractionIL);
0390     theAluminiumDil[theStepN] = (dil * theAluminiumFractionIL);
0391     theHGC_G10_FR4Dil[theStepN] = (dil * theHGC_G10_FR4FractionIL);
0392     theSiliconDil[theStepN] = (dil * theSiliconFractionIL);
0393     theStainlessSteelDil[theStepN] = (dil * theStainlessSteelFractionIL);
0394     theWCuDil[theStepN] = (dil * theWCuFractionIL);
0395     thePolystyreneDil[theStepN] = (dil * thePolystyreneFractionIL);
0396     theHGC_EEConnectorDil[theStepN] = (dil * theHGC_EEConnectorFractionIL);
0397     theHGC_HEConnectorDil[theStepN] = (dil * theHGC_HEConnectorFractionIL);
0398 
0399     theInitialX[theStepN] = prePos.x();
0400     theInitialY[theStepN] = prePos.y();
0401     theInitialZ[theStepN] = prePos.z();
0402     theFinalX[theStepN] = postPos.x();
0403     theFinalY[theStepN] = postPos.y();
0404     theFinalZ[theStepN] = postPos.z();
0405     theVolumeID[theStepN] = volumeID;
0406     theVolumeName[theStepN] = volumeName;
0407     theVolumeCopy[theStepN] = pv->GetCopyNo();
0408     theVolumeX[theStepN] = objectTranslation.x();
0409     theVolumeY[theStepN] = objectTranslation.y();
0410     theVolumeZ[theStepN] = objectTranslation.z();
0411     theVolumeXaxis1[theStepN] = objectRotation->xx();
0412     theVolumeXaxis2[theStepN] = objectRotation->xy();
0413     theVolumeXaxis3[theStepN] = objectRotation->xz();
0414     theVolumeYaxis1[theStepN] = objectRotation->yx();
0415     theVolumeYaxis2[theStepN] = objectRotation->yy();
0416     theVolumeYaxis3[theStepN] = objectRotation->yz();
0417     theVolumeZaxis1[theStepN] = objectRotation->zx();
0418     theVolumeZaxis2[theStepN] = objectRotation->zy();
0419     theVolumeZaxis3[theStepN] = objectRotation->zz();
0420     theMaterialID[theStepN] = materialID;
0421     theMaterialName[theStepN] = materialName;
0422     theMaterialX0[theStepN] = radlen;
0423     theMaterialLambda0[theStepN] = intlen;
0424     theMaterialDensity[theStepN] = density;
0425     theStepID[theStepN] = track->GetDefinition()->GetPDGEncoding();
0426     theStepInitialPt[theStepN] = prePoint->GetMomentum().perp();
0427     theStepInitialEta[theStepN] = prePoint->GetMomentum().eta();
0428     theStepInitialPhi[theStepN] = prePoint->GetMomentum().phi();
0429     theStepInitialEnergy[theStepN] = prePoint->GetTotalEnergy();
0430     theStepInitialPx[theStepN] = prePoint->GetMomentum().x();
0431     theStepInitialPy[theStepN] = prePoint->GetMomentum().y();
0432     theStepInitialPz[theStepN] = prePoint->GetMomentum().z();
0433     theStepInitialBeta[theStepN] = prePoint->GetBeta();
0434     theStepInitialGamma[theStepN] = prePoint->GetGamma();
0435     theStepInitialMass[theStepN] = prePoint->GetMass();
0436     theStepFinalPt[theStepN] = postPoint->GetMomentum().perp();
0437     theStepFinalEta[theStepN] = postPoint->GetMomentum().eta();
0438     theStepFinalPhi[theStepN] = postPoint->GetMomentum().phi();
0439     theStepFinalEnergy[theStepN] = postPoint->GetTotalEnergy();
0440     theStepFinalPx[theStepN] = postPoint->GetMomentum().x();
0441     theStepFinalPy[theStepN] = postPoint->GetMomentum().y();
0442     theStepFinalPz[theStepN] = postPoint->GetMomentum().z();
0443     theStepFinalBeta[theStepN] = postPoint->GetBeta();
0444     theStepFinalGamma[theStepN] = postPoint->GetGamma();
0445     theStepFinalMass[theStepN] = postPoint->GetMass();
0446     int preProcType = -99;
0447     int postProcType = -99;
0448     if (interactionPre)
0449       preProcType = interactionPre->GetProcessType();
0450     theStepPreProcess[theStepN] = preProcType;
0451     if (interactionPost)
0452       postProcType = interactionPost->GetProcessType();
0453     theStepPostProcess[theStepN] = postProcType;
0454 
0455     LogDebug("MaterialBudget") << "MaterialBudgetData: Step " << theStepN << "\tDelta MB = " << theDmb[theStepN]
0456                                << std::endl
0457                                << " Support " << theSupportDmb[theStepN] << " Sensitive " << theSensitiveDmb[theStepN]
0458                                << " Cables " << theCablesDmb[theStepN] << " Cooling " << theCoolingDmb[theStepN]
0459                                << " Electronics " << theElectronicsDmb[theStepN] << " Other " << theOtherDmb[theStepN]
0460                                << " Air " << theAirDmb[theStepN] << std::endl
0461                                << "\tDelta IL = " << theDil[theStepN] << std::endl
0462                                << " Support " << theSupportDil[theStepN] << " Sensitive " << theSensitiveDil[theStepN]
0463                                << " Cables " << theCablesDil[theStepN] << " Cooling " << theCoolingDil[theStepN]
0464                                << " Electronics " << theElectronicsDil[theStepN] << " Other " << theOtherDil[theStepN]
0465                                << " Air " << theAirDil[theStepN];
0466 
0467     if (interactionPre)
0468       LogDebug("MaterialBudget") << "MaterialBudgetData: Process Pre " << interactionPre->GetProcessName() << " type "
0469                                  << theStepPreProcess[theStepN] << " name "
0470                                  << interactionPre->GetProcessTypeName(G4ProcessType(theStepPreProcess[theStepN]));
0471     if (interactionPost)
0472       LogDebug("MaterialBudget")
0473           << "MaterialBudgetData: Process Post " << interactionPost->GetProcessName() << " type "
0474           << theStepPostProcess[theStepN] << " name "
0475           << interactionPost->GetProcessTypeName(G4ProcessType(theStepPostProcess[theStepN]))
0476           << " Pre x = " << theInitialX[theStepN] << " y = " << theInitialY[theStepN]
0477           << " z = " << theInitialZ[theStepN]
0478           << " Polar Radius = " << sqrt(prePos.x() * prePos.x() + prePos.y() * prePos.y())
0479           << " Pt = " << theStepInitialPt[theStepN] << " Energy = " << theStepInitialEnergy[theStepN] << " Final: "
0480           << " Post x = " << theFinalX[theStepN] << " y = " << theFinalY[theStepN] << " z = " << theFinalZ[theStepN]
0481           << " Polar Radius = " << sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y())
0482           << " Pt = " << theStepFinalPt[theStepN] << " Energy = " << theStepFinalEnergy[theStepN] << std::endl
0483           << " Volume " << volumeID << " name " << theVolumeName[theStepN] << " copy number " << theVolumeCopy[theStepN]
0484           << " material " << theMaterialID[theStepN] << " " << theMaterialName[theStepN]
0485           << " Density = " << theMaterialDensity[theStepN] << " g/cm3"
0486           << " X0 = " << theMaterialX0[theStepN] << " mm"
0487           << " Lambda0 = " << theMaterialLambda0[theStepN] << " mm" << std::endl
0488           << " Particle " << theStepID[theStepN] << " Initial Pt = " << theStepInitialPt[theStepN] << " MeV/c"
0489           << " eta = " << theStepInitialEta[theStepN] << " phi = " << theStepInitialPhi[theStepN]
0490           << " Energy = " << theStepInitialEnergy[theStepN] << " MeV"
0491           << " Mass = " << theStepInitialMass[theStepN] << " MeV/c2"
0492           << " Beta = " << theStepInitialBeta[theStepN] << " Gamma = " << theStepInitialGamma[theStepN] << std::endl
0493           << " Particle " << theStepID[theStepN] << " Final Pt = " << theStepFinalPt[theStepN] << " MeV/c"
0494           << " eta = " << theStepFinalEta[theStepN] << " phi = " << theStepFinalPhi[theStepN]
0495           << " Energy = " << theStepFinalEnergy[theStepN] << " MeV"
0496           << " Mass = " << theStepFinalMass[theStepN] << " MeV/c2"
0497           << " Beta = " << theStepFinalBeta[theStepN] << " Gamma = " << theStepFinalGamma[theStepN] << std::endl
0498           << " Volume Centre x = " << theVolumeX[theStepN] << " y = " << theVolumeY[theStepN]
0499           << " z = " << theVolumeZ[theStepN] << "Polar Radius = "
0500           << sqrt(theVolumeX[theStepN] * theVolumeX[theStepN] + theVolumeY[theStepN] * theVolumeY[theStepN])
0501           << std::endl
0502           << " x axis = (" << theVolumeXaxis1[theStepN] << "," << theVolumeXaxis2[theStepN] << ","
0503           << theVolumeXaxis3[theStepN] << ")" << std::endl
0504           << " y axis = (" << theVolumeYaxis1[theStepN] << "," << theVolumeYaxis2[theStepN] << ","
0505           << theVolumeYaxis3[theStepN] << ")" << std::endl
0506           << " z axis = (" << theVolumeZaxis1[theStepN] << "," << theVolumeZaxis2[theStepN] << ","
0507           << theVolumeZaxis3[theStepN] << ")" << std::endl
0508           << " Secondaries" << std::endl;
0509 
0510     for (G4TrackVector::const_iterator iSec = aStep->GetSecondary()->begin(); iSec != aStep->GetSecondary()->end();
0511          iSec++) {
0512       LogDebug("MaterialBudget") << "MaterialBudgetData: tid " << (*iSec)->GetDefinition()->GetPDGEncoding()
0513                                  << " created through process type " << (*iSec)->GetCreatorProcess()->GetProcessType()
0514                                  << (*iSec)->GetCreatorProcess()->GetProcessTypeName(
0515                                         G4ProcessType((*iSec)->GetCreatorProcess()->GetProcessType()));
0516     }
0517   }
0518 
0519   theTrkLen = aStep->GetTrack()->GetTrackLength();
0520   thePVname = static_cast<std::string>(dd4hep::dd::noNamespace(pv->GetName()));
0521   thePVcopyNo = pv->GetCopyNo();
0522   theRadLen = radlen;
0523   theIntLen = intlen;
0524   theTotalMB += dmb;
0525   theTotalIL += dil;
0526 
0527   theSupportMB += (dmb * theSupportFractionMB);
0528   theSensitiveMB += (dmb * theSensitiveFractionMB);
0529   theCoolingMB += (dmb * theCoolingFractionMB);
0530   theElectronicsMB += (dmb * theElectronicsFractionMB);
0531   theOtherMB += (dmb * theOtherFractionMB);
0532 
0533   //HGCal
0534   theAirMB += (dmb * theAirFractionMB);
0535   theCablesMB += (dmb * theCablesFractionMB);
0536   theCopperMB += (dmb * theCopperFractionMB);
0537   theH_ScintillatorMB += (dmb * theH_ScintillatorFractionMB);
0538   theLeadMB += (dmb * theLeadFractionMB);
0539   theEpoxyMB += (dmb * theEpoxyFractionMB);
0540   theKaptonMB += (dmb * theKaptonFractionMB);
0541   theAluminiumMB += (dmb * theAluminiumFractionMB);
0542   theHGC_G10_FR4MB += (dmb * theHGC_G10_FR4FractionMB);
0543   theSiliconMB += (dmb * theSiliconFractionMB);
0544   theStainlessSteelMB += (dmb * theStainlessSteelFractionMB);
0545   theWCuMB += (dmb * theWCuFractionMB);
0546   thePolystyreneMB += (dmb * thePolystyreneFractionMB);
0547   theHGC_EEConnectorMB += (dmb * theHGC_EEConnectorFractionMB);
0548   theHGC_HEConnectorMB += (dmb * theHGC_HEConnectorFractionMB);
0549 
0550   theSupportIL += (dil * theSupportFractionIL);
0551   theSensitiveIL += (dil * theSensitiveFractionIL);
0552   theCoolingIL += (dil * theCoolingFractionIL);
0553   theElectronicsIL += (dil * theElectronicsFractionIL);
0554   theOtherIL += (dil * theOtherFractionIL);
0555   //HGCal
0556   theAirIL += (dil * theAirFractionIL);
0557   theCablesIL += (dil * theCablesFractionIL);
0558   theCopperIL += (dil * theCopperFractionIL);
0559   theH_ScintillatorIL += (dil * theH_ScintillatorFractionIL);
0560   theLeadIL += (dil * theLeadFractionIL);
0561   theEpoxyIL += (dil * theEpoxyFractionIL);
0562   theKaptonIL += (dil * theKaptonFractionIL);
0563   theAluminiumIL += (dil * theAluminiumFractionIL);
0564   theHGC_G10_FR4IL += (dil * theHGC_G10_FR4FractionIL);
0565   theSiliconIL += (dil * theSiliconFractionIL);
0566   theStainlessSteelIL += (dil * theStainlessSteelFractionIL);
0567   theWCuIL += (dil * theWCuFractionIL);
0568   thePolystyreneIL += (dil * thePolystyreneFractionIL);
0569   theHGC_EEConnectorIL += (dil * theHGC_EEConnectorFractionIL);
0570   theHGC_HEConnectorIL += (dil * theHGC_HEConnectorFractionIL);
0571 
0572   // rr
0573 
0574   theStepN++;
0575 }