File indexing completed on 2025-01-27 02:50:34
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "SimG4CMS/Forward/interface/TimingSD.h"
0011
0012 #include "SimG4Core/Notification/interface/TrackInformation.h"
0013 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0014 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0015
0016 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "G4Step.hh"
0020 #include "G4StepPoint.hh"
0021 #include "G4Track.hh"
0022 #include "G4VPhysicalVolume.hh"
0023 #include "G4SDManager.hh"
0024 #include "G4VProcess.hh"
0025
0026 #include "G4PhysicalConstants.hh"
0027 #include <CLHEP/Units/SystemOfUnits.h>
0028
0029 #include <vector>
0030 #include <iostream>
0031
0032 static const float invgev = 1.0 / CLHEP::GeV;
0033 static const double invns = 1.0 / CLHEP::nanosecond;
0034 static const double invdeg = 1.0 / CLHEP::deg;
0035
0036
0037 TimingSD::TimingSD(const std::string& name, const SensitiveDetectorCatalog& clg, const SimTrackManager* manager)
0038 : SensitiveTkDetector(name, clg),
0039 theManager(manager),
0040 theHC(nullptr),
0041 currentHit(nullptr),
0042 theTrack(nullptr),
0043 preStepPoint(nullptr),
0044 postStepPoint(nullptr),
0045 unitID(0),
0046 previousUnitID(0),
0047 primID(-2),
0048 hcID(-1),
0049 tsID(-2),
0050 primaryID(0),
0051 tSliceID(-1),
0052 timeFactor(1.0),
0053 energyCut(1.e+9),
0054 energyHistoryCut(1.e+9),
0055 hitClassID(0) {
0056 slave = new TrackingSlaveSD(name);
0057 theEnumerator = new G4ProcessTypeEnumerator();
0058 }
0059
0060 TimingSD::~TimingSD() {
0061 delete slave;
0062 delete theEnumerator;
0063 }
0064
0065 void TimingSD::Initialize(G4HCofThisEvent* HCE) {
0066 edm::LogVerbatim("TimingSim") << "TimingSD : Initialize called for " << GetName() << " time slice factor "
0067 << timeFactor << "\n MC truth cuts in are " << energyCut / CLHEP::GeV << " GeV and "
0068 << energyHistoryCut / CLHEP::GeV << " GeV";
0069
0070 theHC = new BscG4HitCollection(GetName(), collectionName[0]);
0071 if (hcID < 0) {
0072 hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0073 }
0074 HCE->AddHitsCollection(hcID, theHC);
0075
0076 tsID = -2;
0077 primID = -2;
0078 }
0079
0080 void TimingSD::setTimeFactor(double val) {
0081 if (val <= 0.0) {
0082 return;
0083 }
0084 timeFactor = val;
0085 #ifdef EDM_ML_DEBUG
0086 edm::LogVerbatim("TimingSim") << "TimingSD : for " << GetName() << " time slice factor is set to " << timeFactor;
0087 #endif
0088 }
0089
0090 void TimingSD::setCuts(double eCut, double historyCut) {
0091 if (eCut > 0.) {
0092 energyCut = eCut;
0093 }
0094 if (historyCut > 0.) {
0095 energyHistoryCut = historyCut;
0096 }
0097 #ifdef EDM_ML_DEBUG
0098 edm::LogVerbatim("TimingSim") << "TimingSD : for " << GetName() << " MC truth cuts in are " << energyCut / CLHEP::GeV
0099 << " GeV and " << energyHistoryCut / CLHEP::GeV << " GeV";
0100 #endif
0101 }
0102
0103 bool TimingSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0104 edeposit = aStep->GetTotalEnergyDeposit();
0105 if (edeposit > 0.f) {
0106 getStepInfo(aStep);
0107
0108 if (-1 <= primaryID && !hitExists(aStep)) {
0109 createNewHit(aStep);
0110 }
0111 }
0112 return true;
0113 }
0114
0115 void TimingSD::getStepInfo(const G4Step* aStep) {
0116 const G4Track* newTrack = aStep->GetTrack();
0117
0118
0119 primaryID = getTrackID(newTrack);
0120 if (primaryID < -1) {
0121 return;
0122 }
0123
0124 preStepPoint = aStep->GetPreStepPoint();
0125 postStepPoint = aStep->GetPostStepPoint();
0126 hitPointExit = postStepPoint->GetPosition();
0127 setToLocal(preStepPoint, hitPointExit, hitPointLocalExit);
0128
0129
0130
0131 if (newTrack->GetDefinition()->GetPDGCharge() == 0.0) {
0132 hitPoint = hitPointExit;
0133 hitPointLocal = hitPointLocalExit;
0134 tof = (float)(postStepPoint->GetGlobalTime() * invns);
0135 } else {
0136 hitPoint = preStepPoint->GetPosition();
0137 setToLocal(preStepPoint, hitPoint, hitPointLocal);
0138 tof = (float)(preStepPoint->GetGlobalTime() * invns);
0139 }
0140
0141 #ifdef EDM_ML_DEBUG
0142 double distGlobal =
0143 std::sqrt(std::pow(hitPoint.x() - hitPointExit.x(), 2) + std::pow(hitPoint.y() - hitPointExit.y(), 2) +
0144 std::pow(hitPoint.z() - hitPointExit.z(), 2));
0145 double distLocal = std::sqrt(std::pow(hitPointLocal.x() - hitPointLocalExit.x(), 2) +
0146 std::pow(hitPointLocal.y() - hitPointLocalExit.y(), 2) +
0147 std::pow(hitPointLocal.z() - hitPointLocalExit.z(), 2));
0148 LogDebug("TimingSim") << "TimingSD:"
0149 << "\n Global entry point: " << hitPoint << "\n Global exit point: " << hitPointExit
0150 << "\n Global step length: " << distGlobal << "\n Local entry point: " << hitPointLocal
0151 << "\n Local exit point: " << hitPointLocalExit << "\n Local step length: " << distLocal;
0152 if (std::fabs(distGlobal - distLocal) > 1.e-6) {
0153 LogDebug("TimingSim") << "DIFFERENCE IN DISTANCE \n";
0154 }
0155 #endif
0156
0157 incidentEnergy = preStepPoint->GetKineticEnergy();
0158
0159
0160 if (newTrack != theTrack) {
0161 theTrack = newTrack;
0162 TrackInformation* info = nullptr;
0163 if (incidentEnergy > energyCut) {
0164 info = cmsTrackInformation(theTrack);
0165 info->setStoreTrack();
0166 info->setIdLastStoredAncestor(theTrack->GetTrackID());
0167 }
0168 if (incidentEnergy > energyHistoryCut) {
0169 if (nullptr == info) {
0170 info = cmsTrackInformation(theTrack);
0171 }
0172 info->putInHistory();
0173 }
0174 #ifdef EDM_ML_DEBUG
0175 if (info != nullptr) {
0176 LogDebug("TimingSim") << "TrackInformation for ID = " << theTrack->GetTrackID();
0177 info->Print();
0178 }
0179 #endif
0180 }
0181
0182 edeposit *= invgev;
0183 if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0184 edepositEM = edeposit;
0185 edepositHAD = 0.f;
0186 } else {
0187 edepositEM = 0.f;
0188 edepositHAD = edeposit;
0189 }
0190
0191 tSlice = timeFactor * preStepPoint->GetGlobalTime() * invns;
0192 tSliceID = (int)tSlice;
0193
0194 setHitClassID(aStep);
0195 unitID = setDetUnitId(aStep);
0196 }
0197
0198 bool TimingSD::hitExists(const G4Step* aStep) {
0199 if (!currentHit) {
0200 return false;
0201 }
0202
0203
0204 if (tSliceID == tsID && unitID == previousUnitID) {
0205 updateHit();
0206 return true;
0207 }
0208
0209
0210
0211
0212 bool found = false;
0213 int thehc_entries = theHC->entries();
0214 for (int j = 0; j < thehc_entries; ++j) {
0215 BscG4Hit* aHit = (*theHC)[j];
0216 if (aHit->getTimeSliceID() == tSliceID && aHit->getUnitID() == unitID) {
0217 if (checkHit(aStep, aHit)) {
0218 currentHit = aHit;
0219 found = true;
0220 break;
0221 }
0222 }
0223 }
0224 if (found) {
0225 updateHit();
0226 }
0227 return found;
0228 }
0229
0230 bool TimingSD::checkHit(const G4Step*, BscG4Hit* hit) {
0231
0232 if (tof < hit->getTof()) {
0233 hit->setTrackID(primaryID);
0234 hit->setIncidentEnergy((float)incidentEnergy);
0235 hit->setPabs(float(preStepPoint->GetMomentum().mag()) * invgev);
0236 hit->setTof(tof);
0237 hit->setParticleType(theTrack->GetDefinition()->GetPDGEncoding());
0238
0239 float ThetaAtEntry = (float)(hitPointLocal.theta() * invdeg);
0240 float PhiAtEntry = (float)(hitPointLocal.phi() * invdeg);
0241
0242 hit->setThetaAtEntry(ThetaAtEntry);
0243 hit->setPhiAtEntry(PhiAtEntry);
0244
0245 hit->setEntry(hitPoint);
0246 hit->setEntryLocalP(hitPointLocal);
0247 hit->setExitLocalP(hitPointLocalExit);
0248
0249 hit->setParentId(theTrack->GetParentID());
0250 hit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess()));
0251
0252 hit->setVertexPosition(theTrack->GetVertexPosition());
0253 }
0254 return true;
0255 }
0256
0257 void TimingSD::storeHit(BscG4Hit* hit) {
0258 if (primID < 0)
0259 return;
0260 if (hit == nullptr) {
0261 edm::LogWarning("BscSim") << "BscSD: hit to be stored is NULL !!";
0262 return;
0263 }
0264
0265 theHC->insert(hit);
0266 }
0267
0268 void TimingSD::createNewHit(const G4Step* aStep) {
0269 #ifdef EDM_ML_DEBUG
0270 const G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
0271 edm::LogVerbatim("TimingSim") << "TimingSD CreateNewHit for " << GetName() << " PV " << currentPV->GetName()
0272 << " PVid = " << currentPV->GetCopyNo() << " Unit " << unitID << "\n primary "
0273 << primaryID << " Tof(ns)= " << tof << " time slice " << tSliceID
0274 << " E(MeV)= " << incidentEnergy << " trackID " << theTrack->GetTrackID() << " "
0275 << theTrack->GetDefinition()->GetParticleName() << " parentID "
0276 << theTrack->GetParentID();
0277
0278 if (theTrack->GetCreatorProcess() != nullptr) {
0279 edm::LogVerbatim("TimingSim") << theTrack->GetCreatorProcess()->GetProcessName();
0280 } else {
0281 edm::LogVerbatim("TimingSim") << " is primary particle";
0282 }
0283 #endif
0284
0285 currentHit = new BscG4Hit;
0286 currentHit->setTrackID(primaryID);
0287 currentHit->setTimeSlice(tSlice);
0288 currentHit->setUnitID(unitID);
0289 currentHit->setIncidentEnergy((float)incidentEnergy);
0290
0291 currentHit->setPabs(float(preStepPoint->GetMomentum().mag() * invgev));
0292 currentHit->setTof(tof);
0293 currentHit->setParticleType(theTrack->GetDefinition()->GetPDGEncoding());
0294
0295 float ThetaAtEntry = hitPointLocal.theta() * invdeg;
0296 float PhiAtEntry = hitPointLocal.phi() * invdeg;
0297
0298 currentHit->setThetaAtEntry(ThetaAtEntry);
0299 currentHit->setPhiAtEntry(PhiAtEntry);
0300
0301 currentHit->setEntry(hitPoint);
0302 currentHit->setEntryLocalP(hitPointLocal);
0303 currentHit->setExitLocalP(hitPointLocalExit);
0304
0305 currentHit->setParentId(theTrack->GetParentID());
0306 currentHit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess()));
0307
0308 currentHit->setVertexPosition(theTrack->GetVertexPosition());
0309
0310 updateHit();
0311 storeHit(currentHit);
0312 }
0313
0314 void TimingSD::updateHit() {
0315 currentHit->addEnergyDeposit(edepositEM, edepositHAD);
0316
0317 #ifdef EDM_ML_DEBUG
0318 edm::LogVerbatim("TimingSim") << "updateHit: " << GetName() << " add eloss(GeV) " << edeposit
0319 << "CurrentHit=" << currentHit << ", PostStepPoint= " << postStepPoint->GetPosition();
0320 #endif
0321
0322
0323 tsID = tSliceID;
0324 primID = primaryID;
0325 previousUnitID = unitID;
0326 }
0327
0328 void TimingSD::setToLocal(const G4StepPoint* stepPoint, const G4ThreeVector& globalPoint, G4ThreeVector& localPoint) {
0329 const G4VTouchable* touch = stepPoint->GetTouchable();
0330 localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(globalPoint);
0331 }
0332
0333 void TimingSD::EndOfEvent(G4HCofThisEvent*) {
0334 int nhits = theHC->entries();
0335 if (0 == nhits) {
0336 return;
0337 }
0338
0339 for (int j = 0; j < nhits; ++j) {
0340 BscG4Hit* aHit = (*theHC)[j];
0341 Local3DPoint locEntryPoint = ConvertToLocal3DPoint(aHit->getEntryLocalP());
0342 Local3DPoint locExitPoint = ConvertToLocal3DPoint(aHit->getExitLocalP());
0343
0344 #ifdef EDM_ML_DEBUG
0345 edm::LogVerbatim("TimingSim") << "TimingSD: Hit for storage \n"
0346 << *aHit << "\n Entry point: " << locEntryPoint << "\n Exit point: " << locExitPoint;
0347 #endif
0348
0349 slave->processHits(PSimHit(locEntryPoint,
0350 locExitPoint,
0351 aHit->getPabs(),
0352 aHit->getTof(),
0353 aHit->getEnergyLoss(),
0354 aHit->getParticleType(),
0355 aHit->getUnitID(),
0356 aHit->getTrackID(),
0357 aHit->getThetaAtEntry(),
0358 aHit->getPhiAtEntry(),
0359 aHit->getProcessId()));
0360 }
0361 }
0362
0363 void TimingSD::PrintAll() {
0364 #ifdef EDM_ML_DEBUG
0365 edm::LogVerbatim("TimingSim") << "TimingSD: Collection " << theHC->GetName();
0366 #endif
0367 theHC->PrintAllHits();
0368 }
0369
0370 void TimingSD::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0371 if (slave->name() == hname) {
0372 cc = slave->hits();
0373 }
0374 }
0375
0376 void TimingSD::update(const BeginOfEvent* i) {
0377 #ifdef EDM_ML_DEBUG
0378 edm::LogVerbatim("TimingSim") << " Dispatched BeginOfEvent for " << GetName();
0379 #endif
0380 clearHits();
0381 }
0382
0383 void TimingSD::clearHits() { slave->Initialize(); }
0384
0385 int TimingSD::getTrackID(const G4Track* aTrack) {
0386 LogDebug("TimingSim") << "primary ID: " << aTrack->GetTrackID();
0387 return aTrack->GetTrackID();
0388 }
0389
0390 void TimingSD::setHitClassID(const G4Step* aStep) { hitClassID = 0; }