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