diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index b054bed..2325215 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -189,6 +189,7 @@ class AnalysisManager { // std::vector fPixelFromPrimaryPizero; // std::vector fPixelFromFSLPizero; std::vector fPixelFromPrimaryLepton; + std::vector fPixelFromPrimaryEMShower; // Truth position of hit in x, y, z std::vector fPixelTruthX; diff --git a/Pinpoint/include/PixelHit.hh b/Pinpoint/include/PixelHit.hh index 590e250..6f84fe5 100644 --- a/Pinpoint/include/PixelHit.hh +++ b/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; diff --git a/Pinpoint/include/TrackInformation.hh b/Pinpoint/include/TrackInformation.hh index 59c4ae5..1a4abbc 100644 --- a/Pinpoint/include/TrackInformation.hh +++ b/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> fHitIndices; // std::vector* fHitIndices; //TODO: Make this a smart pointer? diff --git a/Pinpoint/src/AnalysisManager.cc b/Pinpoint/src/AnalysisManager.cc index b564339..525379a 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/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) // { diff --git a/Pinpoint/src/PixelHit.cc b/Pinpoint/src/PixelHit.cc index 4d3f29e..b5d2735 100644 --- a/Pinpoint/src/PixelHit.cc +++ b/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} { } \ No newline at end of file diff --git a/Pinpoint/src/PixelSD.cc b/Pinpoint/src/PixelSD.cc index 6d0c536..c1013ab 100644 --- a/Pinpoint/src/PixelSD.cc +++ b/Pinpoint/src/PixelSD.cc @@ -86,6 +86,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) const auto* info =static_cast(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) ); } diff --git a/Pinpoint/src/TrackInformation.cc b/Pinpoint/src/TrackInformation.cc index 8f4f8a2..2e2df1f 100644 --- a/Pinpoint/src/TrackInformation.cc +++ b/Pinpoint/src/TrackInformation.cc @@ -9,6 +9,7 @@ TrackInformation::TrackInformation() fFromPrimaryPizero = 0; fFromFSLPizero = 0; fFromPrimaryLepton = 0; + fFromPrimaryEMShower = 0; fHitIndices = std::make_shared>(); } @@ -18,6 +19,7 @@ TrackInformation::TrackInformation(const G4Track* aTrack) fFromPrimaryPizero = 0; fFromFSLPizero = 0; fFromPrimaryLepton = 0; + fFromPrimaryEMShower = 0; fHitIndices = std::make_shared>(); } @@ -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; } diff --git a/Pinpoint/src/TrackingAction.cc b/Pinpoint/src/TrackingAction.cc index 900d7fa..f27147e 100644 --- a/Pinpoint/src/TrackingAction.cc +++ b/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(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(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; iSetTrackIsFromPrimaryEMShower(1); + (*secondaries)[i]->SetUserInformation(info); + } + } + } + } + if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") { // in case of tau decay pizero