replaced THXF to THX in many function prototypes
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix3D.cxx
1 /* $Id$ */
2
3 // ------------------------------------------------------
4 //
5 // Class to handle 3d-corrections.
6 //
7 // ------------------------------------------------------
8 //
9
10 #include <TDirectory.h>
11 #include <TH1F.h>
12 #include <TH2F.h>
13 #include <TH3F.h>
14 #include <TString.h>
15
16 #include <AliLog.h>
17
18 #include "AliCorrectionMatrix3D.h"
19 #include "AliCorrectionMatrix2D.h"
20 #include "AliPWG0Helper.h"
21
22 //____________________________________________________________________
23 ClassImp(AliCorrectionMatrix3D)
24
25 //____________________________________________________________________
26 AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
27   AliCorrectionMatrix()
28 {
29   // default constructor
30 }
31
32 //____________________________________________________________________
33 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
34   : AliCorrectionMatrix(c)
35 {
36   // copy constructor
37   ((AliCorrectionMatrix3D &)c).Copy(*this);
38 }
39
40 //____________________________________________________________________
41 AliCorrectionMatrix3D &AliCorrectionMatrix3D::operator=(const AliCorrectionMatrix3D &c)
42 {
43   // assigment operator
44
45   if (this != &c)
46     ((AliCorrectionMatrix3D &) c).Copy(*this);
47
48   return *this;
49 }
50
51 //____________________________________________________________________
52 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
53               Int_t nBinX, Float_t Xmin, Float_t Xmax,
54               Int_t nBinY, Float_t Ymin, Float_t Ymax,
55               Int_t nBinZ, Float_t Zmin, Float_t Zmax)
56   : AliCorrectionMatrix(name, title)
57 {
58   //
59   // constructor
60   //
61
62   Float_t* binLimitsX = new Float_t[nBinX+1];
63   for (Int_t i=0; i<=nBinX; ++i)
64     binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
65
66   Float_t* binLimitsY = new Float_t[nBinY+1];
67   for (Int_t i=0; i<=nBinY; ++i)
68     binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
69
70   Float_t* binLimitsZ = new Float_t[nBinZ+1];
71   for (Int_t i=0; i<=nBinZ; ++i)
72     binLimitsZ[i] = Zmin + (Zmax - Zmin) / nBinZ * i;
73
74   CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
75
76   delete[] binLimitsX;
77   delete[] binLimitsY;
78   delete[] binLimitsZ;
79 }
80
81 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
82       Int_t nBinX, Float_t Xmin, Float_t Xmax,
83       Int_t nBinY, Float_t Ymin, Float_t Ymax,
84       Int_t nBinZ, const Float_t* zbins)
85   : AliCorrectionMatrix(name, title)
86 {
87   // constructor with variable bin sizes
88
89   Float_t* binLimitsX = new Float_t[nBinX+1];
90   for (Int_t i=0; i<=nBinX; ++i)
91     binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
92
93   Float_t* binLimitsY = new Float_t[nBinY+1];
94   for (Int_t i=0; i<=nBinY; ++i)
95     binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
96
97   CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
98
99   delete[] binLimitsX;
100   delete[] binLimitsY;
101 }
102
103 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3* hBinning)
104   : AliCorrectionMatrix(name, title)
105 {
106   // constructor with variable bin sizes (uses binning of hBinning)
107
108   // do not add this hists to the directory
109   Bool_t oldStatus = TH1::AddDirectoryStatus();
110   TH1::AddDirectory(kFALSE);
111
112   fhMeas = (TH3F*)hBinning->Clone("measured");
113   fhGene = (TH3F*)hBinning->Clone("generated");
114   fhCorr = (TH3F*)hBinning->Clone("correction");
115
116   fhMeas->SetTitle(Form("%s measured", GetTitle()));
117   fhGene->SetTitle(Form("%s generated", GetTitle()));
118   fhCorr->SetTitle(Form("%s correction", GetTitle()));
119
120   fhMeas->Reset();
121   fhGene->Reset();
122   fhCorr->Reset();
123
124   TH1::AddDirectory(oldStatus);
125
126   fhMeas->Sumw2();
127   fhGene->Sumw2();
128   fhCorr->Sumw2();
129 }
130
131 //____________________________________________________________________
132 void AliCorrectionMatrix3D::CreateHists(Int_t nBinX, const Float_t* binLimitsX,
133       Int_t nBinY, const Float_t* binLimitsY,
134       Int_t nBinZ, const Float_t* binLimitsZ)
135 {
136   // create the histograms
137
138   // do not add this hists to the directory
139   Bool_t oldStatus = TH1::AddDirectoryStatus();
140   TH1::AddDirectory(kFALSE);
141
142   fhMeas  = new TH3F("measured",   Form("%s measured",GetTitle()),   nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
143   fhGene  = new TH3F("generated",  Form("%s generated",GetTitle()),  nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
144   fhCorr  = new TH3F("correction", Form("%s correction",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
145
146   fhMeas->Sumw2();
147   fhGene->Sumw2();
148   fhCorr->Sumw2();
149
150   TH1::AddDirectory(oldStatus);
151 }
152
153
154 //____________________________________________________________________
155 AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
156 {
157   //
158   // destructor
159   //
160
161   // histograms already deleted in base class
162 }
163
164 //____________________________________________________________________
165 TH3* AliCorrectionMatrix3D::GetGeneratedHistogram()
166 {
167   // return generated histogram casted to correct type
168   return dynamic_cast<TH3F*> (fhGene);
169 }
170
171 //____________________________________________________________________
172 TH3* AliCorrectionMatrix3D::GetMeasuredHistogram()
173 {
174   // return measured histogram casted to correct type
175   return dynamic_cast<TH3F*> (fhMeas);
176 }
177
178 //____________________________________________________________________
179 TH3* AliCorrectionMatrix3D::GetCorrectionHistogram()
180 {
181   // return correction histogram casted to correct type
182   return dynamic_cast<TH3F*> (fhCorr);
183 }
184
185 //____________________________________________________________________
186 AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Char_t* opt, Float_t aMin, Float_t aMax)
187 {
188   // returns a 2D projection of this correction
189
190   TString option = opt;
191
192   // unzoom
193   fhMeas->GetXaxis()->SetRange(0, 0);
194   fhMeas->GetYaxis()->SetRange(0, 0);
195   fhMeas->GetZaxis()->SetRange(0, 0);
196
197   fhGene->GetXaxis()->SetRange(0, 0);
198   fhGene->GetYaxis()->SetRange(0, 0);
199   fhGene->GetZaxis()->SetRange(0, 0);
200
201   if (aMin<aMax) {
202     if (option.Contains("xy") || option.Contains("yx")) {
203       Int_t bMin = fhMeas->GetZaxis()->FindBin(aMin);
204       Int_t bMax = fhMeas->GetZaxis()->FindBin(aMax);
205       fhGene->GetZaxis()->SetRange(bMin, bMax);
206       fhMeas->GetZaxis()->SetRange(bMin, bMax);
207     }
208     else if (option.Contains("xz") || option.Contains("zx")) {
209       Int_t bMin = fhMeas->GetYaxis()->FindBin(aMin);
210       Int_t bMax = fhMeas->GetYaxis()->FindBin(aMax);
211       fhGene->GetYaxis()->SetRange(bMin, bMax);
212       fhMeas->GetYaxis()->SetRange(bMin, bMax);
213     }
214     else if (option.Contains("yz") || option.Contains("zy")) {
215       Int_t bMin = fhMeas->GetXaxis()->FindBin(aMin);
216       Int_t bMax = fhMeas->GetXaxis()->FindBin(aMax);
217       fhGene->GetXaxis()->SetRange(bMin, bMax);
218       fhMeas->GetXaxis()->SetRange(bMin, bMax);
219     }
220     else {
221       AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
222       return 0;
223     }
224   }
225
226   AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
227
228   TH2F* meas = (TH2F*) ((TH3F*)fhMeas)->Project3D(option)->Clone(Form("%s_meas", corr2D->GetName()));
229   TH2F* gene = (TH2F*) ((TH3F*)fhGene)->Project3D(option)->Clone(Form("%s_gene", corr2D->GetName()));
230
231   TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
232   corr->Reset();
233
234   corr2D->SetGeneratedHistogram(gene);
235   corr2D->SetMeasuredHistogram(meas);
236   corr2D->SetCorrectionHistogram(corr);
237
238   corr2D->Divide();
239
240   // unzoom
241   fhMeas->GetXaxis()->SetRange(0, 0);
242   fhMeas->GetYaxis()->SetRange(0, 0);
243   fhMeas->GetZaxis()->SetRange(0, 0);
244
245   fhGene->GetXaxis()->SetRange(0, 0);
246   fhGene->GetYaxis()->SetRange(0, 0);
247   fhGene->GetZaxis()->SetRange(0, 0);
248
249   return corr2D;
250 }
251
252 //____________________________________________________________________
253 TH1* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Char_t* opt, Float_t aMin1, Float_t aMax1, Float_t aMin2, Float_t aMax2)
254 {
255   // returns a 1D projection of this correction
256   AliDebug(AliLog::kWarning, Form("WARNING: test"));
257   
258   AliCorrectionMatrix2D* corr2D;
259   if (strcmp(opt,"x")==0) {
260     corr2D = Get2DCorrection("yx",aMin1,aMax1);
261     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
262   }
263   if (strcmp(opt,"y")==0) {
264     corr2D = Get2DCorrection("xy",aMin1,aMax1);
265     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
266   }  
267   if (strcmp(opt,"z")==0) {
268     corr2D = Get2DCorrection("yz",aMin1,aMax1);
269     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
270   }  
271   AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));
272   
273   return 0;
274 }
275
276 //____________________________________________________________________
277 void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az)
278 {
279   // add value to measured histogram
280   GetMeasuredHistogram()->Fill(ax, ay, az);
281 }
282
283 //____________________________________________________________________
284 void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az)
285 {
286   // add value to generated histogram
287   GetGeneratedHistogram()->Fill(ax, ay, az);
288 }
289
290 //____________________________________________________________________
291 Float_t AliCorrectionMatrix3D::GetCorrection(Float_t ax, Float_t ay, Float_t az) const
292 {
293   // returns a value of the correction map
294   return fhCorr->GetBinContent(fhCorr->FindBin(ax, ay, az));
295 }
296
297 //____________________________________________________________________
298 //void AliCorrectionMatrix3D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge, Int_t nBinsZedge)
299 void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
300 {
301   // so what do we do here...
302 }
303
304 //____________________________________________________________________
305 void AliCorrectionMatrix3D::SaveHistograms()
306 {
307   //
308   // saves the histograms
309   //
310
311   AliCorrectionMatrix::SaveHistograms();
312
313   if (GetGeneratedHistogram() && GetMeasuredHistogram())
314   {
315     gDirectory->cd(GetName());
316     
317     AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
318     
319     gDirectory->cd("..");
320   }
321 }
322
323 //____________________________________________________________________
324 Int_t AliCorrectionMatrix3D::CheckEmptyBins(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax, Bool_t quiet)
325 {
326   //
327   // counts the number of empty Bins in a given region
328   //
329
330   TH3* hist = GetGeneratedHistogram();
331   if (!hist)
332     return -1;
333
334   Int_t emptyBins = 0;
335   for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
336     for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
337       for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
338         if (hist->GetBinContent(x, y, z) == 0)
339         {
340           if (!quiet)
341             printf("Empty bin in %s at vtx = %f, eta = %f, pt = %f\n", GetName(), hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
342           ++emptyBins;
343         }
344
345   return emptyBins;
346 }
347
348 //____________________________________________________________________
349 TH1* AliCorrectionMatrix3D::PlotBinErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax)
350 {
351   //
352   // makes a 1d plots of the relative bin errors
353   //
354
355   TH3* hist = GetCorrectionHistogram();
356   if (!hist)
357     return 0;
358
359   TH1F* target = new TH1F("relerrors", "relerrors", 100, 0, 10);
360
361   for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
362     for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
363       for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
364         if (hist->GetBinContent(x, y, z) != 0)
365           target->Fill(100 * hist->GetBinError(x, y, z) / hist->GetBinContent(x, y, z));
366
367   return target;
368 }
369