1 #ifndef ALIPOISSONCALCULATOR_H
2 #define ALIPOISSONCALCULATOR_H
10 * A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
11 * using Poisson statistics.
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.
16 * The data is grouped in to regions as defined by the parameters
17 * fXLumping and fYLumping. 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
22 * \langle m\rangle = -\log\left(\frac{e}{t}\right)
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
28 * c &=& \frac{1}{1 - \exp{-\langle m\rangle}}\\ &=&
29 * \frac{1}{1 - \frac{e}{t}}
31 * and the final number in each cell is then
32 * @f$h_i c \langle m\rangle@f$
33 * where @f$h_i@f$ is the number of hits in the cell @f$i@f$
36 class AliPoissonCalculator : public TNamed
42 AliPoissonCalculator();
47 AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/);
51 * @param o Object to copy from
53 AliPoissonCalculator(const AliPoissonCalculator& o);
58 virtual ~AliPoissonCalculator();
62 * @param o Object to assign from
64 * @return Reference to this object
66 AliPoissonCalculator& operator=(const AliPoissonCalculator& o);
68 * Set the number of eta bins to group into a region
70 * @param nx Number of @f$\eta@f$ bins per region
71 * @param ny Number of @f$\phi@f$ bins per region
73 void SetLumping(UShort_t nx, UShort_t ny);
75 * Set the number of X bins to group into a region
77 * @param nx Number of eta bins per region
79 void SetXLumping(UShort_t nx) { SetLumping(nx, fYLumping); } //*MENU*
81 * Set the number of Y bins to group into a region
83 * @param ny Number of eta bins per region
85 void SetYLumping(UShort_t ny) { SetLumping(fYLumping, ny); } //*MENU*
87 * Intialize this object
89 * @param xLumping If larger than 0, set the eta lumping to this
90 * @param yLumping If larger than 0, set the phi lumping to this
92 void Init(Int_t xLumping=-1, Int_t yLumping=-1);
95 * Initialize this object.
97 * Also book the cache histograms
99 * @param xaxis The X-axis
100 * @param yaxis The Y-axis
102 void Define(const TAxis& xaxis, const TAxis& yaxis);
104 * Make output stuff for the passed list
109 * Output stuff to the passed list
111 * @param d List to add output histograms to
113 void Output(TList* d);
115 * Reset the cache histogram
117 * @param base Base histogram
119 void Reset(const TH2D* base);
121 * Fill in an observation
123 * @param strip X axis bin number
124 * @param sec Y axis bin number
125 * @param hit True if hit
126 * @param weight Weight if this
128 void Fill(UShort_t strip, UShort_t sec, Bool_t hit, Double_t weight=1);
130 * Calculate result and store in @a output
132 * @param correct Whether to apply correction or not
134 * @return The result histogram (fBase overwritten)
136 TH2D* Result(Bool_t correct=true);
138 * @return Always true
140 Bool_t IsFolder() const { return kTRUE; }
144 * @param option Not used
146 void Print(const Option_t* option="") const;
150 * @param b Object to browse
152 void Browse(TBrowser* b);
155 * Get the empty versus total histogram
157 * @return Empty versus total
159 TH2D* GetEmptyVsTotal() const { return fEmptyVsTotal; }
161 * Get the histogram of the means
165 TH1D* GetMean() const { return fMean; }
167 * Get the occupancy histogram
169 * @return Occupancy histogram
171 TH1D* GetOccupancy() const { return fOcc; }
173 * Get the correction histogram
175 * @return correction histogram
177 TH2D* GetCorrection() const { return fCorr; }
180 * Get the X bin in the reduced historgam
182 * @param ix X bin in full histogram
184 * @return X bin in reduced histogram
186 Int_t GetReducedXBin(Int_t ix) const;
188 * Get the X bin in the reduced historgam
192 * @return X bin in reduced histogram
194 Int_t GetReducedXBin(Double_t x) const;
196 * Get the Y bin in the reduced historgam
198 * @param iy Y bin in full histogram
200 * @return Y bin in reduced histogram
202 Int_t GetReducedYBin(Int_t iy) const;
204 * Get the Y bin in the reduced historgam
208 * @return Y bin in reduced histogram
210 Int_t GetReducedYBin(Double_t y) const;
214 * check that the lumping parameter makes sense
216 * @param which Which axis
217 * @param nBins Number of bins
218 * @param lumping Lumping
220 * @return The new value of the lumping
222 Int_t CheckLumping(char which, Int_t nBins, Int_t lumping) const;
224 * Clean up allocated space
231 * This is based on the fact that for a Poisson
233 * P(n;\lambda) = \frac{-\lambda^n e^{-\lambda}}{n!}
235 * we have the probability for 0 observation
237 * P(0;\lambda) = e^{-\lambda} = \frac{N_{empty}}{N_{total}}
239 * and so we get that the mean is the defined region is
241 * \lambda = -\log\left(\frac{N_{empty}}{N_{total}}\right)
244 * Note the boundary conditions
245 * - @f$N_{total}=0 \rightarrow\lambda=0@f$
246 * - @f$N_{empty}<\epsilon\rightarrow N_{empty} = \epsilon@f$
248 * @param empty Number of empty bins
249 * @param total Total number of bins
251 * @return The mean in the defined region
253 Double_t CalculateMean(Double_t empty, Double_t total) const;
255 * The mean @f$\lambda@f$ calculated above is not the full story.
256 * In addition it needs to be corrected using the expression
258 * \frac{1}{1-e^{\lambda}} =
259 * \frac{1}{1-\frac{N_{empty}}{N_{total}}}
262 * Note the boundary conditions
263 * - @f$N_{total}=0 \rightarrow\lambda=0@f$
264 * - @f$|N_{total}-N_{empty}|<\epsilon\rightarrow N_{empty} =
265 * N_{total}-\epsilon@f$
267 * @param empty Number of empty bins
268 * @param total Total number of bins
270 * @return The correction to the mean.
272 Double_t CalculateCorrection(Double_t empty, Double_t total) const;
273 UShort_t fXLumping; // Grouping of eta bins
274 UShort_t fYLumping; // Grouping of phi bins
275 TH2D* fTotal; // Total number of strips in a region
276 TH2D* fEmpty; // Total number of strips in a region
277 TH2D* fBasic; // Total number basic hits in a region
278 TH2D* fEmptyVsTotal; // Empty versus total cells
279 TH1D* fMean; // Mean calculated by poisson method
280 TH1D* fOcc; // Histogram of occupancies
281 TH2D* fCorr; // Correction as a function of mean
282 ClassDef(AliPoissonCalculator,2) // Calculate N_ch using Poisson