Geometry/src/MaterialBudgetAction.cc

0001
0002 #include "Validation/Geometry/interface/MaterialBudgetAction.h"
0003 #include "Validation/Geometry/interface/MaterialBudgetData.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0010 #include "SimG4Core/Notification/interface/EndOfRun.h"
0011 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0012 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0013 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0014
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017 #include "G4UImanager.hh"
0018 #include "G4UIterminal.hh"
0019 #include "G4UItcsh.hh"
0020 #include "G4LogicalVolumeStore.hh"
0021 #include "G4TouchableHistory.hh"
0022 #include "G4VPhysicalVolume.hh"
0023 #include "G4VProcess.hh"
0024 #include "G4ParticleTable.hh"
0025 #include "G4ParticleDefinition.hh"
0026 #include "G4ProcessManager.hh"
0027 #include "G4TransportationManager.hh"
0028 #include "DD4hep/Filter.h"
0029
0030 #include <iostream>
0031
0032 //-------------------------------------------------------------------------
0033 MaterialBudgetAction::MaterialBudgetAction(const edm::ParameterSet& iPSet) {
0034   theData = std::make_shared<MaterialBudgetData>();
0035
0036   edm::ParameterSet m_Anal = iPSet.getParameter<edm::ParameterSet>("MaterialBudgetAction");
0037
0038   //---- Accumulate material budget only inside selected volumes
0039   std::string theHistoList = m_Anal.getParameter<std::string>("HistogramList");
0040   std::vector<std::string> volList = m_Anal.getParameter<std::vector<std::string> >("SelectedVolumes");
0041
0042   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: List of the selected volumes:";
0043   for (const auto& it : volList) {
0044     if (it != "None") {
0045       theVolumeList.push_back(it);
0046       edm::LogVerbatim("MaterialBudget") << it;
0047     }
0048   }
0049
0050   // log
0051   if (theHistoList == "Tracker") {
0052     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: running in Tracker Mode";
0053   } else if (theHistoList == "ECAL") {
0054     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: running in Ecal Mode";
0055   } else if (theHistoList == "Mtd") {
0056     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: running in Mtd Mode";
0057   } else if (theHistoList == "HGCal") {
0058     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: running in HGCal Mode";
0059   } else {
0060     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: running in General Mode";
0061   }
0062
0063   //---- Stop track when a process occurs
0064   theProcessToStop = m_Anal.getParameter<std::string>("StopAfterProcess");
0065   LogDebug("MaterialBudget") << "MaterialBudgetAction: stop at process " << theProcessToStop;
0066
0067   //---- Save histos to ROOT file
0068   std::string saveToHistosFile = m_Anal.getParameter<std::string>("HistosFile");
0069   if (saveToHistosFile != "None") {
0070     saveToHistos = true;
0071     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: saving histograms to " << saveToHistosFile;
0072     theHistoMgr = std::make_shared<TestHistoMgr>();
0073     if (theHistoList == "Tracker") {
0074       theHistos = std::make_shared<MaterialBudgetTrackerHistos>(theData, theHistoMgr, saveToHistosFile);
0075     } else if (theHistoList == "ECAL") {
0076       theHistos = std::make_shared<MaterialBudgetEcalHistos>(theData, theHistoMgr, saveToHistosFile);
0077     } else if (theHistoList == "Mtd") {
0078       theHistos = std::make_shared<MaterialBudgetMtdHistos>(theData, theHistoMgr, saveToHistosFile);
0079       theData->setMtdmode(true);
0080     } else if (theHistoList == "HGCal") {
0081       theHistos = std::make_shared<MaterialBudgetHGCalHistos>(theData,
0082                                                               theHistoMgr,
0083                                                               saveToHistosFile,
0084                                                               m_Anal.getParameter<double>("minZ"),
0085                                                               m_Anal.getParameter<double>("maxZ"),
0086                                                               m_Anal.getParameter<int>("nintZ"),
0087                                                               m_Anal.getParameter<double>("rMin"),
0088                                                               m_Anal.getParameter<double>("rMax"),
0089                                                               m_Anal.getParameter<int>("nrbin"),
0090                                                               m_Anal.getParameter<double>("etaMin"),
0091                                                               m_Anal.getParameter<double>("etaMax"),
0092                                                               m_Anal.getParameter<int>("netabin"),
0093                                                               m_Anal.getParameter<double>("phiMin"),
0094                                                               m_Anal.getParameter<double>("phiMax"),
0095                                                               m_Anal.getParameter<int>("nphibin"),
0096                                                               m_Anal.getParameter<double>("RMin"),
0097                                                               m_Anal.getParameter<double>("RMax"),
0098                                                               m_Anal.getParameter<int>("nRbin"));
0099       //In HGCal mode, so tell data class
0100       theData->setHGCalmode(true);
0101     } else {
0102       theHistos = std::make_shared<MaterialBudgetHistos>(theData, theHistoMgr, saveToHistosFile);
0103     }
0104   } else {
0105     edm::LogWarning("MaterialBudget") << "MaterialBudgetAction: No histograms file specified";
0106     saveToHistos = false;
0107   }
0108
0109   //---- Save material budget info to TEXT file
0110   std::string saveToTxtFile = m_Anal.getParameter<std::string>("TextFile");
0111   if (saveToTxtFile != "None") {
0112     saveToTxt = true;
0113     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: saving text info to " << saveToTxtFile;
0114     theTxt = std::make_shared<MaterialBudgetTxt>(theData, saveToTxtFile);
0115   } else {
0116     saveToTxt = false;
0117   }
0118
0119   //---- Compute all the steps even if not stored on file
0120   bool allSteps = m_Anal.getParameter<bool>("AllStepsToTree");
0121   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: all steps are computed " << allSteps;
0122   if (allSteps)
0123     theData->SetAllStepsToTree();
0124
0125   //---- Save tree to ROOT file
0126   std::string saveToTreeFile = m_Anal.getParameter<std::string>("TreeFile");
0127   if (saveToTreeFile != "None") {
0128     saveToTree = true;
0129     theTree = std::make_shared<MaterialBudgetTree>(theData, saveToTreeFile);
0130   } else {
0131     saveToTree = false;
0132   }
0133   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: saving ROOT TTree to " << saveToTreeFile;
0134
0135   //---- Track the first decay products of the main particle
0136   // if their kinetic energy is greater than  Ekin
0137   storeDecay = m_Anal.getUntrackedParameter<bool>("storeDecay", false);
0138   Ekin = m_Anal.getUntrackedParameter<double>("EminDecayProd", 1000.0);  // MeV
0139   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction: decay products steps are stored (" << storeDecay
0140                                      << ") if their kinetic energy is greater than " << Ekin << " MeV";
0141   firstParticle = false;
0142 }
0143
0144 //-------------------------------------------------------------------------
0145 MaterialBudgetAction::~MaterialBudgetAction() {}
0146
0147 //-------------------------------------------------------------------------
0148 void MaterialBudgetAction::update(const BeginOfRun*) {
0149   //----- Check that selected volumes are indeed part of the geometry
0150   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0151
0152   for (const auto& volcite : theVolumeList) {
0153     bool volFound = false;
0154     for (const auto& lvcite : *lvs) {
0155       if (static_cast<std::string>(dd4hep::dd::noNamespace(lvcite->GetName())) == static_cast<std::string>(volcite)) {
0156         volFound = true;
0157         break;
0158       }
0159     }
0160     if (!volFound) {
0161       edm::LogWarning("MaterialBudget") << "MaterialBudgetAction: selected volume not found in geometry " << volcite;
0162     }
0163   }
0164
0165   //----- Check process selected is one of the available ones
0166   bool procFound = false;
0167   if (theProcessToStop == "None") {
0168     procFound = true;
0169   } else {
0170     G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
0171     int siz = partTable->size();
0172     for (int ii = 0; ii < siz; ii++) {
0173       G4ParticleDefinition* particle = partTable->GetParticle(ii);
0174
0175       //--- All processes of this particle
0176       G4ProcessManager* pmanager = particle->GetProcessManager();
0177       G4ProcessVector* pvect = pmanager->GetProcessList();
0178       int sizproc = pvect->size();
0179       for (int jj = 0; jj < sizproc; jj++) {
0180         if ((*pvect)[jj]->GetProcessName() == theProcessToStop) {
0181           procFound = true;
0182           break;
0183         }
0184       }
0185     }
0186   }
0187
0188   if (!procFound) {
0189     edm::LogWarning("MaterialBudget") << "MaterialBudgetAction: selected process to stop tracking not found "
0190                                       << theProcessToStop;
0191   }
0192 }
0193
0194 //-------------------------------------------------------------------------
0195 void MaterialBudgetAction::update(const BeginOfTrack* trk) {
0196   const G4Track* aTrack = (*trk)();  // recover G4 pointer if wanted
0197
0198   // that was a temporary action while we're sorting out
0199   // about # of secondaries (produced if CutsPerRegion=true)
0200
0201   LogDebug("MaterialBudget") << "MaterialBudgetAction: Track ID " << aTrack->GetTrackID() << "Track parent ID "
0202                              << aTrack->GetParentID() << "PDG Id. = " << aTrack->GetDefinition()->GetPDGEncoding()
0203                              << "Ekin = " << aTrack->GetKineticEnergy() << " MeV";
0204
0205   if (aTrack->GetCreatorProcess())
0206     LogDebug("MaterialBudget") << "MaterialBudgetAction: produced through "
0207                                << aTrack->GetCreatorProcess()->GetProcessType();
0208
0209   if (aTrack->GetTrackID() == 1) {
0210     firstParticle = true;
0211   } else {
0212     firstParticle = false;
0213   }
0214
0215   if (storeDecay) {  // if record of the decay is requested
0216     if (aTrack->GetCreatorProcess()) {
0217       if (aTrack->GetParentID() == 1 &&
0218           //      aTrack->GetCreatorProcess()->GetProcessType() == 6
0219           //      &&
0220           aTrack->GetKineticEnergy() > Ekin) {
0221         // continue
0222       } else {
0223         G4Track* aTracknc = const_cast<G4Track*>(aTrack);
0224         aTracknc->SetTrackStatus(fStopAndKill);
0225         return;
0226       }
0227     }  // particles produced from a decay (type=6) of the main particle (ID=1) with Kinetic Energy [MeV] > Ekin
0228   } else {  // kill all the other particles (take only the main one until it disappears) if decay not stored
0229     if (aTrack->GetParentID() != 0) {
0230       G4Track* aTracknc = const_cast<G4Track*>(aTrack);
0231       aTracknc->SetTrackStatus(fStopAndKill);
0232       return;
0233     }
0234   }
0235
0236   theData->dataStartTrack(aTrack);
0237
0238   if (saveToTree)
0239     theTree->fillStartTrack();
0240   if (saveToHistos)
0241     theHistos->fillStartTrack();
0242   if (saveToTxt)
0243     theTxt->fillStartTrack();
0244 }
0245
0246 //-------------------------------------------------------------------------
0247 void MaterialBudgetAction::update(const G4Step* aStep) {
0248   //----- Check it is inside one of the volumes selected
0249   if (!theVolumeList.empty()) {
0250     if (!CheckTouchableInSelectedVolumes(aStep->GetTrack()->GetTouchable()))
0251       return;
0252   }
0253
0254   //---------- each step
0255   theData->dataPerStep(aStep);
0256   if (saveToTree)
0257     theTree->fillPerStep();
0258   if (saveToHistos)
0259     theHistos->fillPerStep();
0260   if (saveToTxt)
0261     theTxt->fillPerStep();
0262
0263   //----- Stop tracking after selected process
0264   if (StopAfterProcess(aStep)) {
0265     G4Track* track = aStep->GetTrack();
0266     track->SetTrackStatus(fStopAndKill);
0267   }
0268
0269   return;
0270 }
0271
0272 //-------------------------------------------------------------------------
0273 std::string MaterialBudgetAction::getSubDetectorName(G4StepPoint* aStepPoint) {
0274   const G4TouchableHistory* theTouchable = (const G4TouchableHistory*)(aStepPoint->GetTouchable());
0275   G4int num_levels = theTouchable->GetHistoryDepth();
0276
0277   if (theTouchable->GetVolume()) {
0278     return static_cast<std::string>(dd4hep::dd::noNamespace(theTouchable->GetVolume(num_levels - 1)->GetName()));
0279   } else {
0280     return "OutOfWorld";
0281   }
0282 }
0283
0284 //-------------------------------------------------------------------------
0285 std::string MaterialBudgetAction::getPartName(G4StepPoint* aStepPoint) {
0286   const G4TouchableHistory* theTouchable = (const G4TouchableHistory*)(aStepPoint->GetTouchable());
0287   G4int num_levels = theTouchable->GetHistoryDepth();
0288
0289   if (theTouchable->GetVolume()) {
0290     return static_cast<std::string>(dd4hep::dd::noNamespace(theTouchable->GetVolume(num_levels - 3)->GetName()));
0291   } else {
0292     return "OutOfWorld";
0293   }
0294 }
0295
0296 //-------------------------------------------------------------------------
0297 void MaterialBudgetAction::update(const EndOfTrack* trk) {
0298   const G4Track* aTrack = (*trk)();  // recover G4 pointer if wanted
0299   theData->dataEndTrack(aTrack);
0300
0301   if (saveToTree)
0302     theTree->fillEndTrack();
0303   if (saveToHistos)
0304     theHistos->fillEndTrack();
0305   if (saveToTxt)
0306     theTxt->fillEndTrack();
0307 }
0308
0309 //-------------------------------------------------------------------------
0310 void MaterialBudgetAction::update(const EndOfRun*) {
0311   // endOfRun calls TestHistoMgr::save() allowing to write
0312   // the ROOT files containing the histograms
0313
0314   if (saveToHistos)
0315     theHistos->endOfRun();
0316   if (saveToTxt)
0317     theHistos->endOfRun();
0318   if (saveToTree)
0319     theTree->endOfRun();
0320
0321   return;
0322 }
0323
0324 //-------------------------------------------------------------------------
0325 bool MaterialBudgetAction::CheckTouchableInSelectedVolumes(const G4VTouchable* touch) {
0326   size_t volh = touch->GetHistoryDepth();
0327   for (int ii = volh; ii >= 0; ii--) {
0328     G4String name = static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(ii)->GetName()));
0329     if (std::find(theVolumeList.begin(), theVolumeList.end(), name) != theVolumeList.end())
0330       return true;
0331   }
0332   return false;
0333 }
0334
0335 //-------------------------------------------------------------------------
0336 bool MaterialBudgetAction::StopAfterProcess(const G4Step* aStep) {
0337   if (theProcessToStop.empty())
0338     return false;
0339
0340   if (aStep->GetPostStepPoint()->GetProcessDefinedStep() == nullptr)
0341     return false;
0342   if (aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == theProcessToStop) {
0343     edm::LogVerbatim("MaterialBudget") << "MaterialBudgetAction :"
0344                                        << aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
0345     return true;
0346   } else {
0347     return false;
0348   }
0349 }