Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-06 04:01:06

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