benchmarks/tracking_performances/doCompare_truth_real_widebins_mom.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_Mom(double etamin, double etamax, double xmin=0., double xmax=0.);
0013 void doCompare_truth_real_widebins_mom(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, double range =0.3, 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 nfiles = 6;
0027    double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,20.0};
0028    std::vector<double> momV_truth, momV_real, momresolV_truth, err_momresolV_truth, momresolV_real, err_momresolV_real;
0029    momV_truth.clear(); momV_real.clear(); momresolV_truth.clear(); err_momresolV_truth.clear(); momresolV_real.clear(); err_momresolV_real.clear();
0030    TString symbolname = "";
0031    if (particle == "pi-") symbolname = "#pi^{-}";
0032    else symbolname = particle;
0033    ofstream outfile;
0034    outfile.open ("Mom_resol.txt",ios_base::app);
0035
0036    TF1 *f1=new TF1("f1","FitMomentumResolution",0.,30.0,2);
0037    f1->SetParLimits(0,0.,0.1);
0038    f1->SetParLimits(1,0.,5.0);
0039
0040    TCanvas *c_mom = new TCanvas("cmom","cmom",1400,1000);
0041    c_mom->SetMargin(0.10, 0.05 ,0.1,0.05);
0042    c_mom->SetGridy();
0043
0044  //Reading the root file
0045   TFile *fmom_truth[nfiles], *fmom_real[nfiles];
0046   TGraphErrors *gr_mom_truth, *gr_mom_real;
0047
0048   TMultiGraph *mgMom;
0049   TLegend *lmom;
0050   mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %");
0051
0052   lmom = new TLegend(0.65,0.80,0.90,0.93);
0053   lmom->SetTextSize(0.03);
0054   lmom->SetBorderSize(0);
0055   lmom->SetHeader(extra_legend.Data(), "C");
0056   lmom->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.Data(), etamin, etamax), "C");
0057
0058   TF1 *func_truth = new TF1("func_truth","gaus",-0.5,0.5);
0059   TF1 *func_real = new TF1("func_real","gaus",-0.5,0.5);
0060
0061   for (int i =0; i<nfiles; ++i){
0062
0063    TCanvas *cp = new TCanvas("cp","cp",1400,1000);
0064    cp->SetMargin(0.10, 0.05 ,0.1,0.07);
0065
0066      fmom_truth[i] = TFile::Open(Form("./truthseed/pi-/mom/Performances_mom_%1.1f_mom_resol_truth_%s.root",mom[i],particle.Data()));
0067      fmom_real[i] = TFile::Open(Form("./realseed/pi-/mom/Performances_mom_%1.1f_mom_resol_realseed_%s.root",mom[i],particle.Data()));
0068
0069      TH1D *hist_truth = (TH1D*) fmom_truth[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax));
0070      hist_truth->Rebin(2);
0071      hist_truth->SetName(Form("truthseed_hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data()));
0072      hist_truth->SetTitle(Form("Truth Seed (%s): Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax));
0073
0074      double mu_truth = hist_truth->GetMean();
0075      double sigma_truth = hist_truth->GetStdDev();
0076          hist_truth->GetXaxis()->SetRangeUser(-1.0*range,1.0*range);
0077      func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth); // fit with in 2 sigma range
0078      hist_truth->Fit(func_truth,"NR+");
0079      mu_truth = func_truth->GetParameter(1);
0080      sigma_truth = func_truth->GetParameter(2);
0081      func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth);
0082      hist_truth->Fit(func_truth,"R+");
0083      float truth_par2 = func_truth->GetParameter(2)*100;
0084      float truth_par2_err = func_truth->GetParError(2)*100;
0085      momV_truth.push_back(mom[i]);
0086      momresolV_truth.push_back(truth_par2);
0087      err_momresolV_truth.push_back(truth_par2_err);
0088
0089      TH1D *hist_real = (TH1D*) fmom_real[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax));
0090      hist_real->Rebin(2);
0091      hist_real->SetName(Form("realseed_hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data()));
0092      hist_real->SetTitle(Form("Realistic Seed (%s): Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax));
0093
0094      double mu_real = hist_real->GetMean();
0095      double sigma_real = hist_real->GetStdDev();
0096          hist_real->GetXaxis()->SetRangeUser(-1.0*range,1.0*range);
0097          func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real); // fit with in 2 sigma range
0098      hist_real->Fit(func_real,"NR+");
0099      mu_real = func_real->GetParameter(1);
0100      sigma_real = func_real->GetParameter(2);
0101      func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real);
0102      hist_real->Fit(func_real,"R+");
0103      float real_par2 = func_real->GetParameter(2)*100;
0104      float real_par2_err = func_real->GetParError(2)*100;
0105      momV_real.push_back(mom[i]);
0106      momresolV_real.push_back(real_par2);
0107      err_momresolV_real.push_back(real_par2_err);
0108      cp->cd();
0109      hist_truth->Draw();
0110      cp->SaveAs(Form("Debug_Plots/truth/%s/mom/truth_mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax));
0111      cp->Clear();
0112      cp->cd();
0113      hist_real->Draw();
0114      cp->SaveAs(Form("Debug_Plots/real/%s/mom/real_mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax));
0115  } // all files
0116
0117     const int size_truth = momV_truth.size();
0118     double p_truth[size_truth], err_p_truth[size_truth], sigma_p_truth[size_truth], err_sigma_p_truth[size_truth];
0119
0120     for (int i=0; i<size_truth; i++){
0121     p_truth[i] = momV_truth.at(i);
0122     sigma_p_truth[i] = momresolV_truth.at(i);
0123     err_sigma_p_truth[i] = err_momresolV_truth.at(i);
0124     err_p_truth[i] = 0.;
0125     }
0126
0127         const int size_real = momV_real.size();
0128     double p_real[size_real], err_p_real[size_real], sigma_p_real[size_real], err_sigma_p_real[size_real];
0129
0130     for (int i=0; i<size_real; i++){
0131     p_real[i] = momV_real.at(i);
0132         sigma_p_real[i] = momresolV_real.at(i);
0133     err_sigma_p_real[i] = err_momresolV_real.at(i);
0134     err_p_real[i] = 0.;
0135     }
0136
0137   TFile *fout = new TFile(Form("Final_Results/%s/mom/mom_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
0138     TGraphErrors *gr1 = new TGraphErrors(size_truth,p_truth,sigma_p_truth,err_p_truth,err_sigma_p_truth);
0139         gr1->SetName("gr_truthseed");
0140     gr1->SetMarkerStyle(25);
0141     gr1->SetMarkerColor(kBlue);
0142     gr1->SetMarkerSize(2.0);
0143     gr1->SetTitle(";p (GeV/c);#sigmap/p");
0144     gr1->GetXaxis()->CenterTitle();
0145     gr1->GetYaxis()->CenterTitle();
0146
0147         TGraphErrors *gr2 = new TGraphErrors(size_real,p_real,sigma_p_real,err_p_real,err_sigma_p_real);
0148         gr2->SetName("gr_realseed");
0149     gr2->SetMarkerStyle(34);
0150     gr2->SetMarkerColor(kRed);
0151     gr2->SetMarkerSize(2.0);
0152     gr2->SetTitle(";p (GeV/c);#sigmap/p");
0153     gr2->GetXaxis()->CenterTitle();
0154     gr2->GetYaxis()->CenterTitle();
0155
0156     mgMom->Add(gr1);
0157     mgMom->Add(gr2);
0158     c_mom->cd();
0159     mgMom->GetXaxis()->SetRangeUser(0.40,20.2);
0160     mgMom->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY())); // 50% more of the maximum value on yaxis
0161     mgMom->Draw("AP");
0162     lmom->AddEntry(gr1,"Truth Seeding");
0163     lmom->AddEntry(gr2,"Realistic Seeding");
0164     lmom->Draw("same");
0165     draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax());
0166     c_mom->SaveAs(Form("Final_Results/%s/mom/mom_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax));
0167
0168     // Write the numbers in output file for comparisons
0169     outfile << extra_legend << endl;
0170   outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"p (GeV/c) \t"<<setw(20)<<"Resol  #mum (Truth)"<<setw(20)<<"Resol #mum (Real)"<<endl;
0171   for (Int_t i = 0; i<gr1->GetN(); ++i){
0172   double x,ytrue, yreal;
0173   gr1->GetPoint(i,x,ytrue);    gr2->GetPoint(i,x,yreal);
0174   outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<ytrue<<setw(20)<<yreal<<endl;
0175   }
0176   outfile.close();
0177
0178     fout->cd();
0179     mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax));
0180     mgMom->Write();
0181     fout->Close();
0182 }
0183
0184 //===Fit Momentum Resolution
0185 float FitMomentumResolution(Double_t *x, Double_t *par)
0186 {
0187   float func = sqrt(par[0]*par[0]*x[0]*x[0]+par[1]*par[1]);
0188   return func;
0189 }
0190
0191 //From Yellow report from section 11.2.2
0192
0193 void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.)
0194 {
0195
0196    TF1 *dd4hep_p;
0197    if (etamin >= -3.5 && etamax <= -2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
0198    else if (etamin >= -2.5 && etamax <= -1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax);
0199    else if (etamin >= -1.0 && etamax <= 1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+0.5^2)",xmin,xmax);
0200    else if (etamin >= 1.0 && etamax <= 2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax);
0201    else if (etamin >= 2.5 && etamax <= 3.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
0202    else return;
0203    dd4hep_p->SetLineStyle(7);
0204    dd4hep_p->SetLineColor(kMagenta);
0205    dd4hep_p->SetLineWidth(3.0);
0206    dd4hep_p->Draw("same");
0207
0208   TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
0209   l->SetTextSize(0.03);
0210   l->SetBorderSize(0);
0211   l->AddEntry(dd4hep_p,"PWGReq","l");
0212   l->Draw("same");
0213  }