File indexing completed on 2023-11-14 23:29:33
0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005
0006 #include "SimDataFormats/CaloHit/interface/MaterialInformation.h"
0007
0008 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0009 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0010 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0011 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0012 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0013 #include "SimG4Core/Notification/interface/Observer.h"
0014 #include "SimG4Core/Watcher/interface/SimProducer.h"
0015
0016 #include "G4LogicalVolumeStore.hh"
0017 #include "G4Step.hh"
0018 #include "G4Track.hh"
0019 #include "G4TransportationManager.hh"
0020 #include "DD4hep/Filter.h"
0021
0022 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0023 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0024
0025 #include <cmath>
0026 #include <iomanip>
0027 #include <iostream>
0028 #include <map>
0029 #include <memory>
0030 #include <string>
0031 #include <vector>
0032 #include <utility>
0033
0034
0035
0036 class MaterialBudgetVolume : public SimProducer,
0037 public Observer<const BeginOfRun*>,
0038 public Observer<const BeginOfEvent*>,
0039 public Observer<const BeginOfTrack*>,
0040 public Observer<const G4Step*>,
0041 public Observer<const EndOfTrack*> {
0042 public:
0043 MaterialBudgetVolume(const edm::ParameterSet& p);
0044 MaterialBudgetVolume(const MaterialBudgetVolume&) = delete;
0045 const MaterialBudgetVolume& operator=(const MaterialBudgetVolume&) = delete;
0046 ~MaterialBudgetVolume() override {}
0047
0048 void produce(edm::Event&, const edm::EventSetup&) override;
0049
0050 struct MatInfo {
0051 double stepL, radL, intL;
0052 MatInfo(double s = 0, double r = 0, double l = 0) : stepL(s), radL(r), intL(l) {}
0053 };
0054
0055 private:
0056
0057 void update(const BeginOfRun* run) override;
0058 void update(const BeginOfEvent* evt) override;
0059 void update(const BeginOfTrack*) override;
0060 void update(const EndOfTrack*) override;
0061 void update(const G4Step* step) override;
0062
0063 void endOfEvent(edm::MaterialInformationContainer& matbg);
0064 bool loadLV();
0065 int findLV(const G4VTouchable*);
0066 bool stopAfter(const G4Step*);
0067
0068 private:
0069 std::vector<std::string> lvNames_;
0070 std::vector<int> lvLevel_;
0071 int iaddLevel_;
0072 double rMax_, zMax_;
0073 bool init_;
0074 std::map<int, std::pair<G4LogicalVolume*, int> > mapLV_;
0075 std::vector<MatInfo> lengths_;
0076 std::vector<MaterialInformation> store_;
0077 };
0078
0079 MaterialBudgetVolume::MaterialBudgetVolume(const edm::ParameterSet& p) : init_(false) {
0080 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudgetVolume");
0081
0082 lvNames_ = m_p.getParameter<std::vector<std::string> >("lvNames");
0083 lvLevel_ = m_p.getParameter<std::vector<int> >("lvLevels");
0084 iaddLevel_ = (m_p.getParameter<bool>("useDD4hep")) ? 1 : 0;
0085 rMax_ = m_p.getParameter<double>("rMax") * CLHEP::m;
0086 zMax_ = m_p.getParameter<double>("zMax") * CLHEP::m;
0087
0088 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: Studies Material budget for " << lvNames_.size()
0089 << " volumes with addLevel " << iaddLevel_ << " and with rMax " << rMax_
0090 << " mm; zMax " << zMax_ << " mm";
0091 std::ostringstream st1;
0092 for (unsigned int k = 0; k < lvNames_.size(); ++k)
0093 st1 << " [" << k << "] " << lvNames_[k] << " at " << lvLevel_[k];
0094 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: Volumes" << st1.str();
0095
0096 produces<edm::MaterialInformationContainer>("MaterialInformation");
0097 #ifdef EDM_ML_DEBUG
0098 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: will produce MaterialInformationContainer";
0099 #endif
0100 }
0101
0102 void MaterialBudgetVolume::produce(edm::Event& e, const edm::EventSetup&) {
0103 std::unique_ptr<edm::MaterialInformationContainer> matbg(new edm::MaterialInformationContainer);
0104 endOfEvent(*matbg);
0105 e.put(std::move(matbg), "MaterialInformation");
0106 }
0107
0108 void MaterialBudgetVolume::update(const BeginOfRun* run) {
0109 init_ = loadLV();
0110
0111 #ifdef EDM_ML_DEBUG
0112 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Finds " << mapLV_.size()
0113 << " logical volumes with return flag " << init_;
0114 #endif
0115 }
0116
0117 void MaterialBudgetVolume::update(const BeginOfEvent* evt) {
0118 #ifdef EDM_ML_DEBUG
0119 int iev = (*evt)()->GetEventID();
0120 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: =====> Begin event = " << iev << std::endl;
0121 #endif
0122 if (!init_)
0123 init_ = loadLV();
0124
0125 store_.clear();
0126 }
0127
0128 void MaterialBudgetVolume::update(const BeginOfTrack* trk) {
0129 lengths_ = std::vector<MatInfo>(mapLV_.size());
0130
0131 #ifdef EDM_ML_DEBUG
0132 const G4Track* aTrack = (*trk)();
0133 const G4ThreeVector& mom = aTrack->GetMomentum();
0134 double theEnergy = aTrack->GetTotalEnergy();
0135 int theID = static_cast<int>(aTrack->GetDefinition()->GetPDGEncoding());
0136 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolumme: Track " << aTrack->GetTrackID() << " Code " << theID
0137 << " Energy " << theEnergy / CLHEP::GeV << " GeV; Momentum " << mom / CLHEP::GeV
0138 << " GeV";
0139 #endif
0140 }
0141
0142 void MaterialBudgetVolume::update(const G4Step* aStep) {
0143 G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0144 double step = aStep->GetStepLength();
0145 double radl = material->GetRadlen();
0146 double intl = material->GetNuclearInterLength();
0147 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0148 int index = findLV(touch);
0149 if (index >= 0) {
0150 lengths_[index].stepL += step;
0151 lengths_[index].radL += (step / radl);
0152 lengths_[index].intL += (step / intl);
0153 }
0154 #ifdef EDM_ML_DEBUG
0155 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Step in "
0156 << touch->GetVolume(0)->GetLogicalVolume()->GetName() << " Index " << index
0157 << " Step " << step << " RadL " << step / radl << " IntL " << step / intl;
0158 #endif
0159
0160
0161 if (stopAfter(aStep)) {
0162 G4Track* track = aStep->GetTrack();
0163 track->SetTrackStatus(fStopAndKill);
0164 }
0165 }
0166
0167 void MaterialBudgetVolume::update(const EndOfTrack* trk) {
0168 const G4Track* aTrack = (*trk)();
0169 int id = aTrack->GetTrackID();
0170 double eta = aTrack->GetMomentumDirection().eta();
0171 double phi = aTrack->GetMomentumDirection().phi();
0172 for (unsigned int k = 0; k < lengths_.size(); ++k) {
0173 MaterialInformation info(lvNames_[k], id, eta, phi, lengths_[k].stepL, lengths_[k].radL, lengths_[k].intL);
0174 store_.emplace_back(info);
0175 #ifdef EDM_ML_DEBUG
0176 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Volume[" << k << "]: " << info;
0177 #endif
0178 }
0179 }
0180
0181 void MaterialBudgetVolume::endOfEvent(edm::MaterialInformationContainer& matbg) {
0182 #ifdef EDM_ML_DEBUG
0183 unsigned int kount(0);
0184 #endif
0185 for (const auto& element : store_) {
0186 MaterialInformation info(element.vname(),
0187 element.id(),
0188 element.trackEta(),
0189 element.trackPhi(),
0190 element.stepLength(),
0191 element.radiationLength(),
0192 element.interactionLength());
0193 matbg.push_back(info);
0194 #ifdef EDM_ML_DEBUG
0195 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume:: Info[" << kount << "] " << info;
0196 ++kount;
0197 #endif
0198 }
0199 }
0200
0201 bool MaterialBudgetVolume::loadLV() {
0202 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0203 bool flag(false);
0204 if (lvs != nullptr) {
0205 std::vector<G4LogicalVolume*>::const_iterator lvcite;
0206 for (unsigned int i = 0; i < lvNames_.size(); i++) {
0207 G4LogicalVolume* lv = nullptr;
0208 std::string name(lvNames_[i]);
0209 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0210 std::string namx(dd4hep::dd::noNamespace((*lvcite)->GetName()));
0211 if (namx == name) {
0212 lv = (*lvcite);
0213 break;
0214 }
0215 }
0216 if (lv != nullptr)
0217 mapLV_[i] = std::make_pair(lv, (lvLevel_[i] + iaddLevel_));
0218 }
0219 flag = true;
0220 #ifdef EDM_ML_DEBUG
0221 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Finds " << mapLV_.size() << " logical volumes";
0222 unsigned int k(0);
0223 for (const auto& lvs : mapLV_) {
0224 edm::LogVerbatim("MaterialBudget") << "Entry[" << k << "] " << lvs.first << ": (" << (lvs.second).first << ", "
0225 << (lvs.second).second << ") : " << lvNames_[lvs.first];
0226 ++k;
0227 }
0228 #endif
0229 }
0230 return flag;
0231 }
0232
0233 int MaterialBudgetVolume::findLV(const G4VTouchable* touch) {
0234 int level(-1);
0235 int levels = ((touch->GetHistoryDepth()) + 1);
0236 for (const auto& lvs : mapLV_) {
0237 if ((lvs.second).second <= levels) {
0238 int ii = levels - (lvs.second).second;
0239 if ((touch->GetVolume(ii)->GetLogicalVolume()) == (lvs.second).first) {
0240 level = lvs.first;
0241 break;
0242 }
0243 }
0244 }
0245 #ifdef EDM_ML_DEBUG
0246 if (level < 0) {
0247 edm::LogVerbatim("MaterialBudget") << "findLV: Gets " << level << " from " << levels << " levels in touchables";
0248 for (int i = 0; i < levels; ++i)
0249 edm::LogVerbatim("MaterialBudget") << "[" << (levels - i) << "] " << touch->GetVolume(i)->GetLogicalVolume()
0250 << " : " << touch->GetVolume(i)->GetLogicalVolume()->GetName();
0251 }
0252 #endif
0253 return level;
0254 }
0255
0256 bool MaterialBudgetVolume::stopAfter(const G4Step* aStep) {
0257 bool flag(false);
0258 if ((rMax_ > 0.0) && (zMax_ > 0.0)) {
0259 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0260 double rr = hitPoint.perp();
0261 double zz = std::abs(hitPoint.z());
0262
0263 if (rr > rMax_ || zz > zMax_) {
0264 #ifdef EDM_ML_DEBUG
0265 edm::LogVerbatim("MaterialBudget") << " MaterialBudgetHcal::StopAfter R = " << rr << " and Z = " << zz;
0266 #endif
0267 flag = true;
0268 }
0269 }
0270 return flag;
0271 }
0272
0273 #include "FWCore/Framework/interface/MakerMacros.h"
0274 #include "FWCore/PluginManager/interface/ModuleDef.h"
0275 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0276
0277 DEFINE_SIMWATCHER(MaterialBudgetVolume);
