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