PWG2/FORWARD/analysis2/AliForwardUtil.h

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