Update for Ds
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix2D.cxx
CommitLineData
bf21645b 1/* $Id$ */
2
3// ------------------------------------------------------
4//
5// Class to handle 2d-corrections.
6//
7// ------------------------------------------------------
8//
9
10#include <TH2F.h>
69b09e3b 11#include <TMath.h>
bf21645b 12
13#include <AliLog.h>
14
15#include "AliCorrectionMatrix2D.h"
16
17//____________________________________________________________________
18ClassImp(AliCorrectionMatrix2D)
19
20//____________________________________________________________________
21AliCorrectionMatrix2D::AliCorrectionMatrix2D() :
22 AliCorrectionMatrix()
23{
24 // default constructor
25}
26
27//____________________________________________________________________
28AliCorrectionMatrix2D::AliCorrectionMatrix2D(const AliCorrectionMatrix2D& c)
29 : AliCorrectionMatrix(c)
30{
31 // copy constructor
32 ((AliCorrectionMatrix2D &)c).Copy(*this);
33}
34
35//____________________________________________________________________
61385583 36AliCorrectionMatrix2D &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//____________________________________________________________________
bf21645b 47AliCorrectionMatrix2D::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
29771dc8 56 // do not add this hists to the directory
57 Bool_t oldStatus = TH1::AddDirectoryStatus();
58 TH1::AddDirectory(kFALSE);
59
82ac9fc4 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);
29771dc8 63
64 TH1::AddDirectory(oldStatus);
bf21645b 65
66 fhMeas->Sumw2();
67 fhGene->Sumw2();
68 fhCorr->Sumw2();
69}
70
71//____________________________________________________________________
72AliCorrectionMatrix2D::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
29771dc8 80 // do not add this hists to the directory
81 Bool_t oldStatus = TH1::AddDirectoryStatus();
82 TH1::AddDirectory(kFALSE);
bf21645b 83
29771dc8 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();
bf21645b 91 fhGene->Sumw2();
92 fhCorr->Sumw2();
93}
94
95//____________________________________________________________________
96AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
97{
98 //
99 // destructor
100 //
101
102 // histograms already deleted in base class
103}
104
5a6310fe 105TH2* AliCorrectionMatrix2D::GetGeneratedHistogram() const
bf21645b 106{
107 // return generated histogram casted to correct type
5a6310fe 108 return dynamic_cast<TH2*> (fhGene);
bf21645b 109}
110
5a6310fe 111TH2* AliCorrectionMatrix2D::GetMeasuredHistogram() const
bf21645b 112{
113 // return measured histogram casted to correct type
5a6310fe 114 return dynamic_cast<TH2*> (fhMeas);
bf21645b 115}
116
117//____________________________________________________________________
69b09e3b 118TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(const Char_t* opt, Float_t min, Float_t max, Bool_t binomialErrors)
bf21645b 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) {
06e4b91b 128 Int_t binMin = fhMeas->GetYaxis()->FindBin(min);
129 Int_t binMax = fhGene->GetYaxis()->FindBin(max);
7584d357 130
06e4b91b 131 if (min>=max) {
132 meas1D = ((TH2F*)fhMeas)->ProjectionX();
133 gene1D = ((TH2F*)fhGene)->ProjectionX();
7584d357 134 }
135 else {
dd367a14 136 Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);
7584d357 137
5a6310fe 138 meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_x_pm", GetName()),binMin,binMax);
139 gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_x_pg", GetName()),binMin,binMax);
7584d357 140 }
bf21645b 141 }
5a6310fe 142 else if (strcmp(opt,"y")==0) {
06e4b91b 143 Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
144 Int_t binMax = fhMeas->GetXaxis()->FindBin(max);
dd367a14 145
06e4b91b 146 if (min>=max) {
147 meas1D = ((TH2F*)fhMeas)->ProjectionY();
148 gene1D = ((TH2F*)fhGene)->ProjectionY();
7584d357 149 }
150 else {
5a6310fe 151 Printf("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax);
dd367a14 152
5a6310fe 153 meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_y_pm", GetName()), binMin, binMax);
154 gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_y_pg", GetName()), binMin, binMax);
7584d357 155 }
bf21645b 156 }
5a6310fe 157 else {
158 Printf("ERROR: Invalid option");
159 return 0;
160 }
161
69b09e3b 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
bf21645b 173 gene1D->SetName(Form("corr_1D_%s",fName.Data()));
174 gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
51f6de65 175
176 TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
177 divided->Reset();
69b09e3b 178
179 divided->Divide(gene1D, meas1D, 1, 1, (binomialErrors) ? "B" : "");
bf21645b 180
51f6de65 181 Printf("%p %p", gene1D, meas1D);
bf21645b 182
51f6de65 183 return (TH1F*)divided;
bf21645b 184}
185
186//____________________________________________________________________
187void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
188{
189 // add value to measured histogram
06e4b91b 190 ((TH2F*)fhMeas)->Fill(ax, ay);
bf21645b 191}
192
193//____________________________________________________________________
194void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
195{
196 // add value to generated histogram
06e4b91b 197 ((TH2F*)fhGene)->Fill(ax, ay);
bf21645b 198}
199
200//____________________________________________________________________
201Float_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//____________________________________________________________________
208void 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
745d6088 268//____________________________________________________________________
269void 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}