Forward/src/TimingSD.cc

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