bug fixed that was considered also under/overflow bins in vertex axis
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix3D.cxx
CommitLineData
bf21645b 1/* $Id$ */
2
3// ------------------------------------------------------
4//
5// Class to handle 3d-corrections.
6//
7// ------------------------------------------------------
8//
9
d51554df 10#include <TDirectory.h>
f10a1859 11#include <TH1F.h>
d51554df 12#include <TH2F.h>
13#include <TH3F.h>
06e4b91b 14#include <TString.h>
bf21645b 15
16#include <AliLog.h>
17
18#include "AliCorrectionMatrix3D.h"
06e4b91b 19#include "AliCorrectionMatrix2D.h"
847489f7 20#include "AliPWG0Helper.h"
bf21645b 21
22//____________________________________________________________________
23ClassImp(AliCorrectionMatrix3D)
24
25//____________________________________________________________________
26AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
27 AliCorrectionMatrix()
28{
29 // default constructor
30}
31
32//____________________________________________________________________
33AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
34 : AliCorrectionMatrix(c)
35{
36 // copy constructor
37 ((AliCorrectionMatrix3D &)c).Copy(*this);
38}
39
40//____________________________________________________________________
61385583 41AliCorrectionMatrix3D &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//____________________________________________________________________
bf21645b 52AliCorrectionMatrix3D::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
083a636e 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;
29771dc8 65
083a636e 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;
29771dc8 69
083a636e 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;
bf21645b 73
083a636e 74 CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
75
76 delete[] binLimitsX;
77 delete[] binLimitsY;
78 delete[] binLimitsZ;
bf21645b 79}
80
1afae8ff 81AliCorrectionMatrix3D::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
083a636e 97 CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
29771dc8 98
083a636e 99 delete[] binLimitsX;
100 delete[] binLimitsY;
101}
29771dc8 102
06e4b91b 103AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3F* hBinning)
083a636e 104 : AliCorrectionMatrix(name, title)
105{
06e4b91b 106 // constructor with variable bin sizes (uses binning of hBinning)
083a636e 107
06e4b91b 108 // do not add this hists to the directory
109 Bool_t oldStatus = TH1::AddDirectoryStatus();
110 TH1::AddDirectory(kFALSE);
083a636e 111
06e4b91b 112 fhMeas = (TH3F*)hBinning->Clone("measured");
113 fhGene = (TH3F*)hBinning->Clone("generated");
114 fhCorr = (TH3F*)hBinning->Clone("correction");
1afae8ff 115
82ac9fc4 116 fhMeas->SetTitle(Form("%s measured", GetTitle()));
117 fhGene->SetTitle(Form("%s generated", GetTitle()));
118 fhCorr->SetTitle(Form("%s correction", GetTitle()));
06e4b91b 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();
083a636e 129}
130
131//____________________________________________________________________
132void 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);
1afae8ff 145
146 fhMeas->Sumw2();
147 fhGene->Sumw2();
148 fhCorr->Sumw2();
083a636e 149
150 TH1::AddDirectory(oldStatus);
1afae8ff 151}
152
083a636e 153
bf21645b 154//____________________________________________________________________
155AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
156{
157 //
158 // destructor
159 //
160
161 // histograms already deleted in base class
162}
163
164//____________________________________________________________________
165TH3F* AliCorrectionMatrix3D::GetGeneratedHistogram()
166{
167 // return generated histogram casted to correct type
168 return dynamic_cast<TH3F*> (fhGene);
169}
170
171//____________________________________________________________________
172TH3F* AliCorrectionMatrix3D::GetMeasuredHistogram()
173{
174 // return measured histogram casted to correct type
175 return dynamic_cast<TH3F*> (fhMeas);
176}
177
178//____________________________________________________________________
179TH3F* AliCorrectionMatrix3D::GetCorrectionHistogram()
180{
181 // return correction histogram casted to correct type
182 return dynamic_cast<TH3F*> (fhCorr);
183}
184
185//____________________________________________________________________
06e4b91b 186AliCorrectionMatrix2D* 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
dd367a14 192 // unzoom
193 fhMeas->GetXaxis()->UnZoom();
194 fhMeas->GetYaxis()->UnZoom();
195 fhMeas->GetZaxis()->UnZoom();
196
197 fhGene->GetXaxis()->UnZoom();
198 fhGene->GetYaxis()->UnZoom();
199 fhGene->GetZaxis()->UnZoom();
200
06e4b91b 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);
dd367a14 205 fhGene->GetZaxis()->SetRange(bMin, bMax);
206 fhMeas->GetZaxis()->SetRange(bMin, bMax);
06e4b91b 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);
dd367a14 211 fhGene->GetYaxis()->SetRange(bMin, bMax);
212 fhMeas->GetYaxis()->SetRange(bMin, bMax);
06e4b91b 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);
dd367a14 217 fhGene->GetXaxis()->SetRange(bMin, bMax);
218 fhMeas->GetXaxis()->SetRange(bMin, bMax);
06e4b91b 219 }
220 else {
221 AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
222 return 0;
223 }
224 }
dd367a14 225
06e4b91b 226 AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
227
dd367a14 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()));
06e4b91b 230
dd367a14 231 TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
06e4b91b 232 corr->Reset();
233
234 corr2D->SetGeneratedHistogram(gene);
235 corr2D->SetMeasuredHistogram(meas);
236 corr2D->SetCorrectionHistogram(corr);
237
238 corr2D->Divide();
239
240 // unzoom
dd367a14 241 fhMeas->GetXaxis()->UnZoom();
242 fhMeas->GetYaxis()->UnZoom();
243 fhMeas->GetZaxis()->UnZoom();
06e4b91b 244
245 fhGene->GetXaxis()->UnZoom();
246 fhGene->GetYaxis()->UnZoom();
247 fhGene->GetZaxis()->UnZoom();
248
249 return corr2D;
250}
251
252//____________________________________________________________________
253TH1F* 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;
dd367a14 259 if (strcmp(opt,"x")==0) {
260 corr2D = Get2DCorrection("yx",aMin1,aMax1);
06e4b91b 261 return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
262 }
263 if (strcmp(opt,"y")==0) {
264 corr2D = Get2DCorrection("xy",aMin1,aMax1);
dd367a14 265 return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
06e4b91b 266 }
267 if (strcmp(opt,"z")==0) {
dd367a14 268 corr2D = Get2DCorrection("yz",aMin1,aMax1);
06e4b91b 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
06e4b91b 276//____________________________________________________________________
bf21645b 277void 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//____________________________________________________________________
284void 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//____________________________________________________________________
291Float_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)
299void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
300{
301 // so what do we do here...
302}
303
304//____________________________________________________________________
305void AliCorrectionMatrix3D::SaveHistograms()
306{
307 //
308 // saves the histograms
309 //
310
311 AliCorrectionMatrix::SaveHistograms();
312
92d2d8ad 313 if (GetGeneratedHistogram() && GetMeasuredHistogram())
29771dc8 314 {
315 gDirectory->cd(GetName());
316
317 AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
318
319 gDirectory->cd("..");
320 }
bf21645b 321}
f10a1859 322
323//____________________________________________________________________
324Int_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 TH3F* 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//____________________________________________________________________
349TH1F* 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 TH3F* 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