Skip to content

Add charm tau trees #26

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
benw22022 merged 6 commits into from
Feb 20, 2026
Merged

Add charm tau trees #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d2f4ecb
Add variable to track if electron originates from primary electron
tobias-boeckh Jan 29, 2026
1f16d17
Merge branch 'master' of github.com:PinpointDetector/Pinpoint_G4
tobias-boeckh Feb 6, 2026
0c3898d
Fix G4Exceptions: remove slow particles in field region which may curl
tobias-boeckh Feb 6, 2026
289836a
Fix vertex z position is fUseFixedZPosition is false
tobias-boeckh Feb 6, 2026
2f75a75
Add trees with decay products of primary tau and charm
tobias-boeckh Feb 20, 2026
2074ad5
Fix pythia event id
tobias-boeckh Feb 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions Pinpoint/include/AnalysisManager.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
#include "AnalysisManagerMessenger.hh"
#include "FPFParticle.hh"

class G4Track;

class AnalysisManager {
public:

Expand Down Expand Up @@ -54,6 +56,8 @@ class AnalysisManager {
void bookHitsTrees();
void bookScintTrees();
void bookFaserTree();
void bookTauTree();
void bookCharmTree();

void FillEventTree(const G4Event* event);
void FillPrimariesTree(const G4Event* event);
Expand All @@ -62,6 +66,7 @@ class AnalysisManager {
void FillHitsOutput();
void FillScintOutput();
void FillFaserOutput();
void FillPythiaTree(const G4Event* event, TTree* tree, std::vector<int> pdg_ids);

float_t GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m);

Expand All @@ -87,6 +92,8 @@ class AnalysisManager {
TTree* fTrk;
TTree* fPrim;
TTree* fGeom;
TTree* fTauTree;
TTree* fCharmTree;

TDirectory* fHits;
TTree* fPixelHitsTree;
Expand Down Expand Up @@ -134,6 +141,12 @@ class AnalysisManager {
double trackKinE;
int trackNPoints;
double trackTheta;
double trackProdX;
double trackProdY;
double trackProdZ;
double trackDecayX;
double trackDecayY;
double trackDecayZ;
// std::vector<double> trackPointX;
// std::vector<double> trackPointY;
// std::vector<double> trackPointZ;
Expand Down Expand Up @@ -228,6 +241,24 @@ class AnalysisManager {
std::vector<float> fFaserE; // Energy
std::vector<float> fFaserEdep; // Energy deposit
std::vector<float> fFaserCharge; // Particle charge

//
//----------------------------------------------------
// Output variables for tau and charm trees
UInt_t fPythiaEventID;
int fPythiaTPDG;
int fPythiaPPDG;
int fPythiaTID;
int fPythiaPID;
double fPythiaProdX;
double fPythiaProdY;
double fPythiaProdZ;
double fPythiaDecayX;
double fPythiaDecayY;
double fPythiaDecayZ;
double fPythiaPx;
double fPythiaPy;
double fPythiaPz;
};

#endif
171 changes: 151 additions & 20 deletions Pinpoint/src/AnalysisManager.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
#include <vector>
#include <functional>
#include <iostream>
#include <string>
Expand Down Expand Up @@ -67,6 +66,8 @@ AnalysisManager::AnalysisManager()
fTrk = nullptr;
fPrim = nullptr;
fPixelHitsTree = nullptr;
fTauTree = nullptr;
fCharmTree = nullptr;
// fActsParticlesTree = nullptr;

fSaveTrack = false;
Expand Down Expand Up @@ -148,9 +149,12 @@ void AnalysisManager::bookTrkTree()
fTrk->Branch("trackKinE", &trackKinE, "trackKinE/D");
fTrk->Branch("trackNPoints", &trackNPoints, "trackNPoints/I");
fTrk->Branch("trackTheta", &trackTheta, "trackTheta/D");
// fTrk->Branch("trackPointX", &trackPointX);
// fTrk->Branch("trackPointY", &trackPointY);
// fTrk->Branch("trackPointZ", &trackPointZ);
fTrk->Branch("trackProdX", &trackProdX, "trackProdX/D");
fTrk->Branch("trackProdY", &trackProdY, "trackProdY/D");
fTrk->Branch("trackProdZ", &trackProdZ, "trackProdZ/D");
fTrk->Branch("trackDecayX", &trackDecayX, "trackDecayX/D");
fTrk->Branch("trackDecayY", &trackDecayY, "trackDecayY/D");
fTrk->Branch("trackDecayZ", &trackDecayZ, "trackDecayZ/D");
}


Expand Down Expand Up @@ -231,7 +235,7 @@ void AnalysisManager::bookScintTrees()
fScintTree->Branch("fromPrimaryLepton", &fScintFromPrimaryLepton);
}


//// --- FASER tracking spectrometer hits ---
void AnalysisManager::bookFaserTree()
{
fFile->cd(fHits->GetName());
Expand All @@ -252,6 +256,46 @@ void AnalysisManager::bookFaserTree()
fFaserHitsTree->Branch("charge", &fFaserCharge);
}

//// --- Children of primary tau decay ---
void AnalysisManager::bookTauTree()
{
fTauTree = new TTree("tau", "Children of tau decay");
fTauTree->Branch("evtID", &fPythiaEventID, "evtID/i");
fTauTree->Branch("PID", &fPythiaPID, "PID/I");
fTauTree->Branch("TID", &fPythiaTID, "TID/I");
fTauTree->Branch("PDG", &fPythiaTPDG, "PDG/I");
fTauTree->Branch("ParentPDG", &fPythiaPPDG, "ParentPDG/I");
fTauTree->Branch("ProdX", &fPythiaProdX, "ProdX/D");
fTauTree->Branch("ProdY", &fPythiaProdY, "ProdY/D");
fTauTree->Branch("ProdZ", &fPythiaProdZ, "ProdZ/D");
fTauTree->Branch("DecayX", &fPythiaDecayX, "DecayX/D");
fTauTree->Branch("DecayY", &fPythiaDecayY, "DecayY/D");
fTauTree->Branch("DecayZ", &fPythiaDecayZ, "DecayZ/D");
fTauTree->Branch("Px", &fPythiaPx, "Px/D");
fTauTree->Branch("Py", &fPythiaPy, "Py/D");
fTauTree->Branch("Pz", &fPythiaPz, "Pz/D");
}

//// --- Children of primary tau decay ---
void AnalysisManager::bookCharmTree()
{
fCharmTree = new TTree("charm", "Children of tau decay");
fCharmTree->Branch("evtID", &fPythiaEventID, "evtID/i");
fCharmTree->Branch("PID", &fPythiaPID, "PID/I");
fCharmTree->Branch("TID", &fPythiaTID, "TID/I");
fCharmTree->Branch("PDG", &fPythiaTPDG, "PDG/I");
fCharmTree->Branch("ParentPDG", &fPythiaPPDG, "ParentPDG/I");
fCharmTree->Branch("ProdX", &fPythiaProdX, "ProdX/D");
fCharmTree->Branch("ProdY", &fPythiaProdY, "ProdY/D");
fCharmTree->Branch("ProdZ", &fPythiaProdZ, "ProdZ/D");
fCharmTree->Branch("DecayX", &fPythiaDecayX, "DecayX/D");
fCharmTree->Branch("DecayY", &fPythiaDecayY, "DecayY/D");
fCharmTree->Branch("DecayZ", &fPythiaDecayZ, "DecayZ/D");
fCharmTree->Branch("Px", &fPythiaPx, "Px/D");
fCharmTree->Branch("Py", &fPythiaPy, "Py/D");
fCharmTree->Branch("Pz", &fPythiaPz, "Pz/D");
}

//---------------------------------------------------------------------
//---------------------------------------------------------------------

Expand All @@ -270,6 +314,8 @@ void AnalysisManager::BeginOfRun()
bookPrimTree();
bookGeomTree();
if (fSaveTrack) bookTrkTree();
bookTauTree();
bookCharmTree();

bookHitsTrees();
bookScintTrees();
Expand All @@ -287,6 +333,8 @@ void AnalysisManager::EndOfRun()
fFile->cd();
fEvt->Write();
fPrim->Write();
fTauTree->Write();
fCharmTree->Write();
FillGeomTree();
fGeom->Write();
if (fSaveTrack) fTrk->Write();
Expand Down Expand Up @@ -370,6 +418,10 @@ void AnalysisManager::EndOfEvent(const G4Event *event)

// FILL PRIMARIES/TRAJECTORIES TREE
FillPrimariesTree(event);

fPythiaEventID = evtID;
FillPythiaTree(event, fTauTree, {15});
FillPythiaTree(event, fCharmTree, {411, 421, 431, 4122, 4112, 4212, 4222, 4132, 4232});
if(fSaveTrack) FillTrajectoriesTree(event);

//-----------------------------------------------------------
Expand Down Expand Up @@ -490,11 +542,11 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event)
primVx, primVy, primVz, primVt,
primPx, primPy, primPz,energy));

G4cout << G4endl;
G4cout << "PrimaryParticleInfo: PDG code " << primPDG << G4endl
<< "Particle unique ID : " << primTrackID << G4endl
<< "Momentum : (" << primPx << ", " << primPy << ", " << primPz << ") MeV" << G4endl
<< "Vertex : (" << primVx << ", " << primVy << ", " << primVz << ") mm" << G4endl;
// G4cout << G4endl;
// G4cout << "PrimaryParticleInfo: PDG code " << primPDG << G4endl
// << "Particle unique ID : " << primTrackID << G4endl
// << "Momentum : (" << primPx << ", " << primPy << ", " << primPz << ") MeV" << G4endl
// << "Vertex : (" << primVx << ", " << primVy << ", " << primVz << ") mm" << G4endl;

fPrim->Fill();
}
Expand Down Expand Up @@ -531,17 +583,34 @@ void AnalysisManager::FillTrajectoriesTree(const G4Event* event)
trackNPoints = trajectory->GetPointEntries();
trackTheta = trajectory->GetInitialMomentum().theta();
count_tracks++;
// for (size_t j = 0; j < trackNPoints; ++j)
// {
// G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition();
// trackPointX.push_back( pos.x() );
// trackPointY.push_back( pos.y() );
// trackPointZ.push_back( pos.z() );
// }

// Find start and decay points based on z values
if (trackNPoints > 0) {
G4double minZ = std::numeric_limits<G4double>::max();
G4double maxZ = std::numeric_limits<G4double>::lowest();
G4ThreeVector prodPos, decayPos;

for (size_t j = 0; j < trackNPoints; ++j)
{
G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition();
if (pos.z() < minZ) {
minZ = pos.z();
prodPos = pos;
}
if (pos.z() > maxZ) {
maxZ = pos.z();
decayPos = pos;
}
}

trackProdX = prodPos.x();
trackProdY = prodPos.y();
trackProdZ = prodPos.z();
trackDecayX = decayPos.x();
trackDecayY = decayPos.y();
trackDecayZ = decayPos.z();
}
fTrk->Fill();
// trackPointX.clear();
// trackPointY.clear();
// trackPointZ.clear();
}
G4cout << "Total number of recorded track: " << count_tracks << std::endl;
}
Expand Down Expand Up @@ -721,6 +790,68 @@ void AnalysisManager::FillFaserOutput()
fFaserHitsTree->Fill();
}


void AnalysisManager::FillPythiaTree(const G4Event* event, TTree *tree, std::vector<int> pdg_ids)
{
auto trajectoryContainer = event->GetTrajectoryContainer();
if (!trajectoryContainer)
{
G4cout << "No tracks found: did you enable their storage with '/tracking/storeTrajectory 1'?" << G4endl;
return;
}
for (size_t parentIdx = 0; parentIdx < trajectoryContainer->entries(); ++parentIdx)
{
// Check for primary tau /charm
G4Trajectory *parent = static_cast<G4Trajectory*>((*trajectoryContainer)[parentIdx]);
int parentPDG = std::abs(parent->GetPDGEncoding());
if ((parent->GetParentID() == 0) && (std::find(pdg_ids.begin(), pdg_ids.end(), std::abs(parentPDG)) != pdg_ids.end())) {
G4int parentID = parent->GetTrackID();
fPythiaPPDG = parent->GetPDGEncoding();
// std::cout << "Found primary tau/charm with pdg " << parentPDG << ", track id " << parentID << std::endl;
// Check for children of this primary tauy lepton
for (size_t trajIdx = 0; trajIdx < trajectoryContainer->entries(); ++trajIdx) {
G4Trajectory *traj = static_cast<G4Trajectory*>((*trajectoryContainer)[trajIdx]);
if (traj->GetParentID() == parentID) {
int trajPDG = traj->GetParticleDefinition()->GetPDGEncoding();
// std::cout << "Found child of primary tau/charm with PDG ID " << trajPDG << std::endl;
// Write out trajecotry information
fPythiaTPDG = traj->GetParticleDefinition()->GetPDGEncoding();
fPythiaTID = traj->GetTrackID();
fPythiaPID = traj->GetParentID();
G4int trackNPoints = traj->GetPointEntries();
// Get production and decay vertex
G4double minZ = std::numeric_limits<G4double>::max();
G4double maxZ = std::numeric_limits<G4double>::lowest();
G4ThreeVector prodPos, decayPos;
for (size_t j = 0; j < trackNPoints; ++j)
{
G4ThreeVector pos = traj->GetPoint(j)->GetPosition();
if (pos.z() < minZ) {
minZ = pos.z();
prodPos = pos;
}
if (pos.z() > maxZ) {
maxZ = pos.z();
decayPos = pos;
}
}
fPythiaProdX = prodPos.x();
fPythiaProdY = prodPos.y();
fPythiaProdZ = prodPos.z();
fPythiaDecayX = decayPos.x();
fPythiaDecayY = decayPos.y();
fPythiaDecayZ = decayPos.z();
G4ThreeVector momentum = traj->GetInitialMomentum();
fPythiaPx = momentum.x();
fPythiaPy = momentum.y();
fPythiaPz = momentum.z();
tree->Fill();
}
}
}
}
}

float_t AnalysisManager::GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m)
{
return TMath::Sqrt(px * px + py * py + pz * pz + m * m);
Expand Down
Loading