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       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     // (primaryID = -2) means ECAL GFlash spot inside MTD
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   // exclude ECAL gflash spots inside MTD detectors,
0118   // which not possible correctly handle
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   // neutral particles deliver energy post step
0130   // charged particle start deliver energy pre step
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   // should MC truth be saved
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   // time slice is defined for the entry point
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   // Update if in the same detector and time-slice
0204   if (tSliceID == tsID && unitID == previousUnitID) {
0205     updateHit();
0206     return true;
0207   }
0208
0209   //look in the HitContainer whether a hit with the same primID, unitID,
0210   //tSliceID already exists:
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   // change hit info to fastest primary particle
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   // buffer for next steps:
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   // here we loop over transient hits and make them persistent
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; }