Add additional info to the simTrack

- is or not from backscattering
- is or not a primary
- the GenParticle ID is always saved if available and returned with a
specific method, while the old genpartIndex() returns -1 if it's not a
primary as before

Author: AuroraPerego

SimDataFormats/Track/interface/SimTrack.h

@@ -34,8 +34,8 @@ class SimTrack : public CoreSimTrack {
   bool noVertex() const { return ivert == -1; }
 
   /// index of the corresponding Generator particle in the Event container (-1 if no Genpart)
-  int genpartIndex() const { return igenpart; }
-  bool noGenpart() const { return igenpart == -1; }
+  int genpartIndex() const { return isPrimary() ? igenpart : -1; }
+  bool noGenpart() const { return isPrimary() ? igenpart == -1 : true; }
 
   const math::XYZVectorD& trackerSurfacePosition() const { return tkposition; }
 
@@ -51,27 +51,41 @@ class SimTrack : public CoreSimTrack {
                               int idAtBoundary,
                               math::XYZTLorentzVectorF positionAtBoundary,
                               math::XYZTLorentzVectorF momentumAtBoundary) {
-    crossedBoundary_ = crossedBoundary;
+    if (crossedBoundary)
+      trackInfo_ |= 1 << 2;
     idAtBoundary_ = idAtBoundary;
     positionAtBoundary_ = positionAtBoundary;
     momentumAtBoundary_ = momentumAtBoundary;
   }
-  bool crossedBoundary() const { return crossedBoundary_; }
+  bool crossedBoundary() const { return (trackInfo_ >> 2) & 1; }
   const math::XYZTLorentzVectorF& getPositionAtBoundary() const { return positionAtBoundary_; }
   const math::XYZTLorentzVectorF& getMomentumAtBoundary() const { return momentumAtBoundary_; }
   int getIDAtBoundary() const { return idAtBoundary_; }
 
+  bool isFromBackScattering() const { return (trackInfo_ >> 0) & 1; }
+  void setFromBackScattering() { trackInfo_ |= 1 << 0; }
+
+  bool isPrimary() const { return (trackInfo_ >> 1) & 1; }
+  void setIsPrimary() { trackInfo_ |= 1 << 1; }
+  void setGenParticleID(const int idx) { igenpart = idx; }
+  int getPrimaryID() const { return igenpart; }
+  uint8_t getTrackInfo() const { return trackInfo_; }
+
 private:
   int ivert;
   int igenpart;
 
   math::XYZVectorD tkposition;
   math::XYZTLorentzVectorD tkmomentum;
 
-  bool crossedBoundary_;
   int idAtBoundary_;
   math::XYZTLorentzVectorF positionAtBoundary_;
   math::XYZTLorentzVectorF momentumAtBoundary_;
+  uint8_t trackInfo_;
+  // explanation of trackInfo bits:
+  // 00000001 = simTrack is from backscattering
+  // 00000010 = simTrack is of a primary particle
+  // 00000100 = simTrack crossed the boundary
 };
 
 #include <iosfwd>

SimDataFormats/Track/src/SimTrack.cc

@@ -1,22 +1,22 @@
 #include "SimDataFormats/Track/interface/SimTrack.h"
 
-SimTrack::SimTrack() : ivert(-1), igenpart(-1), crossedBoundary_(false) {}
+SimTrack::SimTrack() : ivert(-1), igenpart(-1), trackInfo_(0) {}
 
 SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p)
-    : Core(ipart, p), ivert(-1), igenpart(-1), crossedBoundary_(false) {}
+    : Core(ipart, p), ivert(-1), igenpart(-1), trackInfo_(0) {}
 
 SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p, int iv, int ig)
-    : Core(ipart, p), ivert(iv), igenpart(ig), crossedBoundary_(false) {}
+    : Core(ipart, p), ivert(iv), igenpart(ig), trackInfo_(0) {}
 
 SimTrack::SimTrack(int ipart,
                    const math::XYZTLorentzVectorD& p,
                    int iv,
                    int ig,
                    const math::XYZVectorD& tkp,
                    const math::XYZTLorentzVectorD& tkm)
-    : Core(ipart, p), ivert(iv), igenpart(ig), tkposition(tkp), tkmomentum(tkm), crossedBoundary_(false) {}
+    : Core(ipart, p), ivert(iv), igenpart(ig), tkposition(tkp), tkmomentum(tkm), trackInfo_(0) {}
 
-SimTrack::SimTrack(const CoreSimTrack& t, int iv, int ig) : Core(t), ivert(iv), igenpart(ig), crossedBoundary_(false) {}
+SimTrack::SimTrack(const CoreSimTrack& t, int iv, int ig) : Core(t), ivert(iv), igenpart(ig), trackInfo_(0) {}
 
 std::ostream& operator<<(std::ostream& o, const SimTrack& t) {
   return o << (SimTrack::Core)(t) << ", " << t.vertIndex() << ", " << t.genpartIndex();

SimDataFormats/Track/src/classes_def.xml

@@ -4,7 +4,8 @@
   <class name="CoreSimTrack" ClassVersion="10">
    <version ClassVersion="10" checksum="3936841839"/>
   </class>
-  <class name="SimTrack" ClassVersion="12">
+  <class name="SimTrack" ClassVersion="13">
+   <version ClassVersion="13" checksum="1912247222"/>
    <version ClassVersion="12" checksum="3470347245"/>
    <version ClassVersion="11" checksum="1785575744"/>
    <version ClassVersion="10" checksum="1430205451"/>

SimG4Core/Notification/interface/TmpSimTrack.h

@@ -40,7 +40,7 @@ class TmpSimTrack {
   const math::XYZVectorD& momentum() const { return ip_; }
   double energy() const { return ie_; }
   int ivert() const { return ivert_; }
-  int igenpart() const { return igenpart_; }
+  int igenpart() const { return isPrimary_ ? igenpart_ : -1; }
   // parent momentum at interaction
   const math::XYZVectorD& parentMomentum() const { return parentMomentum_; }
   // Information at level of tracker surface
@@ -62,6 +62,12 @@ class TmpSimTrack {
   const math::XYZTLorentzVectorF& getPositionAtBoundary() const { return positionAtBoundary_; }
   const math::XYZTLorentzVectorF& getMomentumAtBoundary() const { return momentumAtBoundary_; }
   int getIDAtBoundary() const { return idAtBoundary_; }
+  bool isFromBackScattering() const { return isFromBackScattering_; }
+  void setFromBackScattering() { isFromBackScattering_ = true; }
+  void setGenParticleID(int i) { igenpart_ = i; }
+  bool isPrimary() const { return isPrimary_; }
+  void setIsPrimary() { isPrimary_ = true; }
+  int getPrimaryID() const { return igenpart_; }
 
 private:
   int id_;
@@ -74,7 +80,9 @@ class TmpSimTrack {
   math::XYZVectorD parentMomentum_{math::XYZVectorD(0., 0., 0.)};
   math::XYZVectorD tkSurfacePosition_{math::XYZVectorD(0., 0., 0.)};
   math::XYZTLorentzVectorD tkSurfaceMomentum_{math::XYZTLorentzVectorD(0., 0., 0., 0.)};
+  bool isFromBackScattering_{false};
   bool crossedBoundary_{false};
+  bool isPrimary_{false};
   int idAtBoundary_{-1};
   math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)};
   math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)};

SimG4Core/Notification/interface/TrackWithHistory.h

@@ -14,7 +14,7 @@ class G4PrimaryParticle;
 
 class TrackWithHistory {
 public:
-  /** The constructor is called at time, 
+  /** The constructor is called at time,
      *  when some of the information may not available yet.
      */
   TrackWithHistory(const G4Track *g4track, int pID);
@@ -27,7 +27,7 @@ class TrackWithHistory {
   int trackID() const { return trackID_; }
   int particleID() const { return pdgID_; }
   int parentID() const { return parentID_; }
-  int genParticleID() const { return genParticleID_; }
+  int genParticleID() const { return isPrimary_ ? genParticleID_ : -1; }
   int vertexID() const { return vertexID_; }
   int processType() const { return procType_; }
   int getIDAtBoundary() const { return idAtBoundary_; }
@@ -67,6 +67,11 @@ class TrackWithHistory {
     tkSurfacePosition_ = pos;
     tkSurfaceMomentum_ = mom;
   }
+  bool isFromBackScattering() const { return isFromBackScattering_; }
+  void setFromBackScattering() { isFromBackScattering_ = true; }
+  bool isPrimary() const { return isPrimary_; }
+  void setIsPrimary() { isPrimary_ = true; }
+  int getPrimaryID() const { return genParticleID_; }
 
 private:
   int trackID_;
@@ -88,6 +93,8 @@ class TrackWithHistory {
   bool storeTrack_{false};
   bool saved_{false};
   bool crossedBoundary_{false};
+  bool isFromBackScattering_{false};
+  bool isPrimary_{false};
 };
 
 extern G4ThreadLocal G4Allocator<TrackWithHistory> *fpTrackWithHistoryAllocator;

SimG4Core/Notification/src/SimTrackManager.cc

@@ -73,6 +73,13 @@ void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, b
   std::pair<int, int> thePair(iTrack->trackID(), iTrack->parentID());
   idsave.push_back(thePair);
   if (inHistory) {
+    auto info = static_cast<const TrackInformation *>(track->GetUserInformation());
+    if (info->isInTrkFromBackscattering())
+      iTrack->setFromBackScattering();
+    // set there *for all the tracks* the genParticle ID associated with the G4Track
+    // in the constructor of TrackWithHistory the info isPrimary is saved and used
+    // to give -1 if the track is not a primary
+    iTrack->setGenParticleID(info->mcTruthID());
     m_trackContainer.push_back(iTrack);
     const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition();
     std::pair<int, math::XYZVectorD> p(iTrack->trackID(),
@@ -132,7 +139,10 @@ void SimTrackManager::reallyStoreTracks() {
     // at this stage there is one vertex per track,
     // so the vertex id of track N is also N
     int iParentID = trkH->parentID();
-    int ig = trkH->genParticleID();
+    int ig = trkH->genParticleID(); // filled only for primary tracks
+    bool isBackScatter = trkH->isFromBackScattering();
+    bool isPrimary = trkH->isPrimary();
+    int primaryGenPartId = trkH->getPrimaryID(); // filled if the G4Track had this info
     int ivertex = getOrCreateVertex(trkH, iParentID);
 
     auto ptr = trkH;
@@ -170,6 +180,11 @@ void SimTrackManager::reallyStoreTracks() {
     TmpSimTrack* g4simtrack =
         new TmpSimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom);
     g4simtrack->copyCrossedBoundaryVars(trkH);
+    if (isBackScatter)
+      g4simtrack->setFromBackScattering();
+    if (isPrimary)
+      g4simtrack->setIsPrimary();
+    g4simtrack->setGenParticleID(primaryGenPartId);
     m_simEvent->addTrack(g4simtrack);
   }
 }

SimG4Core/Notification/src/TmpSimEvent.cc

@@ -35,15 +35,25 @@ void TmpSimEvent::load(edm::SimTrackContainer& c) const {
     int iv = trk->ivert();
     int ig = trk->igenpart();
     int id = trk->id();
+    bool isBackScatter = trk->isFromBackScattering();
+    bool isPrimary = trk->isPrimary();
+    int primaryGenPartId = trk->getPrimaryID(); // filled if the G4Track had this info
     // ip = particle ID as PDG
-    // pp = 4-momentum in GeV
+    // p  = 4-momentum in GeV
     // iv = corresponding TmpSimVertex index
-    // ig = corresponding GenParticle index
+    // ig = corresponding GenParticle index only if is a primary, otherwise is -1
+    // primaryGenPartId = corresponding GenParticle index also if not primary
+    // id = corresponding g4Track Id
     SimTrack t = SimTrack(ip, p, iv, ig, trk->trackerSurfacePosition(), trk->trackerSurfaceMomentum());
     t.setTrackId(id);
     t.setEventId(EncodedEventId(0));
     t.setCrossedBoundaryVars(
         trk->crossedBoundary(), trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary());
+    if (isBackScatter)
+      t.setFromBackScattering();
+    if (isPrimary)
+      t.setIsPrimary();
+    t.setGenParticleID(primaryGenPartId);
     c.push_back(t);
   }
   std::stable_sort(c.begin(), c.end(), IdSort());

SimG4Core/Notification/src/TrackWithHistory.cc

@@ -29,10 +29,15 @@ TrackWithHistory::TrackWithHistory(const G4Track* g4trk, int pID) {
   TrackInformation* trkinfo = static_cast<TrackInformation*>(g4trk->GetUserInformation());
   storeTrack_ = trkinfo->storeTrack();
   auto vgprimary = g4trk->GetDynamicParticle()->GetPrimaryParticle();
+  // GetPrimaryParticle() returns the pointer to the corresponding G4PrimaryParticle object
+  // if this particle is a primary particle OR is defined as a
+  // pre-assigned decay product. Otherwise return nullptr.
+  // Therefore here the genParticleID_ is -1 for not primary particles
   if (vgprimary != nullptr) {
     auto priminfo = static_cast<GenParticleInfo*>(vgprimary->GetUserInformation());
     if (nullptr != priminfo) {
       genParticleID_ = priminfo->id();
+      isPrimary_ = true;
     }
   }
   // V.I. weight is computed in the same way as before
@@ -53,6 +58,7 @@ TrackWithHistory::TrackWithHistory(const G4PrimaryParticle* ptr, int trackID, co
   auto priminfo = static_cast<GenParticleInfo*>(ptr->GetUserInformation());
   if (nullptr != priminfo) {
     genParticleID_ = priminfo->id();
+    isPrimary_ = true;
   }
   weight_ = 10000. * genParticleID_;
 }