Add variable to track if electron originates from primary electron (#22)

Author: tobias-boeckh

Pinpoint/include/AnalysisManager.hh

@@ -189,6 +189,7 @@ class AnalysisManager {
     // std::vector<G4bool> fPixelFromPrimaryPizero;
     // std::vector<G4bool> fPixelFromFSLPizero;
     std::vector<G4bool> fPixelFromPrimaryLepton;
+    std::vector<G4bool> fPixelFromPrimaryEMShower;
 
     // Truth position of hit in x, y, z
     std::vector<Float_t> fPixelTruthX;

Pinpoint/include/PixelHit.hh

@@ -12,7 +12,7 @@ class PixelHit : public G4VHit
 {
 public:
   PixelHit() = default;
-  PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim);
+  PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower);
   PixelHit(const PixelHit&) = default;
   ~PixelHit() override = default;
 
@@ -44,6 +44,7 @@ public:
   // void SetFromPrimaryPizero(G4bool fromPrimaryPizero) { fFromPrimaryPizero = fromPrimaryPizero; }
   // void SetFromFSLPizero(G4bool fromFSLPizero) { fFromFSLPizero = fromFSLPizero; }
   void SetFromPrimaryLepton(G4bool fromPrimaryLepton) { fFromPrimaryLepton = fromPrimaryLepton; }
+  void SetFromPrimaryEMShower(G4bool fromPrimaryEMShower) { fFromPrimaryEMShower = fromPrimaryEMShower; }
   // void SetTruthHitPos(G4ThreeVector pos) { fTruthHitPos = pos; }
 
   G4int GetPDGCode() const { return fPDGCode; }
@@ -65,6 +66,7 @@ public:
   // G4bool GetFromPrimaryPizero() const { return fFromPrimaryPizero; }
   // G4bool GetFromFSLPizero() const { return fFromFSLPizero; }
   G4bool GetFromPrimaryLepton() const { return fFromPrimaryLepton; }
+  G4bool GetFromPrimaryEMShower() const { return fFromPrimaryEMShower; }
 
 private:
   G4int fTrackID = -1;
@@ -76,6 +78,7 @@ private:
   G4int fLayerID = -1;
   // G4int fCharge = 0;
   G4bool fFromPrimaryLepton = false;
+  G4bool fFromPrimaryEMShower = false;
   // G4bool fFromPrimaryPizero = false;
   // G4bool fFromFSLPizero = false;
   // G4ThreeVector fTruthHitPos;

Pinpoint/include/TrackInformation.hh

@@ -24,9 +24,11 @@ class TrackInformation : public G4VUserTrackInformation
     inline void SetTrackIsFromPrimaryPizero(G4int i) {fFromPrimaryPizero = i;}
     inline void SetTrackIsFromFSLPizero(G4int i) {fFromFSLPizero = i;}
     inline void SetTrackIsFromPrimaryLepton(G4int i) {fFromPrimaryLepton = i;}
+    inline void SetTrackIsFromPrimaryEMShower(G4int i) {fFromPrimaryEMShower = i;}
     inline G4int IsTrackFromPrimaryPizero() const {return fFromPrimaryPizero;}
     inline G4int IsTrackFromFSLPizero() const {return fFromFSLPizero;}
     inline G4int IsTrackFromPrimaryLepton() const {return fFromPrimaryLepton;}
+    inline G4int IsTrackFromPrimaryEMShower() const {return fFromPrimaryEMShower;}
 
     inline void InsertHit(G4long hitIndex) { 
       if (!fHitIndices) {
@@ -40,7 +42,7 @@ class TrackInformation : public G4VUserTrackInformation
     G4int fFromPrimaryPizero;
     G4int fFromFSLPizero;
     G4int fFromPrimaryLepton;
-
+    G4int fFromPrimaryEMShower;
     std::shared_ptr<std::vector<G4long>> fHitIndices;
 
     // std::vector<G4long>* fHitIndices; //TODO: Make this a smart pointer?

Pinpoint/src/AnalysisManager.cc

@@ -198,6 +198,7 @@ void AnalysisManager::bookHitsTrees()
   // fPixelHitsTree->Branch("hit_fromPrimaryPizero", &fPixelFromPrimaryPizero);
   // fPixelHitsTree->Branch("hit_fromFSLPizero", &fPixelFromFSLPizero);
   fPixelHitsTree->Branch("hit_fromPrimaryLepton", &fPixelFromPrimaryLepton);
+  fPixelHitsTree->Branch("hit_fromPrimaryEMShower", &fPixelFromPrimaryEMShower);
 
   // if (fSaveTruthHits)
   // {
@@ -310,6 +311,7 @@ void AnalysisManager::BeginOfEvent()
   // fPixelFromPrimaryPizero.clear();
   // fPixelFromFSLPizero.clear();
   fPixelFromPrimaryLepton.clear();
+  fPixelFromPrimaryEMShower.clear();
   // fPixelTruthX.clear();
   // fPixelTruthY.clear();
   // fPixelTruthZ.clear();
@@ -581,6 +583,7 @@ void AnalysisManager::FillHitsOutput()
           // fPixelFromPrimaryPizero.push_back(hit->GetFromPrimaryPizero());
           // fPixelFromFSLPizero.push_back(hit->GetFromFSLPizero());
           fPixelFromPrimaryLepton.push_back(hit->GetFromPrimaryLepton());
+          fPixelFromPrimaryEMShower.push_back(hit->GetFromPrimaryEMShower());
 
           // if (fSaveTruthHits)
           // {

Pinpoint/src/PixelHit.cc

@@ -41,7 +41,7 @@ void PixelHit::Print()
      << G4endl;
 }
 
-PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim)
-: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}
+PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower)
+: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}, fFromPrimaryEMShower{isEMShower}
 {
 }

Pinpoint/src/PixelSD.cc

@@ -86,6 +86,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step)
   const auto* info =static_cast<const TrackInformation*>(track->GetUserInformation());
   G4bool fromPrimaryLepton = info && info->IsTrackFromPrimaryLepton() || parentID ==0;
   fromPrimaryLepton = fromPrimaryLepton && (std::abs(pdgid) == 11 || std::abs(pdgid) == 13 || std::abs(pdgid) == 15);
+  G4bool fromPrimaryEMShower = info && info->IsTrackFromPrimaryEMShower();
 
   assert(rowID < fNPixelsY);
   assert(colID < fNPixelsX);
@@ -106,7 +107,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step)
   } else {
     fPixelHits.push_back(
       new PixelHit(edep, rowID, colID, layerID,
-                  trackID, parentID, pdgid, fromPrimaryLepton)
+                  trackID, parentID, pdgid, fromPrimaryLepton, fromPrimaryEMShower)
     );
   }
 

Pinpoint/src/TrackInformation.cc

@@ -9,6 +9,7 @@ TrackInformation::TrackInformation()
   fFromPrimaryPizero = 0;
   fFromFSLPizero = 0;
   fFromPrimaryLepton = 0;
+  fFromPrimaryEMShower = 0;
   fHitIndices = std::make_shared<std::vector<G4long>>();
 }
 
@@ -18,6 +19,7 @@ TrackInformation::TrackInformation(const G4Track* aTrack)
   fFromPrimaryPizero = 0;
   fFromFSLPizero = 0;
   fFromPrimaryLepton = 0;
+  fFromPrimaryEMShower = 0;
   fHitIndices = std::make_shared<std::vector<G4long>>();
 }
 
@@ -30,7 +32,7 @@ ::operator =(const TrackInformation& aTrackInfo)
   fFromPrimaryPizero = aTrackInfo.fFromPrimaryPizero;
   fFromFSLPizero = aTrackInfo.fFromFSLPizero;
   fFromPrimaryLepton = aTrackInfo.fFromPrimaryLepton;
-
+  fFromPrimaryEMShower = aTrackInfo.fFromPrimaryEMShower;
   return *this;
 }
 
@@ -39,5 +41,6 @@ void TrackInformation::Print() const
     G4cout << "Is from primary pizero " << fFromPrimaryPizero << G4endl;
     G4cout << "Is from final state lepton decay pizero " << fFromFSLPizero << G4endl;
     G4cout << "Is from primary lepton (tau or muon)    " << fFromPrimaryLepton << G4endl;
+    G4cout << "Is from primary EM shower    " << fFromPrimaryEMShower << G4endl;
 }
 

Pinpoint/src/TrackingAction.cc

@@ -71,6 +71,42 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
     }
   }
 
+  if (aTrack->GetTrackID()==1 && abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==11)
+  {
+    const auto* info =static_cast<const TrackInformation*>(aTrack->GetUserInformation());
+    TrackInformation* newInfo =  new TrackInformation();
+    if (info) {
+      newInfo->SetTrackIsFromPrimaryPizero(info->IsTrackFromPrimaryPizero());
+      newInfo->SetTrackIsFromFSLPizero(info->IsTrackFromFSLPizero());
+      newInfo->SetTrackIsFromPrimaryLepton(info->IsTrackFromPrimaryLepton());
+    }
+    newInfo->SetTrackIsFromPrimaryEMShower(1);
+    aTrack->SetUserInformation(newInfo);
+  }
+
+  const auto* info =static_cast<const TrackInformation*>(aTrack->GetUserInformation());
+  G4bool fromPrimaryEMShower = info && info->IsTrackFromPrimaryEMShower();
+
+  if (fromPrimaryEMShower &&
+      (abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==11 ||
+       abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==22))
+  {
+    G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
+    if (secondaries)
+    {
+      size_t nSeco = secondaries->size();
+      if (nSeco>0)
+      {
+        for (size_t i=0; i<nSeco; ++i)
+        {
+          TrackInformation* info =  new TrackInformation();
+          info->SetTrackIsFromPrimaryEMShower(1);
+          (*secondaries)[i]->SetUserInformation(info);
+        }
+      }
+    }
+  }
+
   if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") 
   {
     // in case of tau decay pizero