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