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