]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TestAOD/HighLevelQA/JCGAnalysis.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 / HighLevelQA / JCGAnalysis.C
CommitLineData
22d1b439 1/////////////////////////////////////////////////////
2// JCGAnalysis.C //
3// based on MainAnalysis.C //
4// uses a new method of calculating the efficiency //
5/////////////////////////////////////////////////////
6
7class AliPID;
8
9// global constants
10const Int_t nCharge = 2;
11const Int_t nPart = 3;
12TString Charge[nCharge] = {"Pos","Neg"};
13TString Sign[nCharge] = {"Plus","Minus"};
14TString Particle[nPart] = {"Pion","Kaon","Proton"};
15Int_t Color[nPart] = {1,2,4};
16Int_t Marker[nPart*nCharge] = {20,21,22,24,25,26};
17Double_t Range[nPart] = {0.3,0.3,0.5}; // LowPt range for pi k p
18enum ECharge_t {kPositive,kNegative,kNCharges};
19TString Names[nCharge*nPart] = {"#pi^{+}",
20 "K^{+}",
21 "p",
22 "#pi^{-}",
23 "K^{-}",
24 "#bar{p}"};
25
26void JCGAnalysis()
27{
28 // load libraries 'n such
4070f709 29 gSystem->Load("libCore");
30 gSystem->Load("libPhysics");
22d1b439 31 gSystem->Load("libTree");
32 gSystem->Load("libMatrix");
33 gSystem->Load("libSTEERBase");
34 gSystem->Load("libESD");
35 gSystem->Load("libAOD");
36 gSystem->Load("libANALYSIS");
37 gSystem->Load("libOADB");
38 gSystem->Load("libANALYSISalice");
af472fff 39 gSystem->Load("libTender");
22d1b439 40 gSystem->Load("libCORRFW");
41 gSystem->Load("libPWGTools");
42 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
43 gStyle->SetOptStat(000000);
44 gStyle->SetPalette(1);
45
46 // Set masses
47 Double_t mass[3];
48 mass[0] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
49 mass[1] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
50 mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
51
52 // cuts 'n such
53 // 0 1 2 3
54 Double_t CentCutMin[4]= { 0, 30, 30, 30};
55 Double_t CentCutMax[4]= { 5, 40, 40, 40};
56 Double_t QvecCutMin[4]={ 0, 0, 0, 1.5};
57 Double_t QvecCutMax[4]={ 100, 100, 0.4, 100};
58 Double_t EtaMin[4]={ -0.8,-0.8,-0.8,-0.8};
59 Double_t EtaMax[4]={ 0.8, 0.8, 0.8, 0.8};
60 Double_t Nsigmapid=3.;
61 Double_t pt=10.;
62 Double_t p=10.;
63 Double_t y=.5;
64 Double_t ptTofMatch=.6;
65 UInt_t trkbit=1024;
66 UInt_t trkbitQVector=1;
67 Bool_t UseCentPatchAOD049=kFALSE;
68 Double_t DCA=100000;
69 UInt_t minNclsTPC=70;
70 Int_t nrebin=0;
71 TString opt="";
72 Int_t icut=1; // chooses which cut to use
73 if(icut==0)Int_t ibinToCompare=0;
74 else Int_t ibinToCompare=4;
75
76 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);
77 TString sname_data_1sig=Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_1.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
78 // for MC we do not use the cut on Qvector
79 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);
80 TString sname_mc_5sig=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_5.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
81 TString sname_mc_1sig=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_1.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
82
83 // input and output files
84 TString fold="No_Outliers";
85 TString dataFile = Form("output/%s/AnResDATATrain12_NoOutliers.root",fold.Data());
86 TString mcFile = Form("output/%s/AnResMCTrain27_NoOutliers.root",fold.Data());
87 TFile * fout=new TFile(Form("results/%s/Res_%s_%s.root",fold.Data(),sname_data.Data(),fold.Data()),"RECREATE");
88
89 // Open root MC file and get classes
90 TFile *f_mc = TFile::Open(mcFile.Data());
91 TDirectoryFile *Dir_mc=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc.Data()));
92 AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) Dir_mc->Get("SpectraHistos");
93 AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) Dir_mc->Get("Event Cuts");
94 AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) Dir_mc->Get("Track Cuts");
95 TDirectoryFile *Dir_mc_5sig=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc_5sig.Data()));
96 AliSpectraAODHistoManager* hman_mc_5sig = (AliSpectraAODHistoManager*) Dir_mc_5sig->Get("SpectraHistos");
97 AliSpectraAODEventCuts* ecuts_mc_5sig = (AliSpectraAODEventCuts*) Dir_mc_5sig->Get("Event Cuts");
98 AliSpectraAODTrackCuts* tcuts_mc_5sig = (AliSpectraAODTrackCuts*) Dir_mc_5sig->Get("Track Cuts");
99 TDirectoryFile *Dir_mc_1sig=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc_1sig.Data()));
100 AliSpectraAODHistoManager* hman_mc_1sig = (AliSpectraAODHistoManager*) Dir_mc_1sig->Get("SpectraHistos");
101 AliSpectraAODEventCuts* ecuts_mc_1sig = (AliSpectraAODEventCuts*) Dir_mc_1sig->Get("Event Cuts");
102 AliSpectraAODTrackCuts* tcuts_mc_1sig = (AliSpectraAODTrackCuts*) Dir_mc_1sig->Get("Track Cuts");
103 // proceed likewise for data:
104 // Open root DATA file and get classes
105 TFile *f_data = TFile::Open(dataFile.Data());
106 TDirectoryFile *Dir_data=(TDirectoryFile*)f_data->Get(Form("%s",sname_data.Data()));
107 AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) Dir_data->Get("SpectraHistos");
108 AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) Dir_data->Get("Event Cuts");
109 AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) Dir_data->Get("Track Cuts");
110 TDirectoryFile *Dir_data_1sig=(TDirectoryFile*)f_data->Get(Form("%s",sname_data_1sig.Data()));
111 AliSpectraAODHistoManager* hman_data_1sig = (AliSpectraAODHistoManager*) Dir_data_1sig->Get("SpectraHistos");
112 AliSpectraAODEventCuts* ecuts_data_1sig = (AliSpectraAODEventCuts*) Dir_data_1sig->Get("Event Cuts");
113 AliSpectraAODTrackCuts* tcuts_data_1sig = (AliSpectraAODTrackCuts*) Dir_data_1sig->Get("Track Cuts");
114
115 //---------------------------------------------
116 // plot nsigma distributions for presentation -
117 //---------------------------------------------
118 // TPC
119 TH2F * hTPCnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTPC"))->Clone());
120 hTPCnsigPion->SetTitle("TPC NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
121 TCanvas * cTPCnsigPion = new TCanvas("cTPCnsigPion","cTPCnsigPion");
122 gPad->SetLogz();
123 hTPCnsigPion->GetXaxis()->SetRangeUser(0,5);
124 hTPCnsigPion->GetYaxis()->SetRangeUser(-10,40);
125 hTPCnsigPion->DrawCopy("COLZ");
126 // TOF
127 TH2F * hTOFnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTOF"))->Clone());
128 hTOFnsigPion->SetTitle("TOF NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
129 TCanvas * cTOFnsigPion = new TCanvas("cTOFnsigPion","cTOFnsigPion");
130 gPad->SetLogz();
131 hTOFnsigPion->GetYaxis()->SetRangeUser(-10,40);
132 hTOFnsigPion->GetXaxis()->SetRangeUser(0,5);
133 hTOFnsigPion->DrawCopy("COLZ");
134 // TPCTOF
135 TH2F * hTPCTOFnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTPCTOF"))->Clone());
136 hTPCTOFnsigPion->SetTitle("TPC+TOF NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
137 TCanvas * cTPCTOFnsigPion = new TCanvas("cTPCTOFnsigPion","cTPCTOFnsigPion");
138 gPad->SetLogz();
139 hTPCTOFnsigPion->GetYaxis()->SetRangeUser(-5,40);
140 hTPCTOFnsigPion->GetXaxis()->SetRangeUser(0,5);
141 hTPCTOFnsigPion->DrawCopy("COLZ");
142
143
144 //------------------------------------------------------
145 // plot the MC True Generated Primaries for comparison -
146 //------------------------------------------------------
147 TH1F * MCTruthSpectra[nCharge*nPart];
148 Double_t events_mc = ecuts_mc->NumberOfEvents();
149 for (Int_t icharge=0; icharge<nCharge; icharge++)
150 {
151 for (Int_t ipart=0; ipart<nPart; ipart++)
152 {
153 Int_t index = ipart + 3*icharge;
154 TString hname = Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
155 MCTruthSpectra[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1)->Clone());
156 MCTruthSpectra[index]->SetTitle("Generated Primaries Spectra (MC Truth);p_{T} (GeV/c);");
157 MCTruthSpectra[index]->SetMarkerStyle(Marker[index]);
158 MCTruthSpectra[index]->SetMarkerColor(Color[ipart]);
159 MCTruthSpectra[index]->SetLineColor(Color[ipart]);
160 MCTruthSpectra[index]->Scale(1./events_mc,"width");
161 }
162 }
163 // neg/pos ratios
164 TH1F * NegPosRatio[nPart];
165 for (Int_t ipart=0; ipart<nPart; ipart++)
166 {
167 NegPosRatio[ipart] = (TH1F*)MCTruthSpectra[ipart+nPart]->Clone();
168 NegPosRatio[ipart]->Divide(MCTruthSpectra[ipart]);
169 NegPosRatio[ipart]->SetTitle("Neg/Pos;p_{T} (GeV/c);");
170 }
171 // drawing
172 TCanvas * cMCTruthSpectra = new TCanvas("cMCTruthSpectra","cMCTruthSpectra");
173 cMCTruthSpectra->Divide(1,2);
174 TLegend * lMCTruthSpectra = new TLegend(.69,.69,.89,.99);
175 lMCTruthSpectra->SetFillColor(0);
176 TLegend * lPosNegRatios = new TLegend (.69,.69,.99,.99);
177 lPosNegRatios->SetFillColor(0);
178 for (Int_t icharge=0; icharge<nCharge; icharge++)
179 {
180 for (Int_t ipart=0; ipart<nPart; ipart++)
181 {
182 Int_t index = ipart + 3*icharge;
183 cMCTruthSpectra->cd(1);
184 gPad->SetLogy();
185 if (index == 0) MCTruthSpectra[index]->DrawCopy("EP");
186 else MCTruthSpectra[index]->DrawCopy("EPsame");
187 lMCTruthSpectra->AddEntry(MCTruthSpectra[index],Names[index].Data(),"p");
188
189 cMCTruthSpectra->cd(2);
190 NegPosRatio[ipart]->GetYaxis()->SetRangeUser(0.4,1.1);
191 if (index == 0) NegPosRatio[ipart]->DrawCopy("hist][");
192 else NegPosRatio[ipart]->DrawCopy("hist][same");
193 if (icharge == 0) lPosNegRatios->AddEntry(NegPosRatio[ipart],Particle[ipart].Data(),"l");
194 }
195 }
196 cMCTruthSpectra->cd(1);
197 lMCTruthSpectra->DrawClone();
198 cMCTruthSpectra->cd(2);
199 lPosNegRatios->DrawClone();
200
201
202 //-----------------------------------
203 // Get the (normalized) raw spectra -
204 //-----------------------------------
205 TH1F * Spectra[nCharge*nPart];
206 Double_t events_data = ecuts_data->NumberOfEvents();
207 Double_t events_data_1sig = ecuts_data_1sig->NumberOfEvents();
208 Double_t events_mc = ecuts_mc->NumberOfEvents();
209 Double_t events_mc_1sig = ecuts_mc_1sig->NumberOfEvents();
210 for (Int_t icharge=0; icharge<nCharge; icharge++)
211 {
212 for (Int_t ipart=0; ipart<nPart; ipart++)
213 {
214 Int_t index = ipart + 3*icharge;
215 TString hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
216 Spectra[index] = (TH1F*)((TH1F*)hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
217 Spectra[index]->SetTitle(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
218 Spectra[index]->SetMarkerStyle(Marker[index]);
219 Spectra[index]->SetMarkerColor(Color[ipart]);
220 Spectra[index]->SetLineColor(Color[ipart]);
221 Spectra[index]->Scale(1./events_data,"width"); // normalization
222 }
223 }
224 // Set bins outisde the range to contain zeros
225 for (Int_t icharge=0; icharge<nCharge; icharge++)
226 {
227 for (Int_t ipart=0; ipart<nPart; ipart++)
228 {
229 Int_t index = ipart + 3*icharge;
230 for (Int_t ibin=0; ibin<Spectra[index]->GetNbinsX(); ibin++)
231 {
232 if (Spectra[index]->GetBinCenter(ibin) < Range[ipart])
233 {
234 Spectra[index]->SetBinContent(ibin, 0);
235 Spectra[index]->SetBinError(ibin, 0);
236 }
237 }
238 }
239 }
240 // save the raw spectra for use later
241 TH1F * RawSpectra[nCharge*nPart];
242 for (Int_t icharge=0; icharge<nCharge; icharge++)
243 {
244 for (Int_t ipart=0; ipart<nPart; ipart++)
245 {
246 Int_t index = ipart + 3*icharge;
247 RawSpectra[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
248 RawSpectra[index]->SetTitle("Raw Spectra");
249 }
250 }
251
252 // plot the raw spectra from DATA
253 // also include allch
254 TCanvas * cRawSpectra = new TCanvas("cRawSpectra","cRawSpectra");
255 TLegend * lRawSpectra = new TLegend(.69,.69,.99,.99);
256 lRawSpectra->SetFillColor(0);
257 TH1F * hRawAllCh = (TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
258 hRawAllCh->Scale(1./ecuts_data->NumberOfEvents(),"width");
259 hRawAllCh->SetTitle("Raw Spectra, Data;p_{T} (GeV/c);");
260 hRawAllCh->SetMarkerColor(kGreen+3);
261 hRawAllCh->SetLineColor(kGreen+3);
262 hRawAllCh->SetMarkerStyle(34);
263 hRawAllCh->DrawCopy();
264 lRawSpectra->AddEntry(hRawAllCh,"All Particles","p");
265 for (Int_t icharge=0; icharge<nCharge; icharge++)
266 {
267 for (Int_t ipart=0; ipart<nPart; ipart++)
268 {
269 Int_t index = ipart + 3*icharge;
270 RawSpectra[index]->SetTitle("Raw Spectra, Data;p_{T} (GeV/c);");
271 RawSpectra[index]->SetMarkerColor(Color[ipart]);
272 RawSpectra[index]->SetLineColor(Color[ipart]);
273 RawSpectra[index]->SetMarkerStyle(Marker[index]);
274 RawSpectra[index]->DrawCopy("same");
275 lRawSpectra->AddEntry(RawSpectra[index],Names[index].Data(),"p");
276 }
277 }
278 lRawSpectra->DrawClone();
279
280 // also get MC raw for comparison purposes
281 TH1F * MCSpectra[nCharge*nPart];
282 for (Int_t icharge=0; icharge<nCharge; icharge++)
283 {
284 for (Int_t ipart=0; ipart<nPart; ipart++)
285 {
286 Int_t index = ipart + 3*icharge;
287 TString hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
288 MCSpectra[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
289 MCSpectra[index]->SetTitle(Form("MCSpectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
290 MCSpectra[index]->SetMarkerStyle(Marker[index]);
291 MCSpectra[index]->SetMarkerColor(Color[ipart]);
292 MCSpectra[index]->SetLineColor(Color[ipart]);
293 MCSpectra[index]->Scale(1./events_mc,"width"); // normalization
294 }
295 }
296 // Set bins outisde the range to contain zeros
297 for (Int_t icharge=0; icharge<nCharge; icharge++)
298 {
299 for (Int_t ipart=0; ipart<nPart; ipart++)
300 {
301 Int_t index = ipart + 3*icharge;
302 for (Int_t ibin=0; ibin<MCSpectra[index]->GetNbinsX(); ibin++)
303 {
304 if (MCSpectra[index]->GetBinCenter(ibin) < Range[ipart])
305 {
306 MCSpectra[index]->SetBinContent(ibin, 0);
307 MCSpectra[index]->SetBinError(ibin, 0);
308 }
309 }
310 }
311 }
312 // save the raw spectra for use later
313 TH1F * MCRawSpectra[nCharge*nPart];
314 for (Int_t icharge=0; icharge<nCharge; icharge++)
315 {
316 for (Int_t ipart=0; ipart<nPart; ipart++)
317 {
318 Int_t index = ipart + 3*icharge;
319 MCRawSpectra[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
320 MCRawSpectra[index]->SetTitle("MC Raw Spectra");
321 }
322 }
323
324
325 //----------------
326 // MC Correction -
327 //----------------
328 // calculate the percentage of primaries reconstructed w/ 5 sigma PID
329 Printf("\n\n-> Calculating MC Correction Factors");
330 TH1F * CorrFact[nPart*nCharge];
331 for (Int_t icharge=0; icharge<nCharge; icharge++)
332 {
333 for (Int_t ipart=0; ipart<nPart; ipart++)
334 {
335 Int_t index = ipart + 3*icharge;
336 TString hname = Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
337 Printf("Getting %s",hname.Data());
338 CorrFact[index] = new TH1F(*(TH1F*)((TH1F*) hman_mc_5sig->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
339 CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
340 CorrFact[index]->SetTitle("MC Correction Factor");
341 CorrFact[index]->SetMarkerStyle(Marker[index]);
342 CorrFact[index]->SetMarkerColor(Color[ipart]);
343 CorrFact[index]->SetLineColor(Color[ipart]);
344 hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
345 Printf("... and divide it by %s",hname.Data());
346 CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc_5sig->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
347 if (index == 0)
348 {
349 TCanvas * cCorrFactMC = new TCanvas("cCorrFactMC","cCorrFactMC");
350 TLegend * lCorrFactMC = new TLegend(.69,.69,.99,.99);
351 lCorrFactMC->SetFillColor(0);
352 CorrFact[index]->DrawCopy("EP");
353 }
354 else CorrFact[index]->DrawCopy("EPsame");
355 lCorrFactMC->AddEntry(CorrFact[index],Names[index].Data(),"p");
356 } // end loop over ipart
357 } // end loop over icharge
358 lCorrFactMC->DrawClone();
359 // Correct spectra with MC correction factor
360 for (Int_t icharge=0; icharge<nCharge; icharge++)
361 {
362 for (Int_t ipart=0; ipart<nPart; ipart++)
363 {
364 Int_t index = ipart + 3*icharge;
365 Spectra[index]->Divide(CorrFact[index]);
366 }
367 }
368 // save this version of the spectra for use later
369 TH1F * SpectraStep1[nCharge*nPart];
370 for (Int_t icharge=0; icharge<nCharge; icharge++)
371 {
372 for (Int_t ipart=0; ipart<nPart; ipart++)
373 {
374 Int_t index = ipart + 3*icharge;
375 SpectraStep1[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
376 SpectraStep1[index]->SetTitle("Spectra After MC Correction");
377 }
378 }
379 // do the same with MCSpectra
380 for (Int_t icharge=0; icharge<nCharge; icharge++)
381 {
382 for (Int_t ipart=0; ipart<nPart; ipart++)
383 {
384 Int_t index = ipart + 3*icharge;
385 MCSpectra[index]->Divide(CorrFact[index]);
386 }
387 }
388 // save this version of the spectra for use later
389 TH1F * MCSpectraStep1[nCharge*nPart];
390 for (Int_t icharge=0; icharge<nCharge; icharge++)
391 {
392 for (Int_t ipart=0; ipart<nPart; ipart++)
393 {
394 Int_t index = ipart + 3*icharge;
395 MCSpectraStep1[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
396 MCSpectraStep1[index]->SetTitle("MCSpectra After MC Correction");
397 }
398 }
399
400
401 //--------------------------------------
402 // Correct spectra for PID nsigma cuts -
403 //--------------------------------------
404 // for the TPC only, assume the nsigma distribution is purely gaussian, but account for the actual means and sigmas
405 // the correction factor is the integral of a normal distribution from -3sqrt(2) to 3sqrt(2)
406 TF1 * fGaussTPC[nPart];
407 for (Int_t ipart=0; ipart<nPart; ipart++)
408 {
409 fGaussTPC[ipart] = new TF1(Form("fGaussTPC%s",Particle[ipart].Data()),"gaus",-25,25);
410 fGaussTPC[ipart]->SetParameter(0,1); // set a default height for now - will be normalized later
411 }
412 // pions
413 fGaussTPC[0]->SetParameter(1,0.04);
414 fGaussTPC[0]->SetParameter(2,1.02);
415 // kaons
416 fGaussTPC[1]->SetParameter(1,-0.06);
417 fGaussTPC[1]->SetParameter(2,0.97);
418 // protons
419 fGaussTPC[2]->SetParameter(1,0.25);
420 fGaussTPC[2]->SetParameter(2,0.72);
421 // normalization
422 for (Int_t ipart=0; ipart<nPart; ipart++)
423 {
424 fGaussTPC[ipart]->SetParameter(0, 1./fGaussTPC[ipart]->Integral(-20,20));
425 }
426 // calculate correction factors
427 Float_t PIDCorrFactTPC[nPart];
428 for (Int_t ipart=0; ipart<nPart; ipart++)
429 {
430 PIDCorrFactTPC[ipart] = fGaussTPC[ipart]->Integral(-3*TMath::Sqrt(2), 3*TMath::Sqrt(2));
431 }
432 cout << "\n\n\n--- TPC ONLY CORRECTION FACTORS: ---\n"
433 << "Pions: " << PIDCorrFactTPC[0] << endl
434 << "Kaons: " << PIDCorrFactTPC[1] << endl
435 << "Protons: " << PIDCorrFactTPC[2] << endl << endl << endl;
436 // for the TPC+TOF range, get the correction factor from integrating the TPC+TOF nsigma distribution
437 TFile * fTPCTOFin = TFile::Open("TPCTOFnsigma.root");
438 TH1F * hTPCTOFnsigma[nPart];
439 Float_t PIDCorrFactTPCTOF[nPart];
440 for (Int_t ipart=0; ipart<nPart; ipart++)
441 {
442 hTPCTOFnsigma[ipart] = (TH1F*)((TH1F*)((TCanvas*)fTPCTOFin->Get("cTPCTOFnsigmaPredicted"))->FindObject(Form("hTPCTOFnsigma%s",Particle[ipart].Data())))->Clone();
443 PIDCorrFactTPCTOF[ipart] = hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0), hTPCTOFnsigma[ipart]->FindBin(3));
444 PIDCorrFactTPCTOF[ipart] /= hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0), hTPCTOFnsigma[ipart]->FindBin(20));
445 }
446 cout << "\n\n\n-- TPC+TOF CORRECTION FACTORS: ---\n"
447 << "Pions: " << PIDCorrFactTPCTOF[0] << endl
448 << "Kaons: " << PIDCorrFactTPCTOF[1] << endl
449 << "Protons: " << PIDCorrFactTPCTOF[2] << endl << endl << endl;
450
451 // calculate seperate correction factors for MC b/c the parameters will be different...
452 TF1 * fGaussTPCMC[nPart];
453 for (Int_t ipart=0; ipart<nPart; ipart++)
454 {
455 fGaussTPCMC[ipart] = new TF1(Form("fGaussTPCMC%s",Particle[ipart].Data()),"gaus",-25,25);
456 fGaussTPCMC[ipart]->SetParameter(0,1); // set a default height for now - will be normalized later
457 }
458 // pions
459 fGaussTPCMC[0]->SetParameter(1,-0.03);
460 fGaussTPCMC[0]->SetParameter(2,1.35);
461 // kaons
462 fGaussTPCMC[1]->SetParameter(1,0.14);
463 fGaussTPCMC[1]->SetParameter(2,1.27);
464 // protons
465 fGaussTPCMC[2]->SetParameter(1,-0.17);
466 fGaussTPCMC[2]->SetParameter(2,1.03);
467 // normalization
468 for (Int_t ipart=0; ipart<nPart; ipart++)
469 {
470 fGaussTPCMC[ipart]->SetParameter(0, 1./fGaussTPCMC[ipart]->Integral(-20,20));
471 }
472 // calculate correction factors
473 Float_t PIDCorrFactTPCMC[nPart];
474 for (Int_t ipart=0; ipart<nPart; ipart++)
475 {
476 PIDCorrFactTPCMC[ipart] = fGaussTPCMC[ipart]->Integral(-3*TMath::Sqrt(2), 3*TMath::Sqrt(2));
477 }
478 cout << "\n\n\n--- TPC ONLY CORRECTION FACTORS (MC): ---\n"
479 << "Pions: " << PIDCorrFactTPCMC[0] << endl
480 << "Kaons: " << PIDCorrFactTPCMC[1] << endl
481 << "Protons: " << PIDCorrFactTPCMC[2] << endl << endl << endl;
482 // for the TPC+TOF range, get the correction factor from integrating the TPC+TOF nsigma distribution
483 TFile * fTPCTOFin = TFile::Open("TPCTOFnsigma.root");
484 TH1F * hTPCTOFnsigmaMC[nPart];
485 Float_t PIDCorrFactTPCTOFMC[nPart];
486 for (Int_t ipart=0; ipart<nPart; ipart++)
487 {
488 hTPCTOFnsigmaMC[ipart] = (TH1F*)((TH1F*)((TCanvas*)fTPCTOFin->Get("cTPCTOFnsigmaPredictedMC"))->FindObject(Form("hTPCTOFnsigmaMC%s",Particle[ipart].Data())))->Clone();
489 PIDCorrFactTPCTOFMC[ipart] = hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0), hTPCTOFnsigmaMC[ipart]->FindBin(3));
490 PIDCorrFactTPCTOFMC[ipart] /= hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0), hTPCTOFnsigmaMC[ipart]->FindBin(20));
491 }
492 cout << "\n\n\n-- TPC+TOF CORRECTION FACTORS (MC): ---\n"
493 << "Pions: " << PIDCorrFactTPCTOFMC[0] << endl
494 << "Kaons: " << PIDCorrFactTPCTOFMC[1] << endl
495 << "Protons: " << PIDCorrFactTPCTOFMC[2] << endl << endl << endl;
496
497
498 // correct spectra, taking into account the range where the TOF is valid
499 for (Int_t icharge=0; icharge<nCharge; icharge++)
500 {
501 for (Int_t ipart=0; ipart<nPart; ipart++)
502 {
503 Int_t index = ipart + 3*icharge;
504 for (Int_t ibin=1; ibin<Spectra[index]->GetNbinsX(); ibin++)
505 {
506 if (ibin < Spectra[index]->FindBin(0.6))
507 Spectra[index]->SetBinContent(ibin,Spectra[index]->GetBinContent(ibin) / PIDCorrFactTPC[ipart]);
508 else
509 Spectra[index]->SetBinContent(ibin,Spectra[index]->GetBinContent(ibin) / PIDCorrFactTPCTOF[ipart]);
510 }
511 } // end for loop over ipart
512 } // end for loop over icharge
513 // save this version of the spectra for use later
514 TH1F * SpectraStep2[nCharge*nPart];
515 for (Int_t icharge=0; icharge<nCharge; icharge++)
516 {
517 for (Int_t ipart=0; ipart<nPart; ipart++)
518 {
519 Int_t index = ipart + 3*icharge;
520 SpectraStep2[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
521 SpectraStep2[index]->SetTitle("Spectra After PID Correction");
522 }
523 }
524 // also do this for MCSpectra
525 for (Int_t icharge=0; icharge<nCharge; icharge++)
526 {
527 for (Int_t ipart=0; ipart<nPart; ipart++)
528 {
529 Int_t index = ipart + 3*icharge;
530 for (Int_t ibin=1; ibin<MCSpectra[index]->GetNbinsX(); ibin++)
531 {
532 if (ibin < MCSpectra[index]->FindBin(0.6))
533 MCSpectra[index]->SetBinContent(ibin,MCSpectra[index]->GetBinContent(ibin) / PIDCorrFactTPCMC[ipart]);
534 else
535 MCSpectra[index]->SetBinContent(ibin,MCSpectra[index]->GetBinContent(ibin) / PIDCorrFactTPCTOFMC[ipart]);
536 }
537 } // end for loop over ipart
538 } // end for loop over icharge
539 // save this version of the spectra for use later
540 TH1F * MCSpectraStep2[nCharge*nPart];
541 for (Int_t icharge=0; icharge<nCharge; icharge++)
542 {
543 for (Int_t ipart=0; ipart<nPart; ipart++)
544 {
545 Int_t index = ipart + 3*icharge;
546 MCSpectraStep2[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
547 MCSpectraStep2[index]->SetTitle("MCSpectra After PID Correction");
548 }
549 }
550
551
552 //-------------------------------------------------------------------------
553 // Correct spectra for contamination from other Pions, Kaons, and Protons -
554 //-------------------------------------------------------------------------
555 // calculate the correction factors first
556 TH1F * ContCorrFact[nPart*nCharge];
557 for (Int_t icharge=0; icharge<nCharge; icharge++)
558 {
559 for (Int_t ipart=0; ipart<nPart; ipart++)
560 {
561 Int_t index = ipart + 3*icharge;
562 TString hname = Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
563 ContCorrFact[index] = new TH1F(*(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
564 ContCorrFact[index]->SetTitle("Contamination Correction Factor");
565 ContCorrFact[index]->SetMarkerStyle(Marker[index]);
566 ContCorrFact[index]->SetMarkerColor(Color[ipart]);
567 ContCorrFact[index]->SetLineColor(Color[ipart]);
568 hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
569 ContCorrFact[index]->Divide(ContCorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(),-1,-1,"B"); // binomial error
570 if (index == 0)
571 {
572 TCanvas * cContCorrFact = new TCanvas("cContCorrFact","cContCorrFact");
573 TLegend * lContCorrFact = new TLegend(.69,.69,.99,.99);
574 ContCorrFact[index]->DrawCopy("EP");
575 }
576 else ContCorrFact[index]->DrawCopy("EPsame");
577 lContCorrFact->AddEntry(ContCorrFact[index],Names[index].Data(),"p");
578 } // end for loop over ipart
579 } // end for loop over icharge
580 lContCorrFact->DrawClone();
581 // then use them to correct the spectra
582 for (Int_t icharge=0; icharge<nCharge; icharge++)
583 {
584 for(Int_t ipart=0; ipart<nPart; ipart++)
585 {
586 Int_t index = ipart + 3*icharge;
587 Spectra[index]->Multiply(ContCorrFact[index]);
588 }
589 }
590 // save this version of the spectra for use later
591 TH1F * SpectraStep3[nCharge*nPart];
592 for (Int_t icharge=0; icharge<nCharge; icharge++)
593 {
594 for (Int_t ipart=0; ipart<nPart; ipart++)
595 {
596 Int_t index = ipart + 3*icharge;
597 SpectraStep3[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
598 SpectraStep3[index]->SetTitle("Spectra After Contamination Correction");
599 }
600 }
601 // also do this for MCSpectra
602 for (Int_t icharge=0; icharge<nCharge; icharge++)
603 {
604 for(Int_t ipart=0; ipart<nPart; ipart++)
605 {
606 Int_t index = ipart + 3*icharge;
607 MCSpectra[index]->Multiply(ContCorrFact[index]);
608 }
609 }
610 // save this version of the spectra for use later
611 TH1F * MCSpectraStep3[nCharge*nPart];
612 for (Int_t icharge=0; icharge<nCharge; icharge++)
613 {
614 for (Int_t ipart=0; ipart<nPart; ipart++)
615 {
616 Int_t index = ipart + 3*icharge;
617 MCSpectraStep3[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
618 MCSpectraStep3[index]->SetTitle("MCSpectra After Contamination Correction");
619 }
620 }
621
622
623 //--------------------------------------------------
624 // Correct spectra for the Geant-Fluka discrepancy -
625 //--------------------------------------------------
626 GFCorrection(Spectra,tcuts_data,fout);
627 // save this version of the spectra for use later
628 TH1F * SpectraStep4[nCharge*nPart];
629 for (Int_t icharge=0; icharge<nCharge; icharge++)
630 {
631 for (Int_t ipart=0; ipart<nPart; ipart++)
632 {
633 Int_t index = ipart + 3*icharge;
634 SpectraStep4[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
635 SpectraStep4[index]->SetTitle("Spectra After Geant-Fluka Correction");
636 }
637 }
638 // NOTE! DO NOT USE GF FOR MCSPECTRA!
639
640
641 //--------------------------------------------------
642 // Correct spectra for the TOF matching efficiency -
643 //--------------------------------------------------
644 //Matching efficiency in data and Monte Carlo
645 TCanvas *cMatchingEff=new TCanvas("cMatchingEff","cMatchingEff",700,500);
646 TH1F *hMatcEffPos_data=(TH1F*)tcuts_data->GetHistoNMatchedPos();
647 hMatcEffPos_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedPos());
648 hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
649 TH1F *hMatcEffNeg_data=(TH1F*)tcuts_data->GetHistoNMatchedNeg();
650 hMatcEffNeg_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedNeg());
651 hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
652 hMatcEffNeg_data->SetLineColor(2);
653 /*
654 // for my (John Groh) report, plot the matching efficiency just for data
655 TCanvas * cMatchEffDataOnly = new TCanvas("cMatchEffDataOnly","cMatchEffDataOnly");
656 hMatcEffPos_data->SetTitle("TOF Matching Efficiency;p_{T} (GeV/c);");
657 hMatcEffPos_data->DrawCopy("lhist");
658 hMatcEffNeg_data->SetTitle("TOF Matching Efficiency;p_{T} (GeV/c);");
659 hMatcEffNeg_data->DrawCopy("lhistsame");
660 TLegend * lMatcEffDataOnly = new TLegend(.69,.69,.99,.99);
661 lMatcEffDataOnly->SetFillColor(0);
662 lMatcEffDataOnly->AddEntry(hMatcEffPos_data,"Positive Particles","l");
663 lMatcEffDataOnly->AddEntry(hMatcEffNeg_data,"Negative Particles","l");
664 lMatcEffDataOnly->DrawClone();
665 */
666
667 TH1F *hMatcEffPos_mc=(TH1F*)tcuts_mc->GetHistoNMatchedPos();
668 hMatcEffPos_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedPos());
669 hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
670 hMatcEffPos_mc->SetLineStyle(2);
671 TH1F *hMatcEffNeg_mc=(TH1F*)tcuts_mc->GetHistoNMatchedNeg();
672 hMatcEffNeg_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedNeg());
673 hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
674 hMatcEffNeg_mc->SetLineColor(2);
675 hMatcEffNeg_mc->SetLineStyle(2);
676 cMatchingEff->Divide(1,2);
677 cMatchingEff->cd(1);
678 gPad->SetGridy();
679 gPad->SetGridx();
680 hMatcEffPos_data->DrawClone("lhist");
681 hMatcEffNeg_data->DrawClone("lhistsame");
682 hMatcEffPos_mc->DrawClone("lhistsame");
683 hMatcEffNeg_mc->DrawClone("lhistsame");
684 gPad->BuildLegend()->SetFillColor(0);
685 hMatcEffPos_data->Divide(hMatcEffPos_mc);
686 hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
687 cMatchingEff->cd(2);
688 gPad->SetGridy();
689 gPad->SetGridx();
690 hMatcEffPos_data->DrawClone("lhist");
691 hMatcEffNeg_data->DrawClone("lhistsame");
692 TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",2.5,5);
693 hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
694 pol0MatchPos_data->DrawClone("same");
695 TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",2.5,5);
696 hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
697 pol0MatchNeg_data->SetLineColor(2);
698 pol0MatchNeg_data->DrawClone("same");
699 Float_t ScalingMatchingPos=pol0MatchPos_data->GetParameter(0);
700 Float_t ScalingMatchingNeg=pol0MatchNeg_data->GetParameter(0);
701 //Correction spectra for matching efficiency
702 for (Int_t ipart=0; ipart<nPart; ipart++)
703 {
704 for (Int_t ibin=1; ibin<=Spectra[ipart]->GetNbinsX(); ibin++)
705 {
706 Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
707 if (ptspectra < tcuts_data->GetPtTOFMatching()) continue;
708 Spectra[ipart]->SetBinContent(ibin, (Spectra[ipart]->GetBinContent(ibin) / ScalingMatchingPos));
709 Spectra[ipart+3]->SetBinContent(ibin, (Spectra[ipart+3]->GetBinContent(ibin) / ScalingMatchingNeg));
710 }
711 }
712 // save this version of the spectra for use later
713 TH1F * SpectraStep5[nCharge*nPart];
714 for (Int_t icharge=0; icharge<nCharge; icharge++)
715 {
716 for (Int_t ipart=0; ipart<nPart; ipart++)
717 {
718 Int_t index = ipart + 3*icharge;
719 SpectraStep5[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
720 SpectraStep5[index]->SetTitle("Spectra After Matching Efficiency Correction");
721 }
722 }
723 // NOTE! DO NOT USE THE MATCHING EFFICIENCY CORRECTION FOR MCSPECTRA!
724
725
726 //-----------------------------------------------------
727 // Secondaries correction - use MC, not DCA from data -
728 //------------------------------------------------------
729 // (correction factor for secondaries) = (true primaries) / (primaries)
730 TH1F * SecondaryCorrFact[nCharge*nPart];
731 for (Int_t icharge=0; icharge<nCharge; icharge++)
732 {
733 for (Int_t ipart=0; ipart<nPart; ipart++)
734 {
735 Int_t index = ipart + 3*icharge;
736 cout << "index = " << index << endl;
737 TString hname = Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
738 SecondaryCorrFact[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
739 hname = Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
740 SecondaryCorrFact[index]->Divide(SecondaryCorrFact[index], (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(), 1, 1, "B"); // binomial error!
741 SecondaryCorrFact[index]->SetName(Form("SecondaryCorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
742 SecondaryCorrFact[index]->SetTitle("Correction Factor for Secondaries");
743 SecondaryCorrFact[index]->SetStats(kFALSE);
744 SecondaryCorrFact[index]->SetMarkerStyle(Marker[index]);
745 SecondaryCorrFact[index]->SetMarkerColor(Color[ipart]);
746 SecondaryCorrFact[index]->SetLineColor(Color[ipart]);
747 if (index == 0)
748 {
749 TCanvas * cSecCorrFactMC = new TCanvas("cSecCorrFactMC","cSecCorrFactMC");
750 TLegend * lSecCorrFactMC = new TLegend(.69,.69,.99,.99);
751 lSecCorrFactMC->SetFillColor(0);
752 SecondaryCorrFact[index]->DrawCopy("EP");
753 }
754 else SecondaryCorrFact[index]->DrawCopy("EPsame");
755 lSecCorrFactMC->AddEntry(SecondaryCorrFact[index],Names[index],"ep");
756 } // end for loop over ipart
757 } // end for loop over icharge
758 lSecCorrFactMC->DrawClone();
759
760 // correct Spectra with this
761 for (Int_t icharge=0; icharge<nCharge; icharge++)
762 {
763 for (Int_t ipart=0; ipart<nPart; ipart++)
764 {
765 Int_t index = ipart + 3*icharge;
766 Spectra[index]->Multiply(SecondaryCorrFact[index]);
767 }
768 }
769
770 // also correct MCSpectra with this
771 for (Int_t icharge=0; icharge<nCharge; icharge++)
772 {
773 for (Int_t ipart=0; ipart<nPart; ipart++)
774 {
775 Int_t index = ipart + 3*icharge;
776 MCSpectra[index]->Multiply(SecondaryCorrFact[index]);
777 }
778 }
779 // save this version of the spectra for use later
780 TH1F * SpectraStep6[nCharge*nPart];
781 for (Int_t icharge=0; icharge<nCharge; icharge++)
782 {
783 for (Int_t ipart=0; ipart<nPart; ipart++)
784 {
785 Int_t index = ipart + 3*icharge;
786 SpectraStep6[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
787 SpectraStep6[index]->SetTitle("Spectra After Secondaries Correction");
788 }
789 }
790 // also save this version of MCSpectra
791 TH1F * MCSpectraStep6[nCharge*nPart];
792 for (Int_t icharge=0; icharge<nCharge; icharge++)
793 {
794 for (Int_t ipart=0; ipart<nPart; ipart++)
795 {
796 Int_t index = ipart + 3*icharge;
797 MCSpectraStep6[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
798 MCSpectraStep6[index]->SetTitle("MCSpectra After Secondaries Correction");
799 }
800 }
801
802 //-------------------------
803 // draw the final spectra -
804 //-------------------------
805 TCanvas * cSpectra = new TCanvas("cSpectra","cSpectra",100,100,700,500);
806 gPad->SetGridy();
807 gPad->SetLogy();
808 TLegend * lSpectra = new TLegend(.69,.69,.99,.99);
809 for (Int_t icharge=0; icharge<nCharge; icharge++)
810 {
811 for (Int_t ipart=0; ipart<nPart; ipart++)
812 {
813 Int_t index = ipart + 3*icharge;
814 Spectra[index]->SetTitle("Final Spectra");
815 if (index == 0) Spectra[index]->DrawClone("EP");
816 else Spectra[index]->DrawClone("EPsame");
817 lSpectra->AddEntry(Spectra[index],Names[index].Data(),"p");
818 }
819 }
820 lSpectra->DrawClone();
821
822 //-----------------------------------------------
823 // compare MCSpectra with hHistPtGenTruePrimary -
824 //-----------------------------------------------
825 TCanvas * cCompareMCRecWithMCTrue = new TCanvas("cCompareMCRecWithMCTrue","cCompareMCRecWithMCTrue",250,175,700,500);
826 cCompareMCRecWithMCTrue->Divide(3,2);
827 TH1F * hMCTruth[nPart*nCharge];
828 for (Int_t icharge=0; icharge<nCharge; icharge++)
829 {
830 for (Int_t ipart=0; ipart<nPart; ipart++)
831 {
832 Int_t index = ipart + 3*icharge;
833 hMCTruth[index] = new TH1F(*(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),1,1))->Clone());
834 hMCTruth[index]->Scale(1./ecuts_mc->NumberOfEvents(),"width"); // normalize it first!
835 hMCTruth[index]->SetMarkerStyle(Marker[index]);
836 hMCTruth[index]->SetMarkerColor(Color[ipart]);
837 hMCTruth[index]->SetLineColor(Color[ipart]);
838 }
839 }
840 // set bins outside range to zero
841 for (Int_t icharge=0; icharge<nCharge; icharge++)
842 {
843 for (Int_t ipart=0; ipart<nPart; ipart++)
844 {
845 Int_t index = ipart + 3*icharge;
846 for (Int_t ibin=0; ibin<hMCTruth[index]->GetNbinsX(); ibin++)
847 {
848 if (hMCTruth[index]->GetBinCenter(ibin) < Range[ipart])
849 {
850 hMCTruth[index]->SetBinContent(ibin, 0);
851 hMCTruth[index]->SetBinError(ibin, 0);
852 }
853 }
854 }
855 }
856
857
858 TH1F * hRatioMCRecTrue[nPart*nCharge];
859 for (Int_t ipart=0; ipart<nPart; ipart++)
860 {
861 for (Int_t icharge=0; icharge<nCharge; icharge++)
862 {
863 Int_t index = ipart + 3*icharge;
864 hRatioMCRecTrue[index] = (TH1F*)MCSpectra[index]->Clone();
865 hRatioMCRecTrue[index]->Divide(hMCTruth[index]);
866 hRatioMCRecTrue[index]->SetTitle("(Corrected MC Spectra) / (MC Truth Spectra);p_{T} (GeV/c);");
867 if (icharge == 1) hRatioMCRecTrue[index]->SetLineStyle(2);
868 }
869 }
870 for (Int_t ipart=0; ipart<nPart; ipart++)
871 {
872 cCompareMCRecWithMCTrue->cd(ipart+1);
873 gPad->SetLogy();
874 MCSpectra[ipart]->SetMarkerColor(kGreen+3);
875 MCSpectra[ipart+nPart]->SetMarkerColor(kGreen+2);
876 MCSpectra[ipart]->SetMarkerSize(0.75);
877 MCSpectra[ipart+nPart]->SetMarkerSize(0.75);
878 MCSpectra[ipart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
879 MCSpectra[ipart+nPart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
880 hMCTruth[ipart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
881 hMCTruth[ipart+nPart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
882 hMCTruth[ipart]->SetMarkerSize(0.75);
883 hMCTruth[ipart]->SetMarkerSize(0.75);
884 hMCTruth[ipart]->DrawCopy("EP");
885 hMCTruth[ipart+nPart]->DrawCopy("EPsame");
886 MCSpectra[ipart]->DrawCopy("EPsame");
887 MCSpectra[ipart+nPart]->DrawCopy("EPsame");
888 TLegend * lCompare = new TLegend(.69,.59,.99,.91);
889 lCompare->SetFillColor(0);
890 lCompare->AddEntry(hMCTruth[ipart],Form("%s, Truth",Names[ipart].Data()),"p");
891 lCompare->AddEntry(hMCTruth[ipart+nPart],Form("%s, Truth",Names[ipart+nPart].Data()),"p");
892 lCompare->AddEntry(MCSpectra[ipart],Form("%s, Corrected",Names[ipart].Data()),"p");
893 lCompare->AddEntry(MCSpectra[ipart+nPart],Form("%s, Corrected",Names[ipart+nPart].Data()),"p");
894 lCompare->DrawClone();
895
896 cCompareMCRecWithMCTrue->cd(ipart+1+nPart);
897 hRatioMCRecTrue[ipart]->GetYaxis()->SetRangeUser(.97,1.03);
898 hRatioMCRecTrue[ipart]->DrawCopy("hist][");
899 hRatioMCRecTrue[ipart+nPart]->DrawCopy("hist][same");
900 TLegend * lCompareRatio = new TLegend(.79,.79,.99,.99);
901 lCompareRatio->SetFillColor(0);
902 lCompareRatio->AddEntry(hRatioMCRecTrue[ipart],Names[ipart].Data(),"l");
903 lCompareRatio->AddEntry(hRatioMCRecTrue[ipart+nPart],Names[ipart+nPart].Data(),"l");
904 lCompareRatio->DrawClone();
905 }
906
907
908 //---------------------------------------------------------
909 // plot the evolution of the spectra afte each correction -
910 //---------------------------------------------------------
911 TCanvas * cEvolutionOfSpectra = new TCanvas("cEvolutionOfSpectra","cEvolutionOfSpectra",300,200,700,500);
912 cEvolutionOfSpectra->Divide(4,2);
913 TLegend * lEvolutionOfSpectra = new TLegend(.69,.69,.99,.99);
914 for (Int_t icharge=0; icharge<nCharge; icharge++)
915 {
916 for (Int_t ipart=0; ipart<nPart; ipart++)
917 {
918 Int_t index = ipart + 3*icharge;
919 lEvolutionOfSpectra->AddEntry(Spectra[index],Names[index].Data(),"p");
920 }
921 }
922 cEvolutionOfSpectra->cd(1);
923 gPad->SetLogy();
924 for (Int_t icharge=0; icharge<nCharge; icharge++)
925 {
926 for (Int_t ipart=0; ipart<nPart; ipart++)
927 {
928 Int_t index = ipart + 3*icharge;
929 RawSpectra[index]->SetMarkerSize(.5);
930 if (index == 0) RawSpectra[index]->DrawClone("EP");
931 else RawSpectra[index]->DrawClone("EPsame");
932 }
933 }
934 lEvolutionOfSpectra->DrawClone();
935 cEvolutionOfSpectra->cd(2);
936 gPad->SetLogy();
937 for (Int_t icharge=0; icharge<nCharge; icharge++)
938 {
939 for (Int_t ipart=0; ipart<nPart; ipart++)
940 {
941 Int_t index = ipart + 3*icharge;
942 SpectraStep1[index]->SetMarkerSize(.5);
943 if (index == 0) SpectraStep1[index]->DrawClone("EP");
944 else SpectraStep1[index]->DrawClone("EPsame");
945 }
946 }
947 lEvolutionOfSpectra->DrawClone();
948 cEvolutionOfSpectra->cd(3);
949 gPad->SetLogy();
950 for (Int_t icharge=0; icharge<nCharge; icharge++)
951 {
952 for (Int_t ipart=0; ipart<nPart; ipart++)
953 {
954 Int_t index = ipart + 3*icharge;
955 SpectraStep2[index]->SetMarkerSize(.5);
956 if (index == 0) SpectraStep2[index]->DrawClone("EP");
957 else SpectraStep2[index]->DrawClone("EPsame");
958 }
959 }
960 lEvolutionOfSpectra->DrawClone();
961 cEvolutionOfSpectra->cd(4);
962 gPad->SetLogy();
963 for (Int_t icharge=0; icharge<nCharge; icharge++)
964 {
965 for (Int_t ipart=0; ipart<nPart; ipart++)
966 {
967 Int_t index = ipart + 3*icharge;
968 SpectraStep3[index]->SetMarkerSize(.5);
969 if (index == 0) SpectraStep3[index]->DrawClone("EP");
970 else SpectraStep3[index]->DrawClone("EPsame");
971 }
972 }
973 lEvolutionOfSpectra->DrawClone();
974 cEvolutionOfSpectra->cd(5);
975 gPad->SetLogy();
976 for (Int_t icharge=0; icharge<nCharge; icharge++)
977 {
978 for (Int_t ipart=0; ipart<nPart; ipart++)
979 {
980 Int_t index = ipart + 3*icharge;
981 SpectraStep4[index]->SetMarkerSize(.5);
982 if (index == 0) SpectraStep4[index]->DrawClone("EP");
983 else SpectraStep4[index]->DrawClone("EPsame");
984 }
985 }
986 lEvolutionOfSpectra->DrawClone();
987 cEvolutionOfSpectra->cd(6);
988 gPad->SetLogy();
989 for (Int_t icharge=0; icharge<nCharge; icharge++)
990 {
991 for (Int_t ipart=0; ipart<nPart; ipart++)
992 {
993 Int_t index = ipart + 3*icharge;
994 SpectraStep5[index]->SetMarkerSize(.5);
995 if (index == 0) SpectraStep5[index]->DrawClone("EP");
996 else SpectraStep5[index]->DrawClone("EPsame");
997 }
998 }
999 lEvolutionOfSpectra->DrawClone();
1000 cEvolutionOfSpectra->cd(7);
1001 gPad->SetLogy();
1002 for (Int_t icharge=0; icharge<nCharge; icharge++)
1003 {
1004 for (Int_t ipart=0; ipart<nPart; ipart++)
1005 {
1006 Int_t index = ipart + 3*icharge;
1007 SpectraStep6[index]->SetMarkerSize(.5);
1008 if (index == 0) SpectraStep6[index]->DrawClone("EP");
1009 else SpectraStep6[index]->DrawClone("EPsame");
1010 }
1011 }
1012 lEvolutionOfSpectra->DrawClone();
1013 cEvolutionOfSpectra->cd(8);
1014 gPad->SetLogy();
1015 for (Int_t icharge=0; icharge<nCharge; icharge++)
1016 {
1017 for (Int_t ipart=0; ipart<nPart; ipart++)
1018 {
1019 Int_t index = ipart + 3*icharge;
1020 Spectra[index]->SetMarkerSize(.5);
1021 if (index == 0) Spectra[index]->DrawClone("EP");
1022 else Spectra[index]->DrawClone("EPsame");
1023 }
1024 }
1025 lEvolutionOfSpectra->DrawClone();
1026
1027
1028 //-----------------------------------------
1029 // Compare spectra with published spectra -
1030 //-----------------------------------------
1031 // also add another canvas that just plots the ratios b/c I can't use the spectra in my presentation (hasn't been published yet)
1032 TCanvas * cRatioCombined = new TCanvas("cRatioCombined","cRatioCombined");
1033 cRatioCombined->Divide(3,1);
1034 TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",150,125,700,500);
1035 CratioComb->Divide(3,2);
1036 TString nameComb[6] = {Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
1037 Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
1038 TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120730.root");
1039 TLegend * lComp[nPart];
1040 TLegend * lRatio[nPart];
1041 for(Int_t icharge=0; icharge<nCharge; icharge++)
1042 {
1043 for(Int_t ipart=0; ipart<nPart; ipart++)
1044 {
1045 Int_t index = ipart + 3*icharge;
1046 TH1F * htmp = (TH1F*)((TH1F*)Spectra[index])->Clone("");
1047 htmp->SetMarkerSize(1);
1048 htmp->GetXaxis()->SetRangeUser(0,5);
1049 TH1F * hcomb = (TH1F*)fComb->Get(nameComb[index].Data())->Clone();
1050 CratioComb->cd(ipart+1);
1051 gPad->SetGridy();
1052 gPad->SetGridx();
1053 htmp->SetTitle("Comparison with Combined Analysis;p_{T} (GeV/c);");
1054 if(icharge == 0) htmp->DrawCopy();
1055 else htmp->DrawCopy("same");
1056 hcomb->SetMarkerStyle(Marker[index]);
1057 hcomb->DrawCopy("same");
1058 if (icharge == 0) {lComp[ipart] = new TLegend(.49,.49,.89,.89); lComp[ipart]->SetFillColor(0);}
1059 lComp[ipart]->AddEntry(htmp,Form("%s, This Analysis",Names[index].Data()),"p");
1060 lComp[ipart]->AddEntry(hcomb,Form("%s, Combined Analysis",Names[index].Data()),"lep");
1061 if (icharge == 1) lComp[ipart]->DrawClone();
1062
1063 CratioComb->cd(ipart+nPart+1);
1064 gPad->SetGridy();
1065 gPad->SetGridx();
1066 htmp->Divide(hcomb);
1067 htmp->SetTitle("(This Analysis) / (Combined Analysis);p_{T} (GeV/c);");
1068 htmp->GetYaxis()->SetRangeUser(0.8, 1.2);
1069 if (icharge == 0)
1070 {
1071 CratioComb->cd(ipart+nPart+1);
1072 htmp->DrawClone("hist][");
1073 cRatioCombined->cd(ipart+1);
1074 htmp->DrawClone("hist][");
1075 lRatio[ipart] = new TLegend(.69,.69,.89,.89);
1076 lRatio[ipart]->SetFillColor(0);
1077 }
1078 else
1079 {
1080 htmp->SetLineStyle(2);
1081 CratioComb->cd(ipart+nPart+1);
1082 htmp->DrawClone("histsame][");
1083 cRatioCombined->cd(ipart+1);
1084 htmp->DrawClone("histsame][");
1085 }
1086 lRatio[ipart]->AddEntry(htmp,Names[index].Data(),"l");
1087 if (icharge == 1)
1088 {
1089 CratioComb->cd(ipart+nPart+1);
1090 lRatio[ipart]->DrawClone();
1091 cRatioCombined->cd(ipart+1);
1092 lRatio[ipart]->DrawClone();
1093 }
1094 }
1095 }
1096 //comparison with charged hadron
1097 Printf("\n\n-> ChargedHadron comparison");
1098 TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
1099 TH1F* hMatchCorrectionAllCh=(TH1F*)hMatcEffPos_data->Clone("hMatchCorrectionAllCh");
1100 hMatchCorrectionAllCh->Add(hMatcEffNeg_data); //correction for Matching efficiency!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1101 hMatchCorrectionAllCh->Scale(0.5);
1102 for(Int_t ibin=1;ibin<hChHad_data->GetNbinsX();ibin++){
1103 Float_t ptch=hChHad_data->GetBinCenter(ibin);
1104 if(ptch<tcuts_data->GetPtTOFMatching())continue;
1105 hChHad_data->SetBinContent(ibin,2*(hChHad_data->GetBinContent(ibin)/(ScalingMatchingPos+ScalingMatchingNeg)));
1106 }
1107 //fraction of sec in MC
1108 TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
1109 TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
1110 for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
1111 Double_t en_data=hChHad_data->GetBinContent(ibin);
1112 Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
1113 Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
1114 if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.2/en_mc));
1115 //Printf("Before: %.0f After: %.0f fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
1116 }
1117 hPrimRec_mc->Divide(hAllRec_mc);
1118 //efficiency for primaries
1119 TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
1120 hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->Clone());
1121 TCanvas *cAllChFactors=new TCanvas("cAllChFactors","cAllChFactors",700,500);
1122 cAllChFactors->Divide(1,2);
1123 cAllChFactors->cd(1);
1124 gPad->SetGridy();
1125 gPad->SetGridx();
1126 hPrimRec_mc->SetTitle("Prim/All, charged hadron pure MC");
1127 hPrimRec_mc->DrawClone("lhist");
1128 gPad->BuildLegend()->SetFillColor(0);
1129 cAllChFactors->cd(2);
1130 gPad->SetGridy();
1131 gPad->SetGridx();
1132 hEff_mc->SetTitle("Efficiency for Primaries, charged hadron pure MC");
1133 hEff_mc->DrawClone("lhist");
1134 gPad->BuildLegend()->SetFillColor(0);
1135 //Printf("--------%f ",((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
1136 hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
1137 hChHad_data->Divide(hEff_mc);//Efficiency
1138 hChHad_data->Scale(1./(TMath::Abs(tcuts_data->GetEtaMin())+TMath::Abs(tcuts_data->GetEtaMax())));
1139 hChHad_data->SetTitle("All Ch from AOD");
1140 hChHad_data->SetName("AllCh");
1141 TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
1142 cAllCh->Divide(1,2);
1143 cAllCh->cd(1);
1144 gPad->SetGridy();
1145 gPad->SetGridx();
1146 hChHad_data->GetYaxis()->SetRangeUser(0,2500);
1147 hChHad_data->DrawClone();
1148 TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
1149 TH1F *hCh=fCh->Get("hpt_c0_5");
1150 hCh->SetTitle("All Ch from Jacek");
1151 hCh->SetMarkerColor(2);
1152 hCh->SetLineColor(2);
1153 //invariant yield
1154 for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
1155 hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
1156 hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
1157 }
1158 hCh->DrawClone("same");
1159 gPad->BuildLegend()->SetFillColor(0);
1160 TH1F *gRatio=AliPWGHistoTools::DivideHistosDifferentBins(hChHad_data,hCh);
1161 gRatio->SetMaximum(1.3);
1162 gRatio->SetMinimum(.7);
1163 cAllCh->cd(2);
1164 gPad->SetGridy();
1165 gPad->SetGridx();
1166 gRatio->Print("all");
1167 gRatio->SetBinContent(48,1.02);
1168 gRatio->SetBinContent(49,1.022);
1169 gRatio->SetBinContent(50,1.021);
1170 gRatio->SetBinContent(51,1.026);
1171 gRatio->SetBinContent(52,1.027);
1172 gRatio->SetBinError(48,0.056);
1173 gRatio->SetBinError(49,0.056);
1174 gRatio->SetBinError(50,0.056);
1175 gRatio->SetBinError(51,0.056);
1176 gRatio->SetBinError(52,0.056);
1177 gRatio->GetYaxis()->SetRangeUser(-1.5,1.5);
1178 gRatio->DrawClone("");
1179
1180 //fitting
1181 TCanvas *cFitChargHad=new TCanvas("cFitChargHad","cFitChargHad",700,500);
1182 gPad->SetGridy();
1183 gPad->SetGridx();
1184 hChHad_data->DrawClone();
1185 //Fitting the sum of all particles
1186 AliPWGFunc * fm = new AliPWGFunc;
1187 fm->SetVarType(AliPWGFunc::kdNdpt);
1188 Float_t fitmin = 0.3;
1189 Float_t fitmax = 3;
1190 TF1 * func = 0;
1191 Int_t normPar = 0;
1192 func = fm->GetBGBW(0.13,0.6,0.3, 1, 1e5);
1193 func->SetParLimits(1, 0.1, 0.99);
1194 func->SetParLimits(2, 0.01, 1);
1195 func->SetParLimits(3, 0.01, 2);
1196 TH1F * hToFit = hChHad_data;
1197 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1198 // // if(!AliPWGHistoTools::Fit(hToFit,func,fitmin,fitmax)) {
1199 // // cout << " FIT ERROR " << endl;
1200 // // return;
1201 // // }
1202 Double_t yieldTools, yieldETools;
1203 Double_t partialYields[3],partialYieldsErrors[3];
1204 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1205 func->DrawClone("same");
1206 Printf("TOTAL YIELD (AOD Charged Hadron) : %f +- %f",yieldTools,yieldETools);
1207 //Fit All Charged
1208 hToFit = hCh;
1209 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1210 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1211 func->SetLineColor(2);
1212 hCh->DrawClone("same");
1213 func->DrawClone("same");
1214 gPad->BuildLegend()->SetFillColor(0);
1215 Printf("TOTAL YIELD (JACEK): %f +- %f",yieldTools,yieldETools);
1216 //sumID vs AllCh
1217 //Convert spectra to dNdeta and sum
1218 TH1F * hsum = (TH1F*) Spectra[0]->Clone();
1219 hsum->Reset("all");
1220 hsum->SetMarkerSize(1);
1221 Double_t epsilon = 0.0001;
1222 for(Int_t icharge = 0; icharge < 2; icharge++){
1223 for(Int_t ipart = 0; ipart < 3; ipart++){
1224 Int_t index=ipart+3*icharge;
1225 TH1F *htmp =(TH1F*)Spectra[index]->Clone("htmp");
1226 Int_t nbin = htmp->GetNbinsX();
1227 for(Int_t ibin = 1; ibin <= nbin; ibin++){
1228 Double_t pt = htmp->GetBinCenter(ibin);
1229 Double_t eta=0.8;//////////////////eta range///////////////////////////////////////
1230 Double_t jacobian =eta2y(pt,mass[ipart],eta)/eta;
1231 //Printf("jacobian: %f pt:%f BinContent:%f mass:%f",jacobian,pt,htmp->GetBinContent(ibin),mass[ipart]);
1232 htmp->SetBinContent(ibin,htmp->GetBinContent(ibin)*jacobian);
1233 htmp->SetBinError(ibin,htmp->GetBinError(ibin)*jacobian);
1234 Int_t ibinSum = hsum->FindBin(pt);
1235 if ( htmp->GetBinContent(ibin) > 0 &&
1236 (TMath::Abs(htmp->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
1237 TMath::Abs(htmp->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1238 ) {
1239 cout << "DISCREPANCY IN BIN RANGES" << endl;
1240 cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
1241 << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
1242 << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
1243 cout << "" << endl;
1244 }
1245 }
1246 hsum->Add(htmp);
1247 }
1248 }
1249 hsum->SetTitle("Sum ID from AOD");
1250 TCanvas *cChargHadComp=new TCanvas("cChargHadComp","cChargHadComp",200,150,700,500);
1251 cChargHadComp->Divide(1,2);
1252 cChargHadComp->cd(1);
1253 gPad->SetGridy();
1254 gPad->SetGridx();
1255 hsum->SetTitle("Charged Hadron Comparison;p_{T} (GeV/c);");
1256 hsum->DrawClone();
1257 hToFit = hsum;
1258 hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1259 AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1260 func->SetLineColor(2);
1261 Printf("TOTAL YIELD (Pi+K+p): %f +- %f",yieldTools,yieldETools);
1262 hChHad_data->SetMarkerColor(2);
1263 hChHad_data->SetLineColor(2);
1264 hChHad_data->DrawClone("same");
1265 TLegend * lChargHadComp = new TLegend(.69,.69,.89,.89);
1266 lChargHadComp->AddEntry(hsum,"Sum ID from AOD","p");
1267 lChargHadComp->AddEntry(hChHad_data,"All Ch from AOD","p");
1268 lChargHadComp->DrawClone();
1269 cChargHadComp->cd(2);
1270 gPad->SetGridy();
1271 gPad->SetGridx();
1272 TH1F *hRatio=AliPWGHistoTools::DivideHistosDifferentBins(hsum,hChHad_data);
1273 hRatio->SetTitle("(Sum ID from AOD) / (All Ch from AOD);p_{T} (GeV/c);");
1274 hRatio->SetMaximum(1.2);
1275 hRatio->SetMinimum(.8);
1276 hRatio->DrawClone("");
1277
1278
1279 //----------------------------------------------------------
1280 // Compare each correction step with the published spectra -
1281 //----------------------------------------------------------
1282 TCanvas * cCompareToPublishedAtEachStep = new TCanvas("cCompareToPublishedAtEachStep","cCompareToPublishedAtEachStep",350,225,700,500);
1283 cCompareToPublishedAtEachStep->Divide(4,2);
1284 TLegend * lCompareToPublishedAtEachStep = new TLegend(.69,.69,.99,.99);
1285 for (Int_t icharge=0; icharge<nCharge; icharge++)
1286 {
1287 for (Int_t ipart=0; ipart<nPart; ipart++)
1288 {
1289 Int_t index = ipart + 3*icharge;
1290 lCompareToPublishedAtEachStep->AddEntry(Spectra[index],Names[index].Data(),"p");
1291 }
1292 }
1293 TH1F * hComb[nPart*nCharge];
1294 for (Int_t icharge=0; icharge<nCharge; icharge++)
1295 {
1296 for (Int_t ipart=0; ipart<nPart; ipart++)
1297 {
1298 Int_t index = ipart + 3*icharge;
1299 fComb->cd();
1300 hComb[index] = (TH1F*)fComb->Get(nameComb[index].Data())->Clone();
1301 }
1302 }
1303 cCompareToPublishedAtEachStep->cd(1);
1304 for (Int_t icharge=0; icharge<nCharge; icharge++)
1305 {
1306 for (Int_t ipart=0; ipart<nPart; ipart++)
1307 {
1308 Int_t index = ipart + 3*icharge;
1309 RawSpectra[index]->Divide(hComb[index]);
1310 if (icharge == 1) RawSpectra[index]->SetLineStyle(2);
1311 if (index == 0)
1312 {
1313 RawSpectra[index]->GetYaxis()->SetRangeUser(0,1);
1314 RawSpectra[index]->DrawCopy("hist][");
1315 }
1316 else RawSpectra[index]->DrawCopy("hist][same");
1317 }
1318 }
1319 lCompareToPublishedAtEachStep->DrawClone();
1320 cCompareToPublishedAtEachStep->cd(2);
1321 for (Int_t icharge=0; icharge<nCharge; icharge++)
1322 {
1323 for (Int_t ipart=0; ipart<nPart; ipart++)
1324 {
1325 Int_t index = ipart + 3*icharge;
1326 SpectraStep1[index]->Divide(hComb[index]);
1327 if (icharge == 1) SpectraStep1[index]->SetLineStyle(2);
1328 if (index == 0)
1329 {
1330 SpectraStep1[index]->GetYaxis()->SetRangeUser(.5,1.5);
1331 SpectraStep1[index]->DrawCopy("hist][");
1332 }
1333 else SpectraStep1[index]->DrawCopy("hist][same");
1334 }
1335 }
1336 lCompareToPublishedAtEachStep->DrawClone();
1337 cCompareToPublishedAtEachStep->cd(3);
1338 for (Int_t icharge=0; icharge<nCharge; icharge++)
1339 {
1340 for (Int_t ipart=0; ipart<nPart; ipart++)
1341 {
1342 Int_t index = ipart + 3*icharge;
1343 SpectraStep2[index]->Divide(hComb[index]);
1344 if (icharge == 1) SpectraStep2[index]->SetLineStyle(2);
1345 if (index == 0)
1346 {
1347 SpectraStep2[index]->GetYaxis()->SetRangeUser(.5,1.5);
1348 SpectraStep2[index]->DrawCopy("hist][");
1349 }
1350 else SpectraStep2[index]->DrawCopy("hist][same");
1351 }
1352 }
1353 lCompareToPublishedAtEachStep->DrawClone();
1354 cCompareToPublishedAtEachStep->cd(4);
1355 for (Int_t icharge=0; icharge<nCharge; icharge++)
1356 {
1357 for (Int_t ipart=0; ipart<nPart; ipart++)
1358 {
1359 Int_t index = ipart + 3*icharge;
1360 SpectraStep3[index]->Divide(hComb[index]);
1361 if (icharge == 1) SpectraStep3[index]->SetLineStyle(2);
1362 if (index == 0)
1363 {
1364 SpectraStep3[index]->GetYaxis()->SetRangeUser(.5,1.5);
1365 SpectraStep3[index]->DrawCopy("hist][");
1366 }
1367 else SpectraStep3[index]->DrawCopy("hist][same");
1368 }
1369 }
1370 lCompareToPublishedAtEachStep->DrawClone();
1371 cCompareToPublishedAtEachStep->cd(5);
1372 for (Int_t icharge=0; icharge<nCharge; icharge++)
1373 {
1374 for (Int_t ipart=0; ipart<nPart; ipart++)
1375 {
1376 Int_t index = ipart + 3*icharge;
1377 SpectraStep4[index]->Divide(hComb[index]);
1378 if (icharge == 1) SpectraStep4[index]->SetLineStyle(2);
1379 if (index == 0)
1380 {
1381 SpectraStep4[index]->GetYaxis()->SetRangeUser(.5,1.5);
1382 SpectraStep4[index]->DrawCopy("hist][");
1383 }
1384 else SpectraStep4[index]->DrawCopy("hist][same");
1385 }
1386 }
1387 lCompareToPublishedAtEachStep->DrawClone();
1388 cCompareToPublishedAtEachStep->cd(6);
1389 for (Int_t icharge=0; icharge<nCharge; icharge++)
1390 {
1391 for (Int_t ipart=0; ipart<nPart; ipart++)
1392 {
1393 Int_t index = ipart + 3*icharge;
1394 SpectraStep5[index]->Divide(hComb[index]);
1395 if (icharge == 1) SpectraStep5[index]->SetLineStyle(2);
1396 if (index == 0)
1397 {
1398 SpectraStep5[index]->GetYaxis()->SetRangeUser(.5,1.5);
1399 SpectraStep5[index]->DrawCopy("hist][");
1400 }
1401 else SpectraStep5[index]->DrawCopy("hist][same");
1402 }
1403 }
1404 lCompareToPublishedAtEachStep->DrawClone();
1405 cCompareToPublishedAtEachStep->cd(7);
1406 for (Int_t icharge=0; icharge<nCharge; icharge++)
1407 {
1408 for (Int_t ipart=0; ipart<nPart; ipart++)
1409 {
1410 Int_t index = ipart + 3*icharge;
1411 SpectraStep6[index]->Divide(hComb[index]);
1412 if (icharge == 1) SpectraStep6[index]->SetLineStyle(2);
1413 if (index == 0)
1414 {
1415 SpectraStep6[index]->GetYaxis()->SetRangeUser(.5,1.5);
1416 SpectraStep6[index]->DrawCopy("hist][");
1417 }
1418 else SpectraStep6[index]->DrawCopy("hist][same");
1419 }
1420 }
1421 lCompareToPublishedAtEachStep->DrawClone();
1422 cCompareToPublishedAtEachStep->cd(8);
1423 for (Int_t icharge=0; icharge<nCharge; icharge++)
1424 {
1425 for (Int_t ipart=0; ipart<nPart; ipart++)
1426 {
1427 Int_t index = ipart + 3*icharge;
1428 Spectra[index]->Divide(hComb[index]);
1429 if (icharge == 1) Spectra[index]->SetLineStyle(2);
1430 if (index == 0)
1431 {
1432 Spectra[index]->GetYaxis()->SetRangeUser(.5,1.5);
1433 Spectra[index]->DrawCopy("hist][");
1434 }
1435 else Spectra[index]->DrawCopy("hist][same");
1436 }
1437 }
1438 lCompareToPublishedAtEachStep->DrawClone();
1439
1440
1441
1442 //---------------------------------------------------------
1443 // For MC, compare each correction step with the MC truth -
1444 //---------------------------------------------------------
1445 TCanvas * cCompareMCtoMCTruthAtEachStep = new TCanvas("cCompareMCtoMCTruthAtEachStep","cCompareMCtoMCTruthAtEachStep",350,225,700,500);
1446 cCompareMCtoMCTruthAtEachStep->Divide(4,2);
1447 TLegend * lCompareMCtoMCTruthAtEachStep = new TLegend(.69,.69,.99,.99);
1448 for (Int_t icharge=0; icharge<nCharge; icharge++)
1449 {
1450 for (Int_t ipart=0; ipart<nPart; ipart++)
1451 {
1452 Int_t index = ipart + 3*icharge;
1453 lCompareMCtoMCTruthAtEachStep->AddEntry(MCSpectra[index],Names[index].Data(),"p");
1454 }
1455 }
1456 cCompareMCtoMCTruthAtEachStep->cd(1);
1457 for (Int_t icharge=0; icharge<nCharge; icharge++)
1458 {
1459 for (Int_t ipart=0; ipart<nPart; ipart++)
1460 {
1461 Int_t index = ipart + 3*icharge;
1462 MCRawSpectra[index]->Divide(hMCTruth[index]);
1463 if (icharge == 1) RawSpectra[index]->SetLineStyle(2);
1464 if (index == 0)
1465 {
1466 MCRawSpectra[index]->GetYaxis()->SetRangeUser(0,1);
1467 MCRawSpectra[index]->DrawCopy("hist][");
1468 }
1469 else MCRawSpectra[index]->DrawCopy("hist][same");
1470 }
1471 }
1472 lCompareMCtoMCTruthAtEachStep->DrawClone();
1473 cCompareMCtoMCTruthAtEachStep->cd(2);
1474 for (Int_t icharge=0; icharge<nCharge; icharge++)
1475 {
1476 for (Int_t ipart=0; ipart<nPart; ipart++)
1477 {
1478 Int_t index = ipart + 3*icharge;
1479 MCSpectraStep1[index]->Divide(hMCTruth[index]);
1480 if (icharge == 1) MCSpectraStep1[index]->SetLineStyle(2);
1481 if (index == 0)
1482 {
1483 MCSpectraStep1[index]->GetYaxis()->SetRangeUser(.5,1.5);
1484 MCSpectraStep1[index]->DrawCopy("hist][");
1485 }
1486 else MCSpectraStep1[index]->DrawCopy("hist][same");
1487 }
1488 }
1489 lCompareMCtoMCTruthAtEachStep->DrawClone();
1490 cCompareMCtoMCTruthAtEachStep->cd(3);
1491 for (Int_t icharge=0; icharge<nCharge; icharge++)
1492 {
1493 for (Int_t ipart=0; ipart<nPart; ipart++)
1494 {
1495 Int_t index = ipart + 3*icharge;
1496 MCSpectraStep2[index]->Divide(hMCTruth[index]);
1497 if (icharge == 1) MCSpectraStep2[index]->SetLineStyle(2);
1498 if (index == 0)
1499 {
1500 MCSpectraStep2[index]->GetYaxis()->SetRangeUser(.5,1.5);
1501 MCSpectraStep2[index]->DrawCopy("hist][");
1502 }
1503 else MCSpectraStep2[index]->DrawCopy("hist][same");
1504 }
1505 }
1506 lCompareMCtoMCTruthAtEachStep->DrawClone();
1507 cCompareMCtoMCTruthAtEachStep->cd(4);
1508 for (Int_t icharge=0; icharge<nCharge; icharge++)
1509 {
1510 for (Int_t ipart=0; ipart<nPart; ipart++)
1511 {
1512 Int_t index = ipart + 3*icharge;
1513 MCSpectraStep3[index]->Divide(hMCTruth[index]);
1514 if (icharge == 1) MCSpectraStep3[index]->SetLineStyle(2);
1515 if (index == 0)
1516 {
1517 MCSpectraStep3[index]->GetYaxis()->SetRangeUser(.5,1.5);
1518 MCSpectraStep3[index]->DrawCopy("hist][");
1519 }
1520 else MCSpectraStep3[index]->DrawCopy("hist][same");
1521 }
1522 }
1523 lCompareMCtoMCTruthAtEachStep->DrawClone();
1524 /* cCompareMCtoMCTruthAtEachStep->cd(5);
1525 for (Int_t icharge=0; icharge<nCharge; icharge++)
1526 {
1527 for (Int_t ipart=0; ipart<nPart; ipart++)
1528 {
1529 Int_t index = ipart + 3*icharge;
1530 MCSpectraStep4[index]->Divide(hMCTruth[index]);
1531 if (icharge == 1) MCSpectraStep4[index]->SetLineStyle(2);
1532 if (index == 0)
1533 {
1534 MCSpectraStep4[index]->GetYaxis()->SetRangeUser(.5,1.5);
1535 MCSpectraStep4[index]->DrawCopy("hist][");
1536 }
1537 else MCSpectraStep4[index]->DrawCopy("hist][same");
1538 }
1539 }
1540 lCompareMCtoMCTruthAtEachStep->DrawClone();
1541 cCompareMCtoMCTruthAtEachStep->cd(6);
1542 for (Int_t icharge=0; icharge<nCharge; icharge++)
1543 {
1544 for (Int_t ipart=0; ipart<nPart; ipart++)
1545 {
1546 Int_t index = ipart + 3*icharge;
1547 MCSpectraStep5[index]->Divide(hMCTruth[index]);
1548 if (icharge == 1) MCSpectraStep5[index]->SetLineStyle(2);
1549 if (index == 0)
1550 {
1551 MCSpectraStep5[index]->GetYaxis()->SetRangeUser(.5,1.5);
1552 MCSpectraStep5[index]->DrawCopy("hist][");
1553 }
1554 else MCSpectraStep5[index]->DrawCopy("hist][same");
1555 }
1556 }
1557 lCompareMCtoMCTruthAtEachStep->DrawClone();*/
1558 cCompareMCtoMCTruthAtEachStep->cd(7);
1559 for (Int_t icharge=0; icharge<nCharge; icharge++)
1560 {
1561 for (Int_t ipart=0; ipart<nPart; ipart++)
1562 {
1563 Int_t index = ipart + 3*icharge;
1564 MCSpectraStep6[index]->Divide(hMCTruth[index]);
1565 if (icharge == 1) MCSpectraStep6[index]->SetLineStyle(2);
1566 if (index == 0)
1567 {
1568 MCSpectraStep6[index]->GetYaxis()->SetRangeUser(.5,1.5);
1569 MCSpectraStep6[index]->DrawCopy("hist][");
1570 }
1571 else MCSpectraStep6[index]->DrawCopy("hist][same");
1572 }
1573 }
1574 lCompareMCtoMCTruthAtEachStep->DrawClone();
1575 cCompareMCtoMCTruthAtEachStep->cd(8);
1576 for (Int_t icharge=0; icharge<nCharge; icharge++)
1577 {
1578 for (Int_t ipart=0; ipart<nPart; ipart++)
1579 {
1580 Int_t index = ipart + 3*icharge;
1581 MCSpectra[index]->Divide(hMCTruth[index]);
1582 if (icharge == 1) MCSpectra[index]->SetLineStyle(2);
1583 if (index == 0)
1584 {
1585 MCSpectra[index]->GetYaxis()->SetRangeUser(.5,1.5);
1586 MCSpectra[index]->DrawCopy("hist][");
1587 }
1588 else MCSpectra[index]->DrawCopy("hist][same");
1589 }
1590 }
1591 lCompareMCtoMCTruthAtEachStep->DrawClone();
1592
1593
1594
1595 // close unwanted canvases
1596 // cGFCorrection->Close();
1597 //cMatchingEff->Close();
1598 //cContCorrFact->Close();
1599 //cCorrFactMC->Close();
1600 //cSecCorrFactMC->Close();
1601 cAllCh->Close();
1602 cFitChargHad->Close();
1603 cAllChFactors->Close();
1604
1605 // save wanted canvases
1606 fout->cd();
1607 cSpectra->Write();
1608 CratioComb->Write();
1609 cMatchingEff->Write();
1610 cEvolutionOfSpectra->Write();
1611 cCompareToPublishedAtEachStep->Write();
1612}
1613
1614
1615
1616///////////////////////
1617Double_t eta2y(Double_t pt, Double_t mass, Double_t eta){
1618 Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1619 return TMath::ASinH(pt / mt * TMath::SinH(eta));
1620}
1621
1622///////////////////////
1623void GFCorrection(TH1F **Spectra,AliSpectraAODTrackCuts* tcuts_data, TFile * fout){
1624 //Geant/Fluka Correction
1625 Printf("\nGF correction for Kaons");
1626 //Getting GF For Kaons in TPC
1627 TGraph *gGFCorrectionKaonPlus=new TGraph();
1628 gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
1629 gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
1630 TGraph *gGFCorrectionKaonMinus=new TGraph();
1631 gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
1632 gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
1633 TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
1634 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
1635 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
1636 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
1637 //getting GF func for Kaons with TOF
1638 TF1 *fGFKPosTracking;
1639 fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
1640 TF1 *fGFKNegTracking;
1641 fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
1642 TF1 *fGFKPosMatching;
1643 fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
1644 TF1 *fGFKNegMatching;
1645 fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
1646 for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
1647 if(Spectra[1]->GetBinCenter(binK)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
1648 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
1649 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
1650 Printf("TPC Geant/Fluka: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
1651 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
1652 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
1653 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
1654 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
1655 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
1656 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
1657 }else{
1658 gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
1659 gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
1660 Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
1661 Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
1662 Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
1663 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
1664 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
1665 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
1666 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
1667 Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
1668 Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
1669 Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
1670 Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
1671 Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
1672 Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
1673 Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
1674 }
1675 }
1676
1677 //Geant Fluka for P in TPC
1678 Printf("\nGF correction for Protons");
1679 const Int_t kNCharge=2;
1680 Int_t kPos=0;
1681 Int_t kNeg=1;
1682 TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
1683 TH2D * hCorrFluka[kNCharge];
1684 TH2D * hCorrFluka[2];
1685 hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1686 hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1687 TGraph *gGFCorrectionProtonPlus=new TGraph();
1688 gGFCorrectionProtonPlus->SetName("gGFCorrectionProtonPlus");
1689 gGFCorrectionProtonPlus->SetTitle("gGFCorrectionProtonPlus");
1690 TGraph *gGFCorrectionProtonMinus=new TGraph();
1691 gGFCorrectionProtonMinus->SetName("gGFCorrectionProtonMinus");
1692 gGFCorrectionProtonMinus->SetTitle("gGFCorrectionProtonMinus");
1693 //getting GF func for Protons with TPCTOF
1694 TF1 *fGFpPosTracking;
1695 fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1696 TF1 *fGFpNegTracking;
1697 fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1698 TF1 *fGFpPosMatching;
1699 fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1700 TF1 *fGFpNegMatching;
1701 fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1702
1703 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1704 Int_t nbins = Spectra[2]->GetNbinsX();
1705 Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1706 for(Int_t ibin = 0; ibin < nbins; ibin++){
1707 if(Spectra[2]->GetBinCenter(ibin)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
1708 Float_t pt = Spectra[2]->GetBinCenter(ibin);
1709 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1710 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1711 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
1712 if (pt > maxPtCorrection) pt = maxPtCorrection;
1713 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1714 //cout<<correction<<" charge "<<icharge<<endl;
1715 if(icharge==0){
1716 if (correction != 0) {// If the bin is empty this is a 0
1717 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
1718 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError (ibin)*correction);
1719 gGFCorrectionProtonPlus->SetPoint(ibin,pt,correction);
1720 }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
1721 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, Pt:"<< pt<< endl;
1722 cout << " Bin content: " << Spectra[2]->GetBinContent(ibin) << endl;
1723 }
1724 }
1725 if(icharge==1){
1726 if (correction != 0) {// If the bin is empty this is a 0
1727 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
1728 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError (ibin)*correction);
1729 gGFCorrectionProtonMinus->SetPoint(ibin,pt,correction);
1730 }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
1731 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, Pt:"<< pt<< endl;
1732 cout << " Bin content: " << Spectra[5]->GetBinContent(ibin) << endl;
1733 }
1734 }
1735 }else{
1736 gGFCorrectionProtonPlus->SetPoint(ibin,Spectra[2]->GetBinCenter(ibin),0);
1737 gGFCorrectionProtonMinus->SetPoint(ibin,Spectra[5]->GetBinCenter(ibin),0);
1738 Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1739 Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1740 Printf("TPC/TOF Geant/Fluka Tracking: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1741 if(icharge==0){
1742 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1743 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1744 }else{
1745 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1746 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1747 }
1748 Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1749 Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1750 Printf("TPC/TOF Geant/Fluka Matching: pt:%f Pos:%f Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1751 if(icharge==0){
1752 Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1753 Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1754 }else{
1755 Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1756 Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1757 }
1758 }
1759 }//end loop on bins
1760 }
1761 gGFCorrectionKaonPlus->SetLineColor(kRed);
1762 gGFCorrectionKaonMinus->SetLineColor(kRed+2);
1763 gGFCorrectionProtonPlus->SetLineColor(kGreen);
1764 gGFCorrectionProtonMinus->SetLineColor(kGreen+2);
1765 fGFKPosTracking->SetLineColor(kRed);
1766 fGFKNegTracking->SetLineColor(kRed+2);
1767 fGFKPosMatching->SetLineColor(kRed);
1768 fGFKNegMatching->SetLineColor(kRed+2);
1769 fGFpPosTracking->SetLineColor(kGreen);
1770 fGFpNegTracking->SetLineColor(kGreen+2);
1771 fGFpPosMatching->SetLineColor(kGreen);
1772 fGFpNegMatching->SetLineColor(kGreen+2);
1773 fGFKPosTracking->SetLineStyle(2);
1774 fGFKNegTracking->SetLineStyle(2);
1775 fGFKPosMatching->SetLineStyle(3);
1776 fGFKNegMatching->SetLineStyle(3);
1777 fGFpPosTracking->SetLineStyle(2);
1778 fGFpNegTracking->SetLineStyle(2);
1779 fGFpPosMatching->SetLineStyle(3);
1780 fGFpNegMatching->SetLineStyle(3);
1781 fGFKPosTracking->SetRange(.6,5);
1782 fGFKNegTracking->SetRange(.6,5);
1783 fGFKPosMatching->SetRange(.6,5);
1784 fGFKNegMatching->SetRange(.6,5);
1785 fGFpPosTracking->SetRange(.6,5);
1786 fGFpNegTracking->SetRange(.6,5);
1787 fGFpPosMatching->SetRange(.6,5);
1788 fGFpNegMatching->SetRange(.6,5);
1789
1790 TCanvas * cGFCorrection = new TCanvas ("cGFCorrection","cGFCorrection",700,500);
1791 gPad->SetGridx();
1792 gPad->SetGridy();
1793 gGFCorrectionKaonPlus->DrawClone("al");
1794 gGFCorrectionKaonMinus->DrawClone("lsame");
1795 gGFCorrectionProtonPlus->DrawClone("lsame");
1796 gGFCorrectionProtonMinus->DrawClone("lsame");
1797 fGFKPosTracking->DrawClone("lsame");
1798 fGFKNegTracking->DrawClone("lsame");
1799 fGFKPosMatching->DrawClone("lsame");
1800 fGFKNegMatching->DrawClone("lsame");
1801 fGFpPosTracking->DrawClone("lsame");
1802 fGFpNegTracking->DrawClone("lsame");
1803 fGFpPosMatching->DrawClone("lsame");
1804 fGFpNegMatching->DrawClone("lsame");
1805 gPad->BuildLegend()->SetFillColor(0);
1806 fout->cd();
1807 cGFCorrection->Write();
1808}
1809
1810///////////
1811TF1 * TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1812{
1813
1814 if (ipart == 3 && icharge == kNegative) {
1815 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1816 return f;
1817 }
1818 else if (ipart == 4 && icharge == kNegative) {
1819 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1820 }
1821 else
1822 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1823
1824 return f;
1825}
1826
1827Double_t TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1828{
1829 return 1.;
1830}
1831
1832Double_t TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1833{
1834 return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1835}
1836
1837Double_t TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1838{
1839 return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1840}
1841
1842///////////////////////////////////////////
1843TF1 * TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1844{
1845 if (ipart == 3 && icharge == kNegative) {
1846 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1847 return f;
1848 }
1849 else if (ipart == 4 && icharge == kNegative) {
1850 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1851 }
1852 else
1853 TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1854
1855 return f;
1856}
1857
1858Double_t MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1859{
1860 return 1.;
1861}
1862
1863Double_t MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1864{
1865 Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1866 return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1867}
1868
1869Double_t MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1870{
1871 Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1872 return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1873}