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