Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-28 23:36:05

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HGCPassive.cc
0003 // copied from SimG4HGCalValidation
0004 // Description: Main analysis class for HGCal Validation of G4 Hits
0005 ///////////////////////////////////////////////////////////////////////////////
0006 
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 
0012 // to retreive hits
0013 #include "SimDataFormats/CaloHit/interface/PassiveHit.h"
0014 
0015 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0016 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0017 #include "SimG4Core/Notification/interface/Observer.h"
0018 #include "SimG4Core/Watcher/interface/SimProducer.h"
0019 
0020 #include "G4LogicalVolumeStore.hh"
0021 #include "G4PhysicalVolumeStore.hh"
0022 #include "G4Step.hh"
0023 #include "G4TransportationManager.hh"
0024 #include "G4TouchableHistory.hh"
0025 #include "G4Track.hh"
0026 #include "DD4hep/Filter.h"
0027 
0028 #include <array>
0029 #include <cmath>
0030 #include <map>
0031 #include <string>
0032 #include <vector>
0033 
0034 //#define EDM_ML_DEBUG
0035 
0036 class HGCPassive : public SimProducer,
0037                    public Observer<const BeginOfRun *>,
0038                    public Observer<const BeginOfEvent *>,
0039                    public Observer<const G4Step *> {
0040 public:
0041   HGCPassive(const edm::ParameterSet &p);
0042   HGCPassive(const HGCPassive &) = delete;  // stop default
0043   const HGCPassive &operator=(const HGCPassive &) = delete;
0044   ~HGCPassive() override;
0045 
0046   void produce(edm::Event &, const edm::EventSetup &) override;
0047 
0048 private:
0049   // observer classes
0050   void update(const BeginOfRun *run) override;
0051   void update(const BeginOfEvent *evt) override;
0052   void update(const G4Step *step) override;
0053 
0054   // void endOfEvent(edm::PassiveHitContainer &HGCEEAbsE);
0055   void endOfEvent(edm::PassiveHitContainer &hgcPH, unsigned int k);
0056 
0057   typedef std::map<G4LogicalVolume *, std::pair<unsigned int, std::string>>::iterator volumeIterator;
0058   G4VPhysicalVolume *getTopPV();
0059   volumeIterator findLV(G4LogicalVolume *plv);
0060   void storeInfo(
0061       const volumeIterator itr, G4LogicalVolume *plv, unsigned int copy, double time, double energy, bool flag);
0062 
0063 private:
0064   const edm::ParameterSet m_Passive;
0065   const std::vector<std::string> LVNames_;
0066   const std::string motherName_;
0067   const int addlevel_;
0068   G4VPhysicalVolume *topPV_;
0069   G4LogicalVolume *topLV_;
0070   std::map<G4LogicalVolume *, std::pair<unsigned int, std::string>> mapLV_;
0071 
0072   // some private members for ananlysis
0073   unsigned int count_;
0074   bool init_;
0075   std::map<std::pair<G4LogicalVolume *, unsigned int>, std::array<double, 3>> store_;
0076 };
0077 
0078 HGCPassive::HGCPassive(const edm::ParameterSet &p)
0079     : m_Passive(p.getParameter<edm::ParameterSet>("HGCPassive")),
0080       LVNames_(m_Passive.getParameter<std::vector<std::string>>("LVNames")),
0081       motherName_(m_Passive.getParameter<std::string>("MotherName")),
0082       addlevel_((m_Passive.getParameter<bool>("IfDD4hep")) ? 1 : 0),
0083       topPV_(nullptr),
0084       topLV_(nullptr),
0085       count_(0),
0086       init_(false) {
0087 #ifdef EDM_ML_DEBUG
0088   edm::LogVerbatim("HGCSim") << "Name of the mother volume " << motherName_ << " AddLevel " << addlevel_;
0089   unsigned k(0);
0090 #endif
0091   for (const auto &name : LVNames_) {
0092     produces<edm::PassiveHitContainer>(Form("%sPassiveHits", name.c_str()));
0093 #ifdef EDM_ML_DEBUG
0094     edm::LogVerbatim("HGCSim") << "Collection name[" << k << "] " << name;
0095     ++k;
0096 #endif
0097   }
0098 }
0099 
0100 HGCPassive::~HGCPassive() {}
0101 
0102 void HGCPassive::produce(edm::Event &e, const edm::EventSetup &) {
0103   for (unsigned int k = 0; k < LVNames_.size(); ++k) {
0104     std::unique_ptr<edm::PassiveHitContainer> hgcPH(new edm::PassiveHitContainer);
0105     endOfEvent(*hgcPH, k);
0106     e.put(std::move(hgcPH), Form("%sPassiveHits", LVNames_[k].c_str()));
0107   }
0108 }
0109 
0110 void HGCPassive::update(const BeginOfRun *run) {
0111   topPV_ = getTopPV();
0112   if (topPV_ == nullptr) {
0113     edm::LogWarning("HGCSim") << "Cannot find top level volume\n";
0114   } else {
0115     init_ = true;
0116     const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
0117     for (auto lvcite : *lvs) {
0118       findLV(lvcite);
0119     }
0120 
0121 #ifdef EDM_ML_DEBUG
0122     edm::LogVerbatim("HGCSim") << "HGCPassive::Finds " << mapLV_.size() << " logical volumes";
0123     unsigned int k(0);
0124     for (const auto &lvs : mapLV_) {
0125       edm::LogVerbatim("HGCSim") << "Entry[" << k << "] " << lvs.first << ": (" << (lvs.second).first << ", "
0126                                  << (lvs.second).second << ")";
0127       ++k;
0128     }
0129 #endif
0130   }
0131 }
0132 
0133 //=================================================================== per EVENT
0134 void HGCPassive::update(const BeginOfEvent *evt) {
0135   int iev = (*evt)()->GetEventID();
0136   edm::LogVerbatim("HGCSim") << "HGCPassive: =====> Begin event = " << iev << std::endl;
0137 
0138   ++count_;
0139   store_.clear();
0140 }
0141 
0142 // //=================================================================== each
0143 // STEP
0144 void HGCPassive::update(const G4Step *aStep) {
0145   if (aStep != nullptr) {
0146     G4VSensitiveDetector *curSD = aStep->GetPreStepPoint()->GetSensitiveDetector();
0147     const G4VTouchable *touchable = aStep->GetPreStepPoint()->GetTouchable();
0148 
0149     int level = (touchable->GetHistoryDepth());
0150     if (curSD == nullptr) {
0151       G4LogicalVolume *plv = touchable->GetVolume()->GetLogicalVolume();
0152       auto it = (init_) ? mapLV_.find(plv) : findLV(plv);
0153       double time = aStep->GetTrack()->GetGlobalTime();
0154       double energy = (aStep->GetTotalEnergyDeposit()) / CLHEP::GeV;
0155 
0156       unsigned int copy(0);
0157       if (((aStep->GetPostStepPoint() == nullptr) || (aStep->GetTrack()->GetNextVolume() == nullptr)) &&
0158           (aStep->IsLastStepInVolume())) {
0159 #ifdef EDM_ML_DEBUG
0160         edm::LogVerbatim("HGCSim") << static_cast<std::string>(dd4hep::dd::noNamespace(plv->GetName())) << " F|L Step "
0161                                    << aStep->IsFirstStepInVolume() << ":" << aStep->IsLastStepInVolume() << " Position"
0162                                    << aStep->GetPreStepPoint()->GetPosition() << " Track "
0163                                    << aStep->GetTrack()->GetDefinition()->GetParticleName() << " at"
0164                                    << aStep->GetTrack()->GetPosition() << " Volume " << aStep->GetTrack()->GetVolume()
0165                                    << ":" << aStep->GetTrack()->GetNextVolume() << " Status "
0166                                    << aStep->GetTrack()->GetTrackStatus() << " KE "
0167                                    << aStep->GetTrack()->GetKineticEnergy() << " Deposit "
0168                                    << aStep->GetTotalEnergyDeposit() << " Map " << (it != mapLV_.end());
0169 #endif
0170         energy += (aStep->GetPreStepPoint()->GetKineticEnergy() / CLHEP::GeV);
0171       } else {
0172         time = (aStep->GetPostStepPoint()->GetGlobalTime());
0173         copy = (level < 2)
0174                    ? 0
0175                    : static_cast<unsigned int>(touchable->GetReplicaNumber(0) + 1000 * touchable->GetReplicaNumber(1));
0176       }
0177       if (it != mapLV_.end()) {
0178         storeInfo(it, plv, copy, time, energy, true);
0179       } else if (topLV_ != nullptr) {
0180         auto itr = findLV(topLV_);
0181         if (itr != mapLV_.end()) {
0182           storeInfo(itr, topLV_, copy, time, energy, true);
0183         }
0184       }
0185     }  // if (curSD==NULL)
0186 
0187     // Now for the mother volumes
0188     if (level > 0) {
0189       double energy = (aStep->GetTotalEnergyDeposit()) / CLHEP::GeV;
0190       double time = (aStep->GetTrack()->GetGlobalTime());
0191 
0192       for (int i = level; i > 0; --i) {
0193         G4LogicalVolume *plv = touchable->GetVolume(i)->GetLogicalVolume();
0194         auto it = (init_) ? mapLV_.find(plv) : findLV(plv);
0195 #ifdef EDM_ML_DEBUG
0196         edm::LogVerbatim("HGCSim") << "Level: " << level << ":" << i << " "
0197                                    << static_cast<std::string>(dd4hep::dd::noNamespace(plv->GetName()))
0198                                    << " flag in the List " << (it != mapLV_.end());
0199 #endif
0200         if (it != mapLV_.end()) {
0201           unsigned int copy =
0202               (i == level) ? 0
0203                            : (unsigned int)(touchable->GetReplicaNumber(i) + 1000 * touchable->GetReplicaNumber(i + 1));
0204           storeInfo(it, plv, copy, time, energy, false);
0205         }
0206       }
0207     }
0208   }  // if (aStep != NULL)
0209 
0210 }  // end update aStep
0211 
0212 //================================================================ End of EVENT
0213 
0214 void HGCPassive::endOfEvent(edm::PassiveHitContainer &hgcPH, unsigned int k) {
0215 #ifdef EDM_ML_DEBUG
0216   unsigned int kount(0);
0217 #endif
0218   for (const auto &element : store_) {
0219     G4LogicalVolume *lv = (element.first).first;
0220     auto it = mapLV_.find(lv);
0221     if (it != mapLV_.end()) {
0222       if ((it->second).first == k) {
0223         PassiveHit hit(
0224             (it->second).second, (element.first).second, (element.second)[1], (element.second)[2], (element.second)[0]);
0225         hgcPH.push_back(hit);
0226 #ifdef EDM_ML_DEBUG
0227         edm::LogVerbatim("HGCSim") << "HGCPassive[" << k << "] Hit[" << kount << "] " << hit;
0228         ++kount;
0229 #endif
0230       }
0231     }
0232   }
0233 }
0234 
0235 G4VPhysicalVolume *HGCPassive::getTopPV() {
0236   return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0237 }
0238 
0239 HGCPassive::volumeIterator HGCPassive::findLV(G4LogicalVolume *plv) {
0240   auto itr = mapLV_.find(plv);
0241   if (itr == mapLV_.end()) {
0242     std::string name = static_cast<std::string>(dd4hep::dd::noNamespace(plv->GetName()));
0243     for (unsigned int k = 0; k < LVNames_.size(); ++k) {
0244       if (name.find(LVNames_[k]) != std::string::npos) {
0245         mapLV_[plv] = std::pair<unsigned int, std::string>(k, name);
0246         itr = mapLV_.find(plv);
0247         break;
0248       }
0249     }
0250   }
0251   if (topLV_ == nullptr) {
0252     if (static_cast<std::string>(dd4hep::dd::noNamespace(plv->GetName())) == motherName_)
0253       topLV_ = plv;
0254   }
0255   return itr;
0256 }
0257 
0258 void HGCPassive::storeInfo(const HGCPassive::volumeIterator it,
0259                            G4LogicalVolume *plv,
0260                            unsigned int copy,
0261                            double time,
0262                            double energy,
0263                            bool flag) {
0264   std::pair<G4LogicalVolume *, unsigned int> key(plv, copy);
0265   auto itr = store_.find(key);
0266   double ee = (flag) ? energy : 0;
0267   if (itr == store_.end()) {
0268     store_[key] = {{time, energy, energy}};
0269   } else {
0270     (itr->second)[1] += ee;
0271     (itr->second)[2] += energy;
0272   }
0273 #ifdef EDM_ML_DEBUG
0274   itr = store_.find(key);
0275   edm::LogVerbatim("HGCSim") << "HGCPassive: Element " << (it->second).first << ":" << (it->second).second << ":"
0276                              << copy << " T " << (itr->second)[0] << " E " << (itr->second)[1] << ":"
0277                              << (itr->second)[2];
0278 #endif
0279 }
0280 
0281 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0282 #include "FWCore/PluginManager/interface/ModuleDef.h"
0283 
0284 DEFINE_SIMWATCHER(HGCPassive);