HGCalTestBeam/plugins/HGCalTBPassive.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HGCalTBPassive.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/Geometry/interface/DD4hep2DDDName.h"
0016 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0017 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0018 #include "SimG4Core/Notification/interface/Observer.h"
0019 #include "SimG4Core/Watcher/interface/SimProducer.h"
0020
0021 #include "G4LogicalVolumeStore.hh"
0022 #include "G4PhysicalVolumeStore.hh"
0023 #include "G4Step.hh"
0024 #include "G4TransportationManager.hh"
0025 #include "G4TouchableHistory.hh"
0026 #include "G4Track.hh"
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 HGCalTBPassive : public SimProducer,
0037                        public Observer<const BeginOfRun *>,
0038                        public Observer<const BeginOfEvent *>,
0039                        public Observer<const G4Step *> {
0040 public:
0041   HGCalTBPassive(const edm::ParameterSet &p);
0042   HGCalTBPassive(const HGCalTBPassive &) = delete;  // stop default
0043   const HGCalTBPassive &operator=(const HGCalTBPassive &) = delete;
0044   ~HGCalTBPassive() 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 HGCalTBPassive::HGCalTBPassive(const edm::ParameterSet &p)
0079     : m_Passive(p.getParameter<edm::ParameterSet>("HGCalTBPassive")),
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 HGCalTBPassive::~HGCalTBPassive() {}
0101
0102 void HGCalTBPassive::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 HGCalTBPassive::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") << "HGCalTBPassive::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 HGCalTBPassive::update(const BeginOfEvent *evt) {
0135   int iev = (*evt)()->GetEventID();
0136   edm::LogVerbatim("HGCSim") << "HGCalTBPassive: =====> Begin event = " << iev << std::endl;
0137
0138   ++count_;
0139   store_.clear();
0140 }
0141
0142 // //=================================================================== each
0143 // STEP
0144 void HGCalTBPassive::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") << DD4hep2DDDName::noNameSpace(static_cast<std::string>(plv->GetName()))
0161                                    << " F|L Step " << aStep->IsFirstStepInVolume() << ":" << aStep->IsLastStepInVolume()
0162                                    << " Position" << 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                                    << DD4hep2DDDName::noNameSpace(static_cast<std::string>(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 HGCalTBPassive::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") << "HGCalTBPassive[" << k << "] Hit[" << kount << "] " << hit;
0228         ++kount;
0229 #endif
0230       }
0231     }
0232   }
0233 }
0234
0235 G4VPhysicalVolume *HGCalTBPassive::getTopPV() {
0236   return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0237 }
0238
0239 HGCalTBPassive::volumeIterator HGCalTBPassive::findLV(G4LogicalVolume *plv) {
0240   auto itr = mapLV_.find(plv);
0241   if (itr == mapLV_.end()) {
0242     std::string name = DD4hep2DDDName::noNameSpace(static_cast<std::string>(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 (DD4hep2DDDName::noNameSpace(static_cast<std::string>(plv->GetName())) == motherName_)
0253       topLV_ = plv;
0254   }
0255   return itr;
0256 }
0257
0258 void HGCalTBPassive::storeInfo(const HGCalTBPassive::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") << "HGCalTBPassive: 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(HGCalTBPassive);