More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.h
index 5a6a40cfa15df492ca58032646b928b73aca1b71..938678b42be9adc2975ac4a0f747408b6619a04a 100644 (file)
@@ -1,8 +1,14 @@
 #ifndef ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
 #define ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
 #include <TObject.h>
+#include <TString.h>
+#include <TObjArray.h>
 class TH2D;
+class TH1I;
+class TH1;
+class TF1;
 class TAxis;
+class AliESDEvent;
 
 /** 
  * Utilities used in the forward multiplcity analysis 
@@ -11,7 +17,362 @@ class TAxis;
  */
 class AliForwardUtil : public TObject
 {
- public:
+public:
+  //==================================================================
+  /** 
+   * @{ 
+   * @nane Collision/run parameters 
+   */
+  /**                                          
+   * Defined collision types 
+   */
+  enum ECollisionSystem {
+    kUnknown, 
+    kPP, 
+    kPbPb
+  };
+  //__________________________________________________________________
+  /** 
+   * Parse a collision system spec given in a string.   Known values are 
+   * 
+   *  - "pp", "p-p" which returns kPP 
+   *  - "PbPb", "Pb-Pb", "A-A", which returns kPbPb 
+   *  - Everything else gives kUnknown 
+   * 
+   * @param sys Collision system spec 
+   * 
+   * @return Collision system id 
+   */
+  static UShort_t ParseCollisionSystem(const char* sys);
+  /** 
+   * Get a string representation of the collision system 
+   * 
+   * @param sys  Collision system 
+   * - kPP -> "pp"
+   * - kPbPb -> "PbPb" 
+   * - anything else gives "unknown"
+   * 
+   * @return String representation of the collision system 
+   */
+  static const char* CollisionSystemString(UShort_t sys);
+  //__________________________________________________________________
+  /** 
+   * Parse the center of mass energy given as a float and return known 
+   * values as a unsigned integer
+   * 
+   * @param sys  Collision system (needed for AA)
+   * @param cms  Center of mass energy * total charge 
+   * 
+   * @return Center of mass energy per nucleon
+   */
+  static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
+  /** 
+   * Get a string representation of the center of mass energy per nuclean
+   * 
+   * @param sys  Collision system 
+   * @param sNN  Center of mass energy per nucleon
+   * 
+   * @return String representation of the center of mass energy per nuclean
+   */
+  static const char* CenterOfMassEnergyString(UShort_t cms);
+  //__________________________________________________________________
+  /** 
+   * Parse the magnetic field (in kG) as given by a floating point number
+   * 
+   * @param field  Magnetic field in kG 
+   * 
+   * @return Short integer value of magnetic field in kG 
+   */
+  static Short_t ParseMagneticField(Float_t field);
+  /** 
+   * Get a string representation of the magnetic field
+   * 
+   * @param field Magnetic field in kG
+   * 
+   * @return String representation of the magnetic field
+   */
+  static const char* MagneticFieldString(Short_t field);
+  /* @} */
+
+  /** 
+   * @{ 
+   * @name Energy stragling functions 
+   */
+  //__________________________________________________________________
+  /**
+   * Number of steps to do in the Landau, Gaussiam convolution 
+   */
+  static Int_t fgConvolutionSteps;
+  //------------------------------------------------------------------
+  /** 
+   * How many sigma's of the Gaussian in the Landau, Gaussian
+   * convolution to integrate over
+   */
+  static Double_t fgConvolutionNSigma;
+  //------------------------------------------------------------------
+  /** 
+   * Calculate the shifted Landau
+   * @f[
+   *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
+   * @f]
+   *
+   * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
+   * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
+   * @f$\Delta=0,\xi=1@f$. 
+   *
+   * @param x      Where to evaluate @f$ f'_{L}@f$ 
+   * @param delta  Most probable value 
+   * @param xi     The 'width' of the distribution 
+   *
+   * @return @f$ f'_{L}(x;\Delta,\xi) @f$
+   */
+  static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
+  
+  //------------------------------------------------------------------
+  /** 
+   * Calculate the value of a Landau convolved with a Gaussian 
+   * 
+   * @f[ 
+   * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
+   *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
+   *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
+   * @f]
+   * 
+   * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
+   * energy loss, @f$ \xi@f$ the width of the Landau, and 
+   * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
+   * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
+   * noise in the detector.  
+   *
+   * Note that this function uses the constants fgConvolutionSteps and
+   * fgConvolutionNSigma
+   * 
+   * References: 
+   *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
+   *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+   *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
+   * 
+   * @param x         where to evaluate @f$ f@f$
+   * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
+   * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
+   * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
+   * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
+   * 
+   * @return @f$ f@f$ evaluated at @f$ x@f$.  
+   */
+  static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                            Double_t sigma, Double_t sigma_n);
+
+  //------------------------------------------------------------------
+  /** 
+   * Evaluate 
+   * @f[ 
+   *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
+   * @f] 
+   * corresponding to @f$ i@f$ particles i.e., with the substitutions 
+   * @f[ 
+   *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
+   *    \xi       \rightarrow \xi_i       = i \xi\\
+   *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
+   *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
+   * @f] 
+   * 
+   * @param x        Where to evaluate 
+   * @param delta    @f$ \Delta@f$ 
+   * @param xi       @f$ \xi@f$ 
+   * @param sigma    @f$ \sigma@f$ 
+   * @param sigma_n  @f$ \sigma_n@f$
+   * @param i        @f$ i@f$
+   * 
+   * @return @f$ f_i@f$ evaluated
+   */  
+  static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                             Double_t sigma, Double_t sigma_n, Int_t i);
+
+  //------------------------------------------------------------------
+  /** 
+   * Numerically evaluate 
+   * @f[ 
+   *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
+   * @f] 
+   * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
+   * of the parameters is given by 
+   *
+   * - 0: @f$\Delta@f$ 
+   * - 1: @f$\xi@f$ 
+   * - 2: @f$\sigma@f$ 
+   * - 3: @f$\sigma_n@f$ 
+   *
+   * This is the partial derivative with respect to the parameter of
+   * the response function corresponding to @f$ i@f$ particles i.e.,
+   * with the substitutions
+   * @f[ 
+   *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
+   *    \xi       \rightarrow \xi_i       = i \xi\\
+   *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
+   *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
+   * @f] 
+   * 
+   * @param x        Where to evaluate 
+   * @param ipar     Parameter number 
+   * @param dp       @f$ \esilon\delta p_i@f$ for some value of @f$\epsilon@f$
+   * @param delta    @f$ \Delta@f$ 
+   * @param xi       @f$ \xi@f$ 
+   * @param sigma    @f$ \sigma@f$ 
+   * @param sigma_n  @f$ \sigma_n@f$
+   * @param i        @f$ i@f$
+   * 
+   * @return @f$ f_i@f$ evaluated
+   */  
+  static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
+                                  Double_t delta, Double_t xi, 
+                                  Double_t sigma, Double_t sigma_n, Int_t i);
+
+  //------------------------------------------------------------------
+  /** 
+   * Evaluate 
+   * @f[ 
+   *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
+   * @f] 
+   * 
+   * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
+   * Landau with a Gaussian (see LandauGaus).  Note that 
+   * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
+   * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
+   *  
+   * References: 
+   *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
+   *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+   *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
+   * 
+   * @param x        Where to evaluate @f$ f_N@f$
+   * @param delta    @f$ \Delta_1@f$ 
+   * @param xi       @f$ \xi_1@f$
+   * @param sigma    @f$ \sigma_1@f$ 
+   * @param sigma_n  @f$ \sigma_n@f$ 
+   * @param n        @f$ N@f$ in the sum above.
+   * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
+   *                 @f$ i > 1@f$ 
+   * 
+   * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
+   */
+  static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                             Double_t sigma, Double_t sigma_n, Int_t n, 
+                             Double_t* a);
+  /** 
+   * Generate a TF1 object of @f$ f_I@f$ 
+   * 
+   * @param c        Constant
+   * @param delta    @f$ \Delta@f$ 
+   * @param xi       @f$ \xi_1@f$             
+   * @param sigma    @f$ \sigma_1@f$          
+   * @param sigma_n  @f$ \sigma_n@f$          
+   * @param i       @f$ i@f$ - the number of particles
+   * @param xmin     Least value of range
+   * @param xmax     Largest value of range
+   * 
+   * @return Newly allocated TF1 object
+   */
+  static TF1* MakeILandauGaus(Double_t c, 
+                             Double_t delta, Double_t xi, 
+                             Double_t sigma, Double_t sigma_n,
+                             Int_t    i, 
+                             Double_t xmin,  Double_t  xmax);
+  /** 
+   * Generate a TF1 object of @f$ f_N@f$ 
+   * 
+   * @param c         Constant                        
+   * @param delta     @f$ \Delta@f$                   
+   * @param xi               @f$ \xi_1@f$                     
+   * @param sigma     @f$ \sigma_1@f$                 
+   * @param sigma_n   @f$ \sigma_n@f$                 
+   * @param n        @f$ N@f$ - how many particles to sum to
+   * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
+   *                  @f$ i > 1@f$ 
+   * @param xmin      Least value of range  
+   * @param xmax      Largest value of range
+   * 
+   * @return Newly allocated TF1 object
+   */
+  static TF1* MakeNLandauGaus(Double_t c, 
+                             Double_t delta, Double_t  xi, 
+                             Double_t sigma, Double_t  sigma_n,
+                             Int_t    n,     Double_t* a, 
+                             Double_t xmin,  Double_t  xmax);
+                                                   
+  //__________________________________________________________________
+  /** 
+   * Structure to do fits to the energy loss spectrum 
+   * 
+   */
+  struct ELossFitter 
+  {
+    enum { 
+      kC     = 0,
+      kDelta, 
+      kXi, 
+      kSigma, 
+      kSigmaN, 
+      kN, 
+      kA
+    };
+    /** 
+     * Constructor 
+     * 
+     * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
+     * @param maxRange   Maximum range to fit to 
+     * @param minusBins  The number of bins below maximum to use 
+     */
+    ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
+    virtual ~ELossFitter();
+    /** 
+     * Clear internal arrays 
+     * 
+     */
+    void Clear();
+    /** 
+     * Fit a 1-particle signal to the passed energy loss distribution 
+     * 
+     * Note that this function clears the internal arrays first 
+     *
+     * @param dist    Data to fit the function to 
+     * @param sigman If larger than zero, the initial guess of the
+     *               detector induced noise. If zero or less, then this 
+     *               parameter is ignored in the fit (fixed at 0)
+     * 
+     * @return The function fitted to the data 
+     */
+    TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
+    /** 
+     * Fit a N-particle signal to the passed energy loss distribution 
+     *
+     * If there's no 1-particle fit present, it does that first 
+     *
+     * @param dist   Data to fit the function to 
+     * @param n      Number of particle signals to fit 
+     * @param sigman If larger than zero, the initial guess of the
+     *               detector induced noise. If zero or less, then this 
+     *               parameter is ignored in the fit (fixed at 0)
+     * 
+     * @return The function fitted to the data 
+     */
+    TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
+     
+
+    const Double_t fLowCut;     // Lower cut on data 
+    const Double_t fMaxRange;   // Maximum range to fit 
+    const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
+    TObjArray fFitResults;      // Array of fit results 
+    TObjArray fFunctions;       // Array of functions 
+  };
+  /* @} */
+      
+
+  //==================================================================
+  /** 
+   * @{
+   * @name Convenience containers 
+   */
   /** 
    * Structure to hold histograms 
    *
@@ -85,8 +446,46 @@ class AliForwardUtil : public TObject
     TH2D* fFMD2o; // Histogram for FMD2o
     TH2D* fFMD3i; // Histogram for FMD3i
     TH2D* fFMD3o; // Histogram for FMD3o
+
+    ClassDef(Histos,1) 
   };
 
+  //__________________________________________________________________
+  struct RingHistos : public TObject
+  {
+    RingHistos() : fDet(0), fRing('\0'), fName("") {}
+    RingHistos(UShort_t d, Char_t r) 
+      : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
+    {}
+    RingHistos(const RingHistos& o) 
+      : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
+    {}
+    virtual ~RingHistos() {}
+    RingHistos& operator=(const RingHistos& o) 
+    {
+      TObject::operator=(o);
+      fDet  = o.fDet;
+      fRing = o.fRing;
+      fName = o.fName;
+      return *this;
+    }
+    TList* DefineOutputList(TList* d) const;
+    TList* GetOutputList(TList* d) const;
+    TH1* GetOutputHist(TList* d, const char* name) const;
+    Color_t Color() const 
+    { 
+      return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
+             + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
+    }
+
+    UShort_t fDet; 
+    Char_t   fRing;
+    TString  fName;
+
+    ClassDef(RingHistos,1) 
+  };
+  /* @} */
+    
 };
 
 #endif