benchmarks/lfhcal/doCompare_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_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,  momresolV, err_momresolV;
0029   momV.clear(); momresolV.clear(); err_momresolV.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[nfiles];
0046   TGraphErrors *gr_mom;
0047   TMultiGraph *mgMom;
0048   TLegend *lmom;
0049   mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %");
0050
0051   lmom = new TLegend(0.65,0.80,0.90,0.93);
0052   lmom->SetTextSize(0.03);
0053   lmom->SetBorderSize(0);
0054   lmom->SetHeader(extra_legend.Data(), "C");
0055   lmom->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.Data(), etamin, etamax), "C");
0056
0057   TF1 *func = new TF1("func","gaus",-0.5,0.5);
0058
0059   for (int i =0; i<nfiles; ++i){
0060
0061     TCanvas *cp = new TCanvas("cp","cp",1400,1000);
0062     cp->SetMargin(0.10, 0.05 ,0.1,0.07);
0063
0064     //pi-/mom/lfhcal_mom_20.0_mom_resol_pi-.root
0065     fmom[i] = TFile::Open(Form("./%s/mom/lfhcal_mom_%1.1f_mom_resol_%s.root",particle.Data(),mom[i],particle.Data()));
0066
0067     TH1D *hist = (TH1D*) fmom[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax));
0068     hist->Rebin(2);
0069     hist->SetName(Form("hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data()));
0070     hist->SetTitle(Form("Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax));
0071
0072     double mu = hist->GetMean();
0073     double sigma = hist->GetRMS();
0074     double sigma_err = hist->GetRMSError();
0075     /*
0076     hist->GetXaxis()->SetRangeUser(-1.0*range,1.0*range);
0077     func->SetRange(mu-2.0*sigma,mu+2.0*sigma); // fit with in 2 sigma range
0078     hist->Fit(func,"NR+");
0079     mu = func->GetParameter(1);
0080     sigma = func->GetParameter(2);
0081     func->SetRange(mu-2.0*sigma,mu+2.0*sigma);
0082     hist->Fit(func,"R+");
0083     float par2 = func->GetParameter(2)*100;
0084     float par2_err = func->GetParError(2)*100;
0085     */
0086     momV.push_back(mom[i]);
0087     momresolV.push_back(sigma);
0088     err_momresolV.push_back(sigma_err);
0089
0090     cp->cd();
0091     hist->Draw();
0092     //cp->SaveAs(Form("Debug_Plots/%s/mom/mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax));
0093   } // all files
0094
0095   const int size = momV.size();
0096   double p[size], err_p[size], sigma_p[size], err_sigma_p[size];
0097
0098   for (int i=0; i<size; i++){
0099     p[i] = momV.at(i);
0100     sigma_p[i] = momresolV.at(i);
0101     err_sigma_p[i] = err_momresolV.at(i);
0102     err_p[i] = 0.;
0103   }
0104
0105
0106   TFile *fout = new TFile(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.root",particle.Data(),particle.Data(),etamin,etamax),"recreate");
0107   TGraphErrors *gr1 = new TGraphErrors(size,p,sigma_p,err_p,err_sigma_p);
0108   gr1->SetName("grseed");
0109   gr1->SetMarkerStyle(25);
0110   gr1->SetMarkerColor(kBlue);
0111   gr1->SetMarkerSize(2.0);
0112   gr1->SetTitle(";p (GeV/c);#sigmap/p");
0113   gr1->GetXaxis()->CenterTitle();
0114   gr1->GetYaxis()->CenterTitle();
0115
0116
0117   mgMom->Add(gr1);
0118   c_mom->cd();
0119   mgMom->GetXaxis()->SetRangeUser(0.40,20.2);
0120   mgMom->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr1->GetN(),gr1->GetY())); // 50% more of the maximum value on yaxis
0121   mgMom->Draw("AP");
0122   lmom->AddEntry(gr1,"Nominal");
0123   lmom->Draw("same");
0124   //draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax());
0125   c_mom->SaveAs(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.png",particle.Data(),particle.Data(),etamin,etamax));
0126
0127   // Write the numbers in output file for comparisons
0128   outfile << extra_legend << endl;
0129   outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"p (GeV/c) \t"<<setw(20)<<"Resol  #mum "<<endl;
0130   for (Int_t i = 0; i<gr1->GetN(); ++i){
0131     double x,y;
0132     gr1->GetPoint(i,x,y);
0133     outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<y<<endl;
0134   }
0135   outfile.close();
0136
0137   fout->cd();
0138   mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax));
0139   mgMom->Write();
0140   fout->Close();
0141 }
0142
0143 //===Fit Momentum Resolution
0144 float FitMomentumResolution(Double_t *x, Double_t *par)
0145 {
0146   float func = sqrt(par[0]*par[0]*x[0]*x[0]+par[1]*par[1]);
0147   return func;
0148 }
0149
0150 //From Yellow report from section 11.2.2
0151
0152 // 1.2,1.5,2,2.5,3,3.5
0153
0154 /*
0155 void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.)
0156 {
0157
0158    TF1 *dd4hep_p;
0159    if (etamin >= -3.5 && etamax <= -2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax);
0160    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);
0161    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);
0162    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);
0163    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);
0164    else return;
0165    dd4hep_p->SetLineStyle(7);
0166    dd4hep_p->SetLineColor(kMagenta);
0167    dd4hep_p->SetLineWidth(3.0);
0168    dd4hep_p->Draw("same");
0169
0170   TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
0171   l->SetTextSize(0.03);
0172   l->SetBorderSize(0);
0173   l->AddEntry(dd4hep_p,"PWGReq","l");
0174   l->Draw("same");
0175  }
0176 */