Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-04 01:26:19

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: CaloSD.cc
0003 // Description: Sensitive Detector class for calorimeters
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 // #define EDM_ML_DEBUG
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   //Parameters
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   // Treat fine calorimeters
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   // Class variable to determine whether finecalo rules should apply for this step
0172   doFineCaloThisStep_ = (doFineCalo_ && isItFineCalo(aStep->GetPreStepPoint()->GetTouchable()));
0173 
0174   // apply shower library or parameterisation
0175   // independent on energy deposition at a step
0176   if (isParameterized) {
0177     if (getFromLibrary(aStep)) {
0178       // for parameterized showers the primary track should be killed
0179       // secondary tracks should be killed if they are in the same volume
0180       const_cast<G4Track*>(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   // ignore steps without energy deposit
0195   edepositEM = edepositHAD = 0.f;
0196   if (aStep->GetTotalEnergyDeposit() <= 0.0) {
0197     return false;
0198   }
0199 
0200   // check unitID
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   G4TouchableHistory* touch = (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   // zero edep means hit outside the calorimeter
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     // time of initial track
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     // Update if in the same detector, time-slice and for same track
0300     if (currentID == previousID) {
0301       updateHit(currentHit);
0302     } else {
0303       posGlobal = aSpot->GetEnergySpot()->GetPosition();
0304       // Reset entry point for new primary
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   //This initialization is performed at the beginning of an event
0358   //------------------------------------------------------------
0359   theHC = new CaloG4HitCollection(GetName(), collectionName[0]);
0360 
0361   if (hcID < 0) {
0362     hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0363   }
0364   //theHC ownership is transfered here to HCE
0365   HCE->AddHitsCollection(hcID, theHC);
0366 }
0367 
0368 void CaloSD::EndOfEvent(G4HCofThisEvent*) {
0369   // clean the hits for the last tracks
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   // Update if in the same detector, time-slice and for same track
0413   if (currentID == previousID) {
0414     updateHit(currentHit);
0415     return true;
0416   }
0417 
0418   // Note T. Klijnsma:
0419   // This is a rather strange place to set these class variables.
0420   // The code would be much more readable if all logic for determining
0421   // whether to update a hit or create a new hit is done in one place,
0422   // and only then perform the actual updating or creating of the hit.
0423 
0424   // Reset entry point for new primary
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   //look in the HitContainer whether a hit with the same ID already exists:
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 Takes a vector of ints (representing trackIDs), and returns a formatted string
0465 for debugging purposes
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 /* Very short representation of a CaloHitID */
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 /* As above, but with a hit as input */
0487 std::string CaloSD::shortreprID(const CaloG4Hit* hit) { return shortreprID(hit->getID()); }
0488 
0489 /*
0490 Finds the boundary-crossing parent of a track, and stores it in the CaloSD's map
0491 */
0492 unsigned int CaloSD::findBoundaryCrossingParent(const G4Track* track, bool markAsSaveable) {
0493   TrackInformation* trkInfo = cmsTrackInformation(track);
0494   unsigned int id = track->GetTrackID();
0495   // First see if this track is already in the map
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   // Then see if the track itself crosses the boundary
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->storeTrack(true);
0510     return id;
0511   }
0512   // Else, traverse the history of the track
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     // First check if this ancestor is already in the map
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       // Store this parent also for the rest of the traversed decay chain
0530       for (auto ancestorID : decayChain)
0531         boundaryCrossingParentMap_[ancestorID] = it->second;
0532 #ifdef EDM_ML_DEBUG
0533       // In debug mode, still build the rest of the decay chain for debugging
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     // If not, get this parent from the track manager (expensive)
0544     TrackWithHistory* parentTrack = m_trackManager->getTrackByID(parentID, true);
0545     if (parentTrack->crossedBoundary()) {
0546       if (markAsSaveable)
0547         parentTrack->save();
0548       decayChain.push_back(parentID);
0549       // Record this boundary crossing parent for all traversed ancestors
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     // Next iteration
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   // 'Traditional', non-fine history bookkeeping
0605   if (!doFineCaloThisStep_) {
0606     double etrack = 0;
0607     if (currentID.trackID() == primIDSaved) {  // The track is saved; nothing to be done
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->storeTrack(true);
0617         trkInfo->putInHistory();
0618       }
0619     } else {
0620       TrackWithHistory* trkh = tkMap[currentID.trackID()];
0621 #ifdef EDM_ML_DEBUG
0622       edm::LogVerbatim("CaloSim") << "CaloSD : TrackWithHistory pointer for " << currentID.trackID() << " is " << trkh;
0623 #endif
0624       if (trkh != nullptr) {
0625         etrack = sqrt(trkh->momentum().Mag2());
0626         if (etrack >= energyCut) {
0627           trkh->save();
0628 #ifdef EDM_ML_DEBUG
0629           edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " with Hit";
0630 #endif
0631         }
0632       }
0633     }
0634     primIDSaved = currentID.trackID();
0635   }
0636 
0637   if (useMap)
0638     ++totalHits;
0639   return aHit;
0640 }
0641 
0642 void CaloSD::updateHit(CaloG4Hit* aHit) {
0643   aHit->addEnergyDeposit(edepositEM, edepositHAD);
0644 #ifdef EDM_ML_DEBUG
0645   edm::LogVerbatim("CaloSim") << "CaloSD:" << GetName() << " Add energy deposit in " << currentID
0646                               << " Edep_em(MeV)= " << edepositEM << " Edep_had(MeV)= " << edepositHAD;
0647 #endif
0648 
0649   // buffer for next steps:
0650   previousID = currentID;
0651 }
0652 
0653 void CaloSD::resetForNewPrimary(const G4Step* aStep) {
0654   auto const preStepPoint = aStep->GetPreStepPoint();
0655   entrancePoint = preStepPoint->GetPosition();
0656   entranceLocal = setToLocal(entrancePoint, preStepPoint->GetTouchable());
0657   incidentEnergy = preStepPoint->GetKineticEnergy();
0658 #ifdef EDM_ML_DEBUG
0659   edm::LogVerbatim("CaloSim") << "CaloSD::resetForNewPrimary for " << GetName()
0660                               << " ID= " << aStep->GetTrack()->GetTrackID() << " Ein= " << incidentEnergy / CLHEP::GeV
0661                               << " GeV and"
0662                               << " entrance point global: " << entrancePoint << " local: " << entranceLocal;
0663 #endif
0664 }
0665 
0666 double CaloSD::getAttenuation(const G4Step* aStep, double birk1, double birk2, double birk3) const {
0667   double weight = 1.;
0668   double charge = aStep->GetPreStepPoint()->GetCharge();
0669   double length = aStep->GetStepLength();
0670 
0671   if (charge != 0. && length > 0.) {
0672     double density = aStep->GetPreStepPoint()->GetMaterial()->GetDensity();
0673     double dedx = aStep->GetTotalEnergyDeposit() / length;
0674     double rkb = birk1 / density;
0675     double c = birk2 * rkb * rkb;
0676     if (std::abs(charge) >= 2.)
0677       rkb /= birk3;  // based on alpha particle data
0678     weight = 1. / (1. + rkb * dedx + c * dedx * dedx);
0679 #ifdef EDM_ML_DEBUG
0680     edm::LogVerbatim("CaloSim") << "CaloSD::getAttenuation in " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
0681                                 << " Charge " << charge << " dE/dx " << dedx << " Birk Const " << rkb << ", " << c
0682                                 << " Weight = " << weight << " dE " << aStep->GetTotalEnergyDeposit();
0683 #endif
0684   }
0685   return weight;
0686 }
0687 
0688 void CaloSD::update(const BeginOfRun*) { initRun(); }
0689 
0690 void CaloSD::update(const BeginOfEvent* g4Event) {
0691 #ifdef EDM_ML_DEBUG
0692   edm::LogVerbatim("CaloSim") << "CaloSD: Dispatched BeginOfEvent for " << GetName() << " !";
0693 #endif
0694   clearHits();
0695   initEvent(g4Event);
0696 }
0697 
0698 void CaloSD::update(const EndOfTrack* trk) {
0699   int id = (*trk)()->GetTrackID();
0700   TrackInformation* trkI = cmsTrackInformation((*trk)());
0701   int lastTrackID = -1;
0702   if (trkI)
0703     lastTrackID = trkI->getIDonCaloSurface();
0704   if (id == lastTrackID) {
0705     const TrackContainer* trksForThisEvent = m_trackManager->trackContainer();
0706     if (trksForThisEvent != nullptr) {
0707       int it = (int)(trksForThisEvent->size()) - 1;
0708       if (it >= 0) {
0709         TrackWithHistory* trkH = (*trksForThisEvent)[it];
0710         if (trkH->trackID() == (unsigned int)(id))
0711           tkMap[id] = trkH;
0712 #ifdef EDM_ML_DEBUG
0713         edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size "
0714                                     << trksForThisEvent->size() << " with ID " << trkH->trackID();
0715       } else {
0716         edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size "
0717                                     << trksForThisEvent->size() << " with no ID";
0718 #endif
0719       }
0720     }
0721   }
0722 }
0723 
0724 void CaloSD::update(const ::EndOfEvent*) {
0725   endEvent();
0726   slave.get()->ReserveMemory(theHC->entries());
0727 
0728   int count(0);
0729   int wrong(0);
0730   double eEM(0.0);
0731   double eHAD(0.0);
0732   double eEM2(0.0);
0733   double eHAD2(0.0);
0734   double tt(0.0);
0735   double zloc(0.0);
0736   double zglob(0.0);
0737   double ee(0.0);
0738   int hc_entries = theHC->entries();
0739   for (int i = 0; i < hc_entries; ++i) {
0740     if (!saveHit((*theHC)[i])) {
0741       ++wrong;
0742     }
0743     ++count;
0744     double x = (*theHC)[i]->getEM();
0745     eEM += x;
0746     eEM2 += x * x;
0747     x = (*theHC)[i]->getHadr();
0748     eHAD += x;
0749     eHAD2 += x * x;
0750     tt += (*theHC)[i]->getTimeSlice();
0751     ee += (*theHC)[i]->getIncidentEnergy();
0752     zglob += std::abs((*theHC)[i]->getEntry().z());
0753     zloc += std::abs((*theHC)[i]->getEntryLocal().z());
0754   }
0755 
0756   double norm = (count > 0) ? 1.0 / count : 0.0;
0757   eEM *= norm;
0758   eEM2 *= norm;
0759   eHAD *= norm;
0760   eHAD2 *= norm;
0761   eEM2 = std::sqrt(eEM2 - eEM * eEM);
0762   eHAD2 = std::sqrt(eHAD2 - eHAD * eHAD);
0763   tt *= norm;
0764   ee *= norm;
0765   zglob *= norm;
0766   zloc *= norm;
0767 
0768 #ifdef EDM_ML_DEBUG
0769   edm::LogVerbatim("CaloSim") << "CaloSD: " << GetName() << " store " << count << " hits; " << wrong
0770                               << " track IDs not given properly and " << totalHits - count
0771                               << " hits not passing cuts\n EmeanEM= " << eEM << " ErmsEM= " << eEM2
0772                               << "\n EmeanHAD= " << eHAD << " ErmsHAD= " << eHAD2 << " TimeMean= " << tt
0773                               << " E0mean= " << ee << " Zglob= " << zglob << " Zloc= " << zloc << " ";
0774 #endif
0775   tkMap.erase(tkMap.begin(), tkMap.end());
0776   std::vector<std::unique_ptr<CaloG4Hit>>().swap(reusehit);
0777   if (useMap)
0778     hitMap.erase(hitMap.begin(), hitMap.end());
0779   boundaryCrossingParentMap_.clear();
0780 }
0781 
0782 void CaloSD::clearHits() {
0783   cleanIndex = 0;
0784   previousID.reset();
0785   primIDSaved = -99;
0786 #ifdef EDM_ML_DEBUG
0787   edm::LogVerbatim("CaloSim") << "CaloSD: Clears hit vector for " << GetName()
0788                               << " and initialise slave: " << slave.get()->name();
0789 #endif
0790   slave.get()->Initialize();
0791 }
0792 
0793 void CaloSD::reset() {
0794   if (fpCaloG4HitAllocator) {
0795     fpCaloG4HitAllocator->ResetStorage();
0796   }
0797 }
0798 
0799 void CaloSD::initRun() {}
0800 
0801 void CaloSD::initEvent(const BeginOfEvent*) {}
0802 
0803 void CaloSD::endEvent() {}
0804 
0805 int CaloSD::getTrackID(const G4Track* aTrack) {
0806   int primaryID = 0;
0807   TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0808   if (trkInfo) {
0809     primaryID = trkInfo->getIDonCaloSurface();
0810 #ifdef EDM_ML_DEBUG
0811     edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDonCaloSurface() << ":" << aTrack->GetTrackID() << ":"
0812                                 << primaryID;
0813 #endif
0814   } else {
0815     primaryID = aTrack->GetTrackID();
0816 #ifdef EDM_ML_DEBUG
0817     edm::LogWarning("CaloSim") << "CaloSD: Problem with primaryID **** set by force to TkID **** " << primaryID;
0818 #endif
0819   }
0820   return primaryID;
0821 }
0822 
0823 int CaloSD::setTrackID(const G4Step* aStep) {
0824   auto const theTrack = aStep->GetTrack();
0825   TrackInformation* trkInfo = cmsTrackInformation(theTrack);
0826   int primaryID = trkInfo->getIDonCaloSurface();
0827   if (primaryID <= 0) {
0828     primaryID = theTrack->GetTrackID();
0829   }
0830 #ifdef EDM_ML_DEBUG
0831   edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDonCaloSurface() << ":" << theTrack->GetTrackID() << ":"
0832                               << primaryID;
0833 #endif
0834 
0835   if (primaryID != previousID.trackID()) {
0836     resetForNewPrimary(aStep);
0837   }
0838 #ifdef EDM_ML_DEBUG
0839   edm::LogVerbatim("CaloSim") << "CaloSD::setTrackID for " << GetName()
0840                               << " trackID= " << aStep->GetTrack()->GetTrackID() << " primaryID= " << primaryID;
0841 #endif
0842   return primaryID;
0843 }
0844 
0845 uint16_t CaloSD::getDepth(const G4Step*) { return 0; }
0846 
0847 bool CaloSD::filterHit(CaloG4Hit* hit, double time) {
0848   double emin(eminHit);
0849   if (hit->getDepth() > 0)
0850     emin = eminHitD;
0851 #ifdef EDM_ML_DEBUG
0852   edm::LogVerbatim("CaloSim") << "CaloSD::filterHit(..) Depth " << hit->getDepth() << " Emin = " << emin << " ("
0853                               << eminHit << ", " << eminHitD << ")";
0854 #endif
0855   return ((time <= tmaxHit) && (hit->getEnergyDeposit() > emin));
0856 }
0857 
0858 double CaloSD::getResponseWt(const G4Track* aTrack) {
0859   double wt = 1.0;
0860   if (meanResponse.get()) {
0861     TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0862     wt = meanResponse.get()->getWeight(trkInfo->genParticlePID(), trkInfo->genParticleP());
0863   }
0864   return wt;
0865 }
0866 
0867 void CaloSD::storeHit(CaloG4Hit* hit) {
0868   if (hit == nullptr || previousID.trackID() < 0) {
0869     edm::LogWarning("CaloSim") << "CaloSD: hit to be stored is nullptr !!"
0870                                << " previousID.trackID()= " << previousID.trackID();
0871     return;
0872   }
0873 
0874   theHC->insert(hit);
0875   if (useMap)
0876     hitMap.insert(std::pair<CaloHitID, CaloG4Hit*>(previousID, hit));
0877 }
0878 
0879 bool CaloSD::saveHit(CaloG4Hit* aHit) {
0880   int tkID;
0881   bool ok = true;
0882 
0883   double time = aHit->getTimeSlice();
0884   if (corrTOFBeam)
0885     time += correctT;
0886 
0887   // More strict bookkeeping for finecalo
0888   if (doFineCalo_ && aHit->isFinecaloTrackID()) {
0889 #ifdef EDM_ML_DEBUG
0890     edm::LogVerbatim("DoFineCalo") << "Saving hit " << shortreprID(aHit);
0891 #endif
0892     if (!m_trackManager)
0893       throw cms::Exception("Unknown", "CaloSD") << "m_trackManager not set, needed for finecalo!";
0894     if (!m_trackManager->trackExists(aHit->getTrackID()))
0895       throw cms::Exception("Unknown", "CaloSD")
0896           << "Error on hit " << shortreprID(aHit) << ": Parent track not in track manager";
0897     slave.get()->processHits(aHit->getUnitID(),
0898                              aHit->getEM() / CLHEP::GeV,
0899                              aHit->getHadr() / CLHEP::GeV,
0900                              time,
0901                              aHit->getTrackID(),
0902                              aHit->getDepth());
0903   }
0904   // Regular, not-fine way:
0905   else {
0906     if (m_trackManager) {
0907       tkID = m_trackManager->giveMotherNeeded(aHit->getTrackID());
0908       if (tkID == 0) {
0909         if (m_trackManager->trackExists(aHit->getTrackID()))
0910           tkID = (aHit->getTrackID());
0911         else {
0912           ok = false;
0913         }
0914       }
0915     } else {
0916       tkID = aHit->getTrackID();
0917       ok = false;
0918     }
0919 #ifdef EDM_ML_DEBUG
0920     edm::LogVerbatim("DoFineCalo") << "Saving hit " << shortreprID(aHit) << " with trackID=" << tkID;
0921 #endif
0922     slave.get()->processHits(
0923         aHit->getUnitID(), aHit->getEM() / CLHEP::GeV, aHit->getHadr() / CLHEP::GeV, time, tkID, aHit->getDepth());
0924   }
0925 
0926 #ifdef EDM_ML_DEBUG
0927   if (!ok)
0928     edm::LogWarning("CaloSim") << "CaloSD:Cannot find track ID for " << aHit->getTrackID();
0929   edm::LogVerbatim("CaloSim") << "CalosD: Track ID " << aHit->getTrackID() << " changed to " << tkID
0930                               << " by SimTrackManager Status " << ok;
0931 #endif
0932 
0933 #ifdef EDM_ML_DEBUG
0934   edm::LogVerbatim("CaloSim") << "CaloSD: Store Hit at " << std::hex << aHit->getUnitID() << std::dec << " "
0935                               << aHit->getDepth() << " due to " << tkID << " in time " << time << " of energy "
0936                               << aHit->getEM() / CLHEP::GeV << " GeV (EM) and " << aHit->getHadr() / CLHEP::GeV
0937                               << " GeV (Hadr)";
0938 #endif
0939   return ok;
0940 }
0941 
0942 void CaloSD::update(const BeginOfTrack* trk) {
0943   int primary = -1;
0944   TrackInformation* trkInfo = cmsTrackInformation((*trk)());
0945   if (trkInfo->isPrimary())
0946     primary = (*trk)()->GetTrackID();
0947 
0948 #ifdef EDM_ML_DEBUG
0949   edm::LogVerbatim("CaloSim") << "New track: isPrimary " << trkInfo->isPrimary() << " primary ID = " << primary
0950                               << " primary ancestor ID " << primAncestor;
0951 #endif
0952 
0953   // update the information if a different primary track ID
0954 
0955   if (primary > 0 && primary != primAncestor) {
0956     primAncestor = primary;
0957 
0958     // clean the hits information
0959 
0960     if (theHC->entries() > 0)
0961       cleanHitCollection();
0962   }
0963 }
0964 
0965 void CaloSD::cleanHitCollection() {
0966   std::vector<CaloG4Hit*>* theCollection = theHC->GetVector();
0967 
0968 #ifdef EDM_ML_DEBUG
0969   edm::LogVerbatim("CaloSim") << "CaloSD: collection before merging, size = " << theHC->entries();
0970 #endif
0971   if (reusehit.empty())
0972     reusehit.reserve(theHC->entries() - cleanIndex);
0973 
0974   // if no map used, merge before hits to have the save situation as a map
0975   if (!useMap) {
0976     std::vector<CaloG4Hit*> hitvec;
0977 
0978     hitvec.swap(*theCollection);
0979     sort((hitvec.begin() + cleanIndex), hitvec.end(), CaloG4HitLess());
0980 #ifdef EDM_ML_DEBUG
0981     edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: sort hits in buffer starting from "
0982                                 << "element = " << cleanIndex;
0983     for (unsigned int i = 0; i < hitvec.size(); ++i) {
0984       if (hitvec[i] == nullptr)
0985         edm::LogVerbatim("CaloSim") << i << " has a null pointer";
0986       else
0987         edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i];
0988     }
0989 #endif
0990     CaloG4HitEqual equal;
0991     for (unsigned int i = cleanIndex; i < hitvec.size(); ++i) {
0992       int jump = 0;
0993       for (unsigned int j = i + 1; j < hitvec.size() && equal(hitvec[i], hitvec[j]); ++j) {
0994         ++jump;
0995         // merge j to i
0996         (*hitvec[i]).addEnergyDeposit(*hitvec[j]);
0997         (*hitvec[j]).setEM(0.);
0998         (*hitvec[j]).setHadr(0.);
0999         reusehit.emplace_back(hitvec[j]);
1000         hitvec[j] = nullptr;
1001       }
1002       i += jump;
1003     }
1004 #ifdef EDM_ML_DEBUG
1005     edm::LogVerbatim("CaloSim") << "CaloSD: cleanHitCollection merge the hits in buffer ";
1006     for (unsigned int i = 0; i < hitvec.size(); ++i) {
1007       if (hitvec[i] == nullptr)
1008         edm::LogVerbatim("CaloSim") << i << " has a null pointer";
1009       else
1010         edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i];
1011     }
1012 #endif
1013     //move all nullptr to end of list and then remove them
1014     hitvec.erase(
1015         std::stable_partition(hitvec.begin() + cleanIndex, hitvec.end(), [](CaloG4Hit* p) { return p != nullptr; }),
1016         hitvec.end());
1017 #ifdef EDM_ML_DEBUG
1018     edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: remove the merged hits in buffer,"
1019                                 << " new size = " << hitvec.size();
1020 #endif
1021     hitvec.swap(*theCollection);
1022     totalHits = theHC->entries();
1023   }
1024 
1025 #ifdef EDM_ML_DEBUG
1026   edm::LogVerbatim("CaloSim") << "CaloSD: collection after merging, size= " << theHC->entries()
1027                               << " Size of reusehit= " << reusehit.size()
1028                               << "\n      starting hit selection from index = " << cleanIndex;
1029 #endif
1030 
1031   int addhit = 0;
1032   for (unsigned int i = cleanIndex; i < theCollection->size(); ++i) {
1033     CaloG4Hit* aHit((*theCollection)[i]);
1034 
1035     // selection
1036 
1037     double time = aHit->getTimeSlice();
1038     if (corrTOFBeam)
1039       time += correctT;
1040     if (!filterHit(aHit, time)) {
1041 #ifdef EDM_ML_DEBUG
1042       edm::LogVerbatim("CaloSim") << "CaloSD: dropped CaloG4Hit "
1043                                   << " " << *aHit;
1044 #endif
1045 
1046       // create the list of hits to be reused
1047 
1048       reusehit.emplace_back((*theCollection)[i]);
1049       (*theCollection)[i] = nullptr;
1050       ++addhit;
1051     }
1052   }
1053 
1054 #ifdef EDM_ML_DEBUG
1055   edm::LogVerbatim("CaloSim") << "CaloSD: Size of reusehit after selection = " << reusehit.size()
1056                               << " Number of added hit = " << addhit;
1057 #endif
1058   if (useMap) {
1059     if (addhit > 0) {
1060       int offset = reusehit.size() - addhit;
1061       for (int ii = addhit - 1; ii >= 0; --ii) {
1062         CaloHitID theID = reusehit[offset + ii]->getID();
1063         hitMap.erase(theID);
1064       }
1065     }
1066   }
1067 
1068   //move all nullptr to end of list and then remove them
1069   theCollection->erase(
1070       std::stable_partition(
1071           theCollection->begin() + cleanIndex, theCollection->end(), [](CaloG4Hit* p) { return p != nullptr; }),
1072       theCollection->end());
1073 #ifdef EDM_ML_DEBUG
1074   edm::LogVerbatim("CaloSim") << "CaloSD: hit collection after selection, size = " << theHC->entries();
1075   theHC->PrintAllHits();
1076 #endif
1077 
1078   cleanIndex = theHC->entries();
1079 }
1080 
1081 void CaloSD::printDetectorLevels(const G4VTouchable* touch) const {
1082   //Print name and copy numbers
1083   int level = ((touch->GetHistoryDepth()) + 1);
1084   std::ostringstream st1;
1085   st1 << level << " Levels:";
1086   if (level > 0) {
1087     for (int ii = 0; ii < level; ii++) {
1088       int i = level - ii - 1;
1089       G4VPhysicalVolume* pv = touch->GetVolume(i);
1090       std::string name = (pv != nullptr) ? pv->GetName() : "Unknown";
1091       st1 << " " << name << ":" << touch->GetReplicaNumber(i);
1092     }
1093   }
1094   edm::LogVerbatim("CaloSim") << st1.str();
1095 }