renamed CorrectionMatrix class
[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix2D.cxx
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 //____________________________________________________________________
17 ClassImp(AliCorrectionMatrix2D)
18
19 //____________________________________________________________________
20 AliCorrectionMatrix2D::AliCorrectionMatrix2D() :
21   AliCorrectionMatrix()
22 {
23   // default constructor
24 }
25
26 //____________________________________________________________________
27 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const AliCorrectionMatrix2D& c)
28   : AliCorrectionMatrix(c)
29 {
30   // copy constructor
31   ((AliCorrectionMatrix2D &)c).Copy(*this);
32 }
33
34 //____________________________________________________________________
35 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
36                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
37                                        Int_t nBinY, Float_t Ymin, Float_t Ymax) 
38   : AliCorrectionMatrix(name, title)
39 {
40   //
41   // constructor
42   //
43
44   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
45   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
46   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
47
48   fhMeas->Sumw2();
49   fhGene->Sumw2();
50   fhCorr->Sumw2();
51 }
52
53 //____________________________________________________________________
54 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
55                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
56   : AliCorrectionMatrix(name, title)
57 {
58   //
59   // constructor
60   //
61
62   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
63   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
64   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);
65
66   fhMeas->Sumw2();
67   fhGene->Sumw2();
68   fhCorr->Sumw2();
69 }
70
71 //____________________________________________________________________
72 AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
73 {
74   //
75   // destructor
76   //
77
78   // histograms already deleted in base class
79 }
80
81 TH2F* AliCorrectionMatrix2D::GetGeneratedHistogram() const
82 {
83   // return generated histogram casted to correct type
84   return dynamic_cast<TH2F*> (fhGene);
85 }
86
87 TH2F* AliCorrectionMatrix2D::GetMeasuredHistogram() const
88 {
89   // return measured histogram casted to correct type
90   return dynamic_cast<TH2F*> (fhMeas);
91 }
92
93 //____________________________________________________________________
94 TH1F* AliCorrectionMatrix2D::Get1DCorrection(Char_t* opt)
95 {
96   //
97   // integrate the correction over one variable 
98   // 
99
100   TH1D* meas1D = 0;
101   TH1D* gene1D = 0; 
102
103   if (strcmp(opt,"x")==0) {
104     meas1D = GetMeasuredHistogram()->ProjectionX();
105     gene1D = GetGeneratedHistogram()->ProjectionX();
106   }
107   if (strcmp(opt,"y")==0) {
108     meas1D = GetMeasuredHistogram()->ProjectionY();
109     gene1D = GetGeneratedHistogram()->ProjectionY();
110   }
111   gene1D->Sumw2();
112
113   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
114   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
115
116   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
117   
118   return (TH1F*)gene1D;   
119 }
120
121 //____________________________________________________________________
122 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
123 {
124   // add value to measured histogram
125   GetMeasuredHistogram()->Fill(ax, ay);
126 }
127
128 //____________________________________________________________________
129 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
130 {
131   // add value to generated histogram
132   GetGeneratedHistogram()->Fill(ax, ay);
133 }
134
135 //____________________________________________________________________
136 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
137 {
138   // returns a value of the correction map
139   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
140 }
141
142 //____________________________________________________________________
143 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
144 {
145   // remove edges of correction histogram by removing 
146   // - bins with content less than cut
147   // - bins next to bins with zero bin content
148   
149   Int_t nBinsX = fhCorr->GetNbinsX();
150   Int_t nBinsY = fhCorr->GetNbinsY();
151
152   // set bin content to zero for bins with content smaller cut
153   for (Int_t bx=0; bx<=nBinsX; bx++) {
154     for (Int_t by=0; by<=nBinsY; by++) {
155       if (fhCorr->GetBinContent(bx,by)>cut) {
156           fhCorr->SetBinContent(bx,by,0);
157           fhCorr->SetBinError(bx,by,0);
158       }
159     }
160   }
161
162   // set bin content to zero for bins next to bins with zero
163   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
164   tmp->Reset();
165   
166   Bool_t done = kFALSE;
167   Int_t nBinsXCount = 0;
168   Int_t nBinsYCount = 0;
169   while (!done) {    
170     if (nBinsXCount<nBinsXedge) 
171       for (Int_t bx=0; bx<=nBinsX; bx++) {
172         for (Int_t by=0; by<=nBinsY; by++) {
173           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
174               (fhCorr->GetBinContent(bx-1,by)==0))
175             tmp->SetBinContent(bx,by,1);        
176           
177         }
178       }
179     if (nBinsYCount<nBinsYedge) 
180       for (Int_t bx=0; bx<=nBinsX; bx++) {
181         for (Int_t by=0; by<=nBinsY; by++) {
182           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
183               (fhCorr->GetBinContent(bx,by-1)==0))
184             tmp->SetBinContent(bx,by,1);        
185         }
186       }    
187     for (Int_t bx=0; bx<=nBinsX; bx++) {
188       for (Int_t by=0; by<=nBinsY; by++) {
189         if (tmp->GetBinContent(bx,by)==1) {
190           fhCorr->SetBinContent(bx,by,0);
191           fhCorr->SetBinError(bx,by,0);
192         }
193       }
194     }
195     nBinsXCount++;
196     nBinsYCount++;
197     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
198   }
199   tmp->Delete();  
200
201 }
202