removing checking of vertex resolution in getvertex
[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   TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
167   divided->Reset();
168  
169   divided->Divide(gene1D, meas1D, 1, 1, "B");
170
171   Printf("%p %p", gene1D, meas1D);
172   
173   return (TH1F*)divided;   
174 }
175
176 //____________________________________________________________________
177 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
178 {
179   // add value to measured histogram
180   ((TH2F*)fhMeas)->Fill(ax, ay);
181 }
182
183 //____________________________________________________________________
184 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
185 {
186   // add value to generated histogram
187   ((TH2F*)fhGene)->Fill(ax, ay);
188 }
189
190 //____________________________________________________________________
191 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
192 {
193   // returns a value of the correction map
194   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
195 }
196
197 //____________________________________________________________________
198 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
199 {
200   // remove edges of correction histogram by removing 
201   // - bins with content less than cut
202   // - bins next to bins with zero bin content
203   
204   Int_t nBinsX = fhCorr->GetNbinsX();
205   Int_t nBinsY = fhCorr->GetNbinsY();
206
207   // set bin content to zero for bins with content smaller cut
208   for (Int_t bx=0; bx<=nBinsX; bx++) {
209     for (Int_t by=0; by<=nBinsY; by++) {
210       if (fhCorr->GetBinContent(bx,by)>cut) {
211           fhCorr->SetBinContent(bx,by,0);
212           fhCorr->SetBinError(bx,by,0);
213       }
214     }
215   }
216
217   // set bin content to zero for bins next to bins with zero
218   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
219   tmp->Reset();
220   
221   Bool_t done = kFALSE;
222   Int_t nBinsXCount = 0;
223   Int_t nBinsYCount = 0;
224   while (!done) {    
225     if (nBinsXCount<nBinsXedge) 
226       for (Int_t bx=0; bx<=nBinsX; bx++) {
227         for (Int_t by=0; by<=nBinsY; by++) {
228           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
229               (fhCorr->GetBinContent(bx-1,by)==0))
230             tmp->SetBinContent(bx,by,1);        
231           
232         }
233       }
234     if (nBinsYCount<nBinsYedge) 
235       for (Int_t bx=0; bx<=nBinsX; bx++) {
236         for (Int_t by=0; by<=nBinsY; by++) {
237           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
238               (fhCorr->GetBinContent(bx,by-1)==0))
239             tmp->SetBinContent(bx,by,1);        
240         }
241       }    
242     for (Int_t bx=0; bx<=nBinsX; bx++) {
243       for (Int_t by=0; by<=nBinsY; by++) {
244         if (tmp->GetBinContent(bx,by)==1) {
245           fhCorr->SetBinContent(bx,by,0);
246           fhCorr->SetBinError(bx,by,0);
247         }
248       }
249     }
250     nBinsXCount++;
251     nBinsYCount++;
252     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
253   }
254   tmp->Delete();  
255
256 }
257
258 //____________________________________________________________________
259 void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
260 {
261         // rebins the histograms, recalculates the correction
262
263         GetGeneratedHistogram()->Rebin2D(x, y);
264         GetMeasuredHistogram()->Rebin2D(x, y);
265         GetCorrectionHistogram()->Rebin2D(x, y);
266         Divide();
267 }