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