Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-06 04:27:13 

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 #ifdef EDM_ML_DEBUG
0436     edm::LogVerbatim("CaloSim") << "CaloSD: Collection " << theHC[k]->GetName();
0437 #endif
0438     theHC[k]->PrintAllHits();
0439   }
0440 }
0441 
0442 void CaloSD::fillHits(edm::PCaloHitContainer& cc, const std::string& hname) {
0443   for (int k = 0; k < nHC_; ++k) {
0444 #ifdef EDM_ML_DEBUG
0445     edm::LogVerbatim("CaloSim") << "CaloSD: Tries to transfer " << slave[k].get()->hits().size() << " hits for "
0446                                 << slave[k].get()->name() << "   " << hname;
0447 #endif
0448     if (slave[k].get()->name() == hname) {
0449       cc = slave[k].get()->hits();
0450       slave[k].get()->Clean();
0451     }
0452   }
0453 }
0454 
0455 G4ThreeVector CaloSD::setToLocal(const G4ThreeVector& global, const G4VTouchable* touch) const {
0456   return touch->GetHistory()->GetTopTransform().TransformPoint(global);
0457 }
0458 
0459 G4ThreeVector CaloSD::setToGlobal(const G4ThreeVector& local, const G4VTouchable* touch) const {
0460   return touch->GetHistory()->GetTopTransform().Inverse().TransformPoint(local);
0461 }
0462 
0463 bool CaloSD::hitExists(const G4Step* aStep, int k) {
0464 #ifdef EDM_ML_DEBUG
0465   edm::LogVerbatim("CaloSim") << "CaloSD: hitExists for " << k;
0466 #endif
0467   // Update if in the same detector, time-slice and for same track
0468   if (currentID[k] == previousID[k]) {
0469     updateHit(currentHit[k], k);
0470     return true;
0471   }
0472 
0473   // Note T. Klijnsma:
0474   // This is a rather strange place to set these class variables.
0475   // The code would be much more readable if all logic for determining
0476   // whether to update a hit or create a new hit is done in one place,
0477   // and only then perform the actual updating or creating of the hit.
0478 
0479   // Reset entry point for new primary (only for the primary hit)
0480   if (k == 0) {
0481     posGlobal = aStep->GetPreStepPoint()->GetPosition();
0482     if (currentID[k].trackID() != previousID[k].trackID()) {
0483       resetForNewPrimary(aStep);
0484     }
0485   }
0486   return checkHit(k);
0487 }
0488 
0489 bool CaloSD::checkHit(int k) {
0490 #ifdef EDM_ML_DEBUG
0491   edm::LogVerbatim("CaloSim") << "CaloSD: checkHit for " << k << " for map " << useMap << ":" << &hitMap[k] << " Nhits "
0492                               << nCheckedHits[k] << " HC " << theHC[k] << " ID " << currentID[k];
0493 #endif
0494   //look in the HitContainer whether a hit with the same ID already exists:
0495   bool found = false;
0496   if (useMap) {
0497     std::map<CaloHitID, CaloG4Hit*>::const_iterator it = hitMap[k].find(currentID[k]);
0498     if (it != hitMap[k].end()) {
0499       currentHit[k] = it->second;
0500       found = true;
0501     }
0502   } else if (nCheckedHits[k] > 0) {
0503     int nhits = theHC[k]->entries();
0504     int minhit = std::max(nhits - nCheckedHits[k], 0);
0505     int maxhit = nhits - 1;
0506 
0507     for (int j = maxhit; j > minhit; --j) {
0508       if ((*theHC[k])[j]->getID() == currentID[k]) {
0509         currentHit[k] = (*theHC[k])[j];
0510         found = true;
0511         break;
0512       }
0513     }
0514   }
0515 
0516   if (found) {
0517     updateHit(currentHit[k], k);
0518   }
0519   return found;
0520 }
0521 
0522 int CaloSD::getNumberOfHits(int k) { return theHC[k]->entries(); }
0523 
0524 /*
0525 Takes a vector of ints (representing trackIDs), and returns a formatted string
0526 for debugging purposes
0527 */
0528 std::string CaloSD::printableDecayChain(const std::vector<unsigned int>& decayChain) {
0529   std::stringstream ss;
0530   for (long unsigned int i = 0; i < decayChain.size(); i++) {
0531     if (i > 0)
0532       ss << " <- ";
0533     ss << decayChain[i];
0534   }
0535   return ss.str();
0536 }
0537 
0538 /* Very short representation of a CaloHitID */
0539 std::string CaloSD::shortreprID(const CaloHitID& ID) {
0540   std::stringstream ss;
0541   ss << GetName() << "/" << ID.unitID() << "/trk" << ID.trackID() << "/d" << ID.depth() << "/time" << ID.timeSliceID();
0542   if (ID.isFinecaloTrackID())
0543     ss << "/FC";
0544   return ss.str();
0545 }
0546 
0547 /* As above, but with a hit as input */
0548 std::string CaloSD::shortreprID(const CaloG4Hit* hit) { return shortreprID(hit->getID()); }
0549 
0550 /*
0551 Finds the boundary-crossing parent of a track, and stores it in the CaloSD's map
0552 */
0553 unsigned int CaloSD::findBoundaryCrossingParent(const G4Track* track, bool markAsSaveable) {
0554   TrackInformation* trkInfo = cmsTrackInformation(track);
0555   unsigned int id = track->GetTrackID();
0556   // First see if this track is already in the map
0557   auto it = boundaryCrossingParentMap_.find(id);
0558   if (it != boundaryCrossingParentMap_.end()) {
0559 #ifdef EDM_ML_DEBUG
0560     edm::LogVerbatim("DoFineCalo") << "Track " << id << " parent already cached: " << it->second;
0561 #endif
0562     return it->second;
0563   }
0564   // Then see if the track itself crosses the boundary
0565   else if (trkInfo->crossedBoundary()) {
0566 #ifdef EDM_ML_DEBUG
0567     edm::LogVerbatim("DoFineCalo") << "Track " << id << " crosses boundary itself";
0568 #endif
0569     boundaryCrossingParentMap_[id] = id;
0570     trkInfo->setStoreTrack();
0571     return id;
0572   }
0573   // Else, traverse the history of the track
0574   std::vector<unsigned int> decayChain{id};
0575 #ifdef EDM_ML_DEBUG
0576   edm::LogVerbatim("DoFineCalo") << "Track " << id << ": Traversing history to find boundary-crossing parent";
0577 #endif
0578   unsigned int parentID = track->GetParentID();
0579   while (true) {
0580     if (parentID == 0)
0581       throw cms::Exception("Unknown", "CaloSD")
0582           << "Hit end of parentage for track " << id << " without finding a boundary-crossing parent";
0583     // First check if this ancestor is already in the map
0584     auto it = boundaryCrossingParentMap_.find(parentID);
0585     if (it != boundaryCrossingParentMap_.end()) {
0586 #ifdef EDM_ML_DEBUG
0587       edm::LogVerbatim("DoFineCalo") << "  Track " << parentID
0588                                      << " boundary-crossing parent already cached: " << it->second;
0589 #endif
0590       // Store this parent also for the rest of the traversed decay chain
0591       for (auto ancestorID : decayChain)
0592         boundaryCrossingParentMap_[ancestorID] = it->second;
0593 #ifdef EDM_ML_DEBUG
0594       // In debug mode, still build the rest of the decay chain for debugging
0595       decayChain.push_back(parentID);
0596       while (parentID != it->second) {
0597         parentID = m_trackManager->getTrackByID(parentID, true)->parentID();
0598         decayChain.push_back(parentID);
0599       }
0600       edm::LogVerbatim("DoFineCalo") << "  Full decay chain: " << printableDecayChain(decayChain);
0601 #endif
0602       return it->second;
0603     }
0604     // If not, get this parent from the track manager (expensive)
0605     TrackWithHistory* parentTrack = m_trackManager->getTrackByID(parentID, true);
0606     if (parentTrack->crossedBoundary()) {
0607       if (markAsSaveable)
0608         parentTrack->setToBeSaved();
0609       decayChain.push_back(parentID);
0610       // Record this boundary crossing parent for all traversed ancestors
0611       for (auto ancestorID : decayChain)
0612         boundaryCrossingParentMap_[ancestorID] = parentID;
0613 #ifdef EDM_ML_DEBUG
0614       edm::LogVerbatim("DoFineCalo") << "  Found boundary-crossing ancestor " << parentID << " for track " << id
0615                                      << "; decay chain: " << printableDecayChain(decayChain);
0616 #endif
0617       return parentID;
0618     }
0619     // Next iteration
0620     decayChain.push_back(parentID);
0621     parentID = parentTrack->parentID();
0622   }
0623 }
0624 
0625 CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack, int k) {
0626 #ifdef EDM_ML_DEBUG
0627   edm::LogVerbatim("CaloSim") << "CaloSD::CreateNewHit " << getNumberOfHits() << " for " << GetName()
0628                               << " Unit:" << currentID[k].unitID() << " " << currentID[k].depth()
0629                               << " Edep= " << edepositEM << " " << edepositHAD
0630                               << " primaryID= " << currentID[k].trackID()
0631                               << " timeSlice= " << currentID[k].timeSliceID() << " ID= " << theTrack->GetTrackID()
0632                               << " " << theTrack->GetDefinition()->GetParticleName()
0633                               << " E(GeV)= " << theTrack->GetKineticEnergy() / CLHEP::GeV
0634                               << " parentID= " << theTrack->GetParentID() << "\n Ein= " << incidentEnergy
0635                               << " entranceGlobal: " << entrancePoint << " entranceLocal: " << entranceLocal
0636                               << " posGlobal: " << posGlobal;
0637 #endif
0638 
0639   CaloG4Hit* aHit;
0640   if (!reusehit[k].empty()) {
0641     aHit = reusehit[k].back().release();
0642     aHit->setEM(0.f);
0643     aHit->setHadr(0.f);
0644     reusehit[k].pop_back();
0645   } else {
0646     aHit = new CaloG4Hit;
0647   }
0648 
0649   aHit->setID(currentID[k]);
0650   aHit->setEntry(entrancePoint.x(), entrancePoint.y(), entrancePoint.z());
0651   aHit->setEntryLocal(entranceLocal.x(), entranceLocal.y(), entranceLocal.z());
0652   aHit->setPosition(posGlobal.x(), posGlobal.y(), posGlobal.z());
0653   aHit->setIncidentEnergy(incidentEnergy);
0654   updateHit(aHit, k);
0655 
0656   storeHit(aHit, k);
0657   TrackInformation* trkInfo = cmsTrackInformation(theTrack);
0658 
0659 #ifdef EDM_ML_DEBUG
0660   if (doFineCaloThisStep_)
0661     edm::LogVerbatim("DoFineCalo") << "New hit " << shortreprID(aHit) << " using finecalo;"
0662                                    << " isItFineCalo(post)=" << isItFineCalo(aStep->GetPostStepPoint()->GetTouchable())
0663                                    << " isItFineCalo(pre)=" << isItFineCalo(aStep->GetPreStepPoint()->GetTouchable());
0664 #endif
0665 
0666   // 'Traditional', non-fine history bookkeeping
0667   if (!doFineCaloThisStep_) {
0668     double etrack = 0;
0669     if (currentID[k].trackID() == primIDSaved[k]) {  // The track is saved; nothing to be done
0670     } else if (currentID[k].trackID() == theTrack->GetTrackID()) {
0671       etrack = theTrack->GetKineticEnergy();
0672 #ifdef EDM_ML_DEBUG
0673       edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID[k].trackID() << " etrack " << etrack
0674                                   << " eCut " << energyCut << " force: " << forceSave
0675                                   << " save: " << (etrack >= energyCut || forceSave);
0676 #endif
0677       if (etrack >= energyCut || forceSave) {
0678         trkInfo->setStoreTrack();
0679       }
0680     } else {
0681       TrackWithHistory* trkh = tkMap[currentID[k].trackID()];
0682 #ifdef EDM_ML_DEBUG
0683       edm::LogVerbatim("CaloSim") << "CaloSD : TrackWithHistory pointer for " << currentID[k].trackID() << " is "
0684                                   << trkh;
0685 #endif
0686       if (trkh != nullptr) {
0687         etrack = sqrt(trkh->momentum().Mag2());
0688         if (etrack >= energyCut) {
0689           trkh->setToBeSaved();
0690 #ifdef EDM_ML_DEBUG
0691           edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID[k].trackID() << " with Hit";
0692 #endif
0693         }
0694       }
0695     }
0696     primIDSaved[k] = currentID[k].trackID();
0697   }
0698 
0699   if (useMap)
0700     ++totalHits[k];
0701   return aHit;
0702 }
0703 
0704 void CaloSD::updateHit(CaloG4Hit* aHit, int k) {
0705   aHit->addEnergyDeposit(edepositEM, edepositHAD);
0706 #ifdef EDM_ML_DEBUG
0707   edm::LogVerbatim("CaloSim") << "CaloSD:" << GetName() << " Add energy deposit in " << currentID[k]
0708                               << " Edep_em(MeV)= " << edepositEM << " Edep_had(MeV)= " << edepositHAD;
0709 #endif
0710 
0711   // buffer for next steps:
0712   previousID[k] = currentID[k];
0713 }
0714 
0715 void CaloSD::resetForNewPrimary(const G4Step* aStep) {
0716   auto const preStepPoint = aStep->GetPreStepPoint();
0717   entrancePoint = preStepPoint->GetPosition();
0718   entranceLocal = setToLocal(entrancePoint, preStepPoint->GetTouchable());
0719   incidentEnergy = preStepPoint->GetKineticEnergy();
0720 #ifdef EDM_ML_DEBUG
0721   edm::LogVerbatim("CaloSim") << "CaloSD::resetForNewPrimary for " << GetName()
0722                               << " ID= " << aStep->GetTrack()->GetTrackID() << " Ein= " << incidentEnergy / CLHEP::GeV
0723                               << " GeV and"
0724                               << " entrance point global: " << entrancePoint << " local: " << entranceLocal;
0725 #endif
0726 }
0727 
0728 double CaloSD::getAttenuation(const G4Step* aStep, double birk1, double birk2, double birk3) const {
0729   double weight = 1.;
0730   double charge = aStep->GetPreStepPoint()->GetCharge();
0731   double length = aStep->GetStepLength();
0732 
0733   if (charge != 0. && length > 0.) {
0734     double density = aStep->GetPreStepPoint()->GetMaterial()->GetDensity();
0735     double dedx = aStep->GetTotalEnergyDeposit() / length;
0736     double rkb = birk1 / density;
0737     double c = birk2 * rkb * rkb;
0738     if (std::abs(charge) >= 2.)
0739       rkb /= birk3;  // based on alpha particle data
0740     weight = 1. / (1. + rkb * dedx + c * dedx * dedx);
0741 #ifdef EDM_ML_DEBUG
0742     edm::LogVerbatim("CaloSim") << "CaloSD::getAttenuation in " << aStep->GetPreStepPoint()->GetMaterial()->GetName()
0743                                 << " Charge " << charge << " dE/dx " << dedx << " Birk Const " << rkb << ", " << c
0744                                 << " Weight = " << weight << " dE " << aStep->GetTotalEnergyDeposit();
0745 #endif
0746   }
0747   return weight;
0748 }
0749 
0750 void CaloSD::update(const BeginOfRun*) { initRun(); }
0751 
0752 void CaloSD::update(const BeginOfEvent* g4Event) {
0753 #ifdef EDM_ML_DEBUG
0754   edm::LogVerbatim("CaloSim") << "CaloSD: Dispatched BeginOfEvent for " << GetName() << " !";
0755 #endif
0756   clearHits();
0757   initEvent(g4Event);
0758 }
0759 
0760 void CaloSD::update(const EndOfTrack* trk) {
0761   int id = (*trk)()->GetTrackID();
0762   TrackInformation* trkI = cmsTrackInformation((*trk)());
0763   int lastTrackID = -1;
0764   if (trkI)
0765     lastTrackID = trkI->getIDonCaloSurface();
0766   if (id == lastTrackID) {
0767     auto trksForThisEvent = m_trackManager->trackContainer();
0768     if (!trksForThisEvent->empty()) {
0769       TrackWithHistory* trkH = trksForThisEvent->back();
0770       if (trkH->trackID() == id) {
0771         tkMap[id] = trkH;
0772 #ifdef EDM_ML_DEBUG
0773         edm::LogVerbatim("CaloSim") << "CaloSD: get track " << id << " from Container of size "
0774                                     << trksForThisEvent->size() << " with ID " << trkH->trackID();
0775 #endif
0776       }
0777     }
0778   }
0779 }
0780 
0781 void CaloSD::update(const ::EndOfEvent*) {
0782   endEvent();
0783   for (int k = 0; k < nHC_; ++k) {
0784     slave[k].get()->ReserveMemory(theHC[k]->entries());
0785 
0786     int count(0);
0787     double eEM(0.0);
0788     double eHAD(0.0);
0789     double eEM2(0.0);
0790     double eHAD2(0.0);
0791 #ifdef EDM_ML_DEBUG
0792     int wrong(0);
0793     double tt(0.0);
0794     double zloc(0.0);
0795     double zglob(0.0);
0796     double ee(0.0);
0797 #endif
0798     int hc_entries = theHC[k]->entries();
0799     for (int i = 0; i < hc_entries; ++i) {
0800 #ifdef EDM_ML_DEBUG
0801       if (!saveHit((*theHC[k])[i], k)) {
0802         ++wrong;
0803       }
0804 #else
0805       saveHit((*theHC[k])[i], k);
0806 #endif
0807 
0808       ++count;
0809       double x = (*theHC[k])[i]->getEM();
0810       eEM += x;
0811       eEM2 += x * x;
0812       x = (*theHC[k])[i]->getHadr();
0813       eHAD += x;
0814       eHAD2 += x * x;
0815 #ifdef EDM_ML_DEBUG
0816       tt += (*theHC[k])[i]->getTimeSlice();
0817       ee += (*theHC[k])[i]->getIncidentEnergy();
0818       zglob += std::abs((*theHC[k])[i]->getEntry().z());
0819       zloc += std::abs((*theHC[k])[i]->getEntryLocal().z());
0820 #endif
0821     }
0822 
0823     double norm = (count > 0) ? 1.0 / count : 0.0;
0824     eEM *= norm;
0825     eEM2 *= norm;
0826     eHAD *= norm;
0827     eHAD2 *= norm;
0828     eEM2 = std::sqrt(eEM2 - eEM * eEM);
0829     eHAD2 = std::sqrt(eHAD2 - eHAD * eHAD);
0830 #ifdef EDM_ML_DEBUG
0831     tt *= norm;
0832     ee *= norm;
0833     zglob *= norm;
0834     zloc *= norm;
0835     edm::LogVerbatim("CaloSim") << "CaloSD: " << GetName() << " store " << count << " hits; " << wrong
0836                                 << " track IDs not given properly and " << totalHits - count
0837                                 << " hits not passing cuts\n EmeanEM= " << eEM << " ErmsEM= " << eEM2
0838                                 << "\n EmeanHAD= " << eHAD << " ErmsHAD= " << eHAD2 << " TimeMean= " << tt
0839                                 << " E0mean= " << ee << " Zglob= " << zglob << " Zloc= " << zloc << " ";
0840 #endif
0841     std::vector<std::unique_ptr<CaloG4Hit>>().swap(reusehit[k]);
0842     if (useMap)
0843       hitMap[k].erase(hitMap[k].begin(), hitMap[k].end());
0844   }
0845   tkMap.erase(tkMap.begin(), tkMap.end());
0846   boundaryCrossingParentMap_.clear();
0847 }
0848 
0849 void CaloSD::clearHits() {
0850   for (int k = 0; k < nHC_; ++k) {
0851     previousID[k].reset();
0852     primIDSaved[k] = -99;
0853     cleanIndex[k] = 0;
0854 #ifdef EDM_ML_DEBUG
0855     edm::LogVerbatim("CaloSim") << "CaloSD: Clears hit vector for " << GetName()
0856                                 << " and initialise slave: " << slave[k].get()->name();
0857 #endif
0858     slave[k].get()->Initialize();
0859   }
0860 }
0861 
0862 void CaloSD::reset() {
0863   if (fpCaloG4HitAllocator) {
0864     fpCaloG4HitAllocator->ResetStorage();
0865   }
0866 }
0867 
0868 void CaloSD::initRun() {}
0869 
0870 void CaloSD::initEvent(const BeginOfEvent*) {}
0871 
0872 void CaloSD::endEvent() {}
0873 
0874 int CaloSD::getTrackID(const G4Track* aTrack) {
0875   int primaryID = 0;
0876   TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0877   if (trkInfo) {
0878     primaryID = trkInfo->getIDonCaloSurface();
0879 #ifdef EDM_ML_DEBUG
0880     edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDonCaloSurface() << ":" << aTrack->GetTrackID() << ":"
0881                                 << primaryID;
0882 #endif
0883   } else {
0884     primaryID = aTrack->GetTrackID();
0885 #ifdef EDM_ML_DEBUG
0886     edm::LogWarning("CaloSim") << "CaloSD: Problem with primaryID **** set by force to TkID **** " << primaryID;
0887 #endif
0888   }
0889   return primaryID;
0890 }
0891 
0892 int CaloSD::setTrackID(const G4Step* aStep) {
0893   auto const theTrack = aStep->GetTrack();
0894   TrackInformation* trkInfo = cmsTrackInformation(theTrack);
0895   int primaryID = trkInfo->getIDonCaloSurface();
0896   if (primaryID <= 0) {
0897     primaryID = theTrack->GetTrackID();
0898   }
0899 #ifdef EDM_ML_DEBUG
0900   edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDonCaloSurface() << ":" << theTrack->GetTrackID() << ":"
0901                               << primaryID;
0902 #endif
0903 
0904   if (primaryID != previousID[0].trackID()) {
0905     resetForNewPrimary(aStep);
0906   }
0907 #ifdef EDM_ML_DEBUG
0908   edm::LogVerbatim("CaloSim") << "CaloSD::setTrackID for " << GetName()
0909                               << " trackID= " << aStep->GetTrack()->GetTrackID() << " primaryID= " << primaryID;
0910 #endif
0911   return primaryID;
0912 }
0913 
0914 uint16_t CaloSD::getDepth(const G4Step*) { return 0; }
0915 
0916 void CaloSD::processSecondHit(const G4Step*, const G4Track*) {}
0917 
0918 bool CaloSD::filterHit(CaloG4Hit* hit, double time) {
0919   double emin(eminHit);
0920   if (hit->getDepth() > 0)
0921     emin = eminHitD;
0922 #ifdef EDM_ML_DEBUG
0923   edm::LogVerbatim("CaloSim") << "CaloSD::filterHit(..) Depth " << hit->getDepth() << " Emin = " << emin << " ("
0924                               << eminHit << ", " << eminHitD << ")";
0925 #endif
0926   return ((time <= tmaxHit) && (hit->getEnergyDeposit() > emin));
0927 }
0928 
0929 double CaloSD::getResponseWt(const G4Track* aTrack, int k) {
0930   double wt = 1.0;
0931   if (meanResponse[k].get()) {
0932     TrackInformation* trkInfo = cmsTrackInformation(aTrack);
0933     wt = meanResponse[k].get()->getWeight(trkInfo->genParticlePID(), trkInfo->genParticleP());
0934   }
0935   return wt;
0936 }
0937 
0938 void CaloSD::storeHit(CaloG4Hit* hit, int k) {
0939   if (hit == nullptr || previousID[k].trackID() < 0) {
0940     edm::LogWarning("CaloSim") << "CaloSD: hit to be stored is nullptr !!"
0941                                << " previousID.trackID()= " << previousID[k].trackID();
0942     return;
0943   }
0944 
0945   theHC[k]->insert(hit);
0946   if (useMap)
0947     hitMap[k].insert(std::pair<CaloHitID, CaloG4Hit*>(previousID[k], hit));
0948 }
0949 
0950 bool CaloSD::saveHit(CaloG4Hit* aHit, int k) {
0951   int tkID;
0952   bool ok = true;
0953 
0954   double time = aHit->getTimeSlice();
0955   if (corrTOFBeam)
0956     time += correctT;
0957 
0958   // More strict bookkeeping for finecalo
0959   if (doFineCalo_ && aHit->isFinecaloTrackID()) {
0960 #ifdef EDM_ML_DEBUG
0961     edm::LogVerbatim("DoFineCalo") << "Saving hit " << shortreprID(aHit);
0962 #endif
0963     if (!m_trackManager)
0964       throw cms::Exception("Unknown", "CaloSD") << "m_trackManager not set, needed for finecalo!";
0965     if (!m_trackManager->trackExists(aHit->getTrackID()))
0966       throw cms::Exception("Unknown", "CaloSD")
0967           << "Error on hit " << shortreprID(aHit) << ": Parent track not in track manager";
0968     slave[k].get()->processHits(aHit->getUnitID(),
0969                                 aHit->getEM() / CLHEP::GeV,
0970                                 aHit->getHadr() / CLHEP::GeV,
0971                                 time,
0972                                 aHit->getTrackID(),
0973                                 aHit->getDepth());
0974   }
0975   // Regular, not-fine way:
0976   else {
0977     if (m_trackManager) {
0978       tkID = m_trackManager->giveMotherNeeded(aHit->getTrackID());
0979       if (tkID == 0) {
0980         if (m_trackManager->trackExists(aHit->getTrackID()))
0981           tkID = (aHit->getTrackID());
0982         else {
0983           ok = false;
0984         }
0985       }
0986     } else {
0987       tkID = aHit->getTrackID();
0988       ok = false;
0989     }
0990 #ifdef EDM_ML_DEBUG
0991     edm::LogVerbatim("DoFineCalo") << "Saving hit " << shortreprID(aHit) << " with trackID=" << tkID;
0992 #endif
0993 
0994     slave[k].get()->processHits(
0995         aHit->getUnitID(), aHit->getEM() / CLHEP::GeV, aHit->getHadr() / CLHEP::GeV, time, tkID, aHit->getDepth());
0996   }
0997 
0998 #ifdef EDM_ML_DEBUG
0999   if (!ok)
1000     edm::LogWarning("CaloSim") << "CaloSD:Cannot find track ID for " << aHit->getTrackID();
1001   edm::LogVerbatim("CaloSim") << "CalosD: Track ID " << aHit->getTrackID() << " changed to " << tkID
1002                               << " by SimTrackManager Status " << ok;
1003 #endif
1004 
1005 #ifdef EDM_ML_DEBUG
1006   edm::LogVerbatim("CaloSim") << "CaloSD: Store Hit at " << std::hex << aHit->getUnitID() << std::dec << " "
1007                               << aHit->getDepth() << " due to " << tkID << " in time " << time << " of energy "
1008                               << aHit->getEM() / CLHEP::GeV << " GeV (EM) and " << aHit->getHadr() / CLHEP::GeV
1009                               << " GeV (Hadr)";
1010 #endif
1011   return ok;
1012 }
1013 
1014 void CaloSD::update(const BeginOfTrack* trk) {
1015   int primary = -1;
1016   TrackInformation* trkInfo = cmsTrackInformation((*trk)());
1017   if (trkInfo->isPrimary())
1018     primary = (*trk)()->GetTrackID();
1019 
1020 #ifdef EDM_ML_DEBUG
1021   edm::LogVerbatim("CaloSim") << "New track: isPrimary " << trkInfo->isPrimary() << " primary ID = " << primary
1022                               << " primary ancestor ID " << primAncestor;
1023 #endif
1024 
1025   // update the information if a different primary track ID
1026 
1027   if (primary > 0 && primary != primAncestor) {
1028     primAncestor = primary;
1029 
1030     // clean the hits information
1031 
1032     bool clean(false);
1033     for (int k = 0; k < nHC_; ++k)
1034       if (theHC[k]->entries() > 0)
1035         clean = true;
1036     if (clean)
1037       cleanHitCollection();
1038   }
1039 }
1040 
1041 void CaloSD::cleanHitCollection() {
1042   for (int k = 0; k < nHC_; ++k) {
1043     std::vector<CaloG4Hit*>* theCollection = theHC[k]->GetVector();
1044 
1045 #ifdef EDM_ML_DEBUG
1046     edm::LogVerbatim("CaloSim") << "CaloSD: collection before merging, size = " << theHC[k]->entries();
1047 #endif
1048     if (reusehit[k].empty())
1049       reusehit[k].reserve(theHC[k]->entries() - cleanIndex[k]);
1050 
1051     // if no map used, merge before hits to have the save situation as a map
1052     if (!useMap) {
1053       std::vector<CaloG4Hit*> hitvec;
1054 
1055       hitvec.swap(*theCollection);
1056       sort((hitvec.begin() + cleanIndex[k]), hitvec.end(), CaloG4HitLess());
1057 #ifdef EDM_ML_DEBUG
1058       edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: sort hits in buffer starting from "
1059                                   << "element = " << cleanIndex[k];
1060       for (unsigned int i = 0; i < hitvec.size(); ++i) {
1061         if (hitvec[i] == nullptr)
1062           edm::LogVerbatim("CaloSim") << i << " has a null pointer";
1063         else
1064           edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i];
1065       }
1066 #endif
1067       CaloG4HitEqual equal;
1068       for (unsigned int i = cleanIndex[k]; i < hitvec.size(); ++i) {
1069         int jump = 0;
1070         for (unsigned int j = i + 1; j < hitvec.size() && equal(hitvec[i], hitvec[j]); ++j) {
1071           ++jump;
1072           // merge j to i
1073           (*hitvec[i]).addEnergyDeposit(*hitvec[j]);
1074           (*hitvec[j]).setEM(0.);
1075           (*hitvec[j]).setHadr(0.);
1076           reusehit[k].emplace_back(hitvec[j]);
1077           hitvec[j] = nullptr;
1078         }
1079         i += jump;
1080       }
1081 #ifdef EDM_ML_DEBUG
1082       edm::LogVerbatim("CaloSim") << "CaloSD: cleanHitCollection merge the hits in buffer ";
1083       for (unsigned int i = 0; i < hitvec.size(); ++i) {
1084         if (hitvec[i] == nullptr)
1085           edm::LogVerbatim("CaloSim") << i << " has a null pointer";
1086         else
1087           edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i];
1088       }
1089 #endif
1090       //move all nullptr to end of list and then remove them
1091       hitvec.erase(std::stable_partition(
1092                        hitvec.begin() + cleanIndex[k], hitvec.end(), [](CaloG4Hit* p) { return p != nullptr; }),
1093                    hitvec.end());
1094 #ifdef EDM_ML_DEBUG
1095       edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: remove the merged hits in buffer, new size = "
1096                                   << hitvec.size();
1097 #endif
1098       hitvec.swap(*theCollection);
1099       totalHits[k] = theHC[k]->entries();
1100     }
1101 
1102 #ifdef EDM_ML_DEBUG
1103     edm::LogVerbatim("CaloSim") << "CaloSD: collection after merging, size= " << theHC[k]->entries()
1104                                 << " Size of reusehit= " << reusehit[k].size()
1105                                 << "\n      starting hit selection from index = " << cleanIndex[k];
1106 #endif
1107 
1108     int addhit = 0;
1109     for (unsigned int i = cleanIndex[k]; i < theCollection->size(); ++i) {
1110       CaloG4Hit* aHit((*theCollection)[i]);
1111 
1112       // selection
1113 
1114       double time = aHit->getTimeSlice();
1115       if (corrTOFBeam)
1116         time += correctT;
1117       if (!filterHit(aHit, time)) {
1118 #ifdef EDM_ML_DEBUG
1119         edm::LogVerbatim("CaloSim") << "CaloSD: dropped CaloG4Hit "
1120                                     << " " << *aHit;
1121 #endif
1122 
1123         // create the list of hits to be reused
1124 
1125         reusehit[k].emplace_back((*theCollection)[i]);
1126         (*theCollection)[i] = nullptr;
1127         ++addhit;
1128       }
1129     }
1130 
1131 #ifdef EDM_ML_DEBUG
1132     edm::LogVerbatim("CaloSim") << "CaloSD: Size of reusehit after selection = " << reusehit[k].size()
1133                                 << " Number of added hit = " << addhit;
1134 #endif
1135     if (useMap) {
1136       if (addhit > 0) {
1137         int offset = reusehit[k].size() - addhit;
1138         for (int ii = addhit - 1; ii >= 0; --ii) {
1139           CaloHitID theID = reusehit[k][offset + ii]->getID();
1140           hitMap[k].erase(theID);
1141         }
1142       }
1143     }
1144 
1145     //move all nullptr to end of list and then remove them
1146     theCollection->erase(
1147         std::stable_partition(
1148             theCollection->begin() + cleanIndex[k], theCollection->end(), [](CaloG4Hit* p) { return p != nullptr; }),
1149         theCollection->end());
1150 #ifdef EDM_ML_DEBUG
1151     edm::LogVerbatim("CaloSim") << "CaloSD: hit collection after selection, size = " << theHC[k]->entries();
1152     theHC[k]->PrintAllHits();
1153 #endif
1154 
1155     cleanIndex[k] = theHC[k]->entries();
1156   }
1157 }
1158 
1159 void CaloSD::printDetectorLevels(const G4VTouchable* touch) const {
1160   //Print name and copy numbers
1161   int level = ((touch->GetHistoryDepth()) + 1);
1162   std::ostringstream st1;
1163   st1 << level << " Levels:";
1164   if (level > 0) {
1165     for (int ii = 0; ii < level; ii++) {
1166       int i = level - ii - 1;
1167       G4VPhysicalVolume* pv = touch->GetVolume(i);
1168       std::string name = (pv != nullptr) ? pv->GetName() : "Unknown";
1169       st1 << " " << name << ":" << touch->GetReplicaNumber(i);
1170     }
1171   }
1172   edm::LogVerbatim("CaloSim") << st1.str();
1173 }
[ 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!