Lcp_Reconstruction/Pull_distributions/CheckPull_Vertex.C

0001 // Code to draw the pull distribution of the primary vertex
0002 // Shyam Kumar; INFN Bari; shyam.kumar@ba.infn.it; shyam055119@gmail.com
0003
0004 void CheckPull_Vertex(TString filelist = "test_D0.list"){
0005
0006     gStyle->SetPalette(kRainBow);
0007     gStyle->SetTitleSize(0.045,"XY");
0008     gStyle->SetTitleSize(0.04,"XY");
0009     gStyle->SetLabelSize(0.04,"XY");
0010     gStyle->SetTitleOffset(1.0,"XY");
0011     gStyle->SetOptStat(1);
0012     gStyle->SetOptFit(1);
0013     gStyle->SetOptTitle(0);
0014     gStyle->SetGridColor(kBlack);
0015     gStyle->SetGridWidth(2);
0016     gStyle->SetGridStyle(2);
0017
0018     TH1F *hPullVtxX = new TH1F("hPullVtxX", "Pull x distribution of Primary vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}; Entries (a.u.)", 1000, -10., 10.);
0019      TH1F *hPullVtxY = new TH1F("hPullVtxY", "Pull y distribution of Primary vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}; Entries (a.u.)", 1000, -10., 10.);
0020      TH1F *hPullVtxZ = new TH1F("hPullVtxZ", "Pull z distribution of Primary vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}; Entries (a.u.)", 2000, -20, 20);
0021
0022      TH1F *hPullVtxX_1 = new TH1F("hPullVtxX_1", "Pull x distribution (Ntracks = 1) of Primary vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}; Entries (a.u.)", 1000, -10., 10.);
0023      TH1F *hPullVtxY_1 = new TH1F("hPullVtxY_1", "Pull y distribution (Ntracks = 1) of Primary vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}; Entries (a.u.)", 1000, -10., 10.);
0024      TH1F *hPullVtxZ_1 = new TH1F("hPullVtxZ_1", "Pull z distribution (Ntracks = 1) of Primary vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}; Entries (a.u.)", 2000, -20, 20);
0025
0026      TH1F *hPullVtxX_2 = new TH1F("hPullVtxX_2", "Pull x distribution (Ntracks > 1) of Primary vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}; Entries (a.u.)", 1000, -10., 10.);
0027      TH1F *hPullVtxY_2 = new TH1F("hPullVtxY_2", "Pull y distribution (Ntracks > 1) of Primary vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}; Entries (a.u.)", 1000, -10., 10.);
0028      TH1F *hPullVtxZ_2 = new TH1F("hPullVtxZ_2", "Pull z distribution (Ntracks > 1) of Primary vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}; Entries (a.u.)", 2000, -20, 20);
0029
0030      TH1F *hResVtxX = new TH1F("hResVtxX", "Residual x distribution (Ntracks > 1) of Primary vertex;(Vx_{rec}-Vx_{mc}) (mm); Entries (a.u.)", 1000, -10., 10.);
0031      TH1F *hResVtxY = new TH1F("hResVtxY", "Residual y distribution (Ntracks > 1) of Primary vertex;(Vy_{rec}-Vy_{mc}) (mm); Entries (a.u.)", 1000, -10., 10.);
0032      TH1F *hResVtxZ = new TH1F("hResVtxZ", "Residual z distribution (Ntracks > 1) of Primary vertex;(Vz_{rec}-Vz_{mc}) (mm); Entries (a.u.)", 2000, -5., 5.);
0033      TH1F *hChi2 = new TH1F("hChi2", "Chi2 distribution (Ntracks > 1) of Primary vertex; #chi^{2}; Entries (a.u.)", 2000, 0., 200.);
0034
0035      TFile *fout = new TFile(Form("Pull_distribution_vertex_%s.root",filelist.Data()),"recreate");
0036
0037
0038      // Read all files
0039     int nfiles = 0;
0040      char filename[512];
0041     TChain *chain = new TChain("events");
0042
0043     ifstream *inputstream = new ifstream;
0044      inputstream->open(filelist.Data());
0045      if(!inputstream) printf("[Error] Cannot open file list: %s\n", filelist.Data());
0046
0047     while (inputstream->good())
0048       {
0049       inputstream->getline(filename, 512);
0050       if (inputstream->good())
0051       {
0052         TFile *ftmp = TFile::Open(filename, "READ");
0053         if (!ftmp || !ftmp->IsOpen() || !ftmp->GetNkeys())
0054         {
0055          printf("[e] Skipping bad file: %s\n", filename);
0056          if (ftmp) { ftmp->Close(); delete ftmp; }
0057          continue;
0058         }
0059         cout << "[i] Add " << nfiles << "th file: " << filename << endl;
0060         chain->Add(filename);
0061         nfiles++;
0062
0063         ftmp->Close();
0064         delete ftmp;
0065       }
0066     }
0067
0068     inputstream->close();
0069     printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0070
0071     TTreeReader treereader(chain);
0072     TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0073      TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0074      TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0075     TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0076     TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0077     TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0078
0079     TTreeReaderArray<float> CTVx = {treereader, "CentralTrackVertices.position.x"};
0080     TTreeReaderArray<float> CTVy = {treereader, "CentralTrackVertices.position.y"};
0081     TTreeReaderArray<float> CTVz = {treereader, "CentralTrackVertices.position.z"};
0082     TTreeReaderArray<float> CTVerr_xx = {treereader, "CentralTrackVertices.positionError.xx"};
0083     TTreeReaderArray<float> CTVerr_yy = {treereader, "CentralTrackVertices.positionError.yy"};
0084     TTreeReaderArray<float> CTVerr_zz = {treereader, "CentralTrackVertices.positionError.zz"};
0085     TTreeReaderArray<float> CTVchi2 = {treereader, "CentralTrackVertices.chi2"};
0086     TTreeReaderArray<int> prim_vtx_index = {treereader, "PrimaryVertices_objIdx.index"};
0087      TTreeReaderArray<float> rcCharge = {treereader, "ReconstructedChargedParticles.charge"};
0088
0089     int nevents = 0;
0090      while(treereader.Next())
0091     {
0092       if(nevents%1000==0) printf("\n[i] Events %d\n",nevents);
0093       // find MC primary vertex
0094       int nMCPart = mcPartMass.GetSize();
0095
0096       TVector3 vertex_mc(-999., -999., -999.);
0097       for(int imc=0; imc<nMCPart; imc++)
0098     {
0099       if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0100         {
0101           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0102           break;
0103         }
0104     }
0105       // get RC primary vertex and it's error
0106      TVector3 vertex_rc(-999., -999., -999.);
0107      TVector3 err_vertex_rc(-999., -999., -999.);
0108      double chi2 = 0.;
0109      if(prim_vtx_index.GetSize()>0)
0110     {
0111       int rc_vtx_index = prim_vtx_index[0];
0112       vertex_rc.SetXYZ(CTVx[rc_vtx_index], CTVy[rc_vtx_index], CTVz[rc_vtx_index]);
0113       err_vertex_rc.SetXYZ(sqrt(CTVerr_xx[rc_vtx_index]), sqrt(CTVerr_yy[rc_vtx_index]), sqrt(CTVerr_zz[rc_vtx_index]));
0114       chi2 = CTVchi2[rc_vtx_index];
0115     }
0116     hPullVtxX->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x());
0117     hPullVtxY->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y());
0118      hPullVtxZ->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());
0119
0120      // Separate two cases one with one track while other more than one track
0121      if (rcCharge.GetSize()==1){
0122     hPullVtxX_1->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x());
0123     hPullVtxY_1->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y());
0124      hPullVtxZ_1->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());
0125      }
0126      else{
0127     hPullVtxX_2->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x());
0128     hPullVtxY_2->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y());
0129      hPullVtxZ_2->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());
0130
0131      hResVtxX->Fill((vertex_rc.x()-vertex_mc.x()));
0132     hResVtxY->Fill((vertex_rc.y()-vertex_mc.y()));
0133      hResVtxZ->Fill((vertex_rc.z()-vertex_mc.z()));
0134      hChi2->Fill(chi2);
0135
0136      }
0137
0138      nevents++;
0139     }
0140
0141
0142     fout->cd();
0143     hPullVtxX->Write();
0144     hPullVtxY->Write();
0145     hPullVtxZ->Write();
0146     hPullVtxX_1->Write();
0147     hPullVtxY_1->Write();
0148     hPullVtxZ_1->Write();
0149     hPullVtxX_2->Write();
0150     hPullVtxY_2->Write();
0151     hPullVtxZ_2->Write();
0152     hResVtxX->Write();
0153     hResVtxY->Write();
0154     hResVtxZ->Write();
0155     hChi2->Write();
0156     fout->Close();
0157     delete inputstream;
0158     delete chain;
0159 }