Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxSpectra.cxx
1 #include "AliHighPtDeDxSpectra.h"
2
3 #ifndef __IOSTREAM__
4 #include <iostream>
5 #endif
6
7 using namespace std;
8
9 ClassImp(AliHighPtDeDxSpectra);
10
11 //
12 // AliHighPtDeDxSpectra class
13 //
14 // This class contains the AliHighPtDeDxSpectra information 
15 //
16
17 //_________________________________________________________
18 AliHighPtDeDxSpectra::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 //_________________________________________________________
39 AliHighPtDeDxSpectra::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 //_________________________________________________________
60 AliHighPtDeDxSpectra::~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 //_________________________________________________________
75 TH1D* 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 //_________________________________________________________
103 TH1D* 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 //_________________________________________________________
122 TH1D* 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 //_________________________________________________________
184 TH1D* 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 //_________________________________________________________
195 TF1* 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 //_________________________________________________________
205 void 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 //_________________________________________________________
265 void 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 // }