Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-11-16 09:02:22 

0001 // Code to draw Basic Tracking Performances
0002 // Shyam Kumar:INFN Bari, shyam.kumar@ba.infn.it; shyam055119@gmail.com
0003 
0004 #include <TGraphErrors.h>
0005 #include <TF1.h>
0006 #include <TRandom.h>
0007 #include <TCanvas.h>
0008 #include <TLegend.h>
0009 #include <TMath.h>
0010 
0011 void draw_Performance(Int_t nevent=-1)
0012 {
0013 
0014 //==========Style of the plot============
0015    gStyle->SetPalette(1);
0016    gStyle->SetOptTitle(1);
0017    gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
0018    gStyle->SetTitleSize(.04,"X");gStyle->SetTitleSize(.04,"Y");
0019    gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
0020    gStyle->SetHistLineWidth(2);
0021    gStyle->SetOptFit(1);
0022    gStyle->SetOptStat(0);
0023 
0024 //=======Reading the root file DD4HEP===========
0025    TFile* file = new TFile("tracking_test_gun.edm4eic.root"); // Tree with tracks and hits
0026                 // Create the tree reader and its data containers
0027    TTreeReader myReader("events", file); // name of tree and file
0028 
0029    TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge");
0030    TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x");
0031    TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y");
0032    TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
0033    TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x");
0034    TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y");
0035    TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z");
0036    TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus");
0037    TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");
0038 
0039    TTreeReaderArray<Float_t> px_rec(myReader, "ReconstructedChargedParticles.momentum.x");
0040    TTreeReaderArray<Float_t> py_rec(myReader, "ReconstructedChargedParticles.momentum.y");
0041    TTreeReaderArray<Float_t> pz_rec(myReader, "ReconstructedChargedParticles.momentum.z");
0042 
0043  const int ngraph = 7;
0044  TCanvas * c[ngraph];
0045  for (int i =0; i<ngraph; ++i){
0046  c[i] = new TCanvas(Form("c%d",i),Form("c%d",i),1200,1000);
0047  c[i]->SetMargin(0.09, 0.1 ,0.1,0.06);
0048  }
0049 
0050  // X-Y Hits
0051 
0052  Int_t nbins = 200;
0053  Double_t x= 4.0;
0054  TH2D *hetavspmc = new TH2D("hetavspmc","hetavspmc",400,0.,40.,nbins,-x,x);
0055  TH2D *hetavsprec = new TH2D("hetavsprec","hetavsprec",400,0.,40.,nbins,-x,x);
0056  TH2D *heff_pnum = new TH2D("heff_pnum","heff_p",400,0.,40.,nbins,-x,x);
0057  TH2D *heff_pden = new TH2D("heff_pden","heff_p",400,0.,40.,nbins,-x,x);
0058 
0059  TH2D *hetavsptmc = new TH2D("hetavsptmc","hetavsptmc",400,0.,40.,nbins,-x,x);
0060  TH2D *hetavsptrec = new TH2D("hetavsptrec","hetavsptrec",400,0.,40.,nbins,-x,x);
0061  TH2D *heff_ptnum = new TH2D("heff_ptnum","heff_p",400,0.,40.,nbins,-x,x);
0062  TH2D *heff_ptden = new TH2D("heff_ptden","heff_p",400,0.,40.,nbins,-x,x);
0063 
0064  TH1D *hpresol = new TH1D("hpresol","hpresol",200,-0.5,0.5);
0065 
0066  hetavspmc->GetXaxis()->SetTitle("p_{mc} (GeV/c)");
0067  hetavspmc->GetYaxis()->SetTitle("#eta_{mc}");
0068  hetavspmc->GetXaxis()->CenterTitle();
0069  hetavspmc->GetYaxis()->CenterTitle();
0070 
0071  hetavsptmc->GetXaxis()->SetTitle("pt_{mc} (GeV/c)");
0072  hetavsptmc->GetYaxis()->SetTitle("#eta_{mc}");
0073  hetavsptmc->GetXaxis()->CenterTitle();
0074  hetavsptmc->GetYaxis()->CenterTitle();
0075 
0076  hetavsprec->GetXaxis()->SetTitle("p_{rec} (GeV/c)");
0077  hetavsprec->GetYaxis()->SetTitle("#eta_{rec}");
0078  hetavsprec->GetXaxis()->CenterTitle();
0079  hetavsprec->GetYaxis()->CenterTitle();
0080 
0081  hetavsptrec->GetXaxis()->SetTitle("pt_{rec} (GeV/c)");
0082  hetavsptrec->GetYaxis()->SetTitle("#eta_{rec}");
0083  hetavsptrec->GetXaxis()->CenterTitle();
0084  hetavsptrec->GetYaxis()->CenterTitle();
0085 
0086  hpresol->GetXaxis()->SetTitle("dp/p");
0087  hpresol->GetYaxis()->SetTitle("Entries (a.u.)");
0088  hpresol->GetXaxis()->CenterTitle();
0089  hpresol->GetYaxis()->CenterTitle();
0090 
0091 
0092 //////////////////////////////////////////////////////////////////////
0093  int count = 0;
0094   while (myReader.Next())
0095   {
0096    if (nevent>0 && count>nevent) continue;
0097  //  cout<<"=====Event No. "<<count<<"============="<<endl;
0098      Double_t pmc = 0;  Double_t etamc = 0; Double_t ptmc = 0;
0099    // MC Particle
0100      for (int iParticle = 0; iParticle < charge.GetSize(); ++iParticle){
0101    //  cout<<" PDG: "<<pdg[iParticle]<<" Status: "<<status[iParticle]<<" Pt: "<<sqrt(px_mc[iParticle]*px_mc[iParticle]+py_mc[iParticle]*py_mc[iParticle])<<endl;
0102      if (status[iParticle] ==1){
0103      pmc = sqrt(px_mc[iParticle]*px_mc[iParticle]+py_mc[iParticle]*py_mc[iParticle]+pz_mc[iParticle]*pz_mc[iParticle]);
0104      Double_t ptmc = sqrt(px_mc[iParticle]*px_mc[iParticle]+py_mc[iParticle]*py_mc[iParticle]);
0105      Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pz_mc[iParticle]/pmc))/2));
0106      hetavspmc->Fill(pmc,etamc);
0107      hetavsptmc->Fill(ptmc,etamc);
0108      heff_pden->Fill(pmc,etamc);
0109      heff_ptden->Fill(ptmc,etamc);
0110      if (px_rec.GetSize()==1){
0111      Double_t prec = sqrt(px_rec[iParticle]*px_rec[iParticle]+py_rec[iParticle]*py_rec[iParticle]+pz_rec[iParticle]*pz_rec[iParticle]);
0112      Double_t ptrec = sqrt(px_rec[iParticle]*px_rec[iParticle]+py_rec[iParticle]*py_rec[iParticle]);
0113      Double_t etarec = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pz_rec[iParticle]/prec))/2));
0114      hetavsprec->Fill(prec,etarec);
0115      hetavsptrec->Fill(ptrec,etarec);
0116      hpresol->Fill((prec-pmc)/pmc);
0117      heff_pnum->Fill(pmc,etamc);
0118      heff_ptnum->Fill(ptmc,etamc);
0119      }
0120      }
0121      }
0122 
0123     count++;
0124 
0125   }
0126 
0127       heff_pnum->Divide(heff_pden);
0128       heff_ptnum->Divide(heff_ptden);
0129       heff_pnum->GetYaxis()->SetTitle("Acceptance");
0130       heff_pnum->GetXaxis()->SetTitle("p_mc");
0131       heff_pnum->GetZaxis()->SetRangeUser(0.,1.1);
0132 
0133       heff_ptnum->GetYaxis()->SetTitle("Acceptance");
0134       heff_ptnum->GetXaxis()->SetTitle("pt_mc");
0135       heff_ptnum->GetZaxis()->SetRangeUser(0.,1.1);
0136 
0137      c[0]->cd();
0138      hetavspmc->Draw("colz");
0139      c[0]->SaveAs("eta_mcvspmc.png");
0140      c[1]->cd();
0141      hetavsprec->Draw("colz");
0142      c[1]->SaveAs("eta_mcvsprec.png");
0143      c[2]->cd();
0144      hetavsptmc->Draw("colz");
0145      c[2]->SaveAs("eta_mcvsptmc.png");
0146      c[3]->cd();
0147      hetavsptrec->Draw("colz");
0148      c[3]->SaveAs("eta_mcvsptrec.png");
0149      c[4]->cd();
0150      hpresol->Draw("hist");
0151      c[4]->SaveAs("ptresol.png");
0152      c[5]->cd();
0153      heff_pnum->Draw("colz");
0154      c[5]->SaveAs("effvsp_2D.png");
0155      c[6]->cd();
0156      heff_ptnum->Draw("colz");
0157      c[6]->SaveAs("effvspt_2D.png");
0158 
0159 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!