Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:32:37

0001 #include "Validation/Geometry/interface/MaterialBudgetForward.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 
0016 #include <iostream>
0017 //#define DebugLog
0018 const int MaterialBudgetForward::maxSet;
0019 
0020 MaterialBudgetForward::MaterialBudgetForward(const edm::ParameterSet& p) {
0021   edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudgetForward");
0022   detTypes = m_p.getParameter<std::vector<std::string> >("DetectorTypes");
0023   constituents = m_p.getParameter<std::vector<int> >("Constituents");
0024   stackOrder = m_p.getParameter<std::vector<int> >("StackOrder");
0025   detNames = m_p.getParameter<std::vector<std::string> >("DetectorNames");
0026   detLevels = m_p.getParameter<std::vector<int> >("DetectorLevels");
0027   etaRegions = m_p.getParameter<std::vector<double> >("EtaBoundaries");
0028   regionTypes = m_p.getParameter<std::vector<int> >("RegionTypes");
0029   boundaries = m_p.getParameter<std::vector<double> >("Boundaries");
0030   edm::LogInfo("MaterialBudget") << "MaterialBudgetForward initialized for " << detTypes.size() << " detector types";
0031   unsigned int nc = 0;
0032   for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
0033     edm::LogInfo("MaterialBudget") << "Type[" << ii << "] : " << detTypes[ii] << " with " << constituents[ii]
0034                                    << ", order " << stackOrder[ii] << " constituents --> ";
0035     for (int kk = 0; kk < constituents[ii]; ++kk) {
0036       std::string name = "Unknown";
0037       int level = -1;
0038       if (nc < detNames.size()) {
0039         name = detNames[nc];
0040         level = detLevels[nc];
0041         ++nc;
0042       }
0043       edm::LogInfo("MaterialBudget") << "    Constituent[" << kk << "]: " << name << " at level " << level;
0044     }
0045   }
0046   edm::LogInfo("MaterialBudget") << "MaterialBudgetForward Stop condition for " << etaRegions.size() << " eta regions";
0047   for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
0048     edm::LogInfo("MaterialBudget") << "Region[" << ii << "] : Eta < " << etaRegions[ii] << " boundary type "
0049                                    << regionTypes[ii] << " limit " << boundaries[ii];
0050   }
0051   book(m_p);
0052 }
0053 
0054 MaterialBudgetForward::~MaterialBudgetForward() {}
0055 
0056 void MaterialBudgetForward::update(const BeginOfRun*) {
0057   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0058   std::vector<G4LogicalVolume*>::const_iterator lvcite;
0059 
0060   unsigned int kount = detNames.size();
0061   for (unsigned int ii = 0; ii < kount; ++ii)
0062     logVolumes.push_back(nullptr);
0063 
0064   for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0065     for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
0066       if (strcmp(detNames[ii].c_str(), (*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") << "MaterialBudgetForward::Finds " << detNames.size() - kount << " out of "
0076                                  << detNames.size() << " LV addresses";
0077   for (unsigned int ii = 0; ii < detNames.size(); ++ii) {
0078     std::string name("Unknown");
0079     if (logVolumes[ii] != nullptr)
0080       name = 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 MaterialBudgetForward::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 DebugLog
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 MaterialBudgetForward::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 DebugLog
0149   edm::LogInfo("MaterialBudget") << "MaterialBudgetForward::Step in "
0150                                  << touch->GetVolume(0)->GetLogicalVolume()->GetName() << " Index " << indx << " Step "
0151                                  << step << " RadL " << step / radl << " IntL " << step / intl;
0152 #endif
0153   //----- Stop tracking after selected position
0154   if (stopAfter(aStep)) {
0155     G4Track* track = aStep->GetTrack();
0156     track->SetTrackStatus(fStopAndKill);
0157   }
0158 }
0159 
0160 void MaterialBudgetForward::update(const EndOfTrack* trk) {
0161   for (unsigned int ii = 0; ii < detTypes.size(); ++ii) {
0162     for (unsigned int jj = 0; jj <= detTypes.size(); ++jj) {
0163       if (stackOrder[jj] == (int)(ii + 1)) {
0164         for (unsigned int kk = 0; kk <= detTypes.size(); ++kk) {
0165           if (stackOrder[kk] == (int)(ii)) {
0166             radLen[jj] += radLen[kk];
0167             intLen[jj] += intLen[kk];
0168 #ifdef DebugLog
0169             edm::LogInfo("MaterialBudget")
0170                 << "MaterialBudgetForward::Add " << kk << ":" << stackOrder[kk] << " to " << jj << ":" << stackOrder[jj]
0171                 << " RadL " << radLen[kk] << " : " << radLen[jj] << " IntL " << intLen[kk] << " : " << intLen[jj]
0172                 << " StepL " << stepLen[kk] << " : " << stepLen[jj];
0173 #endif
0174             break;
0175           }
0176         }
0177         break;
0178       }
0179     }
0180   }
0181 
0182   for (unsigned int ii = 0; ii <= detTypes.size(); ++ii) {
0183     me100[ii]->Fill(eta_, radLen[ii]);
0184     me200[ii]->Fill(eta_, intLen[ii]);
0185     me300[ii]->Fill(eta_, stepLen[ii]);
0186     me400[ii]->Fill(eta_);
0187     me500[ii]->Fill(eta_, phi_, radLen[ii]);
0188     me600[ii]->Fill(eta_, phi_, intLen[ii]);
0189     me700[ii]->Fill(eta_, phi_, stepLen[ii]);
0190     me800[ii]->Fill(eta_, phi_);
0191 
0192     std::string name("Unknown");
0193     if (ii < detTypes.size())
0194       name = detTypes[ii];
0195 #ifdef DebugLog
0196     edm::LogInfo("MaterialBudget") << "MaterialBudgetForward::Volume[" << ii << "]: " << name << " eta " << eta_
0197                                    << " == Step " << stepLen[ii] << " RadL " << radLen[ii] << " IntL " << intLen[ii];
0198 #endif
0199   }
0200 }
0201 
0202 void MaterialBudgetForward::book(const edm::ParameterSet& m_p) {
0203   // Book histograms
0204   edm::Service<TFileService> tfile;
0205 
0206   if (!tfile.isAvailable())
0207     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0208                                       << "please add it to config file";
0209 
0210   int binEta = m_p.getUntrackedParameter<int>("NBinEta", 320);
0211   int binPhi = m_p.getUntrackedParameter<int>("NBinPhi", 180);
0212   double minEta = m_p.getUntrackedParameter<double>("MinEta", -8.0);
0213   double maxEta = m_p.getUntrackedParameter<double>("MaxEta", 8.0);
0214   double maxPhi = CLHEP::pi;
0215   edm::LogInfo("MaterialBudget") << "MaterialBudgetForward: Booking user "
0216                                  << "histos === with " << binEta << " bins "
0217                                  << "in eta from " << minEta << " to " << maxEta << " and " << binPhi << " bins "
0218                                  << "in phi from " << -maxPhi << " to " << maxPhi;
0219 
0220   char name[20], title[80];
0221   std::string named;
0222   for (int i = 0; i < std::min((int)(detTypes.size() + 1), maxSet); i++) {
0223     if (i >= (int)(detTypes.size()))
0224       named = "Unknown";
0225     else
0226       named = detTypes[i];
0227     sprintf(name, "%d", i + 100);
0228     sprintf(title, "MB(X0) prof Eta in %s", named.c_str());
0229     me100[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
0230     sprintf(name, "%d", i + 200);
0231     sprintf(title, "MB(L0) prof Eta in %s", named.c_str());
0232     me200[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
0233     sprintf(name, "%d", i + 300);
0234     sprintf(title, "MB(Step) prof Eta in %s", named.c_str());
0235     me300[i] = tfile->make<TProfile>(name, title, binEta, minEta, maxEta);
0236     sprintf(name, "%d", i + 400);
0237     sprintf(title, "Eta in %s", named.c_str());
0238     me400[i] = tfile->make<TH1F>(name, title, binEta, minEta, maxEta);
0239     sprintf(name, "%d", i + 500);
0240     sprintf(title, "MB(X0) prof Eta Phi in %s", named.c_str());
0241     me500[i] = tfile->make<TProfile2D>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
0242     sprintf(name, "%d", i + 600);
0243     sprintf(title, "MB(L0) prof Eta Phi in %s", named.c_str());
0244     me600[i] = tfile->make<TProfile2D>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
0245     sprintf(name, "%d", i + 700);
0246     sprintf(title, "MB(Step) prof Eta Phi in %s", named.c_str());
0247     me700[i] = tfile->make<TProfile2D>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
0248     sprintf(name, "%d", i + 800);
0249     sprintf(title, "Eta vs Phi in %s", named.c_str());
0250     me800[i] = tfile->make<TH2F>(name, title, binEta / 2, minEta, maxEta, binPhi / 2, -maxPhi, maxPhi);
0251   }
0252 
0253   edm::LogInfo("MaterialBudget") << "MaterialBudgetForward: Booking user "
0254                                  << "histos done ===";
0255 }
0256 
0257 bool MaterialBudgetForward::stopAfter(const G4Step* aStep) {
0258   G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0259   if (aStep->GetPostStepPoint() != nullptr)
0260     hitPoint = aStep->GetPostStepPoint()->GetPosition();
0261   double rr = hitPoint.perp();
0262   double zz = std::abs(hitPoint.z());
0263 
0264   bool flag(false);
0265   for (unsigned int ii = 0; ii < etaRegions.size(); ++ii) {
0266 #ifdef DebugLog
0267     edm::LogInfo("MaterialBudget") << " MaterialBudgetForward::Eta " << eta_ << " in Region[" << ii << "] with "
0268                                    << etaRegions[ii] << " type " << regionTypes[ii] << "|" << boundaries[ii];
0269 #endif
0270     if (fabs(eta_) < etaRegions[ii]) {
0271       if (regionTypes[ii] == 0) {
0272         if (rr >= boundaries[ii] - 0.001)
0273           flag = true;
0274       } else {
0275         if (zz >= boundaries[ii] - 0.001)
0276           flag = true;
0277       }
0278 #ifdef DebugLog
0279       if (flag)
0280         edm::LogInfo("MaterialBudget") << " MaterialBudgetForward::Stop after R = " << rr << " and Z = " << zz;
0281 #endif
0282       break;
0283     }
0284   }
0285 #ifdef DebugLog
0286   edm::LogInfo("MaterialBudget") << " MaterialBudgetForward:: R = " << rr << " and Z = " << zz << " Flag " << flag;
0287 #endif
0288   return flag;
0289 }