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 TArrayD;
  26 class AliESDEvent;
  27 class AliAODEvent;
  28 class AliAnalysisTaskSE;
  29
  30 /**
  31  * Utilities used in the forward multiplcity analysis
  32  *
  33  * @ingroup pwglf_forward
  34  */
  35 class AliForwardUtil : public TObject
  36 {
  37 public:
  38   enum {
  39     kSkipRing = (1 << 19) // Bit for skipping a histogram
  40   };
  41   /**
  42    * Get the standard color for a ring
  43    *
  44    * @param d Detector
  45    * @param r Ring
  46    *
  47    * @return
  48    */
  49   static Color_t RingColor(UShort_t d, Char_t r)
  50   {
  51     return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
  52             + ((r == 'I' || r == 'i') ? 2 : -3));
  53   }
  54   //==================================================================
  55   /**
  56    * @{
  57    * @name AliROOT version
  58    */
  59   /**
  60    * Get the revision number of AliROOT
  61    *
  62    * @return Subversion revision number of AliROOT used
  63    */
  64   static ULong_t AliROOTRevision();
  65   /**
  66    * Get the branch identifier of AliROOT.  In case of trunk, return
  67    * 0xFFFFFFFF, while for @b vM-N-R{-S}, we get
  68    *
  69    * @code
  70    *    ((M & 0xFF) << 12 | (N & 0xFF) << 8 | (R & 0xFF) << 3 | (X))
  71    * @endcode
  72    * where @c X is 0xAA if @b S is specified (e.g., analysis tag).
  73    *
  74    * @return branch identifer encoded in bits
  75    */
  76   static ULong_t AliROOTBranch();
  77   //==================================================================
  78   /**
  79    * @{
  80    * @name Collision/run parameters
  81    */
  82   /**
  83    * Defined collision types
  84    */
  85   enum ECollisionSystem {
  86     kUnknown,
  87     kPP,
  88     kPbPb,
  89     kPPb
  90   };
  91   //__________________________________________________________________
  92   /**
  93    * Calculate the beam rapidity.
  94    *
  95    * @b Note: The beam energy is given in GeV/charge
  96    *
  97    * @param beam Beam energy in GeV/charge
  98    * @param z    Charge number of projectile
  99    * @param a    Mass number of projectile
 100    *
 101    * @return The rapidity of the beam
 102    */
 103   static Float_t BeamRapidity(Float_t beam, UShort_t z, UShort_t a);
 104   /**
 105    * Calculate the center of mass energy from the beam energy per
 106    * charge and the nucleus numbers.
 107    *
 108    * @param beam Beam energy in GeV/charge
 109    * @param z1   Charge number of projectile
 110    * @param a1   Mass number of projectile
 111    * @param z2   Charge number of projectile
 112    * @param a2   Mass number of projectile
 113    *
 114    * @return The center of mass energy in GeV/nucleon
 115    */
 116   static Float_t CenterOfMassEnergy(Float_t beam, UShort_t z1, UShort_t a1,
 117                                     Short_t z2=-1, Short_t a2=-1);
 118   /**
 119    * Calculate the center of mass rapidity (shift)
 120    *
 121    * @param z1 Charge number of projectile
 122    * @param a1 Mass number of projectile
 123    * @param z2 Charge number of projectile
 124    * @param a2 Mass number of projectile
 125    * x
 126    * @return Rapidity of the center of mass
 127    */
 128   static Float_t CenterOfMassRapidity(UShort_t z1, UShort_t a1,
 129                                       Short_t z2=-1, Short_t a2=-1);
 130   /**
 131    * Parse a collision system spec given in a string.   Known values are
 132    *
 133    *  - "pp", "p-p" which returns kPP
 134    *  - "PbPb", "Pb-Pb", "A-A", which returns kPbPb
 135    *  - "pPb", "p-Pb", "pA", p-A" which returns kPPb
 136    *  - Everything else gives kUnknown
 137    *
 138    * @param sys Collision system spec
 139    *
 140    * @return Collision system id
 141    */
 142   static UShort_t ParseCollisionSystem(const char* sys);
 143   /**
 144    * Get a string representation of the collision system
 145    *
 146    * @param sys  Collision system
 147    * - kPP -> "pp"
 148    * - kPbPb -> "PbPb"
 149    * - kPPb -> "pPb"
 150    * - anything else gives "unknown"
 151    *
 152    * @return String representation of the collision system
 153    */
 154   static const char* CollisionSystemString(UShort_t sys);
 155   //__________________________________________________________________
 156   /**
 157    * Parse the center of mass energy given as a float and return known
 158    * values as a unsigned integer
 159    *
 160    * @param sys  Collision system (needed for AA)
 161    * @param cms  Center of mass energy * total charge
 162    *
 163    * @return Center of mass energy per nucleon
 164    */
 165   static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
 166   /**
 167    * Get a string representation of the center of mass energy per nuclean
 168    *
 169    * @param cms  Center of mass energy per nucleon
 170    *
 171    * @return String representation of the center of mass energy per nuclean
 172    */
 173   static const char* CenterOfMassEnergyString(UShort_t cms);
 174   //__________________________________________________________________
 175   /**
 176    * Parse the magnetic field (in kG) as given by a floating point number
 177    *
 178    * @param field  Magnetic field in kG
 179    *
 180    * @return Short integer value of magnetic field in kG
 181    */
 182   static Short_t ParseMagneticField(Float_t field);
 183   /**
 184    * Get a string representation of the magnetic field
 185    *
 186    * @param field Magnetic field in kG
 187    *
 188    * @return String representation of the magnetic field
 189    */
 190   static const char* MagneticFieldString(Short_t field);
 191   /* @} */
 192
 193   //==================================================================
 194   /**
 195    * @{
 196    * @name Recalculate @f$\eta@f$, @f$\phi@f$, etc.
 197    */
 198   /**
 199    * Get the radius of a strip.
 200    *
 201    * @param ring  Ring identifier 'I' or 'O'
 202    * @param strip Strip number
 203    *
 204    * @return Radial distance from beam of the strip
 205    */
 206   static Double_t GetStripR(Char_t ring, UShort_t strip);
 207   /**
 208    * Get eta from strip
 209    *
 210    * @param det, ring, sec, strip, zvtx
 211    *
 212    * @return eta
 213    */
 214   static Double_t GetEtaFromStrip(UShort_t det, Char_t ring,
 215                                   UShort_t sec, UShort_t strip, Double_t zvtx);
 216   /**
 217    * Get the azimuthal angle of a strip
 218    *
 219    * @param ring  Ring identifier 'I' or 'O'
 220    * @param strip Strip number
 221    * @param phi   Straight forward strip phi
 222    * @param xvtx  Ip X coordinate
 223    * @param yvtx  Ip Y coordinate
 224    *
 225    * @return The phi angle correctef for (X,Y) off set.
 226    */
 227   static Double_t GetPhiFromStrip(Char_t ring, UShort_t strip,
 228                                   Double_t phi, Double_t xvtx, Double_t yvtx);
 229   /* @} */
 230
 231   //==================================================================
 232   /**
 233    * @{
 234    * @name Manager related tasks
 235    */
 236   /**
 237    * Get the AOD event - either from the input (AOD analysis) or the
 238    * output (ESD analysis)
 239    *
 240    * @param task Task to do the investigation for
 241    *
 242    * @return Found AOD event or null
 243    */
 244   static AliAODEvent* GetAODEvent(AliAnalysisTaskSE* task);
 245   /**
 246    * Check if we have something that will provide and AOD event
 247    *
 248    * @return 0 if there's nothing that provide an AOD event, 1 if it
 249    * is provided on the input (AOD analysis) or 2 if it is provided on
 250    * the output (ESD analysis)
 251    */
 252   static UShort_t CheckForAOD();
 253   /**
 254    * Check if we have a particular (kind) of task in our train
 255    *
 256    * @param clsOrName  Class name or name of task
 257    * @param cls If true, look for a task of a particular class -
 258    * otherwise search for a speficially name task
 259    *
 260    * @return true if the needed task was found
 261    */
 262   static Bool_t CheckForTask(const char* clsOrName, Bool_t cls=true);
 263   /* @} */
 264
 265   //==================================================================
 266   /**
 267    * @{
 268    * @name Member functions to store and retrieve analysis parameters
 269    */
 270   static TObject* MakeParameter(const char* name, UShort_t value);
 271   static TObject* MakeParameter(const char* name, Int_t value);
 272   static TObject* MakeParameter(const char* name, Double_t value);
 273   static TObject* MakeParameter(const char* name, Bool_t value);
 274   static TObject* MakeParameter(const char* name, ULong_t value);
 275   static void GetParameter(TObject* o, UShort_t& value);
 276   static void GetParameter(TObject* o, Int_t& value);
 277   static void GetParameter(TObject* o, Double_t& value);
 278   static void GetParameter(TObject* o, Bool_t& value);
 279   static void GetParameter(TObject* o, ULong_t& value);
 280   /* @} */
 281
 282   //==================================================================
 283   /**
 284    * @{
 285    * @name Axis functions
 286    */
 287   static TAxis* MakeFullIpZAxis(Int_t nCenter=20);
 288   static void MakeFullIpZAxis(Int_t nCenter, TArrayD& bins);
 289   static void MakeLogScale(Int_t nBins, Int_t minOrder, Int_t maxOrder, TArrayD& bins);
 290   static void PrintTask(const TObject& o);
 291   static void PrintName(const char* name);
 292   static void PrintField(const char* name, const char* value, ...);
 293   /* @} */
 294
 295   //==================================================================
 296 #if 0 // Moved to separate classes
 297   /**
 298    * @{
 299    * @name Energy stragling functions
 300    */
 301   //__________________________________________________________________
 302   /**
 303    * Number of steps to do in the Landau, Gaussiam convolution
 304    */
 305   static Int_t fgConvolutionSteps; // Number of convolution steps
 306   //------------------------------------------------------------------
 307   /**
 308    * How many sigma's of the Gaussian in the Landau, Gaussian
 309    * convolution to integrate over
 310    */
 311   static Double_t fgConvolutionNSigma; // Number of convolution sigmas
 312   //------------------------------------------------------------------
 313   /**
 314    * Calculate the shifted Landau
 315    * @f[
 316    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 317    * @f]
 318    *
 319    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 320    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 321    * @f$\Delta=0,\xi=1@f$.
 322    *
 323    * @param x      Where to evaluate @f$ f'_{L}@f$
 324    * @param delta  Most probable value
 325    * @param xi     The 'width' of the distribution
 326    *
 327    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
 328    */
 329   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
 330
 331   //------------------------------------------------------------------
 332   /**
 333    * Calculate the value of a Landau convolved with a Gaussian
 334    *
 335    * @f[
 336    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 337    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 338    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 339    * @f]
 340    *
 341    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 342    * energy loss, @f$ \xi@f$ the width of the Landau, and
 343    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 344    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 345    * noise in the detector.
 346    *
 347    * Note that this function uses the constants fgConvolutionSteps and
 348    * fgConvolutionNSigma
 349    *
 350    * References:
 351    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 352    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 353    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 354    *
 355    * @param x         where to evaluate @f$ f@f$
 356    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 357    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 358    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 359    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 360    *
 361    * @return @f$ f@f$ evaluated at @f$ x@f$.
 362    */
 363   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi,
 364                              Double_t sigma, Double_t sigma_n);
 365
 366   //------------------------------------------------------------------
 367   /**
 368    * Evaluate
 369    * @f[
 370    *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 371    * @f]
 372    * corresponding to @f$ i@f$ particles i.e., with the substitutions
 373    * @f{eqnarray*}{
 374    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
 375    *    \xi       \rightarrow \xi_i       &=& i \xi\\
 376    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
 377    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 378    * @f}
 379    *
 380    * @param x        Where to evaluate
 381    * @param delta    @f$ \Delta@f$
 382    * @param xi       @f$ \xi@f$
 383    * @param sigma    @f$ \sigma@f$
 384    * @param sigma_n  @f$ \sigma_n@f$
 385    * @param i        @f$ i @f$
 386    *
 387    * @return @f$ f_i @f$ evaluated
 388    */
 389   static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 390                               Double_t sigma, Double_t sigma_n, Int_t i);
 391
 392   //------------------------------------------------------------------
 393   /**
 394    * Numerically evaluate
 395    * @f[
 396    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 397    * @f]
 398    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 399    * of the parameters is given by
 400    *
 401    * - 0: @f$\Delta@f$
 402    * - 1: @f$\xi@f$
 403    * - 2: @f$\sigma@f$
 404    * - 3: @f$\sigma_n@f$
 405    *
 406    * This is the partial derivative with respect to the parameter of
 407    * the response function corresponding to @f$ i@f$ particles i.e.,
 408    * with the substitutions
 409    * @f[
 410    *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
 411    *    \xi       \rightarrow \xi_i       = i \xi\\
 412    *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
 413    *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 414    * @f]
 415    *
 416    * @param x        Where to evaluate
 417    * @param ipar     Parameter number
 418    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 419    * @param delta    @f$ \Delta@f$
 420    * @param xi       @f$ \xi@f$
 421    * @param sigma    @f$ \sigma@f$
 422    * @param sigma_n  @f$ \sigma_n@f$
 423    * @param i        @f$ i@f$
 424    *
 425    * @return @f$ f_i@f$ evaluated
 426    */
 427   static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
 428                                    Double_t delta, Double_t xi,
 429                                    Double_t sigma, Double_t sigma_n, Int_t i);
 430
 431   //------------------------------------------------------------------
 432   /**
 433    * Evaluate
 434    * @f[
 435    *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 436    * @f]
 437    *
 438    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 439    * Landau with a Gaussian (see LandauGaus).  Note that
 440    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 441    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 442    *
 443    * References:
 444    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 445    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 446    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 447    *
 448    * @param x        Where to evaluate @f$ f_N@f$
 449    * @param delta    @f$ \Delta_1@f$
 450    * @param xi       @f$ \xi_1@f$
 451    * @param sigma    @f$ \sigma_1@f$
 452    * @param sigma_n  @f$ \sigma_n@f$
 453    * @param n        @f$ N@f$ in the sum above.
 454    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 455    *                 @f$ i > 1@f$
 456    *
 457    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
 458    */
 459   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 460                               Double_t sigma, Double_t sigma_n, Int_t n,
 461                               const Double_t* a);
 462   /**
 463    * Generate a TF1 object of @f$ f_I@f$
 464    *
 465    * @param c        Constant
 466    * @param delta    @f$ \Delta@f$
 467    * @param xi       @f$ \xi_1@f$
 468    * @param sigma    @f$ \sigma_1@f$
 469    * @param sigma_n  @f$ \sigma_n@f$
 470    * @param i        @f$ i@f$ - the number of particles
 471    * @param xmin     Least value of range
 472    * @param xmax     Largest value of range
 473    *
 474    * @return Newly allocated TF1 object
 475    */
 476   static TF1* MakeILandauGaus(Double_t c,
 477                               Double_t delta, Double_t xi,
 478                               Double_t sigma, Double_t sigma_n,
 479                               Int_t    i,
 480                               Double_t xmin,  Double_t  xmax);
 481   /**
 482    * Generate a TF1 object of @f$ f_N@f$
 483    *
 484    * @param c         Constant
 485    * @param delta     @f$ \Delta@f$
 486    * @param xi        @f$ \xi_1@f$
 487    * @param sigma     @f$ \sigma_1@f$
 488    * @param sigma_n   @f$ \sigma_n@f$
 489    * @param n         @f$ N@f$ - how many particles to sum to
 490    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 491    *                  @f$ i > 1@f$
 492    * @param xmin      Least value of range
 493    * @param xmax      Largest value of range
 494    *
 495    * @return Newly allocated TF1 object
 496    */
 497   static TF1* MakeNLandauGaus(Double_t c,
 498                               Double_t delta, Double_t  xi,
 499                               Double_t sigma, Double_t  sigma_n,
 500                               Int_t    n,     const Double_t* a,
 501                               Double_t xmin,  Double_t  xmax);
 502
 503   //__________________________________________________________________
 504   /**
 505    * Structure to do fits to the energy loss spectrum
 506    *
 507    * @ingroup pwglf_forward
 508    */
 509   struct ELossFitter
 510   {
 511     enum {
 512       kC     = 0,
 513       kDelta,
 514       kXi,
 515       kSigma,
 516       kSigmaN,
 517       kN,
 518       kA
 519     };
 520     /**
 521      * Constructor
 522      *
 523      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
 524      * @param maxRange   Maximum range to fit to
 525      * @param minusBins  The number of bins below maximum to use
 526      */
 527     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins);
 528     /**
 529      * Destructor
 530      *
 531      */
 532     virtual ~ELossFitter();
 533     void SetDebug(Bool_t debug=true) { fDebug = debug; }
 534     /**
 535      * Clear internal arrays
 536      *
 537      */
 538     void Clear();
 539     /**
 540      * Fit a 1-particle signal to the passed energy loss distribution
 541      *
 542      * Note that this function clears the internal arrays first
 543      *
 544      * @param dist    Data to fit the function to
 545      * @param sigman If larger than zero, the initial guess of the
 546      *               detector induced noise. If zero or less, then this
 547      *               parameter is ignored in the fit (fixed at 0)
 548      *
 549      * @return The function fitted to the data
 550      */
 551     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
 552     /**
 553      * Fit a N-particle signal to the passed energy loss distribution
 554      *
 555      * If there's no 1-particle fit present, it does that first
 556      *
 557      * @param dist   Data to fit the function to
 558      * @param n      Number of particle signals to fit
 559      * @param sigman If larger than zero, the initial guess of the
 560      *               detector induced noise. If zero or less, then this
 561      *               parameter is ignored in the fit (fixed at 0)
 562      *
 563      * @return The function fitted to the data
 564      */
 565     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
 566     /**
 567      * Fit a composite distribution of energy loss from both primaries
 568      * and secondaries
 569      *
 570      * @param dist   Distribution
 571      * @param sigman If larger than zero, the initial guess of the
 572      *                detector included noise.  If zero or less this
 573      *                parameter is fixed to 0.
 574      *
 575      * @return Function fitted to the data
 576      */
 577     TF1* FitComposite(TH1* dist, Double_t sigman);
 578     /**
 579      * Get Lower cut on data
 580      *
 581      * @return Lower cut on data
 582      */
 583     Double_t GetLowCut() const { return fLowCut; }
 584     /**
 585      * Get Maximum range to fit
 586      *
 587      * @return Maximum range to fit
 588      */
 589     Double_t GetMaxRange() const { return fMaxRange; }
 590     /**
 591      * Get Number of bins from maximum to fit 1st peak
 592      *
 593      * @return Number of bins from maximum to fit 1st peak
 594      */
 595     UShort_t GetMinusBins() const { return fMinusBins; }
 596     /**
 597      * Get Array of fit results
 598      *
 599      * @return Array of fit results
 600      */
 601     const TObjArray& GetFitResults() const { return fFitResults; }
 602     /**
 603      * Get Array of fit results
 604      *
 605      * @return Array of fit results
 606      */
 607     TObjArray& GetFitResults() { return fFitResults; }
 608     /**
 609      * Get Array of functions
 610      *
 611      * @return Array of functions
 612      */
 613     const TObjArray& GetFunctions() const { return fFunctions; }
 614     /**
 615      * Get Array of functions
 616      *
 617      * @return Array of functions
 618      */
 619     TObjArray& GetFunctions() { return fFunctions; }
 620   private:
 621     const Double_t fLowCut;     // Lower cut on data
 622     const Double_t fMaxRange;   // Maximum range to fit
 623     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
 624     TObjArray fFitResults;      // Array of fit results
 625     TObjArray fFunctions;       // Array of functions
 626     Bool_t    fDebug;
 627   };
 628   /* @} */
 629 #endif
 630
 631   //==================================================================
 632   /**
 633    * @{
 634    * @name Convenience containers
 635    */
 636   /**
 637    * Structure to hold histograms
 638    *
 639    * @ingroup pwglf_forward
 640    */
 641   struct Histos : public TObject
 642   {
 643     /**
 644      * Constructor
 645      *
 646      *
 647      */
 648     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
 649     /**
 650      * Copy constructor
 651      *
 652      * @param o Object to copy from
 653      */
 654     Histos(const Histos& o)
 655       : TObject(o),
 656         fFMD1i(o.fFMD1i),
 657         fFMD2i(o.fFMD2i),
 658         fFMD2o(o.fFMD2o),
 659         fFMD3i(o.fFMD3i),
 660         fFMD3o(o.fFMD3o)
 661     {}
 662     /**
 663      * Assignement operator
 664      *
 665      * @return Reference to this
 666      */
 667     Histos& operator=(const Histos&) { return *this;}
 668     /**
 669      * Destructor.  This does not delete the interally allocated
 670      * memory.  Use the member function Delete for that.
 671      */
 672     ~Histos();
 673     /**
 674      * Clear internal memory.  Note, if the internal histograms are
 675      * added to an output container, then we must not free this
 676      * memory.
 677      */
 678     void Delete(Option_t* opt="");
 679     /**
 680      * Initialize the object
 681      *
 682      * @param etaAxis Eta axis to use
 683      */
 684     void Init(const TAxis& etaAxis);
 685     /**
 686      * Re-initialize the object with new @f$\eta@f$ axis
 687      *
 688      * @param etaAxis Eta axis to use
 689      */
 690     void ReInit(const TAxis& etaAxis);
 691     /**
 692      * Make a histogram
 693      *
 694      * @param d        Detector
 695      * @param r        Ring
 696      * @param etaAxis  Eta axis to use
 697      *
 698      * @return Newly allocated histogram
 699      */
 700     static TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis);
 701     /**
 702      * Set the @f$\eta@f$ axis
 703      *
 704      * @param hist    Histogram
 705      * @param etaAxis @f$\eta@f$ axis to use
 706      */
 707     static void RebinEta(TH2D* hist, const TAxis& etaAxis);
 708     /**
 709      * Clear data
 710      *
 711      * @param option Not used
 712      */
 713     void  Clear(Option_t* option="");
 714     // const TH2D* Get(UShort_t d, Char_t r) const;
 715     /**
 716      * Get the histogram for a particular detector,ring
 717      *
 718      * @param d Detector
 719      * @param r Ring
 720      *
 721      * @return Histogram for detector,ring or nul
 722      */
 723     TH2D* Get(UShort_t d, Char_t r) const;
 724     TH2D* fFMD1i; // Histogram for FMD1i
 725     TH2D* fFMD2i; // Histogram for FMD2i
 726     TH2D* fFMD2o; // Histogram for FMD2o
 727     TH2D* fFMD3i; // Histogram for FMD3i
 728     TH2D* fFMD3o; // Histogram for FMD3o
 729
 730     ClassDef(Histos,2)
 731   };
 732
 733   //__________________________________________________________________
 734   /**
 735    * Base class for structure holding ring specific histograms
 736    *
 737    * @ingroup pwglf_forward
 738    */
 739   struct RingHistos : public TObject
 740   {
 741     /**
 742      * Constructor
 743      *
 744      */
 745     RingHistos() : fDet(0), fRing('\0'), fName(""), fkNSector(0), fkNStrip(0) {}
 746     /**
 747      *
 748      *
 749      * @param d Detector
 750      * @param r Ring
 751      */
 752     RingHistos(UShort_t d, Char_t r)
 753       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)),
 754         fkNSector(r == 'i' || r == 'I' ? 20 : 40),
 755         fkNStrip(r == 'i' || r == 'I' ? 512 : 256)
 756     {}
 757     /**
 758      * Copy constructor
 759      *
 760      * @param o Object to copy from
 761      */
 762     RingHistos(const RingHistos& o)
 763       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName),
 764         fkNSector(o.fkNSector), fkNStrip(o.fkNStrip)
 765     {}
 766     /**
 767      *
 768      */
 769     virtual ~RingHistos() {}
 770     /**
 771      * Assignement operator
 772      *
 773      * @param o Object to assign from
 774      *
 775      * @return Reference to this
 776      */
 777     RingHistos& operator=(const RingHistos& o)
 778     {
 779       if (&o == this) return *this;
 780       TObject::operator=(o);
 781       fDet  = o.fDet;
 782       fRing = o.fRing;
 783       fName = o.fName;
 784       fkNSector = o.fkNSector;
 785       fkNStrip  = o.fkNStrip;
 786       return *this;
 787     }
 788     /**
 789      * Define the outout list in @a d
 790      *
 791      * @param d Where to put the output list
 792      *
 793      * @return Newly allocated TList object or null
 794      */
 795     TList* DefineOutputList(TList* d) const;
 796     /**
 797      * Get our output list from the container @a d
 798      *
 799      * @param d where to get the output list from
 800      *
 801      * @return The found TList or null
 802      */
 803     TList* GetOutputList(const TList* d) const;
 804     /**
 805      * Find a specific histogram in the source list @a d
 806      *
 807      * @param d     (top)-container
 808      * @param name  Name of histogram
 809      *
 810      * @return Found histogram or null
 811      */
 812     TH1* GetOutputHist(const TList* d, const char* name) const;
 813     /**
 814      * Get the colour of this ring
 815      *
 816      *
 817      * @return
 818      */
 819     Color_t Color() const
 820     {
 821       return AliForwardUtil::RingColor(fDet, fRing);
 822     }
 823     /**
 824      * The name of this ring
 825      *
 826      * @return Name of this ring
 827      */
 828     const char* GetName() const { return fName.Data(); }
 829     /**
 830      * Get number of sectors
 831      */
 832     const UShort_t& NSector() const { return fkNSector; }
 833     /**
 834      * Get number of strips
 835      */
 836     const UShort_t& NStrip() const { return fkNStrip; }
 837     UShort_t fDet;   // Detector
 838     Char_t   fRing;  // Ring
 839     TString  fName;  // Name
 840     UShort_t fkNSector; // Number of sectors
 841     UShort_t fkNStrip;  // Number of strips
 842
 843     ClassDef(RingHistos,1)
 844   };
 845   /* @} */
 846
 847   //__________________________________________________________________
 848   /**
 849    * A guard idom for producing debug output
 850    *
 851    */
 852   struct DebugGuard
 853   {
 854     /**
 855      * Constructor
 856      *
 857      * @param lvl       Current level
 858      * @param msgLvl    Target level
 859      * @param format    @c printf -like format
 860      *
 861      * @return
 862      */
 863     DebugGuard(Int_t lvl, Int_t msgLvl, const char* format, ...);
 864     /**
 865      * Destructor
 866      */
 867     ~DebugGuard();
 868     /**
 869      * Make a message
 870      *
 871      * @param lvl    Current level
 872      * @param msgLvl Target level
 873      * @param format @c printf -like format
 874      */
 875     static void Message(Int_t lvl, Int_t msgLvl, const char* format, ...);
 876   private:
 877     /**
 878      * Output the message
 879      *
 880      * @param in    Direction
 881      * @param msg   Message
 882      */
 883     static void Output(int in, TString& msg);
 884     /**
 885      * Format a message
 886      *
 887      * @param out     Output is stored here
 888      * @param format  @c printf -like format
 889      * @param ap      List of arguments
 890      */
 891     static void Format(TString& out, const char* format, va_list ap);
 892     TString fMsg;
 893   };
 894 private:
 895   /**
 896    * Constructor
 897    */
 898   AliForwardUtil() {}
 899   /**
 900    * Copy constructor
 901    *
 902    * @param o Object to copy from
 903    */
 904   AliForwardUtil(const AliForwardUtil& o) : TObject(o) {}
 905   /**
 906    * Assingment operator
 907    *
 908    *
 909    * @return Reference to this object
 910    */
 911   AliForwardUtil& operator=(const AliForwardUtil&) { return *this; }
 912   /**
 913    * Destructor
 914    */
 915   ~AliForwardUtil() {}
 916
 917
 918   ClassDef(AliForwardUtil,1) // Utilities - do not make object
 919 };
 920
 921 // #ifdef LOG_NO_DEBUG
 922 // # define DGUARD(L,N,F,...) do {} while(false)
 923 // #else
 924 /**
 925  * Macro to declare a DebugGuard
 926  *
 927  * @param L Current debug level
 928  * @param N Target debug level
 929  * @param F @c printf -like Format
 930  */
 931 # define DGUARD(L,N,F,...)                                      \
 932   AliForwardUtil::DebugGuard _GUARD(L,N,F, ## __VA_ARGS__)
 933 /**
 934  * Macro to make a debug message, using DebugGuard::Message
 935  *
 936  * @param L Current debug level
 937  * @param N Target debug level
 938  * @param F @c printf -like Format
 939  */
 940 # define DMSG(L,N,F,...)                                        \
 941   AliForwardUtil::DebugGuard::Message(L,N,F, ## __VA_ARGS__)
 942 // #endif
 943 #endif
 944 // Local Variables:
 945 //  mode: C++
 946 // End:
 947