New classes for finding multiple vertices (in case of pile-up). They will be used...
[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//____________________________________________________________________
a6e0ebfe 117TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(const 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()));
51f6de65 165
166 TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
167 divided->Reset();
168
169 divided->Divide(gene1D, meas1D, 1, 1, "B");
bf21645b 170
51f6de65 171 Printf("%p %p", gene1D, meas1D);
bf21645b 172
51f6de65 173 return (TH1F*)divided;
bf21645b 174}
175
176//____________________________________________________________________
177void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
178{
179 // add value to measured histogram
06e4b91b 180 ((TH2F*)fhMeas)->Fill(ax, ay);
bf21645b 181}
182
183//____________________________________________________________________
184void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
185{
186 // add value to generated histogram
06e4b91b 187 ((TH2F*)fhGene)->Fill(ax, ay);
bf21645b 188}
189
190//____________________________________________________________________
191Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
192{
193 // returns a value of the correction map
194 return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
195}
196
197//____________________________________________________________________
198void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
199{
200 // remove edges of correction histogram by removing
201 // - bins with content less than cut
202 // - bins next to bins with zero bin content
203
204 Int_t nBinsX = fhCorr->GetNbinsX();
205 Int_t nBinsY = fhCorr->GetNbinsY();
206
207 // set bin content to zero for bins with content smaller cut
208 for (Int_t bx=0; bx<=nBinsX; bx++) {
209 for (Int_t by=0; by<=nBinsY; by++) {
210 if (fhCorr->GetBinContent(bx,by)>cut) {
211 fhCorr->SetBinContent(bx,by,0);
212 fhCorr->SetBinError(bx,by,0);
213 }
214 }
215 }
216
217 // set bin content to zero for bins next to bins with zero
218 TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
219 tmp->Reset();
220
221 Bool_t done = kFALSE;
222 Int_t nBinsXCount = 0;
223 Int_t nBinsYCount = 0;
224 while (!done) {
225 if (nBinsXCount<nBinsXedge)
226 for (Int_t bx=0; bx<=nBinsX; bx++) {
227 for (Int_t by=0; by<=nBinsY; by++) {
228 if ((fhCorr->GetBinContent(bx+1,by)==0)||
229 (fhCorr->GetBinContent(bx-1,by)==0))
230 tmp->SetBinContent(bx,by,1);
231
232 }
233 }
234 if (nBinsYCount<nBinsYedge)
235 for (Int_t bx=0; bx<=nBinsX; bx++) {
236 for (Int_t by=0; by<=nBinsY; by++) {
237 if ((fhCorr->GetBinContent(bx,by+1)==0)||
238 (fhCorr->GetBinContent(bx,by-1)==0))
239 tmp->SetBinContent(bx,by,1);
240 }
241 }
242 for (Int_t bx=0; bx<=nBinsX; bx++) {
243 for (Int_t by=0; by<=nBinsY; by++) {
244 if (tmp->GetBinContent(bx,by)==1) {
245 fhCorr->SetBinContent(bx,by,0);
246 fhCorr->SetBinError(bx,by,0);
247 }
248 }
249 }
250 nBinsXCount++;
251 nBinsYCount++;
252 if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
253 }
254 tmp->Delete();
255
256}
257
745d6088 258//____________________________________________________________________
259void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
260{
261 // rebins the histograms, recalculates the correction
262
263 GetGeneratedHistogram()->Rebin2D(x, y);
264 GetMeasuredHistogram()->Rebin2D(x, y);
265 GetCorrectionHistogram()->Rebin2D(x, y);
266 Divide();
267}