File indexing completed on 2021-12-06 04:01:06
0001
0002
0003
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
0036
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
0075
0076
0077
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
0153 setParameterized(true);
0154
0155 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0156
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
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
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
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
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
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
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())) {
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
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
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) {
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
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
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 }
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
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 }
