Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-08 03:34:56 

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