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