o) adding log tags to all files
[u/mrichter/AliRoot.git] / PWG0 / CorrectionMatrix2D.cxx
1 /* $Id$ */
2
3 #include "CorrectionMatrix2D.h"
4
5 //____________________________________________________________________
6 ClassImp(CorrectionMatrix2D);
7
8 //____________________________________________________________________
9 CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name, Char_t* title,
10                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
11                                        Int_t nBinY, Float_t Ymin, Float_t Ymax) 
12   : TNamed(name, title)
13 {
14
15   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
16   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
17   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
18   fhRatio = new TH2F(Form("ratio_%s",name),Form("ratio_%s",title), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
19 }
20
21 //____________________________________________________________________
22 CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name,Char_t* title, 
23                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
24   : TNamed(name, title)
25 {
26   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
27   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
28   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);
29   fhRatio = new TH2F(Form("ratio_%s",name),Form("ratio_%s",title), nBinX, X, nBinY, Y);
30 }
31
32
33 //____________________________________________________________________
34 CorrectionMatrix2D::~CorrectionMatrix2D() {
35   // Destructor
36   //
37   if (fhMeas)  delete fhMeas;
38   if (fhGene)  delete fhCorr;
39   if (fhRatio) delete fhRatio;
40   if (fhCorr)  delete fhCorr;
41 }
42
43
44 //________________________________________________________________________
45 void CorrectionMatrix2D::SetAxisTitles(Char_t* titleX, Char_t* titleY) 
46
47   fhMeas ->SetXTitle(titleX);  fhMeas ->SetYTitle(titleY);
48   fhGene ->SetXTitle(titleX);  fhGene ->SetYTitle(titleY);
49   fhCorr ->SetXTitle(titleX);  fhCorr ->SetYTitle(titleY);
50   fhRatio->SetXTitle(titleX);  fhRatio->SetYTitle(titleY);
51 }
52
53
54 //____________________________________________________________________
55 void CorrectionMatrix2D::Finish() {  
56
57 if (!fhMeas || !fhGene)  return;
58  
59
60   fhRatio->Divide(fhMeas, fhGene, 1,1,"B");
61   fhCorr->Divide(fhGene, fhMeas, 1,1,"B");
62
63   Int_t nBinsVtx = fhCorr->GetNbinsX();
64   Int_t nBinsEta = fhCorr->GetNbinsY();
65
66   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
67
68   // cut at 0.2
69   for (Int_t bx=0; bx<=nBinsVtx; bx++) {
70     for (Int_t by=0; by<=nBinsEta; by++) {
71       if (tmp->GetBinContent(bx,by)<0.2) {
72         fhCorr->SetBinContent(bx,by,0);
73         fhCorr->SetBinError(bx,by,0);
74         
75         tmp->SetBinContent(bx,by,0);
76       }
77       else 
78         tmp->SetBinContent(bx,by,1);
79     }
80   }
81   
82 }
83
84 //____________________________________________________________________
85 void
86 CorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge) {
87
88   // remove edges of correction histogram by removing 
89   // - bins with content less than cut
90   // - bins next to bins with zero bin content
91   
92   Int_t nBinsX = fhCorr->GetNbinsX();
93   Int_t nBinsY = fhCorr->GetNbinsY();
94
95   // set bin content to zero for bins with content smaller cut
96   for (Int_t bx=0; bx<=nBinsX; bx++) {
97     for (Int_t by=0; by<=nBinsY; by++) {
98       if (fhCorr->GetBinContent(bx,by)>cut) {
99           fhCorr->SetBinContent(bx,by,0);
100           fhCorr->SetBinError(bx,by,0);
101       }
102     }
103   }
104
105   // set bin content to zero for bins next to bins with zero
106   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
107   tmp->Reset();
108   
109   Bool_t done = kFALSE;
110   Int_t nBinsXCount = 0;
111   Int_t nBinsYCount = 0;
112   while (!done) {    
113     if (nBinsXCount<nBinsXedge) 
114       for (Int_t bx=0; bx<=nBinsX; bx++) {
115         for (Int_t by=0; by<=nBinsY; by++) {
116           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
117               (fhCorr->GetBinContent(bx-1,by)==0))
118             tmp->SetBinContent(bx,by,1);        
119           
120         }
121       }
122     if (nBinsYCount<nBinsYedge) 
123       for (Int_t bx=0; bx<=nBinsX; bx++) {
124         for (Int_t by=0; by<=nBinsY; by++) {
125           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
126               (fhCorr->GetBinContent(bx,by-1)==0))
127             tmp->SetBinContent(bx,by,1);        
128         }
129       }    
130     for (Int_t bx=0; bx<=nBinsX; bx++) {
131       for (Int_t by=0; by<=nBinsY; by++) {
132         if (tmp->GetBinContent(bx,by)==1) {
133           fhCorr->SetBinContent(bx,by,0);
134           fhCorr->SetBinError(bx,by,0);
135         }
136       }
137     }
138     nBinsXCount++;
139     nBinsYCount++;
140     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
141   }
142   tmp->Delete();  
143
144 }
145
146 //____________________________________________________________________
147 Bool_t CorrectionMatrix2D::LoadHistograms(Char_t* fileName, Char_t* dir) {
148
149   TFile* fin = TFile::Open(fileName);  
150   
151   if(!fin) 
152   {
153     //Info("LoadHistograms",Form(" %s file does not exist",fileName));
154   return kFALSE;
155   }
156   
157   if(fhGene)  {delete fhGene;  fhGene=0;}
158   if(fhCorr)  {delete fhCorr;  fhCorr=0;}
159   if(fhRatio) {delete fhRatio; fhRatio=0;}
160   if(fhMeas)  {delete fhMeas;  fhMeas=0;}
161   
162   fhMeas  = (TH2F*)fin->Get(Form("%s/meas",dir));
163       if(!fhMeas)  Info("LoadHistograms","No meas  hist available");
164   fhGene  = (TH2F*)fin->Get(Form("%s/gene",dir));
165       if(!fhGene)  Info("LoadHistograms","No gene  hist available");
166   fhRatio = (TH2F*)fin->Get(Form("%s/ratio",dir));
167       if(!fhRatio) Info("LoadHistograms","No ratio hist available");
168   fhCorr  = (TH2F*)fin->Get(Form("%s/corr",dir));
169       if(!fhCorr) {Info("LoadHistograms","No corr  hist available");
170       return kFALSE;}
171       
172   return kTRUE;
173 }
174
175
176 //____________________________________________________________________
177 void
178 CorrectionMatrix2D::SaveHistograms() {
179
180   gDirectory->mkdir(fName.Data());
181   gDirectory->cd(fName.Data());
182   
183   fhMeas ->Write();
184   fhGene ->Write();
185
186   if (fhCorr)
187     fhCorr->Write();
188
189   if (fhRatio)
190     fhRatio->Write();
191
192   gDirectory->cd("../");
193 }
194
195 //____________________________________________________________________
196 void CorrectionMatrix2D::DrawHistograms()
197 {
198   TCanvas* canvas = new TCanvas("Correction", "Correction", 800, 800);
199   canvas->Divide(2, 2);
200
201   canvas->cd(1);
202   if (fhMeas)
203     fhMeas->Draw("COLZ");
204
205   canvas->cd(2);
206   if (fhGene)
207     fhGene->Draw("COLZ");
208
209   canvas->cd(3);
210   if (fhRatio)
211     fhRatio->Draw("COLZ");
212
213   canvas->cd(4);
214   if (fhCorr)
215     fhCorr->Draw("COLZ");
216 }
217
218