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
0556 ## Install library
0557 # Could also use include(GNUInstallDirs)
0558 # and then destinations of the form ${CMAKE_INSTALL_INCLUDEDIR}
0559 install(TARGETS helloroot
0560   EXPORT helloroot-export
0561   LIBRARY DESTINATION lib
0562   ARCHIVE DESTINATION lib
0563   )
0564
0565 ## Install headers
0566 install (DIRECTORY ${CMAKE_SOURCE_DIR}/include/helloroot
0567   DESTINATION  include/helloroot
0568   )
0569
0570 ## Generate configuration file - this allows you to use cmake in another project
0571 ## to find and link the installed helloroot library
0572 install(EXPORT helloroot-export
0573   FILE
0574   hellorootConfig.cmake
0575   NAMESPACE
0576     helloroot::
0577   DESTINATION
0578   cmake
0579   )
0580
0581 ## Final message
0582 message( " Done!")
0583 ```
0584
0585 The contents of helloroot.hh are -
0586
0587 ```c++
0588 #ifndef HELLO_ROOT_H
0589 #define HELLO_ROOT_H
0590
0591 void HelloRoot();
0592
0593 #endif // HELLO_ROOT_H
0594 ```
0595
0596 The contents of helloexec.cxx are -
0597
0598 ```c++
0599 #include<helloroot/helloroot.hh>
0600
0601 #include<iostream>
0602 #include<string>
0603
0604 int main()
0605 {
0606     std::cout << "Hello from main " << std::endl;
0607     HelloRoot();
0608
0609     return 0;
0610 }
0611 ```
0612
0613 And finally, the contents of helloroot.cxx are -
0614
0615 ```c++
0616 #include<helloroot/helloroot.hh>
0617
0618 #include<iostream>
0619 #include<string>
0620
0621 #include<TH1D.h>
0622 #include<TPad.h>
0623
0624 void HelloRoot()
0625 {
0626   std::cout << "Hello from HelloRoot" << std::endl;
0627
0628   // do something with root
0629   TH1D h("h", "h", 100, -5, 5);
0630   h.FillRandom("gaus", 1000);
0631   h.Draw();
0632   gPad->SaveAs("hello.png");
0633
0634   return;
0635 }
0636 ```
0637 Please consult the README and script comments for further instructions.
0638
0639 ## Python Uproot Scripts
0640
0641 ### EfficiencyAnalysis.py
0642
0643 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`.
0644
0645 ```python
0646 #! /usr/bin/python
0647
0648 #Import relevant packages
0649 import ROOT, math, array
0650 from ROOT import TH1F, TH2F, TMath, TTree, TVector3, TVector2
0651 import uproot as up
0652
0653 #Define and open files
0654 infile="PATH_TO_INPUT_FILE"
0655 ofile=ROOT.TFile.Open("EfficiencyAnalysis_OutPy.root", "RECREATE")
0656
0657 # Open input file and define branches we want to look at with uproot
0658 events_tree = up.open(infile)["events"]
0659
0660 # Get particle information
0661 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0662 partMomX = events_tree["MCParticles.momentum.x"].array()
0663 partMomY = events_tree["MCParticles.momentum.y"].array()
0664 partMomZ = events_tree["MCParticles.momentum.z"].array()
0665 partPdg = events_tree["MCParticles.PDG"].array()
0666
0667 # Get reconstructed track information
0668 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0669 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0670 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0671
0672 # Get assocations between MCParticles and ReconstructedChargedParticles
0673 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0674 simuAssoc = events_tree["ReconstructedChargedParticleAssociations.simID"].array()
0675
0676 # Define histograms below
0677 partEta = ROOT.TH1D("partEta","Eta of Thrown Charged Particles;Eta",100, -5 ,5 )
0678 matchedPartEta = ROOT.TH1D("matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track", 100, -5 ,5);
0679 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charge Particle", 5000, 0, 5);
0680
0681 # Add main analysis loop(s) below
0682 for i in range(0, len(events_tree)): # Loop over all events
0683     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0684         if partGenStat[i][j] == 1: # Select stable particles
0685             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0686             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0687                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0688                 trueEta = trueMom.PseudoRapidity()
0689                 truePhi = trueMom.Phi()
0690
0691                 partEta.Fill(trueEta)
0692                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0693                     if (simuAssoc[i][k] == j):
0694                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0695                         deltaEta = trueEta - recMom.PseudoRapidity()
0696                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
0697                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
0698
0699                         matchedPartEta.Fill(trueEta)
0700                         matchedPartTrackDeltaR.Fill(deltaR)
0701
0702 # Write output histograms to file below
0703 partEta.Write()
0704 matchedPartEta.Write()
0705 matchedPartTrackDeltaR.Write()
0706
0707 # Close files
0708 ofile.Close()
0709 ```
0710 <!--
0711
0712 A "solution" version of the script for the exercise is included below -
0713
0714 ```python
0715 #! /usr/bin/python
0716
0717 #Import relevant packages
0718 import ROOT, math, array
0719 from ROOT import TCanvas, TColor, TGaxis, TH1F, TH2F, TPad, TStyle, gStyle, gPad, TGaxis, TLine, TMath, TPaveText, TTree, TVector3, TVector2
0720 import uproot as up
0721
0722 #Define and open files
0723 infile="PATH_TO_FILE"
0724 ofile=ROOT.TFile.Open("EfficiencyAnalysis_Exercise_OutPy.root", "RECREATE")
0725
0726 # Open input file and define branches we want to look at with uproot
0727 events_tree = up.open(infile)["events"]
0728
0729 # Get particle information
0730 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0731 partMomX = events_tree["MCParticles.momentum.x"].array()
0732 partMomY = events_tree["MCParticles.momentum.y"].array()
0733 partMomZ = events_tree["MCParticles.momentum.z"].array()
0734 partPdg = events_tree["MCParticles.PDG"].array()
0735
0736 # Get reconstructed track information
0737 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0738 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0739 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0740
0741 # Get assocations between MCParticles and ReconstructedChargedParticles
0742 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0743 simuAssoc = events_tree["ReconstructedChargedParticleAssociations.simID"].array()
0744
0745 # Define histograms below
0746 partEta = ROOT.TH1D("partEta","#eta of Thrown Charged Particles; #eta", 120, -6, 6)
0747 matchedPartEta = ROOT.TH1D("matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track; #eta", 120, -6, 6)
0748 partMom = ROOT.TH1D("partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150)
0749 matchedPartMom = ROOT.TH1D("matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150)
0750 partPhi = ROOT.TH1D("partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2)
0751 matchedPartPhi = ROOT.TH1D("matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2)
0752
0753 partPEta = ROOT.TH2D("partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6)
0754 matchedPartPEta = ROOT.TH2D("matchedPartPEta", "P vs #eta of Thrown Charged Particles, with matching track; P(GeV/C); #eta", 150, 0, 150, 120, -6, 6)
0755 partPhiEta = ROOT.TH2D("partPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
0756 matchedPartPhiEta = ROOT.TH2D("matchedPartPhiEta", "#phi vs #eta of Thrown Charged Particles; #phi(rad); #eta", 160, -3.2, 3.2, 120, -6, 6)
0757
0758 matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
0759 matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
0760 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
0761 matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
0762
0763 # Define some histograms for our efficiencies
0764 TrackEff_Eta = ROOT.TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6)
0765 TrackEffMom = ROOT.TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150)
0766 TrackEffPhi = ROOT.TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2)
0767
0768 # 2D Efficiencies
0769 TrackEff_PEta = ROOT.TH2D("TrackEff_PEta", "Tracking efficiency as fn of P and #eta; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6)
0770 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)
0771
0772 # All charged particle histos
0773 ChargedEta = ROOT.TH1D("ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6)
0774 ChargedPhi = ROOT.TH1D("ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2)
0775 ChargedP = ROOT.TH1D("ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150)
0776
0777 # Add main analysis loop(s) below
0778 for i in range(0, len(events_tree)): # Loop over all events
0779     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0780         if partGenStat[i][j] == 1: # Select stable particles
0781             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0782             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0783                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0784                 trueEta = trueMom.PseudoRapidity()
0785                 truePhi = trueMom.Phi()
0786
0787                 partEta.Fill(trueEta)
0788                 partPhi.Fill(truePhi)
0789                 partMom.Fill(trueMom.Mag())
0790                 partPEta.Fill(trueMom.Mag(), trueEta)
0791                 partPhiEta.Fill(truePhi, trueEta)
0792                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0793                     if (simuAssoc[i][k] == j):
0794                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0795                         deltaEta = trueEta - recMom.PseudoRapidity()
0796                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
0797                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
0798                         deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
0799
0800                         matchedPartTrackDeltaEta.Fill(deltaEta)
0801                         matchedPartTrackDeltaPhi.Fill(deltaPhi)
0802                         matchedPartTrackDeltaR.Fill(deltaR)
0803                         matchedPartTrackDeltaMom.Fill(deltaMom)
0804
0805                         matchedPartEta.Fill(trueEta)
0806                         matchedPartPhi.Fill(truePhi)
0807                         matchedPartMom.Fill(trueMom.Mag())
0808                         matchedPartPEta.Fill(trueMom.Mag(), trueEta)
0809                         matchedPartPhiEta.Fill(truePhi, trueEta)
0810     for x in range (0, len(trackMomX[i])): # Loop over all charged particles, thrown or not
0811         CPartMom = ROOT.TVector3(trackMomX[i][x], trackMomY[i][x], trackMomZ[i][x])
0812         CPartEta = CPartMom.PseudoRapidity()
0813         CPartPhi = CPartMom.Phi()
0814
0815         ChargedEta.Fill(CPartEta)
0816         ChargedPhi.Fill(CPartPhi)
0817         ChargedP.Fill(CPartMom.Mag())
0818
0819 # Write output histograms to file below
0820 partEta.Write()
0821 matchedPartEta.Write()
0822 partMom.Write()
0823 matchedPartMom.Write()
0824 partPhi.Write()
0825 matchedPartPhi.Write()
0826 partPEta.Write()
0827 matchedPartPEta.Write()
0828 partPhiEta.Write()
0829 matchedPartPhiEta.Write()
0830 matchedPartTrackDeltaEta.Write()
0831 matchedPartTrackDeltaPhi.Write()
0832 matchedPartTrackDeltaR.Write()
0833 matchedPartTrackDeltaMom.Write()
0834 ChargedEta.Write()
0835 ChargedPhi.Write()
0836 ChargedP.Write()
0837
0838 TrackEff_Eta.Divide(matchedPartEta, partEta, 1, 1, "b")
0839 TrackEffMom.Divide(matchedPartMom, partMom, 1, 1, "b")
0840 TrackEffPhi.Divide(matchedPartPhi, partPhi, 1, 1, "b")
0841 TrackEff_PEta.Divide(matchedPartPEta, partPEta, 1, 1, "b")
0842 TrackEff_PhiEta.Divide(matchedPartPhiEta, partPhiEta, 1, 1, "b")
0843
0844 TrackEff_Eta.Write()
0845 TrackEffMom.Write()
0846 TrackEffPhi.Write()
0847 TrackEff_PEta.Write()
0848 TrackEff_PhiEta.Write()
0849
0850 # Close files
0851 ofile.Close()
0852 ```
0853 Insert your input file path and execute as the example code above.
0854 -->
0855 ### ResolutionAnalysis.py
0856
0857 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`.
0858
0859 ```python
0860 #! /usr/bin/python
0861
0862 #Import relevant packages
0863 import ROOT, math, array
0864 from ROOT import TH1F, TH2F, TMath, TTree, TVector3, TVector2
0865 import uproot as up
0866
0867 #Define and open files
0868 infile="PATH_TO_INPUT_FILE"
0869 ofile=ROOT.TFile.Open("ResolutionAnalysis_OutPy.root", "RECREATE")
0870
0871 # Open input file and define branches we want to look at with uproot
0872 events_tree = up.open(infile)["events"]
0873
0874 # Get particle information
0875 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0876 partMomX = events_tree["MCParticles.momentum.x"].array()
0877 partMomY = events_tree["MCParticles.momentum.y"].array()
0878 partMomZ = events_tree["MCParticles.momentum.z"].array()
0879 partPdg = events_tree["MCParticles.PDG"].array()
0880
0881 # Get reconstructed track information
0882 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0883 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0884 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0885
0886 # Get assocations between MCParticles and ReconstructedChargedParticles
0887 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0888 simuAssoc = events_tree["ReconstructedChargedParticleAssociations.simID"].array()
0889
0890 # Define histograms below
0891 trackMomentumRes = ROOT.TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2)
0892
0893 matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
0894 matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
0895 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
0896 matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
0897
0898 # Add main analysis loop(s) below
0899 for i in range(0, len(events_tree)): # Loop over all events
0900     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
0901         if partGenStat[i][j] == 1: # Select stable particles
0902             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
0903             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
0904                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
0905                 trueEta = trueMom.PseudoRapidity()
0906                 truePhi = trueMom.Phi()
0907
0908                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
0909                     if (simuAssoc[i][k] == j):
0910                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
0911                         deltaEta = trueEta - recMom.PseudoRapidity()
0912                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
0913                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
0914                         deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
0915
0916                         momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag()
0917
0918                         matchedPartTrackDeltaEta.Fill(deltaEta)
0919                         matchedPartTrackDeltaPhi.Fill(deltaPhi)
0920                         matchedPartTrackDeltaR.Fill(deltaR)
0921                         matchedPartTrackDeltaMom.Fill(deltaMom)
0922
0923                         trackMomentumRes.Fill(momRes)
0924
0925 # Write output histograms to file below
0926 trackMomentumRes.Write()
0927 matchedPartTrackDeltaEta.Write()
0928 matchedPartTrackDeltaPhi.Write()
0929 matchedPartTrackDeltaR.Write()
0930 matchedPartTrackDeltaMom.Write()
0931
0932 # Close files
0933 ofile.Close()
0934 ```
0935 <!--
0936
0937 A "solution" version of the script for the exercise is included below -
0938
0939 ```python
0940 #! /usr/bin/python
0941
0942 #Import relevant packages
0943 import ROOT, math, array
0944 from ROOT import TCanvas, TColor, TGaxis, TH1F, TH2F, TPad, TStyle, gStyle, gPad, TGaxis, TLine, TMath, TPaveText, TTree, TVector3, TVector2
0945 import uproot as up
0946
0947 #Define and open files
0948 infile="PATH_TO_FILE"
0949 ofile=ROOT.TFile.Open("ResolutionAnalysis_Exercise_OutPy.root", "RECREATE")
0950
0951 # Open input file and define branches we want to look at with uproot
0952 events_tree = up.open(infile)["events"]
0953
0954 # Get particle information
0955 partGenStat = events_tree["MCParticles.generatorStatus"].array()
0956 partMomX = events_tree["MCParticles.momentum.x"].array()
0957 partMomY = events_tree["MCParticles.momentum.y"].array()
0958 partMomZ = events_tree["MCParticles.momentum.z"].array()
0959 partPdg = events_tree["MCParticles.PDG"].array()
0960
0961 # Get reconstructed track information
0962 trackMomX = events_tree["ReconstructedChargedParticles.momentum.x"].array()
0963 trackMomY = events_tree["ReconstructedChargedParticles.momentum.y"].array()
0964 trackMomZ = events_tree["ReconstructedChargedParticles.momentum.z"].array()
0965
0966 # Get assocations between MCParticles and ReconstructedChargedParticles
0967 recoAssoc = events_tree["ReconstructedChargedParticleAssociations.recID"].array()
0968 simuAssoc = events_tree.["ReconstructedChargedParticleAssociations.simID"].array()
0969
0970 # Define histograms below
0971 trackMomentumRes = ROOT.TH1D("trackMomentumRes","Track Momentum Resolution", 400, -2, 2)
0972 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);
0973 trackMomResEta = ROOT.TH2D("trackMomResEta", "Track Momentum Resolution vs #eta; (P_{rec} - P_{MC})/P_{MC}; #eta_{MC}", 400, -2, 2, 120, -6, 6);
0974
0975 trackMomentumRes_e = ROOT.TH1D("trackMomentumRes_e","e^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0976 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);
0977 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);
0978
0979 trackMomentumRes_mu = ROOT.TH1D("trackMomentumRes_mu","#mu^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0980 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);
0981 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);
0982
0983 trackMomentumRes_pi = ROOT.TH1D("trackMomentumRes_pi","#pi^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0984 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);
0985 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);
0986
0987 trackMomentumRes_K = ROOT.TH1D("trackMomentumRes_K","K^{#pm} Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0988 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);
0989 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);
0990
0991 trackMomentumRes_p = ROOT.TH1D("trackMomentumRes_p","p Track Momentum Resolution; (P_{rec} - P_{MC})/P_{MC}", 400, -2, 2);
0992 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);
0993 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);
0994
0995 matchedPartTrackDeltaEta = ROOT.TH1D("matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructe Charged Particle; #Delta#eta", 100, -0.25, 0.25)
0996 matchedPartTrackDeltaPhi = ROOT.TH1D("matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2)
0997 matchedPartTrackDeltaR = ROOT.TH1D("matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle; #Delta R", 300, 0, 0.3)
0998 matchedPartTrackDeltaMom = ROOT.TH1D("matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10)
0999
1000 # Add main analysis loop(s) below
1001 for i in range(0, len(events_tree)): # Loop over all events
1002     for j in range(0, len(partGenStat[i])): # Loop over all thrown particles
1003         if partGenStat[i][j] == 1: # Select stable particles
1004             pdg = abs(partPdg[i][j]) # Get PDG for each stable particle
1005             if(pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212):
1006                 trueMom = ROOT.TVector3(partMomX[i][j], partMomY[i][j], partMomZ[i][j])
1007                 trueEta = trueMom.PseudoRapidity()
1008                 truePhi = trueMom.Phi()
1009
1010                 for k in range(0,len(simuAssoc[i])): # Loop over associations to find matching ReconstructedChargedParticle
1011                     if (simuAssoc[i][k] == j):
1012                         recMom = ROOT.TVector3(trackMomX[i][recoAssoc[i][k]], trackMomY[i][recoAssoc[i][k]], trackMomZ[i][recoAssoc[i][k]])
1013                         deltaEta = trueEta - recMom.PseudoRapidity()
1014                         deltaPhi = TVector2. Phi_mpi_pi(truePhi - recMom.Phi())
1015                         deltaR = math.sqrt((deltaEta*deltaEta) + (deltaPhi*deltaPhi))
1016                         deltaMom = ((trueMom.Mag()) - (recMom.Mag()))
1017
1018                         momRes = (recMom.Mag() - trueMom.Mag())/trueMom.Mag()
1019
1020                         trackMomentumRes.Fill(momRes)
1021                         trackMomResP.Fill(momRes, trueMom.Mag())
1022                         trackMomResEta.Fill(momRes, trueEta)
1023
1024                         if( pdg == 11):
1025                             trackMomentumRes_e.Fill(momRes)
1026                             trackMomResP_e.Fill(momRes, trueMom.Mag())
1027                             trackMomResEta_e.Fill(momRes, trueEta)
1028                         elif( pdg == 13):
1029                             trackMomentumRes_mu.Fill(momRes)
1030                             trackMomResP_mu.Fill(momRes, trueMom.Mag())
1031                             trackMomResEta_mu.Fill(momRes, trueEta)
1032                         elif( pdg == 211):
1033                             trackMomentumRes_pi.Fill(momRes)
1034                             trackMomResP_pi.Fill(momRes, trueMom.Mag())
1035                             trackMomResEta_pi.Fill(momRes, trueEta)
1036                         elif( pdg == 321):
1037                             trackMomentumRes_K.Fill(momRes)
1038                             trackMomResP_K.Fill(momRes, trueMom.Mag())
1039                             trackMomResEta_K.Fill(momRes, trueEta)
1040                         elif( pdg == 2212):
1041                             trackMomentumRes_p.Fill(momRes)
1042                             trackMomResP_p.Fill(momRes, trueMom.Mag())
1043                             trackMomResEta_p.Fill(momRes, trueEta)
1044
1045                         matchedPartTrackDeltaEta.Fill(deltaEta)
1046                         matchedPartTrackDeltaPhi.Fill(deltaPhi)
1047                         matchedPartTrackDeltaR.Fill(deltaR)
1048                         matchedPartTrackDeltaMom.Fill(deltaMom)
1049
1050 # Write output histograms to file below
1051 trackMomentumRes.Write()
1052 trackMomResP.Write()
1053 trackMomResEta.Write()
1054 trackMomentumRes_e.Write()
1055 trackMomResP_e.Write()
1056 trackMomResEta_e.Write()
1057 trackMomentumRes_mu.Write()
1058 trackMomResP_mu.Write()
1059 trackMomResEta_mu.Write()
1060 trackMomentumRes_pi.Write()
1061 trackMomResP_pi.Write()
1062 trackMomResEta_pi.Write()
1063 trackMomentumRes_K.Write()
1064 trackMomResP_K.Write()
1065 trackMomResEta_K.Write()
1066 trackMomentumRes_p.Write()
1067 trackMomResP_p.Write()
1068 trackMomResEta_p.Write()
1069
1070 matchedPartTrackDeltaEta.Write()
1071 matchedPartTrackDeltaPhi.Write()
1072 matchedPartTrackDeltaR.Write()
1073 matchedPartTrackDeltaMom.Write()
1074
1075 # Close files
1076 ofile.Close()
1077 ```
1078 Insert your input file path and execute as the example code above.
1079 -->
1080 ## RDataFrames Example
1081
1082 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).
1083
1084 ### EfficiencyAnalysisRDF.C
1085
1086 Create a file called `EfficiencyAnalysisRDF.C` and paste the code below in. Remember to change the file path.
1087
1088 Execute this script via - `root -l -q EfficiencyAnalysisRDF.C++`. Do this within eic-shell or somewhere else with the correct EDM4hep/EDM4eic libraries installed.
1089
1090 ```c++
1091 #include <edm4hep/utils/vector_utils.h>
1092 #include <edm4hep/MCParticle.h>
1093 #include <edm4eic/ReconstructedParticle.h>
1094 #include <ROOT/RDataFrame.hxx>
1095 #include <ROOT/RVec.hxx>
1096 #include <TFile.h>
1097
1098 // Define aliases for the data types
1099 using MCP = edm4hep::MCParticleData;
1100 using RecoP = edm4eic::ReconstructedParticleData;
1101
1102 // Define function to vectorize the edm4hep::utils methods
1103 template <typename T>
1104 auto getEta = [](ROOT::VecOps::RVec<T> momenta) {
1105   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::eta(p.momentum); });
1106 };
1107
1108 template <typename T>
1109 auto getPhi = [](ROOT::VecOps::RVec<T> momenta) {
1110   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::angleAzimuthal(p.momentum); });
1111 };
1112
1113 // Define the function to perform the efficiency analysis
1114 void EfficiencyAnalysisRDF(TString infile="PATH_TO_FILE"){
1115
1116   // Set up input file
1117   ROOT::RDataFrame df("events", infile);
1118
1119   // Define new dataframe node with additional columns
1120   auto df1 =  df.Define("statusFilter",  "MCParticles.generatorStatus == 1"    )
1121                 .Define("absPDG",        "abs(MCParticles.PDG)"                )
1122                 .Define("pdgFilter",     "absPDG == 11 || absPDG == 13 || absPDG == 211 || absPDG == 321 || absPDG == 2212")
1123                 .Define("particleFilter","statusFilter && pdgFilter"           )
1124                 .Define("filtMCParts",   "MCParticles[particleFilter]"         )
1125                 .Define("assoFilter",    "Take(particleFilter,ReconstructedChargedParticleAssociations.simID)") // Incase any of the associated particles happen to not be charged
1126                 .Define("assoMCParts",   "Take(MCParticles,ReconstructedChargedParticleAssociations.simID)[assoFilter]")
1127                 .Define("assoRecParts",  "Take(ReconstructedChargedParticles,ReconstructedChargedParticleAssociations.recID)[assoFilter]")
1128                 .Define("filtMCEta",     getEta<MCP>   , {"filtMCParts"} )
1129                 .Define("filtMCPhi",     getPhi<MCP>   , {"filtMCParts"} )
1130                 .Define("accoMCEta",     getEta<MCP>   , {"assoMCParts"} )
1131                 .Define("accoMCPhi",     getPhi<MCP>   , {"assoMCParts"} )
1132                 .Define("assoRecEta",    getEta<RecoP> , {"assoRecParts"})
1133                 .Define("assoRecPhi",    getPhi<RecoP> , {"assoRecParts"})
1134                 .Define("deltaR",        "ROOT::VecOps::DeltaR(assoRecEta, accoMCEta, assoRecPhi, accoMCPhi)");
1135
1136   // Define histograms
1137   auto partEta                = df1.Histo1D({"partEta","Eta of Thrown Charged Particles;Eta",100,-5.,5.},"filtMCEta");
1138   auto matchedPartEta         = df1.Histo1D({"matchedPartEta","Eta of Thrown Charged Particles That Have Matching Track",100,-5.,5.},"accoMCEta");
1139   auto matchedPartTrackDeltaR = df1.Histo1D({"matchedPartTrackDeltaR","Delta R Between Matching Thrown and Reconstructed Charged Particle",5000,0.,5.},"deltaR");
1140
1141   // Write histograms to file
1142   TFile *ofile = TFile::Open("EfficiencyAnalysis_Out_RDF.root","RECREATE");
1143
1144   // Booked Define and Histo1D lazy actions are only performed here
1145   partEta->Write();
1146   matchedPartEta->Write();
1147   matchedPartTrackDeltaR->Write();
1148
1149   ofile->Close(); // Close output file
1150 }
1151 ```
1152 <!--
1153
1154 A "solution" using RDataFrames is included below, please see the notes following this script for some of my thoughts on RDataFrames -
1155
1156 ```c++
1157 #include <edm4hep/utils/vector_utils.h>
1158 #include <edm4hep/MCParticle.h>
1159 #include <edm4eic/ReconstructedParticle.h>
1160 #include <ROOT/RDataFrame.hxx>
1161 #include <ROOT/RVec.hxx>
1162 #include <TFile.h>
1163
1164 // Define aliases for the data types
1165 using MCP = edm4hep::MCParticleData;
1166 using RecoP = edm4eic::ReconstructedParticleData;
1167
1168 // Define function to vectorize the edm4hep::utils methods
1169 template <typename T>
1170 auto getEta = [](ROOT::VecOps::RVec<T> momenta) {
1171   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::eta(p.momentum); });
1172 };
1173
1174 template <typename T>
1175 auto getPhi = [](ROOT::VecOps::RVec<T> momenta) {
1176   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::angleAzimuthal(p.momentum); });
1177 };
1178
1179 template <typename T>
1180 auto getP = [](ROOT::VecOps::RVec<T> momenta) {
1181   //return ROOT::VecOps::Map(momenta, [](const T& p) { return (p.momentum); }); // This is a vector3f
1182   return ROOT::VecOps::Map(momenta, [](const T& p) { return edm4hep::utils::magnitude(p.momentum); }); // This is a the magnitude of that vector3f
1183 };
1184
1185 // Define the function to perform the efficiency analysis
1186 void EfficiencyAnalysisRDF_Exercise(TString infile="PATH_TO_INPUT_FILE"){
1187
1188   // Set up input file
1189   ROOT::RDataFrame df("events", infile);
1190
1191   // Define new dataframe node with additional columns
1192   auto df1 =  df.Define("statusFilter",  "MCParticles.generatorStatus == 1"    )
1193     .Define("absPDG",        "abs(MCParticles.PDG)"                )
1194     .Define("pdgFilter",     "absPDG == 11 || absPDG == 13 || absPDG == 211 || absPDG == 321 || absPDG == 2212")
1195     .Define("particleFilter","statusFilter && pdgFilter"           )
1196     .Define("filtMCParts",   "MCParticles[particleFilter]"         )
1197     .Define("assoFilter",    "Take(particleFilter,ReconstructedChargedParticleAssociations.simID)") // Incase any of the associated particles happen to not be charged
1198     .Define("assoMCParts",   "Take(MCParticles,ReconstructedChargedParticleAssociations.simID)[assoFilter]")
1199     .Define("assoRecParts",  "Take(ReconstructedChargedParticles,ReconstructedChargedParticleAssociations.recID)[assoFilter]")
1200     .Define("filtMCEta",     getEta<MCP>   , {"filtMCParts"} )
1201     .Define("filtMCPhi",     getPhi<MCP>   , {"filtMCParts"} )
1202     .Define("filtMCp",       getP<MCP>     , {"filtMCParts"} )
1203     .Define("assoMCEta",     getEta<MCP>   , {"assoMCParts"} )
1204     .Define("assoMCPhi",     getPhi<MCP>   , {"assoMCParts"} )
1205     .Define("assoMCp",       getP<MCP>     , {"assoMCParts"} )
1206     .Define("assoRecEta",    getEta<RecoP> , {"assoRecParts"})
1207     .Define("assoRecPhi",    getPhi<RecoP> , {"assoRecParts"})
1208     .Define("assoRecp",      getP<RecoP>   , {"assoRecParts"})
1209     .Define("deltaEta",      "assoMCEta - assoRecEta"        )
1210     .Define("deltaPhi",      "ROOT::VecOps::DeltaPhi(assoRecPhi, assoMCPhi)")
1211     .Define("deltaR",        "ROOT::VecOps::DeltaR(assoRecEta, assoMCEta, assoRecPhi, assoMCPhi)")
1212     .Define("deltaMom",      "assoMCp - assoRecp")
1213     .Define("recoEta",       getEta<RecoP>,  {"ReconstructedChargedParticles"})
1214     .Define("recoPhi",       getPhi<RecoP>,  {"ReconstructedChargedParticles"})
1215     .Define("recoP",         getP<RecoP>,    {"ReconstructedChargedParticles"});
1216
1217   // 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
1218   auto partEta                = df1.Histo1D({"partEta","#eta of Thrown Charged Particles;#eta",120,-6.,6.},"filtMCEta");
1219   auto matchedPartEta         = df1.Histo1D({"matchedPartEta","#eta of Thrown Charged Particles That Have Matching Track;#eta",120,-6.,6.},"assoMCEta");
1220   auto partMom = df1.Histo1D({"partMom", "Momentum of Thrown Charged Particles (truth); P(GeV/c)", 150, 0, 150}, "filtMCp");
1221   auto matchedPartMom = df1.Histo1D({"matchedPartMom", "Momentum of Thrown Charged Particles (truth), with matching track; P(GeV/c)", 150, 0, 150}, "assoMCp");
1222   auto partPhi = df1.Histo1D({"partPhi", "#phi of Thrown Charged Particles (truth); #phi(rad)", 320, -3.2, 3.2},"filtMCPhi");
1223   auto matchedPartPhi = df1.Histo1D({"matchedPartPhi", "#phi of Thrown Charged Particles (truth), with matching track; #phi(rad)", 320, -3.2, 3.2}, "assoMCPhi");
1224
1225   auto partPEta = df1.Histo2D({"partPEta", "P vs #eta of Thrown Charged Particles; P(GeV/c); #eta", 150, 0, 150, 120, -6, 6}, "filtMCp", "filtMCEta");
1226   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");
1227   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");
1228   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");
1229
1230   auto matchedPartTrackDeltaEta = df1.Histo1D({"matchedPartTrackDeltaEta","#Delta#eta Between Matching Thrown and Reconstructed Charged Particle; #Delta#eta", 100, -0.25, 0.25}, "deltaEta");
1231   auto matchedPartTrackDeltaPhi = df1.Histo1D({"matchedPartTrackDeltaPhi","#Detla #phi Between Matching Thrown and Reconstructed Charged Particle; #Delta#phi", 200, -0.2, 0.2}, "deltaPhi");
1232   auto matchedPartTrackDeltaR = df1.Histo1D({"matchedPartTrackDeltaR","#Delta R Between Matching Thrown and Reconstructed Charged Particle;#Delta R",300,0.,0.3}, "deltaR");
1233   auto matchedPartTrackDeltaMom = df1.Histo1D({"matchedPartTrackDeltaMom","#Delta P Between Matching Thrown and Reconstructed Charged Particle; #Delta P", 200, -10, 10}, "deltaMom");
1234
1235   // Define some histograms for our efficiencies - Done "old school" root style - Maybe the division can be done direct from a DF?
1236   TH1D *TrackEff_Eta = new TH1D("TrackEff_Eta", "Tracking efficiency as fn of #eta; #eta; Eff(%)", 120, -6, 6);
1237   TH1D *TrackEff_Mom = new TH1D("TrackEff_Mom", "Tracking efficiency as fn of P; P(GeV/c); Eff(%)", 150, 0, 150);
1238   TH1D *TrackEff_Phi = new TH1D("TrackEff_Phi", "Tracking efficiency as fn of #phi; #phi(rad); Eff(%)", 320, -3.2, 3.2);
1239   // 2D Efficiencies
1240   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);
1241   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);
1242
1243   auto ChargedEta = df1.Histo1D({"ChargedEta", "#eta of all charged particles; #eta", 120, -6, 6}, "recoEta");
1244   auto ChargedPhi = df1.Histo1D({"ChargedPhi", "#phi of all charged particles; #phi (rad)", 120, -3.2, 3.2}, "recoPhi");
1245   auto ChargedP = df1.Histo1D({"ChargedP", "P of all charged particles; P(GeV/c)", 150, 0, 150}, "recoP");
1246
1247   // Write histograms to file
1248   TFile *ofile = TFile::Open("EfficiencyAnalysis_Exercise_Out_RDF.root","RECREATE");
1249
1250   // Booked Define and Histo1D lazy actions are only performed here
1251   partEta->Write();
1252   matchedPartEta->Write();
1253   partPhi->Write();
1254   matchedPartPhi->Write();
1255   partMom->Write();
1256   matchedPartMom->Write();
1257   partPEta->Write();
1258   matchedPartPEta->Write();
1259   partPhiEta->Write();
1260   matchedPartPhiEta->Write();
1261   matchedPartTrackDeltaEta->Write();
1262   matchedPartTrackDeltaPhi->Write();
1263   matchedPartTrackDeltaR->Write();
1264   matchedPartTrackDeltaMom->Write();
1265
1266   // Create efficiency histograms by dividing appropriately. Note we must actually get the pointer explicitly.
1267   TrackEff_Eta->Divide(matchedPartEta.GetPtr(), partEta.GetPtr(), 1, 1, "b");
1268   TrackEff_Mom->Divide(matchedPartMom.GetPtr(), partMom.GetPtr(), 1, 1, "b");
1269   TrackEff_Phi->Divide(matchedPartPhi.GetPtr(), partPhi.GetPtr(), 1, 1, "b");
1270   TrackEff_PEta->Divide(matchedPartPEta.GetPtr(), partPEta.GetPtr(), 1, 1, "b");
1271   TrackEff_PhiEta->Divide(matchedPartPhiEta.GetPtr(), partPhiEta.GetPtr(), 1, 1, "b");
1272
1273   TrackEff_Eta->Write();
1274   TrackEff_Mom->Write();
1275   TrackEff_Phi->Write();
1276   TrackEff_PEta->Write();
1277   TrackEff_PhiEta->Write();
1278
1279   ChargedEta->Write();
1280   ChargedPhi->Write();
1281   ChargedP->Write();
1282
1283   ofile->Close(); // Close output file
1284 }
1285 ```
1286 Insert your input file path and execute as the example code above.
1287
1288 > Note:
1289 > - Before writing this example and solution, I have never used RDataFrames before. Some thoughts below.
1290 >   - I find them very difficult to work with. The processes that are "actually" going on are very obscured with RDataFrames.
1291 >   - Finding resources and examples online is much more difficult due to how new RDataFrames are.
1292 >     - This further complicates writing and working with them.
1293 >   - 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.
1294 >   - The script itself seems to be very slow when running.
1295 >   - 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.
1296 >  - In summary, I personally would *not* recommend using RDataFrames for your analysis (at this point in time).
1297 >  - Note that this just my (SKay) thoughts and opinions though.
1298 {: .callout}
1299
1300 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.
1301 -->
1302
1303 {% include links.md %}