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