new way of calculation Nsigma correction
[u/mrichter/AliRoot.git] / PWG0 / CorrectionMatrix2D.cxx
CommitLineData
1caed034 1// ------------------------------------------------------
2//
3// Class to handle 2d-corrections.
4//
5// ------------------------------------------------------
6//
7
dc740de4 8/* $Id$ */
9
ea2ef7cc 10#include <TFile.h>
11#include <TCanvas.h>
12
13#include <AliLog.h>
14
ddec1c19 15#include "CorrectionMatrix2D.h"
16
17//____________________________________________________________________
c56385d9 18ClassImp(CorrectionMatrix2D)
ddec1c19 19
20//____________________________________________________________________
1caed034 21CorrectionMatrix2D::CorrectionMatrix2D(const CorrectionMatrix2D& c)
22 : TNamed(c)
23{
24 // copy constructor
25 ((CorrectionMatrix2D &)c).Copy(*this);
26}
27
28//____________________________________________________________________
b6d3306b 29CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name, Char_t* title,
30 Int_t nBinX, Float_t Xmin, Float_t Xmax,
31 Int_t nBinY, Float_t Ymin, Float_t Ymax)
32 : TNamed(name, title)
ddec1c19 33{
1caed034 34 //
35 // constructor
36 //
37
ddec1c19 38
ea2ef7cc 39 fhMeas = new TH2F(Form("meas_%s",name), Form("meas_%s",title), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
40 fhGene = new TH2F(Form("gene_%s",name), Form("gene_%s",title), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
41 fhCorr = new TH2F(Form("corr_%s",name), Form("corr_%s",title), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
3a282793 42
43 fhMeas->Sumw2();
44 fhGene->Sumw2();
45 fhCorr->Sumw2();
ddec1c19 46}
47
b6d3306b 48//____________________________________________________________________
49CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name,Char_t* title,
50 Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y)
1caed034 51 : TNamed(name, title)
ddec1c19 52{
1caed034 53 //
54 // constructor
55 //
56
ea2ef7cc 57 fhMeas = new TH2F(Form("meas_%s",name), Form("meas_%s",title), nBinX, X, nBinY, Y);
58 fhGene = new TH2F(Form("gene_%s",name), Form("gene_%s",title), nBinX, X, nBinY, Y);
59 fhCorr = new TH2F(Form("corr_%s",name), Form("corr_%s",title), nBinX, X, nBinY, Y);
3a282793 60
61 fhMeas->Sumw2();
62 fhGene->Sumw2();
63 fhCorr->Sumw2();
ddec1c19 64}
65
66
b6d3306b 67//____________________________________________________________________
68CorrectionMatrix2D::~CorrectionMatrix2D() {
1caed034 69 //
70 // destructor
b6d3306b 71 //
72 if (fhMeas) delete fhMeas;
1caed034 73 if (fhGene) delete fhGene;
b6d3306b 74 if (fhCorr) delete fhCorr;
1caed034 75}
76
77//____________________________________________________________________
78CorrectionMatrix2D &CorrectionMatrix2D::operator=(const CorrectionMatrix2D &c)
79{
80 // assigment operator
81
82 if (this != &c)
83 ((CorrectionMatrix2D &) c).Copy(*this);
84
85 return *this;
86}
87
88//____________________________________________________________________
a571621d 89TH1F* CorrectionMatrix2D::Get1DCorrection(Char_t* opt) {
90 //
91 // integrate the correction over one variable
92 //
93
fcf2fb36 94 TH1D* meas1D = 0;
95 TH1D* gene1D = 0;
a571621d 96
97 if (strcmp(opt,"x")==0) {
98 meas1D = fhMeas->ProjectionX();
99 gene1D = fhGene->ProjectionX();
100 }
101 if (strcmp(opt,"y")==0) {
102 meas1D = fhMeas->ProjectionY();
103 gene1D = fhGene->ProjectionY();
104 }
105 gene1D->Sumw2();
106
107 gene1D->SetName(Form("corr_1D_%s",fName.Data()));
108 gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
109
110 gene1D->Divide(gene1D, meas1D, 1, 1, "B");
111
112 return (TH1F*)gene1D;
113}
114
115
116//____________________________________________________________________
1caed034 117void
118CorrectionMatrix2D::Copy(TObject& c) const
119{
120 // copy function
121
122 CorrectionMatrix2D& target = (CorrectionMatrix2D &) c;
123
124 target.fhMeas = fhMeas;
125 target.fhGene = fhGene;
126 target.fhCorr = fhCorr;
b6d3306b 127}
ddec1c19 128
ddec1c19 129
b6d3306b 130//________________________________________________________________________
131void CorrectionMatrix2D::SetAxisTitles(Char_t* titleX, Char_t* titleY)
132{
1caed034 133 //
134 // method for setting the axis titles of the histograms
135 //
136
b6d3306b 137 fhMeas ->SetXTitle(titleX); fhMeas ->SetYTitle(titleY);
138 fhGene ->SetXTitle(titleX); fhGene ->SetYTitle(titleY);
139 fhCorr ->SetXTitle(titleX); fhCorr ->SetYTitle(titleY);
ea2ef7cc 140}
141
142//____________________________________________________________________
143Long64_t CorrectionMatrix2D::Merge(TCollection* list) {
144 // Merge a list of CorrectionMatrix2D objects with this (needed for
145 // PROOF).
146 // Returns the number of merged objects (including this).
147
148 if (!list)
149 return 0;
150
151 if (list->IsEmpty())
152 return 1;
153
154 TIterator* iter = list->MakeIterator();
155 TObject* obj;
156
157 // collections of measured and generated histograms
158 TList* collectionMeas = new TList;
159 TList* collectionGene = new TList;
160
161 Int_t count = 0;
162 while ((obj = iter->Next())) {
163
164 CorrectionMatrix2D* entry = dynamic_cast<CorrectionMatrix2D*> (obj);
165 if (entry == 0)
166 continue;
167
168 collectionMeas->Add(entry->GetMeasuredHistogram());
169 collectionGene->Add(entry->GetGeneratedHistogram());
170
171 count++;
172 }
173 fhMeas->Merge(collectionMeas);
174 fhGene->Merge(collectionGene);
175
176 // is this really faster than just adding the histograms in the list???
177 delete collectionMeas;
178 delete collectionGene;
179
180
181 return count+1;
ddec1c19 182}
183
184
185//____________________________________________________________________
ea2ef7cc 186void CorrectionMatrix2D::Divide() {
1caed034 187 //
ea2ef7cc 188 // divide the histograms to get the correction
1caed034 189 //
ddec1c19 190
1caed034 191 if (!fhMeas || !fhGene) return;
ddec1c19 192
b6d3306b 193 fhCorr->Divide(fhGene, fhMeas, 1,1,"B");
ddec1c19 194
195}
196
197//____________________________________________________________________
198void
1caed034 199CorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
200{
ddec1c19 201 // remove edges of correction histogram by removing
202 // - bins with content less than cut
203 // - bins next to bins with zero bin content
204
b6d3306b 205 Int_t nBinsX = fhCorr->GetNbinsX();
206 Int_t nBinsY = fhCorr->GetNbinsY();
ddec1c19 207
208 // set bin content to zero for bins with content smaller cut
209 for (Int_t bx=0; bx<=nBinsX; bx++) {
210 for (Int_t by=0; by<=nBinsY; by++) {
b6d3306b 211 if (fhCorr->GetBinContent(bx,by)>cut) {
212 fhCorr->SetBinContent(bx,by,0);
213 fhCorr->SetBinError(bx,by,0);
ddec1c19 214 }
215 }
216 }
217
218 // set bin content to zero for bins next to bins with zero
b6d3306b 219 TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
ddec1c19 220 tmp->Reset();
221
222 Bool_t done = kFALSE;
223 Int_t nBinsXCount = 0;
224 Int_t nBinsYCount = 0;
225 while (!done) {
226 if (nBinsXCount<nBinsXedge)
227 for (Int_t bx=0; bx<=nBinsX; bx++) {
228 for (Int_t by=0; by<=nBinsY; by++) {
b6d3306b 229 if ((fhCorr->GetBinContent(bx+1,by)==0)||
230 (fhCorr->GetBinContent(bx-1,by)==0))
ddec1c19 231 tmp->SetBinContent(bx,by,1);
232
233 }
234 }
235 if (nBinsYCount<nBinsYedge)
236 for (Int_t bx=0; bx<=nBinsX; bx++) {
237 for (Int_t by=0; by<=nBinsY; by++) {
b6d3306b 238 if ((fhCorr->GetBinContent(bx,by+1)==0)||
239 (fhCorr->GetBinContent(bx,by-1)==0))
ddec1c19 240 tmp->SetBinContent(bx,by,1);
241 }
242 }
243 for (Int_t bx=0; bx<=nBinsX; bx++) {
244 for (Int_t by=0; by<=nBinsY; by++) {
245 if (tmp->GetBinContent(bx,by)==1) {
b6d3306b 246 fhCorr->SetBinContent(bx,by,0);
247 fhCorr->SetBinError(bx,by,0);
ddec1c19 248 }
249 }
250 }
251 nBinsXCount++;
252 nBinsYCount++;
253 if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
254 }
255 tmp->Delete();
256
257}
258
259//____________________________________________________________________
260Bool_t CorrectionMatrix2D::LoadHistograms(Char_t* fileName, Char_t* dir) {
1caed034 261 //
262 // loads the histograms from a file
263 //
264
ddec1c19 265 TFile* fin = TFile::Open(fileName);
266
1caed034 267 if(!fin) {
ddec1c19 268 //Info("LoadHistograms",Form(" %s file does not exist",fileName));
1caed034 269 return kFALSE;
ddec1c19 270 }
271
b6d3306b 272 if(fhGene) {delete fhGene; fhGene=0;}
273 if(fhCorr) {delete fhCorr; fhCorr=0;}
b6d3306b 274 if(fhMeas) {delete fhMeas; fhMeas=0;}
ddec1c19 275
f2cb0edf 276 fhMeas = (TH2F*)fin->Get(Form("%s/meas_%s", dir,GetName()));
277 if(!fhMeas) Info("LoadHistograms","No meas hist available");
278
279 fhGene = (TH2F*)fin->Get(Form("%s/gene_%s",dir, GetName()));
280 if(!fhGene) Info("LoadHistograms","No gene hist available");
281
282 fhCorr = (TH2F*)fin->Get(Form("%s/corr_%s",dir, GetName()));
283 if(!fhCorr)
284 {
285 Info("LoadHistograms","No corr hist available");
286 return kFALSE;
287 }
ddec1c19 288
289 return kTRUE;
290}
291
292
293//____________________________________________________________________
294void
295CorrectionMatrix2D::SaveHistograms() {
1caed034 296 //
297 // saves the histograms
298 //
a571621d 299
b6d3306b 300 fhMeas ->Write();
301 fhGene ->Write();
ddec1c19 302
b6d3306b 303 if (fhCorr)
304 fhCorr->Write();
ddec1c19 305}
306
307//____________________________________________________________________
308void CorrectionMatrix2D::DrawHistograms()
309{
1caed034 310 //
311 // draws all the four histograms on one TCanvas
312 //
313
a571621d 314 TCanvas* canvas = new TCanvas(Form("correction_%s",fName.Data()),
315 Form("correction_%s",fName.Data()), 800, 800);
ddec1c19 316 canvas->Divide(2, 2);
a571621d 317
ddec1c19 318 canvas->cd(1);
b6d3306b 319 if (fhMeas)
320 fhMeas->Draw("COLZ");
a571621d 321
ddec1c19 322 canvas->cd(2);
b6d3306b 323 if (fhGene)
324 fhGene->Draw("COLZ");
ddec1c19 325
326 canvas->cd(3);
b6d3306b 327 if (fhCorr)
328 fhCorr->Draw("COLZ");
ea2ef7cc 329
a571621d 330 canvas->cd(4);
331
ea2ef7cc 332 // add: draw here the stat. errors of the correction histogram
333
ddec1c19 334}
335
336