new way of calculation Nsigma correction
[u/mrichter/AliRoot.git] / PWG0 / CorrectionMatrix2D.cxx
1 // ------------------------------------------------------
2 //
3 // Class to handle 2d-corrections. 
4 //
5 // ------------------------------------------------------
6 //
7
8 /* $Id$ */
9
10 #include <TFile.h>
11 #include <TCanvas.h>
12
13 #include <AliLog.h>
14
15 #include "CorrectionMatrix2D.h"
16
17 //____________________________________________________________________
18 ClassImp(CorrectionMatrix2D)
19
20 //____________________________________________________________________
21 CorrectionMatrix2D::CorrectionMatrix2D(const CorrectionMatrix2D& c) 
22   : TNamed(c)
23 {
24   // copy constructor
25   ((CorrectionMatrix2D &)c).Copy(*this);
26 }
27
28 //____________________________________________________________________
29 CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name, Char_t* title,
30                                        Int_t nBinX, Float_t Xmin, Float_t Xmax,
31                                        Int_t nBinY, Float_t Ymin, Float_t Ymax) 
32   : TNamed(name, title)
33 {
34   //
35   // constructor
36   //
37
38
39   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
40   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
41   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
42
43   fhMeas->Sumw2();
44   fhGene->Sumw2();
45   fhCorr->Sumw2();
46 }
47
48 //____________________________________________________________________
49 CorrectionMatrix2D::CorrectionMatrix2D(Char_t* name,Char_t* title, 
50                                        Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
51   : TNamed(name, title) 
52 {
53   //
54   // constructor
55   //
56
57   fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
58   fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
59   fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);
60
61   fhMeas->Sumw2();
62   fhGene->Sumw2();
63   fhCorr->Sumw2();
64 }
65
66
67 //____________________________________________________________________
68 CorrectionMatrix2D::~CorrectionMatrix2D() {
69   //
70   // destructor
71   //
72   if (fhMeas)  delete fhMeas;
73   if (fhGene)  delete fhGene;
74   if (fhCorr)  delete fhCorr;
75 }
76
77 //____________________________________________________________________
78 CorrectionMatrix2D &CorrectionMatrix2D::operator=(const CorrectionMatrix2D &c)
79 {
80   // assigment operator
81
82   if (this != &c) 
83     ((CorrectionMatrix2D &) c).Copy(*this);
84
85   return *this;
86 }
87
88 //____________________________________________________________________
89 TH1F* CorrectionMatrix2D::Get1DCorrection(Char_t* opt) {
90   //
91   // integrate the correction over one variable 
92   // 
93
94   TH1D* meas1D = 0;
95   TH1D* gene1D = 0; 
96
97   if (strcmp(opt,"x")==0) {
98     meas1D = fhMeas->ProjectionX();
99     gene1D = fhGene->ProjectionX();      
100   }
101   if (strcmp(opt,"y")==0) {
102     meas1D = fhMeas->ProjectionY();
103     gene1D = fhGene->ProjectionY();      
104   }
105   gene1D->Sumw2();
106
107   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
108   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
109
110   gene1D->Divide(gene1D, meas1D, 1, 1, "B");
111   
112   return (TH1F*)gene1D;   
113 }
114
115
116 //____________________________________________________________________
117 void
118 CorrectionMatrix2D::Copy(TObject& c) const 
119 {
120   // copy function
121
122   CorrectionMatrix2D& target = (CorrectionMatrix2D &) c;
123
124   target.fhMeas  = fhMeas;
125   target.fhGene  = fhGene;
126   target.fhCorr  = fhCorr;
127 }
128
129
130 //________________________________________________________________________
131 void CorrectionMatrix2D::SetAxisTitles(Char_t* titleX, Char_t* titleY) 
132
133   //
134   // method for setting the axis titles of the histograms
135   //
136
137   fhMeas ->SetXTitle(titleX);  fhMeas ->SetYTitle(titleY);
138   fhGene ->SetXTitle(titleX);  fhGene ->SetYTitle(titleY);
139   fhCorr ->SetXTitle(titleX);  fhCorr ->SetYTitle(titleY);
140 }
141
142 //____________________________________________________________________
143 Long64_t CorrectionMatrix2D::Merge(TCollection* list) {
144   // Merge a list of CorrectionMatrix2D objects with this (needed for
145   // PROOF). 
146   // Returns the number of merged objects (including this).
147
148   if (!list)
149     return 0;
150   
151   if (list->IsEmpty())
152     return 1;
153
154   TIterator* iter = list->MakeIterator();
155   TObject* obj;
156
157   // collections of measured and generated histograms
158   TList* collectionMeas = new TList;
159   TList* collectionGene = new TList;
160   
161   Int_t count = 0;
162   while ((obj = iter->Next())) {
163     
164     CorrectionMatrix2D* entry = dynamic_cast<CorrectionMatrix2D*> (obj);
165     if (entry == 0) 
166       continue;
167
168     collectionMeas->Add(entry->GetMeasuredHistogram());
169     collectionGene->Add(entry->GetGeneratedHistogram());
170
171     count++;
172   }
173   fhMeas->Merge(collectionMeas);
174   fhGene->Merge(collectionGene);
175
176   // is this really faster than just adding the histograms in the list???
177   delete collectionMeas;
178   delete collectionGene;
179
180
181   return count+1;
182 }
183
184
185 //____________________________________________________________________
186 void CorrectionMatrix2D::Divide() {  
187   //
188   // divide the histograms to get the correction
189   // 
190
191   if (!fhMeas || !fhGene)  return; 
192
193   fhCorr->Divide(fhGene, fhMeas, 1,1,"B");
194   
195 }
196
197 //____________________________________________________________________
198 void
199 CorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge) 
200 {
201   // remove edges of correction histogram by removing 
202   // - bins with content less than cut
203   // - bins next to bins with zero bin content
204   
205   Int_t nBinsX = fhCorr->GetNbinsX();
206   Int_t nBinsY = fhCorr->GetNbinsY();
207
208   // set bin content to zero for bins with content smaller cut
209   for (Int_t bx=0; bx<=nBinsX; bx++) {
210     for (Int_t by=0; by<=nBinsY; by++) {
211       if (fhCorr->GetBinContent(bx,by)>cut) {
212           fhCorr->SetBinContent(bx,by,0);
213           fhCorr->SetBinError(bx,by,0);
214       }
215     }
216   }
217
218   // set bin content to zero for bins next to bins with zero
219   TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
220   tmp->Reset();
221   
222   Bool_t done = kFALSE;
223   Int_t nBinsXCount = 0;
224   Int_t nBinsYCount = 0;
225   while (!done) {    
226     if (nBinsXCount<nBinsXedge) 
227       for (Int_t bx=0; bx<=nBinsX; bx++) {
228         for (Int_t by=0; by<=nBinsY; by++) {
229           if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
230               (fhCorr->GetBinContent(bx-1,by)==0))
231             tmp->SetBinContent(bx,by,1);        
232           
233         }
234       }
235     if (nBinsYCount<nBinsYedge) 
236       for (Int_t bx=0; bx<=nBinsX; bx++) {
237         for (Int_t by=0; by<=nBinsY; by++) {
238           if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
239               (fhCorr->GetBinContent(bx,by-1)==0))
240             tmp->SetBinContent(bx,by,1);        
241         }
242       }    
243     for (Int_t bx=0; bx<=nBinsX; bx++) {
244       for (Int_t by=0; by<=nBinsY; by++) {
245         if (tmp->GetBinContent(bx,by)==1) {
246           fhCorr->SetBinContent(bx,by,0);
247           fhCorr->SetBinError(bx,by,0);
248         }
249       }
250     }
251     nBinsXCount++;
252     nBinsYCount++;
253     if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
254   }
255   tmp->Delete();  
256
257 }
258
259 //____________________________________________________________________
260 Bool_t CorrectionMatrix2D::LoadHistograms(Char_t* fileName, Char_t* dir) {
261   //
262   // loads the histograms from a file
263   //
264   
265   TFile* fin = TFile::Open(fileName);  
266   
267   if(!fin) {
268     //Info("LoadHistograms",Form(" %s file does not exist",fileName));
269     return kFALSE;
270   }
271   
272   if(fhGene)  {delete fhGene;  fhGene=0;}
273   if(fhCorr)  {delete fhCorr;  fhCorr=0;}
274   if(fhMeas)  {delete fhMeas;  fhMeas=0;}
275   
276   fhMeas  = (TH2F*)fin->Get(Form("%s/meas_%s", dir,GetName()));
277   if(!fhMeas)  Info("LoadHistograms","No meas  hist available");
278   
279   fhGene  = (TH2F*)fin->Get(Form("%s/gene_%s",dir, GetName()));
280   if(!fhGene)  Info("LoadHistograms","No gene  hist available");
281   
282   fhCorr  = (TH2F*)fin->Get(Form("%s/corr_%s",dir, GetName()));
283   if(!fhCorr) 
284   {
285     Info("LoadHistograms","No corr  hist available");
286     return kFALSE;
287   }
288       
289   return kTRUE;
290 }
291
292
293 //____________________________________________________________________
294 void
295 CorrectionMatrix2D::SaveHistograms() {
296   //
297   // saves the histograms 
298   //
299   
300   fhMeas ->Write();
301   fhGene ->Write();
302
303   if (fhCorr)
304     fhCorr->Write();
305 }
306
307 //____________________________________________________________________
308 void CorrectionMatrix2D::DrawHistograms()
309 {
310   //
311   // draws all the four histograms on one TCanvas
312   //
313
314   TCanvas* canvas = new TCanvas(Form("correction_%s",fName.Data()), 
315                                 Form("correction_%s",fName.Data()), 800, 800);
316   canvas->Divide(2, 2);
317   
318   canvas->cd(1);
319   if (fhMeas)
320     fhMeas->Draw("COLZ");
321   
322   canvas->cd(2);
323   if (fhGene)
324     fhGene->Draw("COLZ");
325
326   canvas->cd(3);
327   if (fhCorr)
328     fhCorr->Draw("COLZ");
329
330   canvas->cd(4);
331
332   // add: draw here the stat. errors of the correction histogram
333   
334 }
335
336