File indexing completed on 2024-09-27 07:02:38
0001
0002
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
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 epic ="", TString eicrecon = "")
0014 {
0015
0016
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;
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 ofstream outfile;
0038 outfile.open ("DCAZ_resol.txt",ios_base::app);
0039
0040 TF1 *f1=new TF1("f1","FitPointingAngle",0.,30.0,2);
0041 f1->SetParLimits(0,0.,50000);
0042 f1->SetParLimits(1,0.,50000);
0043
0044
0045 TCanvas *c_dcaz = new TCanvas("c_dcaz","c_dcaz",1400,1000);
0046 c_dcaz->SetMargin(0.10, 0.05 ,0.1,0.05);
0047 c_dcaz->SetGridy();
0048
0049
0050 TFile *fDCA_truth, *fDCA_real;
0051
0052 TMultiGraph *mgDCAZ;
0053 TLegend *lDCAZ;
0054 mgDCAZ = new TMultiGraph("mgDCAZ",";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)");
0055 lDCAZ = new TLegend(0.70,0.80,0.90,0.93);
0056 lDCAZ->SetTextSize(0.03);
0057 lDCAZ->SetBorderSize(0);
0058 lDCAZ->SetHeader(Form("%s ePIC(%s/%s): %1.1f < #eta < %1.1f",symbolname.Data(),epic.Data(),eicrecon.Data(),etamin,etamax),"C");
0059
0060 fDCA_truth = TFile::Open(Form("./truthseed/%s/dca/final_hist_dca_truthseed.root",particle.Data()));
0061 fDCA_real = TFile::Open(Form("./realseed/%s/dca/final_hist_dca_realseed.root",particle.Data()));
0062
0063
0064 TH3D *hist_d0z_truth = (TH3D*) fDCA_truth->Get("h_d0z_3d");
0065 TH3D *hist_d0z_real = (TH3D*) fDCA_real->Get("h_d0z_3d");
0066
0067
0068 Int_t etamin_bin = hist_d0z_truth->GetYaxis()->FindBin(etamin+0.0001);
0069 Int_t etamax_bin = hist_d0z_truth->GetYaxis()->FindBin(etamax-0.0001);
0070
0071 TF1 *func_truth = new TF1("func_truth","gaus",-0.5,0.5);
0072 TF1 *func_real = new TF1("func_real","gaus",-0.5,0.5);
0073
0074 for(int iptbin=0; iptbin<nptbins; ++iptbin){
0075
0076 TCanvas *cp = new TCanvas("cp","cp",1400,1000);
0077 cp->SetMargin(0.10, 0.05 ,0.1,0.07);
0078
0079 double ptmin = (1.0-variation)*pt[iptbin];
0080 double ptmax = (1.0+variation)*pt[iptbin];
0081
0082 Int_t ptmin_bin = hist_d0z_truth->GetZaxis()->FindBin(ptmin+0.0001);
0083 Int_t ptmax_bin = hist_d0z_truth->GetZaxis()->FindBin(ptmax-0.0001);
0084
0085 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");
0086 histd0z_truth_1d->SetTitle(Form("d0_{z} (truth): %1.1f <#eta< %1.1f && %1.2f <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));
0087 histd0z_truth_1d->SetName(Form("eta_%1.1f_%1.1f_d0z_truth_pt_%1.1f",etamin,etamax,pt[iptbin]));
0088 if (histd0z_truth_1d->GetEntries()<100) continue;
0089 double mu_truth = histd0z_truth_1d->GetMean();
0090 double sigma_truth = histd0z_truth_1d->GetStdDev();
0091 func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth);
0092 histd0z_truth_1d->Fit(func_truth,"NR+");
0093 mu_truth = func_truth->GetParameter(2);
0094 sigma_truth = func_truth->GetParError(2);
0095 func_truth->SetRange(mu_truth-2.0*sigma_truth,mu_truth+2.0*sigma_truth);
0096 histd0z_truth_1d->Fit(func_truth,"R+");
0097 float truth_par2 = func_truth->GetParameter(2)*10000;
0098 float truth_par2_err = func_truth->GetParError(2)*10000;
0099 momV_truth.push_back(pt[iptbin]);
0100 dcaZresolV_truth.push_back(truth_par2);
0101 err_dcaZresolV_truth.push_back(truth_par2_err);
0102
0103 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");
0104 histd0z_real_1d->SetTitle(Form("d0_{z} (real): %1.1f <#eta< %1.1f && %1.2f <p_{T}< %1.2f",etamin,etamax,ptmin,ptmax));
0105 histd0z_real_1d->SetName(Form("eta_%1.1f_%1.1f_d0z_real_pt_%1.1f",etamin,etamax,pt[iptbin]));
0106
0107 if (histd0z_real_1d->GetEntries()<100) continue;
0108 double mu_real = histd0z_real_1d->GetMean();
0109 double sigma_real = histd0z_real_1d->GetStdDev();
0110 func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real);
0111 histd0z_real_1d->Fit(func_real,"NR+");
0112 mu_real = func_real->GetParameter(1);
0113 sigma_real = func_real->GetParameter(2);
0114 func_real->SetRange(mu_real-2.0*sigma_real,mu_real+2.0*sigma_real);
0115 histd0z_real_1d->Fit(func_real,"R+");
0116 float real_par2 = func_real->GetParameter(2)*10000;
0117 float real_par2_err = func_real->GetParError(2)*10000;
0118 momV_real.push_back(pt[iptbin]);
0119 dcaZresolV_real.push_back(real_par2);
0120 err_dcaZresolV_real.push_back(real_par2_err);
0121
0122 cp->cd();
0123 histd0z_truth_1d->Draw();
0124 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));
0125 cp->Clear();
0126 cp->cd();
0127 histd0z_real_1d->Draw();
0128 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));
0129 }
0130
0131 const int size_truth = momV_truth.size();
0132 double pt_truth[size_truth], err_pt_truth[size_truth], sigma_dcaz_truth[size_truth], err_sigma_dcaz_truth[size_truth];
0133
0134 for (int i=0; i<size_truth; i++){
0135 pt_truth[i] = momV_truth.at(i);
0136 sigma_dcaz_truth[i] = dcaZresolV_truth.at(i);
0137 err_sigma_dcaz_truth[i] = err_dcaZresolV_truth.at(i);
0138 err_pt_truth[i] = 0.;
0139 }
0140
0141 const int size_real = momV_real.size();
0142 double pt_real[size_real], err_pt_real[size_real], sigma_dcaz_real[size_real], err_sigma_dcaz_real[size_real];
0143
0144 for (int i=0; i<size_real; i++){
0145 pt_real[i] = momV_real.at(i);
0146 sigma_dcaz_real[i] = dcaZresolV_real.at(i);
0147 err_sigma_dcaz_real[i] = err_dcaZresolV_real.at(i);
0148 err_pt_real[i] = 0.;
0149 }
0150
0151 TFile *fout = new TFile(Form("Final_Results/%s/dca/dcaz_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate");
0152 TGraphErrors *gr1 = new TGraphErrors(size_truth,pt_truth,sigma_dcaz_truth,err_pt_truth,err_sigma_dcaz_truth);
0153 gr1->SetName("gr_truthseed");
0154 gr1->SetMarkerStyle(25);
0155 gr1->SetMarkerColor(kMagenta);
0156 gr1->SetMarkerSize(2.0);
0157 gr1->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)");
0158 gr1->GetXaxis()->CenterTitle();
0159 gr1->GetYaxis()->CenterTitle();
0160
0161 TGraphErrors *gr2 = new TGraphErrors(size_real,pt_real,sigma_dcaz_real,err_pt_real,err_sigma_dcaz_real);
0162 gr2->SetName("gr_realseed");
0163 gr2->SetMarkerStyle(34);
0164 gr2->SetMarkerColor(kRed);
0165 gr2->SetMarkerSize(2.0);
0166 gr2->SetTitle(";p_{T} (GeV/c); #sigma_{DCA_{Z}} (#mum)");
0167 gr2->GetXaxis()->CenterTitle();
0168 gr2->GetYaxis()->CenterTitle();
0169
0170
0171 mgDCAZ->Add(gr2);
0172 c_dcaz->cd();
0173
0174 mgDCAZ->GetXaxis()->SetRangeUser(0.18,mgDCAZ->GetXaxis()->GetXmax());
0175 mgDCAZ->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr2->GetN(),gr2->GetY()));
0176 mgDCAZ->Draw("AP");
0177
0178 lDCAZ->AddEntry(gr2,"Realistic Seeding");
0179 lDCAZ->Draw("same");
0180
0181 c_dcaz->SaveAs(Form("Final_Results/%s/dca/dcaz_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax));
0182
0183
0184 outfile<<"ePIC"<<setw(20)<<epic.Data()<<setw(20)<<"EICRecon"<<setw(20)<<eicrecon.Data()<<endl;
0185 outfile<<"Etamin"<<setw(20)<<"Etamax"<<setw(20)<<"Pt (GeV/c) \t"<<setw(20)<<"DCAz Resol #mum (Real)"<<endl;
0186 for (Int_t i = 0; i<gr1->GetN(); ++i){
0187 double x,ytrue, yreal;
0188 gr2->GetPoint(i,x,yreal);
0189 outfile<<etamin<<setw(20)<<etamax<<setw(20)<<x<<setw(20)<<yreal<<endl;
0190 }
0191 outfile.close();
0192 fout->cd();
0193 mgDCAZ->SetName(Form("dcaz_resol_%1.1f_eta_%1.1f",etamin,etamax));
0194 mgDCAZ->Write();
0195 fout->Close();
0196 }
0197
0198
0199
0200 float FitPointingAngle(Double_t *x, Double_t *par)
0201 {
0202 float func = sqrt((par[0]*par[0])/(x[0]*x[0])+par[1]*par[1]);
0203 return func;
0204 }
0205
0206 void draw_req_DCA(double etamin, double etamax, double xmin=0., double xmax=0.)
0207 {
0208 TF1 *PWGReq_DCA2D;
0209 if (etamin >= -3.5 && etamax <= -2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax);
0210 else if (etamin >= -2.5 && etamax <= -1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax);
0211 else if (etamin >= -1.0 && etamax <= 1.0) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((20./x)^2+5.0^2)",xmin,xmax);
0212 else if (etamin >= 1.0 && etamax <= 2.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+20.0^2)",xmin,xmax);
0213 else if (etamin >= 2.5 && etamax <= 3.5) PWGReq_DCA2D = new TF1("PWGReq_DCA2D", "TMath::Sqrt((30./x)^2+40.0^2)",xmin,xmax);
0214 else return;
0215 PWGReq_DCA2D->SetLineStyle(7);
0216 PWGReq_DCA2D->SetLineWidth(3.0);
0217 PWGReq_DCA2D->SetLineColor(kBlue);
0218 PWGReq_DCA2D->Draw("same");
0219
0220 TLegend *l= new TLegend(0.70,0.75,0.90,0.80);
0221 l->SetTextSize(0.03);
0222 l->SetBorderSize(0);
0223 l->AddEntry(PWGReq_DCA2D,"PWGReq","l");
0224 l->Draw("same");
0225 }