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