writing empty Digits too
[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//____________________________________________________________________
81void
82CorrectionMatrix2D::Copy(TObject& c) const
83{
84 // copy function
85
86 CorrectionMatrix2D& target = (CorrectionMatrix2D &) c;
87
88 target.fhMeas = fhMeas;
89 target.fhGene = fhGene;
90 target.fhCorr = fhCorr;
b6d3306b 91}
ddec1c19 92
ddec1c19 93
b6d3306b 94//________________________________________________________________________
95void CorrectionMatrix2D::SetAxisTitles(Char_t* titleX, Char_t* titleY)
96{
1caed034 97 //
98 // method for setting the axis titles of the histograms
99 //
100
b6d3306b 101 fhMeas ->SetXTitle(titleX); fhMeas ->SetYTitle(titleY);
102 fhGene ->SetXTitle(titleX); fhGene ->SetYTitle(titleY);
103 fhCorr ->SetXTitle(titleX); fhCorr ->SetYTitle(titleY);
ea2ef7cc 104}
105
106//____________________________________________________________________
107Long64_t CorrectionMatrix2D::Merge(TCollection* list) {
108 // Merge a list of CorrectionMatrix2D objects with this (needed for
109 // PROOF).
110 // Returns the number of merged objects (including this).
111
112 if (!list)
113 return 0;
114
115 if (list->IsEmpty())
116 return 1;
117
118 TIterator* iter = list->MakeIterator();
119 TObject* obj;
120
121 // collections of measured and generated histograms
122 TList* collectionMeas = new TList;
123 TList* collectionGene = new TList;
124
125 Int_t count = 0;
126 while ((obj = iter->Next())) {
127
128 CorrectionMatrix2D* entry = dynamic_cast<CorrectionMatrix2D*> (obj);
129 if (entry == 0)
130 continue;
131
132 collectionMeas->Add(entry->GetMeasuredHistogram());
133 collectionGene->Add(entry->GetGeneratedHistogram());
134
135 count++;
136 }
137 fhMeas->Merge(collectionMeas);
138 fhGene->Merge(collectionGene);
139
140 // is this really faster than just adding the histograms in the list???
141 delete collectionMeas;
142 delete collectionGene;
143
144
145 return count+1;
ddec1c19 146}
147
148
149//____________________________________________________________________
ea2ef7cc 150void CorrectionMatrix2D::Divide() {
1caed034 151 //
ea2ef7cc 152 // divide the histograms to get the correction
1caed034 153 //
ddec1c19 154
1caed034 155 if (!fhMeas || !fhGene) return;
ddec1c19 156
b6d3306b 157 fhCorr->Divide(fhGene, fhMeas, 1,1,"B");
ddec1c19 158
159}
160
161//____________________________________________________________________
162void
1caed034 163CorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
164{
ddec1c19 165 // remove edges of correction histogram by removing
166 // - bins with content less than cut
167 // - bins next to bins with zero bin content
168
b6d3306b 169 Int_t nBinsX = fhCorr->GetNbinsX();
170 Int_t nBinsY = fhCorr->GetNbinsY();
ddec1c19 171
172 // set bin content to zero for bins with content smaller cut
173 for (Int_t bx=0; bx<=nBinsX; bx++) {
174 for (Int_t by=0; by<=nBinsY; by++) {
b6d3306b 175 if (fhCorr->GetBinContent(bx,by)>cut) {
176 fhCorr->SetBinContent(bx,by,0);
177 fhCorr->SetBinError(bx,by,0);
ddec1c19 178 }
179 }
180 }
181
182 // set bin content to zero for bins next to bins with zero
b6d3306b 183 TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
ddec1c19 184 tmp->Reset();
185
186 Bool_t done = kFALSE;
187 Int_t nBinsXCount = 0;
188 Int_t nBinsYCount = 0;
189 while (!done) {
190 if (nBinsXCount<nBinsXedge)
191 for (Int_t bx=0; bx<=nBinsX; bx++) {
192 for (Int_t by=0; by<=nBinsY; by++) {
b6d3306b 193 if ((fhCorr->GetBinContent(bx+1,by)==0)||
194 (fhCorr->GetBinContent(bx-1,by)==0))
ddec1c19 195 tmp->SetBinContent(bx,by,1);
196
197 }
198 }
199 if (nBinsYCount<nBinsYedge)
200 for (Int_t bx=0; bx<=nBinsX; bx++) {
201 for (Int_t by=0; by<=nBinsY; by++) {
b6d3306b 202 if ((fhCorr->GetBinContent(bx,by+1)==0)||
203 (fhCorr->GetBinContent(bx,by-1)==0))
ddec1c19 204 tmp->SetBinContent(bx,by,1);
205 }
206 }
207 for (Int_t bx=0; bx<=nBinsX; bx++) {
208 for (Int_t by=0; by<=nBinsY; by++) {
209 if (tmp->GetBinContent(bx,by)==1) {
b6d3306b 210 fhCorr->SetBinContent(bx,by,0);
211 fhCorr->SetBinError(bx,by,0);
ddec1c19 212 }
213 }
214 }
215 nBinsXCount++;
216 nBinsYCount++;
217 if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
218 }
219 tmp->Delete();
220
221}
222
223//____________________________________________________________________
224Bool_t CorrectionMatrix2D::LoadHistograms(Char_t* fileName, Char_t* dir) {
1caed034 225 //
226 // loads the histograms from a file
227 //
228
ddec1c19 229 TFile* fin = TFile::Open(fileName);
230
1caed034 231 if(!fin) {
ddec1c19 232 //Info("LoadHistograms",Form(" %s file does not exist",fileName));
1caed034 233 return kFALSE;
ddec1c19 234 }
235
b6d3306b 236 if(fhGene) {delete fhGene; fhGene=0;}
237 if(fhCorr) {delete fhCorr; fhCorr=0;}
b6d3306b 238 if(fhMeas) {delete fhMeas; fhMeas=0;}
ddec1c19 239
ea2ef7cc 240 fhMeas = (TH2F*)fin->Get(Form("%s/meas_%s",dir,fName.Data()));
b6d3306b 241 if(!fhMeas) Info("LoadHistograms","No meas hist available");
ea2ef7cc 242 fhGene = (TH2F*)fin->Get(Form("%s/gene_%s",dir,fName.Data()));
b6d3306b 243 if(!fhGene) Info("LoadHistograms","No gene hist available");
ea2ef7cc 244 fhCorr = (TH2F*)fin->Get(Form("%s/corr_%s",dir,fName.Data()));
b6d3306b 245 if(!fhCorr) {Info("LoadHistograms","No corr hist available");
ddec1c19 246 return kFALSE;}
247
248 return kTRUE;
249}
250
251
252//____________________________________________________________________
253void
254CorrectionMatrix2D::SaveHistograms() {
1caed034 255 //
256 // saves the histograms
257 //
ddec1c19 258
259 gDirectory->mkdir(fName.Data());
260 gDirectory->cd(fName.Data());
261
b6d3306b 262 fhMeas ->Write();
263 fhGene ->Write();
ddec1c19 264
b6d3306b 265 if (fhCorr)
266 fhCorr->Write();
ddec1c19 267
ddec1c19 268 gDirectory->cd("../");
269}
270
271//____________________________________________________________________
272void CorrectionMatrix2D::DrawHistograms()
273{
1caed034 274 //
275 // draws all the four histograms on one TCanvas
276 //
277
ddec1c19 278 TCanvas* canvas = new TCanvas("Correction", "Correction", 800, 800);
279 canvas->Divide(2, 2);
280
281 canvas->cd(1);
b6d3306b 282 if (fhMeas)
283 fhMeas->Draw("COLZ");
ddec1c19 284
285 canvas->cd(2);
b6d3306b 286 if (fhGene)
287 fhGene->Draw("COLZ");
ddec1c19 288
289 canvas->cd(3);
b6d3306b 290 if (fhCorr)
291 fhCorr->Draw("COLZ");
ea2ef7cc 292
293 // add: draw here the stat. errors of the correction histogram
294
ddec1c19 295}
296
297