Merge pull request #42455 from fabiocos/fc-mtdsimhistory

MTD simulation: update trackId for BTL PSimHits

SimG4CMS/Forward/interface/MtdSD.h

@@ -28,6 +28,10 @@ class MtdSD : public TimingSD {
   int getTrackID(const G4Track *) override;
 
 private:
+  static constexpr unsigned int k_idsecOffset = 100000000;
+  static constexpr unsigned int k_idloopOffset = 200000000;
+  static constexpr unsigned int k_idFromCaloOffset = 300000000;
+
   double energyCut;
   double energyHistoryCut;
 

SimG4CMS/Forward/src/MtdSD.cc

@@ -106,16 +106,23 @@ int MtdSD::getTrackID(const G4Track* aTrack) {
     trkInfo->Print();
 #endif
     if (rname == "FastTimerRegionSensBTL") {
-      if (trkInfo->isBTLdaughter()) {
-        theID = trkInfo->idAtBTLentrance();
+      theID = trkInfo->mcTruthID();
+      if (trkInfo->isExtSecondary() && !trkInfo->isInTrkFromBackscattering()) {
+        theID += k_idsecOffset;
+      } else if (trkInfo->isInTrkFromBackscattering()) {
+        theID += k_idFromCaloOffset;
+      } else if (trkInfo->isBTLlooper()) {
+        theID += k_idloopOffset;
       }
 #ifdef EDM_ML_DEBUG
-      edm::LogVerbatim("MtdSim") << "BTL Track ID: " << trkInfo->idAtBTLentrance() << ":" << theID;
+      edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID()
+                                 << " BTL Track ID: " << trkInfo->mcTruthID() << ":" << theID;
 #endif
     } else if (rname == "FastTimerRegionSensETL") {
       theID = trkInfo->getIDonCaloSurface();
 #ifdef EDM_ML_DEBUG
-      edm::LogVerbatim("MtdSim") << "ETL Track ID: " << trkInfo->getIDonCaloSurface() << ":" << theID;
+      edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID()
+                                 << " ETL Track ID: " << trkInfo->mcTruthID() << ":" << theID;
 #endif
     } else {
       throw cms::Exception("MtdSDError") << "MtdSD called in incorrect region " << rname;

SimG4Core/Application/interface/Phase2SteppingAction.h

@@ -35,6 +35,7 @@ class Phase2SteppingAction : public G4UserSteppingAction {
   bool isLowEnergy(const G4LogicalVolume*, const G4Track*) const;
   void PrintKilledTrack(const G4Track*, const TrackStatus&) const;
 
+  const G4VPhysicalVolume* cmse{nullptr};
   const G4VPhysicalVolume* tracker{nullptr};
   const G4VPhysicalVolume* calo{nullptr};
   const G4VPhysicalVolume* btl{nullptr};

SimG4Core/Application/src/Phase2SteppingAction.cc

@@ -168,17 +168,14 @@ void Phase2SteppingAction::UserSteppingAction(const G4Step* aStep) {
       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();
+                                         << " IdAtBTLentrance = " << trkinfo->mcTruthID();
 #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) {
@@ -196,6 +193,16 @@ void Phase2SteppingAction::UserSteppingAction(const G4Step* aStep) {
       if (!trkinfo->crossedBoundary()) {
         trkinfo->setCrossedBoundary(theTrack);
       }
+    } else if (preStep->GetPhysicalVolume() == calo && postStep->GetPhysicalVolume() == cmse) {
+      // store transition calo -> cmse to tag backscattering
+      TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
+      if (!trkinfo->isInTrkFromBackscattering()) {
+        trkinfo->setInTrkFromBackscattering();
+#ifdef DebugLog
+        LogDebug("SimG4CoreApplication") << "Setting flag for backscattering from CALO "
+                                         << trkinfo->isInTrkFromBackscattering();
+#endif
+      }
     }
   } else {
     theTrack->SetTrackStatus(fStopAndKill);
@@ -236,8 +243,10 @@ bool Phase2SteppingAction::initPointer() {
       calo = pvcite;
     } else if (pvname == "BarrelTimingLayer" || pvname == "btl:BarrelTimingLayer_1") {
       btl = pvcite;
+    } else if (pvname == "CMSE" || pvname == "cms:CMSE_1") {
+      cmse = pvcite;
     }
-    if (tracker && calo && btl)
+    if (tracker && calo && btl && cmse)
       break;
   }
   edm::LogVerbatim("SimG4CoreApplication")

SimG4Core/Notification/interface/TrackInformation.h

@@ -90,13 +90,12 @@ class TrackInformation : public G4VUserTrackInformation {
   bool isFromTtoBTL() const { return (mtdStatus_ >> 0) & 1; }
   void setFromBTLtoT() { mtdStatus_ |= 1 << 1; }  // 2nd bit
   bool isFromBTLtoT() const { return (mtdStatus_ >> 1) & 1; }
-  void setBTLdaughter() { mtdStatus_ |= 1 << 2; }  // 3rd bit
-  bool isBTLdaughter() const { return (mtdStatus_ >> 2) & 1; }
-  void setBTLlooper() { mtdStatus_ |= 1 << 3; }  // 4th bit
-  bool isBTLlooper() const { return (mtdStatus_ >> 3) & 1; }
-
-  int idAtBTLentrance() const { return idAtBTLentrance_; }
-  void setIdAtBTLentrance(int id) { idAtBTLentrance_ = id; }
+  void setBTLlooper() { mtdStatus_ |= 1 << 2; }  // 3th bit
+  bool isBTLlooper() const { return (mtdStatus_ >> 2) & 1; }
+  void setInTrkFromBackscattering() { mtdStatus_ |= 1 << 3; }  // 4th bit
+  bool isInTrkFromBackscattering() const { return (mtdStatus_ >> 3) & 1; }
+  void setExtSecondary() { mtdStatus_ |= 1 << 4; }  //5th bit
+  bool isExtSecondary() const { return (mtdStatus_ >> 4) & 1; }
 
   void Print() const override;
 
@@ -119,7 +118,6 @@ class TrackInformation : public G4VUserTrackInformation {
   int mcTruthID_{-1};
   int caloSurfaceParticlePID_{0};
   int castorHitPID_{0};
-  int idAtBTLentrance_{0};
   uint8_t mtdStatus_{0};
   double genParticleP_{0.};
   double caloSurfaceParticleP_{0.};

SimG4Core/Notification/src/MCTruthUtil.cc

@@ -1,4 +1,3 @@
-
 #include "SimG4Core/Notification/interface/MCTruthUtil.h"
 #include "SimG4Core/Notification/interface/TrackInformation.h"
 
@@ -18,8 +17,6 @@ void MCTruthUtil::primary(G4Track *aTrack) {
 void MCTruthUtil::secondary(G4Track *aTrack, const G4Track &mother, int flag) {
   auto motherInfo = static_cast<const TrackInformation *>(mother.GetUserInformation());
   auto trkInfo = new TrackInformation();
-  LogDebug("SimG4CoreApplication") << "MCTruthUtil called for " << aTrack->GetTrackID() << " mother "
-                                   << motherInfo->isPrimary() << " flag " << flag;
 
   // Take care of cascade decays
   if (flag == 1) {
@@ -57,17 +54,26 @@ void MCTruthUtil::secondary(G4Track *aTrack, const G4Track &mother, int flag) {
     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();
-    }
+  // for MTD
+  if (!trkInfo->isPrimary() && !isInBTL(aTrack)) {
+    trkInfo->setExtSecondary();
+  }
+  if (motherInfo->isExtSecondary()) {
+    trkInfo->setExtSecondary();
+  }
+  if (motherInfo->isBTLlooper()) {
+    trkInfo->setBTLlooper();
+  }
+  if (motherInfo->isInTrkFromBackscattering()) {
+    trkInfo->setInTrkFromBackscattering();
   }
 
   aTrack->SetUserInformation(trkInfo);
+#ifdef EDM_ML_DEBUG
+  LogTrace("SimG4CoreApplication") << "MCTruthUtil called for " << aTrack->GetTrackID() << " mother "
+                                   << motherInfo->isPrimary() << " flag " << flag;
+  trkInfo->Print();
+#endif
 }
 
 bool MCTruthUtil::isInBTL(const G4Track *aTrack) {

SimG4Core/Notification/src/TrackInformation.cc

@@ -26,14 +26,16 @@ void TrackInformation::Print() const {
                                << "                    hasHits = " << hasHits_ << "\n"
                                << "                    isPrimary = " << isPrimary_ << "\n"
                                << "                    isGeneratedSecondary = " << isGeneratedSecondary_ << "\n"
+                               << "                    mcTruthID = " << mcTruthID_ << "\n"
                                << "                    isInHistory = " << isInHistory_ << "\n"
                                << "                    idOnCaloSurface = " << getIDonCaloSurface() << "\n"
                                << "                    caloIDChecked = " << caloIDChecked() << "\n"
                                << "                    idCaloVolume = " << idCaloVolume_ << "\n"
                                << "                    idLastVolume = " << idLastVolume_ << "\n"
                                << "                    isFromTtoBTL = " << isFromTtoBTL() << "\n"
                                << "                    isFromBTLtoT = " << isFromBTLtoT() << "\n"
-                               << "                    isBTLdaughter = " << isBTLdaughter() << "\n"
                                << "                    isBTLlooper = " << isBTLlooper() << "\n"
-                               << "                    idAtBTLentrance = " << idAtBTLentrance_;
+                               << "                    isInTrkFromBackscattering = " << isInTrkFromBackscattering()
+                               << "\n"
+                               << "                    isExtSecondary = " << isExtSecondary();
 }