Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-27 02:50:34

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; }