benchmarks/tracking_performances/doCompare_truth_real_widebins_dcaz.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 #define mpi 0.139  // 1.864 GeV/c^2
0011
0012 void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.);
0013 void doCompare_truth_real_widebins_dcaz(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, Bool_t drawreq=1, TString extra_legend = "") // name = p, pt for getting p or pt dependence fitted results
0014 {
0015
0016 //=== style of the plot=========
0017    gStyle->SetPalette(1);
0018    gStyle->SetOptTitle(1);
0019    gStyle->SetTitleOffset(1.0,"XY");
0020    gStyle->SetTitleSize(.04,"XY");
0021    gStyle->SetLabelSize(.04,"XY");
0022    gStyle->SetHistLineWidth(2);
0023    gStyle->SetOptFit(1);
0024    gStyle->SetOptStat(1);
0025
0026    const Int_t nptbins = 10;
0027    double pt[nptbins] ={0.2, 0.3, 0.5,1.0, 1.5, 2.0, 5.0, 8.0, 10., 15.0};
0028    Double_t variation = 0.1; // 10 % variation
0029    std::vector<double> momV_truth, momV_real;
0030    std::vector<double> dcaZresolV_truth, err_dcaZresolV_truth, dcaZresolV_real, err_dcaZresolV_real;
0031    momV_truth.clear(); momV_real.clear(); dcaZresolV_truth.clear(); err_dcaZresolV_truth.clear(); dcaZresolV_real.clear(); err_dcaZresolV_real.clear();
0032
0033    TString symbolname = "";
0034    if (particle == "pi-") symbolname = "#pi^{-}";
0035    else symbolname = particle;
0036
0037    TF1 *f1=new TF1("f1","FitPointingAngle",0.,30.0,2);
0038    f1->SetParLimits(0,0.,50000);
0039    f1->SetParLimits(1,0.,50000);
0040
0041
0042        TCanvas *c_dcaz = new TCanvas("c_dcaz","c_dcaz",1400,1000);
0043        c_dcaz->SetMargin(0.10, 0.05 ,0.1,0.05);
0044        c_dcaz->SetGridy();
0045
0046      //Reading the root file
0047       TFile *fDCA_truth, *fDCA_real;
0048
0049       TMultiGraph *mgDCAZ;
0050       TLegend *lDCAZ;
0051       mgDCAZ = new TMultiGraph("mgDCAZ",";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)");
0052       lDCAZ = new TLegend(0.70,0.80,0.90,0.93);
0053       lDCAZ->SetTextSize(0.03);
0054       lDCAZ->SetBorderSize(0);
0055     lDCAZ->SetHeader(extra_legend.Data(), "C");
0056     lDCAZ->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.Data(), etamin, etamax), "C");
0057
0058     fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data()));
0059       fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data()));
0060
0061      // Truth seeding histograms
0062      TH3D *hist_d0z_truth = (TH3D*) fDCA_truth->Get("h_d0z_3d");
0063    TH3D *hist_d0z_real = (TH3D*) fDCA_real->Get("h_d0z_3d");
0064
0065       // d0z calculation for truth/real (binning are same in both cases)
0066   Int_t etamin_bin = hist_d0z_truth->GetYaxis()->FindBin(etamin+0.0001);
0067     Int_t etamax_bin = hist_d0z_truth->GetYaxis()->FindBin(etamax-0.0001);
0068
0069     TF1 *func_truth = new TF1("func_truth","gaus",-0.5,0.5);
0070      TF1 *func_real = new TF1("func_real","gaus",-0.5,0.5);
0071
0072     for(int iptbin=0; iptbin<nptbins; ++iptbin){
0073
0074    TCanvas *cp = new TCanvas("cp","cp",1400,1000);
0075    cp->SetMargin(0.10, 0.05 ,0.1,0.07);
0076
0077     double ptmin = (1.0-variation)*pt[iptbin]; // 10% range
0078     double ptmax = (1.0+variation)*pt[iptbin]; // 10% range
0079
0080     Int_t ptmin_bin = hist_d0z_truth->GetZaxis()->FindBin(ptmin+0.0001);
0081     Int_t ptmax_bin = hist_d0z_truth->GetZaxis()->FindBin(ptmax-0.0001);
0082
0083     TH1D *histd0z_truth_1d = (TH1D*)hist_d0z_truth->ProjectionX(Form("histd0z_truth_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o");
0084     histd0z_truth_1d->SetTitle(Form("d0_{z} (truth): %1.1f <#eta< %1.1f && %1.2f <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));
0085     histd0z_truth_1d->SetName(Form("eta_%1.1f_%1.1f_d0z_truth_pt_%1.1f",etamin,etamax,pt[iptbin]));
0086   // if (histd0z_truth_1d->GetEntries()<100) continue;
0087    double mu_truth = histd0z_truth_1d->GetMean();
0088    double sigma_truth = histd0z_truth_1d->GetStdDev();
0089    func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range
0090    histd0z_truth_1d->Fit(func_truth,"NR+");
0091    mu_truth = func_truth->GetParameter(2);
0092    sigma_truth = func_truth->GetParError(2);
0093    func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range
0094    histd0z_truth_1d->Fit(func_truth,"R+");
0095    float truth_par2 = func_truth->GetParameter(2)*10000; // cm to mum 10000 factor
0096    float truth_par2_err = func_truth->GetParError(2)*10000;
0097    momV_truth.push_back(pt[iptbin]);
0098    dcaZresolV_truth.push_back(truth_par2);
0099    err_dcaZresolV_truth.push_back(truth_par2_err);
0100
0101     TH1D *histd0z_real_1d = (TH1D*)hist_d0z_real->ProjectionX(Form("histd0z_real_eta%1.1f_%1.1f_pt%1.1f_%1.1f",etamin,etamax,ptmin,ptmax),etamin_bin,etamax_bin,ptmin_bin,ptmax_bin,"o");
0102     histd0z_real_1d->SetTitle(Form("d0_{z} (real): %1.1f <#eta< %1.1f && %1.2f <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));
0103     histd0z_real_1d->SetName(Form("eta_%1.1f_%1.1f_d0z_real_pt_%1.1f",etamin,etamax,pt[iptbin]));
0104
0105  //  if (histd0z_real_1d->GetEntries()<100) continue;
0106    double mu_real = histd0z_real_1d->GetMean();
0107    double sigma_real = histd0z_real_1d->GetStdDev();
0108    func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range
0109    histd0z_real_1d->Fit(func_real,"NR+");
0110    mu_real = func_real->GetParameter(1);
0111    sigma_real = func_real->GetParameter(2);
0112    func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range
0113    histd0z_real_1d->Fit(func_real,"R+");
0114    float real_par2 = func_real->GetParameter(2)*10000; // cm to mum 10000 factor
0115    float real_par2_err = func_real->GetParError(2)*10000;
0116    momV_real.push_back(pt[iptbin]);
0117    dcaZresolV_real.push_back(real_par2);
0118    err_dcaZresolV_real.push_back(real_par2_err);
0119
0120    cp->cd();
0121    histd0z_truth_1d->Draw();
0122    cp->SaveAs(Form("Debug_Plots/truth/%s/dca/truth_dcaz_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax));
0123    cp->Clear();
0124    cp->cd();
0125    histd0z_real_1d->Draw();
0126    cp->SaveAs(Form("Debug_Plots/real/%s/dca/real_dcaz_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),pt[iptbin],etamin,etamax));
0127    }   // ptbin
0128
0129     const int size_truth = momV_truth.size();
0130     double pt_truth[size_truth], err_pt_truth[size_truth], sigma_dcaz_truth[size_truth], err_sigma_dcaz_truth[size_truth];
0131
0132     for (int i=0; i<size_truth; i++){
0133     pt_truth[i] = momV_truth.at(i);
0134     sigma_dcaz_truth[i] = dcaZresolV_truth.at(i);
0135     err_sigma_dcaz_truth[i] = err_dcaZresolV_truth.at(i);
0136     err_pt_truth[i] = 0.;
0137     }
0138
0139      const int size_real = momV_real.size();
0140     double pt_real[size_real], err_pt_real[size_real], sigma_dcaz_real[size_real], err_sigma_dcaz_real[size_real];
0141
0142     for (int i=0; i<size_real; i++){
0143     pt_real[i] = momV_real.at(i);
0144      sigma_dcaz_real[i] = dcaZresolV_real.at(i);
0145     err_sigma_dcaz_real[i] = err_dcaZresolV_real.at(i);
0146     err_pt_real[i] = 0.;
0147     }
0148
0149      TFile *fout = new TFile(Form("Final_Results/%s/dca/dcaz_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
0150     TGraphErrors *gr1 = new TGraphErrors(size_truth,pt_truth,sigma_dcaz_truth,err_pt_truth,err_sigma_dcaz_truth);
0151      gr1->SetName("gr_truthseed");
0152     gr1->SetMarkerStyle(25);
0153     gr1->SetMarkerColor(kMagenta);
0154     gr1->SetMarkerSize(2.0);
0155     gr1->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)");
0156     gr1->GetXaxis()->CenterTitle();
0157     gr1->GetYaxis()->CenterTitle();
0158
0159      TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaz_real,err_pt_real,err_sigma_dcaz_real);
0160      gr2->SetName("gr_realseed");
0161     gr2->SetMarkerStyle(34);
0162     gr2->SetMarkerColor(kRed);
0163     gr2->SetMarkerSize(2.0);
0164     gr2->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)");
0165     gr2->GetXaxis()->CenterTitle();
0166     gr2->GetYaxis()->CenterTitle();
0167
0168     //mgDCAZ->Add(gr1);
0169     mgDCAZ->Add(gr2);
0170     c_dcaz->cd();
0171     //c_dcaz->SetLogy();
0172        mgDCAZ->GetXaxis()->SetRangeUser(0.18,mgDCAZ->GetXaxis()->GetXmax());
0173        mgDCAZ->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY()));
0174      mgDCAZ->Draw("AP");
0175     // lDCAZ->AddEntry(gr1,"Truth Seeding");
0176      lDCAZ->AddEntry(gr2,"Realistic Seeding");
0177      lDCAZ->Draw("same");
0178    //  draw_req_DCA(etamin,etamax,0.,mgDCAZ->GetXaxis()->GetXmax());
0179      c_dcaz->SaveAs(Form("Final_Results/%s/dca/dcaz_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax));
0180
0181      fout->cd();
0182      mgDCAZ->SetName(Form("dcaz_resol_%1.1f_eta_%1.1f",etamin,etamax));
0183      mgDCAZ->Write();
0184      fout->Close();
0185 }
0186
0187 //From Yellow report from section 11.2.2
0188 //===Fit Pointing Resolution
0189 float FitPointingAngle(Double_t *x, Double_t *par)
0190 {
0191   float func = sqrt((par[0]*par[0])/(x[0]*x[0])+par[1]*par[1]);
0192   return func;
0193 }
0194
0195 void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.)
0196 {
0197   TF1 *PWGReq_DCA2D;
0198   if (etamin >= -3.5 && etamax <= -2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax);
0199   else if (etamin >= -2.5 && etamax <= -1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax);
0200   else if (etamin >= -1.0 && etamax <= 1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((20./x)^2+5.0^2)",xmin,xmax);
0201   else if (etamin >= 1.0 && etamax <= 2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax);
0202   else if (etamin >= 2.5 && etamax <= 3.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax);
0203   else return;
0204   PWGReq_DCA2D->SetLineStyle(7);
0205   PWGReq_DCA2D->SetLineWidth(3.0);
0206   PWGReq_DCA2D->SetLineColor(kBlue);
0207   PWGReq_DCA2D->Draw("same");
0208
0209   TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
0210   l->SetTextSize(0.03);
0211   l->SetBorderSize(0);
0212   l->AddEntry(PWGReq_DCA2D,"PWGReq","l");
0213   l->Draw("same");
0214 }