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
   8 /**
   9  * A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
  10  * using Poisson statistics.
  11  *
  12  * The input is assumed to be binned in @f$(\eta,\varphi)@f$ as
  13  * described by the 2D histogram passwd to the Reset member function.
  14  *
  15  * The data is grouped in to regions as defined by the parameters
  16  * fXLumping and fYLumping.  The total number of cells and number
  17  * of empty cells is then calculate in each region.  The mean
  18  * multiplicity over the region is then determined as
  19  *
  20  * @f[
  21  * \langle m\rangle = -\log\left(\frac{e}{t}\right)
  22  * @f]
  23  * where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
  24  * total number of cells in the region.  A correction for counting
  25  * statistics, is then applied
  26  * @f{eqnarray*}{
  27  *    c &=& \frac{1}{1 - \exp{-\langle m\rangle}}\\ &=&
  28  *          \frac{1}{1 - \frac{e}{t}}
  29  * @f}
  30  * and the final number in each cell is then
  31  * @f$h_i c \langle m\rangle@f$
  32  * where @f$h_i@f$ is the number of hits in the cell @f$i@f$
  33  *
  34  */
  35 class AliPoissonCalculator : public TNamed
  36 {
  37 public:
  38   /**
  39    * Constructor
  40    */
  41   AliPoissonCalculator();
  42   /**
  43    * Constructor
  44    *
  45    */
  46   AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/);
  47   /**
  48    * Copy constructor
  49    *
  50    * @param o Object to copy from
  51    */
  52   AliPoissonCalculator(const AliPoissonCalculator& o);
  53
  54   /**
  55    * Destructor
  56    */
  57   virtual ~AliPoissonCalculator();
  58   /**
  59    * Assignment operator
  60    *
  61    * @param o Object to assign from
  62    *
  63    * @return Reference to this object
  64    */
  65   AliPoissonCalculator& operator=(const AliPoissonCalculator& o);
  66   /**
  67    * Set the number of eta bins to group into a region
  68    *
  69    * @param n Number of eta bins per region
  70    */
  71   void SetLumping(UShort_t nx, UShort_t ny);
  72   /**
  73    * Set the number of X bins to group into a region
  74    *
  75    * @param n Number of eta bins per region
  76    */
  77   void SetXLumping(UShort_t nx) { SetLumping(nx, fYLumping); } //*MENU*
  78   /**
  79    * Set the number of Y bins to group into a region
  80    *
  81    * @param n Number of eta bins per region
  82    */
  83   void SetYLumping(UShort_t ny) { SetLumping(fYLumping, ny); } //*MENU*
  84   /**
  85    * Intialize this object
  86    *
  87    * @param xLumping If larger than 0, set the eta lumping to this
  88    * @param yLumping If larger than 0, set the phi lumping to this
  89    */
  90   void Init(Int_t xLumping=-1, Int_t yLumping=-1);
  91   /**
  92    * Make output stuff for the passed list
  93    *
  94    * @param none
  95    */
  96   void MakeOutput();
  97   /**
  98    * Output stuff to the passed list
  99    *
 100    * @param d List to add output histograms to
 101    */
 102   void Output(TList* d);
 103   /**
 104    * Reset the cache histogram
 105    *
 106    * @param xChannels number of channels in x
 107    * @param yChannels number of channels in y
 108    */
 109   void Reset(UShort_t xChannels, UShort_t yChannels);
 110   /**
 111    * Fill in an observation
 112    *
 113    * @param eta     Eta value
 114    * @param phi     Phi value
 115    * @param hit     True if hit
 116    * @param weight  Weight if this
 117    */
 118   void Fill(UShort_t strip, UShort_t sec, Bool_t hit, Double_t weight=1);
 119   /**
 120    * Calculate result and store in @a output
 121    *
 122    * @return The result histogram (fBase overwritten)
 123    */
 124   TH2D* Result(Bool_t correct=true);
 125   /**
 126    * @return Always true
 127    */
 128   Bool_t IsFolder() const { return kTRUE; }
 129   /**
 130    * Print information
 131    *
 132    * @param option Not used
 133    */
 134   void Print(const Option_t* option="") const;
 135   /**
 136    * Browse this object
 137    *
 138    * @param b Object to browse
 139    */
 140   void Browse(TBrowser* b);
 141
 142   /**
 143    * Get the empty versus total histogram
 144    *
 145    * @return Empty versus total
 146    */
 147   TH2D* GetEmptyVsTotal() const { return fEmptyVsTotal; }
 148   /**
 149    * Get the histogram of the means
 150    *
 151    * @return Means
 152    */
 153   TH1D* GetMean() const { return fMean; }
 154   /**
 155    * Get the occupancy histogram
 156    *
 157    * @return Occupancy histogram
 158    */
 159   TH1D* GetOccupancy() const { return fOcc; }
 160   /**
 161    * Get the correction histogram
 162    *
 163    * @return correction histogram
 164    */
 165   TH2D* GetCorrection() const { return fCorr; }
 166
 167   /**
 168    * Get the X bin in the reduced historgam
 169    *
 170    * @param ix X bin in full histogram
 171    *
 172    * @return X bin in reduced histogram
 173    */
 174   Int_t GetReducedXBin(Int_t ix) const;
 175   /**
 176    * Get the X bin in the reduced historgam
 177    *
 178    * @param x X value
 179    *
 180    * @return X bin in reduced histogram
 181    */
 182   Int_t GetReducedXBin(Double_t x) const;
 183   /**
 184    * Get the Y bin in the reduced historgam
 185    *
 186    * @param iy Y bin in full histogram
 187    *
 188    * @return Y bin in reduced histogram
 189    */
 190   Int_t GetReducedYBin(Int_t iy) const;
 191   /**
 192    * Get the Y bin in the reduced historgam
 193    *
 194    * @param y Y value
 195    *
 196    * @return Y bin in reduced histogram
 197    */
 198   Int_t GetReducedYBin(Double_t y) const;
 199
 200 protected:
 201   /**
 202    * Clean up allocated space
 203    *
 204    */
 205   void CleanUp();
 206   /**
 207    * Calculate the mean
 208    *
 209    * This is based on the fact that for a Poisson
 210    * @f[
 211    *   P(n;\lambda) = \frac{-\lambda^n e^{-\lambda}}{n!}
 212    * @f]
 213    * we have the probability for 0 observation
 214    * @f[
 215    *   P(0;\lambda) = e^{-\lambda} = \frac{N_{empty}}{N_{total}}
 216    * @f]
 217    * and so we get that the mean is the defined region is
 218    * @f[
 219    *   \lambda = -\log\left(\frac{N_{empty}}{N_{total}}\right)
 220    * @f]
 221    *
 222    * Note the boundary conditions
 223    * - @f$N_{total}=0 \rightarrow\lambda=0@f$
 224    * - @f$N_{empty}<\epsilon\rightarrow N_{empty} = \epsilon@f$
 225    *
 226    * @param empty Number of empty bins
 227    * @param total Total number of bins
 228    *
 229    * @return The mean in the defined region
 230    */
 231   Double_t CalculateMean(Double_t empty, Double_t total) const;
 232   /**
 233    * The mean @f$\lambda@f$ calculated above is not the full story.
 234    * In addition it needs to be corrected using the expression
 235    * @f[
 236    *   \frac{1}{1-e^{\lambda}} =
 237    *     \frac{1}{1-\frac{N_{empty}}{N_{total}}}
 238    * @f]
 239    *
 240    * Note the boundary conditions
 241    * - @f$N_{total}=0 \rightarrow\lambda=0@f$
 242    * - @f$|N_{total}-N_{empty}|<\epsilon\rightarrow N_{empty} =
 243    * N_{total}-\epsilon@f$
 244    *
 245    * @param empty Number of empty bins
 246    * @param total Total number of bins
 247    *
 248    * @return The correction to the mean.
 249    */
 250   Double_t CalculateCorrection(Double_t empty, Double_t total) const;
 251   UShort_t fXLumping;   // Grouping of eta bins
 252   UShort_t fYLumping;   // Grouping of phi bins
 253   TH2D*    fTotal;        // Total number of strips in a region
 254   TH2D*    fEmpty;        // Total number of strips in a region
 255   TH2D*    fBasic;        // Total number basic hits in a region
 256   TH2D*    fEmptyVsTotal; // Empty versus total cells
 257   TH1D*    fMean;         // Mean calculated by poisson method
 258   TH1D*    fOcc;          // Histogram of occupancies
 259   TH2D*    fCorr;         // Correction as a function of mean
 260   ClassDef(AliPoissonCalculator,2) // Calculate N_ch using Poisson
 261 };
 262
 263 #endif
 264 // Local Variables:
 265 //   mode: C++
 266 // End:
 267