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