Revision of the PSimHit type storage, revert trackId to its basic meaning

RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc

@@ -709,9 +709,9 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup&
               });
               // Check if TP has direct or other sim cluster for BTL
               for (const auto& simClusterRef : simClustersRefs) {
-                if (simClusterRef->trackIdOffset() == 0) {
+                if (simClusterRef->hitProdType() == 0) {
                   isTPmtdDirectETL = true;
-                } else if (simClusterRef->trackIdOffset() != 0) {
+                } else if (simClusterRef->hitProdType() != 0) {
                   isTPmtdOtherETL = true;
                 }
               }
@@ -724,9 +724,9 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup&
               });
               // Check if TP has direct or other sim cluster for BTL
               for (const auto& simClusterRef : simClustersRefs) {
-                if (simClusterRef->trackIdOffset() == 0) {
+                if (simClusterRef->hitProdType() == 0) {
                   isTPmtdDirectBTL = true;
-                } else if (simClusterRef->trackIdOffset() != 0) {
+                } else if (simClusterRef->hitProdType() != 0) {
                   isTPmtdOtherBTL = true;
                 }
               }

SimDataFormats/CaloAnalysis/interface/MtdSimCluster.h

@@ -116,9 +116,9 @@ class MtdSimCluster : public SimCluster {
     ++nsimhits_;
   }
 
-  void setTrackIdOffset(unsigned int offset) { idOffset_ = offset; }
+  void setHitProdType(unsigned int offset) { idOffset_ = offset; }
 
-  unsigned int trackIdOffset() const { return idOffset_; }
+  unsigned int hitProdType() const { return idOffset_; }
 
 protected:
   std::vector<uint64_t> mtdHits_;

SimDataFormats/CaloAnalysis/src/MtdSimLayerCluster.cc

@@ -27,7 +27,7 @@ MtdSimLayerCluster::~MtdSimLayerCluster() {}
 std::ostream &operator<<(std::ostream &s, MtdSimLayerCluster const &tp) {
   s << "CP momentum, q, ID, & Event #: " << tp.p4() << " " << tp.charge() << " " << tp.pdgId() << " "
     << tp.eventId().bunchCrossing() << "." << tp.eventId().event() << std::endl;
-  s << " Offset " << tp.trackIdOffset() << " "
+  s << " Offset " << tp.hitProdType() << " "
     << " LC time " << tp.simLCTime() << " LC energy " << tp.simLCEnergy() << std::endl;
 
   for (MtdSimLayerCluster::genp_iterator hepT = tp.genParticle_begin(); hepT != tp.genParticle_end(); ++hepT) {

SimDataFormats/TrackingHit/interface/PSimHit.h

@@ -14,8 +14,6 @@ class TrackingSlaveSD;  // for friend declaration only
 
 class PSimHit {
 public:
-  static constexpr unsigned int k_tidOffset = 200000000;
-
   PSimHit() : theDetUnitId(0) {}
 
   PSimHit(const Local3DPoint& entry,
@@ -107,15 +105,6 @@ class PSimHit {
    */
   unsigned int trackId() const { return theTrackId; }
 
-  /** In case te SimTrack ID is incremented by the k_tidOffset for hit category definition, this
-   * methods returns the original theTrackId value directly.
-   */
-  unsigned int originalTrackId() const { return theTrackId % k_tidOffset; }
-
-  unsigned int offsetTrackId() const { return theTrackId / k_tidOffset; }
-
-  static unsigned int addTrackIdOffset(unsigned int tId, unsigned int offset) { return offset * k_tidOffset + tId; }
-
   EncodedEventId eventId() const { return theEventId; }
 
   void setEventId(EncodedEventId e) { theEventId = e; }
@@ -128,10 +117,20 @@ class PSimHit {
    *  value with special significance is zero (for "undefined"), so zero should
    *  not be the ID of any process.
    */
-  unsigned short processType() const { return theProcessType; }
+
+  // use 9 bits (up to 511) for process id, reserve the rest for hit production mechanism id
+  // 7 bits field available in PSimHit processType, i.e. up 127, to store processes
+  unsigned short processType() const { return theProcessType & kProcidMask; }
+
+  unsigned short hitProdType() const { return (theProcessType >> kHitidShift) & kHitidMask; }
+  void setHitProdType(unsigned int hitId) { theProcessType |= hitId << kHitidShift; }
 
   void setTof(float tof) { theTof = tof; }
 
+  static constexpr unsigned int kProcidMask = 0x1FF;
+  static constexpr unsigned int kHitidMask = 0x7F;
+  static constexpr unsigned int kHitidShift = 9;
+
 protected:
   // properties
   Local3DPoint theEntryPoint;  // position at entry

SimDataFormats/TrackingHit/interface/SimHitCategory.h

@@ -0,0 +1,22 @@
+#ifndef SimHitCategory_h
+#define SimHitCategory_h
+
+namespace SimHitCategory {
+
+  // Identification code of sim hit production mechanism, subdetector dependent
+
+  // MTD
+
+  static constexpr unsigned int nCategoriesMTD = 5;
+
+  static constexpr unsigned int prodTypeMTD[nCategoriesMTD] = {
+      0,  // direct hit from particle coming from tracker
+      1,  // BTL hit from secondary particle not saved in history
+      2,  // BTL hit from identified looper
+      3,  // BTL hit from back-scattering from CALO volume
+      4   // ETL hit entering from a rear face of disks
+  };
+
+};  // namespace SimHitCategory
+
+#endif

SimDataFormats/TrackingHit/src/classes_def.xml

@@ -1,5 +1,6 @@
 <lcgdict>
-  <class name="PSimHit" ClassVersion="10">
+  <class name="PSimHit" ClassVersion="11">
+   <version ClassVersion="11" checksum="1181355938"/>
    <version ClassVersion="10" checksum="1181355938"/>
   </class>
   <class name="std::vector<PSimHit>"/>

SimDataFormats/TrackingHit/test/BuildFile.xml

@@ -0,0 +1,4 @@
+<use name="SimDataFormats/TrackingHit"/>
+
+<bin file="testPSimHit.cc">
+</bin>

SimDataFormats/TrackingHit/test/testPSimHit.cc

@@ -0,0 +1,24 @@
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+
+#include <iostream>
+#include <iomanip>
+#include <cassert>
+
+int main() {
+  Local3DPoint dummy(0, 0, 0);
+
+  for (int procId = 1; procId < 404; procId++) {
+    for (unsigned short itype = 0; itype < 8; itype++) {
+      PSimHit testHit(dummy, dummy, 0., 0., 0., 11, 0, 0, 0., 0., procId);
+      testHit.setHitProdType(itype);
+      std::cout << " hit procId = " << std::fixed << std::setw(8) << procId << " sim procId = " << std::setw(8)
+                << testHit.processType() << " hit type = " << std::setw(8) << itype << " sim type = " << std::setw(8)
+                << testHit.hitProdType() << std::endl;
+
+      assert(procId == testHit.processType());
+      assert(itype == testHit.hitProdType());
+    }
+  }
+
+  return 0;
+}

SimG4CMS/Forward/interface/BscG4Hit.h

@@ -14,6 +14,8 @@
 #include <cstdint>
 #include <iostream>
 
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+
 class BscG4Hit : public G4VHit {
 public:
   BscG4Hit();
@@ -92,6 +94,8 @@ class BscG4Hit : public G4VHit {
 
   int getParentId() const { return theParentId; };
   int getProcessId() const { return theProcessId; };
+  int getProdType() const { return (theProcessId >> PSimHit::kHitidShift) & PSimHit::kHitidMask; }
+  void setHitProdType(unsigned int hitId) { theProcessId |= hitId << PSimHit::kHitidShift; }
   float getVx() const { return theVx; };
   float getVy() const { return theVy; };
   float getVz() const { return theVz; };

SimG4CMS/Forward/interface/MtdSD.h

@@ -2,7 +2,7 @@
 #define SimG4CMSForward_MtdSD_h
 
 #include "SimG4CMS/Forward/interface/TimingSD.h"
-#include "SimG4CMS/Forward/interface/MtdHitCategory.h"
+#include "SimDataFormats/TrackingHit/interface/SimHitCategory.h"
 
 #include "FWCore/Framework/interface/Frameworkfwd.h"
 #include "FWCore/ParameterSet/interface/ParameterSetfwd.h"

SimG4CMS/Forward/src/BSCG4Hit.cc

@@ -5,6 +5,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "SimG4CMS/Forward/interface/BscG4Hit.h"
+#include "SimDataFormats/TrackingHit/interface/SimHitCategory.h"
 #include <iostream>
 
 BscG4Hit::BscG4Hit() : entry(0., 0., 0.), entrylp(0., 0., 0.), exitlp(0., 0., 0.) {
@@ -141,8 +142,8 @@ std::ostream& operator<<(std::ostream& os, const BscG4Hit& hit) {
      << " EnergyLoss = " << hit.getEnergyLoss() << std::endl
      << " ParticleType = " << hit.getParticleType() << std::endl
      << " Pabs = " << hit.getPabs() << std::endl
-     << " Energy of primary particle (ID = " << hit.getTrackID() << ") = " << hit.getIncidentEnergy() << " (MeV)"
-     << std::endl
+     << " Energy of primary particle (ID = " << hit.getTrackID() << " hit type = " << hit.getProdType()
+     << ") = " << hit.getIncidentEnergy() << " (MeV)" << std::endl
      << " Entry point in Bsc unit number " << hit.getUnitID() << " is: " << hit.getEntry() << " (mm)" << std::endl;
   os << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
   return os;

SimG4CMS/Forward/src/MtdSD.cc

@@ -16,8 +16,6 @@
 
 #include <iostream>
 
-using namespace MtdHitCategory;
-
 //-------------------------------------------------------------------
 MtdSD::MtdSD(const std::string& name,
              const SensitiveDetectorCatalog& clg,
@@ -111,32 +109,12 @@ int MtdSD::getTrackID(const G4Track* aTrack) {
     if (!trkInfo->storeTrack()) {
       theID = trkInfo->idLastStoredAncestor();
     }
-    if (theID >= static_cast<int>(PSimHit::k_tidOffset)) {
-      edm::LogError("MtdSim") << " SimTrack ID " << theID << " exceeds maximum allowed by PSimHit identifier"
-                              << PSimHit::k_tidOffset << " unreliable MTD hit type";
-    }
-    if (rname == "FastTimerRegionSensBTL") {
-      if (trkInfo->isInTrkFromBackscattering()) {
-        theID = PSimHit::addTrackIdOffset(theID, k_idFromCaloOffset);
-      } else if (trkInfo->isExtSecondary() && !trkInfo->isInTrkFromBackscattering() && !trkInfo->storeTrack()) {
-        theID = PSimHit::addTrackIdOffset(theID, k_idsecOffset);
-      } else if (trkInfo->isBTLlooper()) {
-        theID = PSimHit::addTrackIdOffset(theID, k_idloopOffset);
-      }
-#ifdef EDM_ML_DEBUG
-      edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID()
-                                 << " BTL Track ID: " << trkInfo->mcTruthID() << ":" << theID;
-#endif
-    } else if (rname == "FastTimerRegionSensETL") {
-      if (hitClassID == k_idETLfromBack) {
-        theID = PSimHit::addTrackIdOffset(theID, k_idETLfromBack);
-      }
 #ifdef EDM_ML_DEBUG
-      edm::LogVerbatim("MtdSim") << "MtdSD: Track ID: " << aTrack->GetTrackID()
-                                 << " ETL Track ID: " << trkInfo->mcTruthID() << ":" << theID;
+    edm::LogVerbatim("MtdSim") << "MtdSD: current Track ID: " << aTrack->GetTrackID()
+                               << " stored Track ID: " << trkInfo->mcTruthID() << ":" << theID;
 #endif
-      // In the case of ECAL GFlash fast spot may be inside MTD and should be ignored
-    } else {
+    // In the case of ECAL GFlash fast spot may be inside MTD and should be ignored
+    if (rname != "FastTimerRegionSensBTL" && rname != "FastTimerRegionSensETL") {
       throw cms::Exception("MtdSDError") << "MtdSD called in incorrect region " << rname;
     }
   } else {
@@ -148,20 +126,31 @@ int MtdSD::getTrackID(const G4Track* aTrack) {
 }
 
 void MtdSD::setHitClassID(const G4Step* aStep) {
+  hitClassID = 0;
   TrackInformation* trkInfo = cmsTrackInformation(aStep->GetTrack());
   if (nullptr == trkInfo) {
     return;
   }
   const G4String& rname = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetRegion()->GetName();
-  if (rname == "FastTimerRegionSensETL") {
+  if (rname == "FastTimerRegionSensBTL") {
+    if (trkInfo->isInTrkFromBackscattering()) {
+      hitClassID = SimHitCategory::prodTypeMTD[3];
+    } else if (trkInfo->isExtSecondary() && !trkInfo->isInTrkFromBackscattering() && !trkInfo->storeTrack()) {
+      hitClassID = SimHitCategory::prodTypeMTD[1];
+    } else if (trkInfo->isBTLlooper()) {
+      hitClassID = SimHitCategory::prodTypeMTD[2];
+    }
+  } else if (rname == "FastTimerRegionSensETL") {
     double zin = std::abs(aStep->GetPreStepPoint()->GetPosition().z());
     double zout = std::abs(aStep->GetPostStepPoint()->GetPosition().z());
     if (zout - zin < 0.) {
-      hitClassID = k_idETLfromBack;
+      hitClassID = SimHitCategory::prodTypeMTD[4];
       trkInfo->setETLfromBack();
     } else {
-      hitClassID = 0;
       trkInfo->setETLfromFront();
     }
   }
+#ifdef EDM_ML_DEBUG
+  edm::LogVerbatim("MtdSim") << "MtdSD: process type = " << hitClassID;
+#endif
 }

SimG4CMS/Forward/src/TimingSD.cc

@@ -238,6 +238,7 @@ bool TimingSD::checkHit(const G4Step*, BscG4Hit* hit) {
 
     hit->setParentId(theTrack->GetParentID());
     hit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess()));
+    hit->setHitProdType(hitClassID);
 
     hit->setVertexPosition(theTrack->GetVertexPosition());
   }
@@ -294,6 +295,7 @@ void TimingSD::createNewHit(const G4Step* aStep) {
 
   currentHit->setParentId(theTrack->GetParentID());
   currentHit->setProcessId(theEnumerator->processId(theTrack->GetCreatorProcess()));
+  currentHit->setHitProdType(hitClassID);
 
   currentHit->setVertexPosition(theTrack->GetVertexPosition());
 

SimG4Core/Application/test/SimHitCaloHitDumper.cc

@@ -373,8 +373,7 @@ void SimHitCaloHitDumper::analyze(const edm::Event& iEvent, const edm::EventSetu
     for (int ihit = 0; ihit < (*icoll).first; ++ihit) {
       edm::LogPrint("SimHitCaloHitDumper")
           << theMTDHits[nhit] << " Energy = " << theMTDHits[nhit].energyLoss()
-          << " tid orig/offset= " << theMTDHits[nhit].originalTrackId() << " " << theMTDHits[nhit].offsetTrackId()
-          << " Track Id = " << theMTDHits[nhit].trackId();
+          << " Track Id = " << theMTDHits[nhit].trackId() << " hit type = " << theMTDHits[nhit].hitProdType();
       nhit++;
     }
   }

SimG4Core/Notification/src/SimTrackManager.cc

@@ -175,10 +175,6 @@ void SimTrackManager::reallyStoreTracks() {
       }
     }
 
-    if (id >= static_cast<int>(PSimHit::k_tidOffset)) {
-      edm::LogWarning("SimTrackManager::reallyStoreTracks")
-          << " SimTrack ID " << id << " exceeds maximum allowed by PSimHit identifier" << PSimHit::k_tidOffset;
-    }
     TmpSimTrack* g4simtrack =
         new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom);
     g4simtrack->copyCrossedBoundaryVars(trkH);