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