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