adding macro that applies the corrections and visualizes the result
[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 TH2F* AliCorrectionMatrix2D::GetGeneratedHistogram() const
105 {
106   // return generated histogram casted to correct type
107   return dynamic_cast<TH2F*> (fhGene);
108 }
109
110 TH2F* AliCorrectionMatrix2D::GetMeasuredHistogram() const
111 {
112   // return measured histogram casted to correct type
113   return dynamic_cast<TH2F*> (fhMeas);
114 }
115
116 //____________________________________________________________________
117 TH1F* 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_pm", GetName()),binMin,binMax);
138       gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_pg", GetName()),binMin,binMax);
139     }
140   }
141   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       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax));
151
152       meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_pm", GetName()), binMin, binMax);
153       gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_pg", GetName()), binMin, binMax);
154     }
155   }
156   gene1D->Sumw2();
157
158   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
159   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
160
161   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
162   
163   return (TH1F*)gene1D;   
164 }
165
166 //____________________________________________________________________
167 void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
168 {
169   // add value to measured histogram
170   ((TH2F*)fhMeas)->Fill(ax, ay);
171 }
172
173 //____________________________________________________________________
174 void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
175 {
176   // add value to generated histogram
177   ((TH2F*)fhGene)->Fill(ax, ay);
178 }
179
180 //____________________________________________________________________
181 Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
182 {
183   // returns a value of the correction map
184   return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
185 }
186
187 //____________________________________________________________________
188 void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
189 {
190   // remove edges of correction histogram by removing 
191   // - bins with content less than cut
192   // - bins next to bins with zero bin content
193   
194   Int_t nBinsX = fhCorr->GetNbinsX();
195   Int_t nBinsY = fhCorr->GetNbinsY();
196
197   // set bin content to zero for bins with content smaller cut
198   for (Int_t bx=0; bx<=nBinsX; bx++) {
199     for (Int_t by=0; by<=nBinsY; by++) {
200       if (fhCorr->GetBinContent(bx,by)>cut) {
201           fhCorr->SetBinContent(bx,by,0);
202           fhCorr->SetBinError(bx,by,0);
203       }
204     }
205   }
206
207   // set bin content to zero for bins next to bins with zero
208   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
209   tmp->Reset();
210   
211   Bool_t done = kFALSE;
212   Int_t nBinsXCount = 0;
213   Int_t nBinsYCount = 0;
214   while (!done) {    
215     if (nBinsXCount<nBinsXedge) 
216       for (Int_t bx=0; bx<=nBinsX; bx++) {
217         for (Int_t by=0; by<=nBinsY; by++) {
218           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
219               (fhCorr->GetBinContent(bx-1,by)==0))
220             tmp->SetBinContent(bx,by,1);        
221           
222         }
223       }
224     if (nBinsYCount<nBinsYedge) 
225       for (Int_t bx=0; bx<=nBinsX; bx++) {
226         for (Int_t by=0; by<=nBinsY; by++) {
227           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
228               (fhCorr->GetBinContent(bx,by-1)==0))
229             tmp->SetBinContent(bx,by,1);        
230         }
231       }    
232     for (Int_t bx=0; bx<=nBinsX; bx++) {
233       for (Int_t by=0; by<=nBinsY; by++) {
234         if (tmp->GetBinContent(bx,by)==1) {
235           fhCorr->SetBinContent(bx,by,0);
236           fhCorr->SetBinError(bx,by,0);
237         }
238       }
239     }
240     nBinsXCount++;
241     nBinsYCount++;
242     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
243   }
244   tmp->Delete();  
245
246 }
247