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