Merge pull request #41318 from fabiocos/fc-g4ForBTL

MTD simulation: update Track ID stored in MTD SimHits

Geometry/MTDNumberingBuilder/src/CmsMTDStringToEnum.cc

@@ -4,6 +4,8 @@ const CmsMTDStringToEnum::Impl CmsMTDStringToEnum::m_impl;
 
 CmsMTDStringToEnum::Impl::Impl() {
   map_.emplace("FastTimerRegion", GeometricTimingDet::MTD);
+  map_.emplace("FastTimerRegionBTL", GeometricTimingDet::MTD);
+  map_.emplace("FastTimerRegionETL", GeometricTimingDet::MTD);
   map_.emplace("BarrelTimingLayer", GeometricTimingDet::BTL);
   map_.emplace("Layer1", GeometricTimingDet::BTLLayer);
   map_.emplace("Layer1Timing", GeometricTimingDet::BTLLayer);

Geometry/MTDSimData/data/v2/mtdProdCuts.xml

@@ -3,10 +3,17 @@
 <DDDefinition xmlns="http://www.cern.ch/cms/DDL" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.cern.ch/cms/DDL ../../../DetectorDescription/Schema/DDLSchema.xsd">
 
   <SpecParSection label="mtdProdCuts.xml" eval="true">
-    <SpecPar name="fasttime">
+    <SpecPar name="fasttime-btl">
       <PartSelector path="//BarrelTimingLayer"/>
+      <Parameter name="CMSCutsRegion" value="FastTimerRegionBTL" eval="false"/>
+      <Parameter name="ProdCutsForElectrons" value="0.1*mm"/>
+      <Parameter name="ProdCutsForPositrons" value="0.1*mm"/>
+      <Parameter name="ProdCutsForGamma" value="0.1*mm"/>
+      <Parameter name="ProdCutsForProtons" value="0.1*mm"/>
+    </SpecPar>
+    <SpecPar name="fasttime-etl">
       <PartSelector path="//EndcapTimingLayer"/>
-      <Parameter name="CMSCutsRegion" value="FastTimerRegion" eval="false"/>
+      <Parameter name="CMSCutsRegion" value="FastTimerRegionETL" eval="false"/>
       <Parameter name="ProdCutsForElectrons" value="0.1*mm"/>
       <Parameter name="ProdCutsForPositrons" value="0.1*mm"/>
       <Parameter name="ProdCutsForGamma" value="0.1*mm"/>

Geometry/MTDSimData/data/v4/mtdProdCuts.xml

@@ -3,11 +3,18 @@
 <DDDefinition>
 
   <SpecParSection label="mtdProdCuts.xml" eval="true">
-    <SpecPar name="fasttime">
+    <SpecPar name="fasttime-btl">
       <PartSelector path="//BarrelTimingLayer"/>
+      <Parameter name="CMSCutsRegion" value="FastTimerRegionBTL" eval="false"/>
+      <Parameter name="ProdCutsForElectrons" value="0.1*mm"/>
+      <Parameter name="ProdCutsForPositrons" value="0.1*mm"/>
+      <Parameter name="ProdCutsForGamma" value="0.1*mm"/>
+      <Parameter name="ProdCutsForProtons" value="0.1*mm"/>
+    </SpecPar>
+    <SpecPar name="fasttime-etl">
       <PartSelector path="//EndcapTimingLayer"/>
       <PartSelector path="//CALOECTSFront"/>
-      <Parameter name="CMSCutsRegion" value="FastTimerRegion" eval="false"/>
+      <Parameter name="CMSCutsRegion" value="FastTimerRegionETL" eval="false"/>
       <Parameter name="ProdCutsForElectrons" value="0.1*mm"/>
       <Parameter name="ProdCutsForPositrons" value="0.1*mm"/>
       <Parameter name="ProdCutsForGamma" value="0.1*mm"/>

SimG4CMS/Forward/interface/MtdSD.h

@@ -24,7 +24,13 @@ class MtdSD : public TimingSD {
 
   uint32_t setDetUnitId(const G4Step *) override;
 
+protected:
+  int getTrackID(const G4Track *) override;
+
 private:
+  double energyCut;
+  double energyHistoryCut;
+
   void setNumberingScheme(MTDNumberingScheme *);
   void getBaseNumber(const G4Step *);
 

SimG4CMS/Forward/interface/TimingSD.h

@@ -62,6 +62,9 @@ class TimingSD : public SensitiveTkDetector, public Observer<const BeginOfEvent*
   const G4ThreeVector& getLocalEntryPoint() const { return hitPointLocal; };
   const G4ThreeVector& getGlobalEntryPoint() const { return hitPoint; };
 
+  // general method to assign track ID to be stored in hits
+  virtual int getTrackID(const G4Track*);
+
 private:
   void getStepInfo(const G4Step*);
   bool hitExists(const G4Step*);

SimG4CMS/Forward/src/MtdSD.cc

@@ -1,3 +1,5 @@
+//#define EDM_ML_DEBUG
+
 #include "SimG4CMS/Forward/interface/MtdSD.h"
 
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
@@ -13,7 +15,6 @@
 
 #include <iostream>
 
-//#define EDM_ML_DEBUG
 //-------------------------------------------------------------------
 MtdSD::MtdSD(const std::string& name,
              const SensitiveDetectorCatalog& clg,
@@ -40,6 +41,11 @@ MtdSD::MtdSD(const std::string& name,
   if (scheme)
     setNumberingScheme(scheme);
 
+  energyCut = m_p.getParameter<double>("EnergyThresholdForPersistencyInGeV") * CLHEP::GeV;     //default must be 0.5
+  energyHistoryCut = m_p.getParameter<double>("EnergyThresholdForHistoryInGeV") * CLHEP::GeV;  //default must be 0.05
+
+  setCuts(energyCut, energyHistoryCut);
+
   double newTimeFactor = 1. / m_p.getParameter<double>("TimeSliceUnit");
   edm::LogVerbatim("MtdSim") << "New time factor = " << newTimeFactor;
   setTimeFactor(newTimeFactor);
@@ -90,3 +96,34 @@ void MtdSD::getBaseNumber(const G4Step* aStep) {
     }
   }
 }
+
+int MtdSD::getTrackID(const G4Track* aTrack) {
+  int theID = aTrack->GetTrackID();
+  TrackInformation* trkInfo = cmsTrackInformation(aTrack);
+  const G4String& rname = aTrack->GetVolume()->GetLogicalVolume()->GetRegion()->GetName();
+  if (trkInfo != nullptr) {
+#ifdef EDM_ML_DEBUG
+    trkInfo->Print();
+#endif
+    if (rname == "FastTimerRegionSensBTL") {
+      if (trkInfo->isBTLdaughter()) {
+        theID = trkInfo->idAtBTLentrance();
+      }
+#ifdef EDM_ML_DEBUG
+      edm::LogVerbatim("MtdSim") << "BTL Track ID: " << trkInfo->idAtBTLentrance() << ":" << theID;
+#endif
+    } else if (rname == "FastTimerRegionSensETL") {
+      theID = trkInfo->getIDonCaloSurface();
+#ifdef EDM_ML_DEBUG
+      edm::LogVerbatim("MtdSim") << "ETL Track ID: " << trkInfo->getIDonCaloSurface() << ":" << theID;
+#endif
+    } else {
+      throw cms::Exception("MtdSDError") << "MtdSD called in incorrect region " << rname;
+    }
+  } else {
+#ifdef EDM_ML_DEBUG
+    edm::LogWarning("MtdSim") << "MtdSD: Problem with primaryID **** set by force to TkID **** " << theID;
+#endif
+  }
+  return theID;
+}

SimG4CMS/Forward/src/TimingSD.cc

@@ -1,3 +1,5 @@
+//#define EDM_ML_DEBUG
+
 ///////////////////////////////////////////////////////////////////////////////
 // File: TimingSD.cc
 // Date: 02.2006
@@ -27,8 +29,6 @@
 #include <vector>
 #include <iostream>
 
-//#define EDM_ML_DEBUG
-
 static const float invgev = 1.0 / CLHEP::GeV;
 static const double invns = 1.0 / CLHEP::nanosecond;
 static const double invdeg = 1.0 / CLHEP::deg;
@@ -161,6 +161,12 @@ void TimingSD::getStepInfo(const G4Step* aStep) {
       }
       info->putInHistory();
     }
+#ifdef EDM_ML_DEBUG
+    if (info != nullptr) {
+      LogDebug("TimingSim") << "TrackInformation for ID = " << theTrack->GetTrackID();
+      info->Print();
+    }
+#endif
   }
 
   edeposit *= invgev;
@@ -176,7 +182,7 @@ void TimingSD::getStepInfo(const G4Step* aStep) {
   tSliceID = (int)tSlice;
 
   unitID = setDetUnitId(aStep);
-  primaryID = theTrack->GetTrackID();
+  primaryID = getTrackID(theTrack);
 }
 
 bool TimingSD::hitExists(const G4Step* aStep) {
@@ -255,7 +261,7 @@ void TimingSD::createNewHit(const G4Step* aStep) {
   edm::LogVerbatim("TimingSim") << "TimingSD CreateNewHit for " << GetName() << " PV " << currentPV->GetName()
                                 << " PVid = " << currentPV->GetCopyNo() << " Unit " << unitID << "\n primary "
                                 << primaryID << " Tof(ns)= " << tof << " time slice " << tSliceID
-                                << " E(GeV)= " << incidentEnergy << " trackID " << theTrack->GetTrackID() << " "
+                                << " E(MeV)= " << incidentEnergy << " trackID " << theTrack->GetTrackID() << " "
                                 << theTrack->GetDefinition()->GetParticleName() << " parentID "
                                 << theTrack->GetParentID();
 
@@ -365,3 +371,8 @@ void TimingSD::update(const BeginOfEvent* i) {
 }
 
 void TimingSD::clearHits() { slave->Initialize(); }
+
+int TimingSD::getTrackID(const G4Track* aTrack) {
+  LogDebug("TimingSim") << "primary ID: " << aTrack->GetTrackID();
+  return aTrack->GetTrackID();
+}

SimG4Core/Application/interface/SteppingAction.h

@@ -48,6 +48,7 @@ class SteppingAction : public G4UserSteppingAction {
 
   const G4VPhysicalVolume* tracker{nullptr};
   const G4VPhysicalVolume* calo{nullptr};
+  const G4VPhysicalVolume* btl{nullptr};
   const CMSSteppingVerbose* steppingVerbose;
   double theCriticalEnergyForVacuum;
   double theCriticalDensity;

SimG4Core/Application/python/g4SimHits_cfi.py

@@ -518,6 +518,8 @@
         IgnoreTrackID    = cms.bool(False),
         EminHit          = cms.double(0.0),
         CheckID          = cms.untracked.bool(True),
+        EnergyThresholdForPersistencyInGeV = cms.double(1e9), # keep temporarily old behaviour
+        EnergyThresholdForHistoryInGeV = cms.double(1e9) # keep temporarily old behaviour)
     ),
     HGCSD = cms.PSet(
         Verbosity        = cms.untracked.int32(0),

SimG4Core/Application/src/StackingAction.cc

@@ -418,7 +418,8 @@ void StackingAction::initPointer() {
     if (savePDandCinTracker &&
         (rname == "BeamPipe" || rname == "BeamPipeVacuum" || rname == "TrackerPixelSensRegion" ||
          rname == "TrackerPixelDeadRegion" || rname == "TrackerDeadRegion" || rname == "TrackerSensRegion" ||
-         rname == "FastTimerRegion" || rname == "FastTimerRegionSensBTL" || rname == "FastTimerRegionSensETL")) {
+         rname == "FastTimerRegionBTL" || rname == "FastTimerRegionETL" || rname == "FastTimerRegionSensBTL" ||
+         rname == "FastTimerRegionSensETL")) {
       trackerRegions.push_back(reg);
     }
     if (savePDandCinCalo && (rname == "HcalRegion" || rname == "EcalRegion" || rname == "PreshowerSensRegion" ||

SimG4Core/Application/src/SteppingAction.cc

@@ -1,3 +1,5 @@
+//#define EDM_ML_DEBUG
+
 #include "SimG4Core/Application/interface/SteppingAction.h"
 
 #include "SimG4Core/Notification/interface/TrackInformation.h"
@@ -158,10 +160,39 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
       }
     }
   }
-  // check transition tracker/calo
+  // check transition tracker/btl and tracker/calo
   bool isKilled = false;
   if (sAlive == tstat || sVeryForward == tstat) {
-    if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == calo) {
+    // store TrackInformation about transition from one envelope to another
+    if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == btl) {
+      // store transition tracker -> BTL only for tracks entering BTL for the first time
+      TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
+      if (!trkinfo->isFromTtoBTL() && !trkinfo->isFromBTLtoT()) {
+        trkinfo->setFromTtoBTL();
+        trkinfo->setIdAtBTLentrance(theTrack->GetTrackID());
+#ifdef DebugLog
+        LogDebug("SimG4CoreApplication") << "Setting flag for Tracker -> BTL " << trkinfo->isFromTtoBTL()
+                                         << " IdAtBTLentrance = " << trkinfo->idAtBTLentrance();
+#endif
+      } else {
+        trkinfo->setBTLlooper();
+        trkinfo->setIdAtBTLentrance(theTrack->GetTrackID());
+#ifdef DebugLog
+        LogDebug("SimG4CoreApplication") << "Setting flag for BTL looper " << trkinfo->isBTLlooper();
+        trkinfo->Print();
+#endif
+      }
+    } else if (preStep->GetPhysicalVolume() == btl && postStep->GetPhysicalVolume() == tracker) {
+      // store transition BTL -> tracker
+      TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
+      if (!trkinfo->isFromBTLtoT()) {
+        trkinfo->setFromBTLtoT();
+#ifdef DebugLog
+        LogDebug("SimG4CoreApplication") << "Setting flag for BTL -> Tracker " << trkinfo->isFromBTLtoT();
+#endif
+      }
+    } else if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == calo) {
+      // store transition tracker -> calo
       TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
       if (!trkinfo->crossedBoundary()) {
         trkinfo->setCrossedBoundary(theTrack);
@@ -204,12 +235,14 @@ bool SteppingAction::initPointer() {
       tracker = pvcite;
     } else if (pvname == "CALO" || pvname == "caloBase:CALO_1") {
       calo = pvcite;
+    } else if (pvname == "BarrelTimingLayer" || pvname == "btl:BarrelTimingLayer_1") {
+      btl = pvcite;
     }
-    if (tracker && calo)
+    if (tracker && calo && btl)
       break;
   }
   edm::LogVerbatim("SimG4CoreApplication")
-      << "SteppingAction: pointer for Tracker " << tracker << " and for Calo " << calo;
+      << "SteppingAction: pointer for Tracker " << tracker << " and for Calo " << calo << " and for BTL " << btl;
 
   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
   if (numberEkins > 0) {

SimG4Core/Application/src/TrackingAction.cc

@@ -1,3 +1,5 @@
+//#define EDM_ML_DEBUG
+
 #include "SimG4Core/Application/interface/TrackingAction.h"
 
 #include "SimG4Core/Notification/interface/NewTrackAction.h"
@@ -15,8 +17,6 @@
 #include "G4TrackingManager.hh"
 #include "G4SystemOfUnits.hh"
 
-//#define EDM_ML_DEBUG
-
 TrackingAction::TrackingAction(SimTrackManager* stm, CMSSteppingVerbose* sv, const edm::ParameterSet& p)
     : trackManager_(stm),
       steppingVerbose_(sv),

SimG4Core/Notification/interface/NewTrackAction.h

@@ -21,6 +21,8 @@ class NewTrackAction {
 private:
   void addUserInfoToPrimary(G4Track* aTrack) const;
   void addUserInfoToSecondary(G4Track* aTrack, const TrackInformation& motherInfo, int) const;
+
+  bool isInBTL(const G4Track*) const;
 };
 
 #endif

SimG4Core/Notification/interface/TrackInformation.h

@@ -80,6 +80,20 @@ class TrackInformation : public G4VUserTrackInformation {
   }
   int getCastorHitPID() const { return castorHitPID_; }
 
+  // methods for MTD info management
+  //
+  inline void setFromTtoBTL() { mtdStatus_ |= 1 << 0; }  // 1st bit
+  inline bool isFromTtoBTL() const { return (mtdStatus_ >> 0) & 1; }
+  inline void setFromBTLtoT() { mtdStatus_ |= 1 << 1; }  // 2nd bit
+  inline bool isFromBTLtoT() const { return (mtdStatus_ >> 1) & 1; }
+  inline void setBTLdaughter() { mtdStatus_ |= 1 << 2; }  // 3rd bit
+  inline bool isBTLdaughter() const { return (mtdStatus_ >> 2) & 1; }
+  inline void setBTLlooper() { mtdStatus_ |= 1 << 3; }  // 4th bit
+  inline bool isBTLlooper() const { return (mtdStatus_ >> 3) & 1; }
+
+  int idAtBTLentrance() const { return idAtBTLentrance_; }
+  void setIdAtBTLentrance(int id) { idAtBTLentrance_ = id; }
+
   void Print() const override;
 
 private:
@@ -105,6 +119,9 @@ class TrackInformation : public G4VUserTrackInformation {
   bool hasCastorHit_;
   int castorHitPID_;
 
+  uint8_t mtdStatus_;
+  int idAtBTLentrance_;
+
   // Restrict construction to friends
   TrackInformation()
       : G4VUserTrackInformation(),
@@ -126,7 +143,9 @@ class TrackInformation : public G4VUserTrackInformation {
         genParticleP_(0),
         caloSurfaceParticleP_(0),
         hasCastorHit_(false),
-        castorHitPID_(0) {}
+        castorHitPID_(0),
+        mtdStatus_(0),
+        idAtBTLentrance_(0) {}
   friend class NewTrackAction;
 };
 

SimG4Core/Notification/src/NewTrackAction.cc

@@ -1,3 +1,4 @@
+//#define EDM_ML_DEBUG
 
 #include "SimG4Core/Notification/interface/NewTrackAction.h"
 #include "SimG4Core/Notification/interface/TrackInformation.h"
@@ -51,18 +52,41 @@ void NewTrackAction::addUserInfoToSecondary(G4Track *aTrack, const TrackInformat
                                 motherInfo.getIDLastVolume(),
                                 aTrack->GetDefinition()->GetPDGEncoding(),
                                 aTrack->GetMomentum().mag());
+    LogDebug("SimG4CoreApplication") << "NewTrackAction: Id on calo surface " << trkInfo->getIDonCaloSurface()
+                                     << " to be stored " << trkInfo->storeTrack();
   } else {
     // transfer calo ID from mother (to be checked in TrackingAction)
     trkInfo->setIDonCaloSurface(motherInfo.getIDonCaloSurface(),
                                 motherInfo.getIDCaloVolume(),
                                 motherInfo.getIDLastVolume(),
                                 motherInfo.caloSurfaceParticlePID(),
                                 motherInfo.caloSurfaceParticleP());
+    LogDebug("SimG4CoreApplication") << "NewTrackAction: Id on calo surface " << trkInfo->getIDonCaloSurface();
   }
 
   if (motherInfo.hasCastorHit()) {
     trkInfo->setCastorHitPID(motherInfo.getCastorHitPID());
   }
 
+  // manage ID of tracks in BTL to map them to SimTracks to be stored
+  if (isInBTL(aTrack)) {
+    if ((motherInfo.storeTrack() && motherInfo.isFromTtoBTL()) || motherInfo.isBTLdaughter()) {
+      trkInfo->setBTLdaughter();
+      trkInfo->setIdAtBTLentrance(motherInfo.idAtBTLentrance());
+      LogDebug("SimG4CoreApplication") << "NewTrackAction: secondary in BTL " << trkInfo->isBTLdaughter()
+                                       << " from mother ID " << trkInfo->idAtBTLentrance();
+    }
+  }
+
   aTrack->SetUserInformation(trkInfo);
 }
+
+bool NewTrackAction::isInBTL(const G4Track *aTrack) const {
+  bool out = false;
+  G4String tName(aTrack->GetVolume()->GetLogicalVolume()->GetRegion()->GetName());
+  if (tName == "FastTimerRegionBTL" || tName == "FastTimerRegionSensBTL") {
+    out = true;
+  }
+
+  return out;
+}