File indexing completed on 2024-04-06 12:32: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)();
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
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
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
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 }
