TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / corrusingfitsigma.C
CommitLineData
f887e5c3 1//This macro calculates the correction factor for the difference of the n sigma TOF distributions in data and mc
2
3
4class AliSpectraBothHistoManager;
5AliSpectraBothHistoManager* managerdata=0x0;
6AliSpectraBothHistoManager* managermc=0x0;
7
8Double_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)
32Double_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
42Double_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
59Int_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
92void 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");
af472fff 110 gSystem->Load("libTender");
f887e5c3 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}