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