4 TString Charge[]={"Pos","Neg"};
5 TString Sign[]={"Plus","Minus"};
6 TString Particle[]={"Pion","Kaon","Proton"};
7 Int_t Color[3]={1,2,4};
8 Int_t Marker[6]={20,21,22,24,25,26};
9 Double_t projl[3]={0,20,80};
10 Double_t proju[3]={10,40,90};
11 Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
21 gSystem->Load("libCore.so");
22 //gSystem->Load("libGeom.so");
23 gSystem->Load("libPhysics.so");
24 //gSystem->Load("libVMC");
25 gSystem->Load("libTree");
26 //gSystem->Load("libProof");
27 gSystem->Load("libMatrix");
28 gSystem->Load("libSTEERBase");
29 gSystem->Load("libESD");
30 gSystem->Load("libAOD");
31 gSystem->Load("libANALYSIS");
32 gSystem->Load("libOADB");
33 gSystem->Load("libANALYSISalice");
34 gSystem->Load("libTENDER");
35 gSystem->Load("libCORRFW");
36 //gSystem->Load("libPWG0base");
37 //gSystem->Load("libMinuit");
38 gSystem->Load("libPWGTools");
39 gSystem->Load("libPWGLFSPECTRA");
40 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
44 mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
45 mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
46 mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
48 //TString fold="3SigmaPID_AOD086-090_FilterBit10";
49 TString fold="5SigmaPID_AOD046_FilterBit6";
50 Int_t ibinToCompare=-1;
52 TString sname="Cent0to5_QVec0.0to100.0";ibinToCompare=0;
54 TString dataFile = Form("output/%s/Pt.AOD.1._data_ptcut_%s.root",fold.Data(),sname.Data());
55 TString mcFile =Form("output/%s/Pt.AOD.1._MC_%s.root",fold.Data(),sname.Data());
56 gStyle->SetPalette(1);
58 // Open root MC file and get classes
59 cout << "Analysis Macro" << endl;
60 cout << " > Reading MC data" << endl;
61 TFile *_mc = TFile::Open(mcFile.Data());
62 AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
63 AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
64 AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
65 // print info about mc track and Event cuts
66 cout << " -- Info about MC -- "<< endl;
67 ecuts_mc->PrintCuts();
68 // tcuts_mc->PrintCuts();
69 // proceed likewise for data
70 TFile *_data = TFile::Open(dataFile.Data());
71 AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
72 AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
73 AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
74 // print info about track and Event cuts
75 cout << " -- Info about data -- " << endl;
76 ecuts_data->PrintCuts();
77 tcuts_data->PrintCuts();
79 //QAPlots(hman_data,hman_mc);
82 Printf("\n\n-> Calculating MC Correction Factors");
84 TCanvas *ceff=new TCanvas("ceff","ceff",700,500);
85 Bool_t UseMCDCACorrection=kTRUE;
86 for(Int_t icharge=0;icharge<2;icharge++){
87 for(Int_t ipart=0;ipart<3;ipart++){
88 Int_t index=ipart+3*icharge;
89 //BE CAREFUL! depending on the efficiency you choose, you must change the DCA correction (data or data/mc)
91 TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); //MC correction for Prim+Cont+Eff
92 //TString hname=Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());//MC correction for Prim+Eff
93 //TString hname=Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; //MC correction for Cont+Eff
94 //TString hname=Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; // Pure MC efficiency for Prim. BE CAREFUL WITH MUONS!!!
95 Printf("Getting %s",hname.Data());
96 CorrFact[index]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
97 CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
98 CorrFact[index]->SetTitle(Form("CorrFact %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
99 CorrFact[index]->SetMarkerStyle(Marker[index]);
100 CorrFact[index]->SetMarkerColor(Color[ipart]);
101 CorrFact[index]->SetLineColor(Color[ipart]);
102 hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
103 Printf("... and divide it by %s",hname.Data());
104 CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
105 if(index==0)CorrFact[index]->DrawClone();
106 else CorrFact[index]->DrawClone("same");
110 TFile *fESD=new TFile("EffAlex/pionEffPbPb.root");
111 TH1F *hEffESD=(TH1F*)fESD->Get("effMapPionTpcOnlyNeg0");
112 hEffESD->DrawClone("same");
116 printf("\n\n-> Spectra Normalization");
117 AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
118 Double_t events_data = ecuts_data->NumberOfEvents();
119 Printf(": accepted events: %.0f",events_data);
120 Double_t events_mc = ecuts_mc->NumberOfEvents();
121 Printf(": accepted events (MC): %.0f",events_mc);
123 //divide RAW for Correction Factor
124 Printf("\n\n-> Using MC correction factor to correct RAW spectra");
126 for(Int_t icharge=0;icharge<2;icharge++){
127 for(Int_t ipart=0;ipart<3;ipart++){
128 Int_t index=ipart+3*icharge;
129 TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
130 printf("Getting %s",hname.Data());
131 Spectra[index] =(TH1F*)((TH1F*) hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
132 Spectra[index]->SetName(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
133 Spectra[index]->SetTitle(Form("Spectra %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
134 Spectra[index]->SetMarkerStyle(Marker[index]);
135 Spectra[index]->SetMarkerColor(Color[ipart]);
136 Spectra[index]->SetLineColor(Color[ipart]);
137 Printf("... and divide it by %s",hname.Data());
138 Spectra[index]->Divide(CorrFact[index]);
139 Spectra[index]->Scale(1./events_data,"width");//NORMALIZATION
143 //Put Bin Content = 0
144 for(Int_t icharge=0;icharge<2;icharge++){
145 for(Int_t ipart=0;ipart<3;ipart++){
146 Int_t index=ipart+3*icharge;
147 for(Int_t ibin=0;ibin<Spectra[index]->GetNbinsX();ibin++){
148 if(Spectra[index]->GetBinCenter(ibin)<Range[ipart]){
149 Spectra[index]->SetBinContent(ibin,0);
150 Spectra[index]->SetBinError(ibin,0);
155 //DCA Correction with the "right" DCA sample
156 DCACorrection(Spectra,hman_data,hman_mc,UseMCDCACorrection);
158 //DCA Correction forcing loose DCA
159 // TString fold_LooseDCA="5SigmaPID_AOD046_FilterBit5";
160 // TString dataFile_LooseDCA = Form("output/%s/Pt.AOD.1._data_ptcut_%s.root",fold_LooseDCA.Data(),sname.Data());
161 // TString mcFile_LooseDCA =Form("output/%s/Pt.AOD.1._MC_%s.root",fold_LooseDCA.Data(),sname.Data());
162 // TFile *_mc_LooseDCA = TFile::Open(mcFile_LooseDCA.Data());
163 // AliSpectraAODHistoManager* hman_mc_LooseDCA = (AliSpectraAODHistoManager*) _mc_LooseDCA->Get("SpectraHistos");
164 // TFile *_data_LooseDCA = TFile::Open(dataFile_LooseDCA.Data());
165 // AliSpectraAODHistoManager* hman_data_LooseDCA = (AliSpectraAODHistoManager*) _data_LooseDCA->Get("SpectraHistos");
166 // DCACorrection(Spectra,hman_data_LooseDCA,hman_mc_LooseDCA,UseMCDCACorrection);
173 GFCorrection(Spectra,tcuts_data);
176 for(Int_t icharge=0;icharge<2;icharge++){
177 for(Int_t ipart=0;ipart<3;ipart++){
178 Int_t index=ipart+3*icharge;
179 hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
180 MCTruth[index]=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
181 MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
182 MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
183 MCTruth[index]->Scale(1./events_mc,"width");//NORMALIZATION
187 //Drawing Final Spectra
188 TCanvas *cspectra=new TCanvas("cspectra","cspectra",700,500);
192 for(Int_t icharge=0;icharge<2;icharge++){
193 for(Int_t ipart=0;ipart<3;ipart++){
194 Int_t index=ipart+3*icharge;
195 if(index==0)Spectra[index]->DrawClone();
196 else Spectra[index]->DrawClone("same");
202 Printf("\n\n-> Saving spectra in Out%s_%s.root",sname.Data(),fold.Data());
203 TFile * fout=new TFile(Form("results/Out_%s_%s.root",sname.Data(),fold.Data()),"RECREATE");
204 for(Int_t icharge=0;icharge<2;icharge++){
205 for(Int_t ipart=0;ipart<3;ipart++){
206 Int_t index=ipart+3*icharge;
207 Spectra[index]->Write();
208 CorrFact[index]->Write();
209 MCTruth[index]->Write();
213 //if Bin 0-5% with no cut ratio with combined analysis
214 if(ibinToCompare!=-1){
215 TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",700,500);
216 CratioComb->Divide(3,2);
217 TString nameComb[6]={Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
218 Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
219 TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120412.root");
220 TH1F *Spectra_copy[6]=0x0;
221 for(Int_t icharge=0;icharge<2;icharge++){
222 for(Int_t ipart=0;ipart<3;ipart++){
223 Int_t index=ipart+3*icharge;
224 TH1F *htmp=(TH1F*)((TH1F*)Spectra[index])->Clone("");
225 htmp->GetXaxis()->SetRangeUser(0,2.5);
226 TH1F *hcomb=fComb->Get(nameComb[index].Data())->Clone();
227 CratioComb->cd(ipart+1);
230 if(icharge==0)htmp->DrawClone();
231 else htmp->DrawClone("same");
232 hcomb->DrawClone("same");
234 htmp->SetMaximum(1.3);
235 htmp->SetMinimum(0.7);
236 CratioComb->cd(ipart+4);
239 if(icharge==0)htmp->DrawClone();
240 else htmp->DrawClone("same");
245 //comparison with charged hadron
246 Printf("\n\n-> ChargedHadron comparison");
247 TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
249 TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
250 //fraction of sec in MC
251 TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
252 TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
253 for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
254 Double_t en_data=hChHad_data->GetBinContent(ibin);
255 Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
256 Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
257 if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.1/en_mc));
258 Printf("Before: %.0f After: %.0f fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
260 hPrimRec_mc->Divide(hAllRec_mc);
261 //efficiency for primaries
262 TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
263 hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->Clone());
264 //Printf("--------%f ",((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
265 hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
266 hChHad_data->Divide(hEff_mc);//Efficiency
267 hChHad_data->Scale(1./(2*tcuts_data->GetEta()));
268 hChHad_data->SetTitle("All Ch from AOD");
272 hChHad_data->DrawClone();
273 TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
274 TH1F *hCh=fCh->Get("hpt_c0_5");
275 hCh->SetTitle("All Ch from Jacek");
276 hCh->SetMarkerColor(2);
277 hCh->SetLineColor(2);
279 for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
280 hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
281 hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
283 hCh->DrawClone("same");
285 TH1F *gRatio=AliPWGHistoTools::MyDivideHistosDifferentBins(hChHad_data,hCh);
286 gRatio->SetMaximum(1.3);
287 gRatio->SetMinimum(.7);
291 gRatio->DrawClone("");
293 TCanvas *cFitChargHad=new TCanvas("cFitChargHad","cFitChargHad",700,500);
296 hChHad_data->DrawClone();
297 //Fitting the sum of all particles
298 AliPWGFunc * fm = new AliPWGFunc;
299 fm->SetVarType(AliPWGFunc::kdNdpt);
300 Float_t fitmin = 0.3;
304 func = fm->GetBGBW(0.13,0.6,0.3, 1, 1e5);
305 func->SetParLimits(1, 0.1, 0.99);
306 func->SetParLimits(2, 0.01, 1);
307 func->SetParLimits(3, 0.01, 2);
308 TH1F * hToFit = hChHad_data;
309 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
310 // // if(!AliPWGHistoTools::Fit(hToFit,func,fitmin,fitmax)) {
311 // // cout << " FIT ERROR " << endl;
314 Double_t yieldTools, yieldETools;
315 Double_t partialYields[3],partialYieldsErrors[3];
316 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,0, 100, partialYields,partialYieldsErrors);
317 func->DrawClone("same");
318 Printf("TOTAL YIELD (AOD Charged Hadron) : %f +- %f",yieldTools,yieldETools);
321 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
322 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,0, 100, partialYields,partialYieldsErrors);
323 func->SetLineColor(2);
324 hCh->DrawClone("same");
325 func->DrawClone("same");
327 Printf("TOTAL YIELD (JACEK): %f +- %f",yieldTools,yieldETools);
329 //Convert spectra to dNdeta and sum
330 TH1F * hsum = (TH1F*) Spectra[0]->Clone();
332 Double_t epsilon = 0.0001;
333 for(Int_t icharge = 0; icharge < 2; icharge++){
334 for(Int_t ipart = 0; ipart < 3; ipart++){
335 Int_t index=ipart+3*icharge;
336 TH1F *htmp =(TH1F*)Spectra[index]->Clone("htmp");
337 Int_t nbin = htmp->GetNbinsX();
338 for(Int_t ibin = 1; ibin <= nbin; ibin++){
339 Double_t pt = htmp->GetBinCenter(ibin);
340 Double_t eta=0.8;//////////////////eta range///////////////////////////////////////
341 Double_t jacobian =eta2y(pt,mass[ipart],eta)/eta;
342 //Printf("jacobian: %f pt:%f BinContent:%f mass:%f",jacobian,pt,htmp->GetBinContent(ibin),mass[ipart]);
343 htmp->SetBinContent(ibin,htmp->GetBinContent(ibin)*jacobian);
344 htmp->SetBinError(ibin,htmp->GetBinError(ibin)*jacobian);
345 Int_t ibinSum = hsum->FindBin(pt);
346 if ( htmp->GetBinContent(ibin) > 0 &&
347 (TMath::Abs(htmp->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
348 TMath::Abs(htmp->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
350 cout << "DISCREPANCY IN BIN RANGES" << endl;
351 cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
352 << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
353 << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
360 hsum->SetTitle("Sum ID from AOD");
361 TCanvas *cChargHadComp=new TCanvas("cChargHadComp","cChargHadComp",700,500);
362 cChargHadComp->Divide(1,2);
363 cChargHadComp->cd(1);
368 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
369 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,0, 100, partialYields,partialYieldsErrors);
370 func->SetLineColor(2);
371 Printf("TOTAL YIELD (Pi+K+p): %f +- %f",yieldTools,yieldETools);
372 hChHad_data->SetMarkerColor(2);
373 hChHad_data->SetLineColor(2);
374 hChHad_data->DrawClone("same");
376 cChargHadComp->cd(2);
379 hsum->Divide(hChHad_data);
380 hsum->SetMaximum(1.2);
381 hsum->SetMinimum(.8);
385 // //Comparison of efficiency with TPCTOF ESD analysis
386 // Printf("\n\n-> Calculating Efficiency to be compared with ESD analysis");
387 // TH1F *EffTRUEPions;
388 // TH1F *EffSIGMAPions;
389 // TCanvas *cEffESD=new TCanvas("cEffESD","cEffESD",700,500);
390 // cEffESD->Divide(1,2);
394 // TString hname=Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
395 // Printf("Getting %s",hname.Data());
396 // EffTRUEPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
397 // EffTRUEPions->SetName(Form("Eff TRUE_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
398 // EffTRUEPions->SetTitle(Form("Eff TRUE %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
399 // EffTRUEPions->SetMarkerStyle(Marker[icharge]);
400 // EffTRUEPions->SetMarkerColor(Color[ipart]);
401 // EffTRUEPions->SetLineColor(Color[ipart]);
402 // hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
403 // Printf("... and divide it by %s",hname.Data());
404 // EffTRUEPions->Divide(EffTRUEPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
406 // hname=Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
407 // Printf("Getting %s",hname.Data());
408 // EffSIGMAPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
409 // EffSIGMAPions->SetName(Form("Eff SIGMA_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
410 // EffSIGMAPions->SetTitle(Form("Eff SIGMA %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
411 // EffSIGMAPions->SetMarkerStyle(Marker[icharge]);
412 // EffSIGMAPions->SetMarkerColor(Color[ipart+1]);
413 // EffSIGMAPions->SetLineColor(Color[ipart+1]);
414 // hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
415 // Printf("... and divide it by %s",hname.Data());
416 // EffSIGMAPions->Divide(EffSIGMAPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
418 // EffTRUEPions->DrawClone("lhist");
419 // EffSIGMAPions->DrawClone("lhistsame");
420 // hEffESD->DrawClone("lhistsame");
421 // gPad->BuildLegend();
423 // TH1F *hRatioTRUE=AliPWGHistoTools::MyDivideHistosDifferentBins(EffTRUEPions,hEffESD);
424 // hRatioTRUE->DrawClone("lhist");
425 // TH1F *hRatioSIGMA=AliPWGHistoTools::MyDivideHistosDifferentBins(EffSIGMAPions,hEffESD);
426 // hRatioSIGMA->DrawClone("lhistsame");
434 void DCACorrection(TH1F **Spectra, AliSpectraAODHistoManager* hman_data, AliSpectraAODHistoManager* hman_mc,Bool_t UseMCDCACorrection){
435 printf("\n\n-> DCA Correction");
437 Double_t FitRange[2]={-1.,1.};
440 TH1F *hcorrection[2];
441 TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
442 TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500);
443 cRatiocorrection->Divide(2,1);
444 ccorrection->Divide(2,1);
445 TString sample[2]={"data","mc"};
446 for(Int_t icharge=0;icharge<2;icharge++){
447 for(Int_t ipart=0;ipart<3;ipart++){
448 Int_t index=ipart+3*icharge;
449 for(Int_t isample=0;isample<2;isample++){
450 TCanvas *cDCA=new TCanvas(Form("cDCA%d%s",index,sample[isample].Data()),Form("cDCA%d%s",index,sample[isample].Data()),700,500);
451 hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
452 hcorrection[isample]->Reset("all");
454 for(Int_t ibin_data=6;ibin_data<38;ibin_data++){
455 Double_t lowedge=Spectra[index]->GetBinLowEdge(ibin_data);
456 Double_t binwidth=Spectra[index]->GetBinWidth(ibin_data);
457 if(isample==0)TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
458 if(isample==1)TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
459 TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
460 TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
461 TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
462 //Double_t minentries=200;
463 //if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)continue;
467 hToFit->Rebin(nrebin);
468 //Data and MC can have different stat
473 hToFit->Scale(1./hToFit->GetEntries());
474 Double_t normMC=hmc1->GetEntries()+hmc2->GetEntries()+hmc3->GetEntries();
475 hmc1->Scale(1./normMC);
476 hmc2->Scale(1./normMC);
477 hmc3->Scale(1./normMC);
478 cDCA->cd(ibin_data-5);
482 hToFit->DrawClone("lhist");
483 hmc3->DrawClone("lhist");
484 TObjArray *mc = new TObjArray(3); // MC histograms are put in this array
488 TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
489 fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
490 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
491 fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1
492 fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
493 hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
494 hToFit->SetTitle(Form("DCA distr - %s",sample[isample].Data()));
495 Int_t status = fit->Fit(); // perform the fit
496 cout << "fit status: " << status << endl;
497 if (status == 0) { // check on fit status
498 TH1F* result = (TH1F*) fit->GetPlot();
499 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
500 TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
501 TH1F* secMCPred=(TH1F*)fit->GetMCPrediction(2);
503 Double_t v1=0,v2=0,v3=0;
504 Double_t ev1=0,ev2=0,ev3=0;
505 //first method, use directly the fit result
506 fit->GetResult(0,v1,ev1);
507 fit->GetResult(1,v2,ev2);
508 fit->GetResult(2,v3,ev3);
509 result->Scale(hToFit->GetSumOfWeights()/result->GetSumOfWeights());
510 PrimMCPred->Scale(hToFit->GetSumOfWeights()*v1/PrimMCPred->GetSumOfWeights());
511 secStMCPred->Scale(hToFit->GetSumOfWeights()*v2/secStMCPred->GetSumOfWeights());
512 secMCPred->Scale(hToFit->GetSumOfWeights()*v3/secMCPred->GetSumOfWeights());
513 //second method, integrated the MC predisction, it should give the same as the first method
514 // v1=PrimMCPred->Integral();
515 // v2=secStMCPred->Integral();
516 // v3=secMCPred->Integral();
517 //Printf("\n\n\n\n\nv1:%f v2:%f v3:%f ev1:%f ev2:%f ev3:%f ",v1,v2,v3,ev1,ev2,ev3);
518 hcorrection[isample]->SetBinContent(ibin_data,v1/(v1+v2+v3));
519 hcorrection[isample]->SetBinError(ibin_data,0);
521 PrimMCPred->SetLineColor(2);
522 secStMCPred->SetLineColor(6);
523 secMCPred->SetLineColor(4);
524 hToFit->SetMinimum(0.0001);
525 hToFit->DrawClone("lhist");
526 result->SetTitle("Fit result");
527 result->DrawClone("lhistsame");
528 PrimMCPred->DrawClone("lhistsame");
529 secStMCPred->DrawClone("lhistsame");
530 secMCPred->DrawClone("lhistsame");
533 hcorrection[isample]->SetBinContent(ibin_data,1);
534 hcorrection[isample]->SetBinError(ibin_data,0);
540 ccorrection->cd(icharge+1);
543 hcorrection[isample]->SetTitle(Form("DCA corr %s %s %s",Particle[ipart].Data(),Charge[icharge].Data(),sample[isample].Data()));
544 hcorrection[isample]->SetLineWidth(2);
545 hcorrection[isample]->SetLineColor(Color[ipart]);
546 hcorrection[isample]->SetLineStyle(isample+1);
547 hcorrection[isample]->SetMarkerColor(Color[ipart]);
548 hcorrection[isample]->GetXaxis()->SetRangeUser(0.2,2.5);
549 if(ipart==0 && isample==0)hcorrection[isample]->DrawClone("lhist");
550 else hcorrection[isample]->DrawClone("lhistsame");
551 // smooth the DCA correction
552 // TF1 *fitfun = new TF1("fitfun","[0]+[1]*x^[2]+[3]*x^[4]",0.2,2.5);
553 // hcorrection[isample]->Fit(fitfun,"WRN","N",0.35,2);
554 // fitfun->SetLineWidth(1.5);
555 // fitfun->SetLineColor(Color[ipart]);
556 // fitfun->SetLineStyle(isample);
557 // fitfun->DrawClone("same");
558 // for(Int_t ibin=1;ibin<30;ibin++){
559 // hcorrection[isample]->SetBinContent(ibin,fitfun->Eval(hcorrection[isample]->GetBinCenter(ibin)));
562 Spectra[index]->Multiply(hcorrection[0]);//multiply for data
563 Printf("DCACorrection for DATA used: Spectra[index]->Multiply(hcorrection[0])")
564 if(UseMCDCACorrection){
565 Spectra[index]->Divide(hcorrection[1]); //divide by Monte Carlo
566 Printf("DCACorrection for MC used: Spectra[index]->Divide(hcorrection[1]")
568 cRatiocorrection->cd(icharge+1);
571 hcorrection[0]->Divide(hcorrection[1]);
572 if(ipart==0)hcorrection[0]->DrawClone("lhist");
573 else hcorrection[0]->DrawClone("lhistsame");
585 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta){
586 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
587 return TMath::ASinH(pt / mt * TMath::SinH(eta));
589 ///////////////////////
590 void GFCorrection(TH1F **Spectra,AliSpectraAODTrackCuts* tcuts_data){
591 //Geant/Fluka Correction
592 Printf("\nGF correction for Kaons");
593 //Getting GF For Kaons in TPC
594 TGraph *gGFCorrectionKaonPlus=new TGraph();
595 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
596 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
597 TGraph *gGFCorrectionKaonMinus=new TGraph();
598 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
599 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
600 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
601 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
602 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
603 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
604 //getting GF func for Kaons with TOF
605 TF1 *fGFKPosTracking;
606 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
607 TF1 *fGFKNegTracking;
608 fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
609 TF1 *fGFKPosMatching;
610 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
611 TF1 *fGFKNegMatching;
612 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
613 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
614 if(Spectra[1]->GetBinCenter(binK)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
615 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
616 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
617 Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
618 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
619 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
620 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
621 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
622 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
623 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
625 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
626 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
627 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
628 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
629 Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
630 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
631 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
632 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
633 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
634 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
635 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
636 Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
637 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
638 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
639 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
640 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
644 //Geant Fluka for P in TPC
645 Printf("\nGF correction for Protons");
646 const Int_t kNCharge=2;
649 TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
650 TH2D * hCorrFluka[kNCharge];
651 TH2D * hCorrFluka[2];
652 hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
653 hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
654 TGraph *gGFCorrectionProtonPlus=new TGraph();
655 gGFCorrectionProtonPlus->SetName("gGFCorrectionProtonPlus");
656 gGFCorrectionProtonPlus->SetTitle("gGFCorrectionProtonPlus");
657 TGraph *gGFCorrectionProtonMinus=new TGraph();
658 gGFCorrectionProtonMinus->SetName("gGFCorrectionProtonMinus");
659 gGFCorrectionProtonMinus->SetTitle("gGFCorrectionProtonMinus");
660 //getting GF func for Kaons with TPCTOF
661 TF1 *fGFpPosTracking;
662 fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
663 TF1 *fGFpNegTracking;
664 fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
665 TF1 *fGFpPosMatching;
666 fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
667 TF1 *fGFpNegMatching;
668 fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
670 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
671 Int_t nbins = Spectra[2]->GetNbinsX();
672 Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
673 for(Int_t ibin = 0; ibin < nbins; ibin++){
674 if(Spectra[2]->GetBinCenter(ibin)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
675 Float_t pt = Spectra[2]->GetBinCenter(ibin);
676 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
677 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
678 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
679 if (pt > maxPtCorrection) pt = maxPtCorrection;
680 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
681 //cout<<correction<<" charge "<<icharge<<endl;
683 if (correction != 0) {// If the bin is empty this is a 0
684 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
685 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError (ibin)*correction);
686 gGFCorrectionProtonPlus->SetPoint(ibin,pt,correction);
687 }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
688 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, Pt:"<< pt<< endl;
689 cout << " Bin content: " << Spectra[2]->GetBinContent(ibin) << endl;
693 if (correction != 0) {// If the bin is empty this is a 0
694 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
695 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError (ibin)*correction);
696 gGFCorrectionProtonMinus->SetPoint(ibin,pt,correction);
697 }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
698 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, Pt:"<< pt<< endl;
699 cout << " Bin content: " << Spectra[5]->GetBinContent(ibin) << endl;
703 gGFCorrectionProtonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
704 gGFCorrectionProtonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
705 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
706 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
707 Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
708 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrpPosTracking);
709 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrpNegTracking);
710 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrpPosTracking);
711 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrpNegTracking);
712 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
713 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
714 Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
715 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrpPosMatching);
716 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrpNegMatching);
717 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrpPosMatching);
718 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrpNegMatching);
722 gGFCorrectionKaonPlus->SetLineColor(kRed);
723 gGFCorrectionKaonMinus->SetLineColor(kRed+2);
724 gGFCorrectionProtonPlus->SetLineColor(kGreen);
725 gGFCorrectionProtonMinus->SetLineColor(kGreen+2);
726 fGFKPosTracking->SetLineColor(kRed);
727 fGFKNegTracking->SetLineColor(kRed+2);
728 fGFKPosMatching->SetLineColor(kRed);
729 fGFKNegMatching->SetLineColor(kRed+2);
730 fGFpPosTracking->SetLineColor(kGreen);
731 fGFpNegTracking->SetLineColor(kGreen+2);
732 fGFpPosMatching->SetLineColor(kGreen);
733 fGFpNegMatching->SetLineColor(kGreen+2);
734 fGFKPosTracking->SetLineStyle(2);
735 fGFKNegTracking->SetLineStyle(2);
736 fGFKPosMatching->SetLineStyle(3);
737 fGFKNegMatching->SetLineStyle(3);
738 fGFpPosTracking->SetLineStyle(2);
739 fGFpNegTracking->SetLineStyle(2);
740 fGFpPosMatching->SetLineStyle(3);
741 fGFpNegMatching->SetLineStyle(3);
742 fGFKPosTracking->SetRange(.6,5);
743 fGFKNegTracking->SetRange(.6,5);
744 fGFKPosMatching->SetRange(.6,5);
745 fGFKNegMatching->SetRange(.6,5);
746 fGFpPosTracking->SetRange(.6,5);
747 fGFpNegTracking->SetRange(.6,5);
748 fGFpPosMatching->SetRange(.6,5);
749 fGFpNegMatching->SetRange(.6,5);
751 TCanvas * GFCorrection = new TCanvas ("GFCorrection","GFCorrection",700,500);
754 gGFCorrectionKaonPlus->DrawClone("al");
755 gGFCorrectionKaonMinus->DrawClone("lsame");
756 gGFCorrectionProtonPlus->DrawClone("lsame");
757 gGFCorrectionProtonMinus->DrawClone("lsame");
758 fGFKPosTracking->DrawClone("lsame");
759 fGFKNegTracking->DrawClone("lsame");
760 fGFKPosMatching->DrawClone("lsame");
761 fGFKNegMatching->DrawClone("lsame");
762 fGFpPosTracking->DrawClone("lsame");
763 fGFpNegTracking->DrawClone("lsame");
764 fGFpPosMatching->DrawClone("lsame");
765 fGFpNegMatching->DrawClone("lsame");
773 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
776 if (ipart == 3 && icharge == kNegative) {
777 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
780 else if (ipart == 4 && icharge == kNegative) {
781 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
784 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
790 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
796 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
798 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
802 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
804 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
806 ///////////////////////////////////////////
808 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
811 if (ipart == 3 && icharge == kNegative) {
812 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
815 else if (ipart == 4 && icharge == kNegative) {
816 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
819 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
826 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
832 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
834 Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
835 return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
839 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
841 Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
842 return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);