Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-12-06 02:04:14

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HCalSD.cc
0003 // Description: Sensitive Detector class for calorimeters
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "SimG4CMS/Calo/interface/HCalSD.h"
0007 #include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
0008 #include "SimG4CMS/Calo/interface/HcalDumpGeometry.h"
0009 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
0010 #include "SimG4Core/Notification/interface/TrackInformation.h"
0011 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0012 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0013 
0014 #include "CondFormats/GeometryObjects/interface/HcalSimulationParameters.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0018 
0019 #include "G4LogicalVolumeStore.hh"
0020 #include "G4LogicalVolume.hh"
0021 #include "G4Step.hh"
0022 #include "G4Track.hh"
0023 
0024 #include "G4SystemOfUnits.hh"
0025 #include "G4PhysicalConstants.hh"
0026 #include "Randomize.hh"
0027 
0028 #include "DD4hep/Filter.h"
0029 
0030 #include <iostream>
0031 #include <fstream>
0032 #include <iomanip>
0033 #include <sstream>
0034 
0035 //#define EDM_ML_DEBUG
0036 //#define plotDebug
0037 //#define printDebug
0038 
0039 #ifdef plotDebug
0040 #include <TH1F.h>
0041 #endif
0042 
0043 HCalSD::HCalSD(const std::string& name,
0044                const HcalDDDSimConstants* hcns,
0045                const HcalDDDRecConstants* hcnr,
0046                const HcalSimulationConstants* hscs,
0047                const HBHEDarkening* hbd,
0048                const HBHEDarkening* hed,
0049                const SensitiveDetectorCatalog& clg,
0050                edm::ParameterSet const& p,
0051                const SimTrackManager* manager)
0052     : CaloSD(name,
0053              clg,
0054              p,
0055              manager,
0056              (float)(p.getParameter<edm::ParameterSet>("HCalSD").getParameter<double>("TimeSliceUnit")),
0057              p.getParameter<edm::ParameterSet>("HCalSD").getParameter<bool>("IgnoreTrackID")),
0058       hcalConstants_(hcns),
0059       hcalSimConstants_(hscs),
0060       m_HBDarkening(hbd),
0061       m_HEDarkening(hed),
0062       isHF(false),
0063       weight_(1.0),
0064       depth_(1) {
0065   numberingFromDDD.reset(nullptr);
0066   numberingScheme.reset(nullptr);
0067   showerLibrary.reset(nullptr);
0068   hfshower.reset(nullptr);
0069   showerParam.reset(nullptr);
0070   showerPMT.reset(nullptr);
0071   showerBundle.reset(nullptr);
0072   m_HFDarkening.reset(nullptr);
0073   m_HcalTestNS.reset(nullptr);
0074 
0075   //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
0076   //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
0077   //static SimpleConfigurable<double> bk3(1.75,  "HCalSD:BirkC3");
0078   // Values from NIM 80 (1970) 239-244: as implemented in Geant3
0079 
0080   bool dd4hep = p.getParameter<bool>("g4GeometryDD4hepSource");
0081   edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HCalSD");
0082   useBirk = m_HC.getParameter<bool>("UseBirkLaw");
0083   double bunit = (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
0084   birk1 = m_HC.getParameter<double>("BirkC1") * bunit;
0085   birk2 = m_HC.getParameter<double>("BirkC2");
0086   birk3 = m_HC.getParameter<double>("BirkC3");
0087   useShowerLibrary = m_HC.getParameter<bool>("UseShowerLibrary");
0088   useParam = m_HC.getParameter<bool>("UseParametrize");
0089   testNumber = m_HC.getParameter<bool>("TestNumberingScheme");
0090   neutralDensity = m_HC.getParameter<bool>("doNeutralDensityFilter");
0091   usePMTHit = m_HC.getParameter<bool>("UsePMTHits");
0092   betaThr = m_HC.getParameter<double>("BetaThreshold");
0093   eminHitHB = m_HC.getParameter<double>("EminHitHB") * MeV;
0094   eminHitHE = m_HC.getParameter<double>("EminHitHE") * MeV;
0095   eminHitHO = m_HC.getParameter<double>("EminHitHO") * MeV;
0096   eminHitHF = m_HC.getParameter<double>("EminHitHF") * MeV;
0097   useFibreBundle = m_HC.getParameter<bool>("UseFibreBundleHits");
0098   deliveredLumi = m_HC.getParameter<double>("DelivLuminosity");
0099   agingFlagHB = m_HC.getParameter<bool>("HBDarkening");
0100   agingFlagHE = m_HC.getParameter<bool>("HEDarkening");
0101   bool agingFlagHF = m_HC.getParameter<bool>("HFDarkening");
0102   useHF = m_HC.getUntrackedParameter<bool>("UseHF", true);
0103   bool forTBHC = m_HC.getUntrackedParameter<bool>("ForTBHCAL", false);
0104   bool forTBH2 = m_HC.getUntrackedParameter<bool>("ForTBH2", false);
0105   useLayerWt = m_HC.getUntrackedParameter<bool>("UseLayerWt", false);
0106   std::string file = m_HC.getUntrackedParameter<std::string>("WtFile", "None");
0107   testNS_ = m_HC.getUntrackedParameter<bool>("TestNS", false);
0108   edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
0109   applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
0110   bool dumpGeom = m_HC.getUntrackedParameter<bool>("DumpGeometry", false);
0111 
0112 #ifdef EDM_ML_DEBUG
0113   edm::LogVerbatim("HcalSim") << "***************************************************"
0114                               << "\n"
0115                               << "* Constructing a HCalSD  with name " << name << "\n"
0116                               << "\n"
0117                               << "***************************************************";
0118 #endif
0119   edm::LogVerbatim("HcalSim") << "HCalSD:: Use of HF code is set to " << useHF
0120                               << "\n         Use of shower parametrization set to " << useParam
0121                               << "\n         Use of shower library is set to " << useShowerLibrary
0122                               << "\n         Use PMT Hit is set to " << usePMTHit << " with beta Threshold " << betaThr
0123                               << "\n         USe of FibreBundle Hit set to " << useFibreBundle
0124                               << "\n         Use of Birks law is set to " << useBirk
0125                               << " with three constants kB = " << birk1 / bunit << ", C1 = " << birk2
0126                               << ", C2 = " << birk3;
0127   edm::LogVerbatim("HcalSim") << "HCalSD:: Suppression Flag " << suppressHeavy << " protons below " << kmaxProton
0128                               << " MeV,"
0129                               << " neutrons below " << kmaxNeutron << " MeV and"
0130                               << " ions below " << kmaxIon << " MeV\n"
0131                               << "         Threshold for storing hits in HB: " << eminHitHB << " HE: " << eminHitHE
0132                               << " HO: " << eminHitHO << " HF: " << eminHitHF << "\n"
0133                               << "Delivered luminosity for Darkening " << deliveredLumi << " Flag (HE) " << agingFlagHE
0134                               << " Flag (HB) " << agingFlagHB << " Flag (HF) " << agingFlagHF << "\n"
0135                               << "Application of Fiducial Cut " << applyFidCut
0136                               << "Flag for test number|neutral density filter " << testNumber << " " << neutralDensity;
0137 
0138   if (forTBHC) {
0139     useHF = false;
0140     matNames.emplace_back("Scintillator");
0141   } else {
0142     matNames = hcalSimConstants_->hcalsimpar()->hcalMaterialNames_;
0143   }
0144 
0145   HcalNumberingScheme* scheme;
0146   if (testNumber || forTBH2) {
0147     scheme = dynamic_cast<HcalNumberingScheme*>(new HcalTestNumberingScheme(forTBH2));
0148   } else {
0149     scheme = new HcalNumberingScheme();
0150   }
0151   setNumberingScheme(scheme);
0152 
0153   // always call getFromLibrary() method to identify HF region
0154   setParameterized(true);
0155 
0156   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0157   //  std::vector<G4LogicalVolume *>::const_iterator lvcite;
0158   const G4LogicalVolume* lv;
0159   std::string attribute, value;
0160 
0161   if (useHF) {
0162     if (useParam) {
0163       showerParam = std::make_unique<HFShowerParam>(name, hcalConstants_, hcalSimConstants_->hcalsimpar(), p);
0164     } else {
0165       if (useShowerLibrary) {
0166         showerLibrary = std::make_unique<HFShowerLibrary>(name, hcalConstants_, hcalSimConstants_->hcalsimpar(), p);
0167       }
0168       hfshower = std::make_unique<HFShower>(name, hcalConstants_, hcalSimConstants_->hcalsimpar(), p, 0);
0169     }
0170 
0171     // HF volume names
0172     hfNames = hcalSimConstants_->hcalsimpar()->hfNames_;
0173     const std::vector<int>& temp = hcalSimConstants_->hcalsimpar()->hfLevels_;
0174 #ifdef EDM_ML_DEBUG
0175     std::stringstream ss0;
0176     ss0 << "HCalSD: Names to be tested for Volume = HF has " << hfNames.size() << " elements";
0177 #endif
0178     int addlevel = dd4hep ? 1 : 0;
0179     for (unsigned int i = 0; i < hfNames.size(); ++i) {
0180       G4String namv(static_cast<std::string>(dd4hep::dd::noNamespace(hfNames[i])));
0181       lv = nullptr;
0182       for (auto lvol : *lvs) {
0183         if (dd4hep::dd::noNamespace(lvol->GetName()) == namv) {
0184           lv = lvol;
0185           break;
0186         }
0187       }
0188       hfLV.emplace_back(lv);
0189       hfLevels.emplace_back(temp[i] + addlevel);
0190 #ifdef EDM_ML_DEBUG
0191       ss0 << "\n        HF[" << i << "] = " << namv << " LV " << lv << " at level " << (temp[i] + addlevel);
0192 #endif
0193     }
0194 #ifdef EDM_ML_DEBUG
0195     edm::LogVerbatim("HcalSim") << ss0.str();
0196 #endif
0197     // HF Fibre volume names
0198     fibreNames = hcalSimConstants_->hcalsimpar()->hfFibreNames_;
0199     fillLogVolumeVector("HFFibre", fibreNames, fibreLV);
0200     const std::vector<std::string>& pmtNames = hcalSimConstants_->hcalsimpar()->hfPMTNames_;
0201     fillLogVolumeVector("HFPMT", pmtNames, pmtLV);
0202     const std::vector<std::string>& straightNames = hcalSimConstants_->hcalsimpar()->hfFibreStraightNames_;
0203     fillLogVolumeVector("HFFibreBundleStraight", straightNames, fibre1LV);
0204     const std::vector<std::string>& conicalNames = hcalSimConstants_->hcalsimpar()->hfFibreConicalNames_;
0205     fillLogVolumeVector("HFFibreBundleConical", conicalNames, fibre2LV);
0206   }
0207 
0208   //Material list for HB/HE/HO sensitive detectors
0209   const G4MaterialTable* matTab = G4Material::GetMaterialTable();
0210   std::vector<G4Material*>::const_iterator matite;
0211   for (auto const& namx : matNames) {
0212     const G4Material* mat = nullptr;
0213     for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
0214       if (static_cast<std::string>(dd4hep::dd::noNamespace((*matite)->GetName())) == namx) {
0215         mat = (*matite);
0216         break;
0217       }
0218     }
0219     materials.emplace_back(mat);
0220   }
0221 #ifdef EDM_ML_DEBUG
0222   std::stringstream ss1;
0223   for (unsigned int i = 0; i < matNames.size(); ++i) {
0224     if (i / 10 * 10 == i) {
0225       ss1 << "\n";
0226     }
0227     ss1 << "  " << matNames[i];
0228   }
0229   edm::LogVerbatim("HcalSim") << "HCalSD: Material names for HCAL: " << ss1.str();
0230 #endif
0231   if (useLayerWt) {
0232     readWeightFromFile(file);
0233   }
0234   numberingFromDDD = std::make_unique<HcalNumberingFromDDD>(hcalConstants_);
0235 
0236   //Special Geometry parameters
0237   gpar = hcalConstants_->getGparHF();
0238 #ifdef EDM_ML_DEBUG
0239   std::stringstream ss2;
0240   for (unsigned int ig = 0; ig < gpar.size(); ig++) {
0241     ss2 << "\n         gpar[" << ig << "] = " << gpar[ig] / cm << " cm";
0242   }
0243   edm::LogVerbatim("HcalSim") << "Maximum depth for HF " << hcalConstants_->getMaxDepth(2) << gpar.size()
0244                               << " gpar (cm)" << ss2.str();
0245 #endif
0246 
0247   //Test Hcal Numbering Scheme
0248   if (testNS_)
0249     m_HcalTestNS = std::make_unique<HcalTestNS>(hcnr);
0250 
0251   for (int i = 0; i < 9; ++i) {
0252     hit_[i] = time_[i] = dist_[i] = nullptr;
0253   }
0254   hzvem = hzvhad = nullptr;
0255 
0256   if (agingFlagHF) {
0257     m_HFDarkening = std::make_unique<HFDarkening>(m_HC.getParameter<edm::ParameterSet>("HFDarkeningParameterBlock"));
0258   }
0259 #ifdef plotDebug
0260   edm::Service<TFileService> tfile;
0261 
0262   if (tfile.isAvailable()) {
0263     static const char* const labels[] = {"HB",
0264                                          "HE",
0265                                          "HO",
0266                                          "HF Absorber",
0267                                          "HF PMT",
0268                                          "HF Absorber Long",
0269                                          "HF Absorber Short",
0270                                          "HF PMT Long",
0271                                          "HF PMT Short"};
0272     TFileDirectory hcDir = tfile->mkdir("ProfileFromHCalSD");
0273     char name[20], title[60];
0274     for (int i = 0; i < 9; ++i) {
0275       sprintf(title, "Hit energy in %s", labels[i]);
0276       sprintf(name, "HCalSDHit%d", i);
0277       hit_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
0278       sprintf(title, "Energy (MeV)");
0279       hit_[i]->GetXaxis()->SetTitle(title);
0280       hit_[i]->GetYaxis()->SetTitle("Hits");
0281       sprintf(title, "Time of the hit in %s", labels[i]);
0282       sprintf(name, "HCalSDTime%d", i);
0283       time_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
0284       sprintf(title, "Time (ns)");
0285       time_[i]->GetXaxis()->SetTitle(title);
0286       time_[i]->GetYaxis()->SetTitle("Hits");
0287       sprintf(title, "Longitudinal profile in %s", labels[i]);
0288       sprintf(name, "HCalSDDist%d", i);
0289       dist_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
0290       sprintf(title, "Distance (mm)");
0291       dist_[i]->GetXaxis()->SetTitle(title);
0292       dist_[i]->GetYaxis()->SetTitle("Hits");
0293     }
0294     if (useHF && (!useParam)) {
0295       hzvem = hcDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part)", 330, 0.0, 1650.0);
0296       hzvem->GetXaxis()->SetTitle("Longitudinal Profile (EM Part)");
0297       hzvhad = hcDir.make<TH1F>("hzvhad", "Longitudinal Profile (Had Part)", 330, 0.0, 1650.0);
0298       hzvhad->GetXaxis()->SetTitle("Longitudinal Profile (Hadronic Part)");
0299     }
0300   }
0301 #endif
0302   if (dumpGeom) {
0303     std::unique_ptr<HcalNumberingFromDDD> hcn = std::make_unique<HcalNumberingFromDDD>(hcalConstants_);
0304     const auto& lvNames = clg.logicalNames(name);
0305     HcalDumpGeometry geom(lvNames, hcn.get(), testNumber, false);
0306     geom.update();
0307   }
0308 }
0309 
0310 void HCalSD::fillLogVolumeVector(const std::string& value,
0311                                  const std::vector<std::string>& lvnames,
0312                                  std::vector<const G4LogicalVolume*>& lvvec) {
0313   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0314   const G4LogicalVolume* lv;
0315   std::stringstream ss3;
0316   ss3 << "HCalSD: " << lvnames.size() << " names to be tested for Volume <" << value << ">:";
0317   for (unsigned int i = 0; i < lvnames.size(); ++i) {
0318     G4String namv(static_cast<std::string>(dd4hep::dd::noNamespace(lvnames[i])));
0319     lv = nullptr;
0320     for (auto lvol : *lvs) {
0321       if (dd4hep::dd::noNamespace(lvol->GetName()) == namv) {
0322         lv = lvol;
0323         break;
0324       }
0325     }
0326     lvvec.emplace_back(lv);
0327     if (i / 10 * 10 == i) {
0328       ss3 << "\n";
0329     }
0330     ss3 << "  " << namv;
0331   }
0332   edm::LogVerbatim("HcalSim") << ss3.str();
0333 }
0334 
0335 bool HCalSD::getFromLibrary(const G4Step* aStep) {
0336   auto const track = aStep->GetTrack();
0337   depth_ = (aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(0)) % 10;
0338   weight_ = 1.0;
0339   bool kill(false);
0340   isHF = isItHF(aStep);
0341 #ifdef EDM_ML_DEBUG
0342   edm::LogVerbatim("HcalSim") << "GetFromLibrary: isHF " << isHF << " darken " << (m_HFDarkening != nullptr)
0343                               << " useParam " << useParam << " useShowerLibrary " << useShowerLibrary << " Muon? "
0344                               << G4TrackToParticleID::isMuon(track) << " electron? "
0345                               << G4TrackToParticleID::isGammaElectronPositron(track) << " Stable Hadron? "
0346                               << G4TrackToParticleID::isStableHadronIon(track);
0347 #endif
0348   if (isHF) {
0349     if (m_HFDarkening) {
0350       G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0351       const double invcm = 1. / CLHEP::cm;
0352       double r = hitPoint.perp() * invcm;
0353       double z = std::abs(hitPoint.z()) * invcm;
0354       double dose_acquired = 0.;
0355       if (z >= HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
0356         unsigned int hfZLayer = (unsigned int)((z - HFDarkening::lowZLimit) / 5);
0357         if (hfZLayer >= HFDarkening::upperZLimit)
0358           hfZLayer = (HFDarkening::upperZLimit - 1);
0359         float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
0360         for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
0361           dose_acquired = m_HFDarkening->dose(i, r);
0362           weight_ *= m_HFDarkening->degradation(normalized_lumi * dose_acquired);
0363         }
0364       }
0365 #ifdef EDM_ML_DEBUG
0366       edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary: HFLumiDarkening at "
0367                                   << "r= " << r << ", z= " << z << " Dose= " << dose_acquired << " weight= " << weight_;
0368 #endif
0369     }
0370 
0371     if (useParam) {
0372       getFromParam(aStep, kill);
0373 #ifdef EDM_ML_DEBUG
0374       G4String nameVolume = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
0375       edm::LogVerbatim("HcalSim") << "HCalSD: " << getNumberOfHits() << " hits from parametrization in " << nameVolume
0376                                   << " for Track " << track->GetTrackID() << " ("
0377                                   << track->GetDefinition()->GetParticleName() << ")";
0378 #endif
0379     } else if (useShowerLibrary && !G4TrackToParticleID::isMuon(track)) {
0380       if (G4TrackToParticleID::isGammaElectronPositron(track) || G4TrackToParticleID::isStableHadronIon(track)) {
0381 #ifdef EDM_ML_DEBUG
0382         auto nameVolume = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
0383         edm::LogVerbatim("HcalSim") << "HCalSD: Starts shower library from " << nameVolume << " for Track "
0384                                     << track->GetTrackID() << " (" << track->GetDefinition()->GetParticleName() << ")";
0385 
0386 #endif
0387         getFromHFLibrary(aStep, kill);
0388       }
0389     }
0390   }
0391 #ifdef EDM_ML_DEBUG
0392   edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary ID= " << track->GetTrackID() << " ("
0393                               << track->GetDefinition()->GetParticleName() << ") kill= " << kill
0394                               << " weight= " << weight_ << " depth= " << depth_ << " isHF: " << isHF;
0395 #endif
0396   return kill;
0397 }
0398 
0399 double HCalSD::getEnergyDeposit(const G4Step* aStep) {
0400   double destep(0.0);
0401   auto const lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0402   auto const theTrack = aStep->GetTrack();
0403 
0404   if (isHF) {
0405     if (useShowerLibrary && G4TrackToParticleID::isMuon(theTrack) && isItFibre(lv)) {
0406 #ifdef EDM_ML_DEBUG
0407       edm::LogVerbatim("HcalSim") << "HCalSD: Hit at Fibre in LV " << lv->GetName() << " for track "
0408                                   << aStep->GetTrack()->GetTrackID() << " ("
0409                                   << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0410 #endif
0411       hitForFibre(aStep);
0412     }
0413     return destep;
0414   }
0415 
0416   if (isItPMT(lv)) {
0417     if (usePMTHit && showerPMT) {
0418       getHitPMT(aStep);
0419     }
0420 #ifdef EDM_ML_DEBUG
0421     edm::LogVerbatim("HcalSim") << "HCalSD: Hit from PMT parametrization in LV " << lv->GetName() << " for Track "
0422                                 << aStep->GetTrack()->GetTrackID() << " ("
0423                                 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0424 #endif
0425     return destep;
0426 
0427   } else if (isItStraightBundle(lv)) {
0428     if (useFibreBundle && showerBundle) {
0429       getHitFibreBundle(aStep, false);
0430     }
0431 #ifdef EDM_ML_DEBUG
0432     edm::LogVerbatim("HcalSim") << "HCalSD: Hit from straight FibreBundle in LV: " << lv->GetName() << " for track "
0433                                 << aStep->GetTrack()->GetTrackID() << " ("
0434                                 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0435 #endif
0436     return destep;
0437 
0438   } else if (isItConicalBundle(lv)) {
0439     if (useFibreBundle && showerBundle) {
0440       getHitFibreBundle(aStep, true);
0441     }
0442 #ifdef EDM_ML_DEBUG
0443     edm::LogVerbatim("HcalSim") << "HCalSD: Hit from conical FibreBundle PV: " << lv->GetName() << " for track "
0444                                 << aStep->GetTrack()->GetTrackID() << " ("
0445                                 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0446 #endif
0447     return destep;
0448   }
0449 
0450   // normal hit
0451   destep = aStep->GetTotalEnergyDeposit();
0452 
0453   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0454   uint32_t detid = setDetUnitId(aStep);
0455   int det(0), ieta(0), phi(0), z(0), lay, depth(-1);
0456   if (testNumber) {
0457     HcalTestNumbering::unpackHcalIndex(detid, det, z, depth, ieta, phi, lay);
0458     if (z == 0) {
0459       z = -1;
0460     }
0461   } else {
0462     HcalDetId hcid(detid);
0463     det = hcid.subdetId();
0464     ieta = hcid.ietaAbs();
0465     phi = hcid.iphi();
0466     z = hcid.zside();
0467   }
0468   lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
0469 #ifdef EDM_ML_DEBUG
0470   edm::LogVerbatim("HcalSim") << "HCalSD: det: " << det << " ieta: " << ieta << " iphi: " << phi << " zside " << z
0471                               << "  lay: " << lay - 2;
0472 #endif
0473   if (depth_ == 0 && (det == 1 || det == 2) && ((!testNumber) || neutralDensity))
0474     weight_ = hcalConstants_->getLayer0Wt(det, phi, z);
0475   if (useLayerWt) {
0476     G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0477     weight_ = layerWeight(det + 2, hitPoint, depth_, lay);
0478   }
0479 
0480   if (agingFlagHB && m_HBDarkening && det == 1) {
0481     double dweight = m_HBDarkening->degradation(deliveredLumi, ieta, lay);
0482     weight_ *= dweight;
0483 #ifdef EDM_ML_DEBUG
0484     edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi << " coefficient = " << dweight
0485                                 << " Weight= " << weight_;
0486 #endif
0487   }
0488 
0489   if (agingFlagHE && m_HEDarkening && det == 2) {
0490     double dweight = m_HEDarkening->degradation(deliveredLumi, ieta, lay);
0491     weight_ *= dweight;
0492 #ifdef EDM_ML_DEBUG
0493     edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi << " coefficient = " << dweight
0494                                 << " Weight= " << weight_;
0495 #endif
0496   }
0497 
0498   if (suppressHeavy) {
0499     TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0500     if (trkInfo) {
0501       int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0502       if (!(trkInfo->isPrimary())) {  // Only secondary particles
0503         double ke = theTrack->GetKineticEnergy();
0504         if (pdg / 1000000000 == 1 && (pdg / 10000) % 100 > 0 && (pdg / 10) % 100 > 0 && ke < kmaxIon)
0505           weight_ = 0;
0506         if ((pdg == 2212) && (ke < kmaxProton))
0507           weight_ = 0;
0508         if ((pdg == 2112) && (ke < kmaxNeutron))
0509           weight_ = 0;
0510       }
0511     }
0512   }
0513   double wt0(1.0);
0514   if (useBirk) {
0515     const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0516     if (isItScintillator(mat))
0517       wt0 = getAttenuation(aStep, birk1, birk2, birk3);
0518   }
0519   weight_ *= wt0;
0520   double wt1 = getResponseWt(theTrack);
0521   double wt2 = theTrack->GetWeight();
0522   double edep = weight_ * wt1 * destep;
0523   if (wt2 > 0.0) {
0524     edep *= wt2;
0525   }
0526 #ifdef EDM_ML_DEBUG
0527   edm::LogVerbatim("HcalSim") << "HCalSD: edep= " << edep << " Det: " << det + 2 << " depth= " << depth_
0528                               << " weight= " << weight_ << " wt0= " << wt0 << " wt1= " << wt1 << " wt2= " << wt2;
0529 #endif
0530   return edep;
0531 }
0532 
0533 uint32_t HCalSD::setDetUnitId(const G4Step* aStep) {
0534   auto const prePoint = aStep->GetPreStepPoint();
0535   auto const touch = prePoint->GetTouchable();
0536   const G4ThreeVector& hitPoint = prePoint->GetPosition();
0537 
0538   int depth = (touch->GetReplicaNumber(0)) % 10 + 1;
0539   int lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
0540   int det = (touch->GetReplicaNumber(1)) / 1000;
0541 
0542   return setDetUnitId(det, hitPoint, depth, lay);
0543 }
0544 
0545 void HCalSD::setNumberingScheme(HcalNumberingScheme* scheme) {
0546   if (scheme != nullptr) {
0547     edm::LogVerbatim("HcalSim") << "HCalSD: updates numbering scheme for " << GetName();
0548     numberingScheme.reset(scheme);
0549   }
0550 }
0551 
0552 void HCalSD::update(const BeginOfJob* job) {}
0553 
0554 void HCalSD::initRun() {}
0555 
0556 bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
0557   double threshold = 0;
0558   DetId theId(aHit->getUnitID());
0559   switch (theId.subdetId()) {
0560     case HcalBarrel:
0561       threshold = eminHitHB;
0562       break;
0563     case HcalEndcap:
0564       threshold = eminHitHE;
0565       break;
0566     case HcalOuter:
0567       threshold = eminHitHO;
0568       break;
0569     case HcalForward:
0570       threshold = eminHitHF;
0571       break;
0572     default:
0573       break;
0574   }
0575   return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
0576 }
0577 
0578 uint32_t HCalSD::setDetUnitId(int det, const G4ThreeVector& pos, int depth, int lay = 1) {
0579   uint32_t id = 0;
0580   if (det == 0) {
0581 #ifdef printDebug
0582     double eta = std::abs(pos.eta());
0583 #endif
0584     if (std::abs(pos.z()) > maxZ_) {
0585       det = 5;
0586 #ifdef printDebug
0587       if (eta < 2.868)
0588         ++detNull_[2];
0589 #endif
0590     } else if (!(hcalConstants_->isHE())) {
0591       det = 3;
0592 #ifdef printDebug
0593       ++detNull_[0];
0594 #endif
0595     } else {
0596       double minR = minRoff_ + slopeHE_ * std::abs(pos.z());
0597       double maxR = maxRoff_ + slopeHE_ * std::abs(pos.z());
0598       det = ((pos.perp() > minR) && (pos.perp() < maxR)) ? 4 : 3;
0599 #ifdef printDebug
0600       ++detNull_[det - 3];
0601 #endif
0602     }
0603 #ifdef printDEBUG
0604     edm::LogVerbatim("HcalSim") << "Position " << pos.perp() << ":" << std::abs(pos.z()) << " Limits "
0605                                 << !(hcalConstants_->isHE()) << ":" << maxZ_ << " det " << det;
0606   } else {
0607     ++detNull_[3];
0608 #endif
0609   }
0610 
0611   if (numberingFromDDD.get()) {
0612     //get the ID's as eta, phi, depth, ... indices
0613     HcalNumberingFromDDD::HcalID tmp =
0614         numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
0615     id = setDetUnitId(tmp);
0616   }
0617   return id;
0618 }
0619 
0620 uint32_t HCalSD::setDetUnitId(HcalNumberingFromDDD::HcalID& tmp) {
0621   modifyDepth(tmp);
0622   uint32_t id = (numberingScheme.get()) ? numberingScheme->getUnitID(tmp) : 0;
0623   if ((!testNumber) && m_HcalTestNS.get()) {
0624     bool ok = m_HcalTestNS->compare(tmp, id);
0625     if (!ok)
0626       edm::LogWarning("HcalSim") << "Det ID from HCalSD " << HcalDetId(id) << " " << std::hex << id << std::dec
0627                                  << " does not match one from relabller for " << tmp.subdet << ":" << tmp.etaR << ":"
0628                                  << tmp.phi << ":" << tmp.phis << ":" << tmp.depth << ":" << tmp.lay << std::endl;
0629   }
0630   return id;
0631 }
0632 
0633 bool HCalSD::isItHF(const G4Step* aStep) {
0634   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0635   int levels = (touch->GetHistoryDepth()) + 1;
0636   for (unsigned int it = 0; it < hfNames.size(); ++it) {
0637     if (levels >= hfLevels[it]) {
0638       const G4LogicalVolume* lv = touch->GetVolume(levels - hfLevels[it])->GetLogicalVolume();
0639       if (lv == hfLV[it])
0640         return true;
0641     }
0642   }
0643   return false;
0644 }
0645 
0646 bool HCalSD::isItHF(const G4String& name) {
0647   for (const auto& nam : hfNames)
0648     if (name == static_cast<G4String>(nam)) {
0649       return true;
0650     }
0651   return false;
0652 }
0653 
0654 bool HCalSD::isItFibre(const G4LogicalVolume* lv) {
0655   for (auto lvol : fibreLV)
0656     if (lv == lvol) {
0657       return true;
0658     }
0659   return false;
0660 }
0661 
0662 bool HCalSD::isItFibre(const G4String& name) {
0663   for (const auto& nam : fibreNames)
0664     if (name == static_cast<G4String>(nam)) {
0665       return true;
0666     }
0667   return false;
0668 }
0669 
0670 bool HCalSD::isItPMT(const G4LogicalVolume* lv) {
0671   for (auto lvol : pmtLV)
0672     if (lv == lvol) {
0673       return true;
0674     }
0675   return false;
0676 }
0677 
0678 bool HCalSD::isItStraightBundle(const G4LogicalVolume* lv) {
0679   for (auto lvol : fibre1LV)
0680     if (lv == lvol) {
0681       return true;
0682     }
0683   return false;
0684 }
0685 
0686 bool HCalSD::isItConicalBundle(const G4LogicalVolume* lv) {
0687   for (auto lvol : fibre2LV)
0688     if (lv == lvol) {
0689       return true;
0690     }
0691   return false;
0692 }
0693 
0694 bool HCalSD::isItScintillator(const G4Material* mat) {
0695   for (auto amat : materials)
0696     if (amat == mat) {
0697       return true;
0698     }
0699   return false;
0700 }
0701 
0702 bool HCalSD::isItinFidVolume(const G4ThreeVector& hitPoint) {
0703   bool flag = true;
0704   if (applyFidCut) {
0705     int npmt = HFFibreFiducial::PMTNumber(hitPoint);
0706 #ifdef EDM_ML_DEBUG
0707     edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt << " for hit point " << hitPoint;
0708 #endif
0709     if (npmt <= 0)
0710       flag = false;
0711   }
0712 #ifdef EDM_ML_DEBUG
0713   edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint << " return flag " << flag;
0714 #endif
0715   return flag;
0716 }
0717 
0718 void HCalSD::getFromHFLibrary(const G4Step* aStep, bool& isKilled) {
0719   std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, isKilled, weight_, false);
0720   if (!isKilled || hits.empty()) {
0721     return;
0722   }
0723 
0724   int primaryID = setTrackID(aStep);
0725 
0726   // Reset entry point for new primary
0727   resetForNewPrimary(aStep);
0728 
0729   auto const theTrack = aStep->GetTrack();
0730   int det = 5;
0731 
0732   if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0733     edepositEM = 1. * GeV;
0734     edepositHAD = 0.;
0735   } else {
0736     edepositEM = 0.;
0737     edepositHAD = 1. * GeV;
0738   }
0739 #ifdef EDM_ML_DEBUG
0740   edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
0741                               << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0742                               << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV << " GeV";
0743 #endif
0744   for (unsigned int i = 0; i < hits.size(); ++i) {
0745     G4ThreeVector hitPoint = hits[i].position;
0746     if (isItinFidVolume(hitPoint)) {
0747       int depth = hits[i].depth;
0748       double time = hits[i].time;
0749       unsigned int unitID = setDetUnitId(det, hitPoint, depth);
0750       currentID.setID(unitID, time, primaryID, 0);
0751 #ifdef plotDebug
0752       plotProfile(aStep, hitPoint, 1.0 * GeV, time, depth);
0753       bool emType = G4TrackToParticleID::isGammaElectronPositron(theTrack->GetDefinition()->GetPDGEncoding());
0754       plotHF(hitPoint, emType);
0755 #endif
0756       processHit(aStep);
0757     }
0758   }
0759 }
0760 
0761 void HCalSD::hitForFibre(const G4Step* aStep) {  // if not ParamShower
0762 
0763   std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight_);
0764   if (hits.empty()) {
0765     return;
0766   }
0767 
0768   auto const theTrack = aStep->GetTrack();
0769   int primaryID = setTrackID(aStep);
0770   int det = 5;
0771 
0772   if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0773     edepositEM = 1. * GeV;
0774     edepositHAD = 0.;
0775   } else {
0776     edepositEM = 0.;
0777     edepositHAD = 1. * GeV;
0778   }
0779 
0780 #ifdef EDM_ML_DEBUG
0781   edm::LogVerbatim("HcalSim") << "HCalSD::hitForFibre " << hits.size() << " hits for " << GetName() << " of "
0782                               << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0783                               << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV << " GeV in detector type " << det;
0784 #endif
0785 
0786   for (unsigned int i = 0; i < hits.size(); ++i) {
0787     G4ThreeVector hitPoint = hits[i].position;
0788     if (isItinFidVolume(hitPoint)) {
0789       int depth = hits[i].depth;
0790       double time = hits[i].time;
0791       unsigned int unitID = setDetUnitId(det, hitPoint, depth);
0792       currentID.setID(unitID, time, primaryID, 0);
0793 #ifdef plotDebug
0794       plotProfile(aStep, hitPoint, edepositEM, time, depth);
0795       bool emType = (edepositEM > 0.) ? true : false;
0796       plotHF(hitPoint, emType);
0797 #endif
0798       processHit(aStep);
0799     }
0800   }
0801 }
0802 
0803 void HCalSD::getFromParam(const G4Step* aStep, bool& isKilled) {
0804   std::vector<HFShowerParam::Hit> hits = showerParam->getHits(aStep, weight_, isKilled);
0805   if (!isKilled || hits.empty()) {
0806     return;
0807   }
0808 
0809   int primaryID = setTrackID(aStep);
0810   int det = 5;
0811 
0812 #ifdef EDM_ML_DEBUG
0813   edm::LogVerbatim("HcalSim") << "HCalSD::getFromParam " << hits.size() << " hits for " << GetName() << " of "
0814                               << primaryID << " with " << aStep->GetTrack()->GetDefinition()->GetParticleName()
0815                               << " of " << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV
0816                               << " GeV in detector type " << det;
0817 #endif
0818   for (unsigned int i = 0; i < hits.size(); ++i) {
0819     G4ThreeVector hitPoint = hits[i].position;
0820     int depth = hits[i].depth;
0821     double time = hits[i].time;
0822     unsigned int unitID = setDetUnitId(det, hitPoint, depth);
0823     currentID.setID(unitID, time, primaryID, 0);
0824     edepositEM = hits[i].edep * GeV;
0825     edepositHAD = 0.;
0826 #ifdef plotDebug
0827     plotProfile(aStep, hitPoint, edepositEM, time, depth);
0828 #endif
0829     processHit(aStep);
0830   }
0831 }
0832 
0833 void HCalSD::getHitPMT(const G4Step* aStep) {
0834   auto const preStepPoint = aStep->GetPreStepPoint();
0835   auto const theTrack = aStep->GetTrack();
0836   double edep = showerPMT->getHits(aStep);
0837 
0838   if (edep >= 0.) {
0839     edep *= GeV;
0840     double etrack = preStepPoint->GetKineticEnergy();
0841     int primaryID = 0;
0842     if (etrack >= energyCut) {
0843       primaryID = theTrack->GetTrackID();
0844     } else {
0845       primaryID = theTrack->GetParentID();
0846       if (primaryID == 0)
0847         primaryID = theTrack->GetTrackID();
0848     }
0849     // Reset entry point for new primary
0850     resetForNewPrimary(aStep);
0851     //
0852     int det = static_cast<int>(HcalForward);
0853     const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0854     double rr = (hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
0855     double phi = (rr == 0. ? 0. : atan2(hitPoint.y(), hitPoint.x()));
0856     double etaR = showerPMT->getRadius();
0857     int depth = 1;
0858     if (etaR < 0) {
0859       depth = 2;
0860       etaR = -etaR;
0861     }
0862     if (hitPoint.z() < 0)
0863       etaR = -etaR;
0864 #ifdef EDM_ML_DEBUG
0865     edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR " << etaR << " phi " << phi / deg
0866                                 << " depth " << depth;
0867 #endif
0868     double time = (aStep->GetPostStepPoint()->GetGlobalTime());
0869     uint32_t unitID = 0;
0870     if (numberingFromDDD) {
0871       HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, etaR, phi, depth, 1);
0872       unitID = setDetUnitId(tmp);
0873     }
0874     currentID.setID(unitID, time, primaryID, 1);
0875 
0876     edepositHAD = aStep->GetTotalEnergyDeposit();
0877     edepositEM = -edepositHAD + edep;
0878 #ifdef plotDebug
0879     plotProfile(aStep, hitPoint, edep, time, depth);
0880 #endif
0881 #ifdef EDM_ML_DEBUG
0882     double beta = preStepPoint->GetBeta();
0883     edm::LogVerbatim("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName() << " of " << primaryID << " with "
0884                                 << theTrack->GetDefinition()->GetParticleName() << " of "
0885                                 << preStepPoint->GetKineticEnergy() / GeV << " GeV with velocity " << beta << " UnitID "
0886                                 << std::hex << unitID << std::dec;
0887 #endif
0888     processHit(aStep);
0889   }
0890 }
0891 
0892 void HCalSD::getHitFibreBundle(const G4Step* aStep, bool type) {
0893   auto const preStepPoint = aStep->GetPreStepPoint();
0894   auto const theTrack = aStep->GetTrack();
0895   double edep = showerBundle->getHits(aStep, type);
0896 
0897   if (edep >= 0.0) {
0898     edep *= GeV;
0899     double etrack = preStepPoint->GetKineticEnergy();
0900     int primaryID = 0;
0901     if (etrack >= energyCut) {
0902       primaryID = theTrack->GetTrackID();
0903     } else {
0904       primaryID = theTrack->GetParentID();
0905       if (primaryID == 0)
0906         primaryID = theTrack->GetTrackID();
0907     }
0908     // Reset entry point for new primary
0909     resetForNewPrimary(aStep);
0910     //
0911     int det = static_cast<int>(HcalForward);
0912     const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0913     double rr = hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y();
0914     double phi = rr == 0. ? 0. : atan2(hitPoint.y(), hitPoint.x());
0915     double etaR = showerBundle->getRadius();
0916     int depth = 1;
0917     if (etaR < 0.) {
0918       depth = 2;
0919       etaR = -etaR;
0920     }
0921     if (hitPoint.z() < 0.)
0922       etaR = -etaR;
0923 #ifdef EDM_ML_DEBUG
0924     edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR " << etaR << " phi " << phi / deg
0925                                 << " depth " << depth;
0926 #endif
0927     double time = (aStep->GetPostStepPoint()->GetGlobalTime());
0928     uint32_t unitID = 0;
0929     if (numberingFromDDD) {
0930       HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, etaR, phi, depth, 1);
0931       unitID = setDetUnitId(tmp);
0932     }
0933     if (type)
0934       currentID.setID(unitID, time, primaryID, 3);
0935     else
0936       currentID.setID(unitID, time, primaryID, 2);
0937 
0938     edepositHAD = aStep->GetTotalEnergyDeposit();
0939     edepositEM = -edepositHAD + edep;
0940 #ifdef plotDebug
0941     plotProfile(aStep, hitPoint, edep, time, depth);
0942 #endif
0943 #ifdef EDM_ML_DEBUG
0944     double beta = preStepPoint->GetBeta();
0945     edm::LogVerbatim("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for " << GetName() << " of " << primaryID
0946                                 << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0947                                 << preStepPoint->GetKineticEnergy() / GeV << " GeV with velocity " << beta << " UnitID "
0948                                 << std::hex << unitID << std::dec;
0949 #endif
0950     processHit(aStep);
0951   }  // non-zero energy deposit
0952 }
0953 
0954 void HCalSD::readWeightFromFile(const std::string& fName) {
0955   std::ifstream infile;
0956   int entry = 0;
0957   infile.open(fName.c_str(), std::ios::in);
0958   if (infile) {
0959     int det, zside, etaR, phi, lay;
0960     double wt;
0961     while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
0962       uint32_t id = HcalTestNumbering::packHcalIndex(det, zside, 1, etaR, phi, lay);
0963       layerWeights.insert(std::pair<uint32_t, double>(id, wt));
0964       ++entry;
0965 #ifdef EDM_ML_DEBUG
0966       edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry << " ID " << std::hex << id
0967                                   << std::dec << " (" << det << "/" << zside << "/1/" << etaR << "/" << phi << "/"
0968                                   << lay << ") Weight " << wt;
0969 #endif
0970     }
0971     infile.close();
0972   }
0973   edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry << " weights from " << fName;
0974   if (entry <= 0)
0975     useLayerWt = false;
0976 }
0977 
0978 double HCalSD::layerWeight(int det, const G4ThreeVector& pos, int depth, int lay) {
0979   double wt = 1.;
0980   if (numberingFromDDD) {
0981     //get the ID's as eta, phi, depth, ... indices
0982     HcalNumberingFromDDD::HcalID tmp =
0983         numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
0984     modifyDepth(tmp);
0985     uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1, tmp.etaR, tmp.phis, tmp.lay);
0986     std::map<uint32_t, double>::const_iterator ite = layerWeights.find(id);
0987     if (ite != layerWeights.end())
0988       wt = ite->second;
0989 #ifdef EDM_ML_DEBUG
0990     edm::LogVerbatim("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id << std::dec << " (" << tmp.subdet << "/"
0991                                 << tmp.zside << "/1/" << tmp.etaR << "/" << tmp.phis << "/" << tmp.lay << ") Weight "
0992                                 << wt;
0993 #endif
0994   }
0995   return wt;
0996 }
0997 
0998 void HCalSD::plotProfile(const G4Step* aStep, const G4ThreeVector& global, double edep, double time, int id) {
0999   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
1000   static const unsigned int names = 10;
1001   static const G4String modName[names] = {
1002       "HEModule", "HVQF", "HBModule", "MBAT", "MBBT", "MBBTC", "MBBT_R1P", "MBBT_R1M", "MBBT_R1PX", "MBBT_R1MX"};
1003   G4ThreeVector local;
1004   bool found = false;
1005   double depth = -2000;
1006   int idx = 4;
1007   for (int n = 0; n < touch->GetHistoryDepth(); ++n) {
1008     G4String name(static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(n)->GetName())));
1009 #ifdef EDM_ML_DEBUG
1010     edm::LogVerbatim("HcalSim") << "plotProfile Depth " << n << " Name " << name;
1011 #endif
1012     for (unsigned int ii = 0; ii < names; ++ii) {
1013       if (name == modName[ii]) {
1014         found = true;
1015         int dn = touch->GetHistoryDepth() - n;
1016         local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
1017         if (ii == 0) {
1018           depth = local.z() - 4006.5;
1019           idx = 1;
1020         } else if (ii == 1) {
1021           depth = local.z() + 825.0;
1022           idx = 3;
1023         } else if (ii == 2) {
1024           depth = local.x() - 1775.;
1025           idx = 0;
1026         } else {
1027           depth = local.y() + 15.;
1028           idx = 2;
1029         }
1030         break;
1031       }
1032     }
1033     if (found)
1034       break;
1035   }
1036   if (!found)
1037     depth = std::abs(global.z()) - 11500;
1038 #ifdef EDM_ML_DEBUG
1039   edm::LogVerbatim("HcalSim") << "plotProfile Found " << found << " Global " << global << " Local " << local
1040                               << " depth " << depth << " ID " << id << " EDEP " << edep << " Time " << time;
1041 #endif
1042   if (hit_[idx] != nullptr)
1043     hit_[idx]->Fill(edep);
1044   if (time_[idx] != nullptr)
1045     time_[idx]->Fill(time, edep);
1046   if (dist_[idx] != nullptr)
1047     dist_[idx]->Fill(depth, edep);
1048   int jd = 2 * idx + id - 7;
1049   if (jd >= 0 && jd < 4) {
1050     jd += 5;
1051     if (hit_[jd] != nullptr)
1052       hit_[jd]->Fill(edep);
1053     if (time_[jd] != nullptr)
1054       time_[jd]->Fill(time, edep);
1055     if (dist_[jd] != nullptr)
1056       dist_[jd]->Fill(depth, edep);
1057   }
1058 }
1059 
1060 void HCalSD::plotHF(const G4ThreeVector& hitPoint, bool emType) {
1061   double zv = std::abs(hitPoint.z()) - gpar[4];
1062   if (emType) {
1063     if (hzvem != nullptr)
1064       hzvem->Fill(zv);
1065   } else {
1066     if (hzvhad != nullptr)
1067       hzvhad->Fill(zv);
1068   }
1069 }
1070 
1071 void HCalSD::modifyDepth(HcalNumberingFromDDD::HcalID& id) {
1072   if (id.subdet == 4) {
1073     int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
1074     if (hcalConstants_->maxHFDepth(ieta, id.phis) > 2) {
1075       if (id.depth <= 2) {
1076         if (G4UniformRand() > 0.5)
1077           id.depth += 2;
1078       }
1079     }
1080   } else if ((id.subdet == 1 || id.subdet == 2) && testNumber) {
1081     id.depth = (depth_ == 0) ? 1 : 2;
1082   }
1083 }
1084 
1085 void HCalSD::initEvent(const BeginOfEvent*) {
1086 #ifdef printDebug
1087   detNull_ = {0, 0, 0, 0};
1088 #endif
1089 }
1090 
1091 void HCalSD::endEvent() {
1092 #ifdef printDebug
1093   int sum = detNull_[0] + detNull_[1] + detNull_[2];
1094   if (sum > 0)
1095     edm::LogVerbatim("HcalSim") << "NullDets " << detNull_[0] << " " << detNull_[1] << " " << detNull_[2] << " "
1096                                 << detNull_[3] << " " << (static_cast<float>(sum) / (sum + detNull_[3]));
1097 #endif
1098 }