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   /**
 297    * @{
 298    * @name Energy stragling functions
 299    */
 300   //__________________________________________________________________
 301   /**
 302    * Number of steps to do in the Landau, Gaussiam convolution
 303    */
 304   static Int_t fgConvolutionSteps; // Number of convolution steps
 305   //------------------------------------------------------------------
 306   /**
 307    * How many sigma's of the Gaussian in the Landau, Gaussian
 308    * convolution to integrate over
 309    */
 310   static Double_t fgConvolutionNSigma; // Number of convolution sigmas
 311   //------------------------------------------------------------------
 312   /**
 313    * Calculate the shifted Landau
 314    * @f[
 315    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 316    * @f]
 317    *
 318    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 319    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 320    * @f$\Delta=0,\xi=1@f$.
 321    *
 322    * @param x      Where to evaluate @f$ f'_{L}@f$
 323    * @param delta  Most probable value
 324    * @param xi     The 'width' of the distribution
 325    *
 326    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
 327    */
 328   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
 329
 330   //------------------------------------------------------------------
 331   /**
 332    * Calculate the value of a Landau convolved with a Gaussian
 333    *
 334    * @f[
 335    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 336    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 337    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 338    * @f]
 339    *
 340    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 341    * energy loss, @f$ \xi@f$ the width of the Landau, and
 342    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 343    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 344    * noise in the detector.
 345    *
 346    * Note that this function uses the constants fgConvolutionSteps and
 347    * fgConvolutionNSigma
 348    *
 349    * References:
 350    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 351    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 352    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 353    *
 354    * @param x         where to evaluate @f$ f@f$
 355    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 356    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 357    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 358    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 359    *
 360    * @return @f$ f@f$ evaluated at @f$ x@f$.
 361    */
 362   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi,
 363                              Double_t sigma, Double_t sigma_n);
 364
 365   //------------------------------------------------------------------
 366   /**
 367    * Evaluate
 368    * @f[
 369    *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 370    * @f]
 371    * corresponding to @f$ i@f$ particles i.e., with the substitutions
 372    * @f{eqnarray*}{
 373    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
 374    *    \xi       \rightarrow \xi_i       &=& i \xi\\
 375    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
 376    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 377    * @f}
 378    *
 379    * @param x        Where to evaluate
 380    * @param delta    @f$ \Delta@f$
 381    * @param xi       @f$ \xi@f$
 382    * @param sigma    @f$ \sigma@f$
 383    * @param sigma_n  @f$ \sigma_n@f$
 384    * @param i        @f$ i @f$
 385    *
 386    * @return @f$ f_i @f$ evaluated
 387    */
 388   static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 389                               Double_t sigma, Double_t sigma_n, Int_t i);
 390
 391   //------------------------------------------------------------------
 392   /**
 393    * Numerically evaluate
 394    * @f[
 395    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 396    * @f]
 397    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 398    * of the parameters is given by
 399    *
 400    * - 0: @f$\Delta@f$
 401    * - 1: @f$\xi@f$
 402    * - 2: @f$\sigma@f$
 403    * - 3: @f$\sigma_n@f$
 404    *
 405    * This is the partial derivative with respect to the parameter of
 406    * the response function corresponding to @f$ i@f$ particles i.e.,
 407    * with the substitutions
 408    * @f[
 409    *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
 410    *    \xi       \rightarrow \xi_i       = i \xi\\
 411    *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
 412    *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 413    * @f]
 414    *
 415    * @param x        Where to evaluate
 416    * @param ipar     Parameter number
 417    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 418    * @param delta    @f$ \Delta@f$
 419    * @param xi       @f$ \xi@f$
 420    * @param sigma    @f$ \sigma@f$
 421    * @param sigma_n  @f$ \sigma_n@f$
 422    * @param i        @f$ i@f$
 423    *
 424    * @return @f$ f_i@f$ evaluated
 425    */
 426   static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
 427                                    Double_t delta, Double_t xi,
 428                                    Double_t sigma, Double_t sigma_n, Int_t i);
 429
 430   //------------------------------------------------------------------
 431   /**
 432    * Evaluate
 433    * @f[
 434    *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 435    * @f]
 436    *
 437    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 438    * Landau with a Gaussian (see LandauGaus).  Note that
 439    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 440    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 441    *
 442    * References:
 443    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 444    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 445    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 446    *
 447    * @param x        Where to evaluate @f$ f_N@f$
 448    * @param delta    @f$ \Delta_1@f$
 449    * @param xi       @f$ \xi_1@f$
 450    * @param sigma    @f$ \sigma_1@f$
 451    * @param sigma_n  @f$ \sigma_n@f$
 452    * @param n        @f$ N@f$ in the sum above.
 453    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 454    *                 @f$ i > 1@f$
 455    *
 456    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$
 457    */
 458   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 459                               Double_t sigma, Double_t sigma_n, Int_t n,
 460                               const Double_t* a);
 461   /**
 462    * Generate a TF1 object of @f$ f_I@f$
 463    *
 464    * @param c        Constant
 465    * @param delta    @f$ \Delta@f$
 466    * @param xi       @f$ \xi_1@f$
 467    * @param sigma    @f$ \sigma_1@f$
 468    * @param sigma_n  @f$ \sigma_n@f$
 469    * @param i        @f$ i@f$ - the number of particles
 470    * @param xmin     Least value of range
 471    * @param xmax     Largest value of range
 472    *
 473    * @return Newly allocated TF1 object
 474    */
 475   static TF1* MakeILandauGaus(Double_t c,
 476                               Double_t delta, Double_t xi,
 477                               Double_t sigma, Double_t sigma_n,
 478                               Int_t    i,
 479                               Double_t xmin,  Double_t  xmax);
 480   /**
 481    * Generate a TF1 object of @f$ f_N@f$
 482    *
 483    * @param c         Constant
 484    * @param delta     @f$ \Delta@f$
 485    * @param xi        @f$ \xi_1@f$
 486    * @param sigma     @f$ \sigma_1@f$
 487    * @param sigma_n   @f$ \sigma_n@f$
 488    * @param n         @f$ N@f$ - how many particles to sum to
 489    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 490    *                  @f$ i > 1@f$
 491    * @param xmin      Least value of range
 492    * @param xmax      Largest value of range
 493    *
 494    * @return Newly allocated TF1 object
 495    */
 496   static TF1* MakeNLandauGaus(Double_t c,
 497                               Double_t delta, Double_t  xi,
 498                               Double_t sigma, Double_t  sigma_n,
 499                               Int_t    n,     const Double_t* a,
 500                               Double_t xmin,  Double_t  xmax);
 501
 502   //__________________________________________________________________
 503   /**
 504    * Structure to do fits to the energy loss spectrum
 505    *
 506    * @ingroup pwglf_forward
 507    */
 508   struct ELossFitter
 509   {
 510     enum {
 511       kC     = 0,
 512       kDelta,
 513       kXi,
 514       kSigma,
 515       kSigmaN,
 516       kN,
 517       kA
 518     };
 519     /**
 520      * Constructor
 521      *
 522      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
 523      * @param maxRange   Maximum range to fit to
 524      * @param minusBins  The number of bins below maximum to use
 525      */
 526     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins);
 527     /**
 528      * Destructor
 529      *
 530      */
 531     virtual ~ELossFitter();
 532     void SetDebug(Bool_t debug=true) { fDebug = debug; }
 533     /**
 534      * Clear internal arrays
 535      *
 536      */
 537     void Clear();
 538     /**
 539      * Fit a 1-particle signal to the passed energy loss distribution
 540      *
 541      * Note that this function clears the internal arrays first
 542      *
 543      * @param dist    Data to fit the function to
 544      * @param sigman If larger than zero, the initial guess of the
 545      *               detector induced noise. If zero or less, then this
 546      *               parameter is ignored in the fit (fixed at 0)
 547      *
 548      * @return The function fitted to the data
 549      */
 550     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
 551     /**
 552      * Fit a N-particle signal to the passed energy loss distribution
 553      *
 554      * If there's no 1-particle fit present, it does that first
 555      *
 556      * @param dist   Data to fit the function to
 557      * @param n      Number of particle signals to fit
 558      * @param sigman If larger than zero, the initial guess of the
 559      *               detector induced noise. If zero or less, then this
 560      *               parameter is ignored in the fit (fixed at 0)
 561      *
 562      * @return The function fitted to the data
 563      */
 564     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
 565     /**
 566      * Fit a composite distribution of energy loss from both primaries
 567      * and secondaries
 568      *
 569      * @param dist   Distribution
 570      * @param sigman If larger than zero, the initial guess of the
 571      *                detector included noise.  If zero or less this
 572      *                parameter is fixed to 0.
 573      *
 574      * @return Function fitted to the data
 575      */
 576     TF1* FitComposite(TH1* dist, Double_t sigman);
 577     /**
 578      * Get Lower cut on data
 579      *
 580      * @return Lower cut on data
 581      */
 582     Double_t GetLowCut() const { return fLowCut; }
 583     /**
 584      * Get Maximum range to fit
 585      *
 586      * @return Maximum range to fit
 587      */
 588     Double_t GetMaxRange() const { return fMaxRange; }
 589     /**
 590      * Get Number of bins from maximum to fit 1st peak
 591      *
 592      * @return Number of bins from maximum to fit 1st peak
 593      */
 594     UShort_t GetMinusBins() const { return fMinusBins; }
 595     /**
 596      * Get Array of fit results
 597      *
 598      * @return Array of fit results
 599      */
 600     const TObjArray& GetFitResults() const { return fFitResults; }
 601     /**
 602      * Get Array of fit results
 603      *
 604      * @return Array of fit results
 605      */
 606     TObjArray& GetFitResults() { return fFitResults; }
 607     /**
 608      * Get Array of functions
 609      *
 610      * @return Array of functions
 611      */
 612     const TObjArray& GetFunctions() const { return fFunctions; }
 613     /**
 614      * Get Array of functions
 615      *
 616      * @return Array of functions
 617      */
 618     TObjArray& GetFunctions() { return fFunctions; }
 619   private:
 620     const Double_t fLowCut;     // Lower cut on data
 621     const Double_t fMaxRange;   // Maximum range to fit
 622     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
 623     TObjArray fFitResults;      // Array of fit results
 624     TObjArray fFunctions;       // Array of functions
 625     Bool_t    fDebug;
 626   };
 627   /* @} */
 628
 629
 630   //==================================================================
 631   /**
 632    * @{
 633    * @name Convenience containers
 634    */
 635   /**
 636    * Structure to hold histograms
 637    *
 638    * @ingroup pwglf_forward
 639    */
 640   struct Histos : public TObject
 641   {
 642     /**
 643      * Constructor
 644      *
 645      *
 646      */
 647     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
 648     /**
 649      * Copy constructor
 650      *
 651      * @param o Object to copy from
 652      */
 653     Histos(const Histos& o)
 654       : TObject(o),
 655         fFMD1i(o.fFMD1i),
 656         fFMD2i(o.fFMD2i),
 657         fFMD2o(o.fFMD2o),
 658         fFMD3i(o.fFMD3i),
 659         fFMD3o(o.fFMD3o)
 660     {}
 661     /**
 662      * Assignement operator
 663      *
 664      * @return Reference to this
 665      */
 666     Histos& operator=(const Histos&) { return *this;}
 667     /**
 668      * Destructor.  This does not delete the interally allocated
 669      * memory.  Use the member function Delete for that.
 670      */
 671     ~Histos();
 672     /**
 673      * Clear internal memory.  Note, if the internal histograms are
 674      * added to an output container, then we must not free this
 675      * memory.
 676      */
 677     void Delete(Option_t* opt="");
 678     /**
 679      * Initialize the object
 680      *
 681      * @param etaAxis Eta axis to use
 682      */
 683     void Init(const TAxis& etaAxis);
 684     /**
 685      * Re-initialize the object with new @f$\eta@f$ axis
 686      *
 687      * @param etaAxis Eta axis to use
 688      */
 689     void ReInit(const TAxis& etaAxis);
 690     /**
 691      * Make a histogram
 692      *
 693      * @param d        Detector
 694      * @param r        Ring
 695      * @param etaAxis  Eta axis to use
 696      *
 697      * @return Newly allocated histogram
 698      */
 699     static TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis);
 700     /**
 701      * Set the @f$\eta@f$ axis
 702      *
 703      * @param hist    Histogram
 704      * @param etaAxis @f$\eta@f$ axis to use
 705      */
 706     static void RebinEta(TH2D* hist, const TAxis& etaAxis);
 707     /**
 708      * Clear data
 709      *
 710      * @param option Not used
 711      */
 712     void  Clear(Option_t* option="");
 713     // const TH2D* Get(UShort_t d, Char_t r) const;
 714     /**
 715      * Get the histogram for a particular detector,ring
 716      *
 717      * @param d Detector
 718      * @param r Ring
 719      *
 720      * @return Histogram for detector,ring or nul
 721      */
 722     TH2D* Get(UShort_t d, Char_t r) const;
 723     TH2D* fFMD1i; // Histogram for FMD1i
 724     TH2D* fFMD2i; // Histogram for FMD2i
 725     TH2D* fFMD2o; // Histogram for FMD2o
 726     TH2D* fFMD3i; // Histogram for FMD3i
 727     TH2D* fFMD3o; // Histogram for FMD3o
 728
 729     ClassDef(Histos,2)
 730   };
 731
 732   //__________________________________________________________________
 733   /**
 734    * Base class for structure holding ring specific histograms
 735    *
 736    * @ingroup pwglf_forward
 737    */
 738   struct RingHistos : public TObject
 739   {
 740     /**
 741      * Constructor
 742      *
 743      */
 744     RingHistos() : fDet(0), fRing('\0'), fName(""), fkNSector(0), fkNStrip(0) {}
 745     /**
 746      *
 747      *
 748      * @param d Detector
 749      * @param r Ring
 750      */
 751     RingHistos(UShort_t d, Char_t r)
 752       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)),
 753         fkNSector(r == 'i' || r == 'I' ? 20 : 40),
 754         fkNStrip(r == 'i' || r == 'I' ? 512 : 256)
 755     {}
 756     /**
 757      * Copy constructor
 758      *
 759      * @param o Object to copy from
 760      */
 761     RingHistos(const RingHistos& o)
 762       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName),
 763         fkNSector(o.fkNSector), fkNStrip(o.fkNStrip)
 764     {}
 765     /**
 766      *
 767      */
 768     virtual ~RingHistos() {}
 769     /**
 770      * Assignement operator
 771      *
 772      * @param o Object to assign from
 773      *
 774      * @return Reference to this
 775      */
 776     RingHistos& operator=(const RingHistos& o)
 777     {
 778       if (&o == this) return *this;
 779       TObject::operator=(o);
 780       fDet  = o.fDet;
 781       fRing = o.fRing;
 782       fName = o.fName;
 783       fkNSector = o.fkNSector;
 784       fkNStrip  = o.fkNStrip;
 785       return *this;
 786     }
 787     /**
 788      * Define the outout list in @a d
 789      *
 790      * @param d Where to put the output list
 791      *
 792      * @return Newly allocated TList object or null
 793      */
 794     TList* DefineOutputList(TList* d) const;
 795     /**
 796      * Get our output list from the container @a d
 797      *
 798      * @param d where to get the output list from
 799      *
 800      * @return The found TList or null
 801      */
 802     TList* GetOutputList(const TList* d) const;
 803     /**
 804      * Find a specific histogram in the source list @a d
 805      *
 806      * @param d     (top)-container
 807      * @param name  Name of histogram
 808      *
 809      * @return Found histogram or null
 810      */
 811     TH1* GetOutputHist(const TList* d, const char* name) const;
 812     /**
 813      * Get the colour of this ring
 814      *
 815      *
 816      * @return
 817      */
 818     Color_t Color() const
 819     {
 820       return AliForwardUtil::RingColor(fDet, fRing);
 821     }
 822     /**
 823      * The name of this ring
 824      *
 825      * @return Name of this ring
 826      */
 827     const char* GetName() const { return fName.Data(); }
 828     /**
 829      * Get number of sectors
 830      */
 831     const UShort_t& NSector() const { return fkNSector; }
 832     /**
 833      * Get number of strips
 834      */
 835     const UShort_t& NStrip() const { return fkNStrip; }
 836     UShort_t fDet;   // Detector
 837     Char_t   fRing;  // Ring
 838     TString  fName;  // Name
 839     UShort_t fkNSector; // Number of sectors
 840     UShort_t fkNStrip;  // Number of strips
 841
 842     ClassDef(RingHistos,1)
 843   };
 844   /* @} */
 845
 846   //__________________________________________________________________
 847   /**
 848    * A guard idom for producing debug output
 849    *
 850    */
 851   struct DebugGuard
 852   {
 853     /**
 854      * Constructor
 855      *
 856      * @param lvl       Current level
 857      * @param msgLvl    Target level
 858      * @param format    @c printf -like format
 859      *
 860      * @return
 861      */
 862     DebugGuard(Int_t lvl, Int_t msgLvl, const char* format, ...);
 863     /**
 864      * Destructor
 865      */
 866     ~DebugGuard();
 867     /**
 868      * Make a message
 869      *
 870      * @param lvl    Current level
 871      * @param msgLvl Target level
 872      * @param format @c printf -like format
 873      */
 874     static void Message(Int_t lvl, Int_t msgLvl, const char* format, ...);
 875   private:
 876     /**
 877      * Output the message
 878      *
 879      * @param in    Direction
 880      * @param msg   Message
 881      */
 882     static void Output(int in, TString& msg);
 883     /**
 884      * Format a message
 885      *
 886      * @param out     Output is stored here
 887      * @param format  @c printf -like format
 888      * @param ap      List of arguments
 889      */
 890     static void Format(TString& out, const char* format, va_list ap);
 891     TString fMsg;
 892   };
 893 private:
 894   /**
 895    * Constructor
 896    */
 897   AliForwardUtil() {}
 898   /**
 899    * Copy constructor
 900    *
 901    * @param o Object to copy from
 902    */
 903   AliForwardUtil(const AliForwardUtil& o) : TObject(o) {}
 904   /**
 905    * Assingment operator
 906    *
 907    *
 908    * @return Reference to this object
 909    */
 910   AliForwardUtil& operator=(const AliForwardUtil&) { return *this; }
 911   /**
 912    * Destructor
 913    */
 914   ~AliForwardUtil() {}
 915
 916
 917   ClassDef(AliForwardUtil,1) // Utilities - do not make object
 918 };
 919
 920 // #ifdef LOG_NO_DEBUG
 921 // # define DGUARD(L,N,F,...) do {} while(false)
 922 // #else
 923 /**
 924  * Macro to declare a DebugGuard
 925  *
 926  * @param L Current debug level
 927  * @param N Target debug level
 928  * @param F @c printf -like Format
 929  */
 930 # define DGUARD(L,N,F,...)                                      \
 931   AliForwardUtil::DebugGuard _GUARD(L,N,F, ## __VA_ARGS__)
 932 /**
 933  * Macro to make a debug message, using DebugGuard::Message
 934  *
 935  * @param L Current debug level
 936  * @param N Target debug level
 937  * @param F @c printf -like Format
 938  */
 939 # define DMSG(L,N,F,...)                                        \
 940   AliForwardUtil::DebugGuard::Message(L,N,F, ## __VA_ARGS__)
 941 // #endif
 942 #endif
 943 // Local Variables:
 944 //  mode: C++
 945 // End:
 946