Forward/src/TimingSD.cc

0001 //#define EDM_ML_DEBUG
0002
0003 ///////////////////////////////////////////////////////////////////////////////
0004 // File: TimingSD.cc
0005 // Date: 02.2006
0006 // Description: Sensitive Detector class for Timing
0007 // Modifications:
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   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->setStoreTrack();
0157     }
0158     if (incidentEnergy > energyHistoryCut) {
0159       if (nullptr == info) {
0160         info = cmsTrackInformation(theTrack);
0161       }
0162       info->putInHistory();
0163     }
0164 #ifdef EDM_ML_DEBUG
0165     if (info != nullptr) {
0166       LogDebug("TimingSim") << "TrackInformation for ID = " << theTrack->GetTrackID();
0167       info->Print();
0168     }
0169 #endif
0170   }
0171
0172   edeposit *= invgev;
0173   if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
0174     edepositEM = edeposit;
0175     edepositHAD = 0.f;
0176   } else {
0177     edepositEM = 0.f;
0178     edepositHAD = edeposit;
0179   }
0180   // time slice is defined for the entry point
0181   tSlice = timeFactor * preStepPoint->GetGlobalTime() * invns;
0182   tSliceID = (int)tSlice;
0183
0184   unitID = setDetUnitId(aStep);
0185   primaryID = getTrackID(theTrack);
0186 }
0187
0188 bool TimingSD::hitExists(const G4Step* aStep) {
0189   if (!currentHit) {
0190     return false;
0191   }
0192
0193   // Update if in the same detector and time-slice
0194   if (tSliceID == tsID && unitID == previousUnitID) {
0195     updateHit();
0196     return true;
0197   }
0198
0199   //look in the HitContainer whether a hit with the same primID, unitID,
0200   //tSliceID already exists:
0201
0202   bool found = false;
0203   int thehc_entries = theHC->entries();
0204   for (int j = 0; j < thehc_entries; ++j) {
0205     BscG4Hit* aHit = (*theHC)[j];
0206     if (aHit->getTimeSliceID() == tSliceID && aHit->getUnitID() == unitID) {
0207       if (checkHit(aStep, aHit)) {
0208         currentHit = aHit;
0209         found = true;
0210         break;
0211       }
0212     }
0213   }
0214   if (found) {
0215     updateHit();
0216   }
0217   return found;
0218 }
0219
0220 bool TimingSD::checkHit(const G4Step*, BscG4Hit* hit) {
0221   // change hit info to fastest primary particle
0222   if (tof < hit->getTof()) {
0223     hit->setTrackID(primaryID);
0224     hit->setIncidentEnergy((float)incidentEnergy);
0225     hit->setPabs(float(preStepPoint->GetMomentum().mag()) * invgev);
0226     hit->setTof(tof);
0227     hit->setParticleType(theTrack->GetDefinition()->GetPDGEncoding());
0228
0229     float ThetaAtEntry = (float)(hitPointLocal.theta() * invdeg);
0230     float PhiAtEntry = (float)(hitPointLocal.phi() * invdeg);
0231
0232     hit->setThetaAtEntry(ThetaAtEntry);
0233     hit->setPhiAtEntry(PhiAtEntry);
0234
0235     hit->setEntry(hitPoint);
0236     hit->setEntryLocalP(hitPointLocal);
0237     hit->setExitLocalP(hitPointLocalExit);
0238
0239     hit->setParentId(theTrack->GetParentID());
0240     hit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess()));
0241
0242     hit->setVertexPosition(theTrack->GetVertexPosition());
0243   }
0244   return true;
0245 }
0246
0247 void TimingSD::storeHit(BscG4Hit* hit) {
0248   if (primID < 0)
0249     return;
0250   if (hit == nullptr) {
0251     edm::LogWarning("BscSim") << "BscSD: hit to be stored is NULL !!";
0252     return;
0253   }
0254
0255   theHC->insert(hit);
0256 }
0257
0258 void TimingSD::createNewHit(const G4Step* aStep) {
0259 #ifdef EDM_ML_DEBUG
0260   const G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
0261   edm::LogVerbatim("TimingSim") << "TimingSD CreateNewHit for " << GetName() << " PV " << currentPV->GetName()
0262                                 << " PVid = " << currentPV->GetCopyNo() << " Unit " << unitID << "\n primary "
0263                                 << primaryID << " Tof(ns)= " << tof << " time slice " << tSliceID
0264                                 << " E(MeV)= " << incidentEnergy << " trackID " << theTrack->GetTrackID() << " "
0265                                 << theTrack->GetDefinition()->GetParticleName() << " parentID "
0266                                 << theTrack->GetParentID();
0267
0268   if (theTrack->GetCreatorProcess() != nullptr) {
0269     edm::LogVerbatim("TimingSim") << theTrack->GetCreatorProcess()->GetProcessName();
0270   } else {
0271     edm::LogVerbatim("TimingSim") << " is primary particle";
0272   }
0273 #endif
0274
0275   currentHit = new BscG4Hit;
0276   currentHit->setTrackID(primaryID);
0277   currentHit->setTimeSlice(tSlice);
0278   currentHit->setUnitID(unitID);
0279   currentHit->setIncidentEnergy((float)incidentEnergy);
0280
0281   currentHit->setPabs(float(preStepPoint->GetMomentum().mag() * invgev));
0282   currentHit->setTof(tof);
0283   currentHit->setParticleType(theTrack->GetDefinition()->GetPDGEncoding());
0284
0285   float ThetaAtEntry = hitPointLocal.theta() * invdeg;
0286   float PhiAtEntry = hitPointLocal.phi() * invdeg;
0287
0288   currentHit->setThetaAtEntry(ThetaAtEntry);
0289   currentHit->setPhiAtEntry(PhiAtEntry);
0290
0291   currentHit->setEntry(hitPoint);
0292   currentHit->setEntryLocalP(hitPointLocal);
0293   currentHit->setExitLocalP(hitPointLocalExit);
0294
0295   currentHit->setParentId(theTrack->GetParentID());
0296   currentHit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess()));
0297
0298   currentHit->setVertexPosition(theTrack->GetVertexPosition());
0299
0300   updateHit();
0301   storeHit(currentHit);
0302 }
0303
0304 void TimingSD::updateHit() {
0305   currentHit->addEnergyDeposit(edepositEM, edepositHAD);
0306
0307 #ifdef EDM_ML_DEBUG
0308   edm::LogVerbatim("TimingSim") << "updateHit: " << GetName() << " add eloss(GeV) " << edeposit
0309                                 << "CurrentHit=" << currentHit << ", PostStepPoint= " << postStepPoint->GetPosition();
0310 #endif
0311
0312   // buffer for next steps:
0313   tsID = tSliceID;
0314   primID = primaryID;
0315   previousUnitID = unitID;
0316 }
0317
0318 void TimingSD::setToLocal(const G4StepPoint* stepPoint, const G4ThreeVector& globalPoint, G4ThreeVector& localPoint) {
0319   const G4VTouchable* touch = stepPoint->GetTouchable();
0320   localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(globalPoint);
0321 }
0322
0323 void TimingSD::EndOfEvent(G4HCofThisEvent*) {
0324   int nhits = theHC->entries();
0325   if (0 == nhits) {
0326     return;
0327   }
0328   // here we loop over transient hits and make them persistent
0329   for (int j = 0; j < nhits; ++j) {
0330     BscG4Hit* aHit = (*theHC)[j];
0331     Local3DPoint locEntryPoint = ConvertToLocal3DPoint(aHit->getEntryLocalP());
0332     Local3DPoint locExitPoint = ConvertToLocal3DPoint(aHit->getExitLocalP());
0333
0334 #ifdef EDM_ML_DEBUG
0335     edm::LogVerbatim("TimingSim") << "TimingSD: Hit for storage \n"
0336                                   << *aHit << "\n Entry point: " << locEntryPoint << "\n Exit  point: " << locExitPoint;
0337 #endif
0338
0339     slave->processHits(PSimHit(locEntryPoint,
0340                                locExitPoint,
0341                                aHit->getPabs(),
0342                                aHit->getTof(),
0343                                aHit->getEnergyLoss(),
0344                                aHit->getParticleType(),
0345                                aHit->getUnitID(),
0346                                aHit->getTrackID(),
0347                                aHit->getThetaAtEntry(),
0348                                aHit->getPhiAtEntry(),
0349                                aHit->getProcessId()));
0350   }
0351 }
0352
0353 void TimingSD::PrintAll() {
0354 #ifdef EDM_ML_DEBUG
0355   edm::LogVerbatim("TimingSim") << "TimingSD: Collection " << theHC->GetName();
0356 #endif
0357   theHC->PrintAllHits();
0358 }
0359
0360 void TimingSD::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0361   if (slave->name() == hname) {
0362     cc = slave->hits();
0363   }
0364 }
0365
0366 void TimingSD::update(const BeginOfEvent* i) {
0367 #ifdef EDM_ML_DEBUG
0368   edm::LogVerbatim("TimingSim") << " Dispatched BeginOfEvent for " << GetName();
0369 #endif
0370   clearHits();
0371 }
0372
0373 void TimingSD::clearHits() { slave->Initialize(); }
0374
0375 int TimingSD::getTrackID(const G4Track* aTrack) {
0376   LogDebug("TimingSim") << "primary ID: " << aTrack->GetTrackID();
0377   return aTrack->GetTrackID();
0378 }