2f517ffcdf493d725fd7f5871036a03c4f006587
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliPoissonCalculator.cxx
1 #include "AliPoissonCalculator.h"
2 #include <TH2D.h>
3 #include <TBrowser.h>
4 #include <TROOT.h>
5 #include <TMath.h>
6 #include <TList.h>
7 #include <iostream>
8
9 // 
10 // A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
11 // using Poisson statistics. 
12 //
13 // The input is assumed to be binned in @f$(\eta,\varphi)@f$ as
14 // described by the 2D histogram passwd to the Reset member function.  
15 //
16 // The data is grouped in to regions as defined by the parameters
17 // fEtaLumping and fPhiLumping.  The total number of cells and number
18 // of empty cells is then calculate in each region.  The mean
19 // multiplicity over the region is then determined as 
20 //
21 // @f[
22 // \lange m\rangle = -\log\left(\frac{e}{t}\right)
23 // @f]
24 // where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
25 // total number of cells in the region.  A correction for counting
26 // statistics, is then applied 
27 // @f{eqnarray*}
28 //    c &=& \frac{1}{1 - \exp{-\lange m\rangle}}
29 //      &=& \frac{1}{1 - \frac{e}{t}}
30 // @f{eqnarray*}
31 // and the final number in each cell is then 
32 // @f[
33 //   h_i c \lange m\rangle 
34 // @f] 
35 // where @f$h_i@f$ is the number of hits in the cell @f$i@f$ 
36 // 
37 //
38
39 //____________________________________________________________________
40 AliPoissonCalculator::AliPoissonCalculator()
41   : TNamed(),
42     fEtaLumping(32), 
43     fPhiLumping(4), 
44     fTotal(0), 
45     fEmpty(0), 
46     fBasic(0),
47     fEmptyVsTotal(0),
48     fMean(0), 
49     fOcc(0),
50     fCorr(0),
51     fTotalList(),
52     fEmptyList(),
53     fRunningAverage(false)
54 {
55   //
56   // CTOR
57   // 
58 }
59
60 //____________________________________________________________________
61 AliPoissonCalculator::AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/)
62   : TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
63     fEtaLumping(32), 
64     fPhiLumping(4), 
65     fTotal(0), 
66     fEmpty(0), 
67     fBasic(0),
68     fEmptyVsTotal(0),
69     fMean(0), 
70     fOcc(0),
71     fCorr(0),
72     fTotalList(),
73     fEmptyList(),
74     fRunningAverage(false)
75 {
76   //
77   // CTOR
78   //
79   fEmptyList.SetOwner();
80   fTotalList.SetOwner();
81   
82
83 }
84 //____________________________________________________________________
85 AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
86   : TNamed(o),
87     fEtaLumping(o.fEtaLumping),
88     fPhiLumping(o.fPhiLumping),
89     fTotal(0), 
90     fEmpty(0),
91     fBasic(0), 
92     fEmptyVsTotal(0),
93     fMean(0), 
94     fOcc(0),
95     fCorr(0),
96     fTotalList(),
97     fEmptyList(),
98     fRunningAverage(o.fRunningAverage)
99 {
100   Init();
101   Reset(o.fBasic);
102 }
103
104 //____________________________________________________________________
105 AliPoissonCalculator::~AliPoissonCalculator()
106 {
107   // CleanUp();
108 }
109
110 //____________________________________________________________________
111 void
112 AliPoissonCalculator::CleanUp()
113 {
114   if (fTotal)        delete fTotal;        fTotal        = 0;
115   if (fEmpty)        delete fEmpty;        fEmpty        = 0;
116   if (fBasic)        delete fBasic;        fBasic        = 0;
117   if (fEmptyVsTotal) delete fEmptyVsTotal; fEmptyVsTotal = 0;
118   if (fMean)         delete fMean;         fMean         = 0;
119   if (fOcc)          delete fOcc;          fOcc          = 0;
120   if (fCorr)         delete fCorr;         fCorr         = 0;
121 }
122 //____________________________________________________________________
123 AliPoissonCalculator&
124 AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
125 {
126   TNamed::operator=(o);
127   fEtaLumping = o.fEtaLumping;
128   fPhiLumping = o.fPhiLumping;
129   fRunningAverage = o.fRunningAverage;
130   CleanUp();
131   Init();
132   Reset(o.fBasic);
133   return *this;
134 }
135
136 //____________________________________________________________________
137 void
138 AliPoissonCalculator::Init(UShort_t d, Char_t r, Int_t etaLumping, Int_t phiLumping)
139 {
140   // 
141   // Initialize 
142   // 
143   if (etaLumping > 0) SetEtaLumping(etaLumping);
144   if (phiLumping > 0) SetPhiLumping(phiLumping);
145   if(d > 0) {
146
147     Int_t    nEtaF   = (r == 'I' ? 512 : 256);
148     Int_t    nEta    = nEtaF / fEtaLumping;
149     Double_t etaMin  = -0.5;
150     Double_t etaMax  = nEtaF-0.5 ;
151     Int_t    nPhiF   = (r == 'I' ? 20 : 40);
152     Int_t    nPhi    = nPhiF / fPhiLumping;
153     Double_t phiMin  = -0.5;
154     Double_t phiMax  = nPhiF - 0.5;
155     
156     fBasic = new TH2D("basic", "Basic number of hits",
157                       nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
158     fBasic->SetDirectory(0);
159     fBasic->SetXTitle("#eta");
160     fBasic->SetYTitle("#varphi [radians]");
161     fBasic->Sumw2();
162     
163     for(Int_t v = 1 ; v < 11; v++)  { //CHC bins
164       for(Int_t centbin = 0 ; centbin < 13; centbin++) {
165         
166         TH2D* hTotal = new TH2D(Form("totalFMD%d%c_vertex%d_cent%d",d,r,v,centbin),"Total number of bins/region",
167                               nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
168         TH2D* hEmpty = new TH2D(Form("emptyFMD%d%c_vertex%d_cent%d",d,r,v,centbin), "Empty number of bins/region",
169                                 nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
170         hEmpty->Sumw2();
171         hTotal->Sumw2();
172         fEmptyList.Add(hEmpty);
173         fTotalList.Add(hTotal);
174         
175       }
176     }
177   }
178   //Create diagnostics if void
179   if (fEmptyVsTotal) return;
180   
181   Int_t n = fEtaLumping * fPhiLumping + 1;
182   fEmptyVsTotal = new TH2D("emptyVsTotal", 
183                            "# of empty # bins vs total # bins", 
184                            n, -.5, n-.5, n, -.5, n-.5);
185   fEmptyVsTotal->SetXTitle("Total # bins");
186   fEmptyVsTotal->SetYTitle("# empty bins");
187   fEmptyVsTotal->SetZTitle("Correlation");
188   fEmptyVsTotal->SetOption("colz");
189   fEmptyVsTotal->SetDirectory(0);
190
191   n = (fEtaLumping + fPhiLumping);
192   fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
193                    10*n+1, -.1, n+.1);
194   fMean->SetXTitle("#bar{N_{ch}}=#log(empty/total)");
195   fMean->SetYTitle("Events");
196   fMean->SetFillColor(kRed+1);
197   fMean->SetFillStyle(3001);
198   fMean->SetLineColor(kBlack);
199   fMean->SetDirectory(0);
200
201   fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
202                   1000, 0, 100);
203   fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
204   fOcc->SetYTitle("Events");
205   fOcc->SetFillColor(kBlue+1);
206   fOcc->SetFillStyle(3001);
207   fOcc->SetLineColor(kBlack);
208   fOcc->SetDirectory(0);
209
210   fCorr = new TH2D("correction", "Correction as function of mean N_{ch}", 
211                    10*n+1, -.1, n+.1, 100, 0, 10);
212   fCorr->SetXTitle("#bar{N_{ch}}");
213   fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})");
214   fCorr->SetZTitle("Events");
215   fCorr->SetOption("colz");
216   fCorr->SetDirectory(0);
217   
218   
219 }
220 //____________________________________________________________________
221 void AliPoissonCalculator::SetObject(UShort_t d, Char_t r, UShort_t v, Double_t cent) {
222   
223   Int_t centbin = 0;
224   if(cent > 0) {
225     if(cent > 0 && cent <5) centbin = 1;
226     if(cent > 5 && cent <10) centbin = 2;
227     else if (cent>10) centbin = (Int_t)(cent/10.) + 2;
228   }
229   
230   fTotal = static_cast<TH2D*>(fTotalList.FindObject(Form("totalFMD%d%c_vertex%d_cent%d",d,r,v,centbin)));
231   fEmpty = static_cast<TH2D*>(fEmptyList.FindObject(Form("emptyFMD%d%c_vertex%d_cent%d",d,r,v,centbin)));
232   
233   return;
234   
235 }
236 //____________________________________________________________________
237 void
238 AliPoissonCalculator::Output(TList* d)
239 {
240   if (!d) return;
241   if (!fEmptyVsTotal) Init();
242   d->Add(fEmptyVsTotal);
243   d->Add(fMean);
244   d->Add(fOcc);
245   d->Add(fCorr);
246 }
247
248 //____________________________________________________________________
249 void
250 AliPoissonCalculator::Reset(const TH2D* base)
251 {
252   // 
253   // Reset histogram 
254   // 
255   if (!base) return;
256   if (fBasic /* && fTotal && fEmpty*/) {
257     fBasic->Reset();
258     if(!fRunningAverage) {
259       fTotal->Reset();
260       fEmpty->Reset();
261     }
262     return;
263   }
264   /*  
265   Int_t    nEtaF   = base->GetNbinsX();
266   Int_t    nEta    = nEtaF / fEtaLumping;
267   Double_t etaMin  = base->GetXaxis()->GetXmin();
268   Double_t etaMax  = base->GetXaxis()->GetXmax();
269   Int_t    nPhiF   = base->GetNbinsY();
270   Int_t    nPhi    = nPhiF / fPhiLumping;
271   Double_t phiMin  = base->GetYaxis()->GetXmin();
272   Double_t phiMax  = base->GetYaxis()->GetXmax();
273   
274
275   
276   
277   
278   //fTotal = new TH2D("total", "Total number of bins/region",
279   //                nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
280   //fEmpty = new TH2D("empty", "Empty number of bins/region",
281   //nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
282   fBasic = new TH2D("basic", "Basic number of hits",
283                     nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
284     
285   //fTotal->SetDirectory(0);
286   //fEmpty->SetDirectory(0);
287   fBasic->SetDirectory(0);
288   //fTotal->SetXTitle("#eta");
289   //fEmpty->SetXTitle("#eta");
290   fBasic->SetXTitle("#eta");
291   //fTotal->SetYTitle("#varphi [radians]");
292   //fEmpty->SetYTitle("#varphi [radians]");
293   fBasic->SetYTitle("#varphi [radians]");
294   //fTotal->Sumw2();
295   //fEmpty->Sumw2();
296   fBasic->Sumw2();
297   */
298
299 }
300
301 //____________________________________________________________________
302 void
303 AliPoissonCalculator::Fill(UShort_t strip, UShort_t sec, Bool_t hit, 
304                            Double_t weight)
305 {
306   // 
307   // Fill in an observation 
308   // 
309   // Parameters:
310   //    eta     Eta value 
311   //    phi     Phi value
312   //    hit     True if hit 
313   //    weight  Weight if this 
314   //
315   
316   fTotal->Fill(strip, sec);
317   if (hit) fBasic->Fill(strip, sec, weight);
318   else     fEmpty->Fill(strip, sec);
319 }
320
321 //____________________________________________________________________
322 Double_t 
323 AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
324 {
325   if (total <= 0) return 0;
326   if (empty < .001) empty = .001;
327   return -TMath::Log(empty/total);
328 }
329 //____________________________________________________________________
330 Double_t 
331 AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
332 {
333   if (total <= 0) return 0;
334   if (TMath::Abs(empty-total) < .001) empty = total - .001;
335   return 1 / (1 - empty / total);
336 }
337
338 //____________________________________________________________________
339 TH2D*
340 AliPoissonCalculator::Result()
341 {
342   // 
343   // Calculate result and store in @a output
344   // 
345   // Return:
346   //    The result histogram (fBase overwritten)
347   //
348   
349   // Double_t total = fEtaLumping * fPhiLumping;
350   
351   for (Int_t ieta = 1; ieta <= fBasic->GetNbinsX(); ieta++) { 
352     for (Int_t iphi = 1; iphi <= fBasic->GetNbinsY(); iphi++) { 
353       Double_t eta      = fBasic->GetXaxis()->GetBinCenter(ieta);
354       Double_t phi      = fBasic->GetYaxis()->GetBinCenter(iphi);
355       Int_t    jeta     = fEmpty->GetXaxis()->FindBin(eta);
356       Int_t    jphi     = fEmpty->GetYaxis()->FindBin(phi);
357       Double_t empty    = fEmpty->GetBinContent(jeta, jphi);
358       Double_t total    = fTotal->GetBinContent(jeta, jphi);
359       Double_t hits     = fBasic->GetBinContent(ieta,iphi);
360       // Mean in region of interest 
361       Double_t poissonM = CalculateMean(empty, total);
362       Double_t poissonC = CalculateCorrection(empty, total);
363       
364       Double_t poissonV = hits * poissonM * poissonC;
365       Double_t poissonE = TMath::Sqrt(poissonV);
366       if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
367           
368       fBasic->SetBinContent(ieta,iphi,poissonV);
369       fBasic->SetBinError(ieta,iphi,poissonE);
370     }
371   }
372   for (Int_t ieta = 1; ieta <= fEmpty->GetNbinsX(); ieta++) { 
373     for (Int_t iphi = 1; iphi <= fEmpty->GetNbinsY(); iphi++) { 
374       Double_t empty    = fEmpty->GetBinContent(ieta, iphi);
375       Double_t total    = fTotal->GetBinContent(ieta, iphi);
376       Double_t mean     = CalculateMean(empty, total);
377       Double_t corr     = CalculateCorrection(empty, total);
378       fEmptyVsTotal->Fill(total, empty);
379       fMean->Fill(mean);
380       fOcc->Fill(100 * (1 - empty/total));
381       //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
382       fCorr->Fill(mean, corr);
383     }
384   }
385   return fBasic;
386 }
387   
388 //____________________________________________________________________
389 void
390 AliPoissonCalculator::Print(const Option_t*) const
391 {
392   // 
393   // Print information 
394   // 
395   // Parameters:
396   //    option Not used
397   //
398   char ind[gROOT->GetDirLevel()+3];
399   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
400   ind[gROOT->GetDirLevel()] = '\0';
401   std::cout << ind << ClassName() << ": " << GetName() << '\n'
402             << std::boolalpha 
403             << ind << " Eta lumping:            " << fEtaLumping << '\n'
404             << ind << " Phi lumping:            " << fPhiLumping << '\n'
405             << std::noboolalpha << std::endl;
406 }
407 //____________________________________________________________________
408 void
409 AliPoissonCalculator::Browse(TBrowser* b)
410 {
411   // 
412   // Browse this object
413   // 
414   // Parameters:
415   //    b Object to browse 
416   //
417   if (fTotal)        b->Add(fTotal);
418   if (fEmpty)        b->Add(fEmpty);
419   if (fBasic)        b->Add(fBasic);
420   if (fEmptyVsTotal) b->Add(fEmptyVsTotal);
421   if (fMean)         b->Add(fMean);
422   if (fOcc)          b->Add(fOcc);
423   if (fCorr)         b->Add(fCorr);
424 }
425 // 
426 // EOF
427 //