File indexing completed on 2025-04-08 03:34:56
0001
0002
0003
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
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
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
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
0207 doFineCaloThisStep_ = (doFineCalo_ && isItFineCalo(aStep->GetPreStepPoint()->GetTouchable()));
0208
0209
0210
0211 if (isParameterized) {
0212 if (getFromLibrary(aStep)) {
0213
0214
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
0230 edepositEM = edepositHAD = 0.f;
0231 if (aStep->GetTotalEnergyDeposit() <= 0.0) {
0232 return false;
0233 }
0234
0235
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
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
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
0339 if (currentID[0] == previousID[0]) {
0340 updateHit(currentHit[0], 0);
0341 } else {
0342 posGlobal = aSpot->GetEnergySpot()->GetPosition();
0343
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
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
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
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
0474 if (currentID[k] == previousID[k]) {
0475 updateHit(currentHit[k], k);
0476 return true;
0477 }
0478
0479
0480
0481
0482
0483
0484
0485
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
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
0532
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
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
0554 std::string CaloSD::shortreprID(const CaloG4Hit* hit) { return shortreprID(hit->getID()); }
0555
0556
0557
0558
0559 unsigned int CaloSD::findBoundaryCrossingParent(const G4Track* track, bool markAsSaveable) {
0560 TrackInformation* trkInfo = cmsTrackInformation(track);
0561 unsigned int id = track->GetTrackID();
0562
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
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
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
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
0598 for (auto ancestorID : decayChain)
0599 boundaryCrossingParentMap_[ancestorID] = it->second;
0600 #ifdef EDM_ML_DEBUG
0601
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
0612 TrackWithHistory* parentTrack = m_trackManager->getTrackByID(parentID, true);
0613 if (parentTrack->crossedBoundary()) {
0614 if (markAsSaveable)
0615 parentTrack->setToBeSaved();
0616 decayChain.push_back(parentID);
0617
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
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
0674 if (!doFineCaloThisStep_) {
0675 double etrack = 0;
0676 if (currentID[k].trackID() == primIDSaved[k]) {
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
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;
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
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
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
1038
1039 if (primary > 0 && primary != primAncestor) {
1040 primAncestor = primary;
1041
1042
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
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
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
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
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
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
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
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 }