Forward/plugins/SimG4FluxProducer.cc

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/CaloTest/interface/ParticleFlux.h"
0007
0008 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0009 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0010 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0011 #include "SimG4Core/Watcher/interface/SimProducer.h"
0012 #include "SimG4Core/Notification/interface/Observer.h"
0013
0014 #include "G4Step.hh"
0015 #include "G4LogicalVolumeStore.hh"
0016 #include "G4PhysicalVolumeStore.hh"
0017 #include "G4Track.hh"
0018 #include "G4TouchableHistory.hh"
0019 #include "G4TransportationManager.hh"
0020
0021 #include <CLHEP/Units/SystemOfUnits.h>
0022 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0023
0024 #include <cmath>
0025 #include <iostream>
0026 #include <iomanip>
0027 #include <map>
0028 #include <memory>
0029 #include <string>
0030 #include <utility>
0031 #include <vector>
0032
0033 //#define EDM_ML_DEBUG
0034
0035 class SimG4FluxProducer : public SimProducer,
0036                           public Observer<const BeginOfRun *>,
0037                           public Observer<const BeginOfEvent *>,
0038                           public Observer<const G4Step *> {
0039 public:
0040   SimG4FluxProducer(const edm::ParameterSet &p);
0041   SimG4FluxProducer(const SimG4FluxProducer &) = delete;  // stop default
0042   const SimG4FluxProducer &operator=(const SimG4FluxProducer &) = delete;
0043   ~SimG4FluxProducer() override;
0044
0045   void produce(edm::Event &, const edm::EventSetup &) override;
0046
0047 private:
0048   // observer classes
0049   void update(const BeginOfRun *run) override;
0050   void update(const BeginOfEvent *evt) override;
0051   void update(const G4Step *step) override;
0052
0053   void endOfEvent(ParticleFlux &pflx, unsigned int k);
0054   G4VPhysicalVolume *getTopPV();
0055   std::map<G4LogicalVolume *, std::pair<unsigned int, std::string>>::iterator findLV(G4LogicalVolume *plv);
0056
0057   std::vector<std::string> LVNames_;
0058   std::vector<int> LVTypes_;
0059   G4VPhysicalVolume *topPV_;
0060   std::map<G4LogicalVolume *, std::pair<unsigned int, std::string>> mapLV_;
0061
0062   // some private members for ananlysis
0063   unsigned int count_;
0064   bool init_;
0065   std::map<std::pair<G4LogicalVolume *, unsigned int>, ParticleFlux> store_;
0066 };
0067
0068 SimG4FluxProducer::SimG4FluxProducer(const edm::ParameterSet &p) : count_(0), init_(false) {
0069   edm::ParameterSet m_FP = p.getParameter<edm::ParameterSet>("SimG4FluxProducer");
0070   LVNames_ = m_FP.getUntrackedParameter<std::vector<std::string>>("LVNames");
0071   LVTypes_ = m_FP.getUntrackedParameter<std::vector<int>>("LVTypes");
0072
0073   for (unsigned int k = 0; k < LVNames_.size(); ++k) {
0074     produces<ParticleFlux>(Form("%sParticleFlux", LVNames_[k].c_str()));
0075 #ifdef EDM_ML_DEBUG
0076     edm::LogVerbatim("SimG4FluxProducer")
0077         << "SimG4FluxProducer::Collection name[" << k << "] ParticleFlux" << LVNames_[k] << " and type " << LVTypes_[k];
0078 #endif
0079   }
0080 }
0081
0082 SimG4FluxProducer::~SimG4FluxProducer() {}
0083
0084 void SimG4FluxProducer::produce(edm::Event &e, const edm::EventSetup &) {
0085 #ifdef EDM_ML_DEBUG
0086   edm::LogVerbatim("SimG4FluxProducer") << "SimG4FluxProducer: enters produce with " << LVNames_.size() << " LV's";
0087 #endif
0088   for (unsigned int k = 0; k < LVNames_.size(); ++k) {
0089     std::unique_ptr<ParticleFlux> pflux(new ParticleFlux);
0090     endOfEvent(*pflux, k);
0091     std::string name = LVNames_[k] + "ParticleFlux";
0092     e.put(std::move(pflux), name);
0093   }
0094 }
0095
0096 void SimG4FluxProducer::update(const BeginOfRun *run) {
0097 #ifdef EDM_ML_DEBUG
0098   edm::LogVerbatim("SimG4FluxProducer") << "SimG4FluxProducer: enters BeginOfRun";
0099 #endif
0100   topPV_ = getTopPV();
0101   if (topPV_ == nullptr) {
0102     edm::LogWarning("SimG4FluxProducer") << "Cannot find top level volume\n";
0103   } else {
0104     init_ = true;
0105     const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
0106     for (auto lvcite : *lvs) {
0107       findLV(lvcite);
0108     }
0109
0110 #ifdef EDM_ML_DEBUG
0111     edm::LogVerbatim("SimG4FluxProducer") << "SimG4FluxProducer::Finds " << mapLV_.size() << " logical volumes";
0112     unsigned int k(0);
0113     for (const auto &lvs : mapLV_) {
0114       edm::LogVerbatim("SimG4FluxProducer")
0115           << "Entry[" << k << "] " << lvs.first << ": (" << (lvs.second).first << ", " << (lvs.second).second << ")";
0116       ++k;
0117     }
0118 #endif
0119   }
0120 }
0121
0122 void SimG4FluxProducer::update(const BeginOfEvent *evt) {
0123   int iev = (*evt)()->GetEventID();
0124   edm::LogVerbatim("SimG4FluxProducer") << "SimG4FluxProducer: =====> Begin event = " << iev;
0125   ++count_;
0126   store_.clear();
0127 }
0128
0129 void SimG4FluxProducer::update(const G4Step *aStep) {
0130   if (aStep != nullptr) {
0131     G4TouchableHistory *touchable = (G4TouchableHistory *)aStep->GetPreStepPoint()->GetTouchable();
0132     G4LogicalVolume *plv = (G4LogicalVolume *)touchable->GetVolume()->GetLogicalVolume();
0133     auto it = (init_) ? mapLV_.find(plv) : findLV(plv);
0134     //  edm::LogVerbatim("SimG4FluxProducer") << plv->GetName() << " Flag " << (it != mapLV_.end()) << " step " << aStep->IsFirstStepInVolume() << ":"  << aStep->IsLastStepInVolume();
0135     if (it != mapLV_.end()) {
0136       int type = LVTypes_[(it->second).first];
0137       if ((type == 0 && aStep->IsFirstStepInVolume()) || (type == 1 && aStep->IsLastStepInVolume())) {
0138         unsigned int copy = (unsigned int)(touchable->GetReplicaNumber(0));
0139         std::pair<G4LogicalVolume *, unsigned int> key(plv, copy);
0140         auto itr = store_.find(key);
0141         if (itr == store_.end()) {
0142           store_[key] = ParticleFlux((it->second).second, copy);
0143           itr = store_.find(key);
0144         }
0145         G4Track *track = aStep->GetTrack();
0146         int pdgid = track->GetDefinition()->GetPDGEncoding();
0147         int vxtyp = (track->GetCreatorProcess() == nullptr) ? 0 : 1;
0148         double time = (aStep->GetPostStepPoint()->GetGlobalTime());
0149         const double mmTocm(0.1), MeVToGeV(0.001);
0150         ParticleFlux::flux flx(pdgid, vxtyp, time);
0151         flx.vertex = math::GlobalPoint(mmTocm * track->GetVertexPosition().x(),
0152                                        mmTocm * track->GetVertexPosition().y(),
0153                                        mmTocm * track->GetVertexPosition().z());
0154         flx.hitPoint = math::GlobalPoint(
0155             mmTocm * track->GetPosition().x(), mmTocm * track->GetPosition().y(), mmTocm * track->GetPosition().z());
0156         flx.momentum = math::GlobalVector(MeVToGeV * track->GetMomentum().x(),
0157                                           MeVToGeV * track->GetMomentum().y(),
0158                                           MeVToGeV * track->GetMomentum().z());
0159         (itr->second).addFlux(flx);
0160 #ifdef EDM_ML_DEBUG
0161         edm::LogVerbatim("SimG4FluxProducer")
0162             << "SimG4FluxProducer: Element " << (it->second).first << ":" << (it->second).second << ":" << copy
0163             << " ID " << pdgid << " VxType " << vxtyp << " TOF " << time << " Hit Point " << flx.hitPoint << " p "
0164             << flx.momentum << " Vertex " << flx.vertex;
0165 #endif
0166       }  //if(Step ok)
0167     }    //if( it != map.end() )
0168   }      //if (aStep != NULL)
0169 }
0170
0171 void SimG4FluxProducer::endOfEvent(ParticleFlux &flux, unsigned int k) {
0172   bool done(false);
0173   for (const auto &element : store_) {
0174     G4LogicalVolume *lv = (element.first).first;
0175     auto it = mapLV_.find(lv);
0176     if (it != mapLV_.end()) {
0177       if ((it->second).first == k) {
0178         flux = element.second;
0179         done = true;
0180 #ifdef EDM_ML_DEBUG
0181         edm::LogVerbatim("SimG4FluxProducer") << "SimG4FluxProducer[" << k << "] Flux " << flux.getName() << ":"
0182                                               << flux.getId() << " with " << flux.getComponents() << " elements";
0183         std::vector<ParticleFlux::flux> fluxes = flux.getFlux();
0184         unsigned int k(0);
0185         for (auto element : fluxes) {
0186           edm::LogVerbatim("SimG4FluxProducer")
0187               << "Flux[" << k << "] PDGId " << element.pdgId << " VT " << element.vxType << " ToF " << element.tof
0188               << " Vertex " << element.vertex << " Hit " << element.hitPoint << " p " << element.momentum;
0189           ++k;
0190         }
0191 #endif
0192       }
0193     }
0194     if (done)
0195       break;
0196   }
0197 }
0198
0199 G4VPhysicalVolume *SimG4FluxProducer::getTopPV() {
0200   return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0201 }
0202
0203 std::map<G4LogicalVolume *, std::pair<unsigned int, std::string>>::iterator SimG4FluxProducer::findLV(
0204     G4LogicalVolume *plv) {
0205   auto itr = mapLV_.find(plv);
0206   if (itr == mapLV_.end()) {
0207     std::string name = plv->GetName();
0208     for (unsigned int k = 0; k < LVNames_.size(); ++k) {
0209       if (name.find(LVNames_[k]) != std::string::npos) {
0210         mapLV_[plv] = std::pair<unsigned int, std::string>(k, name);
0211         itr = mapLV_.find(plv);
0212         break;
0213       }
0214     }
0215   }
0216   return itr;
0217 }
0218 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0219 #include "FWCore/PluginManager/interface/ModuleDef.h"
0220
0221 DEFINE_SIMWATCHER(SimG4FluxProducer);