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