PWG2/FORWARD/analysis2/AliForwardUtil.cxx

   1 //
   2 // Utilities used in the forward multiplcity analysis
   3 //
   4 //
   5 #include "AliForwardUtil.h"
   6 #include <AliAnalysisManager.h>
   7 #include "AliAODForwardMult.h"
   8 #include <AliLog.h>
   9 #include <AliInputEventHandler.h>
  10 #include <AliESDEvent.h>
  11 #include <AliPhysicsSelection.h>
  12 #include <AliTriggerAnalysis.h>
  13 #include <AliMultiplicity.h>
  14 #include <TH2D.h>
  15 #include <TH1I.h>
  16 #include <TF1.h>
  17 #include <TFitResult.h>
  18 #include <TMath.h>
  19 #include <TError.h>
  20
  21 //====================================================================
  22 UShort_t
  23 AliForwardUtil::ParseCollisionSystem(const char* sys)
  24 {
  25   //
  26   // Parse a collision system spec given in a string.   Known values are
  27   //
  28   //  - "pp", "p-p" which returns kPP
  29   //  - "PbPb", "Pb-Pb", "A-A", which returns kPbPb
  30   //  - Everything else gives kUnknown
  31   //
  32   // Parameters:
  33   //    sys Collision system spec
  34   //
  35   // Return:
  36   //    Collision system id
  37   //
  38   TString s(sys);
  39   s.ToLower();
  40   if (s.Contains("p-p")   || s.Contains("pp"))   return AliForwardUtil::kPP;
  41   if (s.Contains("pb-pb") || s.Contains("pbpb")) return AliForwardUtil::kPbPb;
  42   if (s.Contains("a-a")   || s.Contains("aa"))   return AliForwardUtil::kPbPb;
  43   return AliForwardUtil::kUnknown;
  44 }
  45 //____________________________________________________________________
  46 const char*
  47 AliForwardUtil::CollisionSystemString(UShort_t sys)
  48 {
  49   //
  50   // Get a string representation of the collision system
  51   //
  52   // Parameters:
  53   //    sys  Collision system
  54   // - kPP -> "pp"
  55   // - kPbPb -> "PbPb"
  56   // - anything else gives "unknown"
  57   //
  58   // Return:
  59   //    String representation of the collision system
  60   //
  61   switch (sys) {
  62   case AliForwardUtil::kPP:   return "pp";
  63   case AliForwardUtil::kPbPb: return "PbPb";
  64   }
  65   return "unknown";
  66 }
  67 //____________________________________________________________________
  68 UShort_t
  69 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t v)
  70 {
  71   //
  72   // Parse the center of mass energy given as a float and return known
  73   // values as a unsigned integer
  74   //
  75   // Parameters:
  76   //    sys  Collision system (needed for AA)
  77   //    cms  Center of mass energy * total charge
  78   //
  79   // Return:
  80   //    Center of mass energy per nucleon
  81   //
  82   Float_t energy = v;
  83   if (sys != AliForwardUtil::kPP) energy = energy / 208 * 82;
  84   if (TMath::Abs(energy - 900.)   < 10)  return 900;
  85   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
  86   if (TMath::Abs(energy - 2750.)  < 10)  return 2750;
  87   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
  88   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
  89   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
  90   if (TMath::Abs(energy - 14000.) < 10)  return 14000;
  91   return 0;
  92 }
  93 //____________________________________________________________________
  94 const char*
  95 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
  96 {
  97   //
  98   // Get a string representation of the center of mass energy per nuclean
  99   //
 100   // Parameters:
 101   //    cms  Center of mass energy per nucleon
 102   //
 103   // Return:
 104   //    String representation of the center of mass energy per nuclean
 105   //
 106   return Form("%04dGeV", cms);
 107 }
 108 //____________________________________________________________________
 109 Short_t
 110 AliForwardUtil::ParseMagneticField(Float_t v)
 111 {
 112   //
 113   // Parse the magnetic field (in kG) as given by a floating point number
 114   //
 115   // Parameters:
 116   //    field  Magnetic field in kG
 117   //
 118   // Return:
 119   //    Short integer value of magnetic field in kG
 120   //
 121   if (TMath::Abs(v - 5.) < 1 ) return +5;
 122   if (TMath::Abs(v + 5.) < 1 ) return -5;
 123   if (TMath::Abs(v) < 1)       return 0;
 124   return 999;
 125 }
 126 //____________________________________________________________________
 127 const char*
 128 AliForwardUtil::MagneticFieldString(Short_t f)
 129 {
 130   //
 131   // Get a string representation of the magnetic field
 132   //
 133   // Parameters:
 134   //    field Magnetic field in kG
 135   //
 136   // Return:
 137   //    String representation of the magnetic field
 138   //
 139   return Form("%01dkG", f);
 140 }
 141
 142
 143 //====================================================================
 144 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
 145 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
 146 namespace {
 147   //
 148   // The shift of the most probable value for the ROOT function TMath::Landau
 149   //
 150   const Double_t  mpshift  = -0.22278298;
 151   //
 152   // Integration normalisation
 153   //
 154   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
 155
 156   //
 157   // Utility function to use in TF1 defintition
 158   //
 159   Double_t landauGaus1(Double_t* xp, Double_t* pp)
 160   {
 161     Double_t x        = xp[0];
 162     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
 163     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
 164     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
 165     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
 166     Double_t sigma_n  = pp[AliForwardUtil::ELossFitter::kSigmaN];
 167
 168     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigma_n);
 169   }
 170
 171   //
 172   // Utility function to use in TF1 defintition
 173   //
 174   Double_t landauGausN(Double_t* xp, Double_t* pp)
 175   {
 176     Double_t  x        = xp[0];
 177     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 178     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 179     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 180     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 181     Double_t sigma_n   = pp[AliForwardUtil::ELossFitter::kSigmaN];
 182     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 183     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
 184
 185     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
 186                                                   n, a);
 187   }
 188   //
 189   // Utility function to use in TF1 defintition
 190   //
 191   Double_t landauGausI(Double_t* xp, Double_t* pp)
 192   {
 193     Double_t x         = xp[0];
 194     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 195     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 196     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 197     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 198     Double_t sigma_n   = pp[AliForwardUtil::ELossFitter::kSigmaN];
 199     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 200
 201     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigma_n,i);
 202   }
 203
 204
 205 }
 206 //____________________________________________________________________
 207 Double_t
 208 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
 209 {
 210   //
 211   // Calculate the shifted Landau
 212   // @f[
 213   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 214   // @f]
 215   //
 216   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 217   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 218   // @f$\Delta=0,\xi=1@f$.
 219   //
 220   // Parameters:
 221   //    x      Where to evaluate @f$ f'_{L}@f$
 222   //    delta  Most probable value
 223   //    xi     The 'width' of the distribution
 224   //
 225   // Return:
 226   //    @f$ f'_{L}(x;\Delta,\xi) @f$
 227   //
 228   return TMath::Landau(x, delta - xi * mpshift, xi);
 229 }
 230 //____________________________________________________________________
 231 Double_t
 232 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
 233                            Double_t sigma, Double_t sigma_n)
 234 {
 235   //
 236   // Calculate the value of a Landau convolved with a Gaussian
 237   //
 238   // @f[
 239   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 240   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 241   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 242   // @f]
 243   //
 244   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 245   // energy loss, @f$ \xi@f$ the width of the Landau, and
 246   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 247   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 248   // noise in the detector.
 249   //
 250   // Note that this function uses the constants fgConvolutionSteps and
 251   // fgConvolutionNSigma
 252   //
 253   // References:
 254   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 255   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 256   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 257   //
 258   // Parameters:
 259   //    x         where to evaluate @f$ f@f$
 260   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 261   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 262   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 263   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 264   //
 265   // Return:
 266   //    @f$ f@f$ evaluated at @f$ x@f$.
 267   //
 268   Double_t deltap = delta - xi * mpshift;
 269   Double_t sigma2 = sigma_n*sigma_n + sigma*sigma;
 270   Double_t sigma1 = sigma_n == 0 ? sigma : TMath::Sqrt(sigma2);
 271   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
 272   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
 273   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
 274   Double_t sum    = 0;
 275
 276   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
 277     Double_t x1 = xlow  + (i - .5) * step;
 278     Double_t x2 = xhigh - (i - .5) * step;
 279
 280     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
 281     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
 282   }
 283   return step * sum * invSq2pi / sigma1;
 284 }
 285
 286 //____________________________________________________________________
 287 Double_t
 288 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 289                             Double_t sigma, Double_t sigma_n, Int_t i)
 290 {
 291   //
 292   // Evaluate
 293   // @f[
 294   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 295   // @f]
 296   // corresponding to @f$ i@f$ particles i.e., with the substitutions
 297   // @f{eqnarray*}{
 298   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
 299   //    \xi       \rightarrow \xi_i       &=& i \xi
 300   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
 301   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 302   // @f}
 303   //
 304   // Parameters:
 305   //    x        Where to evaluate
 306   //    delta    @f$ \Delta@f$
 307   //    xi       @f$ \xi@f$
 308   //    sigma    @f$ \sigma@f$
 309   //    sigma_n  @f$ \sigma_n@f$
 310   //    i        @f$ i @f$
 311   //
 312   // Return:
 313   //    @f$ f_i @f$ evaluated
 314   //
 315   Double_t delta_i =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
 316   Double_t xi_i    =  i * xi;
 317   Double_t sigma_i =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
 318   if (sigma_i < 1e-10) {
 319     // Fall back to landau
 320     return AliForwardUtil::Landau(x, delta_i, xi_i);
 321   }
 322   return AliForwardUtil::LandauGaus(x, delta_i, xi_i, sigma_i, sigma_n);
 323 }
 324
 325 //____________________________________________________________________
 326 Double_t
 327 AliForwardUtil::IdLandauGausdPar(Double_t x,
 328                                  UShort_t par,   Double_t dPar,
 329                                  Double_t delta, Double_t xi,
 330                                  Double_t sigma, Double_t sigma_n,
 331                                  Int_t    i)
 332 {
 333   //
 334   // Numerically evaluate
 335   // @f[
 336   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 337   // @f]
 338   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 339   // of the parameters is given by
 340   //
 341   // - 0: @f$\Delta@f$
 342   // - 1: @f$\xi@f$
 343   // - 2: @f$\sigma@f$
 344   // - 3: @f$\sigma_n@f$
 345   //
 346   // This is the partial derivative with respect to the parameter of
 347   // the response function corresponding to @f$ i@f$ particles i.e.,
 348   // with the substitutions
 349   // @f[
 350   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
 351   //    \xi       \rightarrow \xi_i       = i \xi
 352   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
 353   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 354   // @f]
 355   //
 356   // Parameters:
 357   //    x        Where to evaluate
 358   //    ipar     Parameter number
 359   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 360   //    delta    @f$ \Delta@f$
 361   //    xi       @f$ \xi@f$
 362   //    sigma    @f$ \sigma@f$
 363   //    sigma_n  @f$ \sigma_n@f$
 364   //    i        @f$ i@f$
 365   //
 366   // Return:
 367   //    @f$ f_i@f$ evaluated
 368   //
 369   if (dPar == 0) return 0;
 370   Double_t dp      = dPar;
 371   Double_t d2      = dPar / 2;
 372   Double_t delta_i =  i * (delta + xi * TMath::Log(i));
 373   Double_t xi_i    =  i * xi;
 374   Double_t si      =  TMath::Sqrt(Double_t(i));
 375   Double_t sigma_i =  si*sigma;
 376   Double_t y1      = 0;
 377   Double_t y2      = 0;
 378   Double_t y3      = 0;
 379   Double_t y4      = 0;
 380   switch (par) {
 381   case 0:
 382     y1 = ILandauGaus(x, delta_i+i*dp, xi_i, sigma_i, sigma_n, i);
 383     y2 = ILandauGaus(x, delta_i+i*d2, xi_i, sigma_i, sigma_n, i);
 384     y3 = ILandauGaus(x, delta_i-i*d2, xi_i, sigma_i, sigma_n, i);
 385     y4 = ILandauGaus(x, delta_i-i*dp, xi_i, sigma_i, sigma_n, i);
 386     break;
 387   case 1:
 388     y1 = ILandauGaus(x, delta_i, xi_i+i*dp, sigma_i, sigma_n, i);
 389     y2 = ILandauGaus(x, delta_i, xi_i+i*d2, sigma_i, sigma_n, i);
 390     y3 = ILandauGaus(x, delta_i, xi_i-i*d2, sigma_i, sigma_n, i);
 391     y4 = ILandauGaus(x, delta_i, xi_i-i*dp, sigma_i, sigma_n, i);
 392     break;
 393   case 2:
 394     y1 = ILandauGaus(x, delta_i, xi_i, sigma_i+si*dp, sigma_n, i);
 395     y2 = ILandauGaus(x, delta_i, xi_i, sigma_i+si*d2, sigma_n, i);
 396     y3 = ILandauGaus(x, delta_i, xi_i, sigma_i-si*d2, sigma_n, i);
 397     y4 = ILandauGaus(x, delta_i, xi_i, sigma_i-si*dp, sigma_n, i);
 398     break;
 399   case 3:
 400     y1 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n+dp, i);
 401     y2 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n+d2, i);
 402     y3 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n-d2, i);
 403     y4 = ILandauGaus(x, delta_i, xi_i, sigma_i, sigma_n-dp, i);
 404     break;
 405   default:
 406     return 0;
 407   }
 408
 409   Double_t d0  = y1 - y4;
 410   Double_t d1  = 2 * (y2 - y3);
 411
 412   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
 413
 414   return g;
 415 }
 416
 417 //____________________________________________________________________
 418 Double_t
 419 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 420                             Double_t sigma, Double_t sigma_n, Int_t n,
 421                             Double_t* a)
 422 {
 423   //
 424   // Evaluate
 425   // @f[
 426   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 427   // @f]
 428   //
 429   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 430   // Landau with a Gaussian (see LandauGaus).  Note that
 431   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 432   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 433   //
 434   // References:
 435   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 436   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 437   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 438   //
 439   // Parameters:
 440   //    x        Where to evaluate @f$ f_N@f$
 441   //    delta    @f$ \Delta_1@f$
 442   //    xi       @f$ \xi_1@f$
 443   //    sigma    @f$ \sigma_1@f$
 444   //    sigma_n  @f$ \sigma_n@f$
 445   //    n        @f$ N@f$ in the sum above.
 446   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 447   //                 @f$ i > 1@f$
 448   //
 449   // Return:
 450   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$
 451   //
 452   Double_t result = ILandauGaus(x, delta, xi, sigma, sigma_n, 1);
 453   for (Int_t i = 2; i <= n; i++)
 454     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigma_n,i);
 455   return result;
 456 }
 457 namespace {
 458   const Int_t kColors[] = { kRed+1,
 459                             kPink+3,
 460                             kMagenta+2,
 461                             kViolet+2,
 462                             kBlue+1,
 463                             kAzure+3,
 464                             kCyan+1,
 465                             kTeal+2,
 466                             kGreen+2,
 467                             kSpring+3,
 468                             kYellow+2,
 469                             kOrange+2 };
 470 }
 471
 472 //____________________________________________________________________
 473 TF1*
 474 AliForwardUtil::MakeNLandauGaus(Double_t  c,
 475                                 Double_t  delta, Double_t xi,
 476                                 Double_t  sigma, Double_t sigma_n, Int_t n,
 477                                 Double_t* a,
 478                                 Double_t  xmin, Double_t xmax)
 479 {
 480   //
 481   // Generate a TF1 object of @f$ f_N@f$
 482   //
 483   // Parameters:
 484   //    c         Constant
 485   //    delta     @f$ \Delta@f$
 486   //    xi            @f$ \xi_1@f$
 487   //    sigma     @f$ \sigma_1@f$
 488   //    sigma_n   @f$ \sigma_n@f$
 489   //    n             @f$ N@f$ - how many particles to sum to
 490   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 491   //                  @f$ i > 1@f$
 492   //    xmin      Least value of range
 493   //    xmax      Largest value of range
 494   //
 495   // Return:
 496   //    Newly allocated TF1 object
 497   //
 498   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
 499   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
 500   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
 501   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
 502   landaun->SetLineWidth(2);
 503   landaun->SetNpx(500);
 504   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
 505
 506   // Set the initial parameters from the seed fit
 507   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 508   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 509   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 510   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 511   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigma_n);
 512   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);
 513
 514   // Set the range and name of the scale parameters
 515   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 516     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
 517     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
 518   }
 519   return landaun;
 520 }
 521 //____________________________________________________________________
 522 TF1*
 523 AliForwardUtil::MakeILandauGaus(Double_t  c,
 524                                 Double_t  delta, Double_t xi,
 525                                 Double_t  sigma, Double_t sigma_n, Int_t i,
 526                                 Double_t  xmin, Double_t xmax)
 527 {
 528   //
 529   // Generate a TF1 object of @f$ f_I@f$
 530   //
 531   // Parameters:
 532   //    c        Constant
 533   //    delta    @f$ \Delta@f$
 534   //    xi       @f$ \xi_1@f$
 535   //    sigma    @f$ \sigma_1@f$
 536   //    sigma_n  @f$ \sigma_n@f$
 537   //    i            @f$ i@f$ - the number of particles
 538   //    xmin     Least value of range
 539   //    xmax     Largest value of range
 540   //
 541   // Return:
 542   //    Newly allocated TF1 object
 543   //
 544   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
 545   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
 546   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
 547   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
 548   landaui->SetLineWidth(1);
 549   landaui->SetNpx(500);
 550   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
 551
 552   // Set the initial parameters from the seed fit
 553   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 554   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 555   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 556   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 557   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigma_n);
 558   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);
 559
 560   return landaui;
 561 }
 562
 563 //====================================================================
 564 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
 565                                          Double_t maxRange,
 566                                          UShort_t minusBins)
 567   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
 568     fFitResults(0), fFunctions(0)
 569 {
 570   //
 571   // Constructor
 572   //
 573   // Parameters:
 574   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
 575   //    maxRange   Maximum range to fit to
 576   //    minusBins  The number of bins below maximum to use
 577   //
 578   fFitResults.SetOwner();
 579   fFunctions.SetOwner();
 580 }
 581 //____________________________________________________________________
 582 AliForwardUtil::ELossFitter::~ELossFitter()
 583 {
 584   //
 585   // Destructor
 586   //
 587   //
 588   fFitResults.Delete();
 589   fFunctions.Delete();
 590 }
 591 //____________________________________________________________________
 592 void
 593 AliForwardUtil::ELossFitter::Clear()
 594 {
 595   //
 596   // Clear internal arrays
 597   //
 598   //
 599   fFitResults.Clear();
 600   fFunctions.Clear();
 601 }
 602 //____________________________________________________________________
 603 TF1*
 604 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
 605 {
 606   //
 607   // Fit a 1-particle signal to the passed energy loss distribution
 608   //
 609   // Note that this function clears the internal arrays first
 610   //
 611   // Parameters:
 612   //    dist    Data to fit the function to
 613   //    sigman If larger than zero, the initial guess of the
 614   //               detector induced noise. If zero or less, then this
 615   //               parameter is ignored in the fit (fixed at 0)
 616   //
 617   // Return:
 618   //    The function fitted to the data
 619   //
 620
 621   // Clear the cache
 622   Clear();
 623
 624   // Find the fit range
 625   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
 626
 627   // Get the bin with maximum
 628   Int_t    maxBin = dist->GetMaximumBin();
 629   Double_t maxE   = dist->GetBinLowEdge(maxBin);
 630
 631   // Get the low edge
 632   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
 633   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
 634   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
 635   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
 636
 637   // Restore the range
 638   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
 639
 640   // Define the function to fit
 641   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
 642
 643   // Set initial guesses, parameter names, and limits
 644   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
 645   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
 646   landau1->SetNpx(500);
 647   landau1->SetParLimits(kDelta, minE, fMaxRange);
 648   landau1->SetParLimits(kXi,    0.00, fMaxRange);
 649   landau1->SetParLimits(kSigma, 0.01, 0.1);
 650   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
 651   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
 652
 653   // Do the fit, getting the result object
 654   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
 655   landau1->SetRange(minE, fMaxRange);
 656   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
 657   fFunctions.AddAtAndExpand(landau1, 0);
 658
 659   return landau1;
 660 }
 661 //____________________________________________________________________
 662 TF1*
 663 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
 664                                           Double_t sigman)
 665 {
 666   //
 667   // Fit a N-particle signal to the passed energy loss distribution
 668   //
 669   // If there's no 1-particle fit present, it does that first
 670   //
 671   // Parameters:
 672   //    dist   Data to fit the function to
 673   //    n      Number of particle signals to fit
 674   //    sigman If larger than zero, the initial guess of the
 675   //               detector induced noise. If zero or less, then this
 676   //               parameter is ignored in the fit (fixed at 0)
 677   //
 678   // Return:
 679   //    The function fitted to the data
 680   //
 681
 682   // Get the seed fit result
 683   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
 684   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
 685   if (!r || !f) {
 686     f = Fit1Particle(dist, sigman);
 687     r = static_cast<TFitResult*>(fFitResults.At(0));
 688     if (!r || !f) {
 689       ::Warning("FitNLandau", "No first shot at landau fit");
 690       return 0;
 691     }
 692   }
 693
 694   // Get some parameters from seed fit
 695   Double_t delta1  = r->Parameter(kDelta);
 696   Double_t xi1     = r->Parameter(kXi);
 697   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
 698   Double_t minE    = f->GetXmin();
 699
 700   // Array of weights
 701   TArrayD a(n-1);
 702   for (UShort_t i = 2; i <= n; i++)
 703     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
 704   // Make the fit function
 705   TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
 706                                      r->Parameter(kDelta),
 707                                      r->Parameter(kXi),
 708                                      r->Parameter(kSigma),
 709                                      r->Parameter(kSigmaN),
 710                                      n,a.fArray,minE,maxEi);
 711   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
 712   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
 713   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
 714   // Check if we're using the noise sigma
 715   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
 716   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
 717
 718   // Set the range and name of the scale parameters
 719   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 720     landaun->SetParLimits(kA+i-2, 0,1);
 721   }
 722
 723   // Do the fit
 724   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
 725
 726   landaun->SetRange(minE, fMaxRange);
 727   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
 728   fFunctions.AddAtAndExpand(landaun, n-1);
 729
 730   return landaun;
 731 }
 732
 733 //====================================================================
 734 AliForwardUtil::Histos::~Histos()
 735 {
 736   //
 737   // Destructor
 738   //
 739   if (fFMD1i) delete fFMD1i;
 740   if (fFMD2i) delete fFMD2i;
 741   if (fFMD2o) delete fFMD2o;
 742   if (fFMD3i) delete fFMD3i;
 743   if (fFMD3o) delete fFMD3o;
 744 }
 745
 746 //____________________________________________________________________
 747 TH2D*
 748 AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
 749                              const TAxis& etaAxis) const
 750 {
 751   //
 752   // Make a histogram
 753   //
 754   // Parameters:
 755   //    d        Detector
 756   //    r        Ring
 757   //    etaAxis  Eta axis to use
 758   //
 759   // Return:
 760   //    Newly allocated histogram
 761   //
 762   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
 763   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r),
 764                         Form("FMD%d%c cache", d, r),
 765                         etaAxis.GetNbins(), etaAxis.GetXmin(),
 766                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
 767   hist->SetXTitle("#eta");
 768   hist->SetYTitle("#phi [radians]");
 769   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 770   hist->Sumw2();
 771   hist->SetDirectory(0);
 772
 773   return hist;
 774 }
 775 //____________________________________________________________________
 776 void
 777 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
 778 {
 779   //
 780   // Initialize the object
 781   //
 782   // Parameters:
 783   //    etaAxis Eta axis to use
 784   //
 785   fFMD1i = Make(1, 'I', etaAxis);
 786   fFMD2i = Make(2, 'I', etaAxis);
 787   fFMD2o = Make(2, 'O', etaAxis);
 788   fFMD3i = Make(3, 'I', etaAxis);
 789   fFMD3o = Make(3, 'O', etaAxis);
 790 }
 791 //____________________________________________________________________
 792 void
 793 AliForwardUtil::Histos::Clear(Option_t* option)
 794 {
 795   //
 796   // Clear data
 797   //
 798   // Parameters:
 799   //    option Not used
 800   //
 801   fFMD1i->Reset(option);
 802   fFMD2i->Reset(option);
 803   fFMD2o->Reset(option);
 804   fFMD3i->Reset(option);
 805   fFMD3o->Reset(option);
 806 }
 807
 808 //____________________________________________________________________
 809 TH2D*
 810 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
 811 {
 812   //
 813   // Get the histogram for a particular detector,ring
 814   //
 815   // Parameters:
 816   //    d Detector
 817   //    r Ring
 818   //
 819   // Return:
 820   //    Histogram for detector,ring or nul
 821   //
 822   switch (d) {
 823   case 1: return fFMD1i;
 824   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
 825   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
 826   }
 827   return 0;
 828 }
 829 //====================================================================
 830 TList*
 831 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
 832 {
 833   //
 834   // Define the outout list in @a d
 835   //
 836   // Parameters:
 837   //    d Where to put the output list
 838   //
 839   // Return:
 840   //    Newly allocated TList object or null
 841   //
 842   if (!d) return 0;
 843   TList* list = new TList;
 844   list->SetName(fName.Data());
 845   d->Add(list);
 846   return list;
 847 }
 848 //____________________________________________________________________
 849 TList*
 850 AliForwardUtil::RingHistos::GetOutputList(TList* d) const
 851 {
 852   //
 853   // Get our output list from the container @a d
 854   //
 855   // Parameters:
 856   //    d where to get the output list from
 857   //
 858   // Return:
 859   //    The found TList or null
 860   //
 861   if (!d) return 0;
 862   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
 863   return list;
 864 }
 865
 866 //____________________________________________________________________
 867 TH1*
 868 AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
 869 {
 870   //
 871   // Find a specific histogram in the source list @a d
 872   //
 873   // Parameters:
 874   //    d     (top)-container
 875   //    name  Name of histogram
 876   //
 877   // Return:
 878   //    Found histogram or null
 879   //
 880   return static_cast<TH1*>(d->FindObject(name));
 881 }
 882
 883 //
 884 // EOF
 885 //