30f3ec3bd69d85a92446f622b17b6385c64898c9
[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::operator=(const AliCorrectionMatrix2D &c)
36 {
37   // assigment operator
38
39   if (this != &c)
40     ((AliCorrectionMatrix2D &) c).Copy(*this);
41
42   return *this;
43 }
44
45 //____________________________________________________________________
46 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
47                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
48                                        Int_t nBinY, Float_t Ymin, Float_t Ymax) 
49   : AliCorrectionMatrix(name, title)
50 {
51   //
52   // constructor
53   //
54
55   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
56   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
57   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
58
59   fhMeas->Sumw2();
60   fhGene->Sumw2();
61   fhCorr->Sumw2();
62 }
63
64 //____________________________________________________________________
65 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
66                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
67   : AliCorrectionMatrix(name, title)
68 {
69   //
70   // constructor
71   //
72
73   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
74   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
75   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);
76
77   fhMeas->Sumw2();
78   fhGene->Sumw2();
79   fhCorr->Sumw2();
80 }
81
82 //____________________________________________________________________
83 AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
84 {
85   //
86   // destructor
87   //
88
89   // histograms already deleted in base class
90 }
91
92 TH2F* AliCorrectionMatrix2D::GetGeneratedHistogram() const
93 {
94   // return generated histogram casted to correct type
95   return dynamic_cast<TH2F*> (fhGene);
96 }
97
98 TH2F* AliCorrectionMatrix2D::GetMeasuredHistogram() const
99 {
100   // return measured histogram casted to correct type
101   return dynamic_cast<TH2F*> (fhMeas);
102 }
103
104 //____________________________________________________________________
105 TH1F* AliCorrectionMatrix2D::Get1DCorrection(Char_t* opt, Float_t min, Float_t max)
106 {
107   //
108   // integrate the correction over one variable 
109   // 
110
111   TH1D* meas1D = 0;
112   TH1D* gene1D = 0; 
113
114   if (strcmp(opt,"x")==0) {
115     Int_t binMin = GetMeasuredHistogram()->GetYaxis()->FindBin(min);
116     Int_t binMax = GetMeasuredHistogram()->GetYaxis()->FindBin(max);
117
118     if (min==0 && max==0) {
119       meas1D = GetMeasuredHistogram()->ProjectionX();
120       gene1D = GetGeneratedHistogram()->ProjectionX();
121     }
122     else {
123       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including y-bins %d to %d \n", binMin, binMax));
124
125       meas1D = GetMeasuredHistogram()->ProjectionX("pm",binMin,binMax);
126       gene1D = GetGeneratedHistogram()->ProjectionX("pg",binMin,binMax);
127     }
128   }
129   if (strcmp(opt,"y")==0) {
130     Int_t binMin = GetMeasuredHistogram()->GetXaxis()->FindBin(min);
131     Int_t binMax = GetMeasuredHistogram()->GetXaxis()->FindBin(max);
132
133     if (min==0 && max==0) {
134       meas1D = GetMeasuredHistogram()->ProjectionY();
135       gene1D = GetGeneratedHistogram()->ProjectionY();
136     }
137     else {
138       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax));
139
140       meas1D = GetMeasuredHistogram()->ProjectionY("pm", binMin, binMax);
141       gene1D = GetGeneratedHistogram()->ProjectionY("pg", binMin, binMax);
142     }
143   }
144   gene1D->Sumw2();
145
146   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
147   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
148
149   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
150   
151   return (TH1F*)gene1D;   
152 }
153
154 //____________________________________________________________________
155 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
156 {
157   // add value to measured histogram
158   GetMeasuredHistogram()->Fill(ax, ay);
159 }
160
161 //____________________________________________________________________
162 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
163 {
164   // add value to generated histogram
165   GetGeneratedHistogram()->Fill(ax, ay);
166 }
167
168 //____________________________________________________________________
169 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
170 {
171   // returns a value of the correction map
172   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
173 }
174
175 //____________________________________________________________________
176 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
177 {
178   // remove edges of correction histogram by removing 
179   // - bins with content less than cut
180   // - bins next to bins with zero bin content
181   
182   Int_t nBinsX = fhCorr->GetNbinsX();
183   Int_t nBinsY = fhCorr->GetNbinsY();
184
185   // set bin content to zero for bins with content smaller cut
186   for (Int_t bx=0; bx<=nBinsX; bx++) {
187     for (Int_t by=0; by<=nBinsY; by++) {
188       if (fhCorr->GetBinContent(bx,by)>cut) {
189           fhCorr->SetBinContent(bx,by,0);
190           fhCorr->SetBinError(bx,by,0);
191       }
192     }
193   }
194
195   // set bin content to zero for bins next to bins with zero
196   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
197   tmp->Reset();
198   
199   Bool_t done = kFALSE;
200   Int_t nBinsXCount = 0;
201   Int_t nBinsYCount = 0;
202   while (!done) {    
203     if (nBinsXCount<nBinsXedge) 
204       for (Int_t bx=0; bx<=nBinsX; bx++) {
205         for (Int_t by=0; by<=nBinsY; by++) {
206           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
207               (fhCorr->GetBinContent(bx-1,by)==0))
208             tmp->SetBinContent(bx,by,1);        
209           
210         }
211       }
212     if (nBinsYCount<nBinsYedge) 
213       for (Int_t bx=0; bx<=nBinsX; bx++) {
214         for (Int_t by=0; by<=nBinsY; by++) {
215           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
216               (fhCorr->GetBinContent(bx,by-1)==0))
217             tmp->SetBinContent(bx,by,1);        
218         }
219       }    
220     for (Int_t bx=0; bx<=nBinsX; bx++) {
221       for (Int_t by=0; by<=nBinsY; by++) {
222         if (tmp->GetBinContent(bx,by)==1) {
223           fhCorr->SetBinContent(bx,by,0);
224           fhCorr->SetBinError(bx,by,0);
225         }
226       }
227     }
228     nBinsXCount++;
229     nBinsYCount++;
230     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
231   }
232   tmp->Delete();  
233
234 }
235