Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-18 08:23:52

0001 //#define EDM_ML_DEBUG
0002 //#define HcalNumberingTest
0003 
0004 // to make hits in EB/EE/HC
0005 
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0008 #include "DataFormats/Math/interface/Point3D.h"
0009 
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Utilities/interface/isFinite.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 
0017 #include "Geometry/EcalCommonData/interface/EcalBarrelNumberingScheme.h"
0018 #include "Geometry/EcalCommonData/interface/EcalBaseNumber.h"
0019 #include "Geometry/EcalCommonData/interface/EcalEndcapNumberingScheme.h"
0020 #ifdef HcalNumberingTest
0021 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
0022 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0023 #include "Geometry/Records/interface/HcalSimNumberingRecord.h"
0024 #endif
0025 
0026 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0027 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0028 #include "SimDataFormats/CaloHit/interface/PassiveHit.h"
0029 #include "SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h"
0030 
0031 #include "SimG4CMS/Calo/interface/CaloGVHit.h"
0032 #include "SimG4CMS/Calo/interface/HcalNumberingScheme.h"
0033 #include "SimG4CMS/Calo/interface/HcalNumberingFromPS.h"
0034 
0035 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0036 #include "SimG4Core/Notification/interface/Observer.h"
0037 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0038 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0039 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0040 #include "SimG4Core/Watcher/interface/SimProducer.h"
0041 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0042 
0043 #include "G4LogicalVolume.hh"
0044 #include "G4LogicalVolumeStore.hh"
0045 #include "G4NavigationHistory.hh"
0046 #include "G4ParticleTable.hh"
0047 #include "G4PhysicalVolumeStore.hh"
0048 #include "G4Region.hh"
0049 #include "G4RegionStore.hh"
0050 #include "G4Step.hh"
0051 #include "G4SystemOfUnits.hh"
0052 #include "G4Track.hh"
0053 #include "G4Trap.hh"
0054 #include "G4UnitsTable.hh"
0055 #include "G4UserSteppingAction.hh"
0056 #include "G4VPhysicalVolume.hh"
0057 #include "G4VTouchable.hh"
0058 
0059 #include <algorithm>
0060 #include <cmath>
0061 #include <iostream>
0062 #include <iomanip>
0063 #include <string>
0064 #include <utility>
0065 #include <vector>
0066 
0067 class CaloSteppingAction : public SimProducer,
0068                            public Observer<const BeginOfRun*>,
0069                            public Observer<const BeginOfEvent*>,
0070                            public Observer<const EndOfEvent*>,
0071                            public Observer<const G4Step*> {
0072 public:
0073   CaloSteppingAction(const edm::ParameterSet& p);
0074   ~CaloSteppingAction() override;
0075 
0076   void registerConsumes(edm::ConsumesCollector) override;
0077   void produce(edm::Event&, const edm::EventSetup&) override;
0078   void beginRun(edm::EventSetup const&) override;
0079 
0080 private:
0081   void fillHits(edm::PCaloHitContainer& cc, int type);
0082   void fillPassiveHits(edm::PassiveHitContainer& cc);
0083   // observer classes
0084   void update(const BeginOfRun* run) override;
0085   void update(const BeginOfEvent* evt) override;
0086   void update(const G4Step* step) override;
0087   void update(const EndOfEvent* evt) override;
0088 
0089   void NaNTrap(const G4Step*) const;
0090   uint32_t getDetIDHC(int det, int lay, int depth, const math::XYZVectorD& pos) const;
0091   void fillHit(uint32_t id, double dE, double time, int primID, uint16_t depth, double em, int flag);
0092   uint16_t getDepth(bool flag, double crystalDepth, double radl) const;
0093   double curve_LY(double crystalLength, double crystalDepth) const;
0094   double getBirkL3(double dE, double step, double chg, double dens) const;
0095   double getBirkHC(double dE, double step, double chg, double dens) const;
0096   void saveHits(int flag);
0097 
0098   static const int nSD_ = 3;
0099   std::unique_ptr<EcalBarrelNumberingScheme> ebNumberingScheme_;
0100   std::unique_ptr<EcalEndcapNumberingScheme> eeNumberingScheme_;
0101   std::unique_ptr<HcalNumberingFromPS> hcNumberingPS_;
0102 #ifdef HcalNumberingTest
0103   edm::ESGetToken<HcalDDDSimConstants, HcalSimNumberingRecord> ddconsToken_;
0104   std::unique_ptr<HcalNumberingFromDDD> hcNumbering_;
0105 #endif
0106   std::unique_ptr<HcalNumberingScheme> hcNumberingScheme_;
0107   std::unique_ptr<CaloSlaveSD> slave_[nSD_];
0108 
0109   const edm::ParameterSet iC_;
0110   const std::vector<std::string> nameEBSD_, nameEESD_, nameHCSD_;
0111   const std::vector<std::string> nameHitC_;
0112   std::vector<const G4LogicalVolume*> volEBSD_, volEESD_, volHCSD_;
0113   std::map<const G4LogicalVolume*, double> xtalMap_;
0114   std::map<const G4LogicalVolume*, std::string> mapLV_;
0115   const int allSteps_;
0116   const double slopeLY_, birkC1EC_, birkSlopeEC_;
0117   const double birkCutEC_, birkC1HC_, birkC2HC_;
0118   const double birkC3HC_, timeSliceUnit_;
0119   int count_, eventID_;
0120   std::map<std::pair<int, CaloHitID>, CaloGVHit> hitMap_[nSD_];
0121   typedef std::tuple<const G4LogicalVolume*, uint32_t, int, int, double, double, double, double, double, double, double>
0122       PassiveData;
0123   std::vector<PassiveData> store_;
0124 };
0125 
0126 CaloSteppingAction::CaloSteppingAction(const edm::ParameterSet& p)
0127     : iC_(p.getParameter<edm::ParameterSet>("CaloSteppingAction")),
0128       nameEBSD_(iC_.getParameter<std::vector<std::string> >("EBSDNames")),
0129       nameEESD_(iC_.getParameter<std::vector<std::string> >("EESDNames")),
0130       nameHCSD_(iC_.getParameter<std::vector<std::string> >("HCSDNames")),
0131       nameHitC_(iC_.getParameter<std::vector<std::string> >("HitCollNames")),
0132       allSteps_(iC_.getParameter<int>("AllSteps")),
0133       slopeLY_(iC_.getParameter<double>("SlopeLightYield")),
0134       birkC1EC_(iC_.getParameter<double>("BirkC1EC")),
0135       birkSlopeEC_(iC_.getParameter<double>("BirkSlopeEC")),
0136       birkCutEC_(iC_.getParameter<double>("BirkCutEC")),
0137       birkC1HC_(iC_.getParameter<double>("BirkC1HC")),
0138       birkC2HC_(iC_.getParameter<double>("BirkC2HC")),
0139       birkC3HC_(iC_.getParameter<double>("BirkC3HC")),
0140       timeSliceUnit_(iC_.getUntrackedParameter<double>("TimeSliceUnit", 1.0)),
0141       count_(0) {
0142   edm::ParameterSet iC = p.getParameter<edm::ParameterSet>("CaloSteppingAction");
0143 
0144   edm::LogVerbatim("Step") << "CaloSteppingAction:: " << nameEBSD_.size() << " names for EB SD's";
0145   for (unsigned int k = 0; k < nameEBSD_.size(); ++k)
0146     edm::LogVerbatim("Step") << "[" << k << "] " << nameEBSD_[k];
0147   edm::LogVerbatim("Step") << "CaloSteppingAction:: " << nameEESD_.size() << " names for EE SD's";
0148   for (unsigned int k = 0; k < nameEESD_.size(); ++k)
0149     edm::LogVerbatim("Step") << "[" << k << "] " << nameEESD_[k];
0150   edm::LogVerbatim("Step") << "CaloSteppingAction:: " << nameHCSD_.size() << " names for HC SD's";
0151   for (unsigned int k = 0; k < nameHCSD_.size(); ++k)
0152     edm::LogVerbatim("Step") << "[" << k << "] " << nameHCSD_[k];
0153   edm::LogVerbatim("Step") << "CaloSteppingAction::Constants for ECAL: slope " << slopeLY_ << " Birk constants "
0154                            << birkC1EC_ << ":" << birkSlopeEC_ << ":" << birkCutEC_;
0155   edm::LogVerbatim("Step") << "CaloSteppingAction::Constants for HCAL: Birk "
0156                            << "constants " << birkC1HC_ << ":" << birkC2HC_ << ":" << birkC3HC_;
0157   edm::LogVerbatim("Step") << "CaloSteppingAction::Constant for time slice " << timeSliceUnit_;
0158   edm::LogVerbatim("Step") << "CaloSteppingAction:: " << nameHitC_.size() << " hit collections";
0159   for (unsigned int k = 0; k < nameHitC_.size(); ++k)
0160     edm::LogVerbatim("Step") << "[" << k << "] " << nameHitC_[k];
0161 
0162   ebNumberingScheme_ = std::make_unique<EcalBarrelNumberingScheme>();
0163   eeNumberingScheme_ = std::make_unique<EcalEndcapNumberingScheme>();
0164   hcNumberingPS_ = std::make_unique<HcalNumberingFromPS>(iC);
0165   hcNumberingScheme_ = std::make_unique<HcalNumberingScheme>();
0166 #ifdef HcalNumberingTest
0167   hcNumbering_.reset(nullptr);
0168 #endif
0169   for (int k = 0; k < CaloSteppingAction::nSD_; ++k) {
0170     slave_[k] = std::make_unique<CaloSlaveSD>(nameHitC_[k]);
0171     produces<edm::PCaloHitContainer>(nameHitC_[k]);
0172   }
0173   if (allSteps_ != 0)
0174     produces<edm::PassiveHitContainer>("AllPassiveHits");
0175   edm::LogVerbatim("Step") << "CaloSteppingAction:: All Steps Flag " << allSteps_ << " for passive hits";
0176 }
0177 
0178 CaloSteppingAction::~CaloSteppingAction() {
0179   edm::LogVerbatim("Step") << "CaloSteppingAction: -------->  Total number of "
0180                            << "selected entries : " << count_;
0181 }
0182 
0183 void CaloSteppingAction::registerConsumes(edm::ConsumesCollector cc) {
0184 #ifdef HcalNumberingTest
0185   ddconsToken_ = cc.esConsumes<HcalDDDSimConstants, HcalSimNumberingRecord, edm::Transition::BeginRun>();
0186   edm::LogVerbatim("Step") << "CaloSteppingAction::Initialize ESGetToken for HcalDDDSimConstants";
0187 #endif
0188 }
0189 
0190 void CaloSteppingAction::produce(edm::Event& e, const edm::EventSetup&) {
0191   for (int k = 0; k < CaloSteppingAction::nSD_; ++k) {
0192     saveHits(k);
0193     auto product = std::make_unique<edm::PCaloHitContainer>();
0194     fillHits(*product, k);
0195     e.put(std::move(product), nameHitC_[k]);
0196   }
0197   if (allSteps_ != 0) {
0198     std::unique_ptr<edm::PassiveHitContainer> hgcPH(new edm::PassiveHitContainer);
0199     fillPassiveHits(*hgcPH);
0200     e.put(std::move(hgcPH), "AllPassiveHits");
0201   }
0202 }
0203 
0204 void CaloSteppingAction::fillHits(edm::PCaloHitContainer& cc, int type) {
0205   edm::LogVerbatim("Step") << "CaloSteppingAction::fillHits for type " << type << " with "
0206                            << slave_[type].get()->hits().size() << " hits";
0207   cc = slave_[type].get()->hits();
0208   slave_[type].get()->Clean();
0209 }
0210 
0211 void CaloSteppingAction::fillPassiveHits(edm::PassiveHitContainer& cc) {
0212   edm::LogVerbatim("Step") << "CaloSteppingAction::fillPassiveHits with " << store_.size() << " hits";
0213   for (const auto& element : store_) {
0214     auto lv = std::get<0>(element);
0215     auto it = mapLV_.find(lv);
0216     if (it != mapLV_.end()) {
0217       PassiveHit hit(it->second,
0218                      std::get<1>(element),
0219                      std::get<5>(element),
0220                      std::get<6>(element),
0221                      std::get<4>(element),
0222                      std::get<2>(element),
0223                      std::get<3>(element),
0224                      std::get<7>(element),
0225                      std::get<8>(element),
0226                      std::get<9>(element),
0227                      std::get<10>(element));
0228       cc.emplace_back(hit);
0229     }
0230   }
0231 }
0232 
0233 void CaloSteppingAction::beginRun(edm::EventSetup const& es) {
0234   edm::LogVerbatim("Step") << "CaloSteppingAction:: Enter BeginOfRun";
0235 
0236 #ifdef HcalNumberingTest
0237   // Numbering From DDD
0238   const HcalDDDSimConstants* hcons_ = &es.getData(ddconsToken_);
0239   edm::LogVerbatim("Step") << "CaloSteppingAction:: Initialise "
0240                            << "HcalNumberingFromDDD";
0241   hcNumbering_ = std::make_unique<HcalNumberingFromDDD>(hcons_);
0242 #endif
0243 }
0244 
0245 //==================================================================== per RUN
0246 void CaloSteppingAction::update(const BeginOfRun* run) {
0247   int irun = (*run)()->GetRunID();
0248   edm::LogVerbatim("Step") << "CaloSteppingAction:: Begin of Run = " << irun;
0249 
0250   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0251   if (lvs) {
0252     std::map<const std::string, const G4LogicalVolume*> nameMap;
0253     std::map<const std::string, const G4LogicalVolume*>::const_iterator itr;
0254     for (auto lvi = lvs->begin(), lve = lvs->end(); lvi != lve; ++lvi) {
0255       nameMap.emplace((*lvi)->GetName(), *lvi);
0256       if (allSteps_ < 0)
0257         mapLV_[*lvi] = (*lvi)->GetName();
0258     }
0259 
0260     for (auto const& name : nameEBSD_) {
0261       for (itr = nameMap.begin(); itr != nameMap.end(); ++itr) {
0262         const std::string& lvname = itr->first;
0263         if (lvname.find(name) != std::string::npos) {
0264           volEBSD_.emplace_back(itr->second);
0265           int type = (lvname.find("refl") == std::string::npos) ? -1 : 1;
0266           G4Trap* solid = static_cast<G4Trap*>(itr->second->GetSolid());
0267           double dz = 2 * solid->GetZHalfLength() / CLHEP::mm;
0268           xtalMap_.insert(std::pair<const G4LogicalVolume*, double>(itr->second, dz * type));
0269           if ((allSteps_ > 0) && ((allSteps_ % 10) > 0))
0270             mapLV_[itr->second] = itr->first;
0271         }
0272       }
0273     }
0274     for (auto const& name : nameEESD_) {
0275       for (itr = nameMap.begin(); itr != nameMap.end(); ++itr) {
0276         const std::string& lvname = itr->first;
0277         if (lvname.find(name) != std::string::npos) {
0278           volEESD_.emplace_back(itr->second);
0279           int type = (lvname.find("refl") == std::string::npos) ? 1 : -1;
0280           G4Trap* solid = static_cast<G4Trap*>(itr->second->GetSolid());
0281           double dz = 2 * solid->GetZHalfLength() / CLHEP::mm;
0282           xtalMap_.insert(std::pair<const G4LogicalVolume*, double>(itr->second, dz * type));
0283           if ((allSteps_ > 0) && (((allSteps_ / 10) % 10) > 0))
0284             mapLV_[itr->second] = itr->first;
0285         }
0286       }
0287     }
0288 
0289     for (auto const& name : nameHCSD_) {
0290       for (itr = nameMap.begin(); itr != nameMap.end(); ++itr) {
0291         const std::string& lvname = itr->first;
0292         if (lvname.find(name) != std::string::npos) {
0293           volHCSD_.emplace_back(itr->second);
0294           if ((allSteps_ > 0) && (((allSteps_ / 100) % 10) > 0))
0295             mapLV_[itr->second] = itr->first;
0296         }
0297       }
0298     }
0299   }
0300 #ifdef EDM_ML_DEBUG
0301   edm::LogVerbatim("Step") << volEBSD_.size() << " logical volumes for EB SD";
0302   for (unsigned int k = 0; k < volEBSD_.size(); ++k)
0303     edm::LogVerbatim("Step") << "[" << k << "] " << volEBSD_[k];
0304   edm::LogVerbatim("Step") << volEESD_.size() << " logical volumes for EE SD";
0305   for (unsigned int k = 0; k < volEESD_.size(); ++k)
0306     edm::LogVerbatim("Step") << "[" << k << "] " << volEESD_[k];
0307   edm::LogVerbatim("Step") << volHCSD_.size() << " logical volumes for HC SD";
0308   for (unsigned int k = 0; k < volHCSD_.size(); ++k)
0309     edm::LogVerbatim("Step") << "[" << k << "] " << volHCSD_[k];
0310   edm::LogVerbatim("Step") << mapLV_.size() << " logical volumes for Passive hits";
0311   unsigned int k(0);
0312   for (auto itr = mapLV_.begin(); itr != mapLV_.end(); ++itr) {
0313     edm::LogVerbatim("Step") << "[" << k << "] " << itr->second << ":" << itr->first;
0314     ++k;
0315   }
0316 #endif
0317 }
0318 
0319 //=================================================================== per EVENT
0320 void CaloSteppingAction::update(const BeginOfEvent* evt) {
0321   eventID_ = (*evt)()->GetEventID();
0322   edm::LogVerbatim("Step") << "CaloSteppingAction: Begin of event = " << eventID_;
0323   for (int k = 0; k < CaloSteppingAction::nSD_; ++k) {
0324     hitMap_[k].erase(hitMap_[k].begin(), hitMap_[k].end());
0325     slave_[k].get()->Initialize();
0326   }
0327   if (allSteps_ != 0)
0328     store_.clear();
0329 }
0330 
0331 //=================================================================== each STEP
0332 void CaloSteppingAction::update(const G4Step* aStep) {
0333   //  edm::LogVerbatim("Step") <<"CaloSteppingAction: At each Step";
0334   NaNTrap(aStep);
0335   auto lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0336   bool hc = (std::find(volHCSD_.begin(), volHCSD_.end(), lv) != volHCSD_.end());
0337   bool eb = (std::find(volEBSD_.begin(), volEBSD_.end(), lv) != volEBSD_.end());
0338   bool ee = (std::find(volEESD_.begin(), volEESD_.end(), lv) != volEESD_.end());
0339   uint32_t unitID(0);
0340   if (hc || eb || ee) {
0341     double dEStep = aStep->GetTotalEnergyDeposit() / CLHEP::MeV;
0342     auto const theTrack = aStep->GetTrack();
0343     double time = theTrack->GetGlobalTime() / CLHEP::nanosecond;
0344     int primID = theTrack->GetTrackID();
0345     bool em = G4TrackToParticleID::isGammaElectronPositron(theTrack);
0346     auto const touch = aStep->GetPreStepPoint()->GetTouchable();
0347     auto const& hitPoint = aStep->GetPreStepPoint()->GetPosition();
0348     if (hc) {
0349       int depth = (touch->GetReplicaNumber(0)) % 10 + 1;
0350       int lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
0351       int det = (touch->GetReplicaNumber(1)) / 1000;
0352       unitID = getDetIDHC(det, lay, depth, math::XYZVectorD(hitPoint.x(), hitPoint.y(), hitPoint.z()));
0353       if (unitID > 0 && dEStep > 0.0) {
0354         dEStep *= getBirkHC(dEStep,
0355                             (aStep->GetStepLength() / CLHEP::cm),
0356                             aStep->GetPreStepPoint()->GetCharge(),
0357                             (aStep->GetPreStepPoint()->GetMaterial()->GetDensity() / (CLHEP::g / CLHEP::cm3)));
0358         fillHit(unitID, dEStep, time, primID, 0, em, 2);
0359       }
0360     } else {
0361       EcalBaseNumber theBaseNumber;
0362       int size = touch->GetHistoryDepth() + 1;
0363       if (theBaseNumber.getCapacity() < size)
0364         theBaseNumber.setSize(size);
0365       //Get name and copy numbers
0366       if (size > 1) {
0367         for (int ii = 0; ii < size; ii++) {
0368           theBaseNumber.addLevel(touch->GetVolume(ii)->GetName(), touch->GetReplicaNumber(ii));
0369         }
0370       }
0371       unitID = (eb ? (ebNumberingScheme_->getUnitID(theBaseNumber)) : (eeNumberingScheme_->getUnitID(theBaseNumber)));
0372       if (unitID > 0 && dEStep > 0.0) {
0373         auto local = touch->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0374         auto ite = xtalMap_.find(lv);
0375         double crystalLength = ((ite == xtalMap_.end()) ? 230.0 : std::abs(ite->second));
0376         double crystalDepth =
0377             ((ite == xtalMap_.end()) ? 0.0 : (std::abs(0.5 * (ite->second) + (local.z() / CLHEP::mm))));
0378         double radl = aStep->GetPreStepPoint()->GetMaterial()->GetRadlen() / CLHEP::mm;
0379         bool flag = ((ite == xtalMap_.end()) ? true : (((ite->second) >= 0) ? true : false));
0380         auto depth = getDepth(flag, crystalDepth, radl);
0381         dEStep *= (getBirkL3(dEStep,
0382                              (aStep->GetStepLength() / CLHEP::cm),
0383                              aStep->GetPreStepPoint()->GetCharge(),
0384                              (aStep->GetPreStepPoint()->GetMaterial()->GetDensity() / (CLHEP::g / CLHEP::cm3))) *
0385                    curve_LY(crystalLength, crystalDepth));
0386         fillHit(unitID, dEStep, time, primID, depth, em, (eb ? 0 : 1));
0387       }
0388     }
0389   }
0390 
0391   if (allSteps_ != 0) {
0392     auto it = mapLV_.find(lv);
0393     if (it != mapLV_.end()) {
0394       double energy = aStep->GetTotalEnergyDeposit() / CLHEP::MeV;
0395       auto const touch = aStep->GetPreStepPoint()->GetTouchable();
0396       double time = aStep->GetTrack()->GetGlobalTime() / CLHEP::nanosecond;
0397       int trackId = aStep->GetTrack()->GetTrackID();
0398       int pdg = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0399       double stepl = (aStep->GetStepLength() / CLHEP::cm);
0400       double xp = aStep->GetPreStepPoint()->GetPosition().x() / CLHEP::cm;
0401       double yp = aStep->GetPreStepPoint()->GetPosition().y() / CLHEP::cm;
0402       double zp = aStep->GetPreStepPoint()->GetPosition().z() / CLHEP::cm;
0403 #ifdef EDM_ML_DEBUG
0404       edm::LogVerbatim("Step") << "CaloSteppingAction: Volume " << lv->GetName() << " History "
0405                                << touch->GetHistoryDepth() << " Pointers " << aStep->GetPostStepPoint() << ":"
0406                                << aStep->GetTrack()->GetNextVolume() << ":" << aStep->IsLastStepInVolume() << " E "
0407                                << energy << " T " << time << " PDG " << pdg << " step " << stepl << " Position (" << xp
0408                                << ", " << yp << ", " << zp << ")";
0409 #endif
0410       uint32_t copy = (allSteps_ < 0) ? 0 : unitID;
0411       if (((aStep->GetPostStepPoint() == nullptr) || (aStep->GetTrack()->GetNextVolume() == nullptr)) &&
0412           (aStep->IsLastStepInVolume())) {
0413         energy += (aStep->GetPreStepPoint()->GetKineticEnergy() / CLHEP::MeV);
0414       } else {
0415         time = aStep->GetPostStepPoint()->GetGlobalTime() / CLHEP::nanosecond;
0416         if (copy == 0)
0417           copy = (touch->GetHistoryDepth() < 1)
0418                      ? static_cast<uint32_t>(touch->GetReplicaNumber(0))
0419                      : static_cast<uint32_t>(touch->GetReplicaNumber(0) + 1000 * touch->GetReplicaNumber(1));
0420       }
0421       PassiveData key(std::make_tuple(lv, copy, trackId, pdg, time, energy, energy, stepl, xp, yp, zp));
0422       store_.push_back(key);
0423     }
0424   }
0425 }
0426 
0427 //================================================================ End of EVENT
0428 void CaloSteppingAction::update(const EndOfEvent* evt) {
0429   ++count_;
0430   // Fill event input
0431   edm::LogVerbatim("Step") << "CaloSteppingAction: EndOfEvent " << (*evt)()->GetEventID();
0432 }
0433 
0434 void CaloSteppingAction::NaNTrap(const G4Step* aStep) const {
0435   auto currentPos = aStep->GetTrack()->GetPosition();
0436   double xyz = currentPos.x() + currentPos.y() + currentPos.z();
0437   auto currentMom = aStep->GetTrack()->GetMomentum();
0438   xyz += currentMom.x() + currentMom.y() + currentMom.z();
0439 
0440   if (edm::isNotFinite(xyz)) {
0441     auto pCurrentVol = aStep->GetPreStepPoint()->GetPhysicalVolume();
0442     auto& nameOfVol = pCurrentVol->GetName();
0443     throw cms::Exception("Unknown", "CaloSteppingAction")
0444         << " Corrupted Event - NaN detected in volume " << nameOfVol << "\n";
0445   }
0446 }
0447 
0448 uint32_t CaloSteppingAction::getDetIDHC(int det, int lay, int depth, const math::XYZVectorD& pos) const {
0449   HcalNumberingFromDDD::HcalID tmp = hcNumberingPS_.get()->unitID(det, lay, depth, pos);
0450 #ifdef HcalNumberingTest
0451   auto id0 = hcNumberingScheme_.get()->getUnitID(tmp);
0452   HcalNumberingFromDDD::HcalID tmpO = hcNumbering_.get()->unitID(det, pos, depth, lay);
0453   auto idO = hcNumberingScheme_.get()->getUnitID(tmpO);
0454   std::string error = (id0 == idO) ? " ** OK **" : " ** ERROR **";
0455   edm::LogVerbatim("Step") << "Det ID " << HcalDetId(id0) << " Original " << HcalDetId(idO) << error;
0456 #endif
0457   return (hcNumberingScheme_.get()->getUnitID(tmp));
0458 }
0459 
0460 void CaloSteppingAction::fillHit(uint32_t id, double dE, double time, int primID, uint16_t depth, double em, int flag) {
0461   CaloHitID currentID(id, time, primID, depth, timeSliceUnit_);
0462   double edepEM = (em ? dE : 0);
0463   double edepHAD = (em ? 0 : dE);
0464   std::pair<int, CaloHitID> evID = std::make_pair(eventID_, currentID);
0465   auto it = hitMap_[flag].find(evID);
0466   if (it != hitMap_[flag].end()) {
0467     (it->second).addEnergyDeposit(edepEM, edepHAD);
0468   } else {
0469     CaloGVHit aHit;
0470     aHit.setEventID(eventID_);
0471     aHit.setID(currentID);
0472     aHit.addEnergyDeposit(edepEM, edepHAD);
0473     hitMap_[flag][evID] = aHit;
0474   }
0475 }
0476 
0477 uint16_t CaloSteppingAction::getDepth(bool flag, double crystalDepth, double radl) const {
0478   uint16_t depth1 = (flag ? 0 : PCaloHit::kEcalDepthRefz);
0479   uint16_t depth2 = (uint16_t)floor(crystalDepth / radl);
0480   uint16_t depth = (((depth2 & PCaloHit::kEcalDepthMask) << PCaloHit::kEcalDepthOffset) | depth1);
0481 #ifdef EDM_ML_DEBUG
0482   edm::LogVerbatim("Step") << "CaloSteppingAction::getDepth radl " << radl << ":" << crystalDepth << " depth " << depth;
0483 #endif
0484   return depth;
0485 }
0486 
0487 double CaloSteppingAction::curve_LY(double crystalLength, double crystalDepth) const {
0488   double weight = 1.;
0489   double dapd = crystalLength - crystalDepth;
0490   if (dapd >= -0.1 || dapd <= crystalLength + 0.1) {
0491     if (dapd <= 100.)
0492       weight = 1.0 + slopeLY_ - dapd * 0.01 * slopeLY_;
0493 #ifdef EDM_ML_DEBUG
0494     edm::LogVerbatim("Step") << "CaloSteppingAction::curve_LY " << crystalDepth << ":" << crystalLength << ":" << dapd
0495                              << ":" << weight;
0496 #endif
0497   } else {
0498     edm::LogWarning("Step") << "CaloSteppingAction: light coll curve : wrong "
0499                             << "distance to APD " << dapd << " crlength = " << crystalLength
0500                             << " crystal Depth = " << crystalDepth << " weight = " << weight;
0501   }
0502   return weight;
0503 }
0504 
0505 double CaloSteppingAction::getBirkL3(double dEStep, double step, double charge, double density) const {
0506   double weight = 1.;
0507   if (charge != 0. && step > 0.) {
0508     double dedx = dEStep / step;
0509     double rkb = birkC1EC_ / density;
0510     if (dedx > 0) {
0511       weight = 1. - birkSlopeEC_ * log(rkb * dedx);
0512       if (weight < birkCutEC_)
0513         weight = birkCutEC_;
0514       else if (weight > 1.)
0515         weight = 1.;
0516     }
0517 #ifdef EDM_ML_DEBUG
0518     edm::LogVerbatim("Step") << "CaloSteppingAction::getBirkL3 Charge " << charge << " dE/dx " << dedx << " Birk Const "
0519                              << rkb << " Weight = " << weight << " dE " << dEStep << " step " << step;
0520 #endif
0521   }
0522   return weight;
0523 }
0524 
0525 double CaloSteppingAction::getBirkHC(double dEStep, double step, double charge, double density) const {
0526   double weight = 1.;
0527   if (charge != 0. && step > 0.) {
0528     double dedx = dEStep / step;
0529     double rkb = birkC1HC_ / density;
0530     double c = birkC2HC_ * rkb * rkb;
0531     if (std::abs(charge) >= 2.)
0532       rkb /= birkC3HC_;
0533     weight = 1. / (1. + rkb * dedx + c * dedx * dedx);
0534 #ifdef EDM_ML_DEBUG
0535     edm::LogVerbatim("Step") << "CaloSteppingAction::getBirkHC Charge " << charge << " dE/dx " << dedx << " Birk Const "
0536                              << rkb << ", " << c << " Weight = " << weight << " dE " << dEStep;
0537 #endif
0538   }
0539   return weight;
0540 }
0541 
0542 void CaloSteppingAction::saveHits(int type) {
0543   edm::LogVerbatim("Step") << "CaloSteppingAction:: saveHits for type " << type << " with " << hitMap_[type].size()
0544                            << " hits";
0545   slave_[type].get()->ReserveMemory(hitMap_[type].size());
0546   for (auto const& hit : hitMap_[type]) {
0547     slave_[type].get()->processHits(hit.second.getUnitID(),
0548                                     0.001 * hit.second.getEM(),
0549                                     0.001 * hit.second.getHadr(),
0550                                     hit.second.getTimeSlice(),
0551                                     hit.second.getTrackID(),
0552                                     hit.second.getDepth());
0553   }
0554 }
0555 
0556 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0557 #include "FWCore/PluginManager/interface/ModuleDef.h"
0558 
0559 DEFINE_SIMWATCHER(CaloSteppingAction);