File indexing completed on 2023-05-14 22:52:33
0001
0002
0003
0004
0005
0006 #include "SimG4CMS/Calo/interface/CaloSD.h"
0007 #include "SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h"
0008 #include "SimG4Core/Notification/interface/TrackInformation.h"
0009 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0010 #include "SimG4Core/Notification/interface/SimTrackManager.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013
0014 #include "G4EventManager.hh"
0015 #include "G4LogicalVolumeStore.hh"
0016 #include "G4LogicalVolume.hh"
0017 #include "G4SDManager.hh"
0018 #include "G4Step.hh"
0019 #include "G4Track.hh"
0020 #include "G4VProcess.hh"
0021 #include "G4GFlashSpot.hh"
0022 #include "G4ParticleTable.hh"
0023 #include "G4SystemOfUnits.hh"
0024 #include "G4PhysicalConstants.hh"
0025 #include "DD4hep/Filter.h"
0026
0027 #include <fstream>
0028 #include <memory>
0029 #include <sstream>
0030
0031
0032
0033 CaloSD::CaloSD(const std::string& name,
0034 const SensitiveDetectorCatalog& clg,
0035 edm::ParameterSet const& p,
0036 const SimTrackManager* manager,
0037 float timeSliceUnit,
0038 bool ignoreTkID)
0039 : SensitiveCaloDetector(name, clg),
0040 G4VGFlashSensitiveDetector(),
0041 eminHit(0.),
0042 currentHit(nullptr),
0043 m_trackManager(manager),
0044 theHC(nullptr),
0045 ignoreTrackID(ignoreTkID),
0046 hcID(-1),
0047 timeSlice(timeSliceUnit),
0048 eminHitD(0.) {
0049
0050 bool dd4hep = p.getParameter<bool>("g4GeometryDD4hepSource");
0051 int addlevel = dd4hep ? 1 : 0;
0052 edm::ParameterSet m_CaloSD = p.getParameter<edm::ParameterSet>("CaloSD");
0053 energyCut = m_CaloSD.getParameter<double>("EminTrack") * CLHEP::GeV;
0054 tmaxHit = m_CaloSD.getParameter<double>("TmaxHit") * CLHEP::ns;
0055 std::vector<double> eminHits = m_CaloSD.getParameter<std::vector<double>>("EminHits");
0056 std::vector<double> tmaxHits = m_CaloSD.getParameter<std::vector<double>>("TmaxHits");
0057 std::vector<std::string> hcn = m_CaloSD.getParameter<std::vector<std::string>>("HCNames");
0058 std::vector<int> useResMap = m_CaloSD.getParameter<std::vector<int>>("UseResponseTables");
0059 std::vector<double> eminHitX = m_CaloSD.getParameter<std::vector<double>>("EminHitsDepth");
0060 suppressHeavy = m_CaloSD.getParameter<bool>("SuppressHeavy");
0061 kmaxIon = m_CaloSD.getParameter<double>("IonThreshold") * CLHEP::MeV;
0062 kmaxProton = m_CaloSD.getParameter<double>("ProtonThreshold") * CLHEP::MeV;
0063 kmaxNeutron = m_CaloSD.getParameter<double>("NeutronThreshold") * CLHEP::MeV;
0064 nCheckedHits = m_CaloSD.getUntrackedParameter<int>("CheckHits", 25);
0065 useMap = m_CaloSD.getUntrackedParameter<bool>("UseMap", true);
0066 int verbn = m_CaloSD.getUntrackedParameter<int>("Verbosity", 0);
0067 corrTOFBeam = m_CaloSD.getParameter<bool>("CorrectTOFBeam");
0068 double beamZ = m_CaloSD.getParameter<double>("BeamPosition") * CLHEP::cm;
0069 correctT = beamZ / CLHEP::c_light / CLHEP::nanosecond;
0070 doFineCalo_ = m_CaloSD.getParameter<bool>("DoFineCalo");
0071 eMinFine_ = m_CaloSD.getParameter<double>("EminFineTrack") * CLHEP::MeV;
0072 std::vector<std::string> fineNames = m_CaloSD.getParameter<std::vector<std::string>>("FineCaloNames");
0073 std::vector<int> fineLevels = m_CaloSD.getParameter<std::vector<int>>("FineCaloLevels");
0074 std::vector<int> useFines = m_CaloSD.getParameter<std::vector<int>>("UseFineCalo");
0075 for (auto& level : fineLevels)
0076 level += addlevel;
0077
0078 SetVerboseLevel(verbn);
0079 meanResponse.reset(nullptr);
0080 for (unsigned int k = 0; k < hcn.size(); ++k) {
0081 if (name == hcn[k]) {
0082 if (k < eminHits.size())
0083 eminHit = eminHits[k] * CLHEP::MeV;
0084 if (k < eminHitX.size())
0085 eminHitD = eminHitX[k] * CLHEP::MeV;
0086 if (k < tmaxHits.size())
0087 tmaxHit = tmaxHits[k] * CLHEP::ns;
0088 if (k < useResMap.size() && useResMap[k] > 0) {
0089 meanResponse = std::make_unique<CaloMeanResponse>(p);
0090 break;
0091 }
0092 }
0093 }
0094 slave = std::make_unique<CaloSlaveSD>(name);
0095
0096 currentID = CaloHitID(timeSlice, ignoreTrackID);
0097 previousID = CaloHitID(timeSlice, ignoreTrackID);
0098 isParameterized = false;
0099
0100 entrancePoint.set(0., 0., 0.);
0101 entranceLocal.set(0., 0., 0.);
0102 posGlobal.set(0., 0., 0.);
0103 incidentEnergy = edepositEM = edepositHAD = 0.f;
0104
0105 primAncestor = cleanIndex = totalHits = primIDSaved = 0;
0106 forceSave = false;
0107
0108 edm::LogVerbatim("CaloSim") << "CaloSD: Minimum energy of track for saving it " << energyCut / CLHEP::GeV
0109 << " GeV\n Use of HitID Map " << useMap << "\n Check last " << nCheckedHits
0110 << " before saving the hit\n Correct TOF globally by " << correctT
0111 << " ns (Flag =" << corrTOFBeam << ")\n Save hits recorded before " << tmaxHit
0112 << " ns and if energy is above " << eminHit / CLHEP::MeV << " MeV (for depth 0) or "
0113 << eminHitD / CLHEP::MeV << " MeV (for nonzero depths);\n Time Slice Unit "
0114 << timeSlice << "\nIgnore TrackID Flag " << ignoreTrackID << " doFineCalo flag "
0115 << doFineCalo_ << "\nBeam Position " << beamZ / CLHEP::cm << " cm";
0116 if (doFineCalo_)
0117 edm::LogVerbatim("DoFineCalo") << "Using finecalo v2";
0118
0119
0120 edm::LogVerbatim("CaloSim") << "CaloSD: Have a possibility of " << fineNames.size() << " fine calorimeters of which "
0121 << useFines.size() << " are selected";
0122 for (unsigned int k = 0; k < fineNames.size(); ++k)
0123 edm::LogVerbatim("CaloSim") << "[" << k << "] " << fineNames[k] << " at " << fineLevels[k];
0124 std::ostringstream st1;
0125 for (unsigned int k = 0; k < useFines.size(); ++k)
0126 st1 << " [" << k << "] " << useFines[k] << ":" << fineNames[useFines[k]];
0127 edm::LogVerbatim("CaloSim") << "CaloSD used calorimeters" << st1.str();
0128 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0129 std::vector<G4LogicalVolume*>::const_iterator lvcite;
0130 for (unsigned int i = 0; i < useFines.size(); i++) {
0131 G4LogicalVolume* lv = nullptr;
0132 G4String name = static_cast<G4String>(fineNames[useFines[i]]);
0133 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0134 G4String namx(static_cast<std::string>(dd4hep::dd::noNamespace((*lvcite)->GetName())));
0135 if (namx == name) {
0136 lv = (*lvcite);
0137 break;
0138 }
0139 }
0140 if (lv != nullptr) {
0141 CaloSD::Detector detector;
0142 detector.name = name;
0143 detector.lv = lv;
0144 detector.level = fineLevels[useFines[i]];
0145 fineDetectors_.emplace_back(detector);
0146 }
0147 }
0148 #ifdef EDM_ML_DEBUG
0149 edm::LogVerbatim("CaloSim") << "CaloSD::Loads information for " << fineDetectors_.size() << " fine detectors";
0150 unsigned int k(0);
0151 for (const auto& detector : fineDetectors_) {
0152 edm::LogVerbatim("CaloSim") << "Detector[" << k << "] " << detector.name << " at level " << detector.level
0153 << " pointer to LV: " << detector.lv;
0154 }
0155 #endif
0156 }
0157
0158 CaloSD::~CaloSD() {}
0159
0160 G4bool CaloSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0161 NaNTrap(aStep);
0162 ignoreReject = false;
0163
0164 #ifdef EDM_ML_DEBUG
0165 edm::LogVerbatim("CaloSim") << "CaloSD::" << GetName() << " ID= " << aStep->GetTrack()->GetTrackID()
0166 << " prID= " << aStep->GetTrack()->GetParentID()
0167 << " Eprestep= " << aStep->GetPreStepPoint()->GetKineticEnergy()
0168 << " step= " << aStep->GetStepLength() << " Edep= " << aStep->GetTotalEnergyDeposit();
0169 #endif
0170
0171
0172 doFineCaloThisStep_ = (doFineCalo_ && isItFineCalo(aStep->GetPreStepPoint()->GetTouchable()));
0173
0174
0175
0176 if (isParameterized) {
0177 if (getFromLibrary(aStep)) {
0178
0179
0180 (aStep->GetTrack())->SetTrackStatus(fStopAndKill);
0181 if (0 < aStep->GetNumberOfSecondariesInCurrentStep()) {
0182 auto tv = aStep->GetSecondaryInCurrentStep();
0183 auto vol = aStep->GetPreStepPoint()->GetPhysicalVolume();
0184 for (auto& tk : *tv) {
0185 if (tk->GetVolume() == vol) {
0186 const_cast<G4Track*>(tk)->SetTrackStatus(fStopAndKill);
0187 }
0188 }
0189 }
0190 return true;
0191 }
0192 }
0193
0194
0195 edepositEM = edepositHAD = 0.f;
0196 if (aStep->GetTotalEnergyDeposit() <= 0.0) {
0197 return false;
0198 }
0199
0200
0201 unsigned int unitID = setDetUnitId(aStep);
0202 auto const theTrack = aStep->GetTrack();
0203 uint16_t depth = getDepth(aStep);
0204
0205 double time = theTrack->GetGlobalTime() / nanosecond;
0206 int primaryID = getTrackID(theTrack);
0207 if (unitID > 0) {
0208 currentID.setID(unitID, time, primaryID, depth);
0209 } else {
0210 if (!ignoreReject) {
0211 const G4TouchableHistory* touch = static_cast<const G4TouchableHistory*>(theTrack->GetTouchable());
0212 edm::LogVerbatim("CaloSim") << "CaloSD::ProcessHits: unitID= " << unitID << " currUnit= " << currentID.unitID()
0213 << " Detector: " << GetName() << " trackID= " << theTrack->GetTrackID() << " "
0214 << theTrack->GetDefinition()->GetParticleName()
0215 << "\n Edep= " << aStep->GetTotalEnergyDeposit()
0216 << " PV: " << touch->GetVolume(0)->GetName()
0217 << " PVid= " << touch->GetReplicaNumber(0) << " MVid= " << touch->GetReplicaNumber(1);
0218 }
0219 return false;
0220 }
0221 double energy = getEnergyDeposit(aStep);
0222 if (energy <= 0.0) {
0223 return false;
0224 }
0225
0226 if (doFineCaloThisStep_) {
0227 currentID.setID(unitID, time, findBoundaryCrossingParent(theTrack), depth);
0228 currentID.markAsFinecaloTrackID();
0229 }
0230
0231 if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0232 edepositEM = energy;
0233 } else {
0234 edepositHAD = energy;
0235 }
0236 #ifdef EDM_ML_DEBUG
0237 const G4TouchableHistory* touch = static_cast<const G4TouchableHistory*>(theTrack->GetTouchable());
0238 edm::LogVerbatim("CaloSim") << "CaloSD::" << GetName() << " PV:" << touch->GetVolume(0)->GetName()
0239 << " PVid=" << touch->GetReplicaNumber(0) << " MVid=" << touch->GetReplicaNumber(1)
0240 << " Unit:" << std::hex << unitID << std::dec << " Edep=" << edepositEM << " "
0241 << edepositHAD << " ID=" << theTrack->GetTrackID() << " pID=" << theTrack->GetParentID()
0242 << " E=" << theTrack->GetKineticEnergy() << " S=" << aStep->GetStepLength() << "\n "
0243 << theTrack->GetDefinition()->GetParticleName() << " primaryID= " << primaryID
0244 << " currentID= (" << currentID << ") previousID= (" << previousID << ")";
0245 #endif
0246 if (!hitExists(aStep)) {
0247 currentHit = createNewHit(aStep, theTrack);
0248 } else {
0249 #ifdef EDM_ML_DEBUG
0250 edm::LogVerbatim("DoFineCalo") << "Not creating new hit, only updating " << shortreprID(currentHit);
0251 #endif
0252 }
0253 return true;
0254 }
0255
0256 bool CaloSD::ProcessHits(G4GFlashSpot* aSpot, G4TouchableHistory*) {
0257 edepositEM = edepositHAD = 0.f;
0258 const G4Track* track = aSpot->GetOriginatorTrack()->GetPrimaryTrack();
0259
0260 double edep = aSpot->GetEnergySpot()->GetEnergy();
0261 if (edep <= 0.0) {
0262 return false;
0263 }
0264
0265 G4Step fFakeStep;
0266 G4StepPoint* fFakePreStepPoint = fFakeStep.GetPreStepPoint();
0267 G4StepPoint* fFakePostStepPoint = fFakeStep.GetPostStepPoint();
0268 fFakePreStepPoint->SetPosition(aSpot->GetPosition());
0269 fFakePostStepPoint->SetPosition(aSpot->GetPosition());
0270
0271 G4TouchableHandle fTouchableHandle = aSpot->GetTouchableHandle();
0272 fFakePreStepPoint->SetTouchableHandle(fTouchableHandle);
0273 fFakeStep.SetTotalEnergyDeposit(edep);
0274 edep = EnergyCorrected(fFakeStep, track);
0275
0276
0277 if (edep <= 0.0) {
0278 return false;
0279 }
0280
0281 if (G4TrackToParticleID::isGammaElectronPositron(track)) {
0282 edepositEM = edep;
0283 } else {
0284 edepositHAD = edep;
0285 }
0286
0287 unsigned int unitID = setDetUnitId(&fFakeStep);
0288
0289 if (unitID > 0) {
0290
0291 double time = track->GetGlobalTime() / nanosecond;
0292 int primaryID = getTrackID(track);
0293 uint16_t depth = getDepth(&fFakeStep);
0294 currentID.setID(unitID, time, primaryID, depth);
0295 #ifdef EDM_ML_DEBUG
0296 edm::LogVerbatim("CaloSim") << "CaloSD:: GetSpotInfo for Unit 0x" << std::hex << currentID.unitID() << std::dec
0297 << " Edeposit = " << edepositEM << " " << edepositHAD;
0298 #endif
0299
0300 if (currentID == previousID) {
0301 updateHit(currentHit);
0302 } else {
0303 posGlobal = aSpot->GetEnergySpot()->GetPosition();
0304
0305 if (currentID.trackID() != previousID.trackID()) {
0306 entrancePoint = aSpot->GetPosition();
0307 entranceLocal = aSpot->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(entrancePoint);
0308 incidentEnergy = track->GetKineticEnergy();
0309 #ifdef EDM_ML_DEBUG
0310 edm::LogVerbatim("CaloSim") << "CaloSD: Incident energy " << incidentEnergy / CLHEP::GeV << " GeV and"
0311 << " entrance point " << entrancePoint << " (Global) " << entranceLocal
0312 << " (Local)";
0313 #endif
0314 }
0315 if (!checkHit()) {
0316 currentHit = createNewHit(&fFakeStep, track);
0317 }
0318 }
0319 return true;
0320 }
0321 return false;
0322 }
0323
0324 double CaloSD::getEnergyDeposit(const G4Step* aStep) { return aStep->GetTotalEnergyDeposit(); }
0325
0326 double CaloSD::EnergyCorrected(const G4Step& aStep, const G4Track*) { return aStep.GetTotalEnergyDeposit(); }
0327
0328 bool CaloSD::getFromLibrary(const G4Step*) { return false; }
0329
0330 bool CaloSD::isItFineCalo(const G4VTouchable* touch) {
0331 bool ok(false);
0332 int level = ((touch->GetHistoryDepth()) + 1);
0333 for (const auto& detector : fineDetectors_) {
0334 if (level > 0 && level >= detector.level) {
0335 int ii = level - detector.level;
0336 G4LogicalVolume* lv = touch->GetVolume(ii)->GetLogicalVolume();
0337 ok = (lv == detector.lv);
0338 #ifdef EDM_ML_DEBUG
0339 std::string name1 = (lv == 0) ? "Unknown" : lv->GetName();
0340 edm::LogVerbatim("CaloSim") << "CaloSD: volume " << name1 << ":" << detector.name << " at Level "
0341 << detector.level << " Flag " << ok;
0342 #endif
0343 if (ok)
0344 break;
0345 }
0346 }
0347 return ok;
0348 }
0349
0350 void CaloSD::Initialize(G4HCofThisEvent* HCE) {
0351 totalHits = 0;
0352
0353 #ifdef EDM_ML_DEBUG
0354 edm::LogVerbatim("CaloSim") << "CaloSD : Initialize called for " << GetName();
0355 #endif
0356
0357
0358
0359 theHC = new CaloG4HitCollection(GetName(), collectionName[0]);
0360
0361 if (hcID < 0) {
0362 hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0363 }
0364
0365 HCE->AddHitsCollection(hcID, theHC);
0366 }
0367
0368 void CaloSD::EndOfEvent(G4HCofThisEvent*) {
0369
0370
0371 cleanHitCollection();
0372
0373 #ifdef EDM_ML_DEBUG
0374 if (theHC == nullptr)
0375 edm::LogVerbatim("CaloSim") << "CaloSD: EndofEvent entered with no entries";
0376 else
0377 edm::LogVerbatim("CaloSim") << "CaloSD: EndofEvent entered with " << theHC->entries() << " entries";
0378 #endif
0379 }
0380
0381 void CaloSD::clear() {}
0382
0383 void CaloSD::DrawAll() {}
0384
0385 void CaloSD::PrintAll() {
0386 #ifdef EDM_ML_DEBUG
0387 edm::LogVerbatim("CaloSim") << "CaloSD: Collection " << theHC->GetName();
0388 #endif
0389 theHC->PrintAllHits();
0390 }
0391
0392 void CaloSD::fillHits(edm::PCaloHitContainer& cc, const std::string& hname) {
0393 #ifdef EDM_ML_DEBUG
0394 edm::LogVerbatim("CaloSim") << "CaloSD: Tries to transfer " << slave.get()->hits().size() << " hits for "
0395 << slave.get()->name() << " " << hname;
0396 #endif
0397 if (slave.get()->name() == hname) {
0398 cc = slave.get()->hits();
0399 }
0400 slave.get()->Clean();
0401 }
0402
0403 G4ThreeVector CaloSD::setToLocal(const G4ThreeVector& global, const G4VTouchable* touch) const {
0404 return touch->GetHistory()->GetTopTransform().TransformPoint(global);
0405 }
0406
0407 G4ThreeVector CaloSD::setToGlobal(const G4ThreeVector& local, const G4VTouchable* touch) const {
0408 return touch->GetHistory()->GetTopTransform().Inverse().TransformPoint(local);
0409 }
0410
0411 bool CaloSD::hitExists(const G4Step* aStep) {
0412
0413 if (currentID == previousID) {
0414 updateHit(currentHit);
0415 return true;
0416 }
0417
0418
0419
0420
0421
0422
0423
0424
0425 posGlobal = aStep->GetPreStepPoint()->GetPosition();
0426 if (currentID.trackID() != previousID.trackID()) {
0427 resetForNewPrimary(aStep);
0428 }
0429 return checkHit();
0430 }
0431
0432 bool CaloSD::checkHit() {
0433
0434 bool found = false;
0435 if (useMap) {
0436 std::map<CaloHitID, CaloG4Hit*>::const_iterator it = hitMap.find(currentID);
0437 if (it != hitMap.end()) {
0438 currentHit = it->second;
0439 found = true;
0440 }
0441 } else if (nCheckedHits > 0) {
0442 int nhits = theHC->entries();
0443 int minhit = std::max(nhits - nCheckedHits, 0);
0444 int maxhit = nhits - 1;
0445
0446 for (int j = maxhit; j > minhit; --j) {
0447 if ((*theHC)[j]->getID() == currentID) {
0448 currentHit = (*theHC)[j];
0449 found = true;
0450 break;
0451 }
0452 }
0453 }
0454
0455 if (found) {
0456 updateHit(currentHit);
0457 }
0458 return found;
0459 }
0460
0461 int CaloSD::getNumberOfHits() { return theHC->entries(); }
0462
0463
0464
0465
0466
0467 std::string CaloSD::printableDecayChain(const std::vector<unsigned int>& decayChain) {
0468 std::stringstream ss;
0469 for (long unsigned int i = 0; i < decayChain.size(); i++) {
0470 if (i > 0)
0471 ss << " <- ";
0472 ss << decayChain[i];
0473 }
0474 return ss.str();
0475 }
0476
0477
0478 std::string CaloSD::shortreprID(const CaloHitID& ID) {
0479 std::stringstream ss;
0480 ss << GetName() << "/" << ID.unitID() << "/trk" << ID.trackID() << "/d" << ID.depth() << "/time" << ID.timeSliceID();
0481 if (ID.isFinecaloTrackID())
0482 ss << "/FC";
0483 return ss.str();
0484 }
0485
0486
0487 std::string CaloSD::shortreprID(const CaloG4Hit* hit) { return shortreprID(hit->getID()); }
0488
0489
0490
0491
0492 unsigned int CaloSD::findBoundaryCrossingParent(const G4Track* track, bool markAsSaveable) {
0493 TrackInformation* trkInfo = cmsTrackInformation(track);
0494 unsigned int id = track->GetTrackID();
0495
0496 auto it = boundaryCrossingParentMap_.find(id);
0497 if (it != boundaryCrossingParentMap_.end()) {
0498 #ifdef EDM_ML_DEBUG
0499 edm::LogVerbatim("DoFineCalo") << "Track " << id << " parent already cached: " << it->second;
0500 #endif
0501 return it->second;
0502 }
0503
0504 else if (trkInfo->crossedBoundary()) {
0505 #ifdef EDM_ML_DEBUG
0506 edm::LogVerbatim("DoFineCalo") << "Track " << id << " crosses boundary itself";
0507 #endif
0508 boundaryCrossingParentMap_[id] = id;
0509 trkInfo->setStoreTrack();
0510 return id;
0511 }
0512
0513 std::vector<unsigned int> decayChain{id};
0514 #ifdef EDM_ML_DEBUG
0515 edm::LogVerbatim("DoFineCalo") << "Track " << id << ": Traversing history to find boundary-crossing parent";
0516 #endif
0517 unsigned int parentID = track->GetParentID();
0518 while (true) {
0519 if (parentID == 0)
0520 throw cms::Exception("Unknown", "CaloSD")
0521 << "Hit end of parentage for track " << id << " without finding a boundary-crossing parent";
0522
0523 auto it = boundaryCrossingParentMap_.find(parentID);
0524 if (it != boundaryCrossingParentMap_.end()) {
0525 #ifdef EDM_ML_DEBUG
0526 edm::LogVerbatim("DoFineCalo") << " Track " << parentID
0527 << " boundary-crossing parent already cached: " << it->second;
0528 #endif
0529
0530 for (auto ancestorID : decayChain)
0531 boundaryCrossingParentMap_[ancestorID] = it->second;
0532 #ifdef EDM_ML_DEBUG
0533
0534 decayChain.push_back(parentID);
0535 while (parentID != it->second) {
0536 parentID = m_trackManager->getTrackByID(parentID, true)->parentID();
0537 decayChain.push_back(parentID);
0538 }
0539 edm::LogVerbatim("DoFineCalo") << " Full decay chain: " << printableDecayChain(decayChain);
0540 #endif
0541 return it->second;
0542 }
0543
0544 TrackWithHistory* parentTrack = m_trackManager->getTrackByID(parentID, true);
0545 if (parentTrack->crossedBoundary()) {
0546 if (markAsSaveable)
0547 parentTrack->setToBeSaved();
0548 decayChain.push_back(parentID);
0549
0550 for (auto ancestorID : decayChain)
0551 boundaryCrossingParentMap_[ancestorID] = parentID;
0552 #ifdef EDM_ML_DEBUG
0553 edm::LogVerbatim("DoFineCalo") << " Found boundary-crossing ancestor " << parentID << " for track " << id
0554 << "; decay chain: " << printableDecayChain(decayChain);
0555 #endif
0556 return parentID;
0557 }
0558
0559 decayChain.push_back(parentID);
0560 parentID = parentTrack->parentID();
0561 }
0562 }
0563
0564 CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) {
0565 #ifdef EDM_ML_DEBUG
0566 edm::LogVerbatim("CaloSim") << "CaloSD::CreateNewHit " << getNumberOfHits() << " for " << GetName()
0567 << " Unit:" << currentID.unitID() << " " << currentID.depth() << " Edep= " << edepositEM
0568 << " " << edepositHAD << " primaryID= " << currentID.trackID()
0569 << " timeSlice= " << currentID.timeSliceID() << " ID= " << theTrack->GetTrackID() << " "
0570 << theTrack->GetDefinition()->GetParticleName()
0571 << " E(GeV)= " << theTrack->GetKineticEnergy() / CLHEP::GeV
0572 << " parentID= " << theTrack->GetParentID() << "\n Ein= " << incidentEnergy
0573 << " entranceGlobal: " << entrancePoint << " entranceLocal: " << entranceLocal
0574 << " posGlobal: " << posGlobal;
0575 #endif
0576
0577 CaloG4Hit* aHit;
0578 if (!reusehit.empty()) {
0579 aHit = reusehit.back().release();
0580 aHit->setEM(0.f);
0581 aHit->setHadr(0.f);
0582 reusehit.pop_back();
0583 } else {
0584 aHit = new CaloG4Hit;
0585 }
0586
0587 aHit->setID(currentID);
0588 aHit->setEntry(entrancePoint.x(), entrancePoint.y(), entrancePoint.z());
0589 aHit->setEntryLocal(entranceLocal.x(), entranceLocal.y(), entranceLocal.z());
0590 aHit->setPosition(posGlobal.x(), posGlobal.y(), posGlobal.z());
0591 aHit->setIncidentEnergy(incidentEnergy);
0592 updateHit(aHit);
0593
0594 storeHit(aHit);
0595 TrackInformation* trkInfo = cmsTrackInformation(theTrack);
0596
0597 #ifdef EDM_ML_DEBUG
0598 if (doFineCaloThisStep_)
0599 edm::LogVerbatim("DoFineCalo") << "New hit " << shortreprID(aHit) << " using finecalo;"
0600 << " isItFineCalo(post)=" << isItFineCalo(aStep->GetPostStepPoint()->GetTouchable())
0601 << " isItFineCalo(pre)=" << isItFineCalo(aStep->GetPreStepPoint()->GetTouchable());
0602 #endif
0603
0604
0605 if (!doFineCaloThisStep_) {
0606 double etrack = 0;
0607 if (currentID.trackID() == primIDSaved) {
0608 } else if (currentID.trackID() == theTrack->GetTrackID()) {
0609 etrack = theTrack->GetKineticEnergy();
0610 #ifdef EDM_ML_DEBUG
0611 edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " etrack " << etrack
0612 << " eCut " << energyCut << " force: " << forceSave
0613 << " save: " << (etrack >= energyCut || forceSave);
0614 #endif
0615 if (etrack >= energyCut || forceSave) {
0616 trkInfo->setStoreTrack();
0617 }
0618 } else {
0619 TrackWithHistory* trkh = tkMap[currentID.trackID()];
0620 #ifdef EDM_ML_DEBUG
0621 edm::LogVerbatim("CaloSim") << "CaloSD : TrackWithHistory pointer for " << currentID.trackID() << " is " << trkh;
0622 #endif
0623 if (trkh != nullptr) {
0624 etrack = sqrt(trkh->momentum().Mag2());
0625 if (etrack >= energyCut) {
0626 trkh->setToBeSaved();
0627 #ifdef EDM_ML_DEBUG
0628 edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " with Hit";
0629 #endif
0630 }
0631 }
0632 }
0633 primIDSaved = currentID.trackID();
0634 }
0635
0636 if (useMap)
0637 ++totalHits;
0638 return aHit;
0639 }
0640
0641 void CaloSD::updateHit(CaloG4Hit* aHit) {
0642 aHit->addEnergyDeposit(edepositEM, edepositHAD);
0643 #ifdef EDM_ML_DEBUG
0644 edm::LogVerbatim("CaloSim") << "CaloSD:" << GetName() << " Add energy deposit in " << currentID
0645 << " Edep_em(MeV)= " << edepositEM << " Edep_had(MeV)= " << edepositHAD;
0646 #endif
0647
0648
0649 previousID = currentID;
0650 }
0651
0652 void CaloSD::resetForNewPrimary(const G4Step* aStep) {
0653 auto const preStepPoint = aStep->GetPreStepPoint();
0654 entrancePoint = preStepPoint->GetPosition();
0655 entranceLocal = setToLocal(entrancePoint, preStepPoint->GetTouchable());
0656 incidentEnergy = preStepPoint->GetKineticEnergy();
0657 #ifdef EDM_ML_DEBUG
0658 edm::LogVerbatim("CaloSim") << "CaloSD::resetForNewPrimary for " << GetName()
0659 << " ID= " << aStep->GetTrack()->GetTrackID() << " Ein= " << incidentEnergy / CLHEP::GeV
0660 << " GeV and"
0661 << " entrance point global: " << entrancePoint << " local: " << entranceLocal;
0662 #endif
0663 }
0664
0665 double CaloSD::getAttenuation(const G4Step* aStep, double birk1, double birk2, double birk3) const {
0666 double weight = 1.;
0667 double charge = aStep->GetPreStepPoint()->GetCharge();
0668 double length = aStep->GetStepLength();
0669
0670 if (charge != 0. && length > 0.) {
0671 double density = aStep->GetPreStepPoint()->GetMaterial()->GetDensity();
0672 double dedx = aStep->GetTotalEnergyDeposit() / length;
0673 double rkb = birk1 / density;
0674 double c = birk2 * rkb * rkb;
0675 if (std::abs(charge) >= 2.)
0676 rkb /= birk3;
0677 weight = 1. / (1. + rkb * dedx + c * dedx * dedx);
0678 #ifdef EDM_ML_DEBUG
0679 edm::LogVerbatim("CaloSim") << "CaloSD::getAttenuation in " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
0680 << " Charge " << charge << " dE/dx " << dedx << " Birk Const " << rkb << ", " << c
0681 << " Weight = " << weight << " dE " << aStep->GetTotalEnergyDeposit();
0682 #endif
0683 }
0684 return weight;
0685 }
0686
0687 void CaloSD::update(const BeginOfRun*) { initRun(); }
0688
0689 void CaloSD::update(const BeginOfEvent* g4Event) {
0690 #ifdef EDM_ML_DEBUG
0691 edm::LogVerbatim("CaloSim") << "CaloSD: Dispatched BeginOfEvent for " << GetName() << " !";
0692 #endif
0693 clearHits();
0694 initEvent(g4Event);
0695 }
0696
0697 void CaloSD::update(const EndOfTrack* trk) {
0698 int id = (*trk)()->GetTrackID();
0699 TrackInformation* trkI = cmsTrackInformation((*trk)());
0700 int lastTrackID = -1;
0701 if (trkI)
0702 lastTrackID = trkI->getIDonCaloSurface();
0703 if (id == lastTrackID) {
0704 auto trksForThisEvent = m_trackManager->trackContainer();
0705 if (!trksForThisEvent->empty()) {
0706 TrackWithHistory* trkH = trksForThisEvent->back();
0707 if (trkH->trackID() == (unsigned int)(id)) {
0708 tkMap[id] = trkH;
0709 #ifdef EDM_ML_DEBUG
0710 edm::LogVerbatim("CaloSim") << "CaloSD: get track " << id << " from Container of size "
0711 << trksForThisEvent->size() << " with ID " << trkH->trackID();
0712 #endif
0713 }
0714 }
0715 }
0716 }
0717
0718 void CaloSD::update(const ::EndOfEvent*) {
0719 endEvent();
0720 slave.get()->ReserveMemory(theHC->entries());
0721
0722 int count(0);
0723 double eEM(0.0);
0724 double eHAD(0.0);
0725 double eEM2(0.0);
0726 double eHAD2(0.0);
0727 #ifdef EDM_ML_DEBUG
0728 int wrong(0);
0729 double tt(0.0);
0730 double zloc(0.0);
0731 double zglob(0.0);
0732 double ee(0.0);
0733 #endif
0734 int hc_entries = theHC->entries();
0735 for (int i = 0; i < hc_entries; ++i) {
0736 #ifdef EDM_ML_DEBUG
0737 if (!saveHit((*theHC)[i])) {
0738 ++wrong;
0739 }
0740 #else
0741 saveHit((*theHC)[i]);
0742 #endif
0743
0744 ++count;
0745 double x = (*theHC)[i]->getEM();
0746 eEM += x;
0747 eEM2 += x * x;
0748 x = (*theHC)[i]->getHadr();
0749 eHAD += x;
0750 eHAD2 += x * x;
0751 #ifdef EDM_ML_DEBUG
0752 tt += (*theHC)[i]->getTimeSlice();
0753 ee += (*theHC)[i]->getIncidentEnergy();
0754 zglob += std::abs((*theHC)[i]->getEntry().z());
0755 zloc += std::abs((*theHC)[i]->getEntryLocal().z());
0756 #endif
0757 }
0758
0759 double norm = (count > 0) ? 1.0 / count : 0.0;
0760 eEM *= norm;
0761 eEM2 *= norm;
0762 eHAD *= norm;
0763 eHAD2 *= norm;
0764 eEM2 = std::sqrt(eEM2 - eEM * eEM);
0765 eHAD2 = std::sqrt(eHAD2 - eHAD * eHAD);
0766 #ifdef EDM_ML_DEBUG
0767 tt *= norm;
0768 ee *= norm;
0769 zglob *= norm;
0770 zloc *= norm;
0771 edm::LogVerbatim("CaloSim") << "CaloSD: " << GetName() << " store " << count << " hits; " << wrong
0772 << " track IDs not given properly and " << totalHits - count
0773 << " hits not passing cuts\n EmeanEM= " << eEM << " ErmsEM= " << eEM2
0774 << "\n EmeanHAD= " << eHAD << " ErmsHAD= " << eHAD2 << " TimeMean= " << tt
0775 << " E0mean= " << ee << " Zglob= " << zglob << " Zloc= " << zloc << " ";
0776 #endif
0777 tkMap.erase(tkMap.begin(), tkMap.end());
0778 std::vector<std::unique_ptr<CaloG4Hit>>().swap(reusehit);
0779 if (useMap)
0780 hitMap.erase(hitMap.begin(), hitMap.end());
0781 boundaryCrossingParentMap_.clear();
0782 }
0783
0784 void CaloSD::clearHits() {
0785 cleanIndex = 0;
0786 previousID.reset();
0787 primIDSaved = -99;
0788 #ifdef EDM_ML_DEBUG
0789 edm::LogVerbatim("CaloSim") << "CaloSD: Clears hit vector for " << GetName()
0790 << " and initialise slave: " << slave.get()->name();
0791 #endif
0792 slave.get()->Initialize();
0793 }
0794
0795 void CaloSD::reset() {
0796 if (fpCaloG4HitAllocator) {
0797 fpCaloG4HitAllocator->ResetStorage();
0798 }
0799 }
0800
0801 void CaloSD::initRun() {}
0802
0803 void CaloSD::initEvent(const BeginOfEvent*) {}
0804
0805 void CaloSD::endEvent() {}
0806
0807 int CaloSD::getTrackID(const G4Track* aTrack) {
0808 int primaryID = 0;
0809 TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0810 if (trkInfo) {
0811 primaryID = trkInfo->getIDonCaloSurface();
0812 #ifdef EDM_ML_DEBUG
0813 edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDonCaloSurface() << ":" << aTrack->GetTrackID() << ":"
0814 << primaryID;
0815 #endif
0816 } else {
0817 primaryID = aTrack->GetTrackID();
0818 #ifdef EDM_ML_DEBUG
0819 edm::LogWarning("CaloSim") << "CaloSD: Problem with primaryID **** set by force to TkID **** " << primaryID;
0820 #endif
0821 }
0822 return primaryID;
0823 }
0824
0825 int CaloSD::setTrackID(const G4Step* aStep) {
0826 auto const theTrack = aStep->GetTrack();
0827 TrackInformation* trkInfo = cmsTrackInformation(theTrack);
0828 int primaryID = trkInfo->getIDonCaloSurface();
0829 if (primaryID <= 0) {
0830 primaryID = theTrack->GetTrackID();
0831 }
0832 #ifdef EDM_ML_DEBUG
0833 edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDonCaloSurface() << ":" << theTrack->GetTrackID() << ":"
0834 << primaryID;
0835 #endif
0836
0837 if (primaryID != previousID.trackID()) {
0838 resetForNewPrimary(aStep);
0839 }
0840 #ifdef EDM_ML_DEBUG
0841 edm::LogVerbatim("CaloSim") << "CaloSD::setTrackID for " << GetName()
0842 << " trackID= " << aStep->GetTrack()->GetTrackID() << " primaryID= " << primaryID;
0843 #endif
0844 return primaryID;
0845 }
0846
0847 uint16_t CaloSD::getDepth(const G4Step*) { return 0; }
0848
0849 bool CaloSD::filterHit(CaloG4Hit* hit, double time) {
0850 double emin(eminHit);
0851 if (hit->getDepth() > 0)
0852 emin = eminHitD;
0853 #ifdef EDM_ML_DEBUG
0854 edm::LogVerbatim("CaloSim") << "CaloSD::filterHit(..) Depth " << hit->getDepth() << " Emin = " << emin << " ("
0855 << eminHit << ", " << eminHitD << ")";
0856 #endif
0857 return ((time <= tmaxHit) && (hit->getEnergyDeposit() > emin));
0858 }
0859
0860 double CaloSD::getResponseWt(const G4Track* aTrack) {
0861 double wt = 1.0;
0862 if (meanResponse.get()) {
0863 TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0864 wt = meanResponse.get()->getWeight(trkInfo->genParticlePID(), trkInfo->genParticleP());
0865 }
0866 return wt;
0867 }
0868
0869 void CaloSD::storeHit(CaloG4Hit* hit) {
0870 if (hit == nullptr || previousID.trackID() < 0) {
0871 edm::LogWarning("CaloSim") << "CaloSD: hit to be stored is nullptr !!"
0872 << " previousID.trackID()= " << previousID.trackID();
0873 return;
0874 }
0875
0876 theHC->insert(hit);
0877 if (useMap)
0878 hitMap.insert(std::pair<CaloHitID, CaloG4Hit*>(previousID, hit));
0879 }
0880
0881 bool CaloSD::saveHit(CaloG4Hit* aHit) {
0882 int tkID;
0883 bool ok = true;
0884
0885 double time = aHit->getTimeSlice();
0886 if (corrTOFBeam)
0887 time += correctT;
0888
0889
0890 if (doFineCalo_ && aHit->isFinecaloTrackID()) {
0891 #ifdef EDM_ML_DEBUG
0892 edm::LogVerbatim("DoFineCalo") << "Saving hit " << shortreprID(aHit);
0893 #endif
0894 if (!m_trackManager)
0895 throw cms::Exception("Unknown", "CaloSD") << "m_trackManager not set, needed for finecalo!";
0896 if (!m_trackManager->trackExists(aHit->getTrackID()))
0897 throw cms::Exception("Unknown", "CaloSD")
0898 << "Error on hit " << shortreprID(aHit) << ": Parent track not in track manager";
0899 slave.get()->processHits(aHit->getUnitID(),
0900 aHit->getEM() / CLHEP::GeV,
0901 aHit->getHadr() / CLHEP::GeV,
0902 time,
0903 aHit->getTrackID(),
0904 aHit->getDepth());
0905 }
0906
0907 else {
0908 if (m_trackManager) {
0909 tkID = m_trackManager->giveMotherNeeded(aHit->getTrackID());
0910 if (tkID == 0) {
0911 if (m_trackManager->trackExists(aHit->getTrackID()))
0912 tkID = (aHit->getTrackID());
0913 else {
0914 ok = false;
0915 }
0916 }
0917 } else {
0918 tkID = aHit->getTrackID();
0919 ok = false;
0920 }
0921 #ifdef EDM_ML_DEBUG
0922 edm::LogVerbatim("DoFineCalo") << "Saving hit " << shortreprID(aHit) << " with trackID=" << tkID;
0923 #endif
0924 slave.get()->processHits(
0925 aHit->getUnitID(), aHit->getEM() / CLHEP::GeV, aHit->getHadr() / CLHEP::GeV, time, tkID, aHit->getDepth());
0926 }
0927
0928 #ifdef EDM_ML_DEBUG
0929 if (!ok)
0930 edm::LogWarning("CaloSim") << "CaloSD:Cannot find track ID for " << aHit->getTrackID();
0931 edm::LogVerbatim("CaloSim") << "CalosD: Track ID " << aHit->getTrackID() << " changed to " << tkID
0932 << " by SimTrackManager Status " << ok;
0933 #endif
0934
0935 #ifdef EDM_ML_DEBUG
0936 edm::LogVerbatim("CaloSim") << "CaloSD: Store Hit at " << std::hex << aHit->getUnitID() << std::dec << " "
0937 << aHit->getDepth() << " due to " << tkID << " in time " << time << " of energy "
0938 << aHit->getEM() / CLHEP::GeV << " GeV (EM) and " << aHit->getHadr() / CLHEP::GeV
0939 << " GeV (Hadr)";
0940 #endif
0941 return ok;
0942 }
0943
0944 void CaloSD::update(const BeginOfTrack* trk) {
0945 int primary = -1;
0946 TrackInformation* trkInfo = cmsTrackInformation((*trk)());
0947 if (trkInfo->isPrimary())
0948 primary = (*trk)()->GetTrackID();
0949
0950 #ifdef EDM_ML_DEBUG
0951 edm::LogVerbatim("CaloSim") << "New track: isPrimary " << trkInfo->isPrimary() << " primary ID = " << primary
0952 << " primary ancestor ID " << primAncestor;
0953 #endif
0954
0955
0956
0957 if (primary > 0 && primary != primAncestor) {
0958 primAncestor = primary;
0959
0960
0961
0962 if (theHC->entries() > 0)
0963 cleanHitCollection();
0964 }
0965 }
0966
0967 void CaloSD::cleanHitCollection() {
0968 std::vector<CaloG4Hit*>* theCollection = theHC->GetVector();
0969
0970 #ifdef EDM_ML_DEBUG
0971 edm::LogVerbatim("CaloSim") << "CaloSD: collection before merging, size = " << theHC->entries();
0972 #endif
0973 if (reusehit.empty())
0974 reusehit.reserve(theHC->entries() - cleanIndex);
0975
0976
0977 if (!useMap) {
0978 std::vector<CaloG4Hit*> hitvec;
0979
0980 hitvec.swap(*theCollection);
0981 sort((hitvec.begin() + cleanIndex), hitvec.end(), CaloG4HitLess());
0982 #ifdef EDM_ML_DEBUG
0983 edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: sort hits in buffer starting from "
0984 << "element = " << cleanIndex;
0985 for (unsigned int i = 0; i < hitvec.size(); ++i) {
0986 if (hitvec[i] == nullptr)
0987 edm::LogVerbatim("CaloSim") << i << " has a null pointer";
0988 else
0989 edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i];
0990 }
0991 #endif
0992 CaloG4HitEqual equal;
0993 for (unsigned int i = cleanIndex; i < hitvec.size(); ++i) {
0994 int jump = 0;
0995 for (unsigned int j = i + 1; j < hitvec.size() && equal(hitvec[i], hitvec[j]); ++j) {
0996 ++jump;
0997
0998 (*hitvec[i]).addEnergyDeposit(*hitvec[j]);
0999 (*hitvec[j]).setEM(0.);
1000 (*hitvec[j]).setHadr(0.);
1001 reusehit.emplace_back(hitvec[j]);
1002 hitvec[j] = nullptr;
1003 }
1004 i += jump;
1005 }
1006 #ifdef EDM_ML_DEBUG
1007 edm::LogVerbatim("CaloSim") << "CaloSD: cleanHitCollection merge the hits in buffer ";
1008 for (unsigned int i = 0; i < hitvec.size(); ++i) {
1009 if (hitvec[i] == nullptr)
1010 edm::LogVerbatim("CaloSim") << i << " has a null pointer";
1011 else
1012 edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i];
1013 }
1014 #endif
1015
1016 hitvec.erase(
1017 std::stable_partition(hitvec.begin() + cleanIndex, hitvec.end(), [](CaloG4Hit* p) { return p != nullptr; }),
1018 hitvec.end());
1019 #ifdef EDM_ML_DEBUG
1020 edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: remove the merged hits in buffer,"
1021 << " new size = " << hitvec.size();
1022 #endif
1023 hitvec.swap(*theCollection);
1024 totalHits = theHC->entries();
1025 }
1026
1027 #ifdef EDM_ML_DEBUG
1028 edm::LogVerbatim("CaloSim") << "CaloSD: collection after merging, size= " << theHC->entries()
1029 << " Size of reusehit= " << reusehit.size()
1030 << "\n starting hit selection from index = " << cleanIndex;
1031 #endif
1032
1033 int addhit = 0;
1034 for (unsigned int i = cleanIndex; i < theCollection->size(); ++i) {
1035 CaloG4Hit* aHit((*theCollection)[i]);
1036
1037
1038
1039 double time = aHit->getTimeSlice();
1040 if (corrTOFBeam)
1041 time += correctT;
1042 if (!filterHit(aHit, time)) {
1043 #ifdef EDM_ML_DEBUG
1044 edm::LogVerbatim("CaloSim") << "CaloSD: dropped CaloG4Hit "
1045 << " " << *aHit;
1046 #endif
1047
1048
1049
1050 reusehit.emplace_back((*theCollection)[i]);
1051 (*theCollection)[i] = nullptr;
1052 ++addhit;
1053 }
1054 }
1055
1056 #ifdef EDM_ML_DEBUG
1057 edm::LogVerbatim("CaloSim") << "CaloSD: Size of reusehit after selection = " << reusehit.size()
1058 << " Number of added hit = " << addhit;
1059 #endif
1060 if (useMap) {
1061 if (addhit > 0) {
1062 int offset = reusehit.size() - addhit;
1063 for (int ii = addhit - 1; ii >= 0; --ii) {
1064 CaloHitID theID = reusehit[offset + ii]->getID();
1065 hitMap.erase(theID);
1066 }
1067 }
1068 }
1069
1070
1071 theCollection->erase(
1072 std::stable_partition(
1073 theCollection->begin() + cleanIndex, theCollection->end(), [](CaloG4Hit* p) { return p != nullptr; }),
1074 theCollection->end());
1075 #ifdef EDM_ML_DEBUG
1076 edm::LogVerbatim("CaloSim") << "CaloSD: hit collection after selection, size = " << theHC->entries();
1077 theHC->PrintAllHits();
1078 #endif
1079
1080 cleanIndex = theHC->entries();
1081 }
1082
1083 void CaloSD::printDetectorLevels(const G4VTouchable* touch) const {
1084
1085 int level = ((touch->GetHistoryDepth()) + 1);
1086 std::ostringstream st1;
1087 st1 << level << " Levels:";
1088 if (level > 0) {
1089 for (int ii = 0; ii < level; ii++) {
1090 int i = level - ii - 1;
1091 G4VPhysicalVolume* pv = touch->GetVolume(i);
1092 std::string name = (pv != nullptr) ? pv->GetName() : "Unknown";
1093 st1 << " " << name << ":" << touch->GetReplicaNumber(i);
1094 }
1095 }
1096 edm::LogVerbatim("CaloSim") << st1.str();
1097 }