* Changed the trigger bias correction scheme (added MB->INEL, MB->NSD and MB->ND)
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrection.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaCorrection.h"
4
5 #include <TCanvas.h>
6 #include <TH3F.h>
7 #include <TH1D.h>
8
9 //____________________________________________________________________
10 ClassImp(AlidNdEtaCorrection)
11
12 //____________________________________________________________________
13 AlidNdEtaCorrection::AlidNdEtaCorrection()
14   : TNamed(),
15   fTrack2ParticleCorrection(0),
16   fVertexRecoCorrection(0),
17   fTriggerBiasCorrectionMBToINEL(0),
18   fTriggerBiasCorrectionMBToNSD(0),
19   fTriggerBiasCorrectionMBToND(0)
20 {
21   // default constructor
22 }
23
24 //____________________________________________________________________
25 AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title)
26   : TNamed(name, title),
27   fTrack2ParticleCorrection(0),
28   fVertexRecoCorrection(0),
29   fTriggerBiasCorrectionMBToINEL(0),
30   fTriggerBiasCorrectionMBToNSD(0),
31   fTriggerBiasCorrectionMBToND(0)
32 {
33   // constructor
34   //
35
36   Float_t binLimitsPt[] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 10.0, 100.0};
37
38   TString matrixName;
39   matrixName.Form("%s_nTrackToNPart", name);
40
41   fTrack2ParticleCorrection = new AliCorrectionMatrix3D(matrixName, matrixName, 40, -20, 20, 60, -6, 6, 14, binLimitsPt);
42
43   Float_t binLimitsN[]   = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
44                             10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
45   Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
46
47   matrixName.Form("%s_vtxReco", name);
48   fVertexRecoCorrection        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
49
50   matrixName.Form("%s_triggerBias_MBToINEL", name);
51   fTriggerBiasCorrectionMBToINEL       = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
52   matrixName.Form("%s_triggerBias_MBToNSD", name);
53   fTriggerBiasCorrectionMBToNSD        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
54   matrixName.Form("%s_triggerBias_MBToND", name);
55   fTriggerBiasCorrectionMBToND         = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
56   
57   fTrack2ParticleCorrection      ->SetAxisTitles("vtx z [cm]", "#eta", "p_{T} [GeV/c]");
58   fVertexRecoCorrection          ->SetAxisTitles("vtx z [cm]", "Ntracks");
59   fTriggerBiasCorrectionMBToINEL ->SetAxisTitles("vtx z [cm]", "Ntracks");
60   fTriggerBiasCorrectionMBToNSD  ->SetAxisTitles("vtx z [cm]", "Ntracks");
61   fTriggerBiasCorrectionMBToND   ->SetAxisTitles("vtx z [cm]", "Ntracks");
62
63 }
64
65 //____________________________________________________________________
66 AlidNdEtaCorrection::~AlidNdEtaCorrection()
67 {
68   // destructor
69
70   if (fTrack2ParticleCorrection) {
71     delete fTrack2ParticleCorrection;
72     fTrack2ParticleCorrection = 0;
73   }
74
75   if (fVertexRecoCorrection) {
76     delete fVertexRecoCorrection;
77     fVertexRecoCorrection = 0;
78   }
79
80   if (fTriggerBiasCorrectionMBToINEL) {
81     delete fTriggerBiasCorrectionMBToINEL;
82     fTriggerBiasCorrectionMBToINEL = 0;
83   }
84
85   if (fTriggerBiasCorrectionMBToNSD) {
86     delete fTriggerBiasCorrectionMBToNSD;
87     fTriggerBiasCorrectionMBToNSD = 0;
88   }
89
90   if (fTriggerBiasCorrectionMBToND) {
91     delete fTriggerBiasCorrectionMBToND;
92     fTriggerBiasCorrectionMBToND = 0;
93   }
94 }
95
96 //____________________________________________________________________
97 void
98 AlidNdEtaCorrection::Finish() {
99   //
100   // finish method
101   //
102   // divide the histograms in the AliCorrectionMatrix2D objects to get the corrections
103
104   fTrack2ParticleCorrection->Divide();
105
106   TH3F* hist = fTrack2ParticleCorrection->GetCorrectionHistogram();
107   Int_t emptyBins = 0;
108   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
109     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
110       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
111         if (hist->GetBinContent(x, y, z) == 0)
112         {
113           printf("Empty bin in fTrack2ParticleCorrection at vtx = %f, eta = %f, pt = %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
114           ++emptyBins;
115         }
116
117   printf("INFO: In the central region fTrack2ParticleCorrection has %d empty bins\n", emptyBins);
118
119   fVertexRecoCorrection->Divide();
120   fTriggerBiasCorrectionMBToINEL->Divide();
121   fTriggerBiasCorrectionMBToNSD->Divide();
122   fTriggerBiasCorrectionMBToND->Divide();
123 }
124
125 //____________________________________________________________________
126 Long64_t
127 AlidNdEtaCorrection::Merge(TCollection* list) {
128   // Merge a list of dNdEtaCorrection objects with this (needed for
129   // PROOF).
130   // Returns the number of merged objects (including this).
131
132   if (!list)
133     return 0;
134
135   if (list->IsEmpty())
136     return 1;
137
138   TIterator* iter = list->MakeIterator();
139   TObject* obj;
140
141   // collections of measured and generated histograms
142   TList* collectionNtrackToNparticle    = new TList;
143   TList* collectionVertexReco           = new TList;
144   TList* collectionTriggerBiasMBToINEL  = new TList;
145   TList* collectionTriggerBiasMBToNSD   = new TList;
146   TList* collectionTriggerBiasMBToND    = new TList;
147
148   Int_t count = 0;
149   while ((obj = iter->Next())) {
150
151     AlidNdEtaCorrection* entry = dynamic_cast<AlidNdEtaCorrection*> (obj);
152     if (entry == 0)
153       continue;
154
155     collectionNtrackToNparticle  ->Add(entry->GetTrack2ParticleCorrection());
156     collectionVertexReco         ->Add(entry->GetVertexRecoCorrection());
157     collectionTriggerBiasMBToINEL->Add(entry->GetTriggerBiasCorrection("INEL"));
158     collectionTriggerBiasMBToNSD ->Add(entry->GetTriggerBiasCorrection("NSD"));
159     collectionTriggerBiasMBToND  ->Add(entry->GetTriggerBiasCorrection("ND"));
160
161     count++;
162
163     //fNEvents += entry->fNEvents;
164     //fNTriggeredEvents += entry->fNTriggeredEvents;
165   }
166   fTrack2ParticleCorrection      ->Merge(collectionNtrackToNparticle);
167   fVertexRecoCorrection          ->Merge(collectionVertexReco);
168   fTriggerBiasCorrectionMBToINEL ->Merge(collectionTriggerBiasMBToINEL);
169   fTriggerBiasCorrectionMBToNSD  ->Merge(collectionTriggerBiasMBToNSD);
170   fTriggerBiasCorrectionMBToND   ->Merge(collectionTriggerBiasMBToND);
171
172   delete collectionNtrackToNparticle;
173   delete collectionVertexReco;
174   delete collectionTriggerBiasMBToINEL;
175   delete collectionTriggerBiasMBToNSD;
176   delete collectionTriggerBiasMBToND;
177
178   return count+1;
179 }
180
181
182
183 //____________________________________________________________________
184 Bool_t
185 AlidNdEtaCorrection::LoadHistograms(const Char_t* fileName, const Char_t* dir) {
186   //
187   // loads the histograms
188   //
189
190   fTrack2ParticleCorrection      ->LoadHistograms(fileName, dir);
191   fVertexRecoCorrection          ->LoadHistograms(fileName, dir);
192   fTriggerBiasCorrectionMBToINEL ->LoadHistograms(fileName, dir);
193   fTriggerBiasCorrectionMBToNSD  ->LoadHistograms(fileName, dir);
194   fTriggerBiasCorrectionMBToND   ->LoadHistograms(fileName, dir);
195
196   
197
198   return kTRUE;
199 }
200
201 //____________________________________________________________________
202 void
203 AlidNdEtaCorrection::SaveHistograms() {
204   //
205   // save the histograms
206   //
207
208   gDirectory->mkdir(fName.Data());
209   gDirectory->cd(fName.Data());
210
211   fTrack2ParticleCorrection     ->SaveHistograms();
212   fVertexRecoCorrection         ->SaveHistograms();
213   fTriggerBiasCorrectionMBToINEL->SaveHistograms();
214   fTriggerBiasCorrectionMBToNSD ->SaveHistograms();
215   fTriggerBiasCorrectionMBToND  ->SaveHistograms();
216
217   gDirectory->cd("../");
218 }
219
220 //____________________________________________________________________
221 void AlidNdEtaCorrection::DrawHistograms()
222 {
223   //
224   // call the draw histogram method of the two AliCorrectionMatrix2D objects
225
226   fTrack2ParticleCorrection     ->DrawHistograms();
227   fVertexRecoCorrection         ->DrawHistograms();
228   fTriggerBiasCorrectionMBToINEL->DrawHistograms();
229   fTriggerBiasCorrectionMBToNSD ->DrawHistograms();
230   fTriggerBiasCorrectionMBToND  ->DrawHistograms();
231
232 }
233
234 //____________________________________________________________________
235 void AlidNdEtaCorrection::FillEventWithTrigger(Float_t vtx, Float_t n) 
236 {
237   // fill events with trigger.
238   // used to calculate vertex reco and trigger bias corrections
239
240   fVertexRecoCorrection->FillGene(vtx, n);
241   fTriggerBiasCorrectionMBToINEL ->FillMeas(vtx, n);
242   fTriggerBiasCorrectionMBToNSD  ->FillMeas(vtx, n);
243   fTriggerBiasCorrectionMBToND   ->FillMeas(vtx, n);  
244 }
245
246 //____________________________________________________________________
247 void AlidNdEtaCorrection::FillEventAll(Float_t vtx, Float_t n, Char_t* opt)  
248 {
249   // fill all events 
250   // used to calculate trigger bias corrections
251
252   if (strcmp(opt,"INEL")==0) 
253     fTriggerBiasCorrectionMBToINEL->FillGene(vtx, n);
254   else if (strcmp(opt,"NSD")==0)  
255     fTriggerBiasCorrectionMBToNSD->FillGene(vtx, n);
256   else if (strcmp(opt,"ND")==0)   
257     fTriggerBiasCorrectionMBToND->FillGene(vtx, n);
258   else 
259     AliDebug(AliLog::kWarning, Form(" event type %s unknown (use INEL, NSD or ND)",opt));
260 }
261
262 //____________________________________________________________________
263 AliCorrectionMatrix2D* AlidNdEtaCorrection::GetTriggerBiasCorrection(Char_t* opt)  
264 {
265   // returning the trigger bias correction matrix 
266   // option can be used to specify to which collision process (INEL, NSD or ND)  
267   if (strcmp(opt,"INEL")==0) 
268     return fTriggerBiasCorrectionMBToINEL;
269   else if (strcmp(opt,"NSD")==0)  
270     return fTriggerBiasCorrectionMBToNSD;
271   else if (strcmp(opt,"ND")==0)   
272     return fTriggerBiasCorrectionMBToND;
273   else 
274     {
275       AliDebug(AliLog::kWarning, Form(" %s is unknown (use INEL, NSD or ND). returning INEL ",opt)); 
276       return fTriggerBiasCorrectionMBToINEL;
277     }
278 }
279
280 //____________________________________________________________________
281 Float_t AlidNdEtaCorrection::GetTriggerBiasCorrection(Float_t vtx, Float_t n, Char_t* opt) 
282 {
283   // returning the trigger bias correction matrix 
284   // 3rd option can be used to specify to which collision process (INEL, NSD or ND)  
285
286   if (strcmp(opt,"INEL")==0) 
287     return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);
288   else if (strcmp(opt,"NSD")==0)  
289     return fTriggerBiasCorrectionMBToNSD->GetCorrection(vtx, n);
290   else if (strcmp(opt,"ND")==0)   
291     return fTriggerBiasCorrectionMBToND->GetCorrection(vtx, n);
292   else {
293     AliDebug(AliLog::kWarning, Form(" %s unknown (use INEL, NSD or ND). returning corr for INEL",opt)); 
294     return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);
295   }
296 }
297
298
299 //____________________________________________________________________
300 Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
301 {
302   // calculates the fraction of particles measured (some are missed due to the pt cut off)
303   // uses the generated particle histogram from fTrack2ParticleCorrection
304
305   const TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
306
307   // find eta borders, if eta is negative assume -0.8 ... 0.8
308   Int_t etaBegin = 0;
309   Int_t etaEnd = 0;
310   if (eta < -99)
311   {
312     etaBegin = generated->GetYaxis()->FindBin(-0.8);
313     etaEnd = generated->GetYaxis()->FindBin(0.8);
314   }
315   else
316   {
317     etaBegin = generated->GetYaxis()->FindBin(eta);
318     etaEnd = etaBegin;
319   }
320
321   Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
322   Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
323
324   TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", generated->GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
325   ptProj->GetXaxis()->SetTitle(generated->GetZaxis()->GetTitle());
326
327   Int_t ptBin = ptProj->FindBin(ptCutOff);
328   Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
329   Float_t all = ptProj->Integral();
330
331   if (all == 0)
332     return -1;
333
334   Float_t fraction = abovePtCut / all;
335
336   if (debug)
337   {
338     new TCanvas;
339     ptProj->Draw();
340   }
341   else
342     delete ptProj;
343
344   if (debug)
345     printf("AlidNdEtaCorrection::GetMeasuredFraction: pt cut off = %f, eta = %f, => fraction = %f\n", ptCutOff, eta, fraction);
346
347   return fraction;
348 }
349
350 void AlidNdEtaCorrection::ReduceInformation()
351 {
352   // this function deletes the measured and generated histograms from the corrections to reduce the amount of data
353   // in memory
354
355   // these are needed for GetMeasuredFraction(): fTrack2ParticleCorrection->ReduceInformation();
356   fVertexRecoCorrection          ->ReduceInformation();
357   fTriggerBiasCorrectionMBToINEL ->ReduceInformation();
358   fTriggerBiasCorrectionMBToNSD  ->ReduceInformation();
359   fTriggerBiasCorrectionMBToND   ->ReduceInformation();
360 }
361