tutorial-analysis/_extras/exercise_scripts.md

0001 ---
0002 title: "Exercise Scripts"
0003 ---
0004
0005 Included below is a selection of scripts for the exercises in part 3 of this tutorial. Shortly after the tutorial, I will also include "complete" examples for future reference.
0006
0007 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.
0008
0009 ## ROOT TTreeReader Scripts
0010
0011 ### EfficiencyAnalysis.C
0012
0013 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 specifiy your input file path in the first line.
0014
0015 ```c++
0016 void EfficiencyAnalysis(TString infile="PATH_TO_INPUT_FILE"){
0017   // Set output file for the histograms
0018   TFile *ofile = TFile::Open("EfficiencyAnalysis_Out.root","RECREATE");
0019
0020   // Set up input file chain
0021   TChain *mychain = new TChain("events");
0022   mychain->Add(infile);
0023
0024   // Initialize reader
0025   TTreeReader tree_reader(mychain);
0026
0027   // Get Particle Information
0028   TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
0029   TTreeReaderArray<float> partMomX(tree_reader, "MCParticles.momentum.x");
0030   TTreeReaderArray<float> partMomY(tree_reader, "MCParticles.momentum.y");
0031   TTreeReaderArray<float> partMomZ(tree_reader, "MCParticles.momentum.z");
0032   TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
0033
0034   // Get Reconstructed Track Information
0035   TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
0036   TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
0037   TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
0038
0039   // Get Associations Between MCParticles and ReconstructedChargedParticles
0040   TTreeReaderArray<unsigned int> recoAssoc(tree_reader, "ReconstructedChargedParticleAssociations.recID");
0041   TTreeReaderArray<unsigned int> simuAssoc(tree_reader, "ReconstructedChargedParticleAssociations.simID");
0042
0043   // Define Histograms
0044   TH1D *partEta = new TH1D("partEta","Eta of Thrown Charged Particles;Eta",100,-5.,5.);
0045   TH1D *matchedPartEta = new TH1D("matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track",100,-5.,5.);
0046   TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charged Particle",5000,0.,5.);
0047
0048   while(tree_reader.Next()) { // Loop over events
0049     for(unsigned int i=0; i<partGenStat.GetSize(); i++){ // Loop over thrown particles
0050         if(partGenStat[i] == 1){ // Select stable thrown particles
0051             int pdg = TMath::Abs(partPdg[i]);
0052             if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212){ // Look at charged particles (electrons, muons, pions, kaons, protons)
0053                 TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
0054
0055                 float trueEta = trueMom.PseudoRapidity();
0056                 float truePhi = trueMom.Phi();
0057
0058                 partEta->Fill(trueEta);
0059
0060                 for(unsigned int j=0; j<simuAssoc.GetSize(); j++){ // Loop over associations to find matching ReconstructedChargedParticle
0061                     if(simuAssoc[j] == i){ // Find association index matching the index of the thrown particle we are looking at
0062                         TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
0063
0064                         // Check the distance between the thrown and reconstructed particle
0065                         float deltaEta = trueEta - recMom.PseudoRapidity();
0066                         float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
0067                         float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
0068
0069                         matchedPartTrackDeltaR->Fill(deltaR);
0070
0071                         matchedPartEta->Fill(trueEta); // Plot the thrown eta if a matched ReconstructedChargedParticle was found
0072                     }
0073                 } // End loop over associations
0074             } // End PDG check
0075         } // End stable particles condition
0076     } // End loop over thrown particles
0077   } // End loop over events
0078   ofile->Write(); // Write histograms to file
0079   ofile->Close(); // Close output file
0080 }
0081 ```
0082 <!--
0083 A "solution" version of the script for the exercise is included below -
0084
0085 ```c++
0086 void EfficiencyAnalysis_Exercise(TString infile="PATH_TO_FILE"){
0087
0088   // Set output file for the histograms
0089   TFile *ofile = TFile::Open("EfficiencyAnalysis_Exercise_Out.root","RECREATE");
0090
0091   // Analysis code will go here
0092   // Set up input file chain
0093   TChain *mychain = new TChain("events");
0094   mychain->Add(infile);
0095
0096   // Initialize reader
0097   TTreeReader tree_reader(mychain);
0098
0099   // Get Particle Information
0100   TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
0101   TTreeReaderArray<float> partMomX(tree_reader, "MCParticles.momentum.x");
0102   TTreeReaderArray<float> partMomY(tree_reader, "MCParticles.momentum.y");
0103   TTreeReaderArray<float> partMomZ(tree_reader, "MCParticles.momentum.z");
0104   TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
0105
0106   // Get Reconstructed Track Information
0107   TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
0108   TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
0109   TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
0110
0111   // Get Associations Between MCParticles and ReconstructedChargedParticles
0112   TTreeReaderArray<unsigned int> recoAssoc(tree_reader, "ReconstructedChargedParticleAssociations.recID");
0113   TTreeReaderArray<unsigned int> simuAssoc(tree_reader, "ReconstructedChargedParticleAssociations.simID");
0114
0115   // Define Histograms
0116   TH1D *partEta = new TH1D("partEta","#eta of Thrown Charged Particles; #eta", 120, -6, 6);
0117   TH1D *matchedPartEta = new TH1D("matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track; #eta", 120, -6, 6);
0118   TH1D* partMom = new TH1D("partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150);
0119   TH1D* matchedPartMom = new TH1D("matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150);
0120   TH1D* partPhi = new TH1D("partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2);
0121   TH1D* matchedPartPhi = new TH1D("matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2);
0122
0123   TH2D* partPEta = new TH2D("partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6);
0124   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);
0125   TH2D* partPhiEta = new TH2D("partPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6);
0126   TH2D* matchedPartPhiEta = new TH2D("matchedPartPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6);
0127
0128   TH1D *matchedPartTrackDeltaEta = new TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25);
0129   TH1D *matchedPartTrackDeltaPhi = new TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2);
0130   TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3);
0131   TH1D *matchedPartTrackDeltaMom = new TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10);
0132
0133   // Define some histograms for our efficiencies
0134   TH1D *TrackEff_Eta = new TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6);
0135   TH1D *TrackEff_Mom = new TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150);
0136   TH1D *TrackEff_Phi = new TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2);
0137
0138   // 2D Efficiencies
0139   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);
0140   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);
0141
0142   // All charged particle histos
0143   TH1D *ChargedEta = new TH1D("ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6);
0144   TH1D *ChargedPhi = new TH1D("ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2);
0145   TH1D *ChargedP = new TH1D("ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150);
0146
0147   while(tree_reader.Next()) { // Loop over events
0148
0149     for(unsigned int i=0; i<partGenStat.GetSize(); i++) // Loop over thrown particles
0150       {
0151         if(partGenStat[i] == 1) // Select stable thrown particles
0152           {
0153             int pdg = TMath::Abs(partPdg[i]);
0154
0155             if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212) // Look at charged particles (electrons, muons, pions, kaons, protons)
0156               {
0157                 TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
0158
0159                 float trueEta = trueMom.PseudoRapidity();
0160                 float truePhi = trueMom.Phi();
0161
0162                 partEta->Fill(trueEta);
0163                 partPhi->Fill(truePhi);
0164                 partMom->Fill(trueMom.Mag());
0165                 partPEta->Fill(trueMom.Mag(), trueEta);
0166                 partPhiEta->Fill(truePhi, trueEta);
0167
0168                 // Loop over associations to find matching ReconstructedChargedParticle
0169                 for(unsigned int j=0; j<simuAssoc.GetSize(); j++)
0170                   {
0171                     if(simuAssoc[j] == i) // Find association index matching the index of the thrown particle we are looking at
0172                       {
0173                         TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
0174
0175                         // Check the distance between the thrown and reconstructed particle
0176                         float deltaEta = trueEta - recMom.PseudoRapidity();
0177                         float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
0178                         float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
0179                         float deltaMom = ((trueMom.Mag()) - (recMom.Mag()));
0180
0181                         matchedPartTrackDeltaEta->Fill(deltaEta);
0182                         matchedPartTrackDeltaPhi->Fill(deltaPhi);
0183                         matchedPartTrackDeltaR->Fill(deltaR);
0184                         matchedPartTrackDeltaMom->Fill(deltaMom);
0185
0186                         matchedPartEta->Fill(trueEta); // Plot the thrown eta if a matched ReconstructedChargedParticle was found
0187                         matchedPartPhi->Fill(truePhi);
0188                         matchedPartMom->Fill(trueMom.Mag());
0189
0190                         matchedPartPEta->Fill(trueMom.Mag(), trueEta);
0191                         matchedPartPhiEta->Fill(truePhi, trueEta);
0192
0193                       }
0194                   }// End loop over associations
0195               } // End PDG check
0196           } // End stable particles condition
0197       } // End loop over thrown particles
0198     // Loop over all charged particles and fill some histograms of kinematics quantities
0199     for(unsigned int k=0; k<trackMomX.GetSize(); k++){ // Loop over all charged particles, thrown or not
0200
0201       TVector3 CPartMom(trackMomX[k], trackMomY[k], trackMomZ[k]);
0202
0203       float CPartEta = CPartMom.PseudoRapidity();
0204       float CPartPhi = CPartMom.Phi();
0205
0206       ChargedEta->Fill(CPartEta);
0207       ChargedPhi->Fill(CPartPhi);
0208       ChargedP->Fill(CPartMom.Mag());
0209
0210     } // End loop over all charged particles
0211   } // End loop over events
0212
0213   // Take the ratio of the histograms above to get our efficiency plots
0214   TrackEff_Eta->Divide(matchedPartEta, partEta, 1, 1, "b");
0215   TrackEff_Mom->Divide(matchedPartMom, partMom, 1, 1, "b");
0216   TrackEff_Phi->Divide(matchedPartPhi, partPhi, 1, 1, "b");
0217   TrackEff_PEta->Divide(matchedPartPEta, partPEta, 1, 1, "b");
0218   TrackEff_PhiEta->Divide(matchedPartPhiEta, partPhiEta, 1, 1, "b");
0219
0220   ofile->Write(); // Write histograms to file
0221   ofile->Close(); // Close output file
0222 }
0223 ```
0224 Insert your input file path and execute as the example code above.
0225 -->
0226
0227 ### ResolutionAnalysis.C
0228
0229 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 specifiy your input file path in the first line.
0230
0231 ```c++
0232 void ResolutionAnalysis(TString infile="PATH_TO_INPUT_FILE"){
0233   // Set output file for the histograms
0234   TFile *ofile = TFile::Open("ResolutionAnalysis_Out.root","RECREATE");
0235
0236   // Analysis code will go here
0237   // Set up input file chain
0238   TChain *mychain = new TChain("events");
0239   mychain->Add(infile);
0240
0241   // Initialize reader
0242   TTreeReader tree_reader(mychain);
0243
0244   // Get Particle Information
0245   TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
0246   TTreeReaderArray<float> partMomX(tree_reader, "MCParticles.momentum.x");
0247   TTreeReaderArray<float> partMomY(tree_reader, "MCParticles.momentum.y");
0248   TTreeReaderArray<float> partMomZ(tree_reader, "MCParticles.momentum.z");
0249   TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
0250
0251   // Get Reconstructed Track Information
0252   TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
0253   TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
0254   TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
0255
0256   // Get Associations Between MCParticles and ReconstructedChargedParticles
0257   TTreeReaderArray<unsigned int> recoAssoc(tree_reader, "ReconstructedChargedParticleAssociations.recID");
0258   TTreeReaderArray<unsigned int> simuAssoc(tree_reader, "ReconstructedChargedParticleAssociations.simID");
0259
0260   // Define Histograms
0261   TH1D *trackMomentumRes = new TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2);
0262
0263   TH1D *matchedPartTrackDeltaEta = new TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25);
0264   TH1D *matchedPartTrackDeltaPhi = new TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2);
0265   TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3);
0266   TH1D *matchedPartTrackDeltaMom = new TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10);
0267   while(tree_reader.Next()) { // Loop over events
0268     for(unsigned int i=0; i<partGenStat.GetSize(); i++){ // Loop over thrown particles
0269         if(partGenStat[i] == 1){ // Select stable thrown particles
0270             int pdg = TMath::Abs(partPdg[i]);
0271             if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212){ // Look at charged particles (electrons, muons, pions, kaons, protons)
0272                 TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
0273
0274                 float trueEta = trueMom.PseudoRapidity();
0275                 float truePhi = trueMom.Phi();
0276
0277                 for(unsigned int j=0; j<simuAssoc.GetSize(); j++){ // Loop over associations to find matching ReconstructedChargedParticle
0278                     if(simuAssoc[j] == i){ // Find association index matching the index of the thrown particle we are looking at
0279                         TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
0280
0281                         // Check the distance between the thrown and reconstructed particle
0282                         float deltaEta = trueEta - recMom.PseudoRapidity();
0283                         float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
0284                         float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
0285                         float deltaMom = ((trueMom.Mag()) - (recMom.Mag()));
0286                         double momRes = (recMom.Mag()- trueMom.Mag())/trueMom.Mag();
0287
0288                         trackMomentumRes->Fill(momRes);
0289
0290                         matchedPartTrackDeltaEta->Fill(deltaEta);
0291                         matchedPartTrackDeltaPhi->Fill(deltaPhi);
0292                         matchedPartTrackDeltaR->Fill(deltaR);
0293                         matchedPartTrackDeltaMom->Fill(deltaMom);
0294                     }
0295                 } // End loop over associations
0296             } // End PDG check
0297         } // End stable particles condition
0298     } // End loop over thrown particles
0299   } // End loop over events
0300   ofile->Write(); // Write histograms to file
0301   ofile->Close(); // Close output file
0302 }
0303 ```
0304 <!--
0305 A "solution" version of the script for the exercise is included below -
0306
0307 ```c++
0308 void ResolutionAnalysis_Exercise(TString infile="PATH_TO_FILE"){
0309   // Set output file for the histograms
0310   TFile *ofile = TFile::Open("ResolutionAnalysis_Exercise_Out.root","RECREATE");
0311
0312   // Analysis code will go here
0313   // Set up input file chain
0314   TChain *mychain = new TChain("events");
0315   mychain->Add(infile);
0316
0317   // Initialize reader
0318   TTreeReader tree_reader(mychain);
0319
0320   // Get Particle Information
0321   TTreeReaderArray<int> partGenStat(tree_reader, "MCParticles.generatorStatus");
0322   TTreeReaderArray<float> partMomX(tree_reader, "MCParticles.momentum.x");
0323   TTreeReaderArray<float> partMomY(tree_reader, "MCParticles.momentum.y");
0324   TTreeReaderArray<float> partMomZ(tree_reader, "MCParticles.momentum.z");
0325   TTreeReaderArray<int> partPdg(tree_reader, "MCParticles.PDG");
0326
0327   // Get Reconstructed Track Information
0328   TTreeReaderArray<float> trackMomX(tree_reader, "ReconstructedChargedParticles.momentum.x");
0329   TTreeReaderArray<float> trackMomY(tree_reader, "ReconstructedChargedParticles.momentum.y");
0330   TTreeReaderArray<float> trackMomZ(tree_reader, "ReconstructedChargedParticles.momentum.z");
0331
0332   // Get Associations Between MCParticles and ReconstructedChargedParticles
0333   TTreeReaderArray<unsigned int> recoAssoc(tree_reader, "ReconstructedChargedParticleAssociations.recID");
0334   TTreeReaderArray<unsigned int> simuAssoc(tree_reader, "ReconstructedChargedParticleAssociations.simID");
0335
0336   // Define Histograms
0337   TH1D *trackMomentumRes = new TH1D("trackMomentumRes","Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0338   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);
0339   TH2D* trackMomResEta = new TH2D("trackMomResEta", "Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
0340
0341   TH1D *trackMomentumRes_e = new TH1D("trackMomentumRes_e","e^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0342   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);
0343   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);
0344
0345   TH1D *trackMomentumRes_mu = new TH1D("trackMomentumRes_mu","#mu^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0346   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);
0347   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);
0348
0349   TH1D *trackMomentumRes_pi = new TH1D("trackMomentumRes_pi","#pi^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0350   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);
0351   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);
0352
0353   TH1D *trackMomentumRes_K = new TH1D("trackMomentumRes_K","K^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0354   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);
0355   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);
0356
0357   TH1D *trackMomentumRes_p = new TH1D("trackMomentumRes_p","p Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0358   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);
0359   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);
0360
0361   TH1D *matchedPartTrackDeltaEta = new TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25);
0362   TH1D *matchedPartTrackDeltaPhi = new TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2);
0363   TH1D *matchedPartTrackDeltaR = new TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3);
0364   TH1D *matchedPartTrackDeltaMom = new TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10);
0365
0366   while(tree_reader.Next()) { // Loop over events
0367
0368     for(unsigned int i=0; i<partGenStat.GetSize(); i++) // Loop over thrown particles
0369       {
0370         if(partGenStat[i] == 1) // Select stable thrown particles
0371           {
0372             int pdg = TMath::Abs(partPdg[i]);
0373
0374             if(pdg == 11 || pdg == 13 || pdg == 211 || pdg == 321 || pdg == 2212) // Look at charged particles (electrons, muons, pions, kaons, protons)
0375               {
0376                 TVector3 trueMom(partMomX[i],partMomY[i],partMomZ[i]);
0377
0378                 float trueEta = trueMom.PseudoRapidity();
0379                 float truePhi = trueMom.Phi();
0380
0381                 // Loop over associations to find matching ReconstructedChargedParticle
0382                 for(unsigned int j=0; j<simuAssoc.GetSize(); j++)
0383                   {
0384                     if(simuAssoc[j] == i) // Find association index matching the index of the thrown particle we are looking at
0385                       {
0386                         TVector3 recMom(trackMomX[recoAssoc[j]],trackMomY[recoAssoc[j]],trackMomZ[recoAssoc[j]]); // recoAssoc[j] is the index of the matched ReconstructedChargedParticle
0387
0388                         // Check the distance between the thrown and reconstructed particle
0389                         float deltaEta = trueEta - recMom.PseudoRapidity();
0390                         float deltaPhi = TVector2::Phi_mpi_pi(truePhi - recMom.Phi());
0391                         float deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
0392                         float deltaMom = ((trueMom.Mag()) - (recMom.Mag()));
0393
0394                         double momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag();
0395
0396                         trackMomentumRes->Fill(momRes); // Could also multiply by 100 and express as a percentage instead
0397                         trackMomResP->Fill(momRes, trueMom.Mag());
0398                         trackMomResEta->Fill(momRes, trueEta);
0399
0400                         if( pdg == 11){
0401                           trackMomentumRes_e->Fill(momRes);
0402                           trackMomResP_e->Fill(momRes, trueMom.Mag());
0403                           trackMomResEta_e->Fill(momRes, trueEta);
0404                         }
0405                         else if( pdg == 13){
0406                           trackMomentumRes_mu->Fill(momRes);
0407                           trackMomResP_mu->Fill(momRes, trueMom.Mag());
0408                           trackMomResEta_mu->Fill(momRes, trueEta);
0409                         }
0410                         else if( pdg == 211){
0411                           trackMomentumRes_pi->Fill(momRes);
0412                           trackMomResP_pi->Fill(momRes, trueMom.Mag());
0413                           trackMomResEta_pi->Fill(momRes, trueEta);
0414                         }
0415                         else if( pdg == 321){
0416                           trackMomentumRes_K->Fill(momRes);
0417                           trackMomResP_K->Fill(momRes, trueMom.Mag());
0418                           trackMomResEta_K->Fill(momRes, trueEta);
0419                         }
0420                         else if( pdg == 2212){
0421                           trackMomentumRes_p->Fill(momRes);
0422                           trackMomResP_p->Fill(momRes, trueMom.Mag());
0423                           trackMomResEta_p->Fill(momRes, trueEta);
0424                         }
0425
0426                         matchedPartTrackDeltaEta->Fill(deltaEta);
0427                         matchedPartTrackDeltaPhi->Fill(deltaPhi);
0428                         matchedPartTrackDeltaR->Fill(deltaR);
0429                         matchedPartTrackDeltaMom->Fill(deltaMom);
0430
0431                       }
0432                   }// End loop over associations
0433               } // End PDG check
0434           } // End stable particles condition
0435       } // End loop over thrown particles
0436   } // End loop over events
0437
0438   ofile->Write(); // Write histograms to file
0439   ofile->Close(); // Close output file
0440 }
0441 ```
0442
0443 Insert your input file path and execute as the example code above.
0444 -->
0445
0446 ### Compiled ROOT Scripts
0447
0448 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.
0449
0450 Each file is uploaded invidually, but your directory should be structured as follows -
0451
0452 - helloroot
0453     - README.md
0454     - CMakeLists.txt
0455     - include
0456       - helloroot
0457         - helloroot.hh
0458     - src
0459         - helloexec.cxx
0460         - helloroot.cxx
0461
0462 Note that any entry in the above without a file extension is a directory.
0463
0464 The contents of README.md are -
0465
0466 ```
0467 To build using cmake, create a build directory, navigate to it and run cmake. e.g.:
0468
0469 \```
0470 mkdir build
0471 cd build
0472 cmake ..
0473 make
0474 \```
0475 You can specify a number of parallel build threads with the -j flag, e.g.
0476 \```
0477 make -j4
0478 \```
0479
0480 You can specify an install directory to cmake with
0481 -DCMAKE_INSTALL_PREFIX=<path>
0482 then, after building,
0483 \```
0484 make install
0485 \```
0486 to install the headers and libraries under that location.
0487 There is no "make uninstall" but (on Unix-like systems)
0488 you can do
0489 xargs rm < install_manifest.txt
0490 from the cmake build directory.
0491 ```
0492
0493 Note that you should delete the \ characters in this block.
0494
0495 The contents of CMakeLists.txt are -
0496
0497 ```c++
0498 # CMakeLists.txt for helloroot.
0499 # More complicated than needed but demonstrates making and linking your own libraries
0500 # cf. https://cliutils.gitlab.io/modern-cmake/chapters/packages/ROOT.html
0501 # https://root.cern/manual/integrate_root_into_my_cmake_project/
0502
0503 cmake_minimum_required(VERSION 3.10)
0504 project(helloroot VERSION 1.0 LANGUAGES CXX ) # not needed
0505
0506  # Find ROOT. Use at least 6.20 for smoother cmake support
0507  find_package(ROOT 6.20 REQUIRED )
0508
0509  message ( " ROOT Libraries = " ${ROOT_LIBRARIES} )
0510
0511  ##############################################################################################################
0512
0513 # Main target is the libhelloroot library
0514 add_library(
0515   # You can use wildcards but it's cleaner to list the files explicitly
0516    helloroot
0517    SHARED
0518    src/helloroot.cxx
0519    )
0520 ## The particular syntax here is a bit annoying because you have to list all the sub-modules you need
0521 ## but it picks up automatically all the compile options needed for root, e.g. the c++ std version
0522 ## Find all available ROOT modules with `root-config --libs`
0523 target_link_libraries(helloroot PUBLIC ROOT::Core ROOT::RIO ROOT::Rint ROOT::Tree ROOT::EG ROOT::Physics )
0524
0525 ## The above _should_ be true, and it is on most systems. If it's not, uncoment one of the following lines
0526 # target_compile_features(helloroot PUBLIC cxx_std_17)
0527 # target_compile_features(helloroot PUBLIC cxx_std_20)
0528
0529 # include directories - this is also overkill but useful if you want to create dictionaries
0530 # Contact kkauder@gmail.com for that - it's too much for this example
0531 target_include_directories(helloroot
0532 PUBLIC
0533 $<INSTALL_INTERFACE:include>
0534 $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
0535  )
0536
0537 # Can add addtional options here
0538 target_compile_options(helloroot PRIVATE -Wall -Wextra -pedantic -g)
0539
0540 ##############################################################################################################
0541
0542 ## Build executables
0543 add_executable(helloexec src/helloexec.cxx)
0544 # target_compile_options(helloexec PRIVATE -Wall -Wextra -pedantic -g)
0545 target_link_libraries(helloexec helloroot )
0546 target_include_directories(helloexec
0547   PRIVATE
0548   ${ROOT_INCLUDE_DIRS}
0549   )
0550
0551 install(TARGETS helloexec DESTINATION bin)
0552
0553 ##############################################################################################################
0554
0555 ## Install library
0556 # Could also use include(GNUInstallDirs)
0557 # and then destinations of the form ${CMAKE_INSTALL_INCLUDEDIR}
0558 install(TARGETS helloroot
0559   EXPORT helloroot-export
0560   LIBRARY DESTINATION lib
0561   ARCHIVE DESTINATION lib
0562   )
0563
0564 ## Install headers
0565 install (DIRECTORY ${CMAKE_SOURCE_DIR}/include/helloroot
0566   DESTINATION  include/helloroot
0567   )
0568
0569 ## Generate configuration file - this allows you to use cmake in another project
0570 ## to find and link the installed helloroot library
0571 install(EXPORT helloroot-export
0572   FILE
0573   hellorootConfig.cmake
0574   NAMESPACE
0575     helloroot::
0576   DESTINATION
0577   cmake
0578   )
0579
0580 ## Final message
0581 message( " Done!")
0582 ```
0583
0584 The contents of helloroot.hh are -
0585
0586 ```c++
0587 #ifndef HELLO_ROOT_H
0588 #define HELLO_ROOT_H
0589
0590 void HelloRoot();
0591
0592 #endif // HELLO_ROOT_H
0593 ```
0594
0595 The contents of helloexec.cxx are -
0596
0597 ```c++
0598 #include<helloroot/helloroot.hh>
0599
0600 #include<iostream>
0601 #include<string>
0602
0603 int main()
0604 {
0605     std::cout << "Hello from main " << std::endl;
0606     HelloRoot();
0607
0608     return 0;
0609 }
0610 ```
0611
0612 And finally, the contents of helloroot.cxx are -
0613
0614 ```c++
0615 #include<helloroot/helloroot.hh>
0616
0617 #include<iostream>
0618 #include<string>
0619
0620 #include<TH1D.h>
0621 #include<TPad.h>
0622
0623 void HelloRoot()
0624 {
0625   std::cout << "Hello from HelloRoot" << std::endl;
0626
0627   // do something with root
0628   TH1D h("h", "h", 100, -5, 5);
0629   h.FillRandom("gaus", 1000);
0630   h.Draw();
0631   gPad->SaveAs("hello.png");
0632
0633   return;
0634 }
0635 ```
0636 Please consult the README and script comments for further instructions.
0637
0638 ## Python Uproot Scripts
0639
0640 ### EfficiencyAnalysis.py
0641
0642 Create a file called `EfficiencyAnalysis.py` and copy in the code below to get started on the resolution analysis exercise. Note that you will need to correctly specifiy your input file path in the variable `infile`.
0643
0644 ```python
0645 #! /usr/bin/python
0646
0647 #Import relevant packages
0648 import ROOT, math, array
0649 from ROOT import TH1F, TH2F, TMath, TTree, TVector3, TVector2
0650 import uproot as up
0651
0652 #Define and open files
0653 infile="PATH_TO_INPUT_FILE"
0654 ofile=ROOT.TFile.Open("EfficiencyAnalysis_OutPy.root", "RECREATE")
0655
0656 # Open input file and define branches we want to look at with uproot
0657 events_tree = up.open(infile)["events"]
0658
0659 # Get particle information
0660 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0661 partMomX = events_tree["MCParticles.momentum.x"].array()
0662 partMomY = events_tree["MCParticles.momentum.y"].array()
0663 partMomZ = events_tree["MCParticles.momentum.z"].array()
0664 partPdg = events_tree["MCParticles.PDG"].array()
0665
0666 # Get reconstructed track information
0667 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0668 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0669 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0670
0671 # Get assocations between MCParticles and ReconstructedChargedParticles
0672 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0673 simuAssoc = events_tree["ReconstructedChargedParticleAssociations.simID"].array()
0674
0675 # Define histograms below
0676 partEta = ROOT.TH1D("partEta","Eta of Thrown Charged Particles;Eta",100, -5 ,5 )
0677 matchedPartEta = ROOT.TH1D("matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track", 100, -5 ,5);
0678 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charge Particle", 5000, 0, 5);
0679
0680 # Add main analysis loop(s) below
0681 for i in range(0, len(partGenStat)): # Loop over all events
0682     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0683         if partGenStat[i][j] == 1: # Select stable particles
0684             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0685             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0686                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0687                 trueEta = trueMom.PseudoRapidity()
0688                 truePhi = trueMom.Phi()
0689                 partEta.Fill(trueEta)
0690                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0691                     if (simuAssoc[i][k] == j):
0692                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0693                         deltaEta = trueEta - recMom.PseudoRapidity()
0694                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
0695                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
0696                         matchedPartEta.Fill(trueEta)
0697                         matchedPartTrackDeltaR.Fill(deltaR)
0698
0699 # Write output histograms to file below
0700 partEta.Write()
0701 matchedPartEta.Write()
0702 matchedPartTrackDeltaR.Write()
0703
0704 # Close files
0705 ofile.Close()
0706 ```
0707 <!--
0708
0709 A "solution" version of the script for the exercise is included below -
0710
0711 ```python
0712 #! /usr/bin/python
0713
0714 #Import relevant packages
0715 import ROOT, math, array
0716 from ROOT import TCanvas, TColor, TGaxis, TH1F, TH2F, TPad, TStyle, gStyle, gPad, TGaxis, TLine, TMath, TPaveText, TTree, TVector3, TVector2
0717 import uproot as up
0718
0719 #Define and open files
0720 infile="PATH_TO_FILE"
0721 ofile=ROOT.TFile.Open("EfficiencyAnalysis_Exercise_OutPy.root", "RECREATE")
0722
0723 # Open input file and define branches we want to look at with uproot
0724 events_tree = up.open(infile)["events"]
0725
0726 # Get particle information
0727 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0728 partMomX = events_tree["MCParticles.momentum.x"].array()
0729 partMomY = events_tree["MCParticles.momentum.y"].array()
0730 partMomZ = events_tree["MCParticles.momentum.z"].array()
0731 partPdg = events_tree["MCParticles.PDG"].array()
0732
0733 # Get reconstructed track information
0734 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0735 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0736 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0737
0738 # Get assocations between MCParticles and ReconstructedChargedParticles
0739 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0740 simuAssoc = events_tree["ReconstructedChargedParticleAssociations.simID"].array()
0741
0742 # Define histograms below
0743 partEta = ROOT.TH1D("partEta","#eta of Thrown Charged Particles; #eta", 120, -6, 6)
0744 matchedPartEta = ROOT.TH1D("matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track; #eta", 120, -6, 6)
0745 partMom = ROOT.TH1D("partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150)
0746 matchedPartMom = ROOT.TH1D("matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150)
0747 partPhi = ROOT.TH1D("partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2)
0748 matchedPartPhi = ROOT.TH1D("matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2)
0749
0750 partPEta = ROOT.TH2D("partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6)
0751 matchedPartPEta = ROOT.TH2D("matchedPartPEta", "P vs #eta of Thrown Charged Particles, with matching track; P(GeV/C); #eta", 150, 0, 150, 120, -6, 6)
0752 partPhiEta = ROOT.TH2D("partPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
0753 matchedPartPhiEta = ROOT.TH2D("matchedPartPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
0754
0755 matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
0756 matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
0757 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
0758 matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
0759
0760 # Define some histograms for our efficiencies
0761 TrackEff_Eta = ROOT.TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6)
0762 TrackEffMom = ROOT.TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150)
0763 TrackEffPhi = ROOT.TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2)
0764
0765 # 2D Efficiencies
0766 TrackEff_PEta = ROOT.TH2D("TrackEff_PEta", "Tracking efficiency as fn of P and #eta; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6)
0767 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)
0768
0769 # All charged particle histos
0770 ChargedEta = ROOT.TH1D("ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6)
0771 ChargedPhi = ROOT.TH1D("ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2)
0772 ChargedP = ROOT.TH1D("ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150)
0773
0774 # Add main analysis loop(s) below
0775 for i in range(0, len(partGenStat)): # Loop over all events
0776     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0777         if partGenStat[i][j] == 1: # Select stable particles
0778             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0779             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0780                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0781                 trueEta = trueMom.PseudoRapidity()
0782                 truePhi = trueMom.Phi()
0783                 partEta.Fill(trueEta)
0784                 partPhi.Fill(truePhi)
0785                 partMom.Fill(trueMom.Mag())
0786                 partPEta.Fill(trueMom.Mag(), trueEta)
0787                 partPhiEta.Fill(truePhi, trueEta)
0788                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0789                     if (simuAssoc[i][k] == j):
0790                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0791                         deltaEta = trueEta - recMom.PseudoRapidity()
0792                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
0793                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
0794                         deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
0795                         matchedPartTrackDeltaEta.Fill(deltaEta)
0796                         matchedPartTrackDeltaPhi.Fill(deltaPhi)
0797                         matchedPartTrackDeltaR.Fill(deltaR)
0798                         matchedPartTrackDeltaMom.Fill(deltaMom)
0799                         matchedPartEta.Fill(trueEta)
0800                         matchedPartPhi.Fill(truePhi)
0801                         matchedPartMom.Fill(trueMom.Mag())
0802                         matchedPartPEta.Fill(trueMom.Mag(), trueEta)
0803                         matchedPartPhiEta.Fill(truePhi, trueEta)
0804     for x in range (0, len(trackMomX[i])): # Loop over all charged particles, thrown or not
0805         CPartMom = ROOT.TVector3(trackMomX[i][x], trackMomY[i][x], trackMomZ[i][x])
0806         CPartEta = CPartMom.PseudoRapidity()
0807         CPartPhi = CPartMom.Phi()
0808         ChargedEta.Fill(CPartEta)
0809         ChargedPhi.Fill(CPartPhi)
0810         ChargedP.Fill(CPartMom.Mag())
0811
0812 # Write output histograms to file below
0813 partEta.Write()
0814 matchedPartEta.Write()
0815 partMom.Write()
0816 matchedPartMom.Write()
0817 partPhi.Write()
0818 matchedPartPhi.Write()
0819 partPEta.Write()
0820 matchedPartPEta.Write()
0821 partPhiEta.Write()
0822 matchedPartPhiEta.Write()
0823 matchedPartTrackDeltaEta.Write()
0824 matchedPartTrackDeltaPhi.Write()
0825 matchedPartTrackDeltaR.Write()
0826 matchedPartTrackDeltaMom.Write()
0827 ChargedEta.Write()
0828 ChargedPhi.Write()
0829 ChargedP.Write()
0830 TrackEff_Eta.Divide(matchedPartEta, partEta, 1, 1, "b")
0831 TrackEffMom.Divide(matchedPartMom, partMom, 1, 1, "b")
0832 TrackEffPhi.Divide(matchedPartPhi, partPhi, 1, 1, "b")
0833 TrackEff_PEta.Divide(matchedPartPEta, partPEta, 1, 1, "b")
0834 TrackEff_PhiEta.Divide(matchedPartPhiEta, partPhiEta, 1, 1, "b")
0835 TrackEff_Eta.Write()
0836 TrackEffMom.Write()
0837 TrackEffPhi.Write()
0838 TrackEff_PEta.Write()
0839 TrackEff_PhiEta.Write()
0840
0841 # Close files
0842 ofile.Close()
0843 ```
0844 Insert your input file path and execute as the example code above.
0845 -->
0846 ### ResolutionAnalysis.py
0847
0848 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 specifiy your input file path in the variable `infile`.
0849
0850 ```python
0851 #! /usr/bin/python
0852
0853 #Import relevant packages
0854 import ROOT, math, array
0855 from ROOT import TH1F, TH2F, TMath, TTree, TVector3, TVector2
0856 import uproot as up
0857
0858 #Define and open files
0859 infile="PATH_TO_INPUT_FILE"
0860 ofile=ROOT.TFile.Open("ResolutionAnalysis_OutPy.root", "RECREATE")
0861
0862 # Open input file and define branches we want to look at with uproot
0863 events_tree = up.open(infile)["events"]
0864
0865 # Get particle information
0866 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0867 partMomX = events_tree["MCParticles.momentum.x"].array()
0868 partMomY = events_tree["MCParticles.momentum.y"].array()
0869 partMomZ = events_tree["MCParticles.momentum.z"].array()
0870 partPdg = events_tree["MCParticles.PDG"].array()
0871
0872 # Get reconstructed track information
0873 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0874 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0875 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0876
0877 # Get assocations between MCParticles and ReconstructedChargedParticles
0878 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0879 simuAssoc = events_tree["ReconstructedChargedParticleAssociations.simID"].array()
0880
0881 # Define histograms below
0882 trackMomentumRes = ROOT.TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2)
0883
0884 matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
0885 matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
0886 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
0887 matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
0888
0889 # Add main analysis loop(s) below
0890 for i in range(0, len(partGenStat)): # Loop over all events
0891     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0892         if partGenStat[i][j] == 1: # Select stable particles
0893             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0894             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0895                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0896                 trueEta = trueMom.PseudoRapidity()
0897                 truePhi = trueMom.Phi()
0898                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0899                     if (simuAssoc[i][k] == j):
0900                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0901                         deltaEta = trueEta - recMom.PseudoRapidity()
0902                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
0903                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
0904                         deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
0905                         momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag()
0906                         matchedPartTrackDeltaEta.Fill(deltaEta)
0907                         matchedPartTrackDeltaPhi.Fill(deltaPhi)
0908                         matchedPartTrackDeltaR.Fill(deltaR)
0909                         matchedPartTrackDeltaMom.Fill(deltaMom)
0910                         trackMomentumRes.Fill(momRes)
0911
0912 # Write output histograms to file below
0913 trackMomentumRes.Write()
0914 matchedPartTrackDeltaEta.Write()
0915 matchedPartTrackDeltaPhi.Write()
0916 matchedPartTrackDeltaR.Write()
0917 matchedPartTrackDeltaMom.Write()
0918
0919 # Close files
0920 ofile.Close()
0921 ```
0922 <!--
0923
0924 A "solution" version of the script for the exercise is included below -
0925
0926 ```python
0927 #! /usr/bin/python
0928
0929 #Import relevant packages
0930 import ROOT, math, array
0931 from ROOT import TCanvas, TColor, TGaxis, TH1F, TH2F, TPad, TStyle, gStyle, gPad, TGaxis, TLine, TMath, TPaveText, TTree, TVector3, TVector2
0932 import uproot as up
0933
0934 #Define and open files
0935 infile="PATH_TO_FILE"
0936 ofile=ROOT.TFile.Open("ResolutionAnalysis_Exercise_OutPy.root", "RECREATE")
0937
0938 # Open input file and define branches we want to look at with uproot
0939 events_tree = up.open(infile)["events"]
0940
0941 # Get particle information
0942 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0943 partMomX = events_tree["MCParticles.momentum.x"].array()
0944 partMomY = events_tree["MCParticles.momentum.y"].array()
0945 partMomZ = events_tree["MCParticles.momentum.z"].array()
0946 partPdg = events_tree["MCParticles.PDG"].array()
0947
0948 # Get reconstructed track information
0949 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0950 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0951 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0952
0953 # Get assocations between MCParticles and ReconstructedChargedParticles
0954 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0955 simuAssoc = events_tree.["ReconstructedChargedParticleAssociations.simID"].array()
0956
0957 # Define histograms below
0958 trackMomentumRes = ROOT.TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2)
0959 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);
0960 trackMomResEta = ROOT.TH2D("trackMomResEta", "Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
0961
0962 trackMomentumRes_e = ROOT.TH1D("trackMomentumRes_e","e^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0963 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);
0964 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);
0965
0966 trackMomentumRes_mu = ROOT.TH1D("trackMomentumRes_mu","#mu^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0967 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);
0968 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);
0969
0970 trackMomentumRes_pi = ROOT.TH1D("trackMomentumRes_pi","#pi^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0971 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);
0972 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);
0973
0974 trackMomentumRes_K = ROOT.TH1D("trackMomentumRes_K","K^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0975 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);
0976 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);
0977
0978 trackMomentumRes_p = ROOT.TH1D("trackMomentumRes_p","p Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0979 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);
0980 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);
0981
0982 matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
0983 matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
0984 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
0985 matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
0986
0987 # Add main analysis loop(s) below
0988 for i in range(0, len(partGenStat)): # Loop over all events
0989     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0990         if partGenStat[i][j] == 1: # Select stable particles
0991             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0992             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0993                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0994                 trueEta = trueMom.PseudoRapidity()
0995                 truePhi = trueMom.Phi()
0996                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0997                     if (simuAssoc[i][k] == j):
0998                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0999                         deltaEta = trueEta - recMom.PseudoRapidity()
1000                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
1001                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
1002                         deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
1003                         momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag()
1004                         trackMomentumRes.Fill(momRes)
1005                         trackMomResP.Fill(momRes, trueMom.Mag())
1006                         trackMomResEta.Fill(momRes, trueEta)
1007                         if( pdg == 11):
1008                             trackMomentumRes_e.Fill(momRes)
1009                             trackMomResP_e.Fill(momRes, trueMom.Mag())
1010                             trackMomResEta_e.Fill(momRes, trueEta)
1011                         elif( pdg == 13):
1012                             trackMomentumRes_mu.Fill(momRes)
1013                             trackMomResP_mu.Fill(momRes, trueMom.Mag())
1014                             trackMomResEta_mu.Fill(momRes, trueEta)
1015                         elif( pdg == 211):
1016                             trackMomentumRes_pi.Fill(momRes)
1017                             trackMomResP_pi.Fill(momRes, trueMom.Mag())
1018                             trackMomResEta_pi.Fill(momRes, trueEta)
1019                         elif( pdg == 321):
1020                             trackMomentumRes_K.Fill(momRes)
1021                             trackMomResP_K.Fill(momRes, trueMom.Mag())
1022                             trackMomResEta_K.Fill(momRes, trueEta)
1023                         elif( pdg == 2212):
1024                             trackMomentumRes_p.Fill(momRes)
1025                             trackMomResP_p.Fill(momRes, trueMom.Mag())
1026                             trackMomResEta_p.Fill(momRes, trueEta)
1027                         matchedPartTrackDeltaEta.Fill(deltaEta)
1028                         matchedPartTrackDeltaPhi.Fill(deltaPhi)
1029                         matchedPartTrackDeltaR.Fill(deltaR)
1030                         matchedPartTrackDeltaMom.Fill(deltaMom)
1031
1032 # Write output histograms to file below
1033 trackMomentumRes.Write()
1034 trackMomResP.Write()
1035 trackMomResEta.Write()
1036 trackMomentumRes_e.Write()
1037 trackMomResP_e.Write()
1038 trackMomResEta_e.Write()
1039 trackMomentumRes_mu.Write()
1040 trackMomResP_mu.Write()
1041 trackMomResEta_mu.Write()
1042 trackMomentumRes_pi.Write()
1043 trackMomResP_pi.Write()
1044 trackMomResEta_pi.Write()
1045 trackMomentumRes_K.Write()
1046 trackMomResP_K.Write()
1047 trackMomResEta_K.Write()
1048 trackMomentumRes_p.Write()
1049 trackMomResP_p.Write()
1050 trackMomResEta_p.Write()
1051 matchedPartTrackDeltaEta.Write()
1052 matchedPartTrackDeltaPhi.Write()
1053 matchedPartTrackDeltaR.Write()
1054 matchedPartTrackDeltaMom.Write()
1055
1056 # Close files
1057 ofile.Close()
1058 ```
1059 Insert your input file path and execute as the example code above.
1060 -->
1061 ## RDataFrames Example
1062
1063 Note that only the initial stage of the efficiency example is presented here in RDF format. This example was kindly created by [Simon](https://github.com/simonge/EIC_Analysis/blob/main/Analysis-Tutorial/EfficiencyAnalysisRDF.C).
1064
1065 ### EfficiencyAnalysisRDF.C
1066
1067 Create a file called `EfficiencyAnalysisRDF.C` and paste the code below in. Remember to change the file path.
1068
1069 Execute this script via - `root -l -q EfficiencyAnalysisRDF.C++`. Do this within eic-shell or somewhere else with the correct EDM4hep/EDM4eic libraries installed.
1070
1071 ```c++
1072 #include <edm4hep/utils/vector_utils.h>
1073 #include <edm4hep/MCParticle.h>
1074 #include <edm4eic/ReconstructedParticle.h>
1075 #include <ROOT/RDataFrame.hxx>
1076 #include <ROOT/RVec.hxx>
1077 #include <TFile.h>
1078
1079 // Define aliases for the data types
1080 using MCP = edm4hep::MCParticleData;
1081 using RecoP = edm4eic::ReconstructedParticleData;
1082
1083 // Define function to vectorize the edm4hep::utils methods
1084 template <typename T>
1085 auto getEta = [](ROOT::VecOps::RVec<T> momenta) {
1086   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::eta(p.momentum); });
1087 };
1088
1089 template <typename T>
1090 auto getPhi = [](ROOT::VecOps::RVec<T> momenta) {
1091   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::angleAzimuthal(p.momentum); });
1092 };
1093
1094 // Define the function to perform the efficiency analysis
1095 void EfficiencyAnalysisRDF(TString infile="PATH_TO_FILE"){
1096
1097   // Set up input file
1098   ROOT::RDataFrame df("events", infile);
1099
1100   // Define new dataframe node with additional columns
1101   auto df1 =  df.Define("statusFilter",  "MCParticles.generatorStatus == 1"    )
1102                 .Define("absPDG",        "abs(MCParticles.PDG)"                )
1103                 .Define("pdgFilter",     "absPDG == 11 || absPDG == 13 || absPDG == 211 || absPDG == 321 || absPDG == 2212")
1104                 .Define("particleFilter","statusFilter && pdgFilter"           )
1105                 .Define("filtMCParts",   "MCParticles[particleFilter]"         )
1106                 .Define("assoFilter",    "Take(particleFilter,ReconstructedChargedParticleAssociations.simID)") // Incase any of the associated particles happen to not be charged
1107                 .Define("assoMCParts",   "Take(MCParticles,ReconstructedChargedParticleAssociations.simID)[assoFilter]")
1108                 .Define("assoRecParts",  "Take(ReconstructedChargedParticles,ReconstructedChargedParticleAssociations.recID)[assoFilter]")
1109                 .Define("filtMCEta",     getEta<MCP>   , {"filtMCParts"} )
1110                 .Define("filtMCPhi",     getPhi<MCP>   , {"filtMCParts"} )
1111                 .Define("accoMCEta",     getEta<MCP>   , {"assoMCParts"} )
1112                 .Define("accoMCPhi",     getPhi<MCP>   , {"assoMCParts"} )
1113                 .Define("assoRecEta",    getEta<RecoP> , {"assoRecParts"})
1114                 .Define("assoRecPhi",    getPhi<RecoP> , {"assoRecParts"})
1115                 .Define("deltaR",        "ROOT::VecOps::DeltaR(assoRecEta, accoMCEta, assoRecPhi, accoMCPhi)");
1116
1117   // Define histograms
1118   auto partEta                = df1.Histo1D({"partEta","Eta of Thrown Charged Particles;Eta",100,-5.,5.},"filtMCEta");
1119   auto matchedPartEta         = df1.Histo1D({"matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track",100,-5.,5.},"accoMCEta");
1120   auto matchedPartTrackDeltaR = df1.Histo1D({"matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charged Particle",5000,0.,5.},"deltaR");
1121
1122   // Write histograms to file
1123   TFile *ofile = TFile::Open("EfficiencyAnalysis_Out_RDF.root","RECREATE");
1124
1125   // Booked Define and Histo1D lazy actions are only performed here
1126   partEta->Write();
1127   matchedPartEta->Write();
1128   matchedPartTrackDeltaR->Write();
1129
1130   ofile->Close(); // Close output file
1131 }
1132 ```
1133 <!--
1134
1135 A "solution" using RDataFrames is included below, please see the notes following this script for some of my thoughts on RDataFrames -
1136
1137 ```c++
1138 #include <edm4hep/utils/vector_utils.h>
1139 #include <edm4hep/MCParticle.h>
1140 #include <edm4eic/ReconstructedParticle.h>
1141 #include <ROOT/RDataFrame.hxx>
1142 #include <ROOT/RVec.hxx>
1143 #include <TFile.h>
1144
1145 // Define aliases for the data types
1146 using MCP = edm4hep::MCParticleData;
1147 using RecoP = edm4eic::ReconstructedParticleData;
1148
1149 // Define function to vectorize the edm4hep::utils methods
1150 template <typename T>
1151 auto getEta = [](ROOT::VecOps::RVec<T> momenta) {
1152   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::eta(p.momentum); });
1153 };
1154
1155 template <typename T>
1156 auto getPhi = [](ROOT::VecOps::RVec<T> momenta) {
1157   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::angleAzimuthal(p.momentum); });
1158 };
1159
1160 template <typename T>
1161 auto getP = [](ROOT::VecOps::RVec<T> momenta) {
1162   //return ROOT::VecOps::Map(momenta, [](const T& p) { return (p.momentum); }); // This is a vector3f
1163   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::magnitude(p.momentum); }); // This is a the magnitude of that vector3f
1164 };
1165
1166 // Define the function to perform the efficiency analysis
1167 void EfficiencyAnalysisRDF_Exercise(TString infile="PATH_TO_INPUT_FILE"){
1168
1169   // Set up input file
1170   ROOT::RDataFrame df("events", infile);
1171
1172   // Define new dataframe node with additional columns
1173   auto df1 =  df.Define("statusFilter",  "MCParticles.generatorStatus == 1"    )
1174     .Define("absPDG",        "abs(MCParticles.PDG)"                )
1175     .Define("pdgFilter",     "absPDG == 11 || absPDG == 13 || absPDG == 211 || absPDG == 321 || absPDG == 2212")
1176     .Define("particleFilter","statusFilter && pdgFilter"           )
1177     .Define("filtMCParts",   "MCParticles[particleFilter]"         )
1178     .Define("assoFilter",    "Take(particleFilter,ReconstructedChargedParticleAssociations.simID)") // Incase any of the associated particles happen to not be charged
1179     .Define("assoMCParts",   "Take(MCParticles,ReconstructedChargedParticleAssociations.simID)[assoFilter]")
1180     .Define("assoRecParts",  "Take(ReconstructedChargedParticles,ReconstructedChargedParticleAssociations.recID)[assoFilter]")
1181     .Define("filtMCEta",     getEta<MCP>   , {"filtMCParts"} )
1182     .Define("filtMCPhi",     getPhi<MCP>   , {"filtMCParts"} )
1183     .Define("filtMCp",       getP<MCP>     , {"filtMCParts"} )
1184     .Define("assoMCEta",     getEta<MCP>   , {"assoMCParts"} )
1185     .Define("assoMCPhi",     getPhi<MCP>   , {"assoMCParts"} )
1186     .Define("assoMCp",       getP<MCP>     , {"assoMCParts"} )
1187     .Define("assoRecEta",    getEta<RecoP> , {"assoRecParts"})
1188     .Define("assoRecPhi",    getPhi<RecoP> , {"assoRecParts"})
1189     .Define("assoRecp",      getP<RecoP>   , {"assoRecParts"})
1190     .Define("deltaEta",      "assoMCEta - assoRecEta"        )
1191     .Define("deltaPhi",      "ROOT::VecOps::DeltaPhi(assoRecPhi, assoMCPhi)")
1192     .Define("deltaR",        "ROOT::VecOps::DeltaR(assoRecEta, assoMCEta, assoRecPhi, assoMCPhi)")
1193     .Define("deltaMom",      "assoMCp - assoRecp")
1194     .Define("recoEta",       getEta<RecoP>,  {"ReconstructedChargedParticles"})
1195     .Define("recoPhi",       getPhi<RecoP>,  {"ReconstructedChargedParticles"})
1196     .Define("recoP",         getP<RecoP>,    {"ReconstructedChargedParticles"});
1197
1198   // 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
1199   auto partEta                = df1.Histo1D({"partEta","#eta of Thrown Charged Particles;#eta",120,-6.,6.},"filtMCEta");
1200   auto matchedPartEta         = df1.Histo1D({"matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track;#eta",120,-6.,6.},"assoMCEta");
1201   auto partMom = df1.Histo1D({"partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150}, "filtMCp");
1202   auto matchedPartMom = df1.Histo1D({"matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150}, "assoMCp");
1203   auto partPhi = df1.Histo1D({"partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2},"filtMCPhi");
1204   auto matchedPartPhi = df1.Histo1D({"matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2}, "assoMCPhi");
1205
1206   auto partPEta = df1.Histo2D({"partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6}, "filtMCp", "filtMCEta");
1207   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");
1208   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");
1209   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");
1210
1211   auto matchedPartTrackDeltaEta = df1.Histo1D({"matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25}, "deltaEta");
1212   auto matchedPartTrackDeltaPhi = df1.Histo1D({"matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2}, "deltaPhi");
1213   auto matchedPartTrackDeltaR = df1.Histo1D({"matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle;#Delta R",300,0.,0.3}, "deltaR");
1214   auto matchedPartTrackDeltaMom = df1.Histo1D({"matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10}, "deltaMom");
1215
1216   // Define some histograms for our efficiencies - Done "old school" root style - Maybe the division can be done direct from a DF?
1217   TH1D *TrackEff_Eta = new TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6);
1218   TH1D *TrackEff_Mom = new TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150);
1219   TH1D *TrackEff_Phi = new TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2);
1220   // 2D Efficiencies
1221   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);
1222   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);
1223
1224   auto ChargedEta = df1.Histo1D({"ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6}, "recoEta");
1225   auto ChargedPhi = df1.Histo1D({"ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2}, "recoPhi");
1226   auto ChargedP = df1.Histo1D({"ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150}, "recoP");
1227
1228   // Write histograms to file
1229   TFile *ofile = TFile::Open("EfficiencyAnalysis_Exercise_Out_RDF.root","RECREATE");
1230
1231   // Booked Define and Histo1D lazy actions are only performed here
1232   partEta->Write();
1233   matchedPartEta->Write();
1234   partPhi->Write();
1235   matchedPartPhi->Write();
1236   partMom->Write();
1237   matchedPartMom->Write();
1238   partPEta->Write();
1239   matchedPartPEta->Write();
1240   partPhiEta->Write();
1241   matchedPartPhiEta->Write();
1242   matchedPartTrackDeltaEta->Write();
1243   matchedPartTrackDeltaPhi->Write();
1244   matchedPartTrackDeltaR->Write();
1245   matchedPartTrackDeltaMom->Write();
1246
1247   // Create efficiency histograms by dividing appropriately. Note we must actually get the pointer explicitly.
1248   TrackEff_Eta->Divide(matchedPartEta.GetPtr(), partEta.GetPtr(), 1, 1, "b");
1249   TrackEff_Mom->Divide(matchedPartMom.GetPtr(), partMom.GetPtr(), 1, 1, "b");
1250   TrackEff_Phi->Divide(matchedPartPhi.GetPtr(), partPhi.GetPtr(), 1, 1, "b");
1251   TrackEff_PEta->Divide(matchedPartPEta.GetPtr(), partPEta.GetPtr(), 1, 1, "b");
1252   TrackEff_PhiEta->Divide(matchedPartPhiEta.GetPtr(), partPhiEta.GetPtr(), 1, 1, "b");
1253
1254   TrackEff_Eta->Write();
1255   TrackEff_Mom->Write();
1256   TrackEff_Phi->Write();
1257   TrackEff_PEta->Write();
1258   TrackEff_PhiEta->Write();
1259
1260   ChargedEta->Write();
1261   ChargedPhi->Write();
1262   ChargedP->Write();
1263
1264   ofile->Close(); // Close output file
1265 }
1266 ```
1267 Insert your input file path and execute as the example code above.
1268
1269 > Note:
1270 > - Before writing this example and solution, I have never used RDataFrames before. Some thoughts below.
1271 >   - I find them very difficult to work with. The processes that are "actually" going on are very obscured with RDataFrames.
1272 >   - Finding resources and examples online is much more difficult due to how new RDataFrames are.
1273 >     - This further complicates writing and working with them.
1274 >   - Getting this example working and producing the same result as the python and ROOT examples took me as much time as preparing the rest of the tutorial.
1275 >   - The script itself seems to be very slow when running.
1276 >   - I don't think they're a good starting point for someone new to ROOT/Nuclear Physics data analysis. Too much is going on behind the scenes.
1277 >  - In summary, I personally would *not* recommend using RDataFrames for your analysis (at this point in time).
1278 >  - Note that this just my (SKay) thoughts and opinions though.
1279 {: .callout}
1280
1281 Note that due to how much "fun" I had making the efficiency analysis exercise with RDataFrames, I won't be creating a solution for the resolution analysis exercises.
1282 -->
1283
1284 {% include links.md %}