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