TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / corrusingfitsigma.C
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
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
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)));
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);
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         
260                         hist1data->Sumw2();
261                         hist1data->Print("all");
262                         hist1data->Sumw2();
263                         hist1mc->Sumw2();
264                         hist2mc->Sumw2();
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");
289                 corr->Fit("pol0","R0","",minpt,maxpt);
290                 
291                 oufiletxt<<Particle[ipart].Data()<<" "<<((TF1*)corr->GetListOfFunctions()->At(0))->GetParameter(0)<<endl;
292                 lout->Add(corr);
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
305         }
306         TFile* fileout=TFile::Open("TOFcorrPID.root","Recreate");
307         
308         lout->Write("output",TObject::kSingleKey);
309         fileout->Close();
310 }