- fixing mistake
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix.cxx
1 /* $Id$ */
2
3 // ------------------------------------------------------
4 //
5 // Class to handle corrections.
6 //
7 // ------------------------------------------------------
8 //
9
10 #include <TFile.h>
11 #include <TCanvas.h>
12 #include <TH2F.h>
13
14 #include <AliLog.h>
15
16 #include "AliCorrectionMatrix.h"
17
18 //____________________________________________________________________
19 ClassImp(AliCorrectionMatrix)
20
21 //____________________________________________________________________
22 AliCorrectionMatrix::AliCorrectionMatrix() : TNamed(),
23   fhMeas(0),
24   fhGene(0),
25   fhCorr(0)
26 {
27   // default constructor
28 }
29
30 //____________________________________________________________________
31 AliCorrectionMatrix::AliCorrectionMatrix(const Char_t* name, const Char_t* title) : TNamed(name, title),
32   fhMeas(0),
33   fhGene(0),
34   fhCorr(0)
35 {
36   // constructor initializing tnamed
37 }
38
39 //____________________________________________________________________
40 AliCorrectionMatrix::AliCorrectionMatrix(const AliCorrectionMatrix& c) : TNamed(c),
41   fhMeas(0),
42   fhGene(0),
43   fhCorr(0)
44 {
45   // copy constructor
46   ((AliCorrectionMatrix &)c).Copy(*this);
47 }
48
49 //____________________________________________________________________
50 AliCorrectionMatrix::~AliCorrectionMatrix()
51 {
52   //
53   // destructor
54   //
55
56   if (fhMeas)
57   {
58     delete fhMeas;
59     fhMeas = 0;
60   }
61
62   if (fhGene)
63   {
64     delete fhGene;
65     fhGene = 0;
66   }
67
68   if (fhCorr)
69   {
70     delete fhCorr;
71     fhCorr = 0;
72   }
73 }
74
75 //____________________________________________________________________
76 AliCorrectionMatrix &AliCorrectionMatrix::operator=(const AliCorrectionMatrix &c)
77 {
78   // assigment operator
79
80   if (this != &c)
81     ((AliCorrectionMatrix &) c).Copy(*this);
82
83   return *this;
84 }
85
86 //____________________________________________________________________
87 void AliCorrectionMatrix::Copy(TObject& c) const
88 {
89   // copy function
90
91   AliCorrectionMatrix& target = (AliCorrectionMatrix &) c;
92
93   if (fhMeas)
94     target.fhMeas = dynamic_cast<TH1*> (fhMeas->Clone());
95
96   if (fhGene)
97     target.fhGene = dynamic_cast<TH1*> (fhGene->Clone());
98
99   if (fhCorr)
100     target.fhCorr = dynamic_cast<TH1*> (fhCorr->Clone());
101 }
102
103 //________________________________________________________________________
104 void AliCorrectionMatrix::SetAxisTitles(const Char_t* titleX, const Char_t* titleY, const Char_t* titleZ)
105 {
106   //
107   // method for setting the axis titles of the histograms
108   //
109
110   fhMeas ->SetXTitle(titleX);  fhMeas ->SetYTitle(titleY);  fhMeas ->SetZTitle(titleZ);
111   fhGene ->SetXTitle(titleX);  fhGene ->SetYTitle(titleY);  fhGene ->SetZTitle(titleZ);
112   fhCorr ->SetXTitle(titleX);  fhCorr ->SetYTitle(titleY);  fhCorr ->SetZTitle(titleZ);
113 }
114
115 //____________________________________________________________________
116 Long64_t AliCorrectionMatrix::Merge(TCollection* list)
117 {
118   // Merge a list of AliCorrectionMatrix objects with this (needed for
119   // PROOF). 
120   // Returns the number of merged objects (including this).
121
122   if (!list)
123     return 0;
124   
125   if (list->IsEmpty())
126     return 1;
127
128   TIterator* iter = list->MakeIterator();
129   TObject* obj;
130
131   // collections of measured and generated histograms
132   TList* collectionMeas = new TList;
133   TList* collectionGene = new TList;
134   
135   Int_t count = 0;
136   while ((obj = iter->Next())) {
137     
138     AliCorrectionMatrix* entry = dynamic_cast<AliCorrectionMatrix*> (obj);
139     if (entry == 0) 
140       continue;
141
142     collectionMeas->Add(entry->GetMeasuredHistogram());
143     collectionGene->Add(entry->GetGeneratedHistogram());
144
145     count++;
146   }
147   fhMeas->Merge(collectionMeas);
148   fhGene->Merge(collectionGene);
149
150   delete collectionMeas;
151   delete collectionGene;
152
153   return count+1;
154 }
155
156 //____________________________________________________________________
157 void AliCorrectionMatrix::Divide()
158 {
159   //
160   // divides generated by measured to get the correction
161   //
162
163   if (!fhMeas || !fhGene || !fhCorr) {
164     AliDebug(AliLog::kError, "measured or generated histograms not available");
165     return;
166   }
167
168   fhCorr->Divide(fhGene, fhMeas, 1, 1, "B");
169
170   Int_t emptyBins = 0;
171   for (Int_t x=1; x<=fhCorr->GetNbinsX(); ++x)
172     for (Int_t y=1; y<=fhCorr->GetNbinsY(); ++y)
173       for (Int_t z=1; z<=fhCorr->GetNbinsZ(); ++z)
174         if (fhCorr->GetBinContent(x, y, z) == 0)
175           ++emptyBins;
176
177   if (emptyBins > 0)
178     printf("INFO: In %s we have %d empty bins (of %d) in the correction map\n", GetTitle(), emptyBins, fhCorr->GetNbinsX() * fhCorr->GetNbinsY() * fhCorr->GetNbinsZ());
179 }
180
181 //____________________________________________________________________
182 void AliCorrectionMatrix::Multiply()
183 {
184   //
185   // multiplies measured with correction to get the generated
186   //
187
188   if (!fhMeas || !fhGene || !fhCorr)
189     return;
190
191   fhGene->Multiply(fhMeas, fhCorr, 1, 1, "B");
192 }
193
194 //____________________________________________________________________
195 void AliCorrectionMatrix::Add(AliCorrectionMatrix* aMatrixToAdd, Float_t c) {
196   //
197   // adds the measured and generated of aMatrixToAdd to measured and generated of this
198   // 
199   // NB: the correction will naturally stay the same
200   
201   fhMeas->Add(aMatrixToAdd->GetMeasuredHistogram(), c);
202   fhGene->Add(aMatrixToAdd->GetGeneratedHistogram(), c);
203 }
204
205
206 //____________________________________________________________________
207 Bool_t AliCorrectionMatrix::LoadHistograms(const Char_t* dir)
208 {
209   //
210   // loads the histograms from a file
211   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
212   //
213
214   if (!dir)
215     dir = GetName();
216
217   if (!gDirectory->cd(dir))
218     return kFALSE;
219
220   if (fhGene)
221   {
222     delete fhGene;
223     fhGene=0;
224   }
225
226   if (fhCorr)
227   {
228     delete fhCorr;
229     fhCorr=0;
230   }
231
232   if (fhMeas)
233   {
234     delete fhMeas;
235     fhMeas=0;
236   }
237
238   fhMeas  = dynamic_cast<TH1*> (gDirectory->Get("measured"));
239   if (!fhMeas)
240     Info("LoadHistograms", "No measured hist available");
241
242   fhGene  = dynamic_cast<TH1*> (gDirectory->Get("generated"));
243   if (!fhGene)
244     Info("LoadHistograms", "No generated hist available");
245
246   fhCorr  = dynamic_cast<TH1*> (gDirectory->Get("correction"));
247
248   Bool_t success = kTRUE;
249   if (!fhCorr)
250   {
251     Info("LoadHistograms", "No correction hist available");
252     success = kFALSE;
253   }
254
255   gDirectory->cd("..");
256
257   return success;
258 }
259
260 //____________________________________________________________________
261 void AliCorrectionMatrix::SaveHistograms()
262 {
263   //
264   // saves the histograms
265   //
266
267   gDirectory->mkdir(GetName());
268   gDirectory->cd(GetName());
269
270   if (fhMeas)
271     fhMeas ->Write();
272
273   if (fhGene)
274     fhGene ->Write();
275
276   if (fhCorr)
277     fhCorr->Write();
278
279   gDirectory->cd("..");
280 }
281
282 //____________________________________________________________________
283 void AliCorrectionMatrix::DrawHistograms(const Char_t* canvasName)
284 {
285   //
286   // draws all histograms on one TCanvas
287   // if canvasName is 0 the name of this object is taken
288   //
289
290   if (!canvasName)
291     canvasName = Form("%s_canvas", GetName());
292
293   TCanvas* canvas = new TCanvas(canvasName, GetTitle(), 1200, 400);
294   canvas->Divide(3, 1);
295
296   canvas->cd(1);
297   if (fhMeas)
298     fhMeas->Draw("COLZ");
299
300   canvas->cd(2);
301   if (fhGene)
302   {
303     // work around ROOT bug #22011
304     if (fhGene->GetEntries() == 0)
305       fhGene->SetEntries(1);
306     fhGene->Draw("COLZ");
307   }
308
309   canvas->cd(3);
310   if (fhCorr)
311     fhCorr->Draw("COLZ");
312 }
313
314 //____________________________________________________________________
315 void AliCorrectionMatrix::ReduceInformation()
316 {
317   // this function deletes the measured and generated histograms to reduce the amount of data
318   // in memory
319
320   if (fhMeas)
321   {
322     delete fhMeas;
323     fhMeas = 0;
324   }
325
326   if (fhGene)
327   {
328     delete fhGene;
329     fhGene = 0;
330   }
331 }
332
333 //____________________________________________________________________
334 void AliCorrectionMatrix::Reset(Option_t* option)
335 {
336   // resets the histograms
337
338   if (fhGene)
339     fhGene->Reset(option);
340
341   if (fhMeas)
342     fhMeas->Reset(option);
343
344   if (fhCorr)
345     fhCorr->Reset(option);
346 }
347
348 //____________________________________________________________________
349 void AliCorrectionMatrix::SetCorrectionToUnity()
350 {
351   // sets the correction matrix to unity
352
353   if (!fhCorr)
354     return;
355
356   for (Int_t x=0; x<=fhCorr->GetNbinsX()+1; ++x)
357     for (Int_t y=0; y<=fhCorr->GetNbinsY()+1; ++y)
358       for (Int_t z=0; z<=fhCorr->GetNbinsZ()+1; ++z)
359       {
360         fhCorr->SetBinContent(x, y, z, 1);
361         fhCorr->SetBinError(x, y, z, 0);
362       }
363 }
364
365 //____________________________________________________________________
366 void AliCorrectionMatrix::Scale(Double_t factor)
367 {
368   // scales the generated and measured histogram with the given factor
369
370   Printf("Scaling histograms with %f", factor);
371
372   fhMeas->Scale(factor);
373   fhGene->Scale(factor);
374 }
375
376 //____________________________________________________________________
377 void AliCorrectionMatrix::ResetErrorsOnCorrections()
378 {
379   // set the errors on the correction matrix to 0
380
381   if (!fhCorr)
382     return;
383
384   for (Int_t x=0; x<=fhCorr->GetNbinsX()+1; ++x)
385     for (Int_t y=0; y<=fhCorr->GetNbinsY()+1; ++y)
386       for (Int_t z=0; z<=fhCorr->GetNbinsZ()+1; ++z)
387         fhCorr->SetBinError(x, y, z, 0);
388 }