benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaT.C

0001 // Code to compare the tracking performances: Truth seeding vs real seeding
0002 // Shyam Kumar; 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_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.);
0012 void doCompare_truth_real_widebins_dcaT(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, Bool_t drawreq=1, TString epic ="24.06.0", TString eicrecon = "v1.14.0") // name = p, pt for getting p or pt dependence fitted results
0013 {
0014
0015 //=== style of the plot=========
0016    gStyle->SetPalette(1);
0017    gStyle->SetOptTitle(1);
0018    gStyle->SetTitleOffset(1.0,"XY");
0019    gStyle->SetTitleSize(.04,"XY");
0020    gStyle->SetLabelSize(.04,"XY");
0021    gStyle->SetHistLineWidth(2);
0022    gStyle->SetOptFit(1);
0023    gStyle->SetOptStat(1);
0024
0025    const Int_t nptbins = 10;
0026    double pt[nptbins] ={0.2, 0.3, 0.5,1.0, 1.5, 2.0, 5.0, 8.0, 10., 15.0};
0027    Double_t variation = 0.1; // 10 % variation
0028    std::vector<double> momV_truth, momV_real;
0029    std::vector<double> dcaTresolV_truth, err_dcaTresolV_truth, dcaTresolV_real, err_dcaTresolV_real;
0030    momV_truth.clear(); momV_real.clear(); dcaTresolV_truth.clear(); err_dcaTresolV_truth.clear(); dcaTresolV_real.clear(); err_dcaTresolV_real.clear();
0031
0032
0033    TString symbolname = "";
0034    if (particle == "pi-") symbolname = "#pi^{-}";
0035    else symbolname = particle;
0036    // Write the parameters in text file (for comparison)
0037    ofstream outfile;
0038    outfile.open ("DCAT_resol.txt",ios_base::app);
0039
0040
0041    TF1 *f1=new TF1("f1","FitPointingAngle",0.,30.0,2);
0042    f1->SetParLimits(0,0.,50000);
0043    f1->SetParLimits(1,0.,50000);
0044
0045
0046    TCanvas *c_dcaxy = new TCanvas("c_dcaxy","c_dcaxy",1400,1000);
0047    c_dcaxy->SetMargin(0.10, 0.05 ,0.1,0.05);
0048    c_dcaxy->SetGridy();
0049
0050    //Reading the root file
0051     TFile *fDCA_truth, *fDCA_real;
0052     TGraphErrors *gr_dcaT_truth, *gr_dcaT_real, *gr_dcaZ_truth, *gr_dcaZ_real;
0053
0054     TMultiGraph *mgDCAT;
0055     TLegend *lDCAT;
0056      mgDCAT = new TMultiGraph("mgDCAT",";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)");
0057      lDCAT = new TLegend(0.65,0.80,0.90,0.93);
0058      lDCAT->SetTextSize(0.03);
0059      lDCAT->SetBorderSize(0);
0060      lDCAT->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C");
0061
0062       fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data()));
0063       fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data()));
0064
0065      // Truth seeding histograms
0066       TH3D *hist_d0xy_truth = (TH3D*) fDCA_truth->Get("h_d0xy_3d");
0067       TH3D *hist_d0xy_real = (TH3D*) fDCA_real->Get("h_d0xy_3d");
0068
0069       // d0xy calculation for truth/real (binning are same in both cases)
0070      Int_t etamin_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamin+0.0001);
0071      Int_t etamax_bin = hist_d0xy_truth->GetYaxis()->FindBin(etamax-0.0001);
0072
0073      TF1 *func_truth = new TF1("func_truth","gaus",-0.5,0.5);
0074      TF1 *func_real = new TF1("func_real","gaus",-0.5,0.5);
0075
0076     for(int iptbin=0; iptbin<nptbins; ++iptbin){
0077
0078    TCanvas *cp = new TCanvas("cp","cp",1400,1000);
0079    cp->SetMargin(0.10, 0.05 ,0.1,0.07);
0080
0081     double ptmin = (1.0-variation)*pt[iptbin]; // 10% range
0082     double ptmax = (1.0+variation)*pt[iptbin]; // 10% range
0083
0084     Int_t ptmin_bin = hist_d0xy_truth->GetZaxis()->FindBin(ptmin+0.0001);
0085     Int_t ptmax_bin = hist_d0xy_truth->GetZaxis()->FindBin(ptmax-0.0001);
0086
0087     TH1D *histd0xy_truth_1d = (TH1D*)hist_d0xy_truth->ProjectionX(Form("histd0xy_truth_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o");
0088     histd0xy_truth_1d->SetTitle(Form("d0_{xy} (truth): %1.1f <#eta< %1.1f && %1.2f <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));
0089     histd0xy_truth_1d->SetName(Form("eta_%1.1f_%1.1f_d0xy_truth_pt_%1.1f",etamin,etamax,pt[iptbin]));
0090
0091    if (histd0xy_truth_1d->GetEntries()<100) continue;
0092    double mu_truth = histd0xy_truth_1d->GetMean();
0093    double sigma_truth = histd0xy_truth_1d->GetStdDev();
0094    func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range
0095    histd0xy_truth_1d->Fit(func_truth,"NR+");
0096    mu_truth = func_truth->GetParameter(2);
0097    sigma_truth = func_truth->GetParError(2);
0098    func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range
0099    histd0xy_truth_1d->Fit(func_truth,"R+");
0100    float truth_par2 = func_truth->GetParameter(2)*10000; // cm to mum 10000 factor
0101    float truth_par2_err = func_truth->GetParError(2)*10000;
0102    momV_truth.push_back(pt[iptbin]);
0103    dcaTresolV_truth.push_back(truth_par2);
0104    err_dcaTresolV_truth.push_back(truth_par2_err);
0105
0106     TH1D *histd0xy_real_1d = (TH1D*)hist_d0xy_real->ProjectionX(Form("histd0xy_real_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o");
0107     histd0xy_real_1d->SetTitle(Form("d0_{xy} (real): %1.1f <#eta< %1.1f && %1.2f <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));
0108     histd0xy_real_1d->SetName(Form("eta_%1.1f_%1.1f_d0xy_real_pt_%1.1f",etamin,etamax,pt[iptbin]));
0109
0110    if (histd0xy_real_1d->GetEntries()<100) continue;
0111    double mu_real = histd0xy_real_1d->GetMean();
0112    double sigma_real = histd0xy_real_1d->GetStdDev();
0113    func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range
0114    histd0xy_real_1d->Fit(func_real,"NR+");
0115    mu_real = func_real->GetParameter(1);
0116    sigma_real = func_real->GetParameter(2);
0117    func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range
0118    histd0xy_real_1d->Fit(func_real,"R+");
0119    float real_par2 = func_real->GetParameter(2)*10000; // cm to mum 10000 factor
0120    float real_par2_err = func_real->GetParError(2)*10000;
0121    momV_real.push_back(pt[iptbin]);
0122    dcaTresolV_real.push_back(real_par2);
0123    err_dcaTresolV_real.push_back(real_par2_err);
0124
0125    cp->cd();
0126    histd0xy_truth_1d->Draw();
0127    cp->SaveAs(Form("Debug_Plots/truth/%s/dca/truth_dcaxy_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax));
0128    cp->Clear();
0129    cp->cd();
0130    histd0xy_real_1d->Draw();
0131    cp->SaveAs(Form("Debug_Plots/real/%s/dca/real_dcaxy_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax));
0132    }   // ptbin
0133
0134     const int size_truth = momV_truth.size();
0135     double pt_truth[size_truth], err_pt_truth[size_truth], sigma_dcaxy_truth[size_truth], err_sigma_dcaxy_truth[size_truth];
0136
0137     for (int i=0; i<size_truth; i++){
0138     pt_truth[i] = momV_truth.at(i);
0139     sigma_dcaxy_truth[i] = dcaTresolV_truth.at(i);
0140     err_sigma_dcaxy_truth[i] = err_dcaTresolV_truth.at(i);
0141     err_pt_truth[i] = 0.;
0142     }
0143
0144      const int size_real = momV_real.size();
0145     double pt_real[size_real], err_pt_real[size_real], sigma_dcaxy_real[size_real], err_sigma_dcaxy_real[size_real];
0146
0147     for (int i=0; i<size_real; i++){
0148     pt_real[i] = momV_real.at(i);
0149      sigma_dcaxy_real[i] = dcaTresolV_real.at(i);
0150     err_sigma_dcaxy_real[i] = err_dcaTresolV_real.at(i);
0151     err_pt_real[i] = 0.;
0152     }
0153
0154      TFile *fout = new TFile(Form("Final_Results/%s/dca/dcaxy_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
0155     TGraphErrors *gr1 = new TGraphErrors(size_truth,pt_truth,sigma_dcaxy_truth,err_pt_truth,err_sigma_dcaxy_truth);
0156      gr1->SetName("gr_truthseed");
0157     gr1->SetMarkerStyle(25);
0158     gr1->SetMarkerColor(kMagenta);
0159     gr1->SetMarkerSize(2.0);
0160     gr1->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)");
0161     gr1->GetXaxis()->CenterTitle();
0162     gr1->GetYaxis()->CenterTitle();
0163
0164      TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaxy_real,err_pt_real,err_sigma_dcaxy_real);
0165      gr2->SetName("gr_realseed");
0166     gr2->SetMarkerStyle(34);
0167     gr2->SetMarkerColor(kRed);
0168     gr2->SetMarkerSize(2.0);
0169     gr2->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{T}} (#mum)");
0170     gr2->GetXaxis()->CenterTitle();
0171     gr2->GetYaxis()->CenterTitle();
0172
0173 //  mgDCAT->Add(gr1);
0174     mgDCAT->Add(gr2);
0175     c_dcaxy->cd();
0176 //  c_dcaxy->SetLogy();
0177      mgDCAT->GetXaxis()->SetRangeUser(0.18,mgDCAT->GetXaxis()->GetXmax());
0178        mgDCAT->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY()));
0179      mgDCAT->Draw("AP");
0180     // lDCAT->AddEntry(gr1,"Truth Seeding");
0181      lDCAT->AddEntry(gr2,"Realistic Seeding");
0182      lDCAT->Draw("same");
0183      draw_req_DCA(etamin,etamax,0.,mgDCAT->GetXaxis()->GetXmax());
0184      c_dcaxy->SaveAs(Form("Final_Results/%s/dca/dcaxy_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax));
0185
0186    // Write the numbers in output file for comparisons
0187    outfile<<"ePIC"<<setw(20)<<epic.Data()<<setw(20)<<"EICRecon"<<setw(20)<<eicrecon.Data()<<endl;
0188    outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"Pt (GeV/c) \t"<<setw(20)<<"DCAT Resol #mum (Real)"<<endl;
0189    for (Int_t i = 0; i<gr1->GetN(); ++i){
0190    double x,ytrue, yreal;
0191    gr2->GetPoint(i,x,yreal);
0192    outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<yreal<<endl;
0193    }
0194     outfile.close();
0195
0196      fout->cd();
0197      mgDCAT->SetName(Form("dcaxy_resol_%1.1f_eta_%1.1f",etamin,etamax));
0198      mgDCAT->Write();
0199      fout->Close();
0200 }
0201
0202 //From Yellow report from section 11.2.2
0203 //===Fit Pointing Resolution
0204 float FitPointingAngle(Double_t *x, Double_t *par)
0205 {
0206   float func = sqrt((par[0]*par[0])/(x[0]*x[0])+par[1]*par[1]);
0207   return func;
0208 }
0209
0210 void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.)
0211 {
0212   TF1 *PWGReq_DCA2D;
0213   if (etamin >= -3.5 && etamax <= -2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax);
0214   else if (etamin >= -2.5 && etamax <= -1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax);
0215   else if (etamin >= -1.0 && etamax <= 1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((20./x)^2+5.0^2)",xmin,xmax);
0216   else if (etamin >= 1.0 && etamax <= 2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax);
0217   else if (etamin >= 2.5 && etamax <= 3.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax);
0218   else return;
0219   PWGReq_DCA2D->SetLineStyle(7);
0220   PWGReq_DCA2D->SetLineWidth(3.0);
0221   PWGReq_DCA2D->SetLineColor(kBlue);
0222   PWGReq_DCA2D->Draw("same");
0223
0224   TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
0225   l->SetTextSize(0.03);
0226   l->SetBorderSize(0);
0227   l->AddEntry(PWGReq_DCA2D,"PWGReq","l");
0228   l->Draw("same");
0229 }