Add additional information to the SimTracks

SimDataFormats/Track/interface/SimTrack.h

@@ -34,8 +34,9 @@ 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; }
+  bool isPrimary() const { return (trackInfo_ >> 1) & 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 +52,40 @@ 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_ & 1; }
+  void setFromBackScattering() { trackInfo_ |= 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,10 +4,27 @@
   <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"/>
+   <ioread sourceClass = "SimTrack" version="[-12]" targetClass="SimTrack" source="" target="">
+    <![CDATA[
+        SimTrack tmp(newObj->type(), newObj->momentum(), newObj->vertIndex(), newObj->genpartIndex(), newObj->trackerSurfacePosition(), newObj->trackerSurfaceMomentum());
+        tmp.setTrackId(newObj->trackId());
+        tmp.setEventId(newObj->eventId());
+        tmp.setCrossedBoundaryVars(
+            newObj->crossedBoundary(), newObj->getIDAtBoundary(), newObj->getPositionAtBoundary(), newObj->getMomentumAtBoundary());
+        if (newObj->isFromBackScattering()) {
+          tmp.setFromBackScattering();
+        }
+        if (newObj->genpartIndex() != -1) {
+          tmp.setIsPrimary();
+        }
+        *newObj=tmp;
+    ]]>
+   </ioread>
   </class>
   <class name="edm::Wrapper<std::vector<SimTrack> >" />
   <class name="edm::RefProd<std::vector<SimTrack> >"/>

SimG4Core/Application/plugins/OscarMTProducer.cc

@@ -251,7 +251,7 @@ void OscarMTProducer::produce(edm::Event& e, const edm::EventSetup& es) {
 
   if (0 < m_verbose) {
     edm::LogVerbatim("SimG4CoreApplication")
-        << "Produced " << p2->size() << " SimVertecies: position(cm), time(s), parentID, vertexID, processType";
+        << "Produced " << p2->size() << " SimVertices: position(cm), time(s), parentID, vertexID, processType";
     if (1 < m_verbose) {
       int nn = p2->size();
       for (int i = 0; i < nn; ++i) {
@@ -260,12 +260,15 @@ void OscarMTProducer::produce(edm::Event& e, const edm::EventSetup& es) {
     }
     edm::LogVerbatim("SimG4CoreApplication")
         << "Produced " << p1->size()
-        << " SimTracks: pdg, 4-momentum(GeV), vertexID, mcTruthID, flagBoundary, trackID at boundary";
+        << " SimTracks: G4 Id, pdg, 4-momentum(GeV), vertexID, mcTruthID, crossedBoundary -> trackID at boundary, from "
+           "backscattering, isPrimary -> getPrimary";
     if (1 < m_verbose) {
       int nn = p1->size();
       for (int i = 0; i < nn; ++i) {
-        edm::LogVerbatim("Track") << " " << i << ". " << (*p1)[i] << " " << (*p1)[i].crossedBoundary() << " "
-                                  << (*p1)[i].getIDAtBoundary();
+        edm::LogVerbatim("Track") << " " << i << ". " << (*p1)[i].trackId() << ", " << (*p1)[i] << ", "
+                                  << (*p1)[i].crossedBoundary() << "-> " << (*p1)[i].getIDAtBoundary() << ", "
+                                  << (*p1)[i].isFromBackScattering() << ", " << (*p1)[i].isPrimary() << "-> "
+                                  << (*p1)[i].getPrimaryID();
       }
     }
   }

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,15 @@ 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 the *non-primary* tracks the genParticle ID associated with the G4Track
+    // for the primaries this is done in the TrackWithHistory constructor.
+    // In the constructor of TrackWithHistory the info isPrimary is saved and used
+    // to give -1 if the track is not a primary.
+    if (not iTrack->isPrimary())
+      iTrack->setGenParticleID(info->mcTruthID());
     m_trackContainer.push_back(iTrack);
     const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition();
     std::pair<int, math::XYZVectorD> p(iTrack->trackID(),
@@ -124,15 +133,18 @@ void SimTrackManager::storeTracks() {
 void SimTrackManager::reallyStoreTracks() {
   // loop over the (now ordered) vector and really save the tracks
 #ifdef DebugLog
-  edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size();
+  edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer.size();
 #endif
 
   int nn = m_endPoints.size();
   for (auto& trkH : m_trackContainer) {
     // 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 +182,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);
   }
 }
@@ -304,7 +321,7 @@ void SimTrackManager::cleanTracksWithHistory() {
   std::stable_sort(idsave.begin(), idsave.end());
 
 #ifdef DebugLog
-  LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trksForThisEvent->size()
+  LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trackContainer.size()
                               << " tracks with history before branching";
   for (unsigned int it = 0; it < m_trackContainer.size(); it++) {
     LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number "

SimG4Core/Notification/src/TmpSimEvent.cc

@@ -35,15 +35,26 @@ 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();
+    else
+      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_;
 }