Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-21 23:59:15

0001 #include "Validation/Geometry/interface/MaterialBudget.h"
0002 
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0007 
0008 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0009 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0010 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0011 
0012 #include "G4LogicalVolumeStore.hh"
0013 #include "G4Step.hh"
0014 #include "G4Track.hh"
0015 #include "DD4hep/Filter.h"
0016 
0017 #include <iostream>
0018 
0019 MaterialBudget::MaterialBudget(const edm::ParameterSet& p) {
0020   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudget");
0021   detTypes = m_p.getParameter<std::vector<std::string> >("DetectorTypes");
0022   constituents = m_p.getParameter<std::vector<int> >("Constituents");
0023   stackOrder = m_p.getParameter<std::vector<int> >("StackOrder");
0024   detNames = m_p.getParameter<std::vector<std::string> >("DetectorNames");
0025   detLevels = m_p.getParameter<std::vector<int> >("DetectorLevels");
0026   etaRegions = m_p.getParameter<std::vector<double> >("EtaBoundaries");
0027   regionTypes = m_p.getParameter<std::vector<int> >("RegionTypes");
0028   boundaries = m_p.getParameter<std::vector<double> >("Boundaries");
0029   edm::LogInfo("MaterialBudget") << "MaterialBudget initialized for " << detTypes.size() << " detector types";
0030   unsigned int nc = 0;
0031   for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
0032     edm::LogInfo("MaterialBudget") << "Type[" << ii << "] : " << detTypes[ii] << " with " << constituents[ii]
0033                                    << ", order " << stackOrder[ii] << " constituents --> ";
0034     for (int kk = 0; kk < constituents[ii]; ++kk) {
0035       std::string name = "Unknown";
0036       int level = -1;
0037       if (nc < detNames.size()) {
0038         name = detNames[nc];
0039         level = detLevels[nc];
0040         ++nc;
0041       }
0042       edm::LogInfo("MaterialBudget") << "    Constituent[" << kk << "]: " << name << " at level " << level;
0043     }
0044   }
0045   edm::LogInfo("MaterialBudget") << "MaterialBudget Stop condition for " << etaRegions.size() << " eta regions";
0046   for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
0047     edm::LogInfo("MaterialBudget") << "Region[" << ii << "] : Eta < " << etaRegions[ii] << " boundary type "
0048                                    << regionTypes[ii] << " limit " << boundaries[ii];
0049   }
0050   book(m_p);
0051 }
0052 
0053 MaterialBudget::~MaterialBudget() {}
0054 
0055 void MaterialBudget::update(const BeginOfRun*) {
0056   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0057   std::vector<G4LogicalVolume*>::const_iterator lvcite;
0058 
0059   unsigned int kount = detNames.size();
0060   for (unsigned int ii = 0; ii < kount; ++ii)
0061     logVolumes.push_back(nullptr);
0062 
0063   for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0064     for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
0065       if (strcmp(detNames[ii].c_str(),
0066                  static_cast<std::string>(dd4hep::dd::noNamespace((*lvcite)->GetName())).c_str()) == 0) {
0067         logVolumes[ii] = (*lvcite);
0068         kount--;
0069         break;
0070       }
0071     }
0072     if (kount <= 0)
0073       break;
0074   }
0075   edm::LogInfo("MaterialBudget") << "MaterialBudget::Finds " << detNames.size() - kount << " out of " << detNames.size()
0076                                  << " LV addresses";
0077   for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
0078     std::string name("Unknown");
0079     if (logVolumes[ii] != nullptr)
0080       name = dd4hep::dd::noNamespace(logVolumes[ii]->GetName());
0081     edm::LogInfo("MaterialBudget") << "LV[" << ii << "] : " << detNames[ii] << " Address " << logVolumes[ii] << " | "
0082                                    << name;
0083   }
0084 
0085   for (unsigned int ii = 0; ii < (detTypes.size() + 1); ++ii) {
0086     stepLen.push_back(0);
0087     radLen.push_back(0);
0088     intLen.push_back(0);
0089   }
0090   stackOrder.push_back(0);
0091 }
0092 
0093 void MaterialBudget::update(const BeginOfTrack* trk) {
0094   for (unsigned int ii = 0; ii < (detTypes.size() + 1); ++ii) {
0095     stepLen[ii] = 0;
0096     radLen[ii] = 0;
0097     intLen[ii] = 0;
0098   }
0099 
0100   const G4Track* aTrack = (*trk)();  // recover G4 pointer if wanted
0101   const G4ThreeVector& mom = aTrack->GetMomentum();
0102   if (mom.theta() != 0) {
0103     eta_ = mom.eta();
0104   } else {
0105     eta_ = -99;
0106   }
0107   phi_ = mom.phi();
0108   stepT = 0;
0109 
0110 #ifdef EDM_ML_DEBUG
0111   double theEnergy = aTrack->GetTotalEnergy();
0112   int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
0113   edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Track " << aTrack->GetTrackID() << " Code " << theID
0114                                  << " Energy " << theEnergy / CLHEP::GeV << " GeV; Eta " << eta_ << " Phi "
0115                                  << phi_ / CLHEP::deg << " PT " << mom.perp() / CLHEP::GeV << " GeV *****";
0116 #endif
0117 }
0118 
0119 void MaterialBudget::update(const G4Step* aStep) {
0120   //---------- each step
0121   G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0122   double step = aStep->GetStepLength();
0123   double radl = material->GetRadlen();
0124   double intl = material->GetNuclearInterLength();
0125 
0126   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0127   unsigned int indx = detTypes.size();
0128   unsigned int nc = 0;
0129   for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
0130     for (int kk = 0; kk < constituents[ii]; ++kk) {
0131       if (detLevels[nc + kk] <= (int)((touch->GetHistoryDepth()) + 1)) {
0132         int jj = (int)((touch->GetHistoryDepth()) + 1) - detLevels[nc + kk];
0133         if ((touch->GetVolume(jj)->GetLogicalVolume()) == logVolumes[nc + kk]) {
0134           indx = ii;
0135           break;
0136         }
0137       }
0138     }
0139     nc += (unsigned int)(constituents[ii]);
0140     if (indx == ii)
0141       break;
0142   }
0143 
0144   stepT += step;
0145   stepLen[indx] = stepT;
0146   radLen[indx] += (step / radl);
0147   intLen[indx] += (step / intl);
0148 #ifdef EDM_ML_DEBUG
0149   edm::LogInfo("MaterialBudget") << "MaterialBudget::Step in "
0150                                  << dd4hep::dd::noNamespace(touch->GetVolume(0)->GetLogicalVolume()->GetName())
0151                                  << " Index " << indx << " Step " << step << " RadL " << step / radl << " IntL "
0152                                  << step / intl;
0153 #endif
0154   //----- Stop tracking after selected position
0155   if (stopAfter(aStep)) {
0156     G4Track* track = aStep->GetTrack();
0157     track->SetTrackStatus(fStopAndKill);
0158   }
0159 }
0160 
0161 void MaterialBudget::update(const EndOfTrack* trk) {
0162   for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
0163     for (unsigned int jj = 0; jj <= detTypes.size(); ++jj) {
0164       if (stackOrder[jj] == (int)(ii + 1)) {
0165         for (unsigned int kk = 0; kk <= detTypes.size(); ++kk) {
0166           if (stackOrder[kk] == (int)(ii)) {
0167             radLen[jj] += radLen[kk];
0168             intLen[jj] += intLen[kk];
0169 #ifdef EDM_ML_DEBUG
0170             edm::LogInfo("MaterialBudget")
0171                 << "MaterialBudget::Add " << kk << ":" << stackOrder[kk] << " to " << jj << ":" << stackOrder[jj]
0172                 << " RadL " << radLen[kk] << " : " << radLen[jj] << " IntL " << intLen[kk] << " : " << intLen[jj]
0173                 << " StepL " << stepLen[kk] << " : " << stepLen[jj];
0174 #endif
0175             break;
0176           }
0177         }
0178         break;
0179       }
0180     }
0181   }
0182 
0183   for (unsigned int ii = 0; ii <= detTypes.size(); ++ii) {
0184     me100[ii]->Fill(eta_, radLen[ii]);
0185     me200[ii]->Fill(eta_, intLen[ii]);
0186     me300[ii]->Fill(eta_, stepLen[ii]);
0187     me400[ii]->Fill(phi_, radLen[ii]);
0188     me500[ii]->Fill(phi_, intLen[ii]);
0189     me600[ii]->Fill(phi_, stepLen[ii]);
0190 #ifdef EDM_ML_DEBUG
0191     std::string name("Unknown");
0192     if (ii < detTypes.size())
0193       name = detTypes[ii];
0194     edm::LogInfo("MaterialBudget") << "MaterialBudget::Volume[" << ii << "]: " << name << " eta " << eta_ << " == Step "
0195                                    << stepLen[ii] << " RadL " << radLen[ii] << " IntL " << intLen[ii];
0196 #endif
0197   }
0198 }
0199 
0200 void MaterialBudget::book(const edm::ParameterSet& m_p) {
0201   // Book histograms
0202   edm::Service<TFileService> tfile;
0203 
0204   if (!tfile.isAvailable())
0205     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0206                                       << "please add it to config file";
0207 
0208   int binEta = m_p.getUntrackedParameter<int>("NBinEta", 320);
0209   int binPhi = m_p.getUntrackedParameter<int>("NBinPhi", 180);
0210   double minEta = m_p.getUntrackedParameter<double>("MinEta", -8.0);
0211   double maxEta = m_p.getUntrackedParameter<double>("MaxEta", 8.0);
0212   double maxPhi = CLHEP::pi;
0213   edm::LogInfo("MaterialBudget") << "MaterialBudget: Booking user "
0214                                  << "histos === with " << binEta << " bins "
0215                                  << "in eta from " << minEta << " to " << maxEta << " and " << binPhi << " bins "
0216                                  << "in phi from " << -maxPhi << " to " << maxPhi;
0217 
0218   char name[20], title[80];
0219   std::string named;
0220   for (unsigned int i = 0; i <= detTypes.size(); i++) {
0221     if (i >= detTypes.size())
0222       named = "Unknown";
0223     else
0224       named = detTypes[i];
0225     sprintf(name, "%d", i + 100);
0226     sprintf(title, "MB(X0) prof Eta in %s", named.c_str());
0227     me100.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
0228     sprintf(name, "%d", i + 200);
0229     sprintf(title, "MB(L0) prof Eta in %s", named.c_str());
0230     me200.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
0231     sprintf(name, "%d", i + 300);
0232     sprintf(title, "MB(Step) prof Eta in %s", named.c_str());
0233     me300.push_back(tfile->make<TProfile>(name, title, binEta, minEta, maxEta));
0234     sprintf(name, "%d", i + 400);
0235     sprintf(title, "MB(X0) prof Phi in %s", named.c_str());
0236     me400.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
0237     sprintf(name, "%d", i + 500);
0238     sprintf(title, "MB(L0) prof Phi in %s", named.c_str());
0239     me500.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
0240     sprintf(name, "%d", i + 600);
0241     sprintf(title, "MB(Step) prof Phi in %s", named.c_str());
0242     me600.push_back(tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi));
0243   }
0244 
0245   edm::LogInfo("MaterialBudget") << "MaterialBudget: Booking user "
0246                                  << "histos done ===";
0247 }
0248 
0249 bool MaterialBudget::stopAfter(const G4Step* aStep) {
0250   G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0251   if (aStep->GetPostStepPoint() != nullptr)
0252     hitPoint = aStep->GetPostStepPoint()->GetPosition();
0253   double rr = hitPoint.perp();
0254   double zz = std::abs(hitPoint.z());
0255 
0256   bool flag(false);
0257   for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
0258 #ifdef EDM_ML_DEBUG
0259     edm::LogInfo("MaterialBudget") << " MaterialBudget::Eta " << eta_ << " in Region[" << ii << "] with "
0260                                    << etaRegions[ii] << " type " << regionTypes[ii] << "|" << boundaries[ii];
0261 #endif
0262     if (fabs(eta_) < etaRegions[ii]) {
0263       if (regionTypes[ii] == 0) {
0264         if (rr >= boundaries[ii] - 0.001)
0265           flag = true;
0266       } else {
0267         if (zz >= boundaries[ii] - 0.001)
0268           flag = true;
0269       }
0270 #ifdef EDM_ML_DEBUG
0271       if (flag)
0272         edm::LogInfo("MaterialBudget") << " MaterialBudget::Stop after R = " << rr << " and Z = " << zz;
0273 #endif
0274       break;
0275     }
0276   }
0277 #ifdef EDM_ML_DEBUG
0278   edm::LogInfo("MaterialBudget") << " MaterialBudget:: R = " << rr << " and Z = " << zz << " Flag " << flag;
0279 #endif
0280   return flag;
0281 }