7b7f069054a9da68299e9101902c161d35ce32ee
[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   fTriggerCorrection(0),
18   fTriggerBiasCorrection(0)
19   //fNEvents(0),
20   //fNTriggeredEvents(0)
21 {
22   // default constructor
23 }
24
25 //____________________________________________________________________
26 AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title)
27   : TNamed(name, title),
28   fTrack2ParticleCorrection(0),
29   fVertexRecoCorrection(0),
30   fTriggerCorrection(0),
31   fTriggerBiasCorrection(0)
32   //fNEvents(0),
33   //fNTriggeredEvents(0)
34 {
35   // constructor
36   //
37
38   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};
39
40   TString matrixName;
41   matrixName.Form("%s_nTrackToNPart", name);
42
43   fTrack2ParticleCorrection = new AliCorrectionMatrix3D(matrixName, matrixName, 40, -20, 20, 60, -6, 6, 14, binLimitsPt);
44
45   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,
46                             10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
47   Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
48
49   matrixName.Form("%s_vtxReco", name);
50   fVertexRecoCorrection        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
51
52   matrixName.Form("%s_trigger", name);
53   fTriggerCorrection           = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
54
55   matrixName.Form("%s_triggerBias", name);
56   fTriggerBiasCorrection       = new AliCorrectionMatrix2D(matrixName, matrixName, 120,-6,6,100, 0, 10);
57
58   fTrack2ParticleCorrection ->SetAxisTitles("vtx z [cm]", "#eta", "p_{T} [GeV/c]");
59   fVertexRecoCorrection        ->SetAxisTitles("vtx z [cm]", "Ntracks");
60   fTriggerCorrection        ->SetAxisTitles("vtx z [cm]", "Ntracks");
61
62   fTriggerBiasCorrection       ->SetAxisTitles("#eta", "p_{T} [GeV/c]");
63 }
64
65 //____________________________________________________________________
66 AlidNdEtaCorrection::~AlidNdEtaCorrection()
67 {
68   // destructor
69
70   if (fTrack2ParticleCorrection)
71   {
72     delete fTrack2ParticleCorrection;
73     fTrack2ParticleCorrection = 0;
74   }
75
76   if (fVertexRecoCorrection)
77   {
78     delete fVertexRecoCorrection;
79     fVertexRecoCorrection = 0;
80   }
81
82   if (fTriggerCorrection)
83   {
84     delete fTriggerCorrection;
85     fTriggerCorrection = 0;
86   }
87
88   if (fTriggerBiasCorrection)
89   {
90     delete fTriggerBiasCorrection;
91     fTriggerBiasCorrection = 0;
92   }
93 }
94
95 //____________________________________________________________________
96 void
97 AlidNdEtaCorrection::Finish() {
98   //
99   // finish method
100   //
101   // divide the histograms in the AliCorrectionMatrix2D objects to get the corrections
102
103   fTrack2ParticleCorrection->Divide();
104
105   TH3F* hist = fTrack2ParticleCorrection->GetCorrectionHistogram();
106   Int_t emptyBins = 0;
107   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
108     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
109       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
110         if (hist->GetBinContent(x, y, z) == 0)
111         {
112           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));
113           ++emptyBins;
114         }
115
116   printf("INFO: In the central region fTrack2ParticleCorrection has %d empty bins\n", emptyBins);
117
118   fVertexRecoCorrection->Divide();
119   fTriggerCorrection->Divide();
120
121   if (GetNevents() <= 0)
122   {
123     printf("ERROR: GetNevents() returns 0. Cannot scale histogram. Skipping processing of fTriggerBiasCorrection\n");
124     return;
125   }
126   fTriggerBiasCorrection->GetMeasuredHistogram()->Scale(Double_t(GetNtriggeredEvents())/Double_t(GetNevents()));
127   fTriggerBiasCorrection->Divide();
128 }
129
130 //____________________________________________________________________
131 Long64_t
132 AlidNdEtaCorrection::Merge(TCollection* list) {
133   // Merge a list of dNdEtaCorrection objects with this (needed for
134   // PROOF).
135   // Returns the number of merged objects (including this).
136
137   if (!list)
138     return 0;
139
140   if (list->IsEmpty())
141     return 1;
142
143   TIterator* iter = list->MakeIterator();
144   TObject* obj;
145
146   // collections of measured and generated histograms
147   TList* collectionNtrackToNparticle = new TList;
148   TList* collectionVertexReco        = new TList;
149   TList* collectionTriggerBias       = new TList;
150   TList* collectionTrigger           = new TList;
151
152   Int_t count = 0;
153   while ((obj = iter->Next())) {
154
155     AlidNdEtaCorrection* entry = dynamic_cast<AlidNdEtaCorrection*> (obj);
156     if (entry == 0)
157       continue;
158
159     collectionNtrackToNparticle ->Add(entry->GetTrack2ParticleCorrection());
160     collectionVertexReco        ->Add(entry->GetVertexRecoCorrection());
161     collectionTriggerBias        ->Add(entry->GetTriggerBiasCorrection());
162     collectionTrigger        ->Add(entry->GetTriggerCorrection());
163
164     count++;
165
166     //fNEvents += entry->fNEvents;
167     //fNTriggeredEvents += entry->fNTriggeredEvents;
168   }
169   fTrack2ParticleCorrection ->Merge(collectionNtrackToNparticle);
170   fVertexRecoCorrection        ->Merge(collectionVertexReco);
171   fTriggerBiasCorrection        ->Merge(collectionTriggerBias);
172   fTriggerCorrection        ->Merge(collectionTrigger);
173
174   delete collectionNtrackToNparticle;
175   delete collectionVertexReco;
176   delete collectionTriggerBias;
177   delete collectionTrigger;
178
179   return count+1;
180 }
181
182 //____________________________________________________________________
183 Bool_t
184 AlidNdEtaCorrection::LoadHistograms(const Char_t* fileName, const Char_t* dir) {
185   //
186   // loads the histograms
187   //
188
189   fTrack2ParticleCorrection ->LoadHistograms(fileName, dir);
190   fVertexRecoCorrection        ->LoadHistograms(fileName, dir);
191   fTriggerCorrection        ->LoadHistograms(fileName, dir);
192   fTriggerBiasCorrection       ->LoadHistograms(fileName, dir);
193
194   return kTRUE;
195 }
196
197 //____________________________________________________________________
198 void
199 AlidNdEtaCorrection::SaveHistograms() {
200   //
201   // save the histograms
202   //
203
204   gDirectory->mkdir(fName.Data());
205   gDirectory->cd(fName.Data());
206
207   fTrack2ParticleCorrection->SaveHistograms();
208   fVertexRecoCorrection->SaveHistograms();
209   fTriggerCorrection->SaveHistograms();
210   fTriggerBiasCorrection->SaveHistograms();
211
212   gDirectory->cd("../");
213 }
214
215 //____________________________________________________________________
216 void AlidNdEtaCorrection::DrawHistograms()
217 {
218   //
219   // call the draw histogram method of the two AliCorrectionMatrix2D objects
220
221   fTrack2ParticleCorrection ->DrawHistograms();
222   fVertexRecoCorrection        ->DrawHistograms();
223   fTriggerCorrection        ->DrawHistograms();
224   fTriggerBiasCorrection       ->DrawHistograms();
225
226 }
227
228 //____________________________________________________________________
229 Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
230 {
231   // calculates the fraction of particles measured (some are missed due to the pt cut off)
232   // uses the generated particle histogram from fTrack2ParticleCorrection
233
234   const TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
235
236   // find eta borders, if eta is negative assume -0.8 ... 0.8
237   Int_t etaBegin = 0;
238   Int_t etaEnd = 0;
239   if (eta < -99)
240   {
241     etaBegin = generated->GetYaxis()->FindBin(-0.8);
242     etaEnd = generated->GetYaxis()->FindBin(0.8);
243   }
244   else
245   {
246     etaBegin = generated->GetYaxis()->FindBin(eta);
247     etaEnd = etaBegin;
248   }
249
250   Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
251   Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
252
253   TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", generated->GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
254   ptProj->GetXaxis()->SetTitle(generated->GetZaxis()->GetTitle());
255
256   Int_t ptBin = ptProj->FindBin(ptCutOff);
257   Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
258   Float_t all = ptProj->Integral();
259
260   if (all == 0)
261     return -1;
262
263   Float_t fraction = abovePtCut / all;
264
265   if (debug)
266   {
267     new TCanvas;
268     ptProj->Draw();
269   }
270   else
271     delete ptProj;
272
273   if (debug)
274     printf("AlidNdEtaCorrection::GetMeasuredFraction: pt cut off = %f, eta = %f, => fraction = %f\n", ptCutOff, eta, fraction);
275
276   return fraction;
277 }
278
279 void AlidNdEtaCorrection::ReduceInformation()
280 {
281   // this function deletes the measured and generated histograms from the corrections to reduce the amount of data
282   // in memory
283
284   // these are needed for GetMeasuredFraction(): fTrack2ParticleCorrection->ReduceInformation();
285   fVertexRecoCorrection->ReduceInformation();
286   fTriggerCorrection->ReduceInformation();
287   fTriggerBiasCorrection->ReduceInformation();
288 }
289
290 Int_t AlidNdEtaCorrection::GetNevents()
291 {
292   // Return total number of events used for this correction map (taken from fTriggerCorrection)
293
294   if (!fTriggerCorrection)
295     return -1;
296
297   return (Int_t) fTriggerCorrection->GetGeneratedHistogram()->Integral(0, fTriggerCorrection->GetGeneratedHistogram()->GetNbinsX()+1);
298 }
299
300 Int_t AlidNdEtaCorrection::GetNtriggeredEvents()
301 {
302   // Return total number of triggered events used for this correction map (taken from fTriggerCorrection)
303
304   if (!fTriggerCorrection)
305     return -1;
306
307   return (Int_t) fTriggerCorrection->GetMeasuredHistogram()->Integral(0, fTriggerCorrection->GetMeasuredHistogram()->GetNbinsX()+1);
308 }