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