File indexing completed on 2023-10-25 10:04:00
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: "
0343 << (aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0)->GetName())
0344
0345 << " isHF " << isHF << " darken " << (m_HFDarkening != nullptr) << " useParam "
0346 << useParam << " useShowerLibrary " << useShowerLibrary << " Muon? "
0347 << G4TrackToParticleID::isMuon(track) << " electron? "
0348 << G4TrackToParticleID::isGammaElectronPositron(track) << " Stable Hadron? "
0349 << G4TrackToParticleID::isStableHadronIon(track);
0350 #endif
0351 if (isHF) {
0352 if (m_HFDarkening) {
0353 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0354 const double invcm = 1. / CLHEP::cm;
0355 double r = hitPoint.perp() * invcm;
0356 double z = std::abs(hitPoint.z()) * invcm;
0357 double dose_acquired = 0.;
0358 if (z >= HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
0359 unsigned int hfZLayer = (unsigned int)((z - HFDarkening::lowZLimit) / 5);
0360 if (hfZLayer >= HFDarkening::upperZLimit)
0361 hfZLayer = (HFDarkening::upperZLimit - 1);
0362 float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
0363 for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
0364 dose_acquired = m_HFDarkening->dose(i, r);
0365 weight_ *= m_HFDarkening->degradation(normalized_lumi * dose_acquired);
0366 }
0367 }
0368 #ifdef EDM_ML_DEBUG
0369 edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary: HFLumiDarkening at "
0370 << "r= " << r << ", z= " << z << " Dose= " << dose_acquired << " weight= " << weight_;
0371 #endif
0372 }
0373
0374 if (useParam) {
0375 getFromParam(aStep, kill);
0376 #ifdef EDM_ML_DEBUG
0377 G4String nameVolume = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
0378 edm::LogVerbatim("HcalSim") << "HCalSD: " << getNumberOfHits() << " hits from parametrization in " << nameVolume
0379 << " for Track " << track->GetTrackID() << " ("
0380 << track->GetDefinition()->GetParticleName() << ")";
0381 #endif
0382 } else if (useShowerLibrary && !G4TrackToParticleID::isMuon(track)) {
0383 if (G4TrackToParticleID::isGammaElectronPositron(track) || G4TrackToParticleID::isStableHadronIon(track)) {
0384 #ifdef EDM_ML_DEBUG
0385 auto nameVolume = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
0386 edm::LogVerbatim("HcalSim") << "HCalSD: Starts shower library from " << nameVolume << " for Track "
0387 << track->GetTrackID() << " (" << track->GetDefinition()->GetParticleName() << ")";
0388
0389 #endif
0390 getFromHFLibrary(aStep, kill);
0391 }
0392 }
0393 }
0394 #ifdef EDM_ML_DEBUG
0395 edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary ID= " << track->GetTrackID() << " ("
0396 << track->GetDefinition()->GetParticleName() << ") kill= " << kill
0397 << " weight= " << weight_ << " depth= " << depth_ << " isHF: " << isHF;
0398 #endif
0399 return kill;
0400 }
0401
0402 double HCalSD::getEnergyDeposit(const G4Step* aStep) {
0403 double destep(0.0);
0404 auto const lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
0405 auto const theTrack = aStep->GetTrack();
0406
0407 if (isHF) {
0408 if (useShowerLibrary && G4TrackToParticleID::isMuon(theTrack) && isItFibre(lv)) {
0409 #ifdef EDM_ML_DEBUG
0410 edm::LogVerbatim("HcalSim") << "HCalSD: Hit at Fibre in LV " << lv->GetName() << " for track "
0411 << aStep->GetTrack()->GetTrackID() << " ("
0412 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0413 #endif
0414 hitForFibre(aStep);
0415 }
0416 return destep;
0417 }
0418
0419 if (isItPMT(lv)) {
0420 if (usePMTHit && showerPMT) {
0421 getHitPMT(aStep);
0422 }
0423 #ifdef EDM_ML_DEBUG
0424 edm::LogVerbatim("HcalSim") << "HCalSD: Hit from PMT parametrization in LV " << lv->GetName() << " for Track "
0425 << aStep->GetTrack()->GetTrackID() << " ("
0426 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0427 #endif
0428 return destep;
0429
0430 } else if (isItStraightBundle(lv)) {
0431 if (useFibreBundle && showerBundle) {
0432 getHitFibreBundle(aStep, false);
0433 }
0434 #ifdef EDM_ML_DEBUG
0435 edm::LogVerbatim("HcalSim") << "HCalSD: Hit from straight FibreBundle in LV: " << lv->GetName() << " for track "
0436 << aStep->GetTrack()->GetTrackID() << " ("
0437 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0438 #endif
0439 return destep;
0440
0441 } else if (isItConicalBundle(lv)) {
0442 if (useFibreBundle && showerBundle) {
0443 getHitFibreBundle(aStep, true);
0444 }
0445 #ifdef EDM_ML_DEBUG
0446 edm::LogVerbatim("HcalSim") << "HCalSD: Hit from conical FibreBundle PV: " << lv->GetName() << " for track "
0447 << aStep->GetTrack()->GetTrackID() << " ("
0448 << aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
0449 #endif
0450 return destep;
0451 }
0452
0453
0454 destep = aStep->GetTotalEnergyDeposit();
0455
0456 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0457 uint32_t detid = setDetUnitId(aStep);
0458 int det(0), ieta(0), phi(0), z(0), lay, depth(-1);
0459 if (testNumber) {
0460 HcalTestNumbering::unpackHcalIndex(detid, det, z, depth, ieta, phi, lay);
0461 if (z == 0) {
0462 z = -1;
0463 }
0464 } else {
0465 HcalDetId hcid(detid);
0466 det = hcid.subdetId();
0467 ieta = hcid.ietaAbs();
0468 phi = hcid.iphi();
0469 z = hcid.zside();
0470 }
0471 lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
0472 #ifdef EDM_ML_DEBUG
0473 edm::LogVerbatim("HcalSim") << "HCalSD: det: " << det << " ieta: " << ieta << " iphi: " << phi << " zside " << z
0474 << " lay: " << lay - 2;
0475 #endif
0476 if (depth_ == 0 && (det == 1 || det == 2) && ((!testNumber) || neutralDensity))
0477 weight_ = hcalConstants_->getLayer0Wt(det, phi, z);
0478 if (useLayerWt) {
0479 G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
0480 weight_ = layerWeight(det + 2, hitPoint, depth_, lay);
0481 }
0482
0483 if (agingFlagHB && m_HBDarkening && det == 1) {
0484 double dweight = m_HBDarkening->degradation(deliveredLumi, ieta, lay);
0485 weight_ *= dweight;
0486 #ifdef EDM_ML_DEBUG
0487 edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi << " coefficient = " << dweight
0488 << " Weight= " << weight_;
0489 #endif
0490 }
0491
0492 if (agingFlagHE && m_HEDarkening && det == 2) {
0493 double dweight = m_HEDarkening->degradation(deliveredLumi, ieta, lay);
0494 weight_ *= dweight;
0495 #ifdef EDM_ML_DEBUG
0496 edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi << " coefficient = " << dweight
0497 << " Weight= " << weight_;
0498 #endif
0499 }
0500
0501 if (suppressHeavy) {
0502 TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
0503 if (trkInfo) {
0504 int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0505 if (!(trkInfo->isPrimary())) {
0506 double ke = theTrack->GetKineticEnergy();
0507 if (pdg / 1000000000 == 1 && (pdg / 10000) % 100 > 0 && (pdg / 10) % 100 > 0 && ke < kmaxIon)
0508 weight_ = 0;
0509 if ((pdg == 2212) && (ke < kmaxProton))
0510 weight_ = 0;
0511 if ((pdg == 2112) && (ke < kmaxNeutron))
0512 weight_ = 0;
0513 }
0514 }
0515 }
0516 double wt0(1.0);
0517 if (useBirk) {
0518 const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
0519 if (isItScintillator(mat))
0520 wt0 = getAttenuation(aStep, birk1, birk2, birk3);
0521 }
0522 weight_ *= wt0;
0523 double wt1 = getResponseWt(theTrack);
0524 double wt2 = theTrack->GetWeight();
0525 double edep = weight_ * wt1 * destep;
0526 if (wt2 > 0.0) {
0527 edep *= wt2;
0528 }
0529 #ifdef EDM_ML_DEBUG
0530 edm::LogVerbatim("HcalSim") << "HCalSD: edep= " << edep << " Det: " << det + 2 << " depth= " << depth_
0531 << " weight= " << weight_ << " wt0= " << wt0 << " wt1= " << wt1 << " wt2= " << wt2;
0532 #endif
0533 return edep;
0534 }
0535
0536 uint32_t HCalSD::setDetUnitId(const G4Step* aStep) {
0537 auto const prePoint = aStep->GetPreStepPoint();
0538 auto const touch = prePoint->GetTouchable();
0539 const G4ThreeVector& hitPoint = prePoint->GetPosition();
0540
0541 int depth = (touch->GetReplicaNumber(0)) % 10 + 1;
0542 int lay = (touch->GetReplicaNumber(0) / 10) % 100 + 1;
0543 int det = (touch->GetReplicaNumber(1)) / 1000;
0544
0545 uint32_t idx = setDetUnitId(det, hitPoint, depth, lay);
0546 #ifdef EDM_ML_DEBUG
0547 if (depth == 1) {
0548 edm::LogVerbatim("HcalSim") << "HCalSD: Check for " << det << ":" << depth << ":" << lay << " ID " << std::hex
0549 << idx << std::dec;
0550 int det0, z0, depth0, eta0, phi0, lay0(-1);
0551 if (testNumber) {
0552 HcalTestNumbering::unpackHcalIndex(idx, det0, z0, depth0, eta0, phi0, lay0);
0553 } else {
0554 HcalDetId hcid0(idx);
0555 det0 = hcid0.subdetId();
0556 eta0 = hcid0.ietaAbs();
0557 phi0 = hcid0.iphi();
0558 z0 = hcid0.zside();
0559 depth0 = hcid0.depth();
0560 }
0561 edm::LogVerbatim("HcalSim") << "HCalSD: det|z|depth|eta|phi|lay " << det0 << ":" << z0 << ":" << depth0 << ":"
0562 << eta0 << ":" << phi0 << ":" << lay0;
0563 printVolume(touch);
0564 }
0565 #endif
0566 return idx;
0567 }
0568
0569 void HCalSD::setNumberingScheme(HcalNumberingScheme* scheme) {
0570 if (scheme != nullptr) {
0571 edm::LogVerbatim("HcalSim") << "HCalSD: updates numbering scheme for " << GetName();
0572 numberingScheme.reset(scheme);
0573 }
0574 }
0575
0576 void HCalSD::update(const BeginOfJob* job) {}
0577
0578 void HCalSD::initRun() {}
0579
0580 bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
0581 double threshold = 0;
0582 DetId theId(aHit->getUnitID());
0583 switch (theId.subdetId()) {
0584 case HcalBarrel:
0585 threshold = eminHitHB;
0586 break;
0587 case HcalEndcap:
0588 threshold = eminHitHE;
0589 break;
0590 case HcalOuter:
0591 threshold = eminHitHO;
0592 break;
0593 case HcalForward:
0594 threshold = eminHitHF;
0595 break;
0596 default:
0597 break;
0598 }
0599 return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
0600 }
0601
0602 uint32_t HCalSD::setDetUnitId(int det, const G4ThreeVector& pos, int depth, int lay = 1) {
0603 uint32_t id = 0;
0604 if (det == 0) {
0605 #ifdef printDebug
0606 double eta = std::abs(pos.eta());
0607 #endif
0608 if (std::abs(pos.z()) > maxZ_) {
0609 det = 5;
0610 #ifdef printDebug
0611 if (eta < 2.868)
0612 ++detNull_[2];
0613 #endif
0614 } else if (!(hcalConstants_->isHE())) {
0615 det = 3;
0616 #ifdef printDebug
0617 ++detNull_[0];
0618 #endif
0619 } else {
0620 double minR = minRoff_ + slopeHE_ * std::abs(pos.z());
0621 double maxR = maxRoff_ + slopeHE_ * std::abs(pos.z());
0622 det = ((pos.perp() > minR) && (pos.perp() < maxR)) ? 4 : 3;
0623 #ifdef printDebug
0624 ++detNull_[det - 3];
0625 #endif
0626 }
0627 #ifdef printDEBUG
0628 edm::LogVerbatim("HcalSim") << "Position " << pos.perp() << ":" << std::abs(pos.z()) << " Limits "
0629 << !(hcalConstants_->isHE()) << ":" << maxZ_ << " det " << det;
0630 } else {
0631 ++detNull_[3];
0632 #endif
0633 }
0634
0635 if (numberingFromDDD.get()) {
0636
0637 HcalNumberingFromDDD::HcalID tmp =
0638 numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
0639 id = setDetUnitId(tmp);
0640 }
0641 return id;
0642 }
0643
0644 uint32_t HCalSD::setDetUnitId(HcalNumberingFromDDD::HcalID& tmp) {
0645 modifyDepth(tmp);
0646 uint32_t id = (numberingScheme.get()) ? numberingScheme->getUnitID(tmp) : 0;
0647 if ((!testNumber) && m_HcalTestNS.get()) {
0648 bool ok = m_HcalTestNS->compare(tmp, id);
0649 if (!ok)
0650 edm::LogWarning("HcalSim") << "Det ID from HCalSD " << HcalDetId(id) << " " << std::hex << id << std::dec
0651 << " does not match one from relabller for " << tmp.subdet << ":" << tmp.etaR << ":"
0652 << tmp.phi << ":" << tmp.phis << ":" << tmp.depth << ":" << tmp.lay << std::endl;
0653 }
0654 return id;
0655 }
0656
0657 bool HCalSD::isItHF(const G4Step* aStep) {
0658 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0659 int levels = (touch->GetHistoryDepth()) + 1;
0660 for (unsigned int it = 0; it < hfNames.size(); ++it) {
0661 if (levels >= hfLevels[it]) {
0662 const G4LogicalVolume* lv = touch->GetVolume(levels - hfLevels[it])->GetLogicalVolume();
0663 if (lv == hfLV[it])
0664 return true;
0665 }
0666 }
0667 return false;
0668 }
0669
0670 bool HCalSD::isItHF(const G4String& name) {
0671 for (const auto& nam : hfNames)
0672 if (name == static_cast<G4String>(nam)) {
0673 return true;
0674 }
0675 return false;
0676 }
0677
0678 bool HCalSD::isItFibre(const G4LogicalVolume* lv) {
0679 for (auto lvol : fibreLV)
0680 if (lv == lvol) {
0681 return true;
0682 }
0683 return false;
0684 }
0685
0686 bool HCalSD::isItFibre(const G4String& name) {
0687 for (const auto& nam : fibreNames)
0688 if (name == static_cast<G4String>(nam)) {
0689 return true;
0690 }
0691 return false;
0692 }
0693
0694 bool HCalSD::isItPMT(const G4LogicalVolume* lv) {
0695 for (auto lvol : pmtLV)
0696 if (lv == lvol) {
0697 return true;
0698 }
0699 return false;
0700 }
0701
0702 bool HCalSD::isItStraightBundle(const G4LogicalVolume* lv) {
0703 for (auto lvol : fibre1LV)
0704 if (lv == lvol) {
0705 return true;
0706 }
0707 return false;
0708 }
0709
0710 bool HCalSD::isItConicalBundle(const G4LogicalVolume* lv) {
0711 for (auto lvol : fibre2LV)
0712 if (lv == lvol) {
0713 return true;
0714 }
0715 return false;
0716 }
0717
0718 bool HCalSD::isItScintillator(const G4Material* mat) {
0719 for (auto amat : materials)
0720 if (amat == mat) {
0721 return true;
0722 }
0723 return false;
0724 }
0725
0726 bool HCalSD::isItinFidVolume(const G4ThreeVector& hitPoint) {
0727 bool flag = true;
0728 if (applyFidCut) {
0729 int npmt = HFFibreFiducial::PMTNumber(hitPoint);
0730 #ifdef EDM_ML_DEBUG
0731 edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt << " for hit point " << hitPoint;
0732 #endif
0733 if (npmt <= 0)
0734 flag = false;
0735 }
0736 #ifdef EDM_ML_DEBUG
0737 edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint << " return flag " << flag;
0738 #endif
0739 return flag;
0740 }
0741
0742 void HCalSD::getFromHFLibrary(const G4Step* aStep, bool& isKilled) {
0743 std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, isKilled, weight_, false);
0744 if (!isKilled || hits.empty()) {
0745 return;
0746 }
0747
0748 int primaryID = setTrackID(aStep);
0749
0750
0751 resetForNewPrimary(aStep);
0752
0753 auto const theTrack = aStep->GetTrack();
0754 int det = 5;
0755
0756 if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0757 edepositEM = 1. * GeV;
0758 edepositHAD = 0.;
0759 } else {
0760 edepositEM = 0.;
0761 edepositHAD = 1. * GeV;
0762 }
0763 #ifdef EDM_ML_DEBUG
0764 edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary " << hits.size() << " hits for " << GetName() << " of "
0765 << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0766 << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV << " GeV";
0767 #endif
0768 for (unsigned int i = 0; i < hits.size(); ++i) {
0769 G4ThreeVector hitPoint = hits[i].position;
0770 if (isItinFidVolume(hitPoint)) {
0771 int depth = hits[i].depth;
0772 double time = hits[i].time;
0773 unsigned int unitID = setDetUnitId(det, hitPoint, depth);
0774 currentID[0].setID(unitID, time, primaryID, 0);
0775 #ifdef plotDebug
0776 plotProfile(aStep, hitPoint, 1.0 * GeV, time, depth);
0777 bool emType = G4TrackToParticleID::isGammaElectronPositron(theTrack->GetDefinition()->GetPDGEncoding());
0778 plotHF(hitPoint, emType);
0779 #endif
0780 processHit(aStep);
0781 }
0782 }
0783 }
0784
0785 void HCalSD::hitForFibre(const G4Step* aStep) {
0786
0787 std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight_);
0788 if (hits.empty()) {
0789 return;
0790 }
0791
0792 auto const theTrack = aStep->GetTrack();
0793 int primaryID = setTrackID(aStep);
0794 int det = 5;
0795
0796 if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0797 edepositEM = 1. * GeV;
0798 edepositHAD = 0.;
0799 } else {
0800 edepositEM = 0.;
0801 edepositHAD = 1. * GeV;
0802 }
0803
0804 #ifdef EDM_ML_DEBUG
0805 edm::LogVerbatim("HcalSim") << "HCalSD::hitForFibre " << hits.size() << " hits for " << GetName() << " of "
0806 << primaryID << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0807 << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV << " GeV in detector type " << det;
0808 #endif
0809
0810 for (unsigned int i = 0; i < hits.size(); ++i) {
0811 G4ThreeVector hitPoint = hits[i].position;
0812 if (isItinFidVolume(hitPoint)) {
0813 int depth = hits[i].depth;
0814 double time = hits[i].time;
0815 unsigned int unitID = setDetUnitId(det, hitPoint, depth);
0816 currentID[0].setID(unitID, time, primaryID, 0);
0817 #ifdef plotDebug
0818 plotProfile(aStep, hitPoint, edepositEM, time, depth);
0819 bool emType = (edepositEM > 0.) ? true : false;
0820 plotHF(hitPoint, emType);
0821 #endif
0822 processHit(aStep);
0823 }
0824 }
0825 }
0826
0827 void HCalSD::getFromParam(const G4Step* aStep, bool& isKilled) {
0828 std::vector<HFShowerParam::Hit> hits = showerParam->getHits(aStep, weight_, isKilled);
0829 if (!isKilled || hits.empty()) {
0830 return;
0831 }
0832
0833 int primaryID = setTrackID(aStep);
0834 int det = 5;
0835
0836 #ifdef EDM_ML_DEBUG
0837 edm::LogVerbatim("HcalSim") << "HCalSD::getFromParam " << hits.size() << " hits for " << GetName() << " of "
0838 << primaryID << " with " << aStep->GetTrack()->GetDefinition()->GetParticleName()
0839 << " of " << aStep->GetPreStepPoint()->GetKineticEnergy() / GeV
0840 << " GeV in detector type " << det;
0841 #endif
0842 for (unsigned int i = 0; i < hits.size(); ++i) {
0843 G4ThreeVector hitPoint = hits[i].position;
0844 int depth = hits[i].depth;
0845 double time = hits[i].time;
0846 unsigned int unitID = setDetUnitId(det, hitPoint, depth);
0847 currentID[0].setID(unitID, time, primaryID, 0);
0848 edepositEM = hits[i].edep * GeV;
0849 edepositHAD = 0.;
0850 #ifdef plotDebug
0851 plotProfile(aStep, hitPoint, edepositEM, time, depth);
0852 #endif
0853 processHit(aStep);
0854 }
0855 }
0856
0857 void HCalSD::getHitPMT(const G4Step* aStep) {
0858 auto const preStepPoint = aStep->GetPreStepPoint();
0859 auto const theTrack = aStep->GetTrack();
0860 double edep = showerPMT->getHits(aStep);
0861
0862 if (edep >= 0.) {
0863 edep *= GeV;
0864 double etrack = preStepPoint->GetKineticEnergy();
0865 int primaryID = 0;
0866 if (etrack >= energyCut) {
0867 primaryID = theTrack->GetTrackID();
0868 } else {
0869 primaryID = theTrack->GetParentID();
0870 if (primaryID == 0)
0871 primaryID = theTrack->GetTrackID();
0872 }
0873
0874 resetForNewPrimary(aStep);
0875
0876 int det = static_cast<int>(HcalForward);
0877 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0878 double rr = (hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
0879 double phi = (rr == 0. ? 0. : atan2(hitPoint.y(), hitPoint.x()));
0880 double etaR = showerPMT->getRadius();
0881 int depth = 1;
0882 if (etaR < 0) {
0883 depth = 2;
0884 etaR = -etaR;
0885 }
0886 if (hitPoint.z() < 0)
0887 etaR = -etaR;
0888 #ifdef EDM_ML_DEBUG
0889 edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR " << etaR << " phi " << phi / deg
0890 << " depth " << depth;
0891 #endif
0892 double time = (aStep->GetPostStepPoint()->GetGlobalTime());
0893 uint32_t unitID = 0;
0894 if (numberingFromDDD) {
0895 HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, etaR, phi, depth, 1);
0896 unitID = setDetUnitId(tmp);
0897 }
0898 currentID[0].setID(unitID, time, primaryID, 1);
0899
0900 edepositHAD = aStep->GetTotalEnergyDeposit();
0901 edepositEM = -edepositHAD + edep;
0902 #ifdef plotDebug
0903 plotProfile(aStep, hitPoint, edep, time, depth);
0904 #endif
0905 #ifdef EDM_ML_DEBUG
0906 double beta = preStepPoint->GetBeta();
0907 edm::LogVerbatim("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName() << " of " << primaryID << " with "
0908 << theTrack->GetDefinition()->GetParticleName() << " of "
0909 << preStepPoint->GetKineticEnergy() / GeV << " GeV with velocity " << beta << " UnitID "
0910 << std::hex << unitID << std::dec;
0911 #endif
0912 processHit(aStep);
0913 }
0914 }
0915
0916 void HCalSD::getHitFibreBundle(const G4Step* aStep, bool type) {
0917 auto const preStepPoint = aStep->GetPreStepPoint();
0918 auto const theTrack = aStep->GetTrack();
0919 double edep = showerBundle->getHits(aStep, type);
0920
0921 if (edep >= 0.0) {
0922 edep *= GeV;
0923 double etrack = preStepPoint->GetKineticEnergy();
0924 int primaryID = 0;
0925 if (etrack >= energyCut) {
0926 primaryID = theTrack->GetTrackID();
0927 } else {
0928 primaryID = theTrack->GetParentID();
0929 if (primaryID == 0)
0930 primaryID = theTrack->GetTrackID();
0931 }
0932
0933 resetForNewPrimary(aStep);
0934
0935 int det = static_cast<int>(HcalForward);
0936 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0937 double rr = hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y();
0938 double phi = rr == 0. ? 0. : atan2(hitPoint.y(), hitPoint.x());
0939 double etaR = showerBundle->getRadius();
0940 int depth = 1;
0941 if (etaR < 0.) {
0942 depth = 2;
0943 etaR = -etaR;
0944 }
0945 if (hitPoint.z() < 0.)
0946 etaR = -etaR;
0947 #ifdef EDM_ML_DEBUG
0948 edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR " << etaR << " phi " << phi / deg
0949 << " depth " << depth;
0950 #endif
0951 double time = (aStep->GetPostStepPoint()->GetGlobalTime());
0952 uint32_t unitID = 0;
0953 if (numberingFromDDD) {
0954 HcalNumberingFromDDD::HcalID tmp = numberingFromDDD->unitID(det, etaR, phi, depth, 1);
0955 unitID = setDetUnitId(tmp);
0956 }
0957 if (type)
0958 currentID[0].setID(unitID, time, primaryID, 3);
0959 else
0960 currentID[0].setID(unitID, time, primaryID, 2);
0961
0962 edepositHAD = aStep->GetTotalEnergyDeposit();
0963 edepositEM = -edepositHAD + edep;
0964 #ifdef plotDebug
0965 plotProfile(aStep, hitPoint, edep, time, depth);
0966 #endif
0967 #ifdef EDM_ML_DEBUG
0968 double beta = preStepPoint->GetBeta();
0969 edm::LogVerbatim("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for " << GetName() << " of " << primaryID
0970 << " with " << theTrack->GetDefinition()->GetParticleName() << " of "
0971 << preStepPoint->GetKineticEnergy() / GeV << " GeV with velocity " << beta << " UnitID "
0972 << std::hex << unitID << std::dec;
0973 #endif
0974 processHit(aStep);
0975 }
0976 }
0977
0978 void HCalSD::readWeightFromFile(const std::string& fName) {
0979 std::ifstream infile;
0980 int entry = 0;
0981 infile.open(fName.c_str(), std::ios::in);
0982 if (infile) {
0983 int det, zside, etaR, phi, lay;
0984 double wt;
0985 while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
0986 uint32_t id = HcalTestNumbering::packHcalIndex(det, zside, 1, etaR, phi, lay);
0987 layerWeights.insert(std::pair<uint32_t, double>(id, wt));
0988 ++entry;
0989 #ifdef EDM_ML_DEBUG
0990 edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry << " ID " << std::hex << id
0991 << std::dec << " (" << det << "/" << zside << "/1/" << etaR << "/" << phi << "/"
0992 << lay << ") Weight " << wt;
0993 #endif
0994 }
0995 infile.close();
0996 }
0997 edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry << " weights from " << fName;
0998 if (entry <= 0)
0999 useLayerWt = false;
1000 }
1001
1002 double HCalSD::layerWeight(int det, const G4ThreeVector& pos, int depth, int lay) {
1003 double wt = 1.;
1004 if (numberingFromDDD) {
1005
1006 HcalNumberingFromDDD::HcalID tmp =
1007 numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
1008 modifyDepth(tmp);
1009 uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1, tmp.etaR, tmp.phis, tmp.lay);
1010 std::map<uint32_t, double>::const_iterator ite = layerWeights.find(id);
1011 if (ite != layerWeights.end())
1012 wt = ite->second;
1013 #ifdef EDM_ML_DEBUG
1014 edm::LogVerbatim("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id << std::dec << " (" << tmp.subdet << "/"
1015 << tmp.zside << "/1/" << tmp.etaR << "/" << tmp.phis << "/" << tmp.lay << ") Weight "
1016 << wt;
1017 #endif
1018 }
1019 return wt;
1020 }
1021
1022 void HCalSD::plotProfile(const G4Step* aStep, const G4ThreeVector& global, double edep, double time, int id) {
1023 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
1024 static const unsigned int names = 10;
1025 static const G4String modName[names] = {
1026 "HEModule", "HVQF", "HBModule", "MBAT", "MBBT", "MBBTC", "MBBT_R1P", "MBBT_R1M", "MBBT_R1PX", "MBBT_R1MX"};
1027 G4ThreeVector local;
1028 bool found = false;
1029 double depth = -2000;
1030 int idx = 4;
1031 for (int n = 0; n < touch->GetHistoryDepth(); ++n) {
1032 G4String name(static_cast<std::string>(dd4hep::dd::noNamespace(touch->GetVolume(n)->GetName())));
1033 #ifdef EDM_ML_DEBUG
1034 edm::LogVerbatim("HcalSim") << "plotProfile Depth " << n << " Name " << name;
1035 #endif
1036 for (unsigned int ii = 0; ii < names; ++ii) {
1037 if (name == modName[ii]) {
1038 found = true;
1039 int dn = touch->GetHistoryDepth() - n;
1040 local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
1041 if (ii == 0) {
1042 depth = local.z() - 4006.5;
1043 idx = 1;
1044 } else if (ii == 1) {
1045 depth = local.z() + 825.0;
1046 idx = 3;
1047 } else if (ii == 2) {
1048 depth = local.x() - 1775.;
1049 idx = 0;
1050 } else {
1051 depth = local.y() + 15.;
1052 idx = 2;
1053 }
1054 break;
1055 }
1056 }
1057 if (found)
1058 break;
1059 }
1060 if (!found)
1061 depth = std::abs(global.z()) - 11500;
1062 #ifdef EDM_ML_DEBUG
1063 edm::LogVerbatim("HcalSim") << "plotProfile Found " << found << " Global " << global << " Local " << local
1064 << " depth " << depth << " ID " << id << " EDEP " << edep << " Time " << time;
1065 #endif
1066 if (hit_[idx] != nullptr)
1067 hit_[idx]->Fill(edep);
1068 if (time_[idx] != nullptr)
1069 time_[idx]->Fill(time, edep);
1070 if (dist_[idx] != nullptr)
1071 dist_[idx]->Fill(depth, edep);
1072 int jd = 2 * idx + id - 7;
1073 if (jd >= 0 && jd < 4) {
1074 jd += 5;
1075 if (hit_[jd] != nullptr)
1076 hit_[jd]->Fill(edep);
1077 if (time_[jd] != nullptr)
1078 time_[jd]->Fill(time, edep);
1079 if (dist_[jd] != nullptr)
1080 dist_[jd]->Fill(depth, edep);
1081 }
1082 }
1083
1084 void HCalSD::plotHF(const G4ThreeVector& hitPoint, bool emType) {
1085 double zv = std::abs(hitPoint.z()) - gpar[4];
1086 if (emType) {
1087 if (hzvem != nullptr)
1088 hzvem->Fill(zv);
1089 } else {
1090 if (hzvhad != nullptr)
1091 hzvhad->Fill(zv);
1092 }
1093 }
1094
1095 void HCalSD::modifyDepth(HcalNumberingFromDDD::HcalID& id) {
1096 if (id.subdet == 4) {
1097 int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
1098 if (hcalConstants_->maxHFDepth(ieta, id.phis) > 2) {
1099 if (id.depth <= 2) {
1100 if (G4UniformRand() > 0.5)
1101 id.depth += 2;
1102 }
1103 }
1104 } else if ((id.subdet == 1 || id.subdet == 2) && testNumber) {
1105 id.depth = (depth_ == 0) ? 1 : 2;
1106 }
1107 }
1108
1109 void HCalSD::initEvent(const BeginOfEvent*) {
1110 #ifdef printDebug
1111 detNull_ = {0, 0, 0, 0};
1112 #endif
1113 }
1114
1115 void HCalSD::endEvent() {
1116 #ifdef printDebug
1117 int sum = detNull_[0] + detNull_[1] + detNull_[2];
1118 if (sum > 0)
1119 edm::LogVerbatim("HcalSim") << "NullDets " << detNull_[0] << " " << detNull_[1] << " " << detNull_[2] << " "
1120 << detNull_[3] << " " << (static_cast<float>(sum) / (sum + detNull_[3]));
1121 #endif
1122 }
1123
1124 void HCalSD::printVolume(const G4VTouchable* touch) const {
1125 if (touch) {
1126 #ifdef EDM_ML_DEBUG
1127 int level = ((touch->GetHistoryDepth()) + 1);
1128 edm::LogVerbatim("CaloSimX") << "HCalSD::printVolume with " << level << " levels";
1129 static const std::string unknown("Unknown");
1130
1131 for (int ii = 0; ii < level; ii++) {
1132 int i = level - ii - 1;
1133 G4VPhysicalVolume* pv = touch->GetVolume(i);
1134 G4String name = (pv != nullptr) ? pv->GetName() : unknown;
1135 G4int copyno = touch->GetReplicaNumber(i);
1136 edm::LogVerbatim("HcalSim") << "[" << ii << "] " << name << ":" << copyno;
1137 }
1138 #endif
1139 }
1140 }
