Coverity fixes
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix2D.cxx
1 /* $Id$ */
2
3 // ------------------------------------------------------
4 //
5 // Class to handle 2d-corrections.
6 //
7 // ------------------------------------------------------
8 //
9
10 #include <TH2F.h>
11 #include <TMath.h>
12
13 #include <AliLog.h>
14
15 #include "AliCorrectionMatrix2D.h"
16
17 //____________________________________________________________________
18 ClassImp(AliCorrectionMatrix2D)
19
20 //____________________________________________________________________
21 AliCorrectionMatrix2D::AliCorrectionMatrix2D() :
22   AliCorrectionMatrix()
23 {
24   // default constructor
25 }
26
27 //____________________________________________________________________
28 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const AliCorrectionMatrix2D& c)
29   : AliCorrectionMatrix(c)
30 {
31   // copy constructor
32   ((AliCorrectionMatrix2D &)c).Copy(*this);
33 }
34
35 //____________________________________________________________________
36 AliCorrectionMatrix2D &AliCorrectionMatrix2D::operator=(const AliCorrectionMatrix2D &c)
37 {
38   // assigment operator
39
40   if (this != &c)
41     ((AliCorrectionMatrix2D &) c).Copy(*this);
42
43   return *this;
44 }
45
46 //____________________________________________________________________
47 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
48                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
49                                        Int_t nBinY, Float_t Ymin, Float_t Ymax) 
50   : AliCorrectionMatrix(name, title)
51 {
52   //
53   // constructor
54   //
55
56   // do not add this hists to the directory
57   Bool_t oldStatus = TH1::AddDirectoryStatus();
58   TH1::AddDirectory(kFALSE);
59
60   fhMeas  = new TH2F("measured",   Form("%s measured", GetTitle()),   nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
61   fhGene  = new TH2F("generated",  Form("%s generated", GetTitle()),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
62   fhCorr  = new TH2F("correction", Form("%s correction", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
63
64   TH1::AddDirectory(oldStatus);
65
66   fhMeas->Sumw2();
67   fhGene->Sumw2();
68   fhCorr->Sumw2();
69 }
70
71 //____________________________________________________________________
72 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
73                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
74   : AliCorrectionMatrix(name, title)
75 {
76   //
77   // constructor
78   //
79
80   // do not add this hists to the directory
81   Bool_t oldStatus = TH1::AddDirectoryStatus();
82   TH1::AddDirectory(kFALSE);
83
84         fhMeas  = new TH2F("measured",   Form("%s measured",title),   nBinX, X, nBinY, Y);
85   fhGene  = new TH2F("generated",  Form("%s generated",title),  nBinX, X, nBinY, Y);
86   fhCorr  = new TH2F("correction", Form("%s correction",title), nBinX, X, nBinY, Y);
87
88   TH1::AddDirectory(oldStatus);
89
90         fhMeas->Sumw2();
91   fhGene->Sumw2();
92   fhCorr->Sumw2();
93 }
94
95 //____________________________________________________________________
96 AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
97 {
98   //
99   // destructor
100   //
101
102   // histograms already deleted in base class
103 }
104
105 TH2* AliCorrectionMatrix2D::GetGeneratedHistogram() const
106 {
107   // return generated histogram casted to correct type
108   return dynamic_cast<TH2*> (fhGene);
109 }
110
111 TH2* AliCorrectionMatrix2D::GetMeasuredHistogram() const
112 {
113   // return measured histogram casted to correct type
114   return dynamic_cast<TH2*> (fhMeas);
115 }
116
117 //____________________________________________________________________
118 TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(const Char_t* opt, Float_t min, Float_t max, Bool_t binomialErrors)
119 {
120   //
121   // integrate the correction over one variable 
122   // 
123
124   TH1D* meas1D = 0;
125   TH1D* gene1D = 0; 
126
127   if (strcmp(opt,"x")==0) {
128     Int_t binMin = fhMeas->GetYaxis()->FindBin(min);
129     Int_t binMax = fhGene->GetYaxis()->FindBin(max);
130
131     if (min>=max) {
132       meas1D = ((TH2F*)fhMeas)->ProjectionX();
133       gene1D = ((TH2F*)fhGene)->ProjectionX();
134     }
135     else {
136       Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);
137
138       meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_x_pm", GetName()),binMin,binMax);
139       gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_x_pg", GetName()),binMin,binMax);
140     }
141   }
142   else if (strcmp(opt,"y")==0) {
143     Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
144     Int_t binMax = fhMeas->GetXaxis()->FindBin(max);
145
146     if (min>=max) {
147       meas1D = ((TH2F*)fhMeas)->ProjectionY();
148       gene1D = ((TH2F*)fhGene)->ProjectionY();
149     }
150     else {
151       Printf("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax);
152
153       meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_y_pm", GetName()), binMin, binMax);
154       gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_y_pg", GetName()), binMin, binMax);
155     }
156   }
157   else {
158     Printf("ERROR: Invalid option");
159     return 0;
160   }
161
162   if (!binomialErrors)
163   {
164     // set the errors on gene manually, and clear the ones on meas.
165     gene1D->Sumw2();
166     for (Int_t bin=0; bin <= gene1D->GetNbinsX()+1; bin++)
167     {
168       gene1D->SetBinError(bin, TMath::Sqrt(gene1D->GetBinContent(bin)));
169       meas1D->SetBinError(bin, 0);
170     }
171   }
172   
173   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
174   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
175  
176   TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
177   divided->Reset();
178   
179   divided->Divide(gene1D, meas1D, 1, 1, (binomialErrors) ? "B" : "");
180
181   Printf("%p %p", gene1D, meas1D);
182   
183   return (TH1F*)divided;   
184 }
185
186 //____________________________________________________________________
187 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
188 {
189   // add value to measured histogram
190   ((TH2F*)fhMeas)->Fill(ax, ay);
191 }
192
193 //____________________________________________________________________
194 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
195 {
196   // add value to generated histogram
197   ((TH2F*)fhGene)->Fill(ax, ay);
198 }
199
200 //____________________________________________________________________
201 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
202 {
203   // returns a value of the correction map
204   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
205 }
206
207 //____________________________________________________________________
208 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
209 {
210   // remove edges of correction histogram by removing 
211   // - bins with content less than cut
212   // - bins next to bins with zero bin content
213   
214   Int_t nBinsX = fhCorr->GetNbinsX();
215   Int_t nBinsY = fhCorr->GetNbinsY();
216
217   // set bin content to zero for bins with content smaller cut
218   for (Int_t bx=0; bx<=nBinsX; bx++) {
219     for (Int_t by=0; by<=nBinsY; by++) {
220       if (fhCorr->GetBinContent(bx,by)>cut) {
221           fhCorr->SetBinContent(bx,by,0);
222           fhCorr->SetBinError(bx,by,0);
223       }
224     }
225   }
226
227   // set bin content to zero for bins next to bins with zero
228   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
229   tmp->Reset();
230   
231   Bool_t done = kFALSE;
232   Int_t nBinsXCount = 0;
233   Int_t nBinsYCount = 0;
234   while (!done) {    
235     if (nBinsXCount<nBinsXedge) 
236       for (Int_t bx=0; bx<=nBinsX; bx++) {
237         for (Int_t by=0; by<=nBinsY; by++) {
238           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
239               (fhCorr->GetBinContent(bx-1,by)==0))
240             tmp->SetBinContent(bx,by,1);        
241           
242         }
243       }
244     if (nBinsYCount<nBinsYedge) 
245       for (Int_t bx=0; bx<=nBinsX; bx++) {
246         for (Int_t by=0; by<=nBinsY; by++) {
247           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
248               (fhCorr->GetBinContent(bx,by-1)==0))
249             tmp->SetBinContent(bx,by,1);        
250         }
251       }    
252     for (Int_t bx=0; bx<=nBinsX; bx++) {
253       for (Int_t by=0; by<=nBinsY; by++) {
254         if (tmp->GetBinContent(bx,by)==1) {
255           fhCorr->SetBinContent(bx,by,0);
256           fhCorr->SetBinError(bx,by,0);
257         }
258       }
259     }
260     nBinsXCount++;
261     nBinsYCount++;
262     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
263   }
264   tmp->Delete();  
265
266 }
267
268 //____________________________________________________________________
269 void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
270 {
271         // rebins the histograms, recalculates the correction
272
273         GetGeneratedHistogram()->Rebin2D(x, y);
274         GetMeasuredHistogram()->Rebin2D(x, y);
275         GetCorrectionHistogram()->Rebin2D(x, y);
276         Divide();
277 }