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::GetStripR(Char_t ring, UShort_t strip)
 323 {
 324   //Get max R of ring
 325   const Double_t iMinR = 4.5213;
 326   const Double_t iMaxR = 17.2;
 327   const Double_t oMinR = 15.4;
 328   const Double_t oMaxR = 28.0;
 329
 330   Double_t   minR    = (ring == 'I' || ring == 'i') ? iMinR : oMinR;
 331   Double_t   maxR    = (ring == 'I' || ring == 'i') ? iMaxR : oMaxR;
 332   Double_t   nStrips = (ring == 'I' || ring == 'i') ? 512   : 256;
 333   Double_t   rad     =  maxR - minR;
 334   Double_t   segment = rad / nStrips;
 335   Double_t   r       =  minR + segment*strip;
 336
 337   return r;
 338 }
 339
 340 //_____________________________________________________________________
 341 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring,
 342                                          UShort_t sec, UShort_t strip,
 343                                          Double_t zvtx)
 344 {
 345   // Calculate eta from strip with vertex (redundant with
 346   // AliESDFMD::Eta but support displaced vertices)
 347
 348   //Get max R of ring
 349   Double_t   r         = GetStripR(ring, strip);
 350   Int_t      hybrid    = sec / 2;
 351   Bool_t     inner     = (ring == 'I' || ring == 'i');
 352   Double_t   z         = 0;
 353   switch (det) {
 354   case 1: z = 320.266;                     break;
 355   case 2: z = (inner ?  83.666 :  74.966); break;
 356   case 3: z = (inner ? -63.066 : -74.966); break;
 357   default: return -999999;
 358   }
 359   if ((hybrid % 2) == 0) z -= .5;
 360
 361   Double_t   theta = TMath::ATan2(r,z-zvtx);
 362   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 363
 364   return eta;
 365 }
 366
 367 //_____________________________________________________________________
 368 Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip,
 369                                          Double_t phi,
 370                                          Double_t xvtx, Double_t yvtx)
 371 {
 372   // Calculate eta from strip with vertex (redundant with
 373   // AliESDFMD::Eta but support displaced vertices)
 374
 375   // Unknown x,y -> no change
 376   if (yvtx > 999 || xvtx > 999) return phi;
 377
 378   //Get max R of ring
 379   Double_t r   = GetStripR(ring, strip);
 380   Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
 381   Double_t pha = (TMath::Abs(yvtx) < 1e-12  ? 0 : TMath::ATan2(xvtx, yvtx));
 382   Double_t cha = amp * TMath::Cos(phi+pha);
 383   phi += cha;
 384   if (phi < 0)              phi += TMath::TwoPi();
 385   if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
 386   return phi;
 387 }
 388
 389 //====================================================================
 390 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
 391 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
 392 namespace {
 393   //
 394   // The shift of the most probable value for the ROOT function TMath::Landau
 395   //
 396   const Double_t  mpshift  = -0.22278298;
 397   //
 398   // Integration normalisation
 399   //
 400   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
 401
 402   //
 403   // Utility function to use in TF1 defintition
 404   //
 405   Double_t landauGaus1(Double_t* xp, Double_t* pp)
 406   {
 407     Double_t x        = xp[0];
 408     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
 409     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
 410     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
 411     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
 412     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
 413
 414     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
 415   }
 416
 417   //
 418   // Utility function to use in TF1 defintition
 419   //
 420   Double_t landauGausN(Double_t* xp, Double_t* pp)
 421   {
 422     Double_t  x        = xp[0];
 423     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 424     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 425     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 426     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 427     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 428     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 429     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
 430
 431     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
 432                                                   n, a);
 433   }
 434   //
 435   // Utility function to use in TF1 defintition
 436   //
 437   Double_t landauGausI(Double_t* xp, Double_t* pp)
 438   {
 439     Double_t x         = xp[0];
 440     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 441     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 442     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 443     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 444     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 445     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 446
 447     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 448   }
 449
 450
 451 }
 452 //____________________________________________________________________
 453 Double_t
 454 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
 455 {
 456   //
 457   // Calculate the shifted Landau
 458   // @f[
 459   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 460   // @f]
 461   //
 462   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 463   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 464   // @f$\Delta=0,\xi=1@f$.
 465   //
 466   // Parameters:
 467   //    x      Where to evaluate @f$ f'_{L}@f$
 468   //    delta  Most probable value
 469   //    xi     The 'width' of the distribution
 470   //
 471   // Return:
 472   //    @f$ f'_{L}(x;\Delta,\xi) @f$
 473   //
 474   return TMath::Landau(x, delta - xi * mpshift, xi);
 475 }
 476 //____________________________________________________________________
 477 Double_t
 478 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
 479                            Double_t sigma, Double_t sigmaN)
 480 {
 481   //
 482   // Calculate the value of a Landau convolved with a Gaussian
 483   //
 484   // @f[
 485   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 486   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 487   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 488   // @f]
 489   //
 490   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 491   // energy loss, @f$ \xi@f$ the width of the Landau, and
 492   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 493   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 494   // noise in the detector.
 495   //
 496   // Note that this function uses the constants fgConvolutionSteps and
 497   // fgConvolutionNSigma
 498   //
 499   // References:
 500   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 501   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 502   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 503   //
 504   // Parameters:
 505   //    x         where to evaluate @f$ f@f$
 506   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 507   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 508   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 509   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 510   //
 511   // Return:
 512   //    @f$ f@f$ evaluated at @f$ x@f$.
 513   //
 514   Double_t deltap = delta - xi * mpshift;
 515   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
 516   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
 517   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
 518   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
 519   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
 520   Double_t sum    = 0;
 521
 522   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
 523     Double_t x1 = xlow  + (i - .5) * step;
 524     Double_t x2 = xhigh - (i - .5) * step;
 525
 526     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
 527     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
 528   }
 529   return step * sum * invSq2pi / sigma1;
 530 }
 531
 532 //____________________________________________________________________
 533 Double_t
 534 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 535                             Double_t sigma, Double_t sigmaN, Int_t i)
 536 {
 537   //
 538   // Evaluate
 539   // @f[
 540   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 541   // @f]
 542   // corresponding to @f$ i@f$ particles i.e., with the substitutions
 543   // @f{eqnarray*}{
 544   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
 545   //    \xi       \rightarrow \xi_i       &=& i \xi
 546   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
 547   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 548   // @f}
 549   //
 550   // Parameters:
 551   //    x        Where to evaluate
 552   //    delta    @f$ \Delta@f$
 553   //    xi       @f$ \xi@f$
 554   //    sigma    @f$ \sigma@f$
 555   //    sigma_n  @f$ \sigma_n@f$
 556   //    i        @f$ i @f$
 557   //
 558   // Return:
 559   //    @f$ f_i @f$ evaluated
 560   //
 561   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
 562   Double_t xiI    =  i * xi;
 563   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
 564   if (sigmaI < 1e-10) {
 565     // Fall back to landau
 566     return AliForwardUtil::Landau(x, deltaI, xiI);
 567   }
 568   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
 569 }
 570
 571 //____________________________________________________________________
 572 Double_t
 573 AliForwardUtil::IdLandauGausdPar(Double_t x,
 574                                  UShort_t par,   Double_t dPar,
 575                                  Double_t delta, Double_t xi,
 576                                  Double_t sigma, Double_t sigmaN,
 577                                  Int_t    i)
 578 {
 579   //
 580   // Numerically evaluate
 581   // @f[
 582   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 583   // @f]
 584   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 585   // of the parameters is given by
 586   //
 587   // - 0: @f$\Delta@f$
 588   // - 1: @f$\xi@f$
 589   // - 2: @f$\sigma@f$
 590   // - 3: @f$\sigma_n@f$
 591   //
 592   // This is the partial derivative with respect to the parameter of
 593   // the response function corresponding to @f$ i@f$ particles i.e.,
 594   // with the substitutions
 595   // @f[
 596   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
 597   //    \xi       \rightarrow \xi_i       = i \xi
 598   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
 599   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 600   // @f]
 601   //
 602   // Parameters:
 603   //    x        Where to evaluate
 604   //    ipar     Parameter number
 605   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 606   //    delta    @f$ \Delta@f$
 607   //    xi       @f$ \xi@f$
 608   //    sigma    @f$ \sigma@f$
 609   //    sigma_n  @f$ \sigma_n@f$
 610   //    i        @f$ i@f$
 611   //
 612   // Return:
 613   //    @f$ f_i@f$ evaluated
 614   //
 615   if (dPar == 0) return 0;
 616   Double_t dp      = dPar;
 617   Double_t d2      = dPar / 2;
 618   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
 619   Double_t xiI     =  i * xi;
 620   Double_t si      =  TMath::Sqrt(Double_t(i));
 621   Double_t sigmaI  =  si*sigma;
 622   Double_t y1      = 0;
 623   Double_t y2      = 0;
 624   Double_t y3      = 0;
 625   Double_t y4      = 0;
 626   switch (par) {
 627   case 0:
 628     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
 629     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
 630     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
 631     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
 632     break;
 633   case 1:
 634     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
 635     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
 636     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
 637     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
 638     break;
 639   case 2:
 640     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
 641     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
 642     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
 643     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
 644     break;
 645   case 3:
 646     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
 647     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
 648     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
 649     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
 650     break;
 651   default:
 652     return 0;
 653   }
 654
 655   Double_t d0  = y1 - y4;
 656   Double_t d1  = 2 * (y2 - y3);
 657
 658   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
 659
 660   return g;
 661 }
 662
 663 //____________________________________________________________________
 664 Double_t
 665 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 666                             Double_t sigma, Double_t sigmaN, Int_t n,
 667                             const Double_t* a)
 668 {
 669   //
 670   // Evaluate
 671   // @f[
 672   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 673   // @f]
 674   //
 675   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 676   // Landau with a Gaussian (see LandauGaus).  Note that
 677   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 678   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 679   //
 680   // References:
 681   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 682   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 683   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 684   //
 685   // Parameters:
 686   //    x        Where to evaluate @f$ f_N@f$
 687   //    delta    @f$ \Delta_1@f$
 688   //    xi       @f$ \xi_1@f$
 689   //    sigma    @f$ \sigma_1@f$
 690   //    sigma_n  @f$ \sigma_n@f$
 691   //    n        @f$ N@f$ in the sum above.
 692   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 693   //                 @f$ i > 1@f$
 694   //
 695   // Return:
 696   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$
 697   //
 698   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
 699   for (Int_t i = 2; i <= n; i++)
 700     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 701   return result;
 702 }
 703 namespace {
 704   const Int_t kColors[] = { kRed+1,
 705                             kPink+3,
 706                             kMagenta+2,
 707                             kViolet+2,
 708                             kBlue+1,
 709                             kAzure+3,
 710                             kCyan+1,
 711                             kTeal+2,
 712                             kGreen+2,
 713                             kSpring+3,
 714                             kYellow+2,
 715                             kOrange+2 };
 716 }
 717
 718 //____________________________________________________________________
 719 TF1*
 720 AliForwardUtil::MakeNLandauGaus(Double_t  c,
 721                                 Double_t  delta, Double_t xi,
 722                                 Double_t  sigma, Double_t sigmaN, Int_t n,
 723                                 const Double_t* a,
 724                                 Double_t  xmin, Double_t xmax)
 725 {
 726   //
 727   // Generate a TF1 object of @f$ f_N@f$
 728   //
 729   // Parameters:
 730   //    c         Constant
 731   //    delta     @f$ \Delta@f$
 732   //    xi            @f$ \xi_1@f$
 733   //    sigma     @f$ \sigma_1@f$
 734   //    sigma_n   @f$ \sigma_n@f$
 735   //    n             @f$ N@f$ - how many particles to sum to
 736   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 737   //                  @f$ i > 1@f$
 738   //    xmin      Least value of range
 739   //    xmax      Largest value of range
 740   //
 741   // Return:
 742   //    Newly allocated TF1 object
 743   //
 744   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
 745   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
 746   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
 747   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
 748   landaun->SetLineWidth(2);
 749   landaun->SetNpx(500);
 750   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
 751
 752   // Set the initial parameters from the seed fit
 753   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 754   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 755   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 756   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 757   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
 758   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);
 759
 760   // Set the range and name of the scale parameters
 761   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 762     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
 763     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
 764   }
 765   return landaun;
 766 }
 767 //____________________________________________________________________
 768 TF1*
 769 AliForwardUtil::MakeILandauGaus(Double_t  c,
 770                                 Double_t  delta, Double_t xi,
 771                                 Double_t  sigma, Double_t sigmaN, Int_t i,
 772                                 Double_t  xmin, Double_t xmax)
 773 {
 774   //
 775   // Generate a TF1 object of @f$ f_I@f$
 776   //
 777   // Parameters:
 778   //    c        Constant
 779   //    delta    @f$ \Delta@f$
 780   //    xi       @f$ \xi_1@f$
 781   //    sigma    @f$ \sigma_1@f$
 782   //    sigma_n  @f$ \sigma_n@f$
 783   //    i            @f$ i@f$ - the number of particles
 784   //    xmin     Least value of range
 785   //    xmax     Largest value of range
 786   //
 787   // Return:
 788   //    Newly allocated TF1 object
 789   //
 790   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
 791   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
 792   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
 793   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
 794   landaui->SetLineWidth(1);
 795   landaui->SetNpx(500);
 796   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
 797
 798   // Set the initial parameters from the seed fit
 799   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 800   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 801   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 802   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 803   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
 804   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);
 805
 806   return landaui;
 807 }
 808
 809 //====================================================================
 810 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
 811                                          Double_t maxRange,
 812                                          UShort_t minusBins)
 813   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
 814     fFitResults(0), fFunctions(0)
 815 {
 816   //
 817   // Constructor
 818   //
 819   // Parameters:
 820   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
 821   //    maxRange   Maximum range to fit to
 822   //    minusBins  The number of bins below maximum to use
 823   //
 824   fFitResults.SetOwner();
 825   fFunctions.SetOwner();
 826 }
 827 //____________________________________________________________________
 828 AliForwardUtil::ELossFitter::~ELossFitter()
 829 {
 830   //
 831   // Destructor
 832   //
 833   //
 834   fFitResults.Delete();
 835   fFunctions.Delete();
 836 }
 837 //____________________________________________________________________
 838 void
 839 AliForwardUtil::ELossFitter::Clear()
 840 {
 841   //
 842   // Clear internal arrays
 843   //
 844   //
 845   fFitResults.Clear();
 846   fFunctions.Clear();
 847 }
 848 //____________________________________________________________________
 849 TF1*
 850 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
 851 {
 852   //
 853   // Fit a 1-particle signal to the passed energy loss distribution
 854   //
 855   // Note that this function clears the internal arrays first
 856   //
 857   // Parameters:
 858   //    dist    Data to fit the function to
 859   //    sigman If larger than zero, the initial guess of the
 860   //               detector induced noise. If zero or less, then this
 861   //               parameter is ignored in the fit (fixed at 0)
 862   //
 863   // Return:
 864   //    The function fitted to the data
 865   //
 866
 867   // Clear the cache
 868   Clear();
 869
 870   // Find the fit range
 871   dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
 872
 873   // Get the bin with maximum
 874   Int_t    maxBin = dist->GetMaximumBin();
 875   Double_t maxE   = dist->GetBinLowEdge(maxBin);
 876
 877   // Get the low edge
 878   dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
 879   Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
 880   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
 881   Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
 882
 883   // Restore the range
 884   dist->GetXaxis()->SetRangeUser(0, fMaxRange);
 885
 886   // Define the function to fit
 887   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxEE,kSigmaN+1);
 888
 889   // Set initial guesses, parameter names, and limits
 890   landau1->SetParameters(1,0.5,0.07,0.1,sigman);
 891   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
 892   landau1->SetNpx(500);
 893   landau1->SetParLimits(kDelta, minE, fMaxRange);
 894   landau1->SetParLimits(kXi,    0.00, fMaxRange);
 895   landau1->SetParLimits(kSigma, 0.01, 0.1);
 896   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
 897   else              landau1->SetParLimits(kSigmaN, 0, fMaxRange);
 898
 899   // Do the fit, getting the result object
 900   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
 901   landau1->SetRange(minE, fMaxRange);
 902   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
 903   fFunctions.AddAtAndExpand(landau1, 0);
 904
 905   return landau1;
 906 }
 907 //____________________________________________________________________
 908 TF1*
 909 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
 910                                           Double_t sigman)
 911 {
 912   //
 913   // Fit a N-particle signal to the passed energy loss distribution
 914   //
 915   // If there's no 1-particle fit present, it does that first
 916   //
 917   // Parameters:
 918   //    dist   Data to fit the function to
 919   //    n      Number of particle signals to fit
 920   //    sigman If larger than zero, the initial guess of the
 921   //               detector induced noise. If zero or less, then this
 922   //               parameter is ignored in the fit (fixed at 0)
 923   //
 924   // Return:
 925   //    The function fitted to the data
 926   //
 927
 928   // Get the seed fit result
 929   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
 930   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
 931   if (!r || !f) {
 932     f = Fit1Particle(dist, sigman);
 933     r = static_cast<TFitResult*>(fFitResults.At(0));
 934     if (!r || !f) {
 935       ::Warning("FitNLandau", "No first shot at landau fit");
 936       return 0;
 937     }
 938   }
 939
 940   // Get some parameters from seed fit
 941   Double_t delta1  = r->Parameter(kDelta);
 942   Double_t xi1     = r->Parameter(kXi);
 943   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
 944   Double_t minE    = f->GetXmin();
 945
 946   // Array of weights
 947   TArrayD a(n-1);
 948   for (UShort_t i = 2; i <= n; i++)
 949     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
 950   // Make the fit function
 951   TF1* landaun     = MakeNLandauGaus(r->Parameter(kC),
 952                                      r->Parameter(kDelta),
 953                                      r->Parameter(kXi),
 954                                      r->Parameter(kSigma),
 955                                      r->Parameter(kSigmaN),
 956                                      n,a.fArray,minE,maxEi);
 957   landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
 958   landaun->SetParLimits(kXi,     0.00, fMaxRange);       // xi
 959   landaun->SetParLimits(kSigma,  0.01, 1);            // sigma
 960   // Check if we're using the noise sigma
 961   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
 962   else              landaun->SetParLimits(kSigmaN, 0, fMaxRange);
 963
 964   // Set the range and name of the scale parameters
 965   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 966     landaun->SetParLimits(kA+i-2, 0,1);
 967   }
 968
 969   // Do the fit
 970   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
 971
 972   landaun->SetRange(minE, fMaxRange);
 973   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
 974   fFunctions.AddAtAndExpand(landaun, n-1);
 975
 976   return landaun;
 977 }
 978
 979 //====================================================================
 980 AliForwardUtil::Histos::~Histos()
 981 {
 982   //
 983   // Destructor
 984   //
 985 }
 986
 987 //____________________________________________________________________
 988 void
 989 AliForwardUtil::Histos::Delete(Option_t* opt)
 990 {
 991   if (fFMD1i) delete fFMD1i;
 992   if (fFMD2i) delete fFMD2i;
 993   if (fFMD2o) delete fFMD2o;
 994   if (fFMD3i) delete fFMD3i;
 995   if (fFMD3o) delete fFMD3o;
 996   fFMD1i = 0;
 997   fFMD2i = 0;
 998   fFMD2o = 0;
 999   fFMD3i = 0;
1000   fFMD3o = 0;
1001   TObject::Delete(opt);
1002 }
1003
1004 //____________________________________________________________________
1005 TH2D*
1006 AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
1007                              const TAxis& etaAxis) const
1008 {
1009   //
1010   // Make a histogram
1011   //
1012   // Parameters:
1013   //    d        Detector
1014   //    r        Ring
1015   //    etaAxis  Eta axis to use
1016   //
1017   // Return:
1018   //    Newly allocated histogram
1019   //
1020   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
1021   TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r),
1022                         Form("FMD%d%c cache", d, r),
1023                         etaAxis.GetNbins(), etaAxis.GetXmin(),
1024                         etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
1025   hist->SetXTitle("#eta");
1026   hist->SetYTitle("#phi [radians]");
1027   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
1028   hist->Sumw2();
1029   hist->SetDirectory(0);
1030
1031   return hist;
1032 }
1033 //____________________________________________________________________
1034 void
1035 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
1036 {
1037   //
1038   // Initialize the object
1039   //
1040   // Parameters:
1041   //    etaAxis Eta axis to use
1042   //
1043   fFMD1i = Make(1, 'I', etaAxis);
1044   fFMD2i = Make(2, 'I', etaAxis);
1045   fFMD2o = Make(2, 'O', etaAxis);
1046   fFMD3i = Make(3, 'I', etaAxis);
1047   fFMD3o = Make(3, 'O', etaAxis);
1048 }
1049 //____________________________________________________________________
1050 void
1051 AliForwardUtil::Histos::Clear(Option_t* option)
1052 {
1053   //
1054   // Clear data
1055   //
1056   // Parameters:
1057   //    option Not used
1058   //
1059   if (fFMD1i) fFMD1i->Reset(option);
1060   if (fFMD2i) fFMD2i->Reset(option);
1061   if (fFMD2o) fFMD2o->Reset(option);
1062   if (fFMD3i) fFMD3i->Reset(option);
1063   if (fFMD3o) fFMD3o->Reset(option);
1064 }
1065
1066 //____________________________________________________________________
1067 TH2D*
1068 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
1069 {
1070   //
1071   // Get the histogram for a particular detector,ring
1072   //
1073   // Parameters:
1074   //    d Detector
1075   //    r Ring
1076   //
1077   // Return:
1078   //    Histogram for detector,ring or nul
1079   //
1080   switch (d) {
1081   case 1: return fFMD1i;
1082   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1083   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1084   }
1085   return 0;
1086 }
1087 //====================================================================
1088 TList*
1089 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
1090 {
1091   //
1092   // Define the outout list in @a d
1093   //
1094   // Parameters:
1095   //    d Where to put the output list
1096   //
1097   // Return:
1098   //    Newly allocated TList object or null
1099   //
1100   if (!d) return 0;
1101   TList* list = new TList;
1102   list->SetOwner();
1103   list->SetName(fName.Data());
1104   d->Add(list);
1105   return list;
1106 }
1107 //____________________________________________________________________
1108 TList*
1109 AliForwardUtil::RingHistos::GetOutputList(const TList* d) const
1110 {
1111   //
1112   // Get our output list from the container @a d
1113   //
1114   // Parameters:
1115   //    d where to get the output list from
1116   //
1117   // Return:
1118   //    The found TList or null
1119   //
1120   if (!d) return 0;
1121   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1122   return list;
1123 }
1124
1125 //____________________________________________________________________
1126 TH1*
1127 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1128 {
1129   //
1130   // Find a specific histogram in the source list @a d
1131   //
1132   // Parameters:
1133   //    d     (top)-container
1134   //    name  Name of histogram
1135   //
1136   // Return:
1137   //    Found histogram or null
1138   //
1139   return static_cast<TH1*>(d->FindObject(name));
1140 }
1141
1142 //====================================================================
1143 AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl,
1144                                        const char* format, ...)
1145   : fMsg("")
1146 {
1147   if (lvl < msgLvl) return;
1148   va_list ap;
1149   va_start(ap, format);
1150   Format(fMsg, format, ap);
1151   va_end(ap);
1152   Output(+1, fMsg);
1153 }
1154 //____________________________________________________________________
1155 AliForwardUtil::DebugGuard::~DebugGuard()
1156 {
1157   if (fMsg.IsNull()) return;
1158   Output(-1, fMsg);
1159 }
1160 //____________________________________________________________________
1161 void
1162 AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl,
1163                                     const char* format, ...)
1164 {
1165   if (lvl < msgLvl) return;
1166   TString msg;
1167   va_list ap;
1168   va_start(ap, format);
1169   Format(msg, format, ap);
1170   va_end(ap);
1171   Output(0, msg);
1172 }
1173
1174 //____________________________________________________________________
1175 void
1176 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1177 {
1178   static char buf[512];
1179   Int_t n = gROOT->GetDirLevel() + 2;
1180   for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1181   vsnprintf(&(buf[n]), 511-n, format, ap);
1182   buf[511] = '\0';
1183   out = buf;
1184 }
1185 //____________________________________________________________________
1186 void
1187 AliForwardUtil::DebugGuard::Output(int in, TString& msg)
1188 {
1189   msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
1190   AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1191   if      (in > 0) gROOT->IncreaseDirLevel();
1192   else if (in < 0) gROOT->DecreaseDirLevel();
1193 }
1194
1195
1196
1197 //
1198 // EOF
1199 //