]>
Commit | Line | Data |
---|---|---|
f887e5c3 | 1 | //This macro calculates the correction factor for the difference of the n sigma TOF distributions in data and mc |
2 | ||
3 | ||
4 | class AliSpectraBothHistoManager; | |
5 | AliSpectraBothHistoManager* managerdata=0x0; | |
6 | AliSpectraBothHistoManager* managermc=0x0; | |
7 | ||
8 | Double_t Gaus_plus_tail(Double_t* x, Double_t* par ) | |
9 | { | |
10 | ||
11 | Double_t mean=par[0]; | |
12 | Double_t rms=par[1]; | |
13 | Double_t c=par[2]; | |
14 | Double_t slope=par[3]/par[1];//df/dx continues | |
15 | Double_t cut=par[3]; | |
16 | ||
17 | ||
18 | ||
19 | Double_t one_over_sqrt_2pi=1.0/(TMath::Sqrt(2.0*TMath::Pi())); | |
20 | Double_t returnvalue=0.0; | |
21 | ||
22 | Double_t n=0.5*(1.0+TMath::Erf(cut/TMath::Sqrt(2.0)))+TMath::Exp(-cut*cut*0.5)*one_over_sqrt_2pi/(TMath::Abs(rms)*slope); | |
23 | if (x[0]<mean+cut*rms) | |
24 | returnvalue=TMath::Exp(-1.0*(x[0]-mean)*(x[0]-mean)/(2.0*rms*rms))*one_over_sqrt_2pi/TMath::Abs(rms); | |
25 | else | |
26 | returnvalue=TMath::Exp(slope*(mean+cut*rms-x[0]))*TMath::Exp(-cut*cut*0.5)*one_over_sqrt_2pi/TMath::Abs(rms); | |
27 | ||
28 | return c*returnvalue/n; | |
29 | } | |
30 | ||
31 | //Double_t background(Double_t x, Double_t scale, Double_t slope, Double_t startpoint) | |
32 | Double_t background(Double_t* x,Double_t* par) | |
33 | { | |
34 | //at start ponit it is eqaul to scale | |
35 | Double_t scale=par[0]; | |
36 | Double_t slope=par[1]; | |
37 | Double_t startpoint=par[2]; | |
38 | return scale*TMath::Exp((startpoint-x[0])*slope); | |
39 | } | |
40 | ||
41 | ||
42 | Double_t Gaus_on_background(Double_t* x,Double_t* par) | |
43 | { | |
44 | //par[0] peak position | |
45 | //par[1] sigma gaus | |
46 | //par[2] Normalization of signal | |
47 | //par[3] cut | |
48 | //par[4] scale of bg | |
49 | //par[5] slope of bg | |
50 | //par[6] start point of bg const | |
51 | //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<par[6]<<par[7]<<" "<<x[0]<< endl; | |
52 | Double_t returnvalue=0.0; | |
53 | returnvalue+=Gaus_plus_tail(x,par); | |
54 | //returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]); | |
55 | returnvalue+=background(x,&(par[4])); | |
56 | return returnvalue; | |
57 | } | |
58 | ||
59 | Int_t OpenFile(TString nameFile,TString outputname,Int_t mg1) | |
60 | { | |
61 | ||
62 | ||
63 | TFile *file = TFile::Open(nameFile.Data()); | |
64 | if(!file) | |
65 | { | |
66 | cout<<"no file"<<endl; | |
67 | return -1; | |
68 | } | |
69 | TDirectoryFile *dir=(TDirectoryFile*)file->Get(outputname.Data()); | |
70 | if(!dir) | |
71 | { | |
72 | cout<<"no dir "<<outputname.Data()<<endl; | |
73 | return -1; | |
74 | } | |
75 | if(!dir->Get("SpectraHistos")) | |
76 | return -1; | |
77 | if(mg1==1) | |
78 | { | |
79 | managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos"); | |
80 | return managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries(); | |
81 | ||
82 | } | |
83 | else | |
84 | { | |
85 | managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos"); | |
86 | return managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries(); | |
87 | } | |
88 | return -1; | |
89 | } | |
90 | ||
91 | ||
92 | void corrusingfitsigma(TString nameFile1,TString outputname1,TString nameFile2,TString outputname2, Float_t nsigmacut=3.0, Float_t minpt=0.6,Float_t maxpt=1.4, Bool_t saveoption=kFALSE) | |
93 | { | |
94 | ofstream oufiletxt("TOFcorrPID.txt"); | |
95 | Float_t nsigmacut=TMath::Abs(nsigmacut); | |
96 | TString pidmethods[3]={"TPC","TOF","TPCTOF"}; | |
97 | TString Particle[]={"Pion","Kaon","Proton"}; | |
98 | gStyle->SetOptStat(0); | |
99 | TH1::AddDirectory(kFALSE); | |
100 | gSystem->Load("libCore.so"); | |
101 | gSystem->Load("libPhysics.so"); | |
102 | gSystem->Load("libTree"); | |
103 | gSystem->Load("libMatrix"); | |
104 | gSystem->Load("libSTEERBase"); | |
105 | gSystem->Load("libESD"); | |
106 | gSystem->Load("libAOD"); | |
107 | gSystem->Load("libANALYSIS"); | |
108 | gSystem->Load("libOADB"); | |
109 | gSystem->Load("libANALYSISalice"); | |
110 | gSystem->Load("libTENDER"); | |
111 | gSystem->Load("libCORRFW"); | |
112 | gSystem->Load("libPWGTools"); | |
113 | gSystem->Load("libPWGLFspectra"); | |
114 | Int_t ne1= OpenFile(nameFile1,outputname1,1); | |
115 | Int_t ne2= OpenFile(nameFile2,outputname2,2); | |
116 | Int_t imethod=1; | |
117 | oufiletxt<<"file 1 "<<nameFile1.Data()<<endl; | |
118 | oufiletxt<<"outputname 1 "<<outputname1.Data()<<endl; | |
119 | oufiletxt<<"file 2 "<<nameFile2.Data()<<endl; | |
120 | oufiletxt<<"outputname 2 "<<outputname2.Data()<<endl; | |
121 | oufiletxt<<"nsigma "<<nsigmacut<<endl; | |
122 | oufiletxt<<"minpt "<<minpt<<endl; | |
123 | oufiletxt<<"maxpt "<<maxpt<<endl; | |
124 | ||
125 | TList* lout=new TList(); | |
126 | ||
127 | for(Int_t ipart=0;ipart<3;ipart++) | |
128 | { | |
129 | if(!managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data()))) | |
130 | continue; | |
131 | TH2F *nsig_data = (TH2F*)((TH2F*)managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone(); | |
132 | if(!nsig_data) | |
133 | continue; | |
134 | ||
135 | if(!managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data()))) | |
136 | continue; | |
137 | TH2F *nsig_mc = (TH2F*)((TH2F*)managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone(); | |
138 | if(!nsig_mc) | |
139 | continue; | |
140 | Int_t firstbin=nsig_data->GetXaxis()->FindBin(minpt); | |
141 | Int_t lastbin=nsig_data->GetXaxis()->FindBin(maxpt); | |
142 | Float_t ptlow=nsig_data->GetXaxis()->GetBinLowEdge(firstbin); | |
143 | Float_t ptup=nsig_data->GetXaxis()->GetBinUpEdge(lastbin); | |
144 | ||
145 | TH1D* hist1data= new TH1D(Form("datahist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup); | |
146 | TH1D* hist2data= new TH1D(Form("datahist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup); | |
147 | TH1D* sigmadata= new TH1D(Form("sigmadata%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup); | |
148 | ||
149 | hist1data->SetLineColor(kBlack); | |
150 | hist2data->SetLineColor(kBlack); | |
151 | sigmadata->SetMarkerColor(kBlack); | |
152 | sigmadata->SetMarkerStyle(22); | |
153 | hist1data->SetLineStyle(1); | |
154 | hist2data->SetLineStyle(2); | |
155 | hist1data->SetLineWidth(2); | |
156 | hist2data->SetLineWidth(2); | |
157 | ||
158 | ||
159 | TH1D* hist1mc= new TH1D(Form("mchist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup); | |
160 | TH1D* hist2mc= new TH1D(Form("mchist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup); | |
161 | TH1D* sigmamc= new TH1D(Form("sigmamc%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup); | |
162 | ||
163 | hist1mc->SetLineColor(kRed); | |
164 | hist2mc->SetLineColor(kRed); | |
165 | sigmamc->SetMarkerColor(kRed); | |
166 | sigmamc->SetMarkerStyle(23); | |
167 | sigmamc->SetLineColor(kRed); | |
168 | hist1mc->SetLineStyle(1); | |
169 | hist2mc->SetLineStyle(2); | |
170 | hist1mc->SetLineWidth(2); | |
171 | hist2mc->SetLineWidth(2); | |
172 | ||
173 | ||
174 | ||
175 | for(int ibin=firstbin;ibin<lastbin;ibin++) | |
176 | { | |
177 | ||
178 | TH1F *nsig_data_Proj1=(TH1F*)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin); | |
179 | nsig_data_Proj1->Sumw2(); | |
180 | nsig_data_Proj1->GetXaxis()->SetRangeUser(-6,6); | |
181 | nsig_data_Proj1->SetLineColor(kBlack); | |
182 | TH1F *nsig_mc_Proj1=(TH1F*)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin); | |
183 | nsig_mc_Proj1->Sumw2(); | |
184 | nsig_mc_Proj1->GetXaxis()->SetRangeUser(-6,6); | |
185 | nsig_mc_Proj1->SetLineColor(kRed); | |
186 | ||
187 | TF1* gfundata=new TF1("fun_for_fitdata",Gaus_on_background,-5,5,7); | |
188 | gfundata->SetParameter(0,0.0); | |
189 | gfundata->SetParameter(1,1.0); | |
190 | gfundata->SetParameter(2,nsig_data_Proj1->GetBinContent(nsig_data_Proj1->GetXaxis()->FindBin(0.0))); | |
191 | gfundata->FixParameter(3,1.0); | |
192 | gfundata->FixParameter(4,0.0); | |
193 | gfundata->FixParameter(5,1.0); | |
194 | gfundata->FixParameter(6,-5); | |
195 | gfundata->SetLineColor(kBlack); | |
196 | nsig_data_Proj1->Fit(gfundata,"SR+","E"); | |
197 | Float_t sigmadatav=gfundata->GetParameter(1); | |
198 | Float_t meandatav=gfundata->GetParameter(0); | |
199 | sigmadata->SetBinContent(ibin-firstbin+1,gfundata->GetParameter(1)); | |
200 | sigmadata->SetBinError(ibin-firstbin+1,gfundata->GetParError(1)); | |
201 | ||
202 | ||
203 | TF1* gfunmc=new TF1("fun_for_fitmc",Gaus_on_background,-5,5,7); | |
204 | gfunmc->SetParameter(0,0.0); | |
205 | gfunmc->SetParameter(1,1.0); | |
206 | gfunmc->SetParameter(2,nsig_mc_Proj1->GetBinContent(nsig_mc_Proj1->GetXaxis()->FindBin(0.0))); | |
207 | gfunmc->FixParameter(3,1.0); | |
208 | gfunmc->FixParameter(4,0.0); | |
209 | gfunmc->FixParameter(5,1.0); | |
210 | gfunmc->FixParameter(6,-5); | |
211 | gfunmc->SetLineColor(kRed); | |
212 | nsig_mc_Proj1->Fit(gfunmc,"SR+","E"); | |
213 | sigmamc->SetBinContent(ibin-firstbin+1,gfunmc->GetParameter(1)); | |
214 | sigmamc->SetBinError(ibin-firstbin+1,gfunmc->GetParError(1)); | |
215 | Float_t sigmamcv=gfunmc->GetParameter(1); | |
216 | Float_t meanmcv=gfundata->GetParameter(0); | |
217 | ||
218 | ||
176a37e3 | 219 | hist1data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001))); |
220 | Float_t width=nsig_data_Proj1->GetXaxis()->GetBinWidth(-1.0*nsigmacut); | |
221 | Int_t lowbind=nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav); | |
222 | Int_t upbind=nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav); | |
223 | Float_t lowbinpartd=TMath::Abs(nsig_data_Proj1->GetXaxis()->GetBinLowEdge(lowbind)+1.0*nsigmacut*sigmadatav-meandatav); | |
224 | Float_t upbinpartd=TMath::Abs(nsigmacut*sigmadatav+meandatav-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbind)); | |
225 | ||
226 | ||
227 | ||
228 | Float_t corrd=(lowbinpartd*nsig_data_Proj1->GetBinContent(lowbind)+upbinpartd*nsig_data_Proj1->GetBinContent(upbind))/width; | |
229 | ||
230 | hist2data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))-corrd); | |
231 | cout<<corrd/nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))<<endl; | |
232 | ||
f887e5c3 | 233 | hist1mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001))); |
176a37e3 | 234 | Int_t lowbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv); |
235 | Int_t upbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv); | |
236 | Float_t lowbinpartmc=TMath::Abs(nsig_mc_Proj1->GetXaxis()->GetBinLowEdge(lowbinmc)+1.0*nsigmacut*sigmamcv-meanmcv); | |
237 | Float_t upbinpartmc=TMath::Abs(nsigmacut*sigmamcv+meanmcv-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbinmc)); | |
238 | Float_t corrmc=(lowbinpartmc*nsig_mc_Proj1->GetBinContent(lowbinmc)+upbinpartmc*nsig_mc_Proj1->GetBinContent(upbinmc))/width; | |
239 | cout<<corrmc/nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))<<" aaa "<<endl; | |
240 | ||
241 | hist2mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))-corrmc); | |
f887e5c3 | 242 | TCanvas* ctmp=new TCanvas(Form("c%d%s",ibin,Particle[ipart].Data()),Form("c%d%s",ibin,Particle[ipart].Data()),800,800); |
243 | ctmp->cd()->SetLogy(); | |
244 | ||
245 | nsig_data_Proj1->GetXaxis()->SetTitle("#sigma_{TOF}"); | |
246 | nsig_data_Proj1->GetYaxis()->SetTitle("N_{entries}"); | |
247 | nsig_data_Proj1->SetTitle(Form("PionTOF[%.2f,%.2f]",nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin))); | |
248 | nsig_data_Proj1->Draw("E1"); | |
249 | nsig_mc_Proj1->Draw("E1same"); | |
250 | if(saveoption) | |
251 | ctmp->SaveAs(Form("c%d%s.png",ibin,Particle[ipart].Data())); | |
252 | delete ctmp; | |
253 | delete gfundata; | |
254 | delete gfunmc; | |
255 | delete nsig_data_Proj1; | |
256 | delete nsig_mc_Proj1; | |
257 | ||
258 | } | |
259 | ||
176a37e3 | 260 | hist1data->Sumw2(); |
261 | hist1data->Print("all"); | |
262 | hist1data->Sumw2(); | |
263 | hist1mc->Sumw2(); | |
264 | hist2mc->Sumw2(); | |
f887e5c3 | 265 | |
266 | hist2data->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
267 | hist2data->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)"); | |
268 | hist2mc->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
269 | hist2mc->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)"); | |
270 | hist2data->Divide(hist2data,hist1data,1,1,"B"); | |
271 | hist2mc->Divide(hist2mc,hist1mc,1,1,"B"); | |
272 | lout->Add(hist2data); | |
273 | lout->Add(hist2mc); | |
274 | ||
275 | sigmadata->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
276 | sigmadata->GetYaxis()->SetTitle("#sigma_{TOF} fitted"); | |
277 | sigmamc->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
278 | sigmamc->GetYaxis()->SetTitle("#sigma_{TOF} fitted"); | |
279 | ||
280 | ||
281 | ||
282 | lout->Add(sigmadata); | |
283 | lout->Add(sigmamc); | |
284 | ||
285 | TH1D* corr=(TH1D*)hist2data->Clone(Form("corr%s",Particle[ipart].Data())); | |
286 | corr->Divide(hist2mc); | |
287 | corr->SetLineColor(kBlue); | |
288 | corr->GetYaxis()->SetTitle("correction"); | |
176a37e3 | 289 | corr->Fit("pol0","R0","",minpt,maxpt); |
f887e5c3 | 290 | |
291 | oufiletxt<<Particle[ipart].Data()<<" "<<((TF1*)corr->GetListOfFunctions()->At(0))->GetParameter(0)<<endl; | |
292 | lout->Add(corr); | |
176a37e3 | 293 | TCanvas* corrcanvas=new TCanvas(Form("corrcanvas%s",Particle[ipart].Data()), Form("corrcanvas%s",Particle[ipart].Data()),600,600); |
294 | corrcanvas->SetMargin(0.12,0.05,0.15,0.1); | |
295 | corr->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
296 | corr->GetYaxis()->SetTitle("TOF PID mismatch corr."); | |
297 | corr->GetYaxis()->SetTitleOffset(1.5*corr->GetYaxis()->GetTitleOffset()); | |
298 | corr->GetYaxis()->SetRangeUser(0.98,1.04); | |
299 | ||
300 | corrcanvas->cd(); | |
301 | corr->Draw("E1"); | |
302 | corrcanvas->SaveAs(Form("Corr%s.eps",Particle[ipart].Data())); | |
303 | corrcanvas->SaveAs(Form("Corr%s.png",Particle[ipart].Data())); | |
304 | ||
f887e5c3 | 305 | } |
306 | TFile* fileout=TFile::Open("TOFcorrPID.root","Recreate"); | |
307 | ||
308 | lout->Write("output",TObject::kSingleKey); | |
309 | fileout->Close(); | |
310 | } |