File indexing completed on 2021-11-09 23:35:37
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
0067 private:
0068 std::vector<std::string> lvNames_;
0069 std::vector<int> lvLevel_;
0070 int iaddLevel_;
0071 bool init_;
0072 std::map<int, std::pair<G4LogicalVolume*, int> > mapLV_;
0073 std::vector<MatInfo> lengths_;
0074 std::vector<MaterialInformation> store_;
0075 };
0076
0077 MaterialBudgetVolume::MaterialBudgetVolume(const edm::ParameterSet& p) : init_(false) {
0078 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("MaterialBudgetVolume");
0079
0080 lvNames_ = m_p.getParameter<std::vector<std::string> >("lvNames");
0081 lvLevel_ = m_p.getParameter<std::vector<int> >("lvLevels");
0082 iaddLevel_ = (m_p.getParameter<bool>("useDD4hep")) ? 1 : 0;
0083
0084 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: Studies Material budget for " << lvNames_.size()
0085 << " volumes with addLevel " << iaddLevel_;
0086 std::ostringstream st1;
0087 for (unsigned int k = 0; k < lvNames_.size(); ++k)
0088 st1 << " [" << k << "] " << lvNames_[k] << " at " << lvLevel_[k];
0089 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: Volumes" << st1.str();
0090
0091 produces<edm::MaterialInformationContainer>("MaterialInformation");
0092 #ifdef EDM_ML_DEBUG
0093 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: will produce MaterialInformationContainer";
0094 #endif
0095 }
0096
0097 void MaterialBudgetVolume::produce(edm::Event& e, const edm::EventSetup&) {
0098 std::unique_ptr<edm::MaterialInformationContainer> matbg(new edm::MaterialInformationContainer);
0099 endOfEvent(*matbg);
0100 e.put(std::move(matbg), "MaterialInformation");
0101 }
0102
0103 void MaterialBudgetVolume::update(const BeginOfRun* run) {
0104 init_ = loadLV();
0105
0106 #ifdef EDM_ML_DEBUG
0107 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Finds " << mapLV_.size()
0108 << " logical volumes with return flag " << init_;
0109 #endif
0110 }
0111
0112 void MaterialBudgetVolume::update(const BeginOfEvent* evt) {
0113 #ifdef EDM_ML_DEBUG
0114 int iev = (*evt)()->GetEventID();
0115 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: =====> Begin event = " << iev << std::endl;
0116 #endif
0117 if (!init_)
0118 init_ = loadLV();
0119
0120 store_.clear();
0121 }
0122
0123 void MaterialBudgetVolume::update(const BeginOfTrack* trk) {
0124 lengths_ = std::vector<MatInfo>(mapLV_.size());
0125
0126 #ifdef EDM_ML_DEBUG
0127 const G4Track* aTrack = (*trk)();
0128 const G4ThreeVector& mom = aTrack->GetMomentum();
0129 double theEnergy = aTrack->GetTotalEnergy();
0130 int theID = static_cast<int>(aTrack->GetDefinition()->GetPDGEncoding());
0131 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolumme: Track " << aTrack->GetTrackID() << " Code " << theID
0132 << " Energy " << theEnergy / CLHEP::GeV << " GeV; Momentum " << mom / CLHEP::GeV
0133 << " GeV";
0134 #endif
0135 }
0136
0137 void MaterialBudgetVolume::update(const G4Step* aStep) {
0138 G4Material* material = aStep->GetPreStepPoint()->GetMaterial();
0139 double step = aStep->GetStepLength();
0140 double radl = material->GetRadlen();
0141 double intl = material->GetNuclearInterLength();
0142 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0143 int index = findLV(touch);
0144 if (index >= 0) {
0145 lengths_[index].stepL += step;
0146 lengths_[index].radL += (step / radl);
0147 lengths_[index].intL += (step / intl);
0148 }
0149 #ifdef EDM_ML_DEBUG
0150 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Step in "
0151 << touch->GetVolume(0)->GetLogicalVolume()->GetName() << " Index " << index
0152 << " Step " << step << " RadL " << step / radl << " IntL " << step / intl;
0153 #endif
0154 }
0155
0156 void MaterialBudgetVolume::update(const EndOfTrack* trk) {
0157 const G4Track* aTrack = (*trk)();
0158 int id = aTrack->GetTrackID();
0159 double eta = aTrack->GetMomentumDirection().eta();
0160 double phi = aTrack->GetMomentumDirection().phi();
0161 for (unsigned int k = 0; k < lengths_.size(); ++k) {
0162 MaterialInformation info(lvNames_[k], id, eta, phi, lengths_[k].stepL, lengths_[k].radL, lengths_[k].intL);
0163 store_.emplace_back(info);
0164 #ifdef EDM_ML_DEBUG
0165 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Volume[" << k << "]: " << info;
0166 #endif
0167 }
0168 }
0169
0170 void MaterialBudgetVolume::endOfEvent(edm::MaterialInformationContainer& matbg) {
0171 #ifdef EDM_ML_DEBUG
0172 unsigned int kount(0);
0173 #endif
0174 for (const auto& element : store_) {
0175 MaterialInformation info(element.vname(),
0176 element.id(),
0177 element.trackEta(),
0178 element.trackPhi(),
0179 element.stepLength(),
0180 element.radiationLength(),
0181 element.interactionLength());
0182 matbg.push_back(info);
0183 #ifdef EDM_ML_DEBUG
0184 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume:: Info[" << kount << "] " << info;
0185 ++kount;
0186 #endif
0187 }
0188 }
0189
0190 bool MaterialBudgetVolume::loadLV() {
0191 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0192 bool flag(false);
0193 if (lvs != nullptr) {
0194 std::vector<G4LogicalVolume*>::const_iterator lvcite;
0195 for (unsigned int i = 0; i < lvNames_.size(); i++) {
0196 G4LogicalVolume* lv = nullptr;
0197 std::string name(lvNames_[i]);
0198 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0199 std::string namx(dd4hep::dd::noNamespace((*lvcite)->GetName()));
0200 if (namx == name) {
0201 lv = (*lvcite);
0202 break;
0203 }
0204 }
0205 if (lv != nullptr)
0206 mapLV_[i] = std::make_pair(lv, (lvLevel_[i] + iaddLevel_));
0207 }
0208 flag = true;
0209 #ifdef EDM_ML_DEBUG
0210 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume::Finds " << mapLV_.size() << " logical volumes";
0211 unsigned int k(0);
0212 for (const auto& lvs : mapLV_) {
0213 edm::LogVerbatim("MaterialBudget") << "Entry[" << k << "] " << lvs.first << ": (" << (lvs.second).first << ", "
0214 << (lvs.second).second << ") : " << lvNames_[lvs.first];
0215 ++k;
0216 }
0217 #endif
0218 }
0219 return flag;
0220 }
0221
0222 int MaterialBudgetVolume::findLV(const G4VTouchable* touch) {
0223 int level(-1);
0224 int levels = ((touch->GetHistoryDepth()) + 1);
0225 for (const auto& lvs : mapLV_) {
0226 if ((lvs.second).second <= levels) {
0227 int ii = levels - (lvs.second).second;
0228 if ((touch->GetVolume(ii)->GetLogicalVolume()) == (lvs.second).first) {
0229 level = lvs.first;
0230 break;
0231 }
0232 }
0233 }
0234 #ifdef EDM_ML_DEBUG
0235 if (level < 0) {
0236 edm::LogVerbatim("MaterialBudget") << "findLV: Gets " << level << " from " << levels << " levels in touchables";
0237 for (int i = 0; i < levels; ++i)
0238 edm::LogVerbatim("MaterialBudget") << "[" << (levels - i) << "] " << touch->GetVolume(i)->GetLogicalVolume()
0239 << " : " << touch->GetVolume(i)->GetLogicalVolume()->GetName();
0240 }
0241 #endif
0242 return level;
0243 }
0244
0245 #include "FWCore/Framework/interface/MakerMacros.h"
0246 #include "FWCore/PluginManager/interface/ModuleDef.h"
0247 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0248
0249 DEFINE_SIMWATCHER(MaterialBudgetVolume);
