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