Exercise Scripts
Last updated on 2026-07-10 | Edit this page
Included below is a selection of scripts for the exercises in part 3 of this tutorial.
You should be able to copy the code text directly into a new file. The name of the file is included as the title of each script section and in the accompanying descriptive text.
ROOT TTreeReader Scripts
EfficiencyAnalysis.C
Create a file called EfficiencyAnalysis.C and copy in
the code below to get started on the efficiency analysis exercise. Note
that you will need to correctly specify your input file path in the
first line.
CPP
void EfficiencyAnalysis(TString infile="PATH_TO_INPUT_FILE"){
// Set output file for the histograms
TFile *ofile = TFile::Open("EfficiencyAnalysis_Out.root","RECREATE");
// Set up input file chain
TChain *mychain = new TChain("events");
mychain->Add(infile);
// Initialize reader
TTreeReader tree_reader(mychain);
// Get Particle Information
TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
TTreeReaderArray<double> partMomX(tree_reader, "MCParticles.momentum.x");
TTreeReaderArray<double> partMomY(tree_reader, "MCParticles.momentum.y");
TTreeReaderArray<double> partMomZ(tree_reader, "MCParticles.momentum.z");
TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
// Get Reconstructed Track Information
TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
// Get Associations Between MCParticles and ReconstructedChargedParticles
TTreeReaderArray<int> recoAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_rec.index");
TTreeReaderArray<int> simuAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_sim.index");
// Define Histograms
TH1D *partEta = new TH1D("partEta","Eta of Thrown Charged Particles;Eta",100,-5.,5.);
TH1D *matchedPartEta = new TH1D("matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track",100,-5.,5.);
TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charged Particle",5000,0.,5.);
while(tree_reader.Next()) { // Loop over events
for(unsigned int i=0; i<partGenStat.GetSize(); i++){ // Loop over thrown particles
if(partGenStat[i] == 1){ // Select stable thrown particles
int pdg = TMath::Abs(partPdg[i]);
if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212){ // Look at charged particles (electrons, muons, pions, kaons, protons)
TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
float trueEta = trueMom.PseudoRapidity();
float truePhi = trueMom.Phi();
partEta->Fill(trueEta);
for(unsigned int j=0; j<simuAssoc.GetSize(); j++){ // Loop over associations to find matching ReconstructedChargedParticle
if(simuAssoc[j] == i){ // Find association index matching the index of the thrown particle we are looking at
TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
// Check the distance between the thrown and reconstructed particle
float deltaEta = trueEta - recMom.PseudoRapidity();
float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
matchedPartTrackDeltaR->Fill(deltaR);
matchedPartEta->Fill(trueEta); // Plot the thrown eta if a matched ReconstructedChargedParticle was found
}
} // End loop over associations
} // End PDG check
} // End stable particles condition
} // End loop over thrown particles
} // End loop over events
ofile->Write(); // Write histograms to file
ofile->Close(); // Close output file
}A “solution” version of the script for the exercise is included below -
CPP
void EfficiencyAnalysis_Exercise(TString infile="PATH_TO_FILE"){
// Set output file for the histograms
TFile *ofile = TFile::Open("EfficiencyAnalysis_Exercise_Out.root","RECREATE");
// Analysis code will go here
// Set up input file chain
TChain *mychain = new TChain("events");
mychain->Add(infile);
// Initialize reader
TTreeReader tree_reader(mychain);
// Get Particle Information
TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
TTreeReaderArray<double> partMomX(tree_reader, "MCParticles.momentum.x");
TTreeReaderArray<double> partMomY(tree_reader, "MCParticles.momentum.y");
TTreeReaderArray<double> partMomZ(tree_reader, "MCParticles.momentum.z");
TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
// Get Reconstructed Track Information
TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
// Get Associations Between MCParticles and ReconstructedChargedParticles
TTreeReaderArray<int> recoAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_rec.index");
TTreeReaderArray<int> simuAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_sim.index");
// Define Histograms
TH1D *partEta = new TH1D("partEta","#eta of Thrown Charged Particles; #eta", 120, -6, 6);
TH1D *matchedPartEta = new TH1D("matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track; #eta", 120, -6, 6);
TH1D* partMom = new TH1D("partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150);
TH1D* matchedPartMom = new TH1D("matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150);
TH1D* partPhi = new TH1D("partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2);
TH1D* matchedPartPhi = new TH1D("matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2);
TH2D* partPEta = new TH2D("partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6);
TH2D* matchedPartPEta = new TH2D("matchedPartPEta", "P vs #eta of Thrown Charged Particles, with matching track; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6);
TH2D* partPhiEta = new TH2D("partPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6);
TH2D* matchedPartPhiEta = new TH2D("matchedPartPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6);
TH1D *matchedPartTrackDeltaEta = new TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25);
TH1D *matchedPartTrackDeltaPhi = new TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2);
TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3);
TH1D *matchedPartTrackDeltaMom = new TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10);
// Define some histograms for our efficiencies
TH1D *TrackEff_Eta = new TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6);
TH1D *TrackEff_Mom = new TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150);
TH1D *TrackEff_Phi = new TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2);
// 2D Efficiencies
TH2D* TrackEff_PEta = new TH2D("TrackEff_PEta", "Tracking efficiency as fn of P and #eta; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6);
TH2D* TrackEff_PhiEta = new TH2D("TrackEff_PhiEta", "Tracking efficiency as fn of #phi and #eta; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6);
// All charged particle histos
TH1D *ChargedEta = new TH1D("ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6);
TH1D *ChargedPhi = new TH1D("ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2);
TH1D *ChargedP = new TH1D("ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150);
while(tree_reader.Next()) { // Loop over events
for(unsigned int i=0; i<partGenStat.GetSize(); i++) // Loop over thrown particles
{
if(partGenStat[i] == 1) // Select stable thrown particles
{
int pdg = TMath::Abs(partPdg[i]);
if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212) // Look at charged particles (electrons, muons, pions, kaons, protons)
{
TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
float trueEta = trueMom.PseudoRapidity();
float truePhi = trueMom.Phi();
partEta->Fill(trueEta);
partPhi->Fill(truePhi);
partMom->Fill(trueMom.Mag());
partPEta->Fill(trueMom.Mag(), trueEta);
partPhiEta->Fill(truePhi, trueEta);
// Loop over associations to find matching ReconstructedChargedParticle
for(unsigned int j=0; j<simuAssoc.GetSize(); j++)
{
if(simuAssoc[j] == i) // Find association index matching the index of the thrown particle we are looking at
{
TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
// Check the distance between the thrown and reconstructed particle
float deltaEta = trueEta - recMom.PseudoRapidity();
float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
float deltaMom = ((trueMom.Mag()) - (recMom.Mag()));
matchedPartTrackDeltaEta->Fill(deltaEta);
matchedPartTrackDeltaPhi->Fill(deltaPhi);
matchedPartTrackDeltaR->Fill(deltaR);
matchedPartTrackDeltaMom->Fill(deltaMom);
matchedPartEta->Fill(trueEta); // Plot the thrown eta if a matched ReconstructedChargedParticle was found
matchedPartPhi->Fill(truePhi);
matchedPartMom->Fill(trueMom.Mag());
matchedPartPEta->Fill(trueMom.Mag(), trueEta);
matchedPartPhiEta->Fill(truePhi, trueEta);
}
}// End loop over associations
} // End PDG check
} // End stable particles condition
} // End loop over thrown particles
// Loop over all charged particles and fill some histograms of kinematics quantities
for(unsigned int k=0; k<trackMomX.GetSize(); k++){ // Loop over all charged particles, thrown or not
TVector3 CPartMom(trackMomX[k], trackMomY[k], trackMomZ[k]);
float CPartEta = CPartMom.PseudoRapidity();
float CPartPhi = CPartMom.Phi();
ChargedEta->Fill(CPartEta);
ChargedPhi->Fill(CPartPhi);
ChargedP->Fill(CPartMom.Mag());
} // End loop over all charged particles
} // End loop over events
// Take the ratio of the histograms above to get our efficiency plots
TrackEff_Eta->Divide(matchedPartEta, partEta, 1, 1, "b");
TrackEff_Mom->Divide(matchedPartMom, partMom, 1, 1, "b");
TrackEff_Phi->Divide(matchedPartPhi, partPhi, 1, 1, "b");
TrackEff_PEta->Divide(matchedPartPEta, partPEta, 1, 1, "b");
TrackEff_PhiEta->Divide(matchedPartPhiEta, partPhiEta, 1, 1, "b");
ofile->Write(); // Write histograms to file
ofile->Close(); // Close output file
}Insert your input file path and execute as the example code above.
ResolutionAnalysis.C
Create a file called ResolutionAnalysis.C and copy in
the code below to get started on the resolution analysis exercise. Note
that you will need to correctly specify your input file path in the
first line.
CPP
void ResolutionAnalysis(TString infile="PATH_TO_INPUT_FILE"){
// Set output file for the histograms
TFile *ofile = TFile::Open("ResolutionAnalysis_Out.root","RECREATE");
// Analysis code will go here
// Set up input file chain
TChain *mychain = new TChain("events");
mychain->Add(infile);
// Initialize reader
TTreeReader tree_reader(mychain);
// Get Particle Information
TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
TTreeReaderArray<double> partMomX(tree_reader, "MCParticles.momentum.x");
TTreeReaderArray<double> partMomY(tree_reader, "MCParticles.momentum.y");
TTreeReaderArray<double> partMomZ(tree_reader, "MCParticles.momentum.z");
TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
// Get Reconstructed Track Information
TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
// Get Associations Between MCParticles and ReconstructedChargedParticles
TTreeReaderArray<int> recoAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_rec.index");
TTreeReaderArray<int> simuAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_sim.index");
// Define Histograms
TH1D *trackMomentumRes = new TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2);
TH1D *matchedPartTrackDeltaEta = new TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25);
TH1D *matchedPartTrackDeltaPhi = new TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2);
TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3);
TH1D *matchedPartTrackDeltaMom = new TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10);
while(tree_reader.Next()) { // Loop over events
for(unsigned int i=0; i<partGenStat.GetSize(); i++){ // Loop over thrown particles
if(partGenStat[i] == 1){ // Select stable thrown particles
int pdg = TMath::Abs(partPdg[i]);
if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212){ // Look at charged particles (electrons, muons, pions, kaons, protons)
TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
float trueEta = trueMom.PseudoRapidity();
float truePhi = trueMom.Phi();
for(unsigned int j=0; j<simuAssoc.GetSize(); j++){ // Loop over associations to find matching ReconstructedChargedParticle
if(simuAssoc[j] == i){ // Find association index matching the index of the thrown particle we are looking at
TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
// Check the distance between the thrown and reconstructed particle
float deltaEta = trueEta - recMom.PseudoRapidity();
float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
float deltaMom = ((trueMom.Mag()) - (recMom.Mag()));
double momRes = (recMom.Mag()- trueMom.Mag())/trueMom.Mag();
trackMomentumRes->Fill(momRes);
matchedPartTrackDeltaEta->Fill(deltaEta);
matchedPartTrackDeltaPhi->Fill(deltaPhi);
matchedPartTrackDeltaR->Fill(deltaR);
matchedPartTrackDeltaMom->Fill(deltaMom);
}
} // End loop over associations
} // End PDG check
} // End stable particles condition
} // End loop over thrown particles
} // End loop over events
ofile->Write(); // Write histograms to file
ofile->Close(); // Close output file
}A “solution” version of the script for the exercise is included below -
CPP
void ResolutionAnalysis_Exercise(TString infile="PATH_TO_FILE"){
// Set output file for the histograms
TFile *ofile = TFile::Open("ResolutionAnalysis_Exercise_Out.root","RECREATE");
// Analysis code will go here
// Set up input file chain
TChain *mychain = new TChain("events");
mychain->Add(infile);
// Initialize reader
TTreeReader tree_reader(mychain);
// Get Particle Information
TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
TTreeReaderArray<double> partMomX(tree_reader, "MCParticles.momentum.x");
TTreeReaderArray<double> partMomY(tree_reader, "MCParticles.momentum.y");
TTreeReaderArray<double> partMomZ(tree_reader, "MCParticles.momentum.z");
TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
// Get Reconstructed Track Information
TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
// Get Associations Between MCParticles and ReconstructedChargedParticles
TTreeReaderArray<int> recoAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_rec.index");
TTreeReaderArray<int> simuAssoc(tree_reader, "_ReconstructedChargedParticleAssociations_sim.index");
// Define Histograms
TH1D *trackMomentumRes = new TH1D("trackMomentumRes","Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
TH2D* trackMomResP = new TH2D("trackMomResP", "Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
TH2D* trackMomResEta = new TH2D("trackMomResEta", "Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
TH1D *trackMomentumRes_e = new TH1D("trackMomentumRes_e","e^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
TH2D* trackMomResP_e = new TH2D("trackMomResP_e", "e^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 25);
TH2D* trackMomResEta_e = new TH2D("trackMomResEta_e", "e^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
TH1D *trackMomentumRes_mu = new TH1D("trackMomentumRes_mu","#mu^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
TH2D* trackMomResP_mu = new TH2D("trackMomResP_mu", "#mu^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 25);
TH2D* trackMomResEta_mu = new TH2D("trackMomResEta_mu", "#mu^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
TH1D *trackMomentumRes_pi = new TH1D("trackMomentumRes_pi","#pi^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
TH2D* trackMomResP_pi = new TH2D("trackMomResP_pi", "#pi^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
TH2D* trackMomResEta_pi = new TH2D("trackMomResEta_pi", "#pi^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
TH1D *trackMomentumRes_K = new TH1D("trackMomentumRes_K","K^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
TH2D* trackMomResP_K = new TH2D("trackMomResP_K", "K^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
TH2D* trackMomResEta_K = new TH2D("trackMomResEta_K", "K^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
TH1D *trackMomentumRes_p = new TH1D("trackMomentumRes_p","p Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
TH2D* trackMomResP_p = new TH2D("trackMomResP_p", "p Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
TH2D* trackMomResEta_p = new TH2D("trackMomResEta_p", "p Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
TH1D *matchedPartTrackDeltaEta = new TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25);
TH1D *matchedPartTrackDeltaPhi = new TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2);
TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3);
TH1D *matchedPartTrackDeltaMom = new TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10);
while(tree_reader.Next()) { // Loop over events
for(unsigned int i=0; i<partGenStat.GetSize(); i++) // Loop over thrown particles
{
if(partGenStat[i] == 1) // Select stable thrown particles
{
int pdg = TMath::Abs(partPdg[i]);
if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212) // Look at charged particles (electrons, muons, pions, kaons, protons)
{
TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
float trueEta = trueMom.PseudoRapidity();
float truePhi = trueMom.Phi();
// Loop over associations to find matching ReconstructedChargedParticle
for(unsigned int j=0; j<simuAssoc.GetSize(); j++)
{
if(simuAssoc[j] == i) // Find association index matching the index of the thrown particle we are looking at
{
TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
// Check the distance between the thrown and reconstructed particle
float deltaEta = trueEta - recMom.PseudoRapidity();
float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
float deltaMom = ((trueMom.Mag()) - (recMom.Mag()));
double momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag();
trackMomentumRes->Fill(momRes); // Could also multiply by 100 and express as a percentage instead
trackMomResP->Fill(momRes, trueMom.Mag());
trackMomResEta->Fill(momRes, trueEta);
if( pdg == 11){
trackMomentumRes_e->Fill(momRes);
trackMomResP_e->Fill(momRes, trueMom.Mag());
trackMomResEta_e->Fill(momRes, trueEta);
}
else if( pdg == 13){
trackMomentumRes_mu->Fill(momRes);
trackMomResP_mu->Fill(momRes, trueMom.Mag());
trackMomResEta_mu->Fill(momRes, trueEta);
}
else if( pdg == 211){
trackMomentumRes_pi->Fill(momRes);
trackMomResP_pi->Fill(momRes, trueMom.Mag());
trackMomResEta_pi->Fill(momRes, trueEta);
}
else if( pdg == 321){
trackMomentumRes_K->Fill(momRes);
trackMomResP_K->Fill(momRes, trueMom.Mag());
trackMomResEta_K->Fill(momRes, trueEta);
}
else if( pdg == 2212){
trackMomentumRes_p->Fill(momRes);
trackMomResP_p->Fill(momRes, trueMom.Mag());
trackMomResEta_p->Fill(momRes, trueEta);
}
matchedPartTrackDeltaEta->Fill(deltaEta);
matchedPartTrackDeltaPhi->Fill(deltaPhi);
matchedPartTrackDeltaR->Fill(deltaR);
matchedPartTrackDeltaMom->Fill(deltaMom);
}
}// End loop over associations
} // End PDG check
} // End stable particles condition
} // End loop over thrown particles
} // End loop over events
ofile->Write(); // Write histograms to file
ofile->Close(); // Close output file
}Insert your input file path and execute as the example code above.
Compiled ROOT Scripts
As brought up in the tutorial, you may wish to compile your ROOT based scripts for faster processing. Included below are some scripts and a short example of a compiled ROOT macro provided by Kolja Kauder.
Each file is uploaded individually, but your directory should be structured as follows -
- helloroot
- README.md
- CMakeLists.txt
- include
- helloroot
- helloroot.hh
- helloroot
- src
- helloexec.cxx
- helloroot.cxx
Note that any entry in the above without a file extension is a directory.
The contents of README.md are -
To build using cmake, create a build directory, navigate to it and run cmake. e.g.:
\```
mkdir build
cd build
cmake ..
make
\```
You can specify a number of parallel build threads with the -j flag, e.g.
\```
make -j4
\```
You can specify an install directory to cmake with
-DCMAKE_INSTALL_PREFIX=<path>
then, after building,
\```
make install
\```
to install the headers and libraries under that location.
There is no "make uninstall" but (on Unix-like systems)
you can do
xargs rm < install_manifest.txt
from the cmake build directory.Note that you should delete the characters in this block.
The contents of CMakeLists.txt are -
CMAKE
# CMakeLists.txt for helloroot.
# More complicated than needed but demonstrates making and linking your own libraries
# cf. https://cliutils.gitlab.io/modern-cmake/
# https://root.cern/manual/integrate_root_into_my_cmake_project/
cmake_minimum_required(VERSION 3.10)
project(helloroot VERSION 1.0 LANGUAGES CXX ) # not needed
# Find ROOT. Use at least 6.20 for smoother cmake support
find_package(ROOT 6.20 REQUIRED )
message ( " ROOT Libraries = " ${ROOT_LIBRARIES} )
##############################################################################################################
# Main target is the libhelloroot library
add_library(
# You can use wildcards but it's cleaner to list the files explicitly
helloroot
SHARED
src/helloroot.cxx
)
## The particular syntax here is a bit annoying because you have to list all the sub-modules you need
## but it picks up automatically all the compile options needed for root, e.g. the c++ std version
## Find all available ROOT modules with `root-config --libs`
target_link_libraries(helloroot PUBLIC ROOT::Core ROOT::RIO ROOT::Rint ROOT::Tree ROOT::EG ROOT::Physics )
## The above _should_ be true, and it is on most systems. If it's not, uncoment one of the following lines
# target_compile_features(helloroot PUBLIC cxx_std_17)
# target_compile_features(helloroot PUBLIC cxx_std_20)
# include directories - this is also overkill but useful if you want to create dictionaries
# Contact kkauder@gmail.com for that - it's too much for this example
target_include_directories(helloroot
PUBLIC
$<INSTALL_INTERFACE:include>
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
)
# Can add addtional options here
target_compile_options(helloroot PRIVATE -Wall -Wextra -pedantic -g)
##############################################################################################################
## Build executables
add_executable(helloexec src/helloexec.cxx)
# target_compile_options(helloexec PRIVATE -Wall -Wextra -pedantic -g)
target_link_libraries(helloexec helloroot )
target_include_directories(helloexec
PRIVATE
${ROOT_INCLUDE_DIRS}
)
install(TARGETS helloexec DESTINATION bin)
##############################################################################################################
## Install library
# Could also use include(GNUInstallDirs)
# and then destinations of the form ${CMAKE_INSTALL_INCLUDEDIR}
install(TARGETS helloroot
EXPORT helloroot-export
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
)
## Install headers
install (DIRECTORY ${CMAKE_SOURCE_DIR}/include/helloroot
DESTINATION include/helloroot
)
## Generate configuration file - this allows you to use cmake in another project
## to find and link the installed helloroot library
install(EXPORT helloroot-export
FILE
hellorootConfig.cmake
NAMESPACE
helloroot::
DESTINATION
cmake
)
## Final message
message( " Done!")The contents of helloroot.hh are -
The contents of helloexec.cxx are -
CPP
#include<helloroot/helloroot.hh>
#include<iostream>
#include<string>
int main()
{
std::cout << "Hello from main " << std::endl;
HelloRoot();
return 0;
}And finally, the contents of helloroot.cxx are -
CPP
#include<helloroot/helloroot.hh>
#include<iostream>
#include<string>
#include<TH1D.h>
#include<TPad.h>
void HelloRoot()
{
std::cout << "Hello from HelloRoot" << std::endl;
// do something with root
TH1D h("h", "h", 100, -5, 5);
h.FillRandom("gaus", 1000);
h.Draw();
gPad->SaveAs("hello.png");
return;
}Please consult the README and script comments for further instructions.
Python Uproot Scripts - Pythonic Versions
Some template scripts that utilise an python array based approach are included below. For some examples of using uproot to access information in .root files, please consult this notebook which can be run in Google Colab.
Pythonic_EfficiencyAnalysis.py
Create a file called EfficiencyAnalysis.py and copy in
the code below to get started on the efficiency analysis exercise. Note
that you will need to correctly specify your input file path in the
variable fname. Note that some example code to process the
division of two histograms is included as a commented section at the end
of this example.
PYTHON
#! /usr/bin/python
# Import some relevant packages
import uproot as up
import awkward as ak
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.ticker as ticker
import matplotlib.cm as cm
import matplotlib.pylab as plt
import scipy, vector, os
from XRootD import client
from scipy import stats
from matplotlib import pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib import colors as colours
# Set some matplot lib features
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['xaxis.labellocation'] = 'right'
plt.rcParams['yaxis.labellocation'] = 'top'
plt.rcParams["figure.figsize"] = (16,9)
kP6 = ['#5790fc','#f89c20','#e42536','#964a8b','#9c9ca1','#7a21dd'] # Set ROOT kP6 colours - see https://root.cern.ch/doc/v636/classTColor.html
# Open our file
fname = "INPUT_FILE.root"
if os.path.isfile(fname):
file=up.open(fname)
else:
print("Error opening file - ", fname, " check your fname variable!")
# Open the tree
tree = file['events']
# Convert relevant branches to arrays
MCPartBr = tree["MCParticles"].arrays()
RecoAssocRec = tree['_ReconstructedChargedParticleAssociations_rec'].arrays()
RecoAssocSim = tree['_ReconstructedChargedParticleAssociations_sim'].arrays()
ReconChPartBr = tree["ReconstructedChargedParticles"].arrays()
RecID=RecoAssocRec['_ReconstructedChargedParticleAssociations_rec.index'] # Array of reconstructed IDs
SimID=RecoAssocSim['_ReconstructedChargedParticleAssociations_sim.index'] # Array of simulated IDs
# Create some filters, anything with [SimID] or [RecID] will index the event by the associations. This means we will only retain events with a matching truth particle/matching reconstructed particle
BoolMatch=(MCPartBr["MCParticles.PDG"][SimID])==(ReconChPartBr["ReconstructedChargedParticles.PDG"][RecID]) # Use simulated or reconstructed IDs as indices, this checks if the pdg between each array matches
BoolChargeTrack = ((abs(MCPartBr["MCParticles.charge"])!=0) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolChargeTrackMatch = ((abs(MCPartBr["MCParticles.charge"][SimID])!=0) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1))
BoolElec=((abs(MCPartBr["MCParticles.PDG"])==11) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolMuon=((abs(MCPartBr["MCParticles.PDG"])==13) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolPion=((abs(MCPartBr["MCParticles.PDG"])==211) & (MCPartBr["MCParticles.generatorStatus"]==1)) # Use abs to include both positive and negative pions
BoolKaon=((abs(MCPartBr["MCParticles.PDG"])==321) & (MCPartBr["MCParticles.generatorStatus"]==1)) # Use abs to include both positive and negative kaons
BoolProton=((abs(MCPartBr["MCParticles.PDG"])==2212) & (MCPartBr["MCParticles.generatorStatus"]==1)) # Use abs to include both positive and negative protons
# Define some doubles
ElecMass = 511*(10**-6) # Electron mass in GeV
# Convert some branches into arrays of vectors
MC_Parts = vector.zip({'px': MCPartBr["MCParticles.momentum.x"], 'py': MCPartBr["MCParticles.momentum.y"], 'pz': MCPartBr["MCParticles.momentum.z"]})
Rec_Parts = vector.zip({'px': ReconChPartBr["ReconstructedChargedParticles.momentum.x"], 'py': ReconChPartBr["ReconstructedChargedParticles.momentum.y"], 'pz': ReconChPartBr["ReconstructedChargedParticles.momentum.z"], 'E':ReconChPartBr["ReconstructedChargedParticles.energy"]})
# Determine the energy for a few MC particles of specific types
MCEnerElec = np.sqrt(MC_Parts[BoolElec].p**2 + ElecMass**2)
# Calculate some additional quantities which are differences between true and reconstructed values for MC particles with a matched reconstructed track
DeltaEta = MC_Parts[SimID][BoolChargeTrackMatch].eta - Rec_Parts[RecID][BoolChargeTrackMatch].eta
DeltaPhi = MC_Parts[SimID][BoolChargeTrackMatch].phi - Rec_Parts[RecID][BoolChargeTrackMatch].phi
DeltaR = np.sqrt(DeltaEta**2 + DeltaPhi**2)
# Plot some quantities as one image
fig, axs = plt.subplots(2,2, tight_layout=True) # Ironically, this makes things *less* tight
axs[-1, -1].axis('off') # Don't draw any blank subfigs
axs[0,0].hist(ak.flatten(MC_Parts[BoolChargeTrack].eta), bins=100, range=(-5,5),alpha=0.5, color=kP6[1]) # Plot the MC eta values for all charged particles at an MC level
axs[0,0].set_title(r"$\eta_{MC}$ of Charged Particles")
axs[0,0].set(xlabel=r'$\eta_{MC}$', ylabel=r'# Entries / 0.1')
axs[0,1].hist(ak.flatten(MC_Parts[SimID][BoolChargeTrackMatch].eta), bins=100, range=(-5,5),alpha=0.5, color=kP6[1]) # Plot the MC eta values for all charged particles at an MC level that have a matching reconstructed track
axs[0,1].set_title(r"$\eta_{MC}$ of matched Charged Particles")
axs[0,1].set(xlabel=r'$\eta_{MC}$', ylabel=r'# Entries / 0.1')
axs[1,0].hist(ak.flatten(DeltaR), bins=5000, range=(0,5),alpha=0.5, color=kP6[1]) # Plot one of our calculated quantities
axs[1,0].set_title(r"$\Delta R$ of Matched Charged Particles")
axs[1,0].set(xlabel=r'$\Delta R$', ylabel=r'# Entries / 0.001')
plt.savefig("EfficiencyAnalysis_Out.png", dpi = (160))
# Commented out, but to divide histograms we can do the following, just put the array we want to plot as the histo in place of Quantity
#MCHist = np.histogram(ak.flatten(Quantity), bins=100, range=(0,25))
#RecHist = np.histogram(ak.flatten(Quantity), bins=100, range=(0,25))
#with np.errstate(divide='ignore'):
# Division = RecHist[0] / MCHist[0]
#Division = np.nan_to_num(Division,nan=0, posinf = 0)
#Bin_Edges=MCHist[1]
#Bars = 0.5 * (Bin_Edges[1:] + Bin_Edges[:-1])
#BarWidth=Bars[1]-Bars[0]
#plt.bar(Bars, Division, width=BarWidth, alpha=0.5, color='kP6[0]')“Complete” exercise example included below
PYTHON
#! /usr/bin/python
# Import some relevant packages
import uproot as up
import awkward as ak
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.ticker as ticker
import matplotlib.cm as cm
import matplotlib.pylab as plt
import scipy, vector, os
from XRootD import client
from scipy import stats
from matplotlib import pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib import colors as colours
# Set some matplot lib features
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['xaxis.labellocation'] = 'right'
plt.rcParams['yaxis.labellocation'] = 'top'
plt.rcParams["figure.figsize"] = (16,9)
kP6 = ['#5790fc','#f89c20','#e42536','#964a8b','#9c9ca1','#7a21dd'] # Set ROOT kP6 colours - see https://root.cern.ch/doc/v636/classTColor.html
# Open our file
fname = "YOUR_INPUT_FILE"
if os.path.isfile(fname):
file=up.open(fname)
else:
print("Error opening file - ", fname, " check your fname variable!")
# Open the tree
tree = file['events']
# Convert relevant branches to arrays
MCPartBr = tree["MCParticles"].arrays()
RecoAssocRec = tree['_ReconstructedChargedParticleAssociations_rec'].arrays()
RecoAssocSim = tree['_ReconstructedChargedParticleAssociations_sim'].arrays()
ReconChPartBr = tree["ReconstructedChargedParticles"].arrays()
RecID=RecoAssocRec['_ReconstructedChargedParticleAssociations_rec.index'] # Array of reconstructed IDs
SimID=RecoAssocSim['_ReconstructedChargedParticleAssociations_sim.index'] # Array of simulated IDs
# Create some filters, anything with [SimID] or [RecID] will index the event by the associations. This means we will only retain events with a matching truth particle/matching reconstructed particle
BoolMatch=(MCPartBr["MCParticles.PDG"][SimID])==(ReconChPartBr["ReconstructedChargedParticles.PDG"][RecID]) # Use simulated or reconstructed IDs as indices, this checks if the pdg between each array matches
BoolChargeTrack = ((abs(MCPartBr["MCParticles.charge"])!=0) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolChargeTrackMatch = ((abs(MCPartBr["MCParticles.charge"][SimID])!=0) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1))
BoolElec=((abs(MCPartBr["MCParticles.PDG"])==11) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolMuon=((abs(MCPartBr["MCParticles.PDG"])==13) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolPion=((abs(MCPartBr["MCParticles.PDG"])==211) & (MCPartBr["MCParticles.generatorStatus"]==1)) # Use abs to include both positive and negative pions
BoolKaon=((abs(MCPartBr["MCParticles.PDG"])==321) & (MCPartBr["MCParticles.generatorStatus"]==1)) # Use abs to include both positive and negative kaons
BoolProton=((abs(MCPartBr["MCParticles.PDG"])==2212) & (MCPartBr["MCParticles.generatorStatus"]==1)) # Use abs to include both positive and negative protons
BoolElecMatch=((MCPartBr["MCParticles.PDG"][SimID]==11) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1))
BoolPionMatch=((abs(MCPartBr["MCParticles.PDG"][SimID])==211) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1)) # Use abs to include both positive and negative pions
MCStatus = MCPartBr['MCParticles.generatorStatus'] == 1
MCNegCharge = MCPartBr['MCParticles.charge'] == -1
MCScElecPDG = MCPartBr['MCParticles.PDG'] == 11
MCPionPDG = abs(MCPartBr['MCParticles.PDG']) == 211
# We index by the simulation ID to ONLY select events with a matching track
# Define some doubles
ElecMass = 511*(10**-6) # Electron mass in GeV
# Convert some branches into arrays of vectors
MC_Parts = vector.zip({'px': MCPartBr["MCParticles.momentum.x"], 'py': MCPartBr["MCParticles.momentum.y"], 'pz': MCPartBr["MCParticles.momentum.z"]})
Rec_Parts = vector.zip({'px': ReconChPartBr["ReconstructedChargedParticles.momentum.x"], 'py': ReconChPartBr["ReconstructedChargedParticles.momentum.y"], 'pz': ReconChPartBr["ReconstructedChargedParticles.momentum.z"], 'E':ReconChPartBr["ReconstructedChargedParticles.energy"]})
Rec_Vects = vector.zip({'px': MCPartBr["MCParticles.momentum.x"][SimID], 'py': MCPartBr["MCParticles.momentum.y"][SimID], 'pz': MCPartBr["MCParticles.momentum.z"][SimID]})
MC_ScElec = MC_Parts[MCStatus & MCNegCharge & MCScElecPDG]
MC_Pions = MC_Parts[MCStatus & MCPionPDG]
Rec_ScElec = Rec_Vects[BoolElecMatch]
Rec_Pions = Rec_Vects[BoolPionMatch]
# Determine the energy for a few MC particles of specific types
MCEnerElec = np.sqrt(MC_Parts[BoolElec].p**2 + ElecMass**2)
# Calculate some additional quantities which are differences between true and reconstructed values for MC particles with a matched reconstructed track
DeltaEta = MC_Parts[SimID][BoolChargeTrackMatch].eta - Rec_Parts[RecID][BoolChargeTrackMatch].eta
DeltaPhi = MC_Parts[SimID][BoolChargeTrackMatch].phi - Rec_Parts[RecID][BoolChargeTrackMatch].phi
DeltaR = np.sqrt(DeltaEta**2 + DeltaPhi**2)
# Plot some quantities as one image
fig, axs = plt.subplots(2,2, tight_layout=True) # Ironically, this makes things *less* tight
axs[-1, -1].axis('off') # Don't draw any blank subfigs
axs[0,0].hist(ak.flatten(MC_Parts[BoolChargeTrack].eta), bins=100, range=(-5,5),alpha=0.5, color=kP6[1]) # Plot the MC eta values for all charged particles at an MC level
axs[0,0].set_title(r"$\eta_{MC}$ of Charged Particles")
axs[0,0].set(xlabel=r'$\eta_{MC}$', ylabel=r'# Entries / 0.1')
axs[0,1].hist(ak.flatten(MC_Parts[SimID][BoolChargeTrackMatch].eta), bins=100, range=(-5,5),alpha=0.5, color=kP6[1]) # Plot the MC eta values for all charged particles at an MC level that have a matching reconstructed track
axs[0,1].set_title(r"$\eta_{MC}$ of matched Charged Particles")
axs[0,1].set(xlabel=r'$\eta_{MC}$', ylabel=r'# Entries / 0.1')
axs[1,0].hist(ak.flatten(DeltaR), bins=5000, range=(0,5),alpha=0.5, color=kP6[1]) # Plot one of our calculated quantities
axs[1,0].set_title(r"$\Delta R$ of Matched Charged Particles")
axs[1,0].set(xlabel=r'$\Delta R$', ylabel=r'# Entries / 0.001')
plt.savefig("EfficiencyAnalysis_Out.png", dpi = (160))
bins_eta=100
range_eta=(-5,5)
bins_p=100
range_p_elec=(0,25)
range_p_pi=(0,50)
# Make histograms of our scattered electrons and pions - Full truth distributions in P and eta
MC_ScElec_eta = np.histogram(ak.flatten(MC_ScElec.eta), bins = bins_eta, range= range_eta)
MC_Pions_eta = np.histogram(ak.flatten(MC_Pions.eta), bins = bins_eta, range= range_eta)
MC_ScElec_P = np.histogram(ak.flatten(MC_ScElec.p), bins = bins_p, range= range_p_elec)
MC_Pions_P = np.histogram(ak.flatten(MC_Pions.p), bins = bins_p, range= range_p_pi)
# Make histograms of our scattered electrons and pions - Truth distributions for particles that reconstructed in P and eta
Rec_ScElec_eta = np.histogram(ak.flatten(Rec_ScElec.eta), bins = bins_eta, range= range_eta)
Rec_Pions_eta = np.histogram(ak.flatten(Rec_Pions.eta), bins = bins_eta, range= range_eta)
Rec_ScElec_P = np.histogram(ak.flatten(Rec_ScElec.p), bins = bins_p, range= range_p_elec)
Rec_Pions_P = np.histogram(ak.flatten(Rec_Pions.p), bins = bins_p, range= range_p_pi)
# Divide to get efficiency
with np.errstate(divide='ignore'):
Eff_ScElec_eta = Rec_ScElec_eta[0]/MC_ScElec_eta[0]
Eff_Pion_eta = Rec_Pions_eta[0]/MC_Pions_eta[0]
Eff_ScElec_P = Rec_ScElec_P[0]/MC_ScElec_P[0]
Eff_Pion_P = Rec_Pions_P[0]/MC_Pions_P[0]
Eff_ScElec_eta=np.nan_to_num(Eff_ScElec_eta,nan=0,posinf=0)
Eff_Pion_eta=np.nan_to_num(Eff_Pion_eta,nan=0,posinf=0)
Eff_ScElec_P=np.nan_to_num(Eff_ScElec_P,nan=0,posinf=0)
Eff_Pion_P=np.nan_to_num(Eff_Pion_P,nan=0,posinf=0)
fig, axs = plt.subplots(2,2, tight_layout=True) # Ironically, this makes things *less* tight
#axs[-1, -1].axis('off') # Don't draw any blank subfigs
Bin_Edges=Rec_ScElec_eta[1]
Bars=0.5 * (Bin_Edges[1:] + Bin_Edges[:-1])
BarWidth=Bars[1]-Bars[0]
axs[0,0].bar(Bars, Eff_ScElec_eta, width=BarWidth, alpha=0.75, color=kP6[1]) # Plot the MC eta values for all charged particles at an MC level
axs[0,0].set_title(r"Reconstructed $e'$ Efficiency as fn of $\eta$")
axs[0,0].set(xlabel=r'"$\eta$', ylabel=r"$e'$ Effiency")
Bin_Edges=Rec_ScElec_P[1]
Bars=0.5 * (Bin_Edges[1:] + Bin_Edges[:-1])
BarWidth=Bars[1]-Bars[0]
axs[0,1].bar(Bars, Eff_ScElec_P, width=BarWidth, alpha=0.75, color=kP6[1]) # Plot the MC eta values for all charged particles at an MC level that have a matching reconstructed track
axs[0,1].set_title(r"Reconstructed $e'$ Efficiency as fn of $P$")
axs[0,1].set(xlabel=r"$P_{e'}$", ylabel=r"$P_{e'}$")
Bin_Edges=Rec_Pions_eta[1]
Bars=0.5 * (Bin_Edges[1:] + Bin_Edges[:-1])
BarWidth=Bars[1]-Bars[0]
axs[1,0].bar(Bars, Eff_Pion_eta, width=BarWidth, alpha=0.75, color=kP6[1]) # Plot one of our calculated quantities
axs[1,0].set_title(r"Reconstructed $\pi$ Efficiency as fn of $\eta$")
axs[1,0].set(xlabel=r"$\eta$", ylabel=r"$\pi$ Effiency")
Bin_Edges=Rec_Pions_P[1]
Bars=0.5 * (Bin_Edges[1:] + Bin_Edges[:-1])
BarWidth=Bars[1]-Bars[0]
axs[1,1].bar(Bars, Eff_Pion_P, width=BarWidth, alpha=0.75, color=kP6[1]) # Plot one of our calculated quantities
axs[1,1].set_title(r"Reconstructed $\pi$ Efficiency as fn of $P$")
axs[1,1].set(xlabel=r"$P_{\pi}$", ylabel=r"$\pi$ Effiency")
plt.savefig("EfficiencyAnalysis_Exercise_Out.png", dpi = (160))Pythonic_ResolutionAnalysis.py
Create a file called ResolutionAnalysis.py and copy in
the code below to get started on the efficiency analysis exercise. Note
that you will need to correctly specify your input file path in the
variable fname.
PYTHON
#! /usr/bin/python
# Import some relevant packages
import uproot as up
import awkward as ak
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.ticker as ticker
import matplotlib.cm as cm
import matplotlib.pylab as plt
import scipy, vector, os
from XRootD import client
from scipy import stats
from matplotlib import pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib import colors as colours
# Set some matplot lib features
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['xaxis.labellocation'] = 'right'
plt.rcParams['yaxis.labellocation'] = 'top'
plt.rcParams["figure.figsize"] = (16,9)
kP6 = ['#5790fc','#f89c20','#e42536','#964a8b','#9c9ca1','#7a21dd'] # Set ROOT kP6 colours - see https://root.cern.ch/doc/v636/classTColor.html
# Open our file
fname = "INPUT_FILE.root"
if os.path.isfile(fname):
file=up.open(fname)
else:
print("Error opening file - ", fname, " check your fname variable!")
# Open the tree
tree = file['events']
# Convert relevant branches to arrays
MCPartBr = tree["MCParticles"].arrays()
RecoAssocRec = tree['_ReconstructedChargedParticleAssociations_rec'].arrays()
RecoAssocSim = tree['_ReconstructedChargedParticleAssociations_sim'].arrays()
ReconChPartBr = tree["ReconstructedChargedParticles"].arrays()
RecID=RecoAssocRec['_ReconstructedChargedParticleAssociations_rec.index'] # Array of reconstructed IDs
SimID=RecoAssocSim['_ReconstructedChargedParticleAssociations_sim.index'] # Array of simulated IDs
# Create some filters, anything with [SimID] or [RecID] will index the event by the associations. This means we will only retain events with a matching truth particle/matching reconstructed particle
BoolChargeTrack = ((abs(MCPartBr["MCParticles.charge"])!=0) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolChargeTrackMatch = ((abs(MCPartBr["MCParticles.charge"][SimID])!=0) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1))
BoolElec=((abs(MCPartBr["MCParticles.PDG"])==11) & (MCPartBr["MCParticles.generatorStatus"]==1))
# Define some doubles
ElecMass = 511*(10**-6) # Electron mass in GeV
# Convert some branches into arrays of vectors
MC_Parts = vector.zip({'px': MCPartBr["MCParticles.momentum.x"], 'py': MCPartBr["MCParticles.momentum.y"], 'pz': MCPartBr["MCParticles.momentum.z"]})
Rec_Parts = vector.zip({'px': ReconChPartBr["ReconstructedChargedParticles.momentum.x"], 'py': ReconChPartBr["ReconstructedChargedParticles.momentum.y"], 'pz': ReconChPartBr["ReconstructedChargedParticles.momentum.z"], 'E':ReconChPartBr["ReconstructedChargedParticles.energy"]})
# Determine the energy for a few MC particles of specific types
MCEnerElec = np.sqrt(MC_Parts[BoolElec].p**2 + ElecMass**2)
# Calculate some additional quantities which are differences between true and reconstructed values for MC particles with a matched reconstructed track
DeltaEta = MC_Parts[SimID][BoolChargeTrackMatch].eta - Rec_Parts[RecID][BoolChargeTrackMatch].eta
DeltaPhi = MC_Parts[SimID][BoolChargeTrackMatch].phi - Rec_Parts[RecID][BoolChargeTrackMatch].phi
DeltaP = MC_Parts[SimID][BoolChargeTrackMatch].p - Rec_Parts[RecID][BoolChargeTrackMatch].p
DeltaR = np.sqrt(DeltaEta**2 + DeltaPhi**2)
ResP = ((Rec_Parts[RecID][BoolChargeTrackMatch].p - MC_Parts[SimID][BoolChargeTrackMatch].p)/MC_Parts[SimID][BoolChargeTrackMatch].p ) # Momentum resolution as a percentage
# Plot some quantities as one image
fig, axs = plt.subplots(2,3, tight_layout=True) # Ironically, this makes things *less* tight
axs[-1, -1].axis('off') # Don't draw any blank subfigs
axs[0,0].hist(ak.flatten(DeltaEta), bins=100, range=(-0.25,0.25),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed eta
axs[0,0].set_title(r"$\Delta \eta$ of Matched Charged Particles")
axs[0,0].set(xlabel=r'$\Delta \eta$', ylabel=r'# Entries / 0.005')
axs[0,1].hist(ak.flatten(DeltaPhi), bins=200, range=(-0.2,0.2),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed phi
axs[0,1].set_title(r"$\Delta \phi$ of Matched Charged Particles")
axs[0,1].set(xlabel=r'$\Delta \phi$', ylabel='# Entries / 0.002')
axs[0,2].hist(ak.flatten(DeltaR), bins=300, range=(0,0.3),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed R
axs[0,2].set_title(r"$\Delta R$ of Matched Charged Particles")
axs[0,2].set(xlabel=r'$\Delta R$', ylabel=r'# Entries / 0.003')
axs[1,0].hist(ak.flatten(DeltaP), bins=200, range=(-10,10),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed momentum
axs[1,0].set_title(r"$\Delta P$ of Matched Charged Particles")
axs[1,0].set(xlabel=r'$\Delta \eta$', ylabel=r'# Entries / 0.1 GeV/c')
axs[1,1].hist(ak.flatten(ResP), bins=400, range=(-2,2),alpha=0.5, color=kP6[1]) # Plot the momentum resolution
axs[1,1].set_title(r"Momentum Resolution of Matched Charged Particles")
axs[1,1].set(xlabel=r'$(P_{Rec} - P_{MC})/P_{MC}$', ylabel=r'# Entries / 0.01')
plt.savefig("ResolutionAnalysis_Out.png", dpi = (160))“Complete” exercise example included below
PYTHON
#! /usr/bin/python
# Import some relevant packages
import uproot as up
import awkward as ak
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.ticker as ticker
import matplotlib.cm as cm
import matplotlib.pylab as plt
import scipy, vector, os
from XRootD import client
from scipy import stats
from matplotlib import pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib import colors as colours
# Set some matplot lib features
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['xaxis.labellocation'] = 'right'
plt.rcParams['yaxis.labellocation'] = 'top'
plt.rcParams["figure.figsize"] = (16,9)
kP6 = ['#5790fc','#f89c20','#e42536','#964a8b','#9c9ca1','#7a21dd'] # Set ROOT kP6 colours - see https://root.cern.ch/doc/v636/classTColor.html
# Open our file
fname = "YOUR_INPUT_FILE"
if os.path.isfile(fname):
file=up.open(fname)
else:
print("Error opening file - ", fname, " check your fname variable!")
# Open the tree
tree = file['events']
# Convert relevant branches to arrays
MCPartBr = tree["MCParticles"].arrays()
RecoAssocRec = tree['_ReconstructedChargedParticleAssociations_rec'].arrays()
RecoAssocSim = tree['_ReconstructedChargedParticleAssociations_sim'].arrays()
ReconChPartBr = tree["ReconstructedChargedParticles"].arrays()
RecID=RecoAssocRec['_ReconstructedChargedParticleAssociations_rec.index'] # Array of reconstructed IDs
SimID=RecoAssocSim['_ReconstructedChargedParticleAssociations_sim.index'] # Array of simulated IDs
# Create some filters, anything with [SimID] or [RecID] will index the event by the associations. This means we will only retain events with a matching truth particle/matching reconstructed particle
MCStatus = MCPartBr['MCParticles.generatorStatus'] == 1
MCNegCharge = MCPartBr['MCParticles.charge'] == -1
MCScElecPDG = MCPartBr['MCParticles.PDG'] == 11
MCPionPDG = abs(MCPartBr['MCParticles.PDG']) == 211
BoolChargeTrack = ((abs(MCPartBr["MCParticles.charge"])!=0) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolChargeTrackMatch = ((abs(MCPartBr["MCParticles.charge"][SimID])!=0) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1))
BoolElec=((abs(MCPartBr["MCParticles.PDG"])==11) & (MCPartBr["MCParticles.generatorStatus"]==1))
BoolElecMatch=((MCPartBr["MCParticles.PDG"][SimID]==11) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1))
BoolPionMatch=((abs(MCPartBr["MCParticles.PDG"][SimID])==211) & (MCPartBr["MCParticles.generatorStatus"][SimID]==1)) # Use abs to include both positive and negative pions
# Define some doubles
ElecMass = 511*(10**-6) # Electron mass in GeV
# Convert some branches into arrays of vectors
MC_Parts = vector.zip({'px': MCPartBr["MCParticles.momentum.x"], 'py': MCPartBr["MCParticles.momentum.y"], 'pz': MCPartBr["MCParticles.momentum.z"]})
Rec_Parts = vector.zip({'px': ReconChPartBr["ReconstructedChargedParticles.momentum.x"], 'py': ReconChPartBr["ReconstructedChargedParticles.momentum.y"], 'pz': ReconChPartBr["ReconstructedChargedParticles.momentum.z"], 'E':ReconChPartBr["ReconstructedChargedParticles.energy"]})
# We redfine MC vects here to ONLY be the MC particles that have a reconstructed track
MC_Vects = vector.zip({'px': MCPartBr["MCParticles.momentum.x"][SimID], 'py': MCPartBr["MCParticles.momentum.y"][SimID], 'pz': MCPartBr["MCParticles.momentum.z"][SimID]})
MC_ScElec = MC_Vects[BoolElecMatch]
MC_Pions = MC_Vects[BoolPionMatch]
# In this case, we need to access our reconstructed charged particles branch and index it by the reconstructed ID.
RecVects = vector.zip({'px': ReconChPartBr["ReconstructedChargedParticles.momentum.x"][RecID], 'py': ReconChPartBr["ReconstructedChargedParticles.momentum.y"][RecID], 'pz': ReconChPartBr["ReconstructedChargedParticles.momentum.z"][RecID]})
Rec_ScElec = RecVects[BoolElecMatch]
Rec_Pions = RecVects[BoolPionMatch]
ElecMomRes = ((Rec_ScElec.p - MC_ScElec.p)/MC_ScElec.p)*100
ElecEtaRes = ((Rec_ScElec.eta - MC_ScElec.eta)/MC_ScElec.eta)*100
PiMomRes = ((Rec_Pions.p - MC_Pions.p)/MC_Pions.p)*100
PiEtaRes = ((Rec_Pions.eta - MC_Pions.eta)/MC_Pions.eta)*100
# Determine the energy for a few MC particles of specific types
MCEnerElec = np.sqrt(MC_Parts[BoolElec].p**2 + ElecMass**2)
# Calculate some additional quantities which are differences between true and reconstructed values for MC particles with a matched reconstructed track
DeltaEta = MC_Parts[SimID][BoolChargeTrackMatch].eta - Rec_Parts[RecID][BoolChargeTrackMatch].eta
DeltaPhi = MC_Parts[SimID][BoolChargeTrackMatch].phi - Rec_Parts[RecID][BoolChargeTrackMatch].phi
DeltaP = MC_Parts[SimID][BoolChargeTrackMatch].p - Rec_Parts[RecID][BoolChargeTrackMatch].p
DeltaR = np.sqrt(DeltaEta**2 + DeltaPhi**2)
ResP = ((Rec_Parts[RecID][BoolChargeTrackMatch].p - MC_Parts[SimID][BoolChargeTrackMatch].p)/MC_Parts[SimID][BoolChargeTrackMatch].p ) # Momentum resolution as a percentage
# Plot some quantities as one image
fig, axs = plt.subplots(2,3, tight_layout=True) # Ironically, this makes things *less* tight
axs[-1, -1].axis('off') # Don't draw any blank subfigs
axs[0,0].hist(ak.flatten(DeltaEta), bins=100, range=(-0.25,0.25),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed eta
axs[0,0].set_title(r"$\Delta \eta$ of Matched Charged Particles")
axs[0,0].set(xlabel=r'$\Delta \eta$', ylabel=r'# Entries / 0.005')
axs[0,1].hist(ak.flatten(DeltaPhi), bins=200, range=(-0.2,0.2),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed phi
axs[0,1].set_title(r"$\Delta \phi$ of Matched Charged Particles")
axs[0,1].set(xlabel=r'$\Delta \phi$', ylabel='# Entries / 0.002')
axs[0,2].hist(ak.flatten(DeltaR), bins=300, range=(0,0.3),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed R
axs[0,2].set_title(r"$\Delta R$ of Matched Charged Particles")
axs[0,2].set(xlabel=r'$\Delta R$', ylabel=r'# Entries / 0.003')
axs[1,0].hist(ak.flatten(DeltaP), bins=200, range=(-10,10),alpha=0.5, color=kP6[1]) # Plot the difference between true and reconstructed momentum
axs[1,0].set_title(r"$\Delta P$ of Matched Charged Particles")
axs[1,0].set(xlabel=r'$\Delta \eta$', ylabel=r'# Entries / 0.1 GeV/c')
axs[1,1].hist(ak.flatten(ResP), bins=400, range=(-2,2),alpha=0.5, color=kP6[1]) # Plot the momentum resolution
axs[1,1].set_title(r"Momentum Resolution of Matched Charged Particles")
axs[1,1].set(xlabel=r'$(P_{Rec} - P_{MC})/P_{MC}$', ylabel=r'# Entries / 0.01')
plt.savefig("ResolutionAnalysis_Out.png", dpi = (160))
fig, axs = plt.subplots(2,3, tight_layout=True) # Ironically, this makes things *less* tight
axs[0,0].hist(ak.flatten(ElecMomRes), bins=50, range=(-25,25),alpha=0.5, color=kP6[1], label=r"$e'$")
axs[0,0].hist(ak.flatten(PiMomRes), bins=50, range=(-25,25),alpha=0.25, color=kP6[0], label=r"$\pi$")
axs[0,0].set_title(r"Reconstructed $P$ Resolution")
axs[0,0].set(xlabel=r"$P$ Resolution [%]", ylabel=r'# Entries')
axs[0,0].legend(loc='upper right')
axs[0,1].hist(ak.flatten(ElecEtaRes), bins=100, range=(-2,2),alpha=0.5, color=kP6[1], label=r"$e'$")
axs[0,1].hist(ak.flatten(PiEtaRes), bins=100, range=(-2,2),alpha=0.25, color=kP6[0], label=r"$\pi$")
axs[0,1].set_title(r"Reconstructed $\eta$ Resolution")
axs[0,1].set(xlabel=r"$\eta$ Resolution [%]", ylabel=r'# Entries')
axs[0,1].legend(loc='upper right')
# 2D plots
Hist2D1 = axs[0,2].hist2d(np.asarray(ak.flatten(ElecMomRes)), np.asarray(ak.flatten(MC_ScElec.p)), bins=[50,50], range=[[-25,25],[0,20]], cmin=1)
axs[0,2].set_title(r"Reconstructed $P_{e'}$ Resolution as a function of $P_{e'MC}$")
axs[0,2].set(xlabel=r"$P_{e'}$ Resolution [%]", ylabel=r"$P_{e'MC}$")
cb1=plt.colorbar(Hist2D1[3],ax=axs[0,2]) # [3] is the z axis info
cb1.set_label('Counts/bin')
Hist2D2=axs[1,0].hist2d(np.asarray(ak.flatten(PiMomRes)), np.asarray(ak.flatten(MC_Pions.p)), bins=[50,50], range=[[-25,25],[0,20]], cmin=1)
axs[1,0].set_title(r"Reconstructed $P_{\pi}$ Resolution as a function of $P_{\pi MC}$")
axs[1,0].set(xlabel=r"$P_{\pi}$ Resolution [%]", ylabel=r"$P_{\pi'MC}$")
cb2=plt.colorbar(Hist2D2[3],ax=axs[1,0]) # [3] is the z axis info
cb2.set_label('Counts/bin')
Hist2D3=axs[1,1].hist2d(np.asarray(ak.flatten(ElecEtaRes)), np.asarray(ak.flatten(MC_ScElec.eta)), bins=[100,100], range=[[-2,2],[-5,5]], cmin=1)
axs[1,1].set_title(r"Reconstructed $\eta_{e'}$ Resolution as a function of $\eta_{e'MC}$")
axs[1,1].set(xlabel=r"$\eta_{e'}$ Resolution [%]", ylabel=r"$\eta_{e'MC}$")
cb3=plt.colorbar(Hist2D3[3],ax=axs[1,1]) # [3] is the z axis info
cb3.set_label('Counts/bin')
Hist2D4=axs[1,2].hist2d(np.asarray(ak.flatten(PiEtaRes)), np.asarray(ak.flatten(MC_Pions.eta)), bins=[100,100], range=[[-2,2],[-5,5]], cmin=1)
axs[1,2].set_title(r"Reconstructed $\eta_{e'}$ Resolution as a function of $\eta_{e'MC}$")
axs[1,2].set(xlabel=r"$\eta_{\pi}$ Resolution [%]", ylabel=r"$\eta_{\pi MC}$")
cb4=plt.colorbar(Hist2D4[3],ax=axs[1,2]) # [3] is the z axis info
cb4.set_label('Counts/bin')
plt.savefig("ResolutionAnalysis_Exercise_Out.png", dpi = (160))Python Uproot Script - C/ROOT Style (Slow, not recommended!)
EfficiencyAnalysis.py
Create a file called EfficiencyAnalysis.py and copy in
the code below to get started on the efficiency analysis exercise. Note
that you will need to correctly specify your input file path in the
variable infile.
PYTHON
#! /usr/bin/python
#Import relevant packages
import ROOT, math, array
from ROOT import TH1F, TH2F, TMath, TTree, TVector3, TVector2
import uproot as up
#Define and open files
infile="PATH_TO_INPUT_FILE"
ofile=ROOT.TFile.Open("EfficiencyAnalysis_OutPy.root", "RECREATE")
# Open input file and define branches we want to look at with uproot
events_tree = up.open(infile)["events"]
# Get particle information
partGenStat = events_tree["MCParticles.generatorStatus"].array()
partMomX = events_tree["MCParticles.momentum.x"].array()
partMomY = events_tree["MCParticles.momentum.y"].array()
partMomZ = events_tree["MCParticles.momentum.z"].array()
partPdg = events_tree["MCParticles.PDG"].array()
# Get reconstructed track information
trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
# Get assocations between MCParticles and ReconstructedChargedParticles
recoAssoc = events_tree["_ReconstructedChargedParticleAssociations_rec.index"].array()
simuAssoc = events_tree["_ReconstructedChargedParticleAssociations_sim.index"].array()
# Define histograms below
partEta = ROOT.TH1D("partEta","Eta of Thrown Charged Particles;Eta",100, -5 ,5 )
matchedPartEta = ROOT.TH1D("matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track", 100, -5 ,5);
matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charge Particle", 5000, 0, 5);
# Add main analysis loop(s) below
for i in range(0, len(partGenStat)): # Loop over all events
for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
if partGenStat[i][j] == 1: # Select stable particles
pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
trueEta = trueMom.PseudoRapidity()
truePhi = trueMom.Phi()
partEta.Fill(trueEta)
for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
if (simuAssoc[i][k] == j):
recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
deltaEta = trueEta - recMom.PseudoRapidity()
deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
matchedPartEta.Fill(trueEta)
matchedPartTrackDeltaR.Fill(deltaR)
# Write output histograms to file below
partEta.Write()
matchedPartEta.Write()
matchedPartTrackDeltaR.Write()
# Close files
ofile.Close()A “solution” version of the script for the exercise is included below -
PYTHON
#! /usr/bin/python
#Import relevant packages
import ROOT, math, array
from ROOT import TCanvas, TColor, TGaxis, TH1F, TH2F, TPad, TStyle, gStyle, gPad, TGaxis, TLine, TMath, TPaveText, TTree, TVector3, TVector2
import uproot as up
#Define and open files
infile="PATH_TO_FILE"
ofile=ROOT.TFile.Open("EfficiencyAnalysis_Exercise_OutPy.root", "RECREATE")
# Open input file and define branches we want to look at with uproot
events_tree = up.open(infile)["events"]
# Get particle information
partGenStat = events_tree["MCParticles.generatorStatus"].array()
partMomX = events_tree["MCParticles.momentum.x"].array()
partMomY = events_tree["MCParticles.momentum.y"].array()
partMomZ = events_tree["MCParticles.momentum.z"].array()
partPdg = events_tree["MCParticles.PDG"].array()
# Get reconstructed track information
trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
# Get assocations between MCParticles and ReconstructedChargedParticles
recoAssoc = events_tree["_ReconstructedChargedParticleAssociations_rec.index"].array()
simuAssoc = events_tree["_ReconstructedChargedParticleAssociations_sim.index"].array()
# Define histograms below
partEta = ROOT.TH1D("partEta","#eta of Thrown Charged Particles; #eta", 120, -6, 6)
matchedPartEta = ROOT.TH1D("matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track; #eta", 120, -6, 6)
partMom = ROOT.TH1D("partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150)
matchedPartMom = ROOT.TH1D("matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150)
partPhi = ROOT.TH1D("partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2)
matchedPartPhi = ROOT.TH1D("matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2)
partPEta = ROOT.TH2D("partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6)
matchedPartPEta = ROOT.TH2D("matchedPartPEta", "P vs #eta of Thrown Charged Particles, with matching track; P(GeV/C); #eta", 150, 0, 150, 120, -6, 6)
partPhiEta = ROOT.TH2D("partPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
matchedPartPhiEta = ROOT.TH2D("matchedPartPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
# Define some histograms for our efficiencies
TrackEff_Eta = ROOT.TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6)
TrackEffMom = ROOT.TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150)
TrackEffPhi = ROOT.TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2)
# 2D Efficiencies
TrackEff_PEta = ROOT.TH2D("TrackEff_PEta", "Tracking efficiency as fn of P and #eta; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6)
TrackEff_PhiEta = ROOT.TH2D("TrackEff_PhiEta", "Tracking efficiency as fn of #phi and #eta; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
# All charged particle histos
ChargedEta = ROOT.TH1D("ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6)
ChargedPhi = ROOT.TH1D("ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2)
ChargedP = ROOT.TH1D("ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150)
# Add main analysis loop(s) below
for i in range(0, len(partGenStat)): # Loop over all events
for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
if partGenStat[i][j] == 1: # Select stable particles
pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
trueEta = trueMom.PseudoRapidity()
truePhi = trueMom.Phi()
partEta.Fill(trueEta)
partPhi.Fill(truePhi)
partMom.Fill(trueMom.Mag())
partPEta.Fill(trueMom.Mag(), trueEta)
partPhiEta.Fill(truePhi, trueEta)
for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
if (simuAssoc[i][k] == j):
recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
deltaEta = trueEta - recMom.PseudoRapidity()
deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
matchedPartTrackDeltaEta.Fill(deltaEta)
matchedPartTrackDeltaPhi.Fill(deltaPhi)
matchedPartTrackDeltaR.Fill(deltaR)
matchedPartTrackDeltaMom.Fill(deltaMom)
matchedPartEta.Fill(trueEta)
matchedPartPhi.Fill(truePhi)
matchedPartMom.Fill(trueMom.Mag())
matchedPartPEta.Fill(trueMom.Mag(), trueEta)
matchedPartPhiEta.Fill(truePhi, trueEta)
for x in range (0, len(trackMomX[i])): # Loop over all charged particles, thrown or not
CPartMom = ROOT.TVector3(trackMomX[i][x], trackMomY[i][x], trackMomZ[i][x])
CPartEta = CPartMom.PseudoRapidity()
CPartPhi = CPartMom.Phi()
ChargedEta.Fill(CPartEta)
ChargedPhi.Fill(CPartPhi)
ChargedP.Fill(CPartMom.Mag())
# Write output histograms to file below
partEta.Write()
matchedPartEta.Write()
partMom.Write()
matchedPartMom.Write()
partPhi.Write()
matchedPartPhi.Write()
partPEta.Write()
matchedPartPEta.Write()
partPhiEta.Write()
matchedPartPhiEta.Write()
matchedPartTrackDeltaEta.Write()
matchedPartTrackDeltaPhi.Write()
matchedPartTrackDeltaR.Write()
matchedPartTrackDeltaMom.Write()
ChargedEta.Write()
ChargedPhi.Write()
ChargedP.Write()
TrackEff_Eta.Divide(matchedPartEta, partEta, 1, 1, "b")
TrackEffMom.Divide(matchedPartMom, partMom, 1, 1, "b")
TrackEffPhi.Divide(matchedPartPhi, partPhi, 1, 1, "b")
TrackEff_PEta.Divide(matchedPartPEta, partPEta, 1, 1, "b")
TrackEff_PhiEta.Divide(matchedPartPhiEta, partPhiEta, 1, 1, "b")
TrackEff_Eta.Write()
TrackEffMom.Write()
TrackEffPhi.Write()
TrackEff_PEta.Write()
TrackEff_PhiEta.Write()
# Close files
ofile.Close()Insert your input file path and execute as the example code above. ### ResolutionAnalysis.py
Create a file called ResolutionAnalysis.py and copy in
the code below to get started on the resolution analysis exercise. Note
that you will need to correctly specify your input file path in the
variable infile.
PYTHON
#! /usr/bin/python
#Import relevant packages
import ROOT, math, array
from ROOT import TH1F, TH2F, TMath, TTree, TVector3, TVector2
import uproot as up
#Define and open files
infile="PATH_TO_INPUT_FILE"
ofile=ROOT.TFile.Open("ResolutionAnalysis_OutPy.root", "RECREATE")
# Open input file and define branches we want to look at with uproot
events_tree = up.open(infile)["events"]
# Get particle information
partGenStat = events_tree["MCParticles.generatorStatus"].array()
partMomX = events_tree["MCParticles.momentum.x"].array()
partMomY = events_tree["MCParticles.momentum.y"].array()
partMomZ = events_tree["MCParticles.momentum.z"].array()
partPdg = events_tree["MCParticles.PDG"].array()
# Get reconstructed track information
trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
# Get assocations between MCParticles and ReconstructedChargedParticles
recoAssoc = events_tree["_ReconstructedChargedParticleAssociations_rec.index"].array()
simuAssoc = events_tree["_ReconstructedChargedParticleAssociations_sim.index"].array()
# Define histograms below
trackMomentumRes = ROOT.TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2)
matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
# Add main analysis loop(s) below
for i in range(0, len(partGenStat)): # Loop over all events
for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
if partGenStat[i][j] == 1: # Select stable particles
pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
trueEta = trueMom.PseudoRapidity()
truePhi = trueMom.Phi()
for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
if (simuAssoc[i][k] == j):
recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
deltaEta = trueEta - recMom.PseudoRapidity()
deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag()
matchedPartTrackDeltaEta.Fill(deltaEta)
matchedPartTrackDeltaPhi.Fill(deltaPhi)
matchedPartTrackDeltaR.Fill(deltaR)
matchedPartTrackDeltaMom.Fill(deltaMom)
trackMomentumRes.Fill(momRes)
# Write output histograms to file below
trackMomentumRes.Write()
matchedPartTrackDeltaEta.Write()
matchedPartTrackDeltaPhi.Write()
matchedPartTrackDeltaR.Write()
matchedPartTrackDeltaMom.Write()
# Close files
ofile.Close()A “solution” version of the script for the exercise is included below -
PYTHON
#! /usr/bin/python
#Import relevant packages
import ROOT, math, array
from ROOT import TCanvas, TColor, TGaxis, TH1F, TH2F, TPad, TStyle, gStyle, gPad, TGaxis, TLine, TMath, TPaveText, TTree, TVector3, TVector2
import uproot as up
#Define and open files
infile="PATH_TO_FILE"
ofile=ROOT.TFile.Open("ResolutionAnalysis_Exercise_OutPy.root", "RECREATE")
# Open input file and define branches we want to look at with uproot
events_tree = up.open(infile)["events"]
# Get particle information
partGenStat = events_tree["MCParticles.generatorStatus"].array()
partMomX = events_tree["MCParticles.momentum.x"].array()
partMomY = events_tree["MCParticles.momentum.y"].array()
partMomZ = events_tree["MCParticles.momentum.z"].array()
partPdg = events_tree["MCParticles.PDG"].array()
# Get reconstructed track information
trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
# Get assocations between MCParticles and ReconstructedChargedParticles
recoAssoc = events_tree["_ReconstructedChargedParticleAssociations_rec.index"].array()
simuAssoc = events_tree.["_ReconstructedChargedParticleAssociations_sim.index"].array()
# Define histograms below
trackMomentumRes = ROOT.TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2)
trackMomResP = ROOT.TH2D("trackMomResP", "Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
trackMomResEta = ROOT.TH2D("trackMomResEta", "Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
trackMomentumRes_e = ROOT.TH1D("trackMomentumRes_e","e^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
trackMomResP_e = ROOT.TH2D("trackMomResP_e", "e^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 25);
trackMomResEta_e = ROOT.TH2D("trackMomResEta_e", "e^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
trackMomentumRes_mu = ROOT.TH1D("trackMomentumRes_mu","#mu^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
trackMomResP_mu = ROOT.TH2D("trackMomResP_mu", "#mu^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 25);
trackMomResEta_mu = ROOT.TH2D("trackMomResEta_mu", "#mu^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
trackMomentumRes_pi = ROOT.TH1D("trackMomentumRes_pi","#pi^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
trackMomResP_pi = ROOT.TH2D("trackMomResP_pi", "#pi^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
trackMomResEta_pi = ROOT.TH2D("trackMomResEta_pi", "#pi^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
trackMomentumRes_K = ROOT.TH1D("trackMomentumRes_K","K^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
trackMomResP_K = ROOT.TH2D("trackMomResP_K", "K^{#pm} Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
trackMomResEta_K = ROOT.TH2D("trackMomResEta_K", "K^{#pm} Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
trackMomentumRes_p = ROOT.TH1D("trackMomentumRes_p","p Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
trackMomResP_p = ROOT.TH2D("trackMomResP_p", "p Track Momentum Resolution vs P; (P_{rec} - P_{MC})/P_{MC}; P_{MC}(GeV/c)", 400, -2, 2, 150, 0, 150);
trackMomResEta_p = ROOT.TH2D("trackMomResEta_p", "p Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
# Add main analysis loop(s) below
for i in range(0, len(partGenStat)): # Loop over all events
for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
if partGenStat[i][j] == 1: # Select stable particles
pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
trueEta = trueMom.PseudoRapidity()
truePhi = trueMom.Phi()
for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
if (simuAssoc[i][k] == j):
recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
deltaEta = trueEta - recMom.PseudoRapidity()
deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag()
trackMomentumRes.Fill(momRes)
trackMomResP.Fill(momRes, trueMom.Mag())
trackMomResEta.Fill(momRes, trueEta)
if( pdg == 11):
trackMomentumRes_e.Fill(momRes)
trackMomResP_e.Fill(momRes, trueMom.Mag())
trackMomResEta_e.Fill(momRes, trueEta)
elif( pdg == 13):
trackMomentumRes_mu.Fill(momRes)
trackMomResP_mu.Fill(momRes, trueMom.Mag())
trackMomResEta_mu.Fill(momRes, trueEta)
elif( pdg == 211):
trackMomentumRes_pi.Fill(momRes)
trackMomResP_pi.Fill(momRes, trueMom.Mag())
trackMomResEta_pi.Fill(momRes, trueEta)
elif( pdg == 321):
trackMomentumRes_K.Fill(momRes)
trackMomResP_K.Fill(momRes, trueMom.Mag())
trackMomResEta_K.Fill(momRes, trueEta)
elif( pdg == 2212):
trackMomentumRes_p.Fill(momRes)
trackMomResP_p.Fill(momRes, trueMom.Mag())
trackMomResEta_p.Fill(momRes, trueEta)
matchedPartTrackDeltaEta.Fill(deltaEta)
matchedPartTrackDeltaPhi.Fill(deltaPhi)
matchedPartTrackDeltaR.Fill(deltaR)
matchedPartTrackDeltaMom.Fill(deltaMom)
# Write output histograms to file below
trackMomentumRes.Write()
trackMomResP.Write()
trackMomResEta.Write()
trackMomentumRes_e.Write()
trackMomResP_e.Write()
trackMomResEta_e.Write()
trackMomentumRes_mu.Write()
trackMomResP_mu.Write()
trackMomResEta_mu.Write()
trackMomentumRes_pi.Write()
trackMomResP_pi.Write()
trackMomResEta_pi.Write()
trackMomentumRes_K.Write()
trackMomResP_K.Write()
trackMomResEta_K.Write()
trackMomentumRes_p.Write()
trackMomResP_p.Write()
trackMomResEta_p.Write()
matchedPartTrackDeltaEta.Write()
matchedPartTrackDeltaPhi.Write()
matchedPartTrackDeltaR.Write()
matchedPartTrackDeltaMom.Write()
# Close files
ofile.Close()Insert your input file path and execute as the example code above. ## RDataFrames Example
Note that only the initial stage of the efficiency example is presented here in RDF format. This example was kindly created by Simon.
EfficiencyAnalysisRDF.C
Create a file called EfficiencyAnalysisRDF.C and paste
the code below in. Remember to change the file path.
Execute this script via -
root -l -q EfficiencyAnalysisRDF.C++. Do this within
eic-shell or somewhere else with the correct EDM4hep/EDM4eic libraries
installed.
CPP
#include <edm4hep/utils/vector_utils.h>
#include <edm4hep/MCParticle.h>
#include <edm4eic/ReconstructedParticle.h>
#include <ROOT/RDataFrame.hxx>
#include <ROOT/RVec.hxx>
#include <TFile.h>
// Define aliases for the data types
using MCP = edm4hep::MCParticleData;
using RecoP = edm4eic::ReconstructedParticleData;
// Define function to vectorize the edm4hep::utils methods
template <typename T>
auto getEta = [](ROOT::VecOps::RVec<T> momenta) {
return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::eta(p.momentum); });
};
template <typename T>
auto getPhi = [](ROOT::VecOps::RVec<T> momenta) {
return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::angleAzimuthal(p.momentum); });
};
// Define the function to perform the efficiency analysis
void EfficiencyAnalysisRDF(TString infile="PATH_TO_FILE"){
// Set up input file
ROOT::RDataFrame df("events", infile);
// Define new dataframe node with additional columns
auto df1 = df.Define("statusFilter", "MCParticles.generatorStatus == 1" )
.Define("absPDG", "abs(MCParticles.PDG)" )
.Define("pdgFilter", "absPDG == 11 || absPDG == 13 || absPDG == 211 || absPDG == 321 || absPDG == 2212")
.Define("particleFilter","statusFilter && pdgFilter" )
.Define("filtMCParts", "MCParticles[particleFilter]" )
.Define("assoFilter", "Take(particleFilter,ReconstructedChargedParticleAssociations_sim.index)") // Incase any of the associated particles happen to not be charged
.Define("assoMCParts", "Take(MCParticles,ReconstructedChargedParticleAssociations)sim.index)[assoFilter]")
.Define("assoRecParts", "Take(ReconstructedChargedParticles,ReconstructedChargedParticleAssociations._rec.index)[assoFilter]")
.Define("filtMCEta", getEta<MCP> , {"filtMCParts"} )
.Define("filtMCPhi", getPhi<MCP> , {"filtMCParts"} )
.Define("accoMCEta", getEta<MCP> , {"assoMCParts"} )
.Define("accoMCPhi", getPhi<MCP> , {"assoMCParts"} )
.Define("assoRecEta", getEta<RecoP> , {"assoRecParts"})
.Define("assoRecPhi", getPhi<RecoP> , {"assoRecParts"})
.Define("deltaR", "ROOT::VecOps::DeltaR(assoRecEta, accoMCEta, assoRecPhi, accoMCPhi)");
// Define histograms
auto partEta = df1.Histo1D({"partEta","Eta of Thrown Charged Particles;Eta",100,-5.,5.},"filtMCEta");
auto matchedPartEta = df1.Histo1D({"matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track",100,-5.,5.},"accoMCEta");
auto matchedPartTrackDeltaR = df1.Histo1D({"matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charged Particle",5000,0.,5.},"deltaR");
// Write histograms to file
TFile *ofile = TFile::Open("EfficiencyAnalysis_Out_RDF.root","RECREATE");
// Booked Define and Histo1D lazy actions are only performed here
partEta->Write();
matchedPartEta->Write();
matchedPartTrackDeltaR->Write();
ofile->Close(); // Close output file
}A “solution” using RDataFrames is included below,
CPP
#include <edm4hep/utils/vector_utils.h>
#include <edm4hep/MCParticle.h>
#include <edm4eic/ReconstructedParticle.h>
#include <ROOT/RDataFrame.hxx>
#include <ROOT/RVec.hxx>
#include <TFile.h>
// Define aliases for the data types
using MCP = edm4hep::MCParticleData;
using RecoP = edm4eic::ReconstructedParticleData;
// Define function to vectorize the edm4hep::utils methods
template <typename T>
auto getEta = [](ROOT::VecOps::RVec<T> momenta) {
return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::eta(p.momentum); });
};
template <typename T>
auto getPhi = [](ROOT::VecOps::RVec<T> momenta) {
return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::angleAzimuthal(p.momentum); });
};
template <typename T>
auto getP = [](ROOT::VecOps::RVec<T> momenta) {
//return ROOT::VecOps::Map(momenta, [](const T& p) { return (p.momentum); }); // This is a vector3f
return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::magnitude(p.momentum); }); // This is a the magnitude of that vector3f
};
// Define the function to perform the efficiency analysis
void EfficiencyAnalysisRDF_Exercise(TString infile="PATH_TO_INPUT_FILE"){
// Set up input file
ROOT::RDataFrame df("events", infile);
// Define new dataframe node with additional columns
auto df1 = df.Define("statusFilter", "MCParticles.generatorStatus == 1" )
.Define("absPDG", "abs(MCParticles.PDG)" )
.Define("pdgFilter", "absPDG == 11 || absPDG == 13 || absPDG == 211 || absPDG == 321 || absPDG == 2212")
.Define("particleFilter","statusFilter && pdgFilter" )
.Define("filtMCParts", "MCParticles[particleFilter]" )
.Define("assoFilter", "Take(particleFilter,_ReconstructedChargedParticleAssociations_sim.index)") // In case any of the associated particles happen to not be charged
.Define("assoMCParts", "Take(MCParticles,_ReconstructedChargedParticleAssociations_sim.index)[assoFilter]")
.Define("assoRecParts", "Take(ReconstructedChargedParticles,_ReconstructedChargedParticleAssociations_rec.index)[assoFilter]")
.Define("filtMCEta", getEta<MCP> , {"filtMCParts"} )
.Define("filtMCPhi", getPhi<MCP> , {"filtMCParts"} )
.Define("filtMCp", getP<MCP> , {"filtMCParts"} )
.Define("assoMCEta", getEta<MCP> , {"assoMCParts"} )
.Define("assoMCPhi", getPhi<MCP> , {"assoMCParts"} )
.Define("assoMCp", getP<MCP> , {"assoMCParts"} )
.Define("assoRecEta", getEta<RecoP> , {"assoRecParts"})
.Define("assoRecPhi", getPhi<RecoP> , {"assoRecParts"})
.Define("assoRecp", getP<RecoP> , {"assoRecParts"})
.Define("deltaEta", "assoMCEta - assoRecEta" )
.Define("deltaPhi", "ROOT::VecOps::DeltaPhi(assoRecPhi, assoMCPhi)")
.Define("deltaR", "ROOT::VecOps::DeltaR(assoRecEta, assoMCEta, assoRecPhi, assoMCPhi)")
.Define("deltaMom", "assoMCp - assoRecp")
.Define("recoEta", getEta<RecoP>, {"ReconstructedChargedParticles"})
.Define("recoPhi", getPhi<RecoP>, {"ReconstructedChargedParticles"})
.Define("recoP", getP<RecoP>, {"ReconstructedChargedParticles"});
// Define histograms. We create a histogram with the usual naming/titles/bins/range etc, then specify how to fill the histogram based upon things we have defined for our dataframe
auto partEta = df1.Histo1D({"partEta","#eta of Thrown Charged Particles;#eta",120,-6.,6.},"filtMCEta");
auto matchedPartEta = df1.Histo1D({"matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track;#eta",120,-6.,6.},"assoMCEta");
auto partMom = df1.Histo1D({"partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150}, "filtMCp");
auto matchedPartMom = df1.Histo1D({"matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150}, "assoMCp");
auto partPhi = df1.Histo1D({"partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2},"filtMCPhi");
auto matchedPartPhi = df1.Histo1D({"matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2}, "assoMCPhi");
auto partPEta = df1.Histo2D({"partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6}, "filtMCp", "filtMCEta");
auto matchedPartPEta = df1.Histo2D({"matchedPartPEta", "P vs #eta of Thrown Charged Particles, with matching track; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6}, "assoMCp", "assoMCEta");
auto partPhiEta = df1.Histo2D({"partPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6}, "filtMCPhi", "filtMCEta");
auto matchedPartPhiEta = df1.Histo2D({"matchedPartPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6}, "assoMCPhi", "assoMCEta");
auto matchedPartTrackDeltaEta = df1.Histo1D({"matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25}, "deltaEta");
auto matchedPartTrackDeltaPhi = df1.Histo1D({"matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2}, "deltaPhi");
auto matchedPartTrackDeltaR = df1.Histo1D({"matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle;#Delta R",300,0.,0.3}, "deltaR");
auto matchedPartTrackDeltaMom = df1.Histo1D({"matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10}, "deltaMom");
// Define some histograms for our efficiencies - Done "old school" root style - Maybe the division can be done direct from a DF?
TH1D *TrackEff_Eta = new TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6);
TH1D *TrackEff_Mom = new TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150);
TH1D *TrackEff_Phi = new TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2);
// 2D Efficiencies
TH2D* TrackEff_PEta = new TH2D("TrackEff_PEta", "Tracking efficiency as fn of P and #eta; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6);
TH2D* TrackEff_PhiEta = new TH2D("TrackEff_PhiEta", "Tracking efficiency as fn of #phi and #eta; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6);
auto ChargedEta = df1.Histo1D({"ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6}, "recoEta");
auto ChargedPhi = df1.Histo1D({"ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2}, "recoPhi");
auto ChargedP = df1.Histo1D({"ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150}, "recoP");
// Write histograms to file
TFile *ofile = TFile::Open("EfficiencyAnalysis_Exercise_Out_RDF.root","RECREATE");
// Booked Define and Histo1D lazy actions are only performed here
partEta->Write();
matchedPartEta->Write();
partPhi->Write();
matchedPartPhi->Write();
partMom->Write();
matchedPartMom->Write();
partPEta->Write();
matchedPartPEta->Write();
partPhiEta->Write();
matchedPartPhiEta->Write();
matchedPartTrackDeltaEta->Write();
matchedPartTrackDeltaPhi->Write();
matchedPartTrackDeltaR->Write();
matchedPartTrackDeltaMom->Write();
// Create efficiency histograms by dividing appropriately. Note we must actually get the pointer explicitly.
TrackEff_Eta->Divide(matchedPartEta.GetPtr(), partEta.GetPtr(), 1, 1, "b");
TrackEff_Mom->Divide(matchedPartMom.GetPtr(), partMom.GetPtr(), 1, 1, "b");
TrackEff_Phi->Divide(matchedPartPhi.GetPtr(), partPhi.GetPtr(), 1, 1, "b");
TrackEff_PEta->Divide(matchedPartPEta.GetPtr(), partPEta.GetPtr(), 1, 1, "b");
TrackEff_PhiEta->Divide(matchedPartPhiEta.GetPtr(), partPhiEta.GetPtr(), 1, 1, "b");
TrackEff_Eta->Write();
TrackEff_Mom->Write();
TrackEff_Phi->Write();
TrackEff_PEta->Write();
TrackEff_PhiEta->Write();
ChargedEta->Write();
ChargedPhi->Write();
ChargedP->Write();
ofile->Close(); // Close output file
}