Update for Ds
[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 #include <TMath.h>
16
17 #include <AliLog.h>
18
19 #include "AliCorrectionMatrix3D.h"
20 #include "AliCorrectionMatrix2D.h"
21 #include "AliPWG0Helper.h"
22
23 //____________________________________________________________________
24 ClassImp(AliCorrectionMatrix3D)
25
26 //____________________________________________________________________
27 AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
28   AliCorrectionMatrix()
29 {
30   // default constructor
31 }
32
33 //____________________________________________________________________
34 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
35   : AliCorrectionMatrix(c)
36 {
37   // copy constructor
38   ((AliCorrectionMatrix3D &)c).Copy(*this);
39 }
40
41 //____________________________________________________________________
42 AliCorrectionMatrix3D &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 //____________________________________________________________________
53 AliCorrectionMatrix3D::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
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;
66
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;
70
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;
74
75   CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
76
77   delete[] binLimitsX;
78   delete[] binLimitsY;
79   delete[] binLimitsZ;
80 }
81
82 AliCorrectionMatrix3D::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
98   CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
99
100   delete[] binLimitsX;
101   delete[] binLimitsY;
102 }
103
104 AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3* hBinning)
105   : AliCorrectionMatrix(name, title)
106 {
107   // constructor with variable bin sizes (uses binning of hBinning)
108
109   // do not add this hists to the directory
110   Bool_t oldStatus = TH1::AddDirectoryStatus();
111   TH1::AddDirectory(kFALSE);
112
113   fhMeas = (TH3F*)hBinning->Clone("measured");
114   fhGene = (TH3F*)hBinning->Clone("generated");
115   fhCorr = (TH3F*)hBinning->Clone("correction");
116
117   fhMeas->SetTitle(Form("%s measured", GetTitle()));
118   fhGene->SetTitle(Form("%s generated", GetTitle()));
119   fhCorr->SetTitle(Form("%s correction", GetTitle()));
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();
130 }
131
132 //____________________________________________________________________
133 void 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);
146
147   fhMeas->Sumw2();
148   fhGene->Sumw2();
149   fhCorr->Sumw2();
150
151   TH1::AddDirectory(oldStatus);
152 }
153
154
155 //____________________________________________________________________
156 AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
157 {
158   //
159   // destructor
160   //
161
162   // histograms already deleted in base class
163 }
164
165 //____________________________________________________________________
166 TH3* AliCorrectionMatrix3D::GetGeneratedHistogram()
167 {
168   // return generated histogram casted to correct type
169   return dynamic_cast<TH3F*> (fhGene);
170 }
171
172 //____________________________________________________________________
173 TH3* AliCorrectionMatrix3D::GetMeasuredHistogram()
174 {
175   // return measured histogram casted to correct type
176   return dynamic_cast<TH3F*> (fhMeas);
177 }
178
179 //____________________________________________________________________
180 TH3* AliCorrectionMatrix3D::GetCorrectionHistogram()
181 {
182   // return correction histogram casted to correct type
183   return dynamic_cast<TH3F*> (fhCorr);
184 }
185
186 //____________________________________________________________________
187 AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Option_t* opt, Float_t aMin, Float_t aMax)
188 {
189   // returns a 2D projection of this correction
190
191   TString option = opt;
192
193   // unzoom
194   fhMeas->GetXaxis()->SetRange(0, 0);
195   fhMeas->GetYaxis()->SetRange(0, 0);
196   fhMeas->GetZaxis()->SetRange(0, 0);
197
198   fhGene->GetXaxis()->SetRange(0, 0);
199   fhGene->GetYaxis()->SetRange(0, 0);
200   fhGene->GetZaxis()->SetRange(0, 0);
201
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);
206       fhGene->GetZaxis()->SetRange(bMin, bMax);
207       fhMeas->GetZaxis()->SetRange(bMin, bMax);
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);
212       fhGene->GetYaxis()->SetRange(bMin, bMax);
213       fhMeas->GetYaxis()->SetRange(bMin, bMax);
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);
218       fhGene->GetXaxis()->SetRange(bMin, bMax);
219       fhMeas->GetXaxis()->SetRange(bMin, bMax);
220     }
221     else {
222       AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
223       return 0;
224     }
225   }
226
227   AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
228
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()));
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     }
239
240   TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
241   corr->Reset();
242
243   corr2D->SetGeneratedHistogram(gene);
244   corr2D->SetMeasuredHistogram(meas);
245   corr2D->SetCorrectionHistogram(corr);
246
247   corr2D->Divide();
248
249   // unzoom
250   fhMeas->GetXaxis()->SetRange(0, 0);
251   fhMeas->GetYaxis()->SetRange(0, 0);
252   fhMeas->GetZaxis()->SetRange(0, 0);
253
254   fhGene->GetXaxis()->SetRange(0, 0);
255   fhGene->GetYaxis()->SetRange(0, 0);
256   fhGene->GetZaxis()->SetRange(0, 0);
257
258   return corr2D;
259 }
260
261 //____________________________________________________________________
262 TH1* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Option_t* opt, Float_t aMin1, Float_t aMax1, Float_t aMin2, Float_t aMax2)
263 {
264   // returns a 1D projection of this correction
265   
266   AliCorrectionMatrix2D* corr2D = 0;
267   if (strcmp(opt,"x")==0) {
268     corr2D = Get2DCorrection("yx",aMin2,aMax2);
269     return corr2D->Get1DCorrectionHistogram("x",aMin1,aMax1);
270   }
271   if (strcmp(opt,"y")==0) {
272     corr2D = Get2DCorrection("xy",aMin2,aMax2);
273     return corr2D->Get1DCorrectionHistogram("x",aMin1,aMax1);
274   }  
275   if (strcmp(opt,"z")==0) {
276     corr2D = Get2DCorrection("yz",aMin1,aMax1);
277     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
278   }  
279   AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));
280   
281   return 0;
282 }
283
284 //____________________________________________________________________
285 void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az, Double_t weight)
286 {
287   // add value to measured histogram
288   GetMeasuredHistogram()->Fill(ax, ay, az, weight);
289 }
290
291 //____________________________________________________________________
292 void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az, Double_t weight)
293 {
294   // add value to generated histogram
295   GetGeneratedHistogram()->Fill(ax, ay, az, weight);
296 }
297
298 //____________________________________________________________________
299 Float_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)
307 void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
308 {
309   // so what do we do here...
310 }
311
312 //____________________________________________________________________
313 void AliCorrectionMatrix3D::SaveHistograms()
314 {
315   //
316   // saves the histograms
317   //
318
319   AliCorrectionMatrix::SaveHistograms();
320
321   if (GetGeneratedHistogram() && GetMeasuredHistogram())
322   {
323     gDirectory->cd(GetName());
324     
325     AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
326     
327     gDirectory->cd("..");
328   }
329 }
330
331 //____________________________________________________________________
332 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)
333 {
334   //
335   // counts the number of empty Bins in a given region
336   //
337
338   TH3* hist = GetGeneratedHistogram();
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 //____________________________________________________________________
357 TH1* AliCorrectionMatrix3D::PlotBinErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax)
358 {
359   //
360   // makes a 1d plots of the relative bin errors
361   //
362
363   TH3* hist = GetCorrectionHistogram();
364   if (!hist)
365     return 0;
366
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