Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-14 22:52:33 

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       (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   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   // 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->setStoreTrack();
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->setToBeSaved();
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->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   // buffer for next steps:
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;  // based on alpha particle data
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   // More strict bookkeeping for finecalo
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   // Regular, not-fine way:
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   // update the information if a different primary track ID
0956 
0957   if (primary > 0 && primary != primAncestor) {
0958     primAncestor = primary;
0959 
0960     // clean the hits information
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   // if no map used, merge before hits to have the save situation as a map
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         // merge j to i
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     //move all nullptr to end of list and then remove them
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     // selection
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       // create the list of hits to be reused
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   //move all nullptr to end of list and then remove them
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   //Print name and copy numbers
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!