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       info->setIdLastStoredAncestor(gTrack->GetTrackID());
0157     }
0158     //
0159     // Save History?
0160     //
0161     if (gTrack->GetKineticEnergy() > energyHistoryCut) {
0162       if (nullptr == info) {
0163         info = cmsTrackInformation(gTrack);
0164       }
0165       info->putInHistory();
0166       LogDebug("TrackerSimDebug") << " Track inside the tracker selected for HISTORY"
0167                                   << " Track ID= " << gTrack->GetTrackID();
0168     }
0169   }
0170 }
0171
0172 void TkAccumulatingSensitiveDetector::sendHit() {
0173   if (mySimHit == nullptr)
0174     return;
0175   if (printHits) {
0176     TkSimHitPrinter thePrinter("TkHitPositionOSCAR.dat");
0177     thePrinter.startNewSimHit(GetName(),
0178                               oldVolume->GetLogicalVolume()->GetName(),
0179                               mySimHit->detUnitId(),
0180                               mySimHit->trackId(),
0181                               lastTrack,
0182                               eventno);
0183     thePrinter.printLocal(mySimHit->entryPoint(), mySimHit->exitPoint());
0184     thePrinter.printGlobal(globalEntryPoint, globalExitPoint);
0185     thePrinter.printHitData(pname, mySimHit->pabs(), mySimHit->energyLoss(), mySimHit->timeOfFlight());
0186     thePrinter.printGlobalMomentum(px, py, pz);
0187     LogDebug("TrackerSimDebug") << " Storing PSimHit: " << mySimHit->detUnitId() << " " << mySimHit->trackId() << " "
0188                                 << mySimHit->energyLoss() << " " << mySimHit->entryPoint() << " "
0189                                 << mySimHit->exitPoint();
0190   }
0191
0192   if (mySimHit->timeOfFlight() < theTofLimit) {
0193     slaveLowTof.get()->processHits(*mySimHit);  // implicit conversion (slicing) to PSimHit!!!
0194   } else {
0195     slaveHighTof.get()->processHits(*mySimHit);  // implicit conversion (slicing) to PSimHit!!!
0196   }
0197   //
0198   // clean up
0199   delete mySimHit;
0200   mySimHit = nullptr;
0201   lastTrack = 0;
0202   lastId = 0;
0203 }
0204
0205 void TkAccumulatingSensitiveDetector::createHit(const G4Step* aStep) {
0206   // VI: previous hit should be already deleted
0207   //     in past here was a check if a hit is inside a sensitive detector,
0208   //     this is not needed, because call to senstive detector happens
0209   //     only inside the volume
0210   const G4Track* theTrack = aStep->GetTrack();
0211   Local3DPoint theExitPoint = theRotation.get()->transformPoint(LocalPostStepPosition(aStep));
0212   Local3DPoint theEntryPoint;
0213   //
0214   //  Check particle type - for gamma and neutral hadrons energy deposition
0215   //  should be local (VI)
0216   //
0217   if (0.0 == theTrack->GetDefinition()->GetPDGCharge()) {
0218     theEntryPoint = theExitPoint;
0219   } else {
0220     theEntryPoint = theRotation.get()->transformPoint(LocalPreStepPosition(aStep));
0221   }
0222
0223   //
0224   //    This allows to send he skipEvent if it is outside!
0225   //
0226   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0227   float thePabs = preStepPoint->GetMomentum().mag() / CLHEP::GeV;
0228   float theTof = preStepPoint->GetGlobalTime() / CLHEP::nanosecond;
0229   float theEnergyLoss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0230   int theParticleType = G4TrackToParticleID::particleID(theTrack);
0231   uint32_t theDetUnitId = setDetUnitId(aStep);
0232   int theTrackID = theTrack->GetTrackID();
0233   if (theDetUnitId == 0) {
0234     edm::LogWarning("TkAccumulatingSensitiveDetector::createHit") << " theDetUnitId is not valid for " << GetName();
0235     throw cms::Exception("TkAccumulatingSensitiveDetector::createHit")
0236         << "cannot get theDetUnitId for G4Track " << theTrackID;
0237   }
0238
0239   // To whom assign the Hit?
0240   // First iteration: if the track is to be stored, use the current number;
0241   // otherwise, get to the mother
0242   unsigned int theTrackIDInsideTheSimHit = theTrackID;
0243
0244   const TrackInformation* temp = cmsTrackInformation(theTrack);
0245   if (!temp->storeTrack()) {
0246     // Go to the mother!
0247     theTrackIDInsideTheSimHit = theTrack->GetParentID();
0248     LogDebug("TrackerSimDebug") << " TkAccumulatingSensitiveDetector::createHit(): setting the TrackID from "
0249                                 << theTrackIDInsideTheSimHit << " to the mother one " << theTrackIDInsideTheSimHit
0250                                 << " " << theEnergyLoss;
0251   } else {
0252     LogDebug("TrackerSimDebug") << " TkAccumulatingSensitiveDetector:createHit(): leaving the current TrackID "
0253                                 << theTrackIDInsideTheSimHit;
0254   }
0255
0256   const G4ThreeVector& gmd = preStepPoint->GetMomentumDirection();
0257   // convert it to local frame
0258   G4ThreeVector lmd =
0259       ((G4TouchableHistory*)(preStepPoint->GetTouchable()))->GetHistory()->GetTopTransform().TransformAxis(gmd);
0260   Local3DPoint lnmd = theRotation.get()->transformPoint(ConvertToLocal3DPoint(lmd));
0261   float theThetaAtEntry = lnmd.theta();
0262   float thePhiAtEntry = lnmd.phi();
0263
0264   mySimHit = new UpdatablePSimHit(theEntryPoint,
0265                                   theExitPoint,
0266                                   thePabs,
0267                                   theTof,
0268                                   theEnergyLoss,
0269                                   theParticleType,
0270                                   theDetUnitId,
0271                                   theTrackIDInsideTheSimHit,
0272                                   theThetaAtEntry,
0273                                   thePhiAtEntry,
0274                                   theG4ProcTypeEnumerator.get()->processId(theTrack->GetCreatorProcess()));
0275   lastId = theDetUnitId;
0276   lastTrack = theTrackID;
0277
0278   // only for debugging
0279   if (printHits) {
0280     // point on Geant4 unit (mm)
0281     globalEntryPoint = ConvertToLocal3DPoint(preStepPoint->GetPosition());
0282     globalExitPoint = ConvertToLocal3DPoint(aStep->GetPostStepPoint()->GetPosition());
0283     // in CMS unit (GeV)
0284     px = preStepPoint->GetMomentum().x() / CLHEP::GeV;
0285     py = preStepPoint->GetMomentum().y() / CLHEP::GeV;
0286     pz = preStepPoint->GetMomentum().z() / CLHEP::GeV;
0287     oldVolume = preStepPoint->GetPhysicalVolume();
0288     pname = theTrack->GetDefinition()->GetParticleName();
0289     LogDebug("TrackerSimDebug") << " Created PSimHit: " << pname << " " << mySimHit->detUnitId() << " "
0290                                 << mySimHit->trackId() << " " << theTrackID
0291                                 << " p= " << aStep->GetPreStepPoint()->GetMomentum().mag() << " "
0292                                 << mySimHit->energyLoss() << " " << mySimHit->entryPoint() << " "
0293                                 << mySimHit->exitPoint();
0294   }
0295 }
0296
0297 void TkAccumulatingSensitiveDetector::updateHit(const G4Step* aStep) {
0298   // VI: in past here was a check if a hit is inside a sensitive detector,
0299   //     this is not needed, because call to senstive detector happens
0300   //     only inside the volume
0301   Local3DPoint theExitPoint = theRotation.get()->transformPoint(LocalPostStepPosition(aStep));
0302   float theEnergyLoss = aStep->GetTotalEnergyDeposit() / CLHEP::GeV;
0303   mySimHit->setExitPoint(theExitPoint);
0304   mySimHit->addEnergyLoss(theEnergyLoss);
0305   if (printHits) {
0306     globalExitPoint = ConvertToLocal3DPoint(aStep->GetPostStepPoint()->GetPosition());
0307     LogDebug("TrackerSimDebug") << " updateHit: for " << aStep->GetTrack()->GetDefinition()->GetParticleName()
0308                                 << " trackID= " << aStep->GetTrack()->GetTrackID() << " deltaEloss= " << theEnergyLoss
0309                                 << "\n Updated PSimHit: " << mySimHit->detUnitId() << " " << mySimHit->trackId() << " "
0310                                 << mySimHit->energyLoss() << " " << mySimHit->entryPoint() << " "
0311                                 << mySimHit->exitPoint();
0312   }
0313 }
0314
0315 bool TkAccumulatingSensitiveDetector::newHit(const G4Step* aStep) {
0316   const G4Track* theTrack = aStep->GetTrack();
0317
0318   // for neutral particles do not merge hits (V.I.)
0319   if (0.0 == theTrack->GetDefinition()->GetPDGCharge())
0320     return true;
0321
0322   uint32_t theDetUnitId = setDetUnitId(aStep);
0323   int theTrackID = theTrack->GetTrackID();
0324
0325   LogDebug("TrackerSimDebug") << "newHit: OLD(detID,trID) = (" << lastId << "," << lastTrack << "), NEW = ("
0326                               << theDetUnitId << "," << theTrackID << ") Step length(mm)= " << aStep->GetStepLength()
0327                               << " Edep= " << aStep->GetTotalEnergyDeposit()
0328                               << " p= " << aStep->GetPreStepPoint()->GetMomentum().mag();
0329   return ((theTrackID == lastTrack) && (lastId == theDetUnitId) && closeHit(aStep)) ? false : true;
0330 }
0331
0332 bool TkAccumulatingSensitiveDetector::closeHit(const G4Step* aStep) {
0333   const float tolerance2 = 0.0025f;  // (0.5 mm)^2 are allowed between entry and exit
0334   Local3DPoint theEntryPoint = theRotation.get()->transformPoint(LocalPreStepPosition(aStep));
0335   LogDebug("TrackerSimDebug") << " closeHit: distance = " << (mySimHit->exitPoint() - theEntryPoint).mag();
0336
0337   return ((mySimHit->exitPoint() - theEntryPoint).mag2() < tolerance2) ? true : false;
0338 }
0339
0340 void TkAccumulatingSensitiveDetector::EndOfEvent(G4HCofThisEvent*) {
0341   LogDebug("TrackerSimDebug") << " Saving the last hit in a ROU " << GetName();
0342   if (mySimHit != nullptr)
0343     sendHit();
0344 }
0345
0346 void TkAccumulatingSensitiveDetector::update(const BeginOfEvent* i) {
0347   clearHits();
0348   eventno = (*i)()->GetEventID();
0349   delete mySimHit;
0350   mySimHit = nullptr;
0351 }
0352
0353 void TkAccumulatingSensitiveDetector::update(const BeginOfJob* i) { theNumberingScheme = &(numberingScheme(*pDD_)); }
0354
0355 void TkAccumulatingSensitiveDetector::clearHits() {
0356   slaveLowTof.get()->Initialize();
0357   slaveHighTof.get()->Initialize();
0358 }
0359
0360 void TkAccumulatingSensitiveDetector::fillHits(edm::PSimHitContainer& cc, const std::string& hname) {
0361   if (slaveLowTof.get()->name() == hname) {
0362     cc = slaveLowTof.get()->hits();
0363   } else if (slaveHighTof.get()->name() == hname) {
0364     cc = slaveHighTof.get()->hits();
0365   }
0366 }