tracking_performances_dis/analysis/vtx_dis_analysis.cxx

0001 #include <fstream>
0002 #include <iostream>
0003 #include <string>
0004
0005 #include <TChain.h>
0006 #include <TFile.h>
0007 #include <TTreeReader.h>
0008 #include <TTreeReaderArray.h>
0009 #include <TH1D.h>
0010 #include <TH2D.h>
0011 #include <Math/Vector3D.h>
0012 #include <Math/Vector4D.h>
0013 #include <Math/RotationY.h>
0014 #include <TMath.h>
0015 #include <TVector3.h>
0016 #include <TLorentzVector.h>
0017
0018 #include "fmt/color.h"
0019 #include "fmt/core.h"
0020
0021 #include "nlohmann/json.hpp"
0022
0023 void vtx_dis_analysis(const std::string& config_name)
0024 {
0025     // Read our configuration
0026     std::ifstream  config_file{config_name};
0027     nlohmann::json config;
0028     config_file >> config;
0029
0030     const std::string rec_file      = config["rec_file"];
0031     const std::string detector      = config["detector"];
0032     const std::string output_prefix = config["output_prefix"];
0033     const int         ebeam         = config["ebeam"];
0034     const int         pbeam         = config["pbeam"];
0035     const int         Q2_min        = config["Min_Q2"];
0036
0037     fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
0038                 "Running DIS tracking analysis...\n");
0039     fmt::print(" - Detector package: {}\n", detector);
0040     fmt::print(" - input file: {}\n", rec_file);
0041     fmt::print(" - output prefix for histograms: {}\n", output_prefix);
0042     fmt::print(" - ebeam: {}\n", ebeam);
0043     fmt::print(" - pbeam: {}\n", pbeam);
0044     fmt::print(" - Minimum Q2: {}\n", Q2_min);
0045
0046     //--------------------------------------------------------------------------------------------------------------------------------------------
0047
0048     // Set output file for the histograms
0049     std::string output_name_hists = fmt::format("{}.root", output_prefix);
0050     cout << "Output file for histograms = " << output_name_hists << endl;
0051     TFile* ofile = new TFile(output_name_hists.c_str(), "RECREATE");
0052
0053     //--------------------------------------------------------------------------------------------------------------------------------------------
0054
0055     // Set up input file chain
0056     TChain *mychain = new TChain("events");
0057     mychain->Add(rec_file.c_str());
0058
0059     //--------------------------------------------------------------------------------------------------------------------------------------------
0060
0061     // TTreeReader
0062     TTreeReader tree_reader(mychain);
0063
0064     // Reco Vertex
0065     TTreeReaderArray<int> recoVtxType = {tree_reader, "CentralTrackVertices.type"};
0066     TTreeReaderArray<float> recoVtxX = {tree_reader, "CentralTrackVertices.position.x"};
0067     TTreeReaderArray<float> recoVtxY = {tree_reader, "CentralTrackVertices.position.y"};
0068     TTreeReaderArray<float> recoVtxZ = {tree_reader, "CentralTrackVertices.position.z"};
0069
0070     TTreeReaderArray<unsigned int> assoPartBegin = {tree_reader, "CentralTrackVertices.associatedParticles_begin"};
0071     TTreeReaderArray<unsigned int> assoPartEnd = {tree_reader, "CentralTrackVertices.associatedParticles_end"};
0072     TTreeReaderArray<int> assoPartIndex = {tree_reader, "_CentralTrackVertices_associatedParticles.index"};
0073
0074     // MC
0075     TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
0076     TTreeReaderArray<int> mcPDG = {tree_reader, "MCParticles.PDG"};
0077     TTreeReaderArray<float> mcCharge = {tree_reader, "MCParticles.charge"};
0078     TTreeReaderArray<float> mcMomX = {tree_reader, "MCParticles.momentum.x"};
0079     TTreeReaderArray<float> mcMomY = {tree_reader, "MCParticles.momentum.y"};
0080     TTreeReaderArray<float> mcMomZ = {tree_reader, "MCParticles.momentum.z"};
0081
0082     TTreeReaderArray<double> mcVtxX = {tree_reader, "MCParticles.vertex.x"};
0083     TTreeReaderArray<double> mcVtxY = {tree_reader, "MCParticles.vertex.y"};
0084     TTreeReaderArray<double> mcVtxZ = {tree_reader, "MCParticles.vertex.z"};
0085
0086     TTreeReaderArray<unsigned int> mcParentBegin = {tree_reader, "MCParticles.parents_begin"};
0087     TTreeReaderArray<unsigned int> mcParentEnd = {tree_reader, "MCParticles.parents_end"};
0088     TTreeReaderArray<int> mcParentIndex = {tree_reader, "_MCParticles_parents.index"};
0089
0090     // Reco
0091     TTreeReaderArray<int> recoType = {tree_reader, "ReconstructedChargedParticles.type"};
0092
0093
0094     //--------------------------------------------------------------------------------------------------------------------------------------------
0095     // Define Histograms
0096
0097     // Gen
0098     TH1D *numGenTracksHist = new TH1D("numGenTracks","",31,-0.5,30.5);
0099
0100     TH2D *genVtxYvsXHist = new TH2D("genVtxYvsXHist","",200,-1.,1.,200,-1,1);
0101     genVtxYvsXHist->Draw("COLZ");
0102     genVtxYvsXHist->SetTitle("MC Vertex: v_{x} versus v_{y}");
0103     genVtxYvsXHist->GetXaxis()->SetTitle("x-coordinate (in mm)");
0104     genVtxYvsXHist->GetYaxis()->SetTitle("y-coordinate (in mm)");
0105     genVtxYvsXHist->GetXaxis()->CenterTitle(); genVtxYvsXHist->GetYaxis()->CenterTitle();
0106
0107     TH2D *genVtxRvsZHist = new TH2D("genVtxRvsZHist","",200,-100.,100.,100,0,0.5);
0108     genVtxRvsZHist->Draw("COLZ");
0109     genVtxRvsZHist->SetTitle("MC Vertex: v_{r} versus v_{z}");
0110     genVtxRvsZHist->GetXaxis()->SetTitle("z-coordinate (in mm)");
0111     genVtxRvsZHist->GetYaxis()->SetTitle("#sqrt{x^{2} + y^{2}} (in mm)");
0112     genVtxRvsZHist->GetXaxis()->CenterTitle(); genVtxRvsZHist->GetYaxis()->CenterTitle();
0113
0114     // Reco
0115     TH1D *numRecoTracksHist = new TH1D("numRecoTracks","",31,-0.5,30.5);
0116
0117     TH1D *recoVtxEffHist = new TH1D("recoVtxEff","",4,-0.5,3.5);
0118     recoVtxEffHist->Draw("P");
0119     recoVtxEffHist->SetLineColor(kRed); recoVtxEffHist->SetMarkerColor(kRed);
0120     recoVtxEffHist->SetMarkerStyle(8); recoVtxEffHist->SetMarkerSize(1.2);
0121     recoVtxEffHist->SetTitle("Vertexing Efficiency");
0122     recoVtxEffHist->GetXaxis()->SetTitle("Vertex Count"); recoVtxEffHist->GetYaxis()->SetTitle("nEvents/total_events %");
0123     recoVtxEffHist->GetXaxis()->CenterTitle(); recoVtxEffHist->GetYaxis()->CenterTitle();
0124
0125     TH2D *recoVtxYvsXHist = new TH2D("recoVtxYvsX","",200,-1.,1.,200,-1,1);
0126     recoVtxYvsXHist->Draw("COLZ");
0127     recoVtxYvsXHist->SetTitle("Reconstructed Vertex: v_{x} versus v_{y}");
0128     recoVtxYvsXHist->GetXaxis()->SetTitle("x-coordinate (in mm)"); recoVtxYvsXHist->GetYaxis()->SetTitle("y-coordinate (in mm)");
0129     recoVtxYvsXHist->GetXaxis()->CenterTitle(); recoVtxYvsXHist->GetYaxis()->CenterTitle();
0130
0131     TH2D *recoVtxRvsZHist = new TH2D("recoVtxRvsZ","",200,-100.,100.,100,0.,0.8);
0132     recoVtxRvsZHist->Draw("COLZ");
0133     recoVtxRvsZHist->SetTitle("Reconstructed Vertex: v_{r} versus v_{z}");
0134     recoVtxRvsZHist->GetXaxis()->SetTitle("z-coordinate (in mm)");
0135     recoVtxRvsZHist->GetYaxis()->SetTitle("#sqrt{x^{2} + y^{2}} (in mm)");
0136     recoVtxRvsZHist->GetXaxis()->CenterTitle(); recoVtxRvsZHist->GetYaxis()->CenterTitle();
0137
0138     //Reco Vs MC
0139     TH2D *recoVsMCTracksHist = new TH2D("recoVsMCTracks","",31,-0.5,30.5,31,-0.5,30.5);
0140     recoVsMCTracksHist->Draw("COLZ");
0141     recoVsMCTracksHist->SetTitle("Number of particles associated with vertex");
0142     recoVsMCTracksHist->GetXaxis()->SetTitle("N_{MC}"); recoVsMCTracksHist->GetYaxis()->SetTitle("N_{RC}");
0143     recoVsMCTracksHist->GetXaxis()->CenterTitle(); recoVsMCTracksHist->GetYaxis()->CenterTitle();
0144
0145     // Resolution
0146     TH2D *vtxResXvsGenTrkHist = new TH2D("vtxResXvsGenTrk","",31,-0.5,30.5,200,-1,1);
0147     vtxResXvsGenTrkHist->Draw("COLZ");
0148     vtxResXvsGenTrkHist->SetTitle("Vertex Resolution X vs MC Tracks");
0149     vtxResXvsGenTrkHist->GetXaxis()->SetTitle("N_{MC}"); vtxResXvsGenTrkHist->GetYaxis()->SetTitle("recVtx_x - mcVtx_x (in mm)");
0150     vtxResXvsGenTrkHist->GetXaxis()->CenterTitle(); vtxResXvsGenTrkHist->GetYaxis()->CenterTitle();
0151
0152     TH2D *vtxResYvsGenTrkHist = new TH2D("vtxResYvsGenTrk","",31,-0.5,30.5,200,-1,1);
0153     vtxResYvsGenTrkHist->Draw("COLZ");
0154     vtxResYvsGenTrkHist->SetTitle("Vertex Resolution Y vs MC Tracks");
0155     vtxResYvsGenTrkHist->GetXaxis()->SetTitle("N_{MC}");
0156     vtxResYvsGenTrkHist->GetYaxis()->SetTitle("recVtx_y - mcVtx_y (in mm)");
0157     vtxResYvsGenTrkHist->GetXaxis()->CenterTitle(); vtxResYvsGenTrkHist->GetYaxis()->CenterTitle();
0158
0159     TH2D *vtxResZvsGenTrkHist = new TH2D("vtxResZvsGenTrk","",31,-0.5,30.5,200,-1,1);
0160     vtxResZvsGenTrkHist->Draw("COLZ");
0161     vtxResZvsGenTrkHist->SetTitle("Vertex Resolution Z vs MC Tracks");
0162     vtxResZvsGenTrkHist->GetXaxis()->SetTitle("N_{MC}");
0163     vtxResZvsGenTrkHist->GetYaxis()->SetTitle("recVtx_z - mcVtx_z (in mm)");
0164     vtxResZvsGenTrkHist->GetXaxis()->CenterTitle(); vtxResZvsGenTrkHist->GetYaxis()->CenterTitle();
0165
0166     TH2D *vtxResXvsRecoTrkHist = new TH2D("vtxResXvsRecoTrk","",31,-0.5,30.5,200,-1,1);
0167     vtxResXvsRecoTrkHist->Draw("COLZ");
0168     vtxResXvsRecoTrkHist->SetTitle("Vertex Resolution X vs Reconstructed Tracks");
0169     vtxResXvsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}");
0170     vtxResXvsRecoTrkHist->GetYaxis()->SetTitle("recVtx_x - mcVtx_x (in mm)");
0171     vtxResXvsRecoTrkHist->GetXaxis()->CenterTitle(); vtxResXvsRecoTrkHist->GetYaxis()->CenterTitle();
0172
0173     TH2D *vtxResYvsRecoTrkHist = new TH2D("vtxResYvsRecoTrk","",31,-0.5,30.5,200,-1,1);
0174     vtxResYvsRecoTrkHist->Draw("COLZ");
0175     vtxResYvsRecoTrkHist->SetTitle("Vertex Resolution Y vs Reconstructed Tracks");
0176     vtxResYvsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}");
0177     vtxResYvsRecoTrkHist->GetYaxis()->SetTitle("recVtx_y - mcVtx_y (in mm)");
0178     vtxResYvsRecoTrkHist->GetXaxis()->CenterTitle(); vtxResYvsRecoTrkHist->GetYaxis()->CenterTitle();
0179
0180     TH2D *vtxResZvsRecoTrkHist = new TH2D("vtxResZvsRecoTrk","",31,-0.5,30.5,200,-1,1);
0181     vtxResZvsRecoTrkHist->Draw("COLZ");
0182     vtxResZvsRecoTrkHist->SetTitle("Vertex Resolution Z vs Reconstructed Tracks");
0183     vtxResZvsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}");
0184     vtxResZvsRecoTrkHist->GetYaxis()->SetTitle("recVtx_z - mcVtx_z (in mm)");
0185     vtxResZvsRecoTrkHist->GetXaxis()->CenterTitle(); vtxResZvsRecoTrkHist->GetYaxis()->CenterTitle();
0186
0187     TH1D *numGenTrkswithVtxHist = new TH1D("numGenTrkswithVtx","",31,-0.5,30.5);
0188     TH1D *numRecoTrkswithVtxHist = new TH1D("numRecoTrkswithVtx","",31,-0.5,30.5);
0189
0190
0191     int counter(0);
0192
0193     //Loop over events
0194     std::cout<<"Analyzing "<<mychain->GetEntries()<<" events!"<<std::endl;
0195     while (tree_reader.Next()) {
0196
0197         if(counter%100==0) std::cout<<"Analyzing event "<<counter<<std::endl;
0198
0199         counter++;
0200
0201         //////////////////////////////////////////////////////////////////////////
0202         ////////////////////////  Analyze MC Tracks  /////////////////////////
0203         //////////////////////////////////////////////////////////////////////////
0204
0205         //Finding MC vertex using scattered electron
0206         TVector3 mcEvtVtx(-999., -999., -999.);
0207         for(unsigned int i=0; i<mcGenStat.GetSize(); i++)
0208         {
0209         if(mcGenStat[i] != 1) continue;
0210         if(mcPDG[i] != 11) continue;
0211
0212         bool scatEfound = false;
0213         // mcParentBegin and mcParentEnd specify the entries from _MCParticles_parents.index
0214         // _MCParticles_parents.index stores the MCParticle index
0215         for(unsigned int j=mcParentBegin[i]; j<mcParentEnd[i]; j++)
0216         {
0217             int parentPDG = mcPDG[mcParentIndex[j]];
0218             if(parentPDG == 11) scatEfound = true;
0219             }
0220
0221         if(scatEfound == false) continue;
0222         //Scattered electron found
0223         double vtx_mc_x =  mcVtxX[i];
0224         double vtx_mc_y =  mcVtxY[i];
0225         double vtx_mc_z =  mcVtxZ[i];
0226         mcEvtVtx = TVector3(vtx_mc_x, vtx_mc_y, vtx_mc_z);
0227         }
0228         genVtxYvsXHist->Fill(mcEvtVtx.x(), mcEvtVtx.y());
0229         TVector3 mcRadius(mcEvtVtx.x(), mcEvtVtx.y(), 0);
0230         genVtxRvsZHist->Fill(mcEvtVtx.z(), mcRadius.Mag());
0231
0232         //Filtering MC Tracks
0233         int numMCTracks=0;
0234         for(unsigned int i=0; i<mcGenStat.GetSize(); i++)
0235         {
0236         if(mcGenStat[i] != 1) continue;
0237         if(mcCharge[i] == 0) continue;
0238
0239         TVector3 mcPartVtx(mcVtxX[i], mcVtxY[i], mcVtxZ[i]);
0240         TVector3 vtx_diff = mcPartVtx - mcEvtVtx;
0241         if(vtx_diff.Mag() > 1e-4) continue;
0242
0243         TVector3 mcPartMom(mcMomX[i], mcMomY[i], mcMomZ[i]);
0244         if(fabs(mcPartMom.Eta()) > 3.5) continue;
0245
0246         numMCTracks++;
0247         }
0248         numGenTracksHist->Fill(numMCTracks);
0249
0250         //////////////////////////////////////////////////////////////////////////
0251         //////////////////////  Analyze Reconstructed Tracks  //////////////////////
0252         //////////////////////////////////////////////////////////////////////////
0253
0254         numRecoTracksHist->Fill(recoType.GetSize());
0255
0256         //Finding Reconstructed Vertex and Vertexing Efficiency
0257         int nVtx=0;
0258         float diff=999.;
0259         int nAssoPart=0;
0260         TVector3 recoEvtVtx(-999., -999., -999.);
0261         for(unsigned int i=0; i<recoVtxType.GetSize(); i++)
0262         {
0263         nVtx++;
0264
0265         TVector3 recoVtx(recoVtxX[i], recoVtxY[i], recoVtxZ[i]);
0266
0267         //Finding the reconstructed vertex closer to the MC vertex
0268         TVector3 vtx_diff = recoVtx - mcEvtVtx;
0269         if(vtx_diff.Mag() < diff)
0270         {
0271             diff = vtx_diff.Mag();
0272             recoEvtVtx = recoVtx;
0273
0274             for(unsigned int j=assoPartBegin[i]; j<assoPartEnd[i]; j++)
0275             {
0276                 nAssoPart = j;
0277                 }
0278         }
0279         }
0280
0281         recoVtxEffHist->Fill(nVtx);
0282         recoVtxYvsXHist->Fill(recoEvtVtx.x(), recoEvtVtx.y());
0283         TVector3 recoRadius(recoEvtVtx.x(), recoEvtVtx.y(), 0);
0284         recoVtxRvsZHist->Fill(recoEvtVtx.z(), recoRadius.Mag());
0285
0286         vtxResXvsGenTrkHist->Fill(numMCTracks, recoEvtVtx.x() - mcEvtVtx.x());
0287         vtxResYvsGenTrkHist->Fill(numMCTracks, recoEvtVtx.y() - mcEvtVtx.y());
0288         vtxResZvsGenTrkHist->Fill(numMCTracks, recoEvtVtx.z() - mcEvtVtx.z());
0289
0290         vtxResXvsRecoTrkHist->Fill(nAssoPart, recoEvtVtx.x() - mcEvtVtx.x());
0291         vtxResYvsRecoTrkHist->Fill(nAssoPart, recoEvtVtx.y() - mcEvtVtx.y());
0292         vtxResZvsRecoTrkHist->Fill(nAssoPart, recoEvtVtx.z() - mcEvtVtx.z());
0293
0294         if(nVtx !=0) {
0295         numGenTrkswithVtxHist->Fill(numMCTracks);
0296         numRecoTrkswithVtxHist->Fill(recoType.GetSize());
0297
0298         recoVsMCTracksHist->Fill(numMCTracks, nAssoPart);}
0299
0300     }
0301
0302     //--------------------------------------------------------------------------------------------------------------------------------------------
0303     ofile->Write(); // Write histograms to file
0304     ofile->Close(); // Close output file
0305
0306 }
0307