6cfee1835445a13b834a9e5086558b027def1e9f
[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)
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     meas1D = GetMeasuredHistogram()->ProjectionX();
116     gene1D = GetGeneratedHistogram()->ProjectionX();
117   }
118   if (strcmp(opt,"y")==0) {
119     meas1D = GetMeasuredHistogram()->ProjectionY();
120     gene1D = GetGeneratedHistogram()->ProjectionY();
121   }
122   gene1D->Sumw2();
123
124   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
125   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
126
127   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
128   
129   return (TH1F*)gene1D;   
130 }
131
132 //____________________________________________________________________
133 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
134 {
135   // add value to measured histogram
136   GetMeasuredHistogram()->Fill(ax, ay);
137 }
138
139 //____________________________________________________________________
140 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
141 {
142   // add value to generated histogram
143   GetGeneratedHistogram()->Fill(ax, ay);
144 }
145
146 //____________________________________________________________________
147 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
148 {
149   // returns a value of the correction map
150   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
151 }
152
153 //____________________________________________________________________
154 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
155 {
156   // remove edges of correction histogram by removing 
157   // - bins with content less than cut
158   // - bins next to bins with zero bin content
159   
160   Int_t nBinsX = fhCorr->GetNbinsX();
161   Int_t nBinsY = fhCorr->GetNbinsY();
162
163   // set bin content to zero for bins with content smaller cut
164   for (Int_t bx=0; bx<=nBinsX; bx++) {
165     for (Int_t by=0; by<=nBinsY; by++) {
166       if (fhCorr->GetBinContent(bx,by)>cut) {
167           fhCorr->SetBinContent(bx,by,0);
168           fhCorr->SetBinError(bx,by,0);
169       }
170     }
171   }
172
173   // set bin content to zero for bins next to bins with zero
174   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
175   tmp->Reset();
176   
177   Bool_t done = kFALSE;
178   Int_t nBinsXCount = 0;
179   Int_t nBinsYCount = 0;
180   while (!done) {    
181     if (nBinsXCount<nBinsXedge) 
182       for (Int_t bx=0; bx<=nBinsX; bx++) {
183         for (Int_t by=0; by<=nBinsY; by++) {
184           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
185               (fhCorr->GetBinContent(bx-1,by)==0))
186             tmp->SetBinContent(bx,by,1);        
187           
188         }
189       }
190     if (nBinsYCount<nBinsYedge) 
191       for (Int_t bx=0; bx<=nBinsX; bx++) {
192         for (Int_t by=0; by<=nBinsY; by++) {
193           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
194               (fhCorr->GetBinContent(bx,by-1)==0))
195             tmp->SetBinContent(bx,by,1);        
196         }
197       }    
198     for (Int_t bx=0; bx<=nBinsX; bx++) {
199       for (Int_t by=0; by<=nBinsY; by++) {
200         if (tmp->GetBinContent(bx,by)==1) {
201           fhCorr->SetBinContent(bx,by,0);
202           fhCorr->SetBinError(bx,by,0);
203         }
204       }
205     }
206     nBinsXCount++;
207     nBinsYCount++;
208     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
209   }
210   tmp->Delete();  
211
212 }
213