PWG2/FORWARD/analysis2/AliForwardUtil.h

   1 //
   2 // Utilities used in the forward multiplcity analysis
   3 //
   4 //
   5 #ifndef ALIFORWARDUTIL_H
   6 #define ALIFORWARDUTIL_H
   7 #include <TObject.h>
   8 #include <TString.h>
   9 #include <TObjArray.h>
  10 class TH2D;
  11 class TH1I;
  12 class TH1;
  13 class TF1;
  14 class TAxis;
  15 class AliESDEvent;
  16
  17 /**
  18  * Utilities used in the forward multiplcity analysis
  19  *
  20  * @ingroup pwg2_forward
  21  */
  22 class AliForwardUtil : public TObject
  23 {
  24 public:
  25   //==================================================================
  26   /**
  27    * @{
  28    * @name Collision/run parameters
  29    */
  30   /**
  31    * Defined collision types
  32    */
  33   enum ECollisionSystem {
  34     kUnknown,
  35     kPP,
  36     kPbPb
  37   };
  38   //__________________________________________________________________
  39   /**
  40    * Parse a collision system spec given in a string.   Known values are
  41    *
  42    *  - "pp", "p-p" which returns kPP
  43    *  - "PbPb", "Pb-Pb", "A-A", which returns kPbPb
  44    *  - Everything else gives kUnknown
  45    *
  46    * @param sys Collision system spec
  47    *
  48    * @return Collision system id
  49    */
  50   static UShort_t ParseCollisionSystem(const char* sys);
  51   /**
  52    * Get a string representation of the collision system
  53    *
  54    * @param sys  Collision system
  55    * - kPP -> "pp"
  56    * - kPbPb -> "PbPb"
  57    * - anything else gives "unknown"
  58    *
  59    * @return String representation of the collision system
  60    */
  61   static const char* CollisionSystemString(UShort_t sys);
  62   //__________________________________________________________________
  63   /**
  64    * Parse the center of mass energy given as a float and return known
  65    * values as a unsigned integer
  66    *
  67    * @param sys  Collision system (needed for AA)
  68    * @param cms  Center of mass energy * total charge
  69    *
  70    * @return Center of mass energy per nucleon
  71    */
  72   static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
  73   /**
  74    * Get a string representation of the center of mass energy per nuclean
  75    *
  76    * @param cms  Center of mass energy per nucleon
  77    *
  78    * @return String representation of the center of mass energy per nuclean
  79    */
  80   static const char* CenterOfMassEnergyString(UShort_t cms);
  81   //__________________________________________________________________
  82   /**
  83    * Parse the magnetic field (in kG) as given by a floating point number
  84    *
  85    * @param field  Magnetic field in kG
  86    *
  87    * @return Short integer value of magnetic field in kG
  88    */
  89   static Short_t ParseMagneticField(Float_t field);
  90   /**
  91    * Get a string representation of the magnetic field
  92    *
  93    * @param field Magnetic field in kG
  94    *
  95    * @return String representation of the magnetic field
  96    */
  97   static const char* MagneticFieldString(Short_t field);
  98   /* @} */
  99
 100   /**
 101    * @{
 102    * @name Energy stragling functions
 103    */
 104   //__________________________________________________________________
 105   /**
 106    * Number of steps to do in the Landau, Gaussiam convolution
 107    */
 108   static Int_t fgConvolutionSteps;
 109   //------------------------------------------------------------------
 110   /**
 111    * How many sigma's of the Gaussian in the Landau, Gaussian
 112    * convolution to integrate over
 113    */
 114   static Double_t fgConvolutionNSigma;
 115   //------------------------------------------------------------------
 116   /**
 117    * Calculate the shifted Landau
 118    * @f[
 119    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 120    * @f]
 121    *
 122    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 123    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 124    * @f$\Delta=0,\xi=1@f$.
 125    *
 126    * @param x      Where to evaluate @f$ f'_{L}@f$
 127    * @param delta  Most probable value
 128    * @param xi     The 'width' of the distribution
 129    *
 130    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
 131    */
 132   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
 133
 134   //------------------------------------------------------------------
 135   /**
 136    * Calculate the value of a Landau convolved with a Gaussian
 137    *
 138    * @f[
 139    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 140    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 141    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 142    * @f]
 143    *
 144    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 145    * energy loss, @f$ \xi@f$ the width of the Landau, and
 146    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 147    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 148    * noise in the detector.
 149    *
 150    * Note that this function uses the constants fgConvolutionSteps and
 151    * fgConvolutionNSigma
 152    *
 153    * References:
 154    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 155    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 156    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 157    *
 158    * @param x         where to evaluate @f$ f@f$
 159    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 160    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 161    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 162    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 163    *
 164    * @return @f$ f@f$ evaluated at @f$ x@f$.
 165    */
 166   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi,
 167                              Double_t sigma, Double_t sigma_n);
 168
 169   //------------------------------------------------------------------
 170   /**
 171    * Evaluate
 172    * @f[
 173    *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 174    * @f]
 175    * corresponding to @f$ i@f$ particles i.e., with the substitutions
 176    * @f{eqnarray*}{
 177    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
 178    *    \xi       \rightarrow \xi_i       &=& i \xi\\
 179    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
 180    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 181    * @f}
 182    *
 183    * @param x        Where to evaluate
 184    * @param delta    @f$ \Delta@f$
 185    * @param xi       @f$ \xi@f$
 186    * @param sigma    @f$ \sigma@f$
 187    * @param sigma_n  @f$ \sigma_n@f$
 188    * @param i        @f$ i @f$
 189    *
 190    * @return @f$ f_i @f$ evaluated
 191    */
 192   static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 193                               Double_t sigma, Double_t sigma_n, Int_t i);
 194
 195   //------------------------------------------------------------------
 196   /**
 197    * Numerically evaluate
 198    * @f[
 199    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 200    * @f]
 201    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 202    * of the parameters is given by
 203    *
 204    * - 0: @f$\Delta@f$
 205    * - 1: @f$\xi@f$
 206    * - 2: @f$\sigma@f$
 207    * - 3: @f$\sigma_n@f$
 208    *
 209    * This is the partial derivative with respect to the parameter of
 210    * the response function corresponding to @f$ i@f$ particles i.e.,
 211    * with the substitutions
 212    * @f[
 213    *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
 214    *    \xi       \rightarrow \xi_i       = i \xi\\
 215    *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
 216    *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 217    * @f]
 218    *
 219    * @param x        Where to evaluate
 220    * @param ipar     Parameter number
 221    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 222    * @param delta    @f$ \Delta@f$
 223    * @param xi       @f$ \xi@f$
 224    * @param sigma    @f$ \sigma@f$
 225    * @param sigma_n  @f$ \sigma_n@f$
 226    * @param i        @f$ i@f$
 227    *
 228    * @return @f$ f_i@f$ evaluated
 229    */
 230   static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
 231                                    Double_t delta, Double_t xi,
 232                                    Double_t sigma, Double_t sigma_n, Int_t i);
 233
 234   //------------------------------------------------------------------
 235   /**
 236    * Evaluate
 237    * @f[
 238    *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 239    * @f]
 240    *
 241    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 242    * Landau with a Gaussian (see LandauGaus).  Note that
 243    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 244    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 245    *
 246    * References:
 247    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 248    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 249    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 250    *
 251    * @param x        Where to evaluate @f$ f_N@f$
 252    * @param delta    @f$ \Delta_1@f$
 253    * @param xi       @f$ \xi_1@f$
 254    * @param sigma    @f$ \sigma_1@f$
 255    * @param sigma_n  @f$ \sigma_n@f$
 256    * @param n        @f$ N@f$ in the sum above.
 257    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 258    *                 @f$ i > 1@f$
 259    *
 260    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
 261    */
 262   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 263                               Double_t sigma, Double_t sigma_n, Int_t n,
 264                               Double_t* a);
 265   /**
 266    * Generate a TF1 object of @f$ f_I@f$
 267    *
 268    * @param c        Constant
 269    * @param delta    @f$ \Delta@f$
 270    * @param xi       @f$ \xi_1@f$
 271    * @param sigma    @f$ \sigma_1@f$
 272    * @param sigma_n  @f$ \sigma_n@f$
 273    * @param i        @f$ i@f$ - the number of particles
 274    * @param xmin     Least value of range
 275    * @param xmax     Largest value of range
 276    *
 277    * @return Newly allocated TF1 object
 278    */
 279   static TF1* MakeILandauGaus(Double_t c,
 280                               Double_t delta, Double_t xi,
 281                               Double_t sigma, Double_t sigma_n,
 282                               Int_t    i,
 283                               Double_t xmin,  Double_t  xmax);
 284   /**
 285    * Generate a TF1 object of @f$ f_N@f$
 286    *
 287    * @param c         Constant
 288    * @param delta     @f$ \Delta@f$
 289    * @param xi        @f$ \xi_1@f$
 290    * @param sigma     @f$ \sigma_1@f$
 291    * @param sigma_n   @f$ \sigma_n@f$
 292    * @param n         @f$ N@f$ - how many particles to sum to
 293    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 294    *                  @f$ i > 1@f$
 295    * @param xmin      Least value of range
 296    * @param xmax      Largest value of range
 297    *
 298    * @return Newly allocated TF1 object
 299    */
 300   static TF1* MakeNLandauGaus(Double_t c,
 301                               Double_t delta, Double_t  xi,
 302                               Double_t sigma, Double_t  sigma_n,
 303                               Int_t    n,     Double_t* a,
 304                               Double_t xmin,  Double_t  xmax);
 305
 306   //__________________________________________________________________
 307   /**
 308    * Structure to do fits to the energy loss spectrum
 309    *
 310    * @ingroup pwg2_forward
 311    */
 312   struct ELossFitter
 313   {
 314     enum {
 315       kC     = 0,
 316       kDelta,
 317       kXi,
 318       kSigma,
 319       kSigmaN,
 320       kN,
 321       kA
 322     };
 323     /**
 324      * Constructor
 325      *
 326      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
 327      * @param maxRange   Maximum range to fit to
 328      * @param minusBins  The number of bins below maximum to use
 329      */
 330     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins);
 331     /**
 332      * Destructor
 333      *
 334      */
 335     virtual ~ELossFitter();
 336     /**
 337      * Clear internal arrays
 338      *
 339      */
 340     void Clear();
 341     /**
 342      * Fit a 1-particle signal to the passed energy loss distribution
 343      *
 344      * Note that this function clears the internal arrays first
 345      *
 346      * @param dist    Data to fit the function to
 347      * @param sigman If larger than zero, the initial guess of the
 348      *               detector induced noise. If zero or less, then this
 349      *               parameter is ignored in the fit (fixed at 0)
 350      *
 351      * @return The function fitted to the data
 352      */
 353     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
 354     /**
 355      * Fit a N-particle signal to the passed energy loss distribution
 356      *
 357      * If there's no 1-particle fit present, it does that first
 358      *
 359      * @param dist   Data to fit the function to
 360      * @param n      Number of particle signals to fit
 361      * @param sigman If larger than zero, the initial guess of the
 362      *               detector induced noise. If zero or less, then this
 363      *               parameter is ignored in the fit (fixed at 0)
 364      *
 365      * @return The function fitted to the data
 366      */
 367     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
 368
 369
 370     const Double_t fLowCut;     // Lower cut on data
 371     const Double_t fMaxRange;   // Maximum range to fit
 372     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
 373     TObjArray fFitResults;      // Array of fit results
 374     TObjArray fFunctions;       // Array of functions
 375   };
 376   /* @} */
 377
 378
 379   //==================================================================
 380   /**
 381    * @{
 382    * @name Convenience containers
 383    */
 384   /**
 385    * Structure to hold histograms
 386    *
 387    * @ingroup pwg2_forward
 388    */
 389   struct Histos : public TObject
 390   {
 391     /**
 392      * Constructor
 393      *
 394      *
 395      */
 396     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
 397     /**
 398      * Copy constructor
 399      *
 400      * @param o Object to copy from
 401      */
 402     Histos(const Histos& o)
 403       : TObject(o),
 404         fFMD1i(o.fFMD1i),
 405         fFMD2i(o.fFMD2i),
 406         fFMD2o(o.fFMD2o),
 407         fFMD3i(o.fFMD3i),
 408         fFMD3o(o.fFMD3o)
 409     {}
 410     /**
 411      * Assignement operator
 412      *
 413      * @return Reference to this
 414      */
 415     Histos& operator=(const Histos&) { return *this;}
 416     /**
 417      * Destructor
 418      */
 419     ~Histos();
 420     /**
 421      * Initialize the object
 422      *
 423      * @param etaAxis Eta axis to use
 424      */
 425     void Init(const TAxis& etaAxis);
 426     /**
 427      * Make a histogram
 428      *
 429      * @param d        Detector
 430      * @param r        Ring
 431      * @param etaAxis  Eta axis to use
 432      *
 433      * @return Newly allocated histogram
 434      */
 435     TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
 436     /**
 437      * Clear data
 438      *
 439      * @param option Not used
 440      */
 441     void  Clear(Option_t* option="");
 442     // const TH2D* Get(UShort_t d, Char_t r) const;
 443     /**
 444      * Get the histogram for a particular detector,ring
 445      *
 446      * @param d Detector
 447      * @param r Ring
 448      *
 449      * @return Histogram for detector,ring or nul
 450      */
 451     TH2D* Get(UShort_t d, Char_t r) const;
 452     TH2D* fFMD1i; // Histogram for FMD1i
 453     TH2D* fFMD2i; // Histogram for FMD2i
 454     TH2D* fFMD2o; // Histogram for FMD2o
 455     TH2D* fFMD3i; // Histogram for FMD3i
 456     TH2D* fFMD3o; // Histogram for FMD3o
 457
 458     ClassDef(Histos,1)
 459   };
 460
 461   //__________________________________________________________________
 462   /**
 463    * Base class for structure holding ring specific histograms
 464    *
 465    * @ingroup pwg2_forward
 466    */
 467   struct RingHistos : public TObject
 468   {
 469     /**
 470      * Constructor
 471      *
 472      */
 473     RingHistos() : fDet(0), fRing('\0'), fName("") {}
 474     /**
 475      *
 476      *
 477      * @param d Detector
 478      * @param r Ring
 479      */
 480     RingHistos(UShort_t d, Char_t r)
 481       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r))
 482     {}
 483     /**
 484      * Copy constructor
 485      *
 486      * @param o Object to copy from
 487      */
 488     RingHistos(const RingHistos& o)
 489       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
 490     {}
 491     /**
 492      *
 493      */
 494     virtual ~RingHistos() {}
 495     /**
 496      * Assignement operator
 497      *
 498      * @param o Object to assign from
 499      *
 500      * @return Reference to this
 501      */
 502     RingHistos& operator=(const RingHistos& o)
 503     {
 504       TObject::operator=(o);
 505       fDet  = o.fDet;
 506       fRing = o.fRing;
 507       fName = o.fName;
 508       return *this;
 509     }
 510     /**
 511      * Define the outout list in @a d
 512      *
 513      * @param d Where to put the output list
 514      *
 515      * @return Newly allocated TList object or null
 516      */
 517     TList* DefineOutputList(TList* d) const;
 518     /**
 519      * Get our output list from the container @a d
 520      *
 521      * @param d where to get the output list from
 522      *
 523      * @return The found TList or null
 524      */
 525     TList* GetOutputList(TList* d) const;
 526     /**
 527      * Find a specific histogram in the source list @a d
 528      *
 529      * @param d     (top)-container
 530      * @param name  Name of histogram
 531      *
 532      * @return Found histogram or null
 533      */
 534     TH1* GetOutputHist(TList* d, const char* name) const;
 535     /**
 536      *
 537      *
 538      *
 539      * @return
 540      */
 541     Color_t Color() const
 542     {
 543       return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
 544               + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
 545     }
 546     UShort_t fDet;   // Detector
 547     Char_t   fRing;  // Ring
 548     TString  fName;  // Name
 549
 550     ClassDef(RingHistos,1)
 551   };
 552   /* @} */
 553
 554 };
 555
 556 #endif
 557 // Local Variables:
 558 //  mode: C++
 559 // End:
 560