PWGLF/FORWARD/analysis2/AliPoissonCalculator.h

   1 #ifndef ALIPOISSONCALCULATOR_H
   2 #define ALIPOISSONCALCULATOR_H
   3 #include <TNamed.h>
   4 class TH2D;
   5 class TH1D;
   6 class TBrowser;
   7 class TAxis;
   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  * 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
  20  *
  21  * @f[
  22  * \langle 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{-\langle m\rangle}}\\ &=&
  29  *          \frac{1}{1 - \frac{e}{t}}
  30  * @f}
  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$
  34  *
  35  */
  36 class AliPoissonCalculator : public TNamed
  37 {
  38 public:
  39   /**
  40    * Constructor
  41    */
  42   AliPoissonCalculator();
  43   /**
  44    * Constructor
  45    *
  46    */
  47   AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/);
  48   /**
  49    * Copy constructor
  50    *
  51    * @param o Object to copy from
  52    */
  53   AliPoissonCalculator(const AliPoissonCalculator& o);
  54
  55   /**
  56    * Destructor
  57    */
  58   virtual ~AliPoissonCalculator();
  59   /**
  60    * Assignment operator
  61    *
  62    * @param o Object to assign from
  63    *
  64    * @return Reference to this object
  65    */
  66   AliPoissonCalculator& operator=(const AliPoissonCalculator& o);
  67   /**
  68    * Set the number of eta bins to group into a region
  69    *
  70    * @param nx Number of @f$\eta@f$ bins per region
  71    * @param ny Number of @f$\phi@f$ bins per region
  72    */
  73   void SetLumping(UShort_t nx, UShort_t ny);
  74   /**
  75    * Set the number of X bins to group into a region
  76    *
  77    * @param nx Number of eta bins per region
  78    */
  79   void SetXLumping(UShort_t nx) { SetLumping(nx, fYLumping); } //*MENU*
  80   /**
  81    * Set the number of Y bins to group into a region
  82    *
  83    * @param ny Number of eta bins per region
  84    */
  85   void SetYLumping(UShort_t ny) { SetLumping(fYLumping, ny); } //*MENU*
  86   /**
  87    * Intialize this object
  88    *
  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
  91    */
  92   void Init(Int_t xLumping=-1, Int_t yLumping=-1);
  93
  94   /**
  95    * Initialize this object.
  96    *
  97    * Also book the cache histograms
  98    *
  99    * @param xaxis The X-axis
 100    * @param yaxis The Y-axis
 101    */
 102   void Define(const TAxis& xaxis, const TAxis& yaxis);
 103   /**
 104    * Make output stuff for the passed list
 105    *
 106    */
 107   void MakeOutput();
 108   /**
 109    * Output stuff to the passed list
 110    *
 111    * @param d List to add output histograms to
 112    */
 113   void Output(TList* d);
 114   /**
 115    * Reset the cache histogram
 116    *
 117    * @param base Base histogram
 118    */
 119   void Reset(const TH2D* base);
 120   /**
 121    * Fill in an observation
 122    *
 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
 127    */
 128   void Fill(UShort_t strip, UShort_t sec, Bool_t hit, Double_t weight=1);
 129   /**
 130    * Calculate result and store in @a output
 131    *
 132    * @param correct Whether to apply correction or not
 133    *
 134    * @return The result histogram (fBase overwritten)
 135    */
 136   TH2D* Result(Bool_t correct=true);
 137   /**
 138    * After calculating the results, fill the diagnostics histograms
 139    *
 140    */
 141   void FillDiagnostics();
 142   /**
 143    * @return Always true
 144    */
 145   Bool_t IsFolder() const { return kTRUE; }
 146   /**
 147    * Print information
 148    *
 149    * @param option Not used
 150    */
 151   void Print(const Option_t* option="") const;
 152   /**
 153    * Browse this object
 154    *
 155    * @param b Object to browse
 156    */
 157   void Browse(TBrowser* b);
 158
 159   /**
 160    * Get the empty versus total histogram
 161    *
 162    * @return Empty versus total
 163    */
 164   TH2D* GetEmptyVsTotal() const { return fEmptyVsTotal; }
 165   /**
 166    * Get the histogram of the means
 167    *
 168    * @return Means
 169    */
 170   TH1D* GetMean() const { return fMean; }
 171   /**
 172    * Get the occupancy histogram
 173    *
 174    * @return Occupancy histogram
 175    */
 176   TH1D* GetOccupancy() const { return fOcc; }
 177   /**
 178    * Get the correction histogram
 179    *
 180    * @return correction histogram
 181    */
 182   TH2D* GetCorrection() const { return fCorr; }
 183
 184   /**
 185    * Get the X bin in the reduced historgam
 186    *
 187    * @param ix X bin in full histogram
 188    *
 189    * @return X bin in reduced histogram
 190    */
 191   Int_t GetReducedXBin(Int_t ix) const;
 192   /**
 193    * Get the X bin in the reduced historgam
 194    *
 195    * @param x X value
 196    *
 197    * @return X bin in reduced histogram
 198    */
 199   Int_t GetReducedXBin(Double_t x) const;
 200   /**
 201    * Get the Y bin in the reduced historgam
 202    *
 203    * @param iy Y bin in full histogram
 204    *
 205    * @return Y bin in reduced histogram
 206    */
 207   Int_t GetReducedYBin(Int_t iy) const;
 208   /**
 209    * Get the Y bin in the reduced historgam
 210    *
 211    * @param y Y value
 212    *
 213    * @return Y bin in reduced histogram
 214    */
 215   Int_t GetReducedYBin(Double_t y) const;
 216
 217 protected:
 218   /**
 219    * check that the lumping parameter makes sense
 220    *
 221    * @param which   Which axis
 222    * @param nBins   Number of bins
 223    * @param lumping Lumping
 224    *
 225    * @return The new value of the lumping
 226    */
 227   Int_t CheckLumping(char which, Int_t nBins, Int_t lumping) const;
 228   /**
 229    * Clean up allocated space
 230    *
 231    */
 232   void CleanUp();
 233   /**
 234    * Calculate the mean
 235    *
 236    * This is based on the fact that for a Poisson
 237    * @f[
 238    *   P(n;\lambda) = \frac{-\lambda^n e^{-\lambda}}{n!}
 239    * @f]
 240    * we have the probability for 0 observation
 241    * @f[
 242    *   P(0;\lambda) = e^{-\lambda} = \frac{N_{empty}}{N_{total}}
 243    * @f]
 244    * and so we get that the mean is the defined region is
 245    * @f[
 246    *   \lambda = -\log\left(\frac{N_{empty}}{N_{total}}\right)
 247    * @f]
 248    *
 249    * Note the boundary conditions
 250    * - @f$N_{total}=0 \rightarrow\lambda=0@f$
 251    * - @f$N_{empty}<\epsilon\rightarrow N_{empty} = \epsilon@f$
 252    *
 253    * @param empty Number of empty bins
 254    * @param total Total number of bins
 255    *
 256    * @return The mean in the defined region
 257    */
 258   Double_t CalculateMean(Double_t empty, Double_t total) const;
 259   /**
 260    * The mean @f$\lambda@f$ calculated above is not the full story.
 261    * In addition it needs to be corrected using the expression
 262    * @f[
 263    *   \frac{1}{1-e^{\lambda}} =
 264    *     \frac{1}{1-\frac{N_{empty}}{N_{total}}}
 265    * @f]
 266    *
 267    * Note the boundary conditions
 268    * - @f$N_{total}=0 \rightarrow\lambda=0@f$
 269    * - @f$|N_{total}-N_{empty}|<\epsilon\rightarrow N_{empty} =
 270    * N_{total}-\epsilon@f$
 271    *
 272    * @param empty Number of empty bins
 273    * @param total Total number of bins
 274    *
 275    * @return The correction to the mean.
 276    */
 277   Double_t CalculateCorrection(Double_t empty, Double_t total) const;
 278   UShort_t fXLumping;   // Grouping of eta bins
 279   UShort_t fYLumping;   // Grouping of phi bins
 280   TH2D*    fTotal;        // Total number of strips in a region
 281   TH2D*    fEmpty;        // Total number of strips in a region
 282   TH2D*    fBasic;        // Total number basic hits in a region
 283   TH2D*    fEmptyVsTotal; // Empty versus total cells
 284   TH1D*    fMean;         // Mean calculated by poisson method
 285   TH1D*    fOcc;          // Histogram of occupancies
 286   TH2D*    fCorr;         // Correction as a function of mean
 287   ClassDef(AliPoissonCalculator,3) // Calculate N_ch using Poisson
 288 };
 289
 290 #endif
 291 // Local Variables:
 292 //   mode: C++
 293 // End:
 294