]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TestAOD/MainAnalysis.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates in PID code and QA (ChiaraB, Francesco, Pietro, Rossella)
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / MainAnalysis.C
CommitLineData
6da07d74 1class AliPID;
2#include "QAPlots.C"
624a6bb0 3
4TString Charge[]={"Pos","Neg"};
5TString Sign[]={"Plus","Minus"};
6TString Particle[]={"Pion","Kaon","Proton"};
7Int_t Color[3]={1,2,4};
8Int_t Marker[6]={20,21,22,24,25,26};
9Double_t projl[3]={0,20,80};
10Double_t proju[3]={10,40,90};
6da07d74 11Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
12enum ECharge_t {
13 kPositive,
14 kNegative,
15 kNCharges
16};
624a6bb0 17
18
6da07d74 19void MainAnalysis() {
624a6bb0 20
acef4f19 21 gSystem->Load("libCore.so");
6da07d74 22 //gSystem->Load("libGeom.so");
624a6bb0 23 gSystem->Load("libPhysics.so");
6da07d74 24 //gSystem->Load("libVMC");
acef4f19 25 gSystem->Load("libTree");
6da07d74 26 //gSystem->Load("libProof");
acef4f19 27 gSystem->Load("libMatrix");
28 gSystem->Load("libSTEERBase");
29 gSystem->Load("libESD");
30 gSystem->Load("libAOD");
624a6bb0 31 gSystem->Load("libANALYSIS");
acef4f19 32 gSystem->Load("libOADB");
624a6bb0 33 gSystem->Load("libANALYSISalice");
acef4f19 34 gSystem->Load("libTENDER");
35 gSystem->Load("libCORRFW");
36 //gSystem->Load("libPWG0base");
6da07d74 37 //gSystem->Load("libMinuit");
acef4f19 38 gSystem->Load("libPWGTools");
39 gSystem->Load("libPWGLFSPECTRA");
624a6bb0 40 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
6da07d74 41
42 // Set Masses
43 Double_t mass[3];
44 mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
45 mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
46 mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
47
ae0fdd7d 48 //TString fold="3SigmaPID_AOD048-049_FilterBit6";
49 //TString fold="3SigmaPID_AOD048-049_FilterBit5";
50 TString fold="3SigmaPID_AOD048-049_FilterBit5";
51 //TString fold="3SigmaPID_AOD048-049_FilterBit5_test";
6da07d74 52 //TString fold="3SigmaPID_AOD086-090_FilterBit10";
6da07d74 53 Int_t ibinToCompare=-1;
54
55 TString sname="Cent0to5_QVec0.0to100.0";ibinToCompare=0;
56
57 TString dataFile = Form("output/%s/Pt.AOD.1._data_ptcut_%s.root",fold.Data(),sname.Data());
58 TString mcFile =Form("output/%s/Pt.AOD.1._MC_%s.root",fold.Data(),sname.Data());
624a6bb0 59 gStyle->SetPalette(1);
60
61 // Open root MC file and get classes
62 cout << "Analysis Macro" << endl;
63 cout << " > Reading MC data" << endl;
64 TFile *_mc = TFile::Open(mcFile.Data());
65 AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
66 AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
67 AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
68 // print info about mc track and Event cuts
69 cout << " -- Info about MC -- "<< endl;
70 ecuts_mc->PrintCuts();
71 // tcuts_mc->PrintCuts();
72 // proceed likewise for data
73 TFile *_data = TFile::Open(dataFile.Data());
74 AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
75 AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
76 AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
77 // print info about track and Event cuts
78 cout << " -- Info about data -- " << endl;
79 ecuts_data->PrintCuts();
80 tcuts_data->PrintCuts();
81
ae0fdd7d 82 QAPlots(hman_data,hman_mc,ecuts_data,ecuts_mc,tcuts_data,tcuts_mc);
83 return;
624a6bb0 84 //efficiencies
85 Printf("\n\n-> Calculating MC Correction Factors");
86 TH1F *CorrFact[6];
87 TCanvas *ceff=new TCanvas("ceff","ceff",700,500);
ed06106e 88 Bool_t UseMCDCACorrection=kTRUE;
624a6bb0 89 for(Int_t icharge=0;icharge<2;icharge++){
90 for(Int_t ipart=0;ipart<3;ipart++){
91 Int_t index=ipart+3*icharge;
ed06106e 92 //BE CAREFUL! depending on the efficiency you choose, you must change the DCA correction (data or data/mc)
93
6da07d74 94 TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); //MC correction for Prim+Cont+Eff
ae0fdd7d 95 //TString hname=Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());//MC correction for Prim+Eff if (idRec == idGen)
96 //TString hname=Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; //MC correction for Cont+Eff. filled with idRec
97 //TString hname=Form("hHistPtRecPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; //MC correction for Cont+Eff. filled with idGen
98 //TString hname=Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; // Pure MC efficiency for Prim. BE CAREFUL WITH MUONS!!! (idRec == idGen)
624a6bb0 99 Printf("Getting %s",hname.Data());
100 CorrFact[index]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
101 CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
102 CorrFact[index]->SetTitle(Form("CorrFact %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
103 CorrFact[index]->SetMarkerStyle(Marker[index]);
104 CorrFact[index]->SetMarkerColor(Color[ipart]);
105 CorrFact[index]->SetLineColor(Color[ipart]);
6da07d74 106 hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
624a6bb0 107 Printf("... and divide it by %s",hname.Data());
108 CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
109 if(index==0)CorrFact[index]->DrawClone();
110 else CorrFact[index]->DrawClone("same");
111 }
112 }
6da07d74 113
acef4f19 114 TFile *fESD=new TFile("EffAlex/pionEffPbPb.root");
115 TH1F *hEffESD=(TH1F*)fESD->Get("effMapPionTpcOnlyNeg0");
116 hEffESD->DrawClone("same");
624a6bb0 117 gPad->BuildLegend();
ed06106e 118
6da07d74 119 //Normalization
624a6bb0 120 printf("\n\n-> Spectra Normalization");
121 AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
122 Double_t events_data = ecuts_data->NumberOfEvents();
123 Printf(": accepted events: %.0f",events_data);
6da07d74 124 Double_t events_mc = ecuts_mc->NumberOfEvents();
125 Printf(": accepted events (MC): %.0f",events_mc);
624a6bb0 126
624a6bb0 127 //divide RAW for Correction Factor
128 Printf("\n\n-> Using MC correction factor to correct RAW spectra");
129 TH1F *Spectra[6];
130 for(Int_t icharge=0;icharge<2;icharge++){
131 for(Int_t ipart=0;ipart<3;ipart++){
132 Int_t index=ipart+3*icharge;
6da07d74 133 TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
624a6bb0 134 printf("Getting %s",hname.Data());
135 Spectra[index] =(TH1F*)((TH1F*) hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
136 Spectra[index]->SetName(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
137 Spectra[index]->SetTitle(Form("Spectra %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
138 Spectra[index]->SetMarkerStyle(Marker[index]);
139 Spectra[index]->SetMarkerColor(Color[ipart]);
140 Spectra[index]->SetLineColor(Color[ipart]);
141 Printf("... and divide it by %s",hname.Data());
6da07d74 142 Spectra[index]->Divide(CorrFact[index]);
624a6bb0 143 Spectra[index]->Scale(1./events_data,"width");//NORMALIZATION
144 }
145 }
146
ae0fdd7d 147
148
149 //Put Bin Content = 0 for bin below the the range
6da07d74 150 for(Int_t icharge=0;icharge<2;icharge++){
151 for(Int_t ipart=0;ipart<3;ipart++){
152 Int_t index=ipart+3*icharge;
153 for(Int_t ibin=0;ibin<Spectra[index]->GetNbinsX();ibin++){
154 if(Spectra[index]->GetBinCenter(ibin)<Range[ipart]){
155 Spectra[index]->SetBinContent(ibin,0);
156 Spectra[index]->SetBinError(ibin,0);
157 }
624a6bb0 158 }
6da07d74 159 }
624a6bb0 160 }
ae0fdd7d 161
6da07d74 162 //DCA Correction with the "right" DCA sample
ae0fdd7d 163 //DCACorrection(Spectra,hman_data,hman_mc,UseMCDCACorrection);
ed06106e 164
6da07d74 165 //DCA Correction forcing loose DCA
166 // TString fold_LooseDCA="5SigmaPID_AOD046_FilterBit5";
167 // TString dataFile_LooseDCA = Form("output/%s/Pt.AOD.1._data_ptcut_%s.root",fold_LooseDCA.Data(),sname.Data());
168 // TString mcFile_LooseDCA =Form("output/%s/Pt.AOD.1._MC_%s.root",fold_LooseDCA.Data(),sname.Data());
169 // TFile *_mc_LooseDCA = TFile::Open(mcFile_LooseDCA.Data());
170 // AliSpectraAODHistoManager* hman_mc_LooseDCA = (AliSpectraAODHistoManager*) _mc_LooseDCA->Get("SpectraHistos");
171 // TFile *_data_LooseDCA = TFile::Open(dataFile_LooseDCA.Data());
172 // AliSpectraAODHistoManager* hman_data_LooseDCA = (AliSpectraAODHistoManager*) _data_LooseDCA->Get("SpectraHistos");
173 // DCACorrection(Spectra,hman_data_LooseDCA,hman_mc_LooseDCA,UseMCDCACorrection);
174
6da07d74 175 //GFCorrection
176 GFCorrection(Spectra,tcuts_data);
624a6bb0 177
178 TH1F *MCTruth[6];
179 for(Int_t icharge=0;icharge<2;icharge++){
180 for(Int_t ipart=0;ipart<3;ipart++){
181 Int_t index=ipart+3*icharge;
6da07d74 182 hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
183 MCTruth[index]=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
184 MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
185 MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
186 MCTruth[index]->Scale(1./events_mc,"width");//NORMALIZATION
624a6bb0 187 }
188 }
189
ae0fdd7d 190
191 //Matching efficiency in data and Monte Carlo
192 TCanvas *cMatchingEff=new TCanvas("MatchingEff","MatchingEff",700,500);
193 TH1F *hMatcEffPos_data=(TH1F*)tcuts_data->GetHistoNMatchedPos();
194 hMatcEffPos_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedPos());
195 hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
196 TH1F *hMatcEffNeg_data=(TH1F*)tcuts_data->GetHistoNMatchedNeg();
197 hMatcEffNeg_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedNeg());
198 hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
199 hMatcEffNeg_data->SetLineColor(2);
200 TH1F *hMatcEffPos_mc=(TH1F*)tcuts_mc->GetHistoNMatchedPos();
201 hMatcEffPos_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedPos());
202 hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
203 hMatcEffPos_mc->SetLineStyle(2);
204 TH1F *hMatcEffNeg_mc=(TH1F*)tcuts_mc->GetHistoNMatchedNeg();
205 hMatcEffNeg_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedNeg());
206 hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
207 hMatcEffNeg_mc->SetLineColor(2);
208 hMatcEffNeg_mc->SetLineStyle(2);
209 cMatchingEff->Divide(1,2);
210 cMatchingEff->cd(1);
211 gPad->SetGridy();
212 gPad->SetGridx();
213 hMatcEffPos_data->DrawClone("lhist");
214 hMatcEffNeg_data->DrawClone("lhistsame");
215 hMatcEffPos_mc->DrawClone("lhistsame");
216 hMatcEffNeg_mc->DrawClone("lhistsame");
217 gPad->BuildLegend();
218 hMatcEffPos_data->Divide(hMatcEffPos_mc);
219 hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
220 cMatchingEff->cd(2);
221 gPad->SetGridy();
222 gPad->SetGridx();
223 hMatcEffPos_data->DrawClone("lhist");
224 hMatcEffNeg_data->DrawClone("lhistsame");
225 TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",2.5,5);
226 hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
227 pol0MatchPos_data->DrawClone("same");
228 TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",2.5,5);
229 hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
230 pol0MatchNeg_data->SetLineColor(2);
231 pol0MatchNeg_data->DrawClone("same");
232 Float_t ScalingMatchingPos=pol0MatchPos_data->GetParameter(0);
233 Float_t ScalingMatchingNeg=pol0MatchNeg_data->GetParameter(0);
234
235 //Correction spectra for matching efficiency
236 //For the moment I'm using the inclusive correction
237 for(Int_t ipart=0;ipart<3;ipart++){
238 for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++){
239 Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
240 if(ptspectra<tcuts_data->GetPtTOFMatching())continue;
241 //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
242 //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
243 Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/ScalingMatchingPos));
244 Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/ScalingMatchingNeg));
245 }
246 }
247
624a6bb0 248 //Drawing Final Spectra
249 TCanvas *cspectra=new TCanvas("cspectra","cspectra",700,500);
6da07d74 250 gPad->SetGridy();
251 gPad->SetGridy();
252 gPad->SetLogy();
624a6bb0 253 for(Int_t icharge=0;icharge<2;icharge++){
254 for(Int_t ipart=0;ipart<3;ipart++){
255 Int_t index=ipart+3*icharge;
256 if(index==0)Spectra[index]->DrawClone();
257 else Spectra[index]->DrawClone("same");
624a6bb0 258 }
259 }
260 gPad->BuildLegend();
ed06106e 261
624a6bb0 262 //saving spectra
263 Printf("\n\n-> Saving spectra in Out%s_%s.root",sname.Data(),fold.Data());
264 TFile * fout=new TFile(Form("results/Out_%s_%s.root",sname.Data(),fold.Data()),"RECREATE");
265 for(Int_t icharge=0;icharge<2;icharge++){
266 for(Int_t ipart=0;ipart<3;ipart++){
267 Int_t index=ipart+3*icharge;
268 Spectra[index]->Write();
269 CorrFact[index]->Write();
270 MCTruth[index]->Write();
271 }
272 }
ed06106e 273
624a6bb0 274 //if Bin 0-5% with no cut ratio with combined analysis
624a6bb0 275 if(ibinToCompare!=-1){
276 TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",700,500);
277 CratioComb->Divide(3,2);
278 TString nameComb[6]={Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
279 Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
6da07d74 280 TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120412.root");
281 TH1F *Spectra_copy[6]=0x0;
624a6bb0 282 for(Int_t icharge=0;icharge<2;icharge++){
283 for(Int_t ipart=0;ipart<3;ipart++){
6da07d74 284 Int_t index=ipart+3*icharge;
285 TH1F *htmp=(TH1F*)((TH1F*)Spectra[index])->Clone("");
286 htmp->GetXaxis()->SetRangeUser(0,2.5);
624a6bb0 287 TH1F *hcomb=fComb->Get(nameComb[index].Data())->Clone();
288 CratioComb->cd(ipart+1);
289 gPad->SetGridy();
290 gPad->SetGridx();
ae0fdd7d 291 for(Int_t ibin=1;ibin<hcomb->GetNbinsX();ibin++)hcomb->SetBinError(ibin,0);
292
6da07d74 293 if(icharge==0)htmp->DrawClone();
294 else htmp->DrawClone("same");
ae0fdd7d 295 MCTruth[index]->DrawClone("same");
296 hcomb->DrawClone("same");
297 htmp->Divide(hcomb);
6da07d74 298 htmp->SetMaximum(1.3);
299 htmp->SetMinimum(0.7);
624a6bb0 300 CratioComb->cd(ipart+4);
301 gPad->SetGridy();
302 gPad->SetGridx();
6da07d74 303 if(icharge==0)htmp->DrawClone();
304 else htmp->DrawClone("same");
624a6bb0 305 }
306 }
307 }
ed06106e 308
acef4f19 309 //comparison with charged hadron
310 Printf("\n\n-> ChargedHadron comparison");
6da07d74 311 TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
ae0fdd7d 312 TH1F* hMatchCorrectionAllCh=(TH1F*)hMatcEffPos_data->Clone("hMatchCorrectionAllCh");
313 hMatchCorrectionAllCh->Add(hMatcEffNeg_data); //correction for Matching efficiency!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
314 hMatchCorrectionAllCh->Scale(0.5);
315 for(Int_t ibin=1;ibin<hChHad_data->GetNbinsX();ibin++){
316 Float_t ptch=hChHad_data->GetBinCenter(ibin);
317 if(ptch<tcuts_data->GetPtTOFMatching())continue;
318 //hChHad_data->SetBinContent(ibin,(hChHad_data->GetBinContent(ibin)/hMatchCorrectionAllCh->GetBinContent(hMatchCorrectionAllCh->FindBin(ptch))));
319 hChHad_data->SetBinContent(ibin,2*(hChHad_data->GetBinContent(ibin)/(ScalingMatchingPos+ScalingMatchingNeg)));
320 }
acef4f19 321 //fraction of sec in MC
6da07d74 322 TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
323 TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
acef4f19 324 for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
ed06106e 325 Double_t en_data=hChHad_data->GetBinContent(ibin);
326 Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
327 Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
ae0fdd7d 328 if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.2/en_mc));
f6a38178 329 //Printf("Before: %.0f After: %.0f fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
acef4f19 330 }
ed06106e 331 hPrimRec_mc->Divide(hAllRec_mc);
acef4f19 332 //efficiency for primaries
6da07d74 333 TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
334 hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->Clone());
f6a38178 335 TCanvas *cAllChFactors=new TCanvas("cAllChFactors","cAllChFactors",700,500);
336 cAllChFactors->Divide(1,2);
337 cAllChFactors->cd(1);
338 gPad->SetGridy();
339 gPad->SetGridx();
340 hPrimRec_mc->SetTitle("Prim/All, charged hadron pure MC");
341 hPrimRec_mc->DrawClone("lhist");
342 gPad->BuildLegend();
343 cAllChFactors->cd(2);
344 gPad->SetGridy();
345 gPad->SetGridx();
346 hEff_mc->SetTitle("Efficiency for Primaries, charged hadron pure MC");
347 hEff_mc->DrawClone("lhist");
348 gPad->BuildLegend();
6da07d74 349 //Printf("--------%f ",((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
ed06106e 350 hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
351 hChHad_data->Divide(hEff_mc);//Efficiency
352 hChHad_data->Scale(1./(2*tcuts_data->GetEta()));
6da07d74 353 hChHad_data->SetTitle("All Ch from AOD");
f6a38178 354 TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
355 cAllCh->Divide(1,2);
6da07d74 356 cAllCh->cd(1);
ed06106e 357 gPad->SetGridy();
358 gPad->SetGridx();
359 hChHad_data->DrawClone();
360 TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
361 TH1F *hCh=fCh->Get("hpt_c0_5");
6da07d74 362 hCh->SetTitle("All Ch from Jacek");
363 hCh->SetMarkerColor(2);
364 hCh->SetLineColor(2);
ed06106e 365 //invariant yield
366 for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
367 hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
368 hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
369 }
370 hCh->DrawClone("same");
6da07d74 371 gPad->BuildLegend();
ed06106e 372 TH1F *gRatio=AliPWGHistoTools::MyDivideHistosDifferentBins(hChHad_data,hCh);
6da07d74 373 gRatio->SetMaximum(1.3);
374 gRatio->SetMinimum(.7);
375 cAllCh->cd(2);
ed06106e 376 gPad->SetGridy();
377 gPad->SetGridx();
378 gRatio->DrawClone("");
6da07d74 379 //fitting
380 TCanvas *cFitChargHad=new TCanvas("cFitChargHad","cFitChargHad",700,500);
381 gPad->SetGridy();
382 gPad->SetGridx();
383 hChHad_data->DrawClone();
384 //Fitting the sum of all particles
385 AliPWGFunc * fm = new AliPWGFunc;
386 fm->SetVarType(AliPWGFunc::kdNdpt);
387 Float_t fitmin = 0.3;
388 Float_t fitmax = 3;
389 TF1 * func = 0;
390 Int_t normPar = 0;
391 func = fm->GetBGBW(0.13,0.6,0.3, 1, 1e5);
392 func->SetParLimits(1, 0.1, 0.99);
393 func->SetParLimits(2, 0.01, 1);
394 func->SetParLimits(3, 0.01, 2);
395 TH1F * hToFit = hChHad_data;
396 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
397 // // if(!AliPWGHistoTools::Fit(hToFit,func,fitmin,fitmax)) {
398 // // cout << " FIT ERROR " << endl;
399 // // return;
400 // // }
401 Double_t yieldTools, yieldETools;
402 Double_t partialYields[3],partialYieldsErrors[3];
403 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,0, 100, partialYields,partialYieldsErrors);
404 func->DrawClone("same");
405 Printf("TOTAL YIELD (AOD Charged Hadron) : %f +- %f",yieldTools,yieldETools);
406 //Fit All Charged
407 hToFit = hCh;
408 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
409 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,0, 100, partialYields,partialYieldsErrors);
410 func->SetLineColor(2);
411 hCh->DrawClone("same");
412 func->DrawClone("same");
413 gPad->BuildLegend();
414 Printf("TOTAL YIELD (JACEK): %f +- %f",yieldTools,yieldETools);
415 //sumID vs AllCh
416 //Convert spectra to dNdeta and sum
417 TH1F * hsum = (TH1F*) Spectra[0]->Clone();
418 hsum->Reset("all");
419 Double_t epsilon = 0.0001;
420 for(Int_t icharge = 0; icharge < 2; icharge++){
421 for(Int_t ipart = 0; ipart < 3; ipart++){
422 Int_t index=ipart+3*icharge;
423 TH1F *htmp =(TH1F*)Spectra[index]->Clone("htmp");
424 Int_t nbin = htmp->GetNbinsX();
425 for(Int_t ibin = 1; ibin <= nbin; ibin++){
426 Double_t pt = htmp->GetBinCenter(ibin);
427 Double_t eta=0.8;//////////////////eta range///////////////////////////////////////
428 Double_t jacobian =eta2y(pt,mass[ipart],eta)/eta;
429 //Printf("jacobian: %f pt:%f BinContent:%f mass:%f",jacobian,pt,htmp->GetBinContent(ibin),mass[ipart]);
430 htmp->SetBinContent(ibin,htmp->GetBinContent(ibin)*jacobian);
431 htmp->SetBinError(ibin,htmp->GetBinError(ibin)*jacobian);
432 Int_t ibinSum = hsum->FindBin(pt);
433 if ( htmp->GetBinContent(ibin) > 0 &&
434 (TMath::Abs(htmp->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
435 TMath::Abs(htmp->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
436 ) {
437 cout << "DISCREPANCY IN BIN RANGES" << endl;
438 cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
439 << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
440 << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
441 cout << "" << endl;
442 }
443 }
444 hsum->Add(htmp);
445 }
446 }
447 hsum->SetTitle("Sum ID from AOD");
448 TCanvas *cChargHadComp=new TCanvas("cChargHadComp","cChargHadComp",700,500);
449 cChargHadComp->Divide(1,2);
450 cChargHadComp->cd(1);
451 gPad->SetGridy();
452 gPad->SetGridx();
453 hsum->DrawClone();
454 hToFit = hsum;
455 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
456 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,0, 100, partialYields,partialYieldsErrors);
457 func->SetLineColor(2);
458 Printf("TOTAL YIELD (Pi+K+p): %f +- %f",yieldTools,yieldETools);
459 hChHad_data->SetMarkerColor(2);
460 hChHad_data->SetLineColor(2);
461 hChHad_data->DrawClone("same");
acef4f19 462 gPad->BuildLegend();
6da07d74 463 cChargHadComp->cd(2);
464 gPad->SetGridy();
465 gPad->SetGridx();
466 hsum->Divide(hChHad_data);
467 hsum->SetMaximum(1.2);
468 hsum->SetMinimum(.8);
469 hsum->DrawClone("");
ed06106e 470
acef4f19 471
624a6bb0 472}
473
474
475
6da07d74 476
ed06106e 477void DCACorrection(TH1F **Spectra, AliSpectraAODHistoManager* hman_data, AliSpectraAODHistoManager* hman_mc,Bool_t UseMCDCACorrection){
6da07d74 478 printf("\n\n-> DCA Correction");
624a6bb0 479
6da07d74 480 Double_t FitRange[2]={-1.,1.};
481 Int_t nrebin=2;
624a6bb0 482 Printf("\DCACorr");
483 TH1F *hcorrection[2];
484 TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
6da07d74 485 TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500);
486 cRatiocorrection->Divide(2,1);
624a6bb0 487 ccorrection->Divide(2,1);
488 TString sample[2]={"data","mc"};
489 for(Int_t icharge=0;icharge<2;icharge++){
490 for(Int_t ipart=0;ipart<3;ipart++){
491 Int_t index=ipart+3*icharge;
492 for(Int_t isample=0;isample<2;isample++){
493 TCanvas *cDCA=new TCanvas(Form("cDCA%d%s",index,sample[isample].Data()),Form("cDCA%d%s",index,sample[isample].Data()),700,500);
494 hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
495 hcorrection[isample]->Reset("all");
496 cDCA->Divide(8,4);
497 for(Int_t ibin_data=6;ibin_data<38;ibin_data++){
498 Double_t lowedge=Spectra[index]->GetBinLowEdge(ibin_data);
499 Double_t binwidth=Spectra[index]->GetBinWidth(ibin_data);
6da07d74 500 if(isample==0)TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
501 if(isample==1)TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
502 TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
503 TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
504 TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
ae0fdd7d 505 Double_t minentries=1;
6da07d74 506 //if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)continue;
624a6bb0 507 hmc1->Rebin(nrebin);
508 hmc2->Rebin(nrebin);
509 hmc3->Rebin(nrebin);
510 hToFit->Rebin(nrebin);
511 //Data and MC can have different stat
512 hToFit->Sumw2();
513 hmc1->Sumw2();
514 hmc2->Sumw2();
515 hmc3->Sumw2();
516 hToFit->Scale(1./hToFit->GetEntries());
517 Double_t normMC=hmc1->GetEntries()+hmc2->GetEntries()+hmc3->GetEntries();
518 hmc1->Scale(1./normMC);
519 hmc2->Scale(1./normMC);
520 hmc3->Scale(1./normMC);
521 cDCA->cd(ibin_data-5);
522 gPad->SetGridy();
523 gPad->SetGridx();
524 gPad->SetLogy();
525 hToFit->DrawClone("lhist");
526 hmc3->DrawClone("lhist");
527 TObjArray *mc = new TObjArray(3); // MC histograms are put in this array
528 mc->Add(hmc1);
529 mc->Add(hmc2);
530 mc->Add(hmc3);
531 TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
532 fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
533 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
534 fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1
535 fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
536 hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
537 hToFit->SetTitle(Form("DCA distr - %s",sample[isample].Data()));
538 Int_t status = fit->Fit(); // perform the fit
539 cout << "fit status: " << status << endl;
540 if (status == 0) { // check on fit status
541 TH1F* result = (TH1F*) fit->GetPlot();
542 TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
543 TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
544 TH1F* secMCPred=(TH1F*)fit->GetMCPrediction(2);
6da07d74 545
546 Double_t v1=0,v2=0,v3=0;
547 Double_t ev1=0,ev2=0,ev3=0;
548 //first method, use directly the fit result
549 fit->GetResult(0,v1,ev1);
550 fit->GetResult(1,v2,ev2);
551 fit->GetResult(2,v3,ev3);
552 result->Scale(hToFit->GetSumOfWeights()/result->GetSumOfWeights());
553 PrimMCPred->Scale(hToFit->GetSumOfWeights()*v1/PrimMCPred->GetSumOfWeights());
554 secStMCPred->Scale(hToFit->GetSumOfWeights()*v2/secStMCPred->GetSumOfWeights());
555 secMCPred->Scale(hToFit->GetSumOfWeights()*v3/secMCPred->GetSumOfWeights());
556 //second method, integrated the MC predisction, it should give the same as the first method
557 // v1=PrimMCPred->Integral();
558 // v2=secStMCPred->Integral();
559 // v3=secMCPred->Integral();
560 //Printf("\n\n\n\n\nv1:%f v2:%f v3:%f ev1:%f ev2:%f ev3:%f ",v1,v2,v3,ev1,ev2,ev3);
561 hcorrection[isample]->SetBinContent(ibin_data,v1/(v1+v2+v3));
562 hcorrection[isample]->SetBinError(ibin_data,0);
624a6bb0 563 //Drawing section
564 PrimMCPred->SetLineColor(2);
565 secStMCPred->SetLineColor(6);
566 secMCPred->SetLineColor(4);
6da07d74 567 hToFit->SetMinimum(0.0001);
624a6bb0 568 hToFit->DrawClone("lhist");
569 result->SetTitle("Fit result");
570 result->DrawClone("lhistsame");
571 PrimMCPred->DrawClone("lhistsame");
572 secStMCPred->DrawClone("lhistsame");
573 secMCPred->DrawClone("lhistsame");
624a6bb0 574 }
575 else{
576 hcorrection[isample]->SetBinContent(ibin_data,1);
577 hcorrection[isample]->SetBinError(ibin_data,0);
578 }
579 Printf("deleting");
580 delete hToFit;
581 }
ed06106e 582
624a6bb0 583 ccorrection->cd(icharge+1);
584 gPad->SetGridy();
585 gPad->SetGridx();
586 hcorrection[isample]->SetTitle(Form("DCA corr %s %s %s",Particle[ipart].Data(),Charge[icharge].Data(),sample[isample].Data()));
587 hcorrection[isample]->SetLineWidth(2);
588 hcorrection[isample]->SetLineColor(Color[ipart]);
589 hcorrection[isample]->SetLineStyle(isample+1);
590 hcorrection[isample]->SetMarkerColor(Color[ipart]);
591 hcorrection[isample]->GetXaxis()->SetRangeUser(0.2,2.5);
592 if(ipart==0 && isample==0)hcorrection[isample]->DrawClone("lhist");
593 else hcorrection[isample]->DrawClone("lhistsame");
594 // smooth the DCA correction
6da07d74 595 // TF1 *fitfun = new TF1("fitfun","[0]+[1]*x^[2]+[3]*x^[4]",0.2,2.5);
596 // hcorrection[isample]->Fit(fitfun,"WRN","N",0.35,2);
624a6bb0 597 // fitfun->SetLineWidth(1.5);
6da07d74 598 // fitfun->SetLineColor(Color[ipart]);
599 // fitfun->SetLineStyle(isample);
624a6bb0 600 // fitfun->DrawClone("same");
601 // for(Int_t ibin=1;ibin<30;ibin++){
6da07d74 602 // hcorrection[isample]->SetBinContent(ibin,fitfun->Eval(hcorrection[isample]->GetBinCenter(ibin)));
624a6bb0 603 // }
624a6bb0 604 }
605 Spectra[index]->Multiply(hcorrection[0]);//multiply for data
ed06106e 606 Printf("DCACorrection for DATA used: Spectra[index]->Multiply(hcorrection[0])")
607 if(UseMCDCACorrection){
608 Spectra[index]->Divide(hcorrection[1]); //divide by Monte Carlo
609 Printf("DCACorrection for MC used: Spectra[index]->Divide(hcorrection[1]")
610 }
6da07d74 611 cRatiocorrection->cd(icharge+1);
612 gPad->SetGridy();
613 gPad->SetGridx();
614 hcorrection[0]->Divide(hcorrection[1]);
615 if(ipart==0)hcorrection[0]->DrawClone("lhist");
616 else hcorrection[0]->DrawClone("lhistsame");
624a6bb0 617 }
618 }
619 ccorrection->cd(1);
620 gPad->BuildLegend();
621 ccorrection->cd(2);
622 gPad->BuildLegend();
624a6bb0 623}
624
625
6da07d74 626//////////////
627
628Double_t eta2y(Double_t pt, Double_t mass, Double_t eta){
629 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
630 return TMath::ASinH(pt / mt * TMath::SinH(eta));
631}
632///////////////////////
633void GFCorrection(TH1F **Spectra,AliSpectraAODTrackCuts* tcuts_data){
634 //Geant/Fluka Correction
635 Printf("\nGF correction for Kaons");
636 //Getting GF For Kaons in TPC
637 TGraph *gGFCorrectionKaonPlus=new TGraph();
638 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
639 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
640 TGraph *gGFCorrectionKaonMinus=new TGraph();
641 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
642 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
643 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
644 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
645 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
646 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
647 //getting GF func for Kaons with TOF
648 TF1 *fGFKPosTracking;
649 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
650 TF1 *fGFKNegTracking;
651 fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
652 TF1 *fGFKPosMatching;
653 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
654 TF1 *fGFKNegMatching;
655 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
656 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
657 if(Spectra[1]->GetBinCenter(binK)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
658 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
659 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
660 Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
661 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
662 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
663 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
664 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
665 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
666 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
667 }else{
668 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
669 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
670 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
671 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
672 Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
673 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
674 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
675 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
676 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
677 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
678 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
679 Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
680 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
681 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
682 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
683 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
684 }
685 }
686
687 //Geant Fluka for P in TPC
688 Printf("\nGF correction for Protons");
689 const Int_t kNCharge=2;
690 Int_t kPos=0;
691 Int_t kNeg=1;
692 TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
693 TH2D * hCorrFluka[kNCharge];
694 TH2D * hCorrFluka[2];
695 hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
696 hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
697 TGraph *gGFCorrectionProtonPlus=new TGraph();
698 gGFCorrectionProtonPlus->SetName("gGFCorrectionProtonPlus");
699 gGFCorrectionProtonPlus->SetTitle("gGFCorrectionProtonPlus");
700 TGraph *gGFCorrectionProtonMinus=new TGraph();
701 gGFCorrectionProtonMinus->SetName("gGFCorrectionProtonMinus");
702 gGFCorrectionProtonMinus->SetTitle("gGFCorrectionProtonMinus");
703 //getting GF func for Kaons with TPCTOF
704 TF1 *fGFpPosTracking;
705 fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
706 TF1 *fGFpNegTracking;
707 fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
708 TF1 *fGFpPosMatching;
709 fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
710 TF1 *fGFpNegMatching;
711 fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
712
713 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
714 Int_t nbins = Spectra[2]->GetNbinsX();
715 Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
716 for(Int_t ibin = 0; ibin < nbins; ibin++){
717 if(Spectra[2]->GetBinCenter(ibin)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
718 Float_t pt = Spectra[2]->GetBinCenter(ibin);
719 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
720 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
721 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
722 if (pt > maxPtCorrection) pt = maxPtCorrection;
723 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
724 //cout<<correction<<" charge "<<icharge<<endl;
725 if(icharge==0){
726 if (correction != 0) {// If the bin is empty this is a 0
727 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
728 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError (ibin)*correction);
729 gGFCorrectionProtonPlus->SetPoint(ibin,pt,correction);
730 }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
731 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, Pt:"<< pt<< endl;
732 cout << " Bin content: " << Spectra[2]->GetBinContent(ibin) << endl;
733 }
734 }
735 if(icharge==1){
736 if (correction != 0) {// If the bin is empty this is a 0
737 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
738 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError (ibin)*correction);
739 gGFCorrectionProtonMinus->SetPoint(ibin,pt,correction);
740 }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
741 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, Pt:"<< pt<< endl;
742 cout << " Bin content: " << Spectra[5]->GetBinContent(ibin) << endl;
743 }
744 }
745 }else{
ae0fdd7d 746 gGFCorrectionProtonPlus->SetPoint(ibin,Spectra[2]->GetBinCenter(ibin),0);
747 gGFCorrectionProtonMinus->SetPoint(ibin,Spectra[5]->GetBinCenter(ibin),0);
748 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
749 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
750 Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
751 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
752 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
753 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
754 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
755 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
756 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
757 Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
758 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
759 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
760 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
761 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
6da07d74 762 }
763 }//end loop on bins
764 }
765 gGFCorrectionKaonPlus->SetLineColor(kRed);
766 gGFCorrectionKaonMinus->SetLineColor(kRed+2);
767 gGFCorrectionProtonPlus->SetLineColor(kGreen);
768 gGFCorrectionProtonMinus->SetLineColor(kGreen+2);
769 fGFKPosTracking->SetLineColor(kRed);
770 fGFKNegTracking->SetLineColor(kRed+2);
771 fGFKPosMatching->SetLineColor(kRed);
772 fGFKNegMatching->SetLineColor(kRed+2);
773 fGFpPosTracking->SetLineColor(kGreen);
774 fGFpNegTracking->SetLineColor(kGreen+2);
775 fGFpPosMatching->SetLineColor(kGreen);
776 fGFpNegMatching->SetLineColor(kGreen+2);
777 fGFKPosTracking->SetLineStyle(2);
778 fGFKNegTracking->SetLineStyle(2);
779 fGFKPosMatching->SetLineStyle(3);
780 fGFKNegMatching->SetLineStyle(3);
781 fGFpPosTracking->SetLineStyle(2);
782 fGFpNegTracking->SetLineStyle(2);
783 fGFpPosMatching->SetLineStyle(3);
784 fGFpNegMatching->SetLineStyle(3);
785 fGFKPosTracking->SetRange(.6,5);
786 fGFKNegTracking->SetRange(.6,5);
787 fGFKPosMatching->SetRange(.6,5);
788 fGFKNegMatching->SetRange(.6,5);
789 fGFpPosTracking->SetRange(.6,5);
790 fGFpNegTracking->SetRange(.6,5);
791 fGFpPosMatching->SetRange(.6,5);
792 fGFpNegMatching->SetRange(.6,5);
793
794 TCanvas * GFCorrection = new TCanvas ("GFCorrection","GFCorrection",700,500);
795 gPad->SetGridx();
796 gPad->SetGridy();
797 gGFCorrectionKaonPlus->DrawClone("al");
798 gGFCorrectionKaonMinus->DrawClone("lsame");
799 gGFCorrectionProtonPlus->DrawClone("lsame");
800 gGFCorrectionProtonMinus->DrawClone("lsame");
801 fGFKPosTracking->DrawClone("lsame");
802 fGFKNegTracking->DrawClone("lsame");
803 fGFKPosMatching->DrawClone("lsame");
804 fGFKNegMatching->DrawClone("lsame");
805 fGFpPosTracking->DrawClone("lsame");
806 fGFpNegTracking->DrawClone("lsame");
807 fGFpPosMatching->DrawClone("lsame");
808 fGFpNegMatching->DrawClone("lsame");
809 gPad->BuildLegend();
810}
811
812
813
814///////////
815TF1 *
816TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
817{
818
819 if (ipart == 3 && icharge == kNegative) {
820 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
821 return f;
822 }
823 else if (ipart == 4 && icharge == kNegative) {
824 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
825 }
826 else
827 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
828
829 return f;
830}
831
832Double_t
833TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
834{
835 return 1.;
836}
837
838Double_t
839TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
840{
841 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
842}
843
844Double_t
845TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
846{
847 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
848}
849///////////////////////////////////////////
850TF1 *
851TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
852{
853
854 if (ipart == 3 && icharge == kNegative) {
855 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
856 return f;
857 }
858 else if (ipart == 4 && icharge == kNegative) {
859 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
860 }
861 else
862 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
863
864 return f;
865}
866
867
868Double_t
869MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
870{
871 return 1.;
872}
873
874Double_t
875MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
876{
877 Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
878 return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
879}
880
881Double_t
882MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
883{
884 Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
885 return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
886}