Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxSpectra.cxx
CommitLineData
4ebdd20e 1#include "AliHighPtDeDxSpectra.h"
2
3#ifndef __IOSTREAM__
4#include <iostream>
5#endif
6
7using namespace std;
8
9ClassImp(AliHighPtDeDxSpectra);
10
11//
12// AliHighPtDeDxSpectra class
13//
14// This class contains the AliHighPtDeDxSpectra information
15//
16
17//_________________________________________________________
18AliHighPtDeDxSpectra::AliHighPtDeDxSpectra():
19 TNamed(),
20 fDebugLevel(0),
21 fUseMcNoVtxCorrection(kTRUE),
22 fUseFittedEfficiency(kTRUE),
23 fUseBinCorrection(kTRUE),
24 fNevents(0),
25 hPt(0x0),
26 hNevents(0x0),
27 hTriggerEfficiency(0x0),
28 hEfficiency(0x0),
29 hMcNoVtxCorrection(0x0),
30 hBinCorrection(0x0),
31 hMeanPt(0x0),
32 fEfficiency(0x0),
33 fBinFit(0x0)
34{
35 // default constructor - do not use
36}
37
38//_________________________________________________________
39AliHighPtDeDxSpectra::AliHighPtDeDxSpectra(const char* name, const char* title):
40 TNamed(name, title),
41 fDebugLevel(0),
42 fUseMcNoVtxCorrection(kTRUE),
43 fUseFittedEfficiency(kTRUE),
44 fUseBinCorrection(kTRUE),
45 fNevents(0),
46 hPt(0x0),
47 hNevents(0x0),
48 hTriggerEfficiency(0x0),
49 hEfficiency(0x0),
50 hMcNoVtxCorrection(0x0),
51 hBinCorrection(0x0),
52 hMeanPt(0x0),
53 fEfficiency(0x0),
54 fBinFit(0x0)
55{
56 // named constructor
57}
58
59//_________________________________________________________
60AliHighPtDeDxSpectra::~AliHighPtDeDxSpectra()
61{
62 // histograms
63 delete hPt; // pt spectrum for
64 delete hNevents; // events
65 delete hTriggerEfficiency; // trigger efficency
66 delete hEfficiency; // track correction
67 delete hMcNoVtxCorrection;
68 delete hBinCorrection; // bin correction
69 delete hMeanPt; // mean pt of data
70 delete fEfficiency;
71 delete fBinFit;
72}
73
74//_________________________________________________________
75TH1D* AliHighPtDeDxSpectra::ConstructBinCorrection(TH1D* histPt, TProfile* histMeanPt)
76{
77 TH1D* histCorr = dynamic_cast<TH1D*>(histPt->Clone("hBinCorrection"));
78 histCorr->SetDirectory(0);
79 histCorr->Reset();
80
81 fBinFit = new TF1("fBinFit", "[0]*(1.0+x/[1])**(-[2])", 0.0, 50);
82 fBinFit->SetParameters(6, 0.5, 4);
83 histPt->Fit(fBinFit, "0N", "", 3.0, 50.0);
84
85 const Int_t nPtBins = histPt->GetNbinsX();
86 for(Int_t bin = 1; bin < nPtBins; bin++) {
87
88 Float_t meanPt = histMeanPt->GetBinContent(bin);
89 Float_t binPt = histPt->GetXaxis()->GetBinCenter(bin);
90 Float_t correction = 0;
91 if(meanPt>0)
92 correction = fBinFit->Eval(meanPt)/fBinFit->Eval(binPt);
93 if(fDebugLevel>5)
94 cout << "<pT>: " << meanPt << ", bin pT: " << binPt << ", correction: " << correction << endl;
95 histCorr->SetBinContent(bin, correction);
96 histCorr->SetBinError(bin, 0.0);
97 }
98
99 return histCorr;
100}
101
102//_________________________________________________________
103TH1D* AliHighPtDeDxSpectra::ConstructTriggerEfficiency(AliHighPtDeDxMc* mc)
104{
105 TH1D* histEff = dynamic_cast<TH1D*>(mc->GetHistNeventsMcTrig()->Clone("hTriggerEfficiency"));
106 histEff->SetDirectory(0);
107 histEff->Divide(histEff, mc->GetHistNeventsMc(), 1, 1, "B");
108
109 if(fDebugLevel>1) {
110
111 cout << "Trigger efficiency: "
112 << Form("%.1f %% (No vtx), %.1f %% (vtx)",
113 Float_t(histEff->GetBinContent(1)*100.0),
114 Float_t(histEff->GetBinContent(2)*100.0))
115 << endl;
116 }
117
118 return histEff;
119}
120
121//_________________________________________________________
122TH1D* AliHighPtDeDxSpectra::GetEventCount(AliHighPtDeDxData* data, AliHighPtDeDxMc* mc)
123{
124 // Correct
125 TH1D* histEvents = dynamic_cast<TH1D*>(data->GetHistNevents()->Clone("hNevents"));
126
127 Double_t noVtx = histEvents->GetBinContent(1);
128 // Double_t vtx = histEvents->GetBinContent(2);
129 Double_t vtx = histEvents->GetBinContent(2);
130
131 if(fDebugLevel>1) {
132 cout << "BEFORE correction: " << noVtx << " (no vtx), " << vtx << " (vtx)" << endl;
133 }
134
135 if(!fUseMcNoVtxCorrection) {
136
137 TH1D* histVtxStatus = data->GetHistVtxStatus();
138
139 Double_t outsideCut = histVtxStatus->GetBinContent(2);
140 Double_t insideCut = histVtxStatus->GetBinContent(3);
141 Double_t ratio = insideCut/(insideCut+outsideCut);
142
143 if(fDebugLevel>1) {
144 cout << "Inside/Outside/Ratio(I/(I+O)):" << insideCut << "/" << outsideCut << "/" << ratio*100 << " %%" << endl;
145 }
146
147 noVtx *= ratio;
148 histEvents->SetBinContent(1, noVtx);
149
150 if(fDebugLevel>1) {
151 cout << ", AFTER correction: " << noVtx << endl;
152 }
153 } else {
154
155 hMcNoVtxCorrection = dynamic_cast<TH1D*>(mc->GetHistNeventsMcTrig()->Clone("hMcEventCorrection"));
156 hMcNoVtxCorrection->Divide(mc->GetHistNevents());
157
158 Double_t scaleNoVtx = hMcNoVtxCorrection->GetBinContent(1);
159 Double_t scaleVtx = hMcNoVtxCorrection->GetBinContent(2);
160
161 if(fDebugLevel>1) {
162 cout << "Mc event correction: "
163 << Form("%.1f %% (No vtx), %.1f %% (vtx - not corrected)",
164 Float_t(scaleNoVtx*100.0),
165 Float_t(scaleVtx*100.0))
166 << endl;
167 }
168
169 noVtx *= scaleNoVtx;
170 // vtx *= scaleVtx;
171 histEvents->SetBinContent(1, noVtx);
172 // histEvents->SetBinContent(2, vtx);
173 }
174
175 if(fDebugLevel>1) {
176 cout << "AFTER correction: " << noVtx << " (no vtx), " << vtx << " (vtx)" << endl;
177 }
178
179 return histEvents;
180}
181
182
183//_________________________________________________________
184TH1D* AliHighPtDeDxSpectra::ConstructTrackCorrection(AliHighPtDeDxMc* mc)
185{
186 TH1D* histEff = dynamic_cast<TH1D*>(mc->GetHistPt()->Clone("hEfficiency"));
187 histEff->SetDirectory(0);
188 histEff->Divide(histEff, mc->GetHistPtMc(0, 0), 1, 1, "B");
189 // histEff->Divide(histEff, mc->GetHistPtMc());
190
191 return histEff;
192}
193
194//_________________________________________________________
195TF1* AliHighPtDeDxSpectra::FitEfficiency(TH1D* histEff)
196{
197 TF1* func = new TF1("fEfficiency", "[0]*(1-[1]/x)", 0, 50);
198
199 histEff->Fit(func, "0N", "", 3, 50);
200
201 return func;
202}
203
204//_________________________________________________________
205void AliHighPtDeDxSpectra::GetCorrectedSpectra(AliHighPtDeDxData* data,
206 AliHighPtDeDxMc* mc)
207{
208 //
209 // Noralize spectra
210 //
211
212 //
213 // Extract the corrected number of events
214 //
215 Double_t fNeventsVtx = data->GetHistNevents()->GetEntries() - data->GetHistNevents()->GetBinContent(1);
216 // hNevents = GetEventCount(data, mc);
217 // hTriggerEfficiency = ConstructTriggerEfficiency(mc);
218 // hNevents->Divide(hTriggerEfficiency);
219
220 // if(fDebugLevel>1) {
221 // cout << "Events after efficiency: " << hNevents->GetBinContent(1)
222 // << " (no vtx) " << hNevents->GetBinContent(2) << " (vtx)" << endl;
223 // }
224
225 // fNevents = hNevents->GetBinContent(1) + hNevents->GetBinContent(2);
226 fNevents = fNeventsVtx/0.85;
227
228 //
229 // Extract the track correction
230 //
231 hEfficiency = ConstructTrackCorrection(mc);
232
233 hPt = dynamic_cast<TH1D*>(data->GetHistPt()->Clone("hPt"));
234 hPt->SetDirectory(0);
235
236 Double_t etaRange = data->GetEtaHigh() - data->GetEtaLow();
237
238 if(fDebugLevel>1) {
239 cout << "N events (corrected): " << fNevents << endl;
240 cout << "Eta range: " << etaRange << endl;
241
242 cout << endl << "fNeventsVtx(" << fNeventsVtx << ")/fNevents(" << fNevents << ") = "
243 << fNeventsVtx/fNevents*100 << " %" << endl << endl;
244 }
245
246 hPt->Scale(1.0/etaRange);
247 NormalizePt(hPt);
248 hPt->Scale(1.0/fNevents);
249
250 if (fUseFittedEfficiency) {
251
252 fEfficiency = FitEfficiency(hEfficiency);
253 hPt->Divide(fEfficiency);
254 } else
255 hPt->Divide(hEfficiency);
256
257 if (fUseBinCorrection) {
258
259 hBinCorrection = ConstructBinCorrection(hPt, data->GetHistMeanPt());
260 hPt->Divide(hBinCorrection);
261 }
262}
263
264//_________________________________________________________
265void AliHighPtDeDxSpectra::NormalizePt(TH1* hist)
266{
267 const Int_t n = hist->GetNbinsX();
268
269 for(Int_t bin = 1; bin <= n; bin++) {
270
271 const Float_t width = hist->GetXaxis()->GetBinWidth(bin);
272 hist->SetBinError(bin, hist->GetBinError(bin)/width);
273 hist->SetBinContent(bin, hist->GetBinContent(bin)/width);
274 }
275}
276
277
278// //_________________________________________________________
279// void AliHighPtDeDxSpectra::FillTrackInfo(Float_t weight)
280// {
281
282// AliHighPtDeDxBase::FillTrackInfo(weight);
283
284// hPt->Fill(fTrackPt, weight);
285// }