diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index ac2a3af..8c7ba96 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -208,6 +208,7 @@ class AnalysisManager { std::vector fPixelFromPrimaryLepton; std::vector fPixelFromPrimaryEMShower; std::vector fPixelFromCharmedHadron; + std::vector fPixelFromTau; // Truth position of hit in x, y, z std::vector fPixelTruthX; diff --git a/Pinpoint/include/PixelHit.hh b/Pinpoint/include/PixelHit.hh index ca53b2f..28eca09 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, G4bool isEMShower, G4bool isFromCharmed); + PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower, G4bool isFromCharmed, G4bool isFromTau); PixelHit(const PixelHit&) = default; ~PixelHit() override = default; @@ -46,6 +46,7 @@ public: void SetFromPrimaryLepton(G4bool fromPrimaryLepton) { fFromPrimaryLepton = fromPrimaryLepton; } void SetFromPrimaryEMShower(G4bool fromPrimaryEMShower) { fFromPrimaryEMShower = fromPrimaryEMShower; } void SetFromCharmedHadron(G4bool fromCharmedHadron) { fFromCharmedHadron = fromCharmedHadron; } + void SetFromTau(G4bool fromTau) { fFromTau = fromTau; } // void SetTruthHitPos(G4ThreeVector pos) { fTruthHitPos = pos; } G4int GetPDGCode() const { return fPDGCode; } @@ -69,6 +70,7 @@ public: G4bool GetFromPrimaryLepton() const { return fFromPrimaryLepton; } G4bool GetFromPrimaryEMShower() const { return fFromPrimaryEMShower; } G4bool GetFromCharmedHadron() const { return fFromCharmedHadron; } + G4bool GetFromTau() const { return fFromTau; } private: G4int fTrackID = -1; @@ -82,6 +84,7 @@ private: G4bool fFromPrimaryLepton = false; G4bool fFromPrimaryEMShower = false; G4bool fFromCharmedHadron = false; + G4bool fFromTau = false; // G4bool fFromPrimaryPizero = false; // G4bool fFromFSLPizero = false; // G4ThreeVector fTruthHitPos; diff --git a/Pinpoint/include/TrackInformation.hh b/Pinpoint/include/TrackInformation.hh index d4db35d..57b04e3 100644 --- a/Pinpoint/include/TrackInformation.hh +++ b/Pinpoint/include/TrackInformation.hh @@ -26,11 +26,13 @@ class TrackInformation : public G4VUserTrackInformation inline void SetTrackIsFromPrimaryLepton(G4int i) {fFromPrimaryLepton = i;} inline void SetTrackIsFromPrimaryEMShower(G4int i) {fFromPrimaryEMShower = i;} inline void SetTrackIsFromCharmedHadron(G4int i) {fFromCharmedHadron = i;} + inline void SetTrackIsFromTau(G4int i) {fFromTau = 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 G4int IsTrackFromCharmedHadron() const {return fFromCharmedHadron;} + inline G4int IsTrackFromTau() const {return fFromTau;} inline void InsertHit(G4long hitIndex) { if (!fHitIndices) { @@ -46,6 +48,7 @@ class TrackInformation : public G4VUserTrackInformation G4int fFromPrimaryLepton; G4int fFromPrimaryEMShower; G4int fFromCharmedHadron; + G4int fFromTau; 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 d6fc893..36837e5 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/Pinpoint/src/AnalysisManager.cc @@ -206,6 +206,7 @@ void AnalysisManager::bookHitsTrees() fPixelHitsTree->Branch("hit_fromPrimaryLepton", &fPixelFromPrimaryLepton); fPixelHitsTree->Branch("hit_fromPrimaryEMShower", &fPixelFromPrimaryEMShower); fPixelHitsTree->Branch("hit_fromCharmedHadron", &fPixelFromCharmedHadron); + fPixelHitsTree->Branch("hit_fromTau", &fPixelFromTau); // if (fSaveTruthHits) // { @@ -387,6 +388,7 @@ void AnalysisManager::BeginOfEvent() fPixelFromPrimaryLepton.clear(); fPixelFromPrimaryEMShower.clear(); fPixelFromCharmedHadron.clear(); + fPixelFromTau.clear(); // fPixelTruthX.clear(); // fPixelTruthY.clear(); // fPixelTruthZ.clear(); @@ -684,6 +686,7 @@ void AnalysisManager::FillHitsOutput() fPixelFromPrimaryLepton.push_back(hit->GetFromPrimaryLepton()); fPixelFromPrimaryEMShower.push_back(hit->GetFromPrimaryEMShower()); fPixelFromCharmedHadron.push_back(hit->GetFromCharmedHadron()); + fPixelFromTau.push_back(hit->GetFromTau()); // if (fSaveTruthHits) // { diff --git a/Pinpoint/src/PixelHit.cc b/Pinpoint/src/PixelHit.cc index 2120cab..c8f5c3e 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, G4bool isEMShower, G4bool isFromCharmed) -: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}, fFromPrimaryEMShower{isEMShower}, fFromCharmedHadron{isFromCharmed} +PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower, G4bool isFromCharmed, G4bool isFromTau) +: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}, fFromPrimaryEMShower{isEMShower}, fFromCharmedHadron{isFromCharmed}, fFromTau{isFromTau} { } \ No newline at end of file diff --git a/Pinpoint/src/PixelSD.cc b/Pinpoint/src/PixelSD.cc index 2b21a14..9a4ea38 100644 --- a/Pinpoint/src/PixelSD.cc +++ b/Pinpoint/src/PixelSD.cc @@ -88,6 +88,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) fromPrimaryLepton = fromPrimaryLepton && (std::abs(pdgid) == 11 || std::abs(pdgid) == 13 || std::abs(pdgid) == 15); G4bool fromPrimaryEMShower = info && info->IsTrackFromPrimaryEMShower(); G4bool fromCharmedHadron = info && info->IsTrackFromCharmedHadron(); + G4bool fromTau = info && info->IsTrackFromTau(); assert(rowID < fNPixelsY); assert(colID < fNPixelsX); @@ -108,7 +109,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) } else { fPixelHits.push_back( new PixelHit(edep, rowID, colID, layerID, - trackID, parentID, pdgid, fromPrimaryLepton, fromPrimaryEMShower, fromCharmedHadron) + trackID, parentID, pdgid, fromPrimaryLepton, fromPrimaryEMShower, fromCharmedHadron, fromTau) ); } diff --git a/Pinpoint/src/TrackInformation.cc b/Pinpoint/src/TrackInformation.cc index d526025..134fed5 100644 --- a/Pinpoint/src/TrackInformation.cc +++ b/Pinpoint/src/TrackInformation.cc @@ -11,6 +11,7 @@ TrackInformation::TrackInformation() fFromPrimaryLepton = 0; fFromPrimaryEMShower = 0; fFromCharmedHadron = 0; + fFromTau = 0; fHitIndices = std::make_shared>(); } @@ -21,7 +22,7 @@ TrackInformation::TrackInformation(const G4Track* aTrack) fFromFSLPizero = 0; fFromPrimaryLepton = 0; fFromPrimaryEMShower = 0; - fFromCharmedHadron = 0; + fFromTau = 0; fHitIndices = std::make_shared>(); } @@ -35,7 +36,7 @@ ::operator =(const TrackInformation& aTrackInfo) fFromFSLPizero = aTrackInfo.fFromFSLPizero; fFromPrimaryLepton = aTrackInfo.fFromPrimaryLepton; fFromPrimaryEMShower = aTrackInfo.fFromPrimaryEMShower; - fFromCharmedHadron = aTrackInfo.fFromCharmedHadron; + fFromTau = aTrackInfo.fFromTau; return *this; } @@ -46,5 +47,6 @@ void TrackInformation::Print() const G4cout << "Is from primary lepton (tau or muon) " << fFromPrimaryLepton << G4endl; G4cout << "Is from primary EM shower " << fFromPrimaryEMShower << G4endl; G4cout << "Is from charmed hadron " << fFromCharmedHadron << G4endl; + G4cout << "Is from tau " << fFromTau << G4endl; } diff --git a/Pinpoint/src/TrackingAction.cc b/Pinpoint/src/TrackingAction.cc index d04c195..f4a79b1 100644 --- a/Pinpoint/src/TrackingAction.cc +++ b/Pinpoint/src/TrackingAction.cc @@ -157,6 +157,50 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) } } + absPDG = std::abs(aTrack->GetParticleDefinition()->GetPDGEncoding()); + G4bool isTau = absPDG == 15; + G4bool fromTau = info && info->IsTrackFromTau(); + // If this is a primary tau lepton, its decay products are counted as primary particles + if ((isTau) && (aTrack->GetParentID() == 0 || aTrack->GetParentID() == 1)) + { + G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); + if (secondaries) + { + size_t nSeco = secondaries->size(); + if (nSeco>0) + { + for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") + { + TrackInformation* newInfo = new TrackInformation(); + newInfo->SetTrackIsFromTau(1); + (*secondaries)[i]->SetUserInformation(newInfo); + AnalysisManager::GetInstance()->AddOnePrimaryTrack(); + } + } + } + } + } + // If this track comes from a charmed hadron, tag all its decay products + else if (isTau || fromTau) + { + G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); + if (secondaries) + { + size_t nSeco = secondaries->size(); + if (nSeco>0) + { + for (size_t i=0; iSetTrackIsFromTau(1); + (*secondaries)[i]->SetUserInformation(newInfo); + } + } + } + } + if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") { // in case of tau decay pizero