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