Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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;  // stop default
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   // observer classes
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);