891de0b512d9fec6dde391f5046b8865fe8af8ed
[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   // do not add this hists to the directory
56   Bool_t oldStatus = TH1::AddDirectoryStatus();
57   TH1::AddDirectory(kFALSE);
58
59   fhMeas  = new TH2F("measured",   Form("%s measured", GetTitle()),   nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
60   fhGene  = new TH2F("generated",  Form("%s generated", GetTitle()),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
61   fhCorr  = new TH2F("correction", Form("%s correction", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
62
63   TH1::AddDirectory(oldStatus);
64
65   fhMeas->Sumw2();
66   fhGene->Sumw2();
67   fhCorr->Sumw2();
68 }
69
70 //____________________________________________________________________
71 AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
72                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
73   : AliCorrectionMatrix(name, title)
74 {
75   //
76   // constructor
77   //
78
79   // do not add this hists to the directory
80   Bool_t oldStatus = TH1::AddDirectoryStatus();
81   TH1::AddDirectory(kFALSE);
82
83         fhMeas  = new TH2F("measured",   Form("%s measured",title),   nBinX, X, nBinY, Y);
84   fhGene  = new TH2F("generated",  Form("%s generated",title),  nBinX, X, nBinY, Y);
85   fhCorr  = new TH2F("correction", Form("%s correction",title), nBinX, X, nBinY, Y);
86
87   TH1::AddDirectory(oldStatus);
88
89         fhMeas->Sumw2();
90   fhGene->Sumw2();
91   fhCorr->Sumw2();
92 }
93
94 //____________________________________________________________________
95 AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
96 {
97   //
98   // destructor
99   //
100
101   // histograms already deleted in base class
102 }
103
104 TH2* AliCorrectionMatrix2D::GetGeneratedHistogram() const
105 {
106   // return generated histogram casted to correct type
107   return dynamic_cast<TH2*> (fhGene);
108 }
109
110 TH2* AliCorrectionMatrix2D::GetMeasuredHistogram() const
111 {
112   // return measured histogram casted to correct type
113   return dynamic_cast<TH2*> (fhMeas);
114 }
115
116 //____________________________________________________________________
117 TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(Char_t* opt, Float_t min, Float_t max)
118 {
119   //
120   // integrate the correction over one variable 
121   // 
122
123   TH1D* meas1D = 0;
124   TH1D* gene1D = 0; 
125
126   if (strcmp(opt,"x")==0) {
127     Int_t binMin = fhMeas->GetYaxis()->FindBin(min);
128     Int_t binMax = fhGene->GetYaxis()->FindBin(max);
129
130     if (min>=max) {
131       meas1D = ((TH2F*)fhMeas)->ProjectionX();
132       gene1D = ((TH2F*)fhGene)->ProjectionX();
133     }
134     else {
135       Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);
136
137       meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_x_pm", GetName()),binMin,binMax);
138       gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_x_pg", GetName()),binMin,binMax);
139     }
140   }
141   else if (strcmp(opt,"y")==0) {
142     Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
143     Int_t binMax = fhMeas->GetXaxis()->FindBin(max);
144
145     if (min>=max) {
146       meas1D = ((TH2F*)fhMeas)->ProjectionY();
147       gene1D = ((TH2F*)fhGene)->ProjectionY();
148     }
149     else {
150       Printf("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax);
151
152       meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_y_pm", GetName()), binMin, binMax);
153       gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_y_pg", GetName()), binMin, binMax);
154     }
155   }
156   else {
157     Printf("ERROR: Invalid option");
158     return 0;
159   }
160
161   gene1D->Sumw2();
162
163   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
164   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
165
166   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
167   
168   return (TH1F*)gene1D;   
169 }
170
171 //____________________________________________________________________
172 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
173 {
174   // add value to measured histogram
175   ((TH2F*)fhMeas)->Fill(ax, ay);
176 }
177
178 //____________________________________________________________________
179 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
180 {
181   // add value to generated histogram
182   ((TH2F*)fhGene)->Fill(ax, ay);
183 }
184
185 //____________________________________________________________________
186 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
187 {
188   // returns a value of the correction map
189   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
190 }
191
192 //____________________________________________________________________
193 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
194 {
195   // remove edges of correction histogram by removing 
196   // - bins with content less than cut
197   // - bins next to bins with zero bin content
198   
199   Int_t nBinsX = fhCorr->GetNbinsX();
200   Int_t nBinsY = fhCorr->GetNbinsY();
201
202   // set bin content to zero for bins with content smaller cut
203   for (Int_t bx=0; bx<=nBinsX; bx++) {
204     for (Int_t by=0; by<=nBinsY; by++) {
205       if (fhCorr->GetBinContent(bx,by)>cut) {
206           fhCorr->SetBinContent(bx,by,0);
207           fhCorr->SetBinError(bx,by,0);
208       }
209     }
210   }
211
212   // set bin content to zero for bins next to bins with zero
213   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
214   tmp->Reset();
215   
216   Bool_t done = kFALSE;
217   Int_t nBinsXCount = 0;
218   Int_t nBinsYCount = 0;
219   while (!done) {    
220     if (nBinsXCount<nBinsXedge) 
221       for (Int_t bx=0; bx<=nBinsX; bx++) {
222         for (Int_t by=0; by<=nBinsY; by++) {
223           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
224               (fhCorr->GetBinContent(bx-1,by)==0))
225             tmp->SetBinContent(bx,by,1);        
226           
227         }
228       }
229     if (nBinsYCount<nBinsYedge) 
230       for (Int_t bx=0; bx<=nBinsX; bx++) {
231         for (Int_t by=0; by<=nBinsY; by++) {
232           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
233               (fhCorr->GetBinContent(bx,by-1)==0))
234             tmp->SetBinContent(bx,by,1);        
235         }
236       }    
237     for (Int_t bx=0; bx<=nBinsX; bx++) {
238       for (Int_t by=0; by<=nBinsY; by++) {
239         if (tmp->GetBinContent(bx,by)==1) {
240           fhCorr->SetBinContent(bx,by,0);
241           fhCorr->SetBinError(bx,by,0);
242         }
243       }
244     }
245     nBinsXCount++;
246     nBinsYCount++;
247     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
248   }
249   tmp->Delete();  
250
251 }
252
253 //____________________________________________________________________
254 void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
255 {
256         // rebins the histograms, recalculates the correction
257
258         GetGeneratedHistogram()->Rebin2D(x, y);
259         GetMeasuredHistogram()->Rebin2D(x, y);
260         GetCorrectionHistogram()->Rebin2D(x, y);
261         Divide();
262 }