PWGLF/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 <AliAODInputHandler.h>
  11 #include <AliAODHandler.h>
  12 #include <AliAODEvent.h>
  13 #include <AliESDEvent.h>
  14 #include <AliAnalysisTaskSE.h>
  15 #include <AliPhysicsSelection.h>
  16 #include <AliTriggerAnalysis.h>
  17 #include <AliMultiplicity.h>
  18 #include <TParameter.h>
  19 #include <TH2D.h>
  20 #include <TH1I.h>
  21 #include <TF1.h>
  22 #include <TFitResult.h>
  23 #include <TMath.h>
  24 #include <TError.h>
  25 #include <TROOT.h>
  26
  27 //====================================================================
  28 UShort_t
  29 AliForwardUtil::ParseCollisionSystem(const char* sys)
  30 {
  31   //
  32   // Parse a collision system spec given in a string.   Known values are
  33   //
  34   //  - "ppb", "p-pb", "pa", "p-a"  which returns kPPb
  35   //  - "pp", "p-p"                 which returns kPP
  36   //  - "PbPb", "Pb-Pb", "A-A",     which returns kPbPb
  37   //  - Everything else gives kUnknown
  38   //
  39   // Parameters:
  40   //    sys Collision system spec
  41   //
  42   // Return:
  43   //    Collision system id
  44   //
  45   TString s(sys);
  46   s.ToLower();
  47   // we do pA first to avoid pp catch on ppb string (AH)
  48   if (s.Contains("p-pb")  || s.Contains("ppb"))   return AliForwardUtil::kPPb;
  49   if (s.Contains("p-a")   || s.Contains("pa"))    return AliForwardUtil::kPPb;
  50   if (s.Contains("a-p")   || s.Contains("ap"))    return AliForwardUtil::kPPb;
  51   if (s.Contains("p-p")   || s.Contains("pp"))    return AliForwardUtil::kPP;
  52   if (s.Contains("pb-pb") || s.Contains("pbpb"))  return AliForwardUtil::kPbPb;
  53   if (s.Contains("a-a")   || s.Contains("aa"))    return AliForwardUtil::kPbPb;
  54   return AliForwardUtil::kUnknown;
  55 }
  56 //____________________________________________________________________
  57 const char*
  58 AliForwardUtil::CollisionSystemString(UShort_t sys)
  59 {
  60   //
  61   // Get a string representation of the collision system
  62   //
  63   // Parameters:
  64   //    sys  Collision system
  65   // - kPP -> "pp"
  66   // - kPbPb -> "PbPb"
  67   // - anything else gives "unknown"
  68   //
  69   // Return:
  70   //    String representation of the collision system
  71   //
  72   switch (sys) {
  73   case AliForwardUtil::kPP:   return "pp";
  74   case AliForwardUtil::kPbPb: return "PbPb";
  75   case AliForwardUtil::kPPb:  return "pPb";
  76   }
  77   return "unknown";
  78 }
  79 //____________________________________________________________________
  80 Float_t
  81 AliForwardUtil::BeamRapidity(Float_t beam, UShort_t z, UShort_t a)
  82 {
  83   const Double_t pMass = 9.38271999999999995e-01;
  84   const Double_t nMass = 9.39564999999999984e-01;
  85   Double_t       beamE = z * beam / 2;
  86   Double_t       beamM = z * pMass + (a - z) * nMass;
  87   Double_t       beamP = TMath::Sqrt(beamE * beamE - beamM * beamM);
  88   Double_t       beamY = .5* TMath::Log((beamE+beamP) / (beamE-beamP));
  89   return beamY;
  90 }
  91 //____________________________________________________________________
  92 Float_t
  93 AliForwardUtil::CenterOfMassEnergy(Float_t beam,
  94                                    UShort_t z1,
  95                                    UShort_t a1,
  96                                    Short_t z2,
  97                                    Short_t a2)
  98 {
  99   // Calculate the center of mass energy given target/projectile
 100   // mass and charge numbers
 101   if (z2 < 0) z2 = z1;
 102   if (a2 < 0) a2 = a1;
 103   return TMath::Sqrt(Float_t(z1*z2)/a1/a2) * beam;
 104 }
 105 //____________________________________________________________________
 106 Float_t
 107 AliForwardUtil::CenterOfMassRapidity(UShort_t z1,
 108                                      UShort_t a1,
 109                                      Short_t z2,
 110                                      Short_t a2)
 111 {
 112   // Calculate the center of mass rapidity (shift) given target/projectile
 113   // mass and charge numbers
 114   if (z2 < 0) z2 = z1;
 115   if (a2 < 0) a2 = a1;
 116   if (z2 == z1 && a2 == a1) return 0;
 117   return .5 * TMath::Log(Float_t(z1*a2)/z2/a1);
 118 }
 119
 120 //____________________________________________________________________
 121 UShort_t
 122 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t beam)
 123 {
 124   //
 125   // Parse the center of mass energy given as a float and return known
 126   // values as a unsigned integer
 127   //
 128   // Parameters:
 129   //    sys   Collision system (needed for AA)
 130   //    beam  Center of mass energy * total charge
 131   //
 132   // Return:
 133   //    Center of mass energy per nucleon
 134   //
 135   Float_t energy = beam;
 136   // Below no longer needed apparently
 137   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
 138   if (sys == AliForwardUtil::kPPb)
 139     energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
 140   if (TMath::Abs(energy - 900.)   < 10)  return 900;
 141   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
 142   if (TMath::Abs(energy - 2750.)  < 20)  return 2750;
 143   if (TMath::Abs(energy - 4400.)  < 10)  return 4400;
 144   if (TMath::Abs(energy - 5022.)  < 10)  return 5023;
 145   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
 146   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
 147   if (TMath::Abs(energy - 8000.)  < 10)  return 8000;
 148   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
 149   if (TMath::Abs(energy - 14000.) < 10)  return 14000;
 150   return 0;
 151 }
 152 //____________________________________________________________________
 153 const char*
 154 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
 155 {
 156   //
 157   // Get a string representation of the center of mass energy per nuclean
 158   //
 159   // Parameters:
 160   //    cms  Center of mass energy per nucleon
 161   //
 162   // Return:
 163   //    String representation of the center of mass energy per nuclean
 164   //
 165   return Form("%04dGeV", cms);
 166 }
 167 //____________________________________________________________________
 168 Short_t
 169 AliForwardUtil::ParseMagneticField(Float_t v)
 170 {
 171   //
 172   // Parse the magnetic field (in kG) as given by a floating point number
 173   //
 174   // Parameters:
 175   //    field  Magnetic field in kG
 176   //
 177   // Return:
 178   //    Short integer value of magnetic field in kG
 179   //
 180   if (TMath::Abs(v - 5.) < 1 ) return +5;
 181   if (TMath::Abs(v + 5.) < 1 ) return -5;
 182   if (TMath::Abs(v) < 1)       return 0;
 183   return 999;
 184 }
 185 //____________________________________________________________________
 186 const char*
 187 AliForwardUtil::MagneticFieldString(Short_t f)
 188 {
 189   //
 190   // Get a string representation of the magnetic field
 191   //
 192   // Parameters:
 193   //    field Magnetic field in kG
 194   //
 195   // Return:
 196   //    String representation of the magnetic field
 197   //
 198   return Form("%01dkG", f);
 199 }
 200 //_____________________________________________________________________
 201 AliAODEvent* AliForwardUtil::GetAODEvent(AliAnalysisTaskSE* task)
 202 {
 203   // Check if AOD is the output event
 204   if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
 205
 206   AliAODEvent* ret = task->AODEvent();
 207   if (ret) return ret;
 208
 209   // Check if AOD is the input event
 210   ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
 211   if (!ret) ::Warning("GetAODEvent", "No AOD event found");
 212
 213   return ret;
 214 }
 215 //_____________________________________________________________________
 216 UShort_t AliForwardUtil::CheckForAOD()
 217 {
 218   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 219   if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
 220     ::Info("CheckForAOD", "Found AOD Input handler");
 221     return 1;
 222   }
 223   if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
 224     ::Info("CheckForAOD", "Found AOD Output handler");
 225     return 2;
 226   }
 227
 228   ::Warning("CheckForAOD",
 229             "Neither and input nor output AOD handler is specified");
 230   return 0;
 231 }
 232 //_____________________________________________________________________
 233 Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
 234 {
 235   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 236   if (!cls) {
 237     AliAnalysisTask* t = am->GetTask(clsOrName);
 238     if (!t) {
 239       ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
 240       return false;
 241     }
 242     ::Info("CheckForTask", "Found task %s", clsOrName);
 243     return true;
 244   }
 245   TClass* dep = gROOT->GetClass(clsOrName);
 246   if (!dep) {
 247     ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
 248     return false;
 249   }
 250   TIter next(am->GetTasks());
 251   TObject* o = 0;
 252   while ((o = next())) {
 253     if (o->IsA()->InheritsFrom(dep)) {
 254       ::Info("CheckForTask", "Found task of class %s: %s",
 255              clsOrName, o->GetName());
 256       return true;
 257     }
 258   }
 259   ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
 260   return false;
 261 }
 262
 263 //_____________________________________________________________________
 264 TObject* AliForwardUtil::MakeParameter(const Char_t* name, UShort_t value)
 265 {
 266   TParameter<int>* ret = new TParameter<int>(name, value);
 267   ret->SetUniqueID(value);
 268   return ret;
 269 }
 270 //_____________________________________________________________________
 271 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Int_t value)
 272 {
 273   TParameter<int>* ret = new TParameter<int>(name, value);
 274   ret->SetUniqueID(value);
 275   return ret;
 276 }
 277 //_____________________________________________________________________
 278 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Double_t value)
 279 {
 280   TParameter<double>* ret = new TParameter<double>(name, value);
 281   Float_t v = value;
 282   UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
 283   ret->SetUniqueID(*tmp);
 284   return ret;
 285 }
 286 //_____________________________________________________________________
 287 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Bool_t value)
 288 {
 289   TParameter<bool>* ret = new TParameter<bool>(name, value);
 290   ret->SetUniqueID(value);
 291   return ret;
 292 }
 293
 294 //_____________________________________________________________________
 295 void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
 296 {
 297   if (!o) return;
 298   value = o->GetUniqueID();
 299 }
 300 //_____________________________________________________________________
 301 void AliForwardUtil::GetParameter(TObject* o, Int_t& value)
 302 {
 303   if (!o) return;
 304   value = o->GetUniqueID();
 305 }
 306 //_____________________________________________________________________
 307 void AliForwardUtil::GetParameter(TObject* o, Double_t& value)
 308 {
 309   if (!o) return;
 310   UInt_t  i = o->GetUniqueID();
 311   Float_t v = *reinterpret_cast<Float_t*>(&i);
 312   value = v;
 313 }
 314 //_____________________________________________________________________
 315 void AliForwardUtil::GetParameter(TObject* o, Bool_t& value)
 316 {
 317   if (!o) return;
 318   value = o->GetUniqueID();
 319 }
 320
 321 //_____________________________________________________________________
 322 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Double_t zvtx)
 323 {
 324   //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta but support displaced vertices)
 325
 326   //Get max R of ring
 327   Double_t maxR = 0;
 328   Double_t minR = 0;
 329   Bool_t inner = false;
 330   switch (ring) {
 331   case 'i': case 'I': maxR = 17.2; minR =  4.5213; inner = true;  break;
 332   case 'o': case 'O': maxR = 28.0; minR = 15.4;    inner = false; break;
 333   default:
 334     return -99999;
 335   }
 336
 337   Double_t   rad       =  maxR- minR;
 338   Double_t   nStrips   = (ring == 'I' ? 512 : 256);
 339   Double_t   segment   = rad / nStrips;
 340   Double_t   r         =  minR + segment*strip;
 341   Int_t hybrid = sec / 2;
 342
 343   Double_t z = 0;
 344   switch (det) {
 345   case 1: z = 320.266; break;
 346   case 2: z = (inner ? 83.666 : 74.966); break;
 347   case 3: z = (inner ? -63.066 : -74.966); break;
 348   default: return -999999;
 349   }
 350   if ((hybrid % 2) == 0) z -= .5;
 351
 352   Double_t   theta = TMath::ATan2(r,z-zvtx);
 353   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 354
 355   return eta;
 356 }
 357
 358
 359 //====================================================================
 360 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
 361 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
 362 namespace {
 363   //
 364   // The shift of the most probable value for the ROOT function TMath::Landau
 365   //
 366   const Double_t  mpshift  = -0.22278298;
 367   //
 368   // Integration normalisation
 369   //
 370   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
 371
 372   //
 373   // Utility function to use in TF1 defintition
 374   //
 375   Double_t landauGaus1(Double_t* xp, Double_t* pp)
 376   {
 377     Double_t x        = xp[0];
 378     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
 379     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
 380     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
 381     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
 382     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
 383
 384     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
 385   }
 386
 387   //
 388   // Utility function to use in TF1 defintition
 389   //
 390   Double_t landauGausN(Double_t* xp, Double_t* pp)
 391   {
 392     Double_t  x        = xp[0];
 393     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 394     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 395     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 396     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 397     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 398     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 399     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
 400
 401     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
 402                                                   n, a);
 403   }
 404   //
 405   // Utility function to use in TF1 defintition
 406   //
 407   Double_t landauGausI(Double_t* xp, Double_t* pp)
 408   {
 409     Double_t x         = xp[0];
 410     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 411     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 412     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 413     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 414     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 415     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 416
 417     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 418   }
 419
 420
 421 }
 422 //____________________________________________________________________
 423 Double_t
 424 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
 425 {
 426   //
 427   // Calculate the shifted Landau
 428   // @f[
 429   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 430   // @f]
 431   //
 432   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 433   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 434   // @f$\Delta=0,\xi=1@f$.
 435   //
 436   // Parameters:
 437   //    x      Where to evaluate @f$ f'_{L}@f$
 438   //    delta  Most probable value
 439   //    xi     The 'width' of the distribution
 440   //
 441   // Return:
 442   //    @f$ f'_{L}(x;\Delta,\xi) @f$
 443   //
 444   return TMath::Landau(x, delta - xi * mpshift, xi);
 445 }
 446 //____________________________________________________________________
 447 Double_t
 448 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
 449                            Double_t sigma, Double_t sigmaN)
 450 {
 451   //
 452   // Calculate the value of a Landau convolved with a Gaussian
 453   //
 454   // @f[
 455   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 456   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 457   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 458   // @f]
 459   //
 460   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 461   // energy loss, @f$ \xi@f$ the width of the Landau, and
 462   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 463   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 464   // noise in the detector.
 465   //
 466   // Note that this function uses the constants fgConvolutionSteps and
 467   // fgConvolutionNSigma
 468   //
 469   // References:
 470   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 471   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 472   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 473   //
 474   // Parameters:
 475   //    x         where to evaluate @f$ f@f$
 476   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 477   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 478   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 479   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 480   //
 481   // Return:
 482   //    @f$ f@f$ evaluated at @f$ x@f$.
 483   //
 484   Double_t deltap = delta - xi * mpshift;
 485   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
 486   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
 487   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
 488   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
 489   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
 490   Double_t sum    = 0;
 491
 492   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
 493     Double_t x1 = xlow  + (i - .5) * step;
 494     Double_t x2 = xhigh - (i - .5) * step;
 495
 496     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
 497     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
 498   }
 499   return step * sum * invSq2pi / sigma1;
 500 }
 501
 502 //____________________________________________________________________
 503 Double_t
 504 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 505                             Double_t sigma, Double_t sigmaN, Int_t i)
 506 {
 507   //
 508   // Evaluate
 509   // @f[
 510   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 511   // @f]
 512   // corresponding to @f$ i@f$ particles i.e., with the substitutions
 513   // @f{eqnarray*}{
 514   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
 515   //    \xi       \rightarrow \xi_i       &=& i \xi
 516   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
 517   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 518   // @f}
 519   //
 520   // Parameters:
 521   //    x        Where to evaluate
 522   //    delta    @f$ \Delta@f$
 523   //    xi       @f$ \xi@f$
 524   //    sigma    @f$ \sigma@f$
 525   //    sigma_n  @f$ \sigma_n@f$
 526   //    i        @f$ i @f$
 527   //
 528   // Return:
 529   //    @f$ f_i @f$ evaluated
 530   //
 531   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
 532   Double_t xiI    =  i * xi;
 533   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
 534   if (sigmaI < 1e-10) {
 535     // Fall back to landau
 536     return AliForwardUtil::Landau(x, deltaI, xiI);
 537   }
 538   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
 539 }
 540
 541 //____________________________________________________________________
 542 Double_t
 543 AliForwardUtil::IdLandauGausdPar(Double_t x,
 544                                  UShort_t par,   Double_t dPar,
 545                                  Double_t delta, Double_t xi,
 546                                  Double_t sigma, Double_t sigmaN,
 547                                  Int_t    i)
 548 {
 549   //
 550   // Numerically evaluate
 551   // @f[
 552   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 553   // @f]
 554   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 555   // of the parameters is given by
 556   //
 557   // - 0: @f$\Delta@f$
 558   // - 1: @f$\xi@f$
 559   // - 2: @f$\sigma@f$
 560   // - 3: @f$\sigma_n@f$
 561   //
 562   // This is the partial derivative with respect to the parameter of
 563   // the response function corresponding to @f$ i@f$ particles i.e.,
 564   // with the substitutions
 565   // @f[
 566   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
 567   //    \xi       \rightarrow \xi_i       = i \xi
 568   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
 569   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 570   // @f]
 571   //
 572   // Parameters:
 573   //    x        Where to evaluate
 574   //    ipar     Parameter number
 575   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 576   //    delta    @f$ \Delta@f$
 577   //    xi       @f$ \xi@f$
 578   //    sigma    @f$ \sigma@f$
 579   //    sigma_n  @f$ \sigma_n@f$
 580   //    i        @f$ i@f$
 581   //
 582   // Return:
 583   //    @f$ f_i@f$ evaluated
 584   //
 585   if (dPar == 0) return 0;
 586   Double_t dp      = dPar;
 587   Double_t d2      = dPar / 2;
 588   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
 589   Double_t xiI     =  i * xi;
 590   Double_t si      =  TMath::Sqrt(Double_t(i));
 591   Double_t sigmaI  =  si*sigma;
 592   Double_t y1      = 0;
 593   Double_t y2      = 0;
 594   Double_t y3      = 0;
 595   Double_t y4      = 0;
 596   switch (par) {
 597   case 0:
 598     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
 599     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
 600     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
 601     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
 602     break;
 603   case 1:
 604     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
 605     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
 606     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
 607     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
 608     break;
 609   case 2:
 610     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
 611     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
 612     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
 613     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
 614     break;
 615   case 3:
 616     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
 617     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
 618     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
 619     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
 620     break;
 621   default:
 622     return 0;
 623   }
 624
 625   Double_t d0  = y1 - y4;
 626   Double_t d1  = 2 * (y2 - y3);
 627
 628   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
 629
 630   return g;
 631 }
 632
 633 //____________________________________________________________________
 634 Double_t
 635 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 636                             Double_t sigma, Double_t sigmaN, Int_t n,
 637                             const Double_t* a)
 638 {
 639   //
 640   // Evaluate
 641   // @f[
 642   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 643   // @f]
 644   //
 645   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 646   // Landau with a Gaussian (see LandauGaus).  Note that
 647   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 648   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 649   //
 650   // References:
 651   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 652   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 653   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 654   //
 655   // Parameters:
 656   //    x        Where to evaluate @f$ f_N@f$
 657   //    delta    @f$ \Delta_1@f$
 658   //    xi       @f$ \xi_1@f$
 659   //    sigma    @f$ \sigma_1@f$
 660   //    sigma_n  @f$ \sigma_n@f$
 661   //    n        @f$ N@f$ in the sum above.
 662   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 663   //                 @f$ i > 1@f$
 664   //
 665   // Return:
 666   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$
 667   //
 668   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
 669   for (Int_t i = 2; i <= n; i++)
 670     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 671   return result;
 672 }
 673 namespace {
 674   const Int_t kColors[] = { kRed+1,
 675                             kPink+3,
 676                             kMagenta+2,
 677                             kViolet+2,
 678                             kBlue+1,
 679                             kAzure+3,
 680                             kCyan+1,
 681                             kTeal+2,
 682                             kGreen+2,
 683                             kSpring+3,
 684                             kYellow+2,
 685                             kOrange+2 };
 686 }
 687
 688 //____________________________________________________________________
 689 TF1*
 690 AliForwardUtil::MakeNLandauGaus(Double_t  c,
 691                                 Double_t  delta, Double_t xi,
 692                                 Double_t  sigma, Double_t sigmaN, Int_t n,
 693                                 const Double_t* a,
 694                                 Double_t  xmin, Double_t xmax)
 695 {
 696   //
 697   // Generate a TF1 object of @f$ f_N@f$
 698   //
 699   // Parameters:
 700   //    c         Constant
 701   //    delta     @f$ \Delta@f$
 702   //    xi            @f$ \xi_1@f$
 703   //    sigma     @f$ \sigma_1@f$
 704   //    sigma_n   @f$ \sigma_n@f$
 705   //    n             @f$ N@f$ - how many particles to sum to
 706   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 707   //                  @f$ i > 1@f$
 708   //    xmin      Least value of range
 709   //    xmax      Largest value of range
 710   //
 711   // Return:
 712   //    Newly allocated TF1 object
 713   //
 714   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
 715   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
 716   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
 717   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
 718   landaun->SetLineWidth(2);
 719   landaun->SetNpx(500);
 720   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
 721
 722   // Set the initial parameters from the seed fit
 723   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 724   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 725   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 726   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 727   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
 728   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);
 729
 730   // Set the range and name of the scale parameters
 731   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 732     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
 733     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
 734   }
 735   return landaun;
 736 }
 737 //____________________________________________________________________
 738 TF1*
 739 AliForwardUtil::MakeILandauGaus(Double_t  c,
 740                                 Double_t  delta, Double_t xi,
 741                                 Double_t  sigma, Double_t sigmaN, Int_t i,
 742                                 Double_t  xmin, Double_t xmax)
 743 {
 744   //
 745   // Generate a TF1 object of @f$ f_I@f$
 746   //
 747   // Parameters:
 748   //    c        Constant
 749   //    delta    @f$ \Delta@f$
 750   //    xi       @f$ \xi_1@f$
 751   //    sigma    @f$ \sigma_1@f$
 752   //    sigma_n  @f$ \sigma_n@f$
 753   //    i            @f$ i@f$ - the number of particles
 754   //    xmin     Least value of range
 755   //    xmax     Largest value of range
 756   //
 757   // Return:
 758   //    Newly allocated TF1 object
 759   //
 760   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
 761   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
 762   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
 763   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
 764   landaui->SetLineWidth(1);
 765   landaui->SetNpx(500);
 766   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
 767
 768   // Set the initial parameters from the seed fit
 769   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 770   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 771   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 772   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 773   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
 774   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);
 775
 776   return landaui;
 777 }
 778
 779 //====================================================================
 780 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
 781                                          Double_t maxRange,
 782                                          UShort_t minusBins)
 783   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
 784     fFitResults(0), fFunctions(0)
 785 {
 786   //
 787   // Constructor
 788   //
 789   // Parameters:
 790   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
 791   //    maxRange   Maximum range to fit to
 792   //    minusBins  The number of bins below maximum to use
 793   //
 794   fFitResults.SetOwner();
 795   fFunctions.SetOwner();
 796 }
 797 //____________________________________________________________________
 798 AliForwardUtil::ELossFitter::~ELossFitter()
 799 {
 800   //
 801   // Destructor
 802   //
 803   //
 804   fFitResults.Delete();
 805   fFunctions.Delete();
 806 }
 807 //____________________________________________________________________
 808 void
 809 AliForwardUtil::ELossFitter::Clear()
 810 {
 811   //
 812   // Clear internal arrays
 813   //
 814   //
 815   fFitResults.Clear();
 816   fFunctions.Clear();
 817 }
 818 //____________________________________________________________________
 819 TF1*
 820 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
 821 {
 822   //
 823   // Fit a 1-particle signal to the passed energy loss distribution
 824   //
 825   // Note that this function clears the internal arrays first
 826   //
 827   // Parameters:
 828   //    dist    Data to fit the function to
 829   //    sigman If larger than zero, the initial guess of the
 830   //               detector induced noise. If zero or less, then this
 831   //               parameter is ignored in the fit (fixed at 0)
 832   //
 833   // Return:
 834   //    The function fitted to the data
 835   //
 836
 837   // Clear the cache
 838   Clear();
 839
 840   // Find the fit range
 841   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
 842
 843   // Get the bin with maximum
 844   Int_t    maxBin = dist->GetMaximumBin();
 845   Double_t maxE   = dist->GetBinLowEdge(maxBin);
 846
 847   // Get the low edge
 848   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
 849   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
 850   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
 851   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
 852
 853   // Restore the range
 854   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
 855
 856   // Define the function to fit
 857   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
 858
 859   // Set initial guesses, parameter names, and limits
 860   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
 861   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
 862   landau1->SetNpx(500);
 863   landau1->SetParLimits(kDelta, minE, fMaxRange);
 864   landau1->SetParLimits(kXi,    0.00, fMaxRange);
 865   landau1->SetParLimits(kSigma, 0.01, 0.1);
 866   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
 867   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
 868
 869   // Do the fit, getting the result object
 870   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
 871   landau1->SetRange(minE, fMaxRange);
 872   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
 873   fFunctions.AddAtAndExpand(landau1, 0);
 874
 875   return landau1;
 876 }
 877 //____________________________________________________________________
 878 TF1*
 879 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
 880                                           Double_t sigman)
 881 {
 882   //
 883   // Fit a N-particle signal to the passed energy loss distribution
 884   //
 885   // If there's no 1-particle fit present, it does that first
 886   //
 887   // Parameters:
 888   //    dist   Data to fit the function to
 889   //    n      Number of particle signals to fit
 890   //    sigman If larger than zero, the initial guess of the
 891   //               detector induced noise. If zero or less, then this
 892   //               parameter is ignored in the fit (fixed at 0)
 893   //
 894   // Return:
 895   //    The function fitted to the data
 896   //
 897
 898   // Get the seed fit result
 899   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
 900   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
 901   if (!r || !f) {
 902     f = Fit1Particle(dist, sigman);
 903     r = static_cast<TFitResult*>(fFitResults.At(0));
 904     if (!r || !f) {
 905       ::Warning("FitNLandau", "No first shot at landau fit");
 906       return 0;
 907     }
 908   }
 909
 910   // Get some parameters from seed fit
 911   Double_t delta1  = r->Parameter(kDelta);
 912   Double_t xi1     = r->Parameter(kXi);
 913   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
 914   Double_t minE    = f->GetXmin();
 915
 916   // Array of weights
 917   TArrayD a(n-1);
 918   for (UShort_t i = 2; i <= n; i++)
 919     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
 920   // Make the fit function
 921   TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
 922                                      r->Parameter(kDelta),
 923                                      r->Parameter(kXi),
 924                                      r->Parameter(kSigma),
 925                                      r->Parameter(kSigmaN),
 926                                      n,a.fArray,minE,maxEi);
 927   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
 928   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
 929   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
 930   // Check if we're using the noise sigma
 931   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
 932   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
 933
 934   // Set the range and name of the scale parameters
 935   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 936     landaun->SetParLimits(kA+i-2, 0,1);
 937   }
 938
 939   // Do the fit
 940   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
 941
 942   landaun->SetRange(minE, fMaxRange);
 943   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
 944   fFunctions.AddAtAndExpand(landaun, n-1);
 945
 946   return landaun;
 947 }
 948
 949 //====================================================================
 950 AliForwardUtil::Histos::~Histos()
 951 {
 952   //
 953   // Destructor
 954   //
 955 }
 956
 957 //____________________________________________________________________
 958 void
 959 AliForwardUtil::Histos::Delete(Option_t* opt)
 960 {
 961   if (fFMD1i) delete fFMD1i;
 962   if (fFMD2i) delete fFMD2i;
 963   if (fFMD2o) delete fFMD2o;
 964   if (fFMD3i) delete fFMD3i;
 965   if (fFMD3o) delete fFMD3o;
 966   fFMD1i = 0;
 967   fFMD2i = 0;
 968   fFMD2o = 0;
 969   fFMD3i = 0;
 970   fFMD3o = 0;
 971   TObject::Delete(opt);
 972 }
 973
 974 //____________________________________________________________________
 975 TH2D*
 976 AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
 977                              const TAxis& etaAxis) const
 978 {
 979   //
 980   // Make a histogram
 981   //
 982   // Parameters:
 983   //    d        Detector
 984   //    r        Ring
 985   //    etaAxis  Eta axis to use
 986   //
 987   // Return:
 988   //    Newly allocated histogram
 989   //
 990   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
 991   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r),
 992                         Form("FMD%d%c cache", d, r),
 993                         etaAxis.GetNbins(), etaAxis.GetXmin(),
 994                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
 995   hist->SetXTitle("#eta");
 996   hist->SetYTitle("#phi [radians]");
 997   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
 998   hist->Sumw2();
 999   hist->SetDirectory(0);
1000
1001   return hist;
1002 }
1003 //____________________________________________________________________
1004 void
1005 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
1006 {
1007   //
1008   // Initialize the object
1009   //
1010   // Parameters:
1011   //    etaAxis Eta axis to use
1012   //
1013   fFMD1i = Make(1, 'I', etaAxis);
1014   fFMD2i = Make(2, 'I', etaAxis);
1015   fFMD2o = Make(2, 'O', etaAxis);
1016   fFMD3i = Make(3, 'I', etaAxis);
1017   fFMD3o = Make(3, 'O', etaAxis);
1018 }
1019 //____________________________________________________________________
1020 void
1021 AliForwardUtil::Histos::Clear(Option_t* option)
1022 {
1023   //
1024   // Clear data
1025   //
1026   // Parameters:
1027   //    option Not used
1028   //
1029   if (fFMD1i) fFMD1i->Reset(option);
1030   if (fFMD2i) fFMD2i->Reset(option);
1031   if (fFMD2o) fFMD2o->Reset(option);
1032   if (fFMD3i) fFMD3i->Reset(option);
1033   if (fFMD3o) fFMD3o->Reset(option);
1034 }
1035
1036 //____________________________________________________________________
1037 TH2D*
1038 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
1039 {
1040   //
1041   // Get the histogram for a particular detector,ring
1042   //
1043   // Parameters:
1044   //    d Detector
1045   //    r Ring
1046   //
1047   // Return:
1048   //    Histogram for detector,ring or nul
1049   //
1050   switch (d) {
1051   case 1: return fFMD1i;
1052   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1053   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1054   }
1055   return 0;
1056 }
1057 //====================================================================
1058 TList*
1059 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
1060 {
1061   //
1062   // Define the outout list in @a d
1063   //
1064   // Parameters:
1065   //    d Where to put the output list
1066   //
1067   // Return:
1068   //    Newly allocated TList object or null
1069   //
1070   if (!d) return 0;
1071   TList* list = new TList;
1072   list->SetOwner();
1073   list->SetName(fName.Data());
1074   d->Add(list);
1075   return list;
1076 }
1077 //____________________________________________________________________
1078 TList*
1079 AliForwardUtil::RingHistos::GetOutputList(const TList* d) const
1080 {
1081   //
1082   // Get our output list from the container @a d
1083   //
1084   // Parameters:
1085   //    d where to get the output list from
1086   //
1087   // Return:
1088   //    The found TList or null
1089   //
1090   if (!d) return 0;
1091   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1092   return list;
1093 }
1094
1095 //____________________________________________________________________
1096 TH1*
1097 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1098 {
1099   //
1100   // Find a specific histogram in the source list @a d
1101   //
1102   // Parameters:
1103   //    d     (top)-container
1104   //    name  Name of histogram
1105   //
1106   // Return:
1107   //    Found histogram or null
1108   //
1109   return static_cast<TH1*>(d->FindObject(name));
1110 }
1111
1112 //====================================================================
1113 AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl,
1114                                        const char* format, ...)
1115   : fMsg("")
1116 {
1117   if (lvl < msgLvl) return;
1118   va_list ap;
1119   va_start(ap, format);
1120   Format(fMsg, format, ap);
1121   va_end(ap);
1122   Output(+1, fMsg);
1123 }
1124 //____________________________________________________________________
1125 AliForwardUtil::DebugGuard::~DebugGuard()
1126 {
1127   if (fMsg.IsNull()) return;
1128   Output(-1, fMsg);
1129 }
1130 //____________________________________________________________________
1131 void
1132 AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl,
1133                                     const char* format, ...)
1134 {
1135   if (lvl < msgLvl) return;
1136   TString msg;
1137   va_list ap;
1138   va_start(ap, format);
1139   Format(msg, format, ap);
1140   va_end(ap);
1141   Output(0, msg);
1142 }
1143
1144 //____________________________________________________________________
1145 void
1146 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1147 {
1148   static char buf[512];
1149   Int_t n = gROOT->GetDirLevel() + 2;
1150   for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1151   vsnprintf(&(buf[n]), 511-n, format, ap);
1152   buf[511] = '\0';
1153   out = buf;
1154 }
1155 //____________________________________________________________________
1156 void
1157 AliForwardUtil::DebugGuard::Output(int in, TString& msg)
1158 {
1159   msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
1160   AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1161   if      (in > 0) gROOT->IncreaseDirLevel();
1162   else if (in < 0) gROOT->DecreaseDirLevel();
1163 }
1164
1165
1166
1167 //
1168 // EOF
1169 //