Tracker/src/TkAccumulatingSensitiveDetector.cc

0001 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0002
0003 #include "SimG4Core/SensitiveDetector/interface/FrameRotation.h"
0004 #include "SimG4Core/Notification/interface/TrackInformation.h"
0005
0006 #include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h"
0007 #include "SimDataFormats/SimHitMaker/interface/TrackingSlaveSD.h"
0008
0009 #include "SimG4CMS/Tracker/interface/TkAccumulatingSensitiveDetector.h"
0010 #include "SimG4CMS/Tracker/interface/FakeFrameRotation.h"
0011 #include "SimG4CMS/Tracker/interface/TrackerFrameRotation.h"
0012 #include "SimG4CMS/Tracker/interface/TkSimHitPrinter.h"
0013 #include "SimG4CMS/Tracker/interface/TrackerG4SimHitNumberingScheme.h"
0014
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016
0017 #include "SimG4Core/Notification/interface/TrackInformation.h"
0018 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0019 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0020
0021 #include "G4Track.hh"
0022 #include "G4StepPoint.hh"
0023 #include "G4VProcess.hh"
0024
0025 #include <CLHEP/Units/SystemOfUnits.h>
0026
0027 #include <memory>
0028
0029 #include <iostream>
0030 #include <vector>
0031
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033
0034 //#define FAKEFRAMEROTATION
0035
0036 static TrackerG4SimHitNumberingScheme& numberingScheme(const GeometricDet& det) {
0037   static thread_local TrackerG4SimHitNumberingScheme s_scheme(det);
0038   return s_scheme;
0039 }
0040
0041 TkAccumulatingSensitiveDetector::TkAccumulatingSensitiveDetector(const std::string& name,
0042                                                                  const GeometricDet* pDD,
0043                                                                  const SensitiveDetectorCatalog& clg,
0044                                                                  edm::ParameterSet const& p,
0045                                                                  const SimTrackManager* manager)
0046     : SensitiveTkDetector(name, clg),
0047       pDD_(pDD),
0048       theManager(manager),
0049       rTracker(1200. * CLHEP::mm),
0050       zTracker(3000. * CLHEP::mm),
0051       mySimHit(nullptr),
0052       lastId(0),
0053       lastTrack(0),
0054       oldVolume(nullptr),
0055       px(0.0f),
0056       py(0.0f),
0057       pz(0.0f),
0058       eventno(0),
0059       pname("") {
0060   edm::ParameterSet m_TrackerSD = p.getParameter<edm::ParameterSet>("TrackerSD");
0061   allowZeroEnergyLoss = m_TrackerSD.getParameter<bool>("ZeroEnergyLoss");
0062   neverAccumulate = m_TrackerSD.getParameter<bool>("NeverAccumulate");
0063   printHits = m_TrackerSD.getParameter<bool>("PrintHits");
0064   theTofLimit = m_TrackerSD.getParameter<double>("ElectronicSigmaInNanoSeconds") * 3 * CLHEP::ns;  // 3 sigma
0065   energyCut =
0066       m_TrackerSD.getParameter<double>("EnergyThresholdForPersistencyInGeV") * CLHEP::GeV;  //default must be 0.5
0067   energyHistoryCut =
0068       m_TrackerSD.getParameter<double>("EnergyThresholdForHistoryInGeV") * CLHEP::GeV;  //default must be 0.05
0069   rTracker2 = rTracker * rTracker;
0070
0071   // No Rotation given in input, automagically choose one based upon the name
0072   std::string rotType;
0073   theRotation = std::make_unique<TrackerFrameRotation>();
0074   rotType = "TrackerFrameRotation";
0075
0076 #ifdef FAKEFRAMEROTATION
0077   theRotation.reset(new FakeFrameRotation());
0078   rotType = "FakeFrameRotation";
0079 #endif
0080
0081   edm::LogVerbatim("TrackerSim") << " TkAccumulatingSensitiveDetector: "
0082                                  << " Criteria for Saving Tracker SimTracks: \n"
0083                                  << " History: " << energyHistoryCut << " MeV; Persistency: " << energyCut
0084                                  << " MeV;  TofLimit: " << theTofLimit << " ns"
0085                                  << "\n FrameRotation type " << rotType << " rTracker(cm)= " << rTracker / CLHEP::cm
0086                                  << " zTracker(cm)= " << zTracker / CLHEP::cm
0087                                  << " allowZeroEnergyLoss: " << allowZeroEnergyLoss
0088                                  << " neverAccumulate: " << neverAccumulate << " printHits: " << printHits;
0089
0090   slaveLowTof = std::make_unique<TrackingSlaveSD>(name + "LowTof");
0091   slaveHighTof = std::make_unique<TrackingSlaveSD>(name + "HighTof");
0092
0093   std::vector<std::string> temp;
0094   temp.push_back(slaveLowTof.get()->name());
0095   temp.push_back(slaveHighTof.get()->name());
0096   setNames(temp);
0097
0098   theG4ProcTypeEnumerator = std::make_unique<G4ProcessTypeEnumerator>();
0099   theNumberingScheme = nullptr;
0100 }
0101
0102 TkAccumulatingSensitiveDetector::~TkAccumulatingSensitiveDetector() {}
0103
0104 bool TkAccumulatingSensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
0105   LogDebug("TrackerSimDebug") << " Entering a new Step " << aStep->GetTotalEnergyDeposit() << " "
0106                               << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
0107
0108   if (aStep->GetTotalEnergyDeposit() > 0. || allowZeroEnergyLoss) {
0109     if (!mySimHit) {
0110       createHit(aStep);
0111     } else if (neverAccumulate || newHit(aStep)) {
0112       sendHit();
0113       createHit(aStep);
0114     } else {
0115       updateHit(aStep);
0116     }
0117     return true;
0118   }
0119   return false;
0120 }
0121
0122 uint32_t TkAccumulatingSensitiveDetector::setDetUnitId(const G4Step* step) {
0123   return theNumberingScheme->g4ToNumberingScheme(step->GetPreStepPoint()->GetTouchable());
0124 }
0125
0126 void TkAccumulatingSensitiveDetector::update(const BeginOfTrack* bot) {
0127   const G4Track* gTrack = (*bot)();
0128
0129 #ifdef DUMPPROCESSES
0130   if (gTrack->GetCreatorProcess()) {
0131     edm::LogVerbatim("TrackerSim") << " -> PROCESS CREATOR : " << gTrack->GetCreatorProcess()->GetProcessName();
0132   } else {
0133     edm::LogVerbatim("TrackerSim") << " -> No Creator process";
0134   }
0135 #endif
0136
0137   //
0138   //Position
0139   //
0140   const G4ThreeVector& pos = gTrack->GetPosition();
0141   LogDebug("TrackerSimDebug") << " update(..) of " << gTrack->GetDefinition()->GetParticleName()
0142                               << " trackID= " << gTrack->GetTrackID() << " E(MeV)= " << gTrack->GetKineticEnergy()
0143                               << " Ecut= " << energyCut << " R(mm)= " << pos.perp() << " Z(mm)= " << pos.z();
0144
0145   //
0146   // Check if in Tracker Volume
0147   //
0148   if (pos.x() * pos.x() + pos.y() * pos.y() < rTracker2 && std::abs(pos.z()) < zTracker) {
0149     //
0150     // inside the Tracker
0151     //
0152     TrackInformation* info = nullptr;
0153     if (gTrack->GetKineticEnergy() > energyCut) {
0154       info = cmsTrackInformation(gTrack);
0155       info->setStoreTrack();
0156     }
0157     //
0158     // Save History?
0159     //
0160     if (gTrack->GetKineticEnergy() > energyHistoryCut) {
0161       if (nullptr == info) {
0162         info = cmsTrackInformation(gTrack);
0163       }
0164       info->putInHistory();
0165       LogDebug("TrackerSimDebug") << " Track inside the tracker selected for HISTORY"
0166                                   << " Track ID= " << gTrack->GetTrackID();
0167     }
0168   }
0169 }
0170
0171 void TkAccumulatingSensitiveDetector::sendHit() {
0172   if (mySimHit == nullptr)
0173     return;
0174   if (printHits) {
0175     TkSimHitPrinter thePrinter("TkHitPositionOSCAR.dat");
0176     thePrinter.startNewSimHit(GetName(),
0177                               oldVolume->GetLogicalVolume()->GetName(),
0178                               mySimHit->detUnitId(),
0179                               mySimHit->trackId(),
0180                               lastTrack,
0181                               eventno);
0182     thePrinter.printLocal(mySimHit->entryPoint(), mySimHit->exitPoint());
0183     thePrinter.printGlobal(globalEntryPoint, globalExitPoint);
0184     thePrinter.printHitData(pname, mySimHit->pabs(), mySimHit->energyLoss(), mySimHit->timeOfFlight());
0185     thePrinter.printGlobalMomentum(px, py, pz);
0186     LogDebug("TrackerSimDebug") << " Storing PSimHit: " << mySimHit->detUnitId() << " " << mySimHit->trackId() << " "
0187                                 << mySimHit->energyLoss() << " " << mySimHit->entryPoint() << " "
0188                                 << mySimHit->exitPoint();
0189   }
0190
0191   if (mySimHit->timeOfFlight() < theTofLimit) {
0192     slaveLowTof.get()->processHits(*mySimHit);  // implicit conversion (slicing) to PSimHit!!!
0193   } else {
0194     slaveHighTof.get()->processHits(*mySimHit);  // implicit conversion (slicing) to PSimHit!!!
0195   }
0196   //
0197   // clean up
0198   delete mySimHit;
0199   mySimHit = nullptr;
0200   lastTrack = 0;
0201   lastId = 0;
0202 }
0203
0204 void TkAccumulatingSensitiveDetector::createHit(const G4Step* aStep) {
0205   // VI: previous hit should be already deleted
0206   //     in past here was a check if a hit is inside a sensitive detector,
0207   //     this is not needed, because call to senstive detector happens
0208   //     only inside the volume
0209   const G4Track* theTrack = aStep->GetTrack();
0210   Local3DPoint theExitPoint = theRotation.get()->transformPoint(LocalPostStepPosition(aStep));
0211   Local3DPoint theEntryPoint;
0212   //
0213   //  Check particle type - for gamma and neutral hadrons energy deposition
0214   //  should be local (VI)
0215   //
0216   if (0.0 == theTrack->GetDefinition()->GetPDGCharge()) {
0217     theEntryPoint = theExitPoint;
0218   } else {
0219     theEntryPoint = theRotation.get()->transformPoint(LocalPreStepPosition(aStep));
0220   }
0221
0222   //
0223   //    This allows to send he skipEvent if it is outside!
0224   //
0225   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0226   float thePabs = preStepPoint->GetMomentum().mag() / CLHEP::GeV;
0227   float theTof = preStepPoint->GetGlobalTime() / CLHEP::nanosecond;
0228   float theEnergyLoss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0229   int theParticleType = G4TrackToParticleID::particleID(theTrack);
0230   uint32_t theDetUnitId = setDetUnitId(aStep);
0231   int theTrackID = theTrack->GetTrackID();
0232   if (theDetUnitId == 0) {
0233     edm::LogWarning("TkAccumulatingSensitiveDetector::createHit") << " theDetUnitId is not valid for " << GetName();
0234     throw cms::Exception("TkAccumulatingSensitiveDetector::createHit")
0235         << "cannot get theDetUnitId for G4Track " << theTrackID;
0236   }
0237
0238   // To whom assign the Hit?
0239   // First iteration: if the track is to be stored, use the current number;
0240   // otherwise, get to the mother
0241   unsigned int theTrackIDInsideTheSimHit = theTrackID;
0242
0243   const TrackInformation* temp = cmsTrackInformation(theTrack);
0244   if (!temp->storeTrack()) {
0245     // Go to the mother!
0246     theTrackIDInsideTheSimHit = theTrack->GetParentID();
0247     LogDebug("TrackerSimDebug") << " TkAccumulatingSensitiveDetector::createHit(): setting the TrackID from "
0248                                 << theTrackIDInsideTheSimHit << " to the mother one " << theTrackIDInsideTheSimHit
0249                                 << " " << theEnergyLoss;
0250   } else {
0251     LogDebug("TrackerSimDebug") << " TkAccumulatingSensitiveDetector:createHit(): leaving the current TrackID "
0252                                 << theTrackIDInsideTheSimHit;
0253   }
0254
0255   const G4ThreeVector& gmd = preStepPoint->GetMomentumDirection();
0256   // convert it to local frame
0257   G4ThreeVector lmd =
0258       ((G4TouchableHistory*)(preStepPoint->GetTouchable()))->GetHistory()->GetTopTransform().TransformAxis(gmd);
0259   Local3DPoint lnmd = theRotation.get()->transformPoint(ConvertToLocal3DPoint(lmd));
0260   float theThetaAtEntry = lnmd.theta();
0261   float thePhiAtEntry = lnmd.phi();
0262
0263   mySimHit = new UpdatablePSimHit(theEntryPoint,
0264                                   theExitPoint,
0265                                   thePabs,
0266                                   theTof,
0267                                   theEnergyLoss,
0268                                   theParticleType,
0269                                   theDetUnitId,
0270                                   theTrackIDInsideTheSimHit,
0271                                   theThetaAtEntry,
0272                                   thePhiAtEntry,
0273                                   theG4ProcTypeEnumerator.get()->processId(theTrack->GetCreatorProcess()));
0274   lastId = theDetUnitId;
0275   lastTrack = theTrackID;
0276
0277   // only for debugging
0278   if (printHits) {
0279     // point on Geant4 unit (mm)
0280     globalEntryPoint = ConvertToLocal3DPoint(preStepPoint->GetPosition());
0281     globalExitPoint = ConvertToLocal3DPoint(aStep->GetPostStepPoint()->GetPosition());
0282     // in CMS unit (GeV)
0283     px = preStepPoint->GetMomentum().x() / CLHEP::GeV;
0284     py = preStepPoint->GetMomentum().y() / CLHEP::GeV;
0285     pz = preStepPoint->GetMomentum().z() / CLHEP::GeV;
0286     oldVolume = preStepPoint->GetPhysicalVolume();
0287     pname = theTrack->GetDefinition()->GetParticleName();
0288     LogDebug("TrackerSimDebug") << " Created PSimHit: " << pname << " " << mySimHit->detUnitId() << " "
0289                                 << mySimHit->trackId() << " " << theTrackID
0290                                 << " p= " << aStep->GetPreStepPoint()->GetMomentum().mag() << " "
0291                                 << mySimHit->energyLoss() << " " << mySimHit->entryPoint() << " "
0292                                 << mySimHit->exitPoint();
0293   }
0294 }
0295
0296 void TkAccumulatingSensitiveDetector::updateHit(const G4Step* aStep) {
0297   // VI: in past here was a check if a hit is inside a sensitive detector,
0298   //     this is not needed, because call to senstive detector happens
0299   //     only inside the volume
0300   Local3DPoint theExitPoint = theRotation.get()->transformPoint(LocalPostStepPosition(aStep));
0301   float theEnergyLoss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0302   mySimHit->setExitPoint(theExitPoint);
0303   mySimHit->addEnergyLoss(theEnergyLoss);
0304   if (printHits) {
0305     globalExitPoint = ConvertToLocal3DPoint(aStep->GetPostStepPoint()->GetPosition());
0306     LogDebug("TrackerSimDebug") << " updateHit: for " << aStep->GetTrack()->GetDefinition()->GetParticleName()
0307                                 << " trackID= " << aStep->GetTrack()->GetTrackID() << " deltaEloss= " << theEnergyLoss
0308                                 << "\n Updated PSimHit: " << mySimHit->detUnitId() << " " << mySimHit->trackId() << " "
0309                                 << mySimHit->energyLoss() << " " << mySimHit->entryPoint() << " "
0310                                 << mySimHit->exitPoint();
0311   }
0312 }
0313
0314 bool TkAccumulatingSensitiveDetector::newHit(const G4Step* aStep) {
0315   const G4Track* theTrack = aStep->GetTrack();
0316
0317   // for neutral particles do not merge hits (V.I.)
0318   if (0.0 == theTrack->GetDefinition()->GetPDGCharge())
0319     return true;
0320
0321   uint32_t theDetUnitId = setDetUnitId(aStep);
0322   int theTrackID = theTrack->GetTrackID();
0323
0324   LogDebug("TrackerSimDebug") << "newHit: OLD(detID,trID) = (" << lastId << "," << lastTrack << "), NEW = ("
0325                               << theDetUnitId << "," << theTrackID << ") Step length(mm)= " << aStep->GetStepLength()
0326                               << " Edep= " << aStep->GetTotalEnergyDeposit()
0327                               << " p= " << aStep->GetPreStepPoint()->GetMomentum().mag();
0328   return ((theTrackID == lastTrack) && (lastId == theDetUnitId) && closeHit(aStep)) ? false : true;
0329 }
0330
0331 bool TkAccumulatingSensitiveDetector::closeHit(const G4Step* aStep) {
0332   const float tolerance2 = 0.0025f;  // (0.5 mm)^2 are allowed between entry and exit
0333   Local3DPoint theEntryPoint = theRotation.get()->transformPoint(LocalPreStepPosition(aStep));
0334   LogDebug("TrackerSimDebug") << " closeHit: distance = " << (mySimHit->exitPoint() - theEntryPoint).mag();
0335
0336   return ((mySimHit->exitPoint() - theEntryPoint).mag2() < tolerance2) ? true : false;
0337 }
0338
0339 void TkAccumulatingSensitiveDetector::EndOfEvent(G4HCofThisEvent*) {
0340   LogDebug("TrackerSimDebug") << " Saving the last hit in a ROU " << GetName();
0341   if (mySimHit != nullptr)
0342     sendHit();
0343 }
0344
0345 void TkAccumulatingSensitiveDetector::update(const BeginOfEvent* i) {
0346   clearHits();
0347   eventno = (*i)()->GetEventID();
0348   delete mySimHit;
0349   mySimHit = nullptr;
0350 }
0351
0352 void TkAccumulatingSensitiveDetector::update(const BeginOfJob* i) { theNumberingScheme = &(numberingScheme(*pDD_)); }
0353
0354 void TkAccumulatingSensitiveDetector::clearHits() {
0355   slaveLowTof.get()->Initialize();
0356   slaveHighTof.get()->Initialize();
0357 }
0358
0359 void TkAccumulatingSensitiveDetector::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0360   if (slaveLowTof.get()->name() == hname) {
0361     cc = slaveLowTof.get()->hits();
0362   } else if (slaveHighTof.get()->name() == hname) {
0363     cc = slaveHighTof.get()->hits();
0364   }
0365 }