PWGLF/FORWARD/analysis2/AliForwardUtil.cxx

   1 //
   2 // Utilities used in the forward multiplcity analysis
   3 //
   4 //
   5 #include "AliForwardUtil.h"
   6 //#include <ARVersion.h>
   7 #include <AliAnalysisManager.h>
   8 #include "AliAODForwardMult.h"
   9 #include <AliLog.h>
  10 #include <AliInputEventHandler.h>
  11 #include <AliAODInputHandler.h>
  12 #include <AliAODHandler.h>
  13 #include <AliAODEvent.h>
  14 #include <AliESDEvent.h>
  15 #include <AliAnalysisTaskSE.h>
  16 #include <AliPhysicsSelection.h>
  17 #include <AliTriggerAnalysis.h>
  18 #include <AliMultiplicity.h>
  19 #include <TParameter.h>
  20 #include <TH2D.h>
  21 #include <TH1I.h>
  22 #include <TF1.h>
  23 #include <TFitResult.h>
  24 #include <TMath.h>
  25 #include <TError.h>
  26 #include <TROOT.h>
  27
  28 //====================================================================
  29 ULong_t AliForwardUtil::AliROOTRevision()
  30 {
  31 #ifdef ALIROOT_SVN_REVISION
  32   return ALIROOT_SVN_REVISION;
  33 #else
  34   return 0;
  35 #endif
  36 }
  37 //____________________________________________________________________
  38 ULong_t AliForwardUtil::AliROOTBranch()
  39 {
  40 #ifdef ALIROOT_SVN_BRANCH
  41   static ULong_t ret = 0;
  42   if (ret != 0) return ret;
  43
  44   TString str(ALIROOT_SVN_BRANCH);
  45   if (str[0] == 'v') str.Remove(0,1);
  46   if (str.EqualTo("trunk")) return ret = 0xFFFFFFFF;
  47
  48   TObjArray*   tokens = str.Tokenize("-");
  49   TObjString*  pMajor = static_cast<TObjString*>(tokens->At(0));
  50   TObjString*  pMinor = static_cast<TObjString*>(tokens->At(1));
  51   TObjString*  pRelea = static_cast<TObjString*>(tokens->At(2));
  52   TObjString* pAn     = (tokens->GetEntries() > 3 ?
  53     static_cast<TObjString*>(tokens->At(3)) : 0);
  54   TString sMajor = pMajor->String().Strip(TString::kLeading, '0');
  55   TString sMinor = pMinor->String().Strip(TString::kLeading, '0');
  56   TString sRelea = pRelea->String().Strip(TString::kLeading, '0');
  57
  58   ret = (((sMajor.Atoi() & 0xFF) << 12) |
  59     ((sMinor.Atoi() & 0xFF) <<  8) |
  60     ((sRelea.Atoi() & 0xFF) <<  4) |
  61     (pAn ? 0xAA : 0));
  62
  63   return ret;
  64 #else
  65   return 0;
  66 #endif
  67 }
  68
  69 //====================================================================
  70 UShort_t
  71 AliForwardUtil::ParseCollisionSystem(const char* sys)
  72 {
  73   //
  74   // Parse a collision system spec given in a string.   Known values are
  75   //
  76   //  - "ppb", "p-pb", "pa", "p-a"  which returns kPPb
  77   //  - "pp", "p-p"                 which returns kPP
  78   //  - "PbPb", "Pb-Pb", "A-A",     which returns kPbPb
  79   //  - Everything else gives kUnknown
  80   //
  81   // Parameters:
  82   //    sys Collision system spec
  83   //
  84   // Return:
  85   //    Collision system id
  86   //
  87   TString s(sys);
  88   s.ToLower();
  89   // we do pA first to avoid pp catch on ppb string (AH)
  90   if (s.Contains("p-pb")  || s.Contains("ppb"))   return AliForwardUtil::kPPb;
  91   if (s.Contains("p-a")   || s.Contains("pa"))    return AliForwardUtil::kPPb;
  92   if (s.Contains("a-p")   || s.Contains("ap"))    return AliForwardUtil::kPPb;
  93   if (s.Contains("p-p")   || s.Contains("pp"))    return AliForwardUtil::kPP;
  94   if (s.Contains("pb-pb") || s.Contains("pbpb"))  return AliForwardUtil::kPbPb;
  95   if (s.Contains("a-a")   || s.Contains("aa"))    return AliForwardUtil::kPbPb;
  96   return AliForwardUtil::kUnknown;
  97 }
  98 //____________________________________________________________________
  99 const char*
 100 AliForwardUtil::CollisionSystemString(UShort_t sys)
 101 {
 102   //
 103   // Get a string representation of the collision system
 104   //
 105   // Parameters:
 106   //    sys  Collision system
 107   // - kPP -> "pp"
 108   // - kPbPb -> "PbPb"
 109   // - anything else gives "unknown"
 110   //
 111   // Return:
 112   //    String representation of the collision system
 113   //
 114   switch (sys) {
 115   case AliForwardUtil::kPP:   return "pp";
 116   case AliForwardUtil::kPbPb: return "PbPb";
 117   case AliForwardUtil::kPPb:  return "pPb";
 118   }
 119   return "unknown";
 120 }
 121 //____________________________________________________________________
 122 Float_t
 123 AliForwardUtil::BeamRapidity(Float_t beam, UShort_t z, UShort_t a)
 124 {
 125   const Double_t pMass = 9.38271999999999995e-01;
 126   const Double_t nMass = 9.39564999999999984e-01;
 127   Double_t       beamE = z * beam / 2;
 128   Double_t       beamM = z * pMass + (a - z) * nMass;
 129   Double_t       beamP = TMath::Sqrt(beamE * beamE - beamM * beamM);
 130   Double_t       beamY = .5* TMath::Log((beamE+beamP) / (beamE-beamP));
 131   return beamY;
 132 }
 133 //____________________________________________________________________
 134 Float_t
 135 AliForwardUtil::CenterOfMassEnergy(Float_t beam,
 136                                    UShort_t z1,
 137                                    UShort_t a1,
 138                                    Short_t z2,
 139                                    Short_t a2)
 140 {
 141   // Calculate the center of mass energy given target/projectile
 142   // mass and charge numbers
 143   if (z2 < 0) z2 = z1;
 144   if (a2 < 0) a2 = a1;
 145   return TMath::Sqrt(Float_t(z1*z2)/a1/a2) * beam;
 146 }
 147 //____________________________________________________________________
 148 Float_t
 149 AliForwardUtil::CenterOfMassRapidity(UShort_t z1,
 150                                      UShort_t a1,
 151                                      Short_t z2,
 152                                      Short_t a2)
 153 {
 154   // Calculate the center of mass rapidity (shift) given target/projectile
 155   // mass and charge numbers
 156   if (z2 < 0) z2 = z1;
 157   if (a2 < 0) a2 = a1;
 158   if (z2 == z1 && a2 == a1) return 0;
 159   return .5 * TMath::Log(Float_t(z1*a2)/z2/a1);
 160 }
 161
 162 //____________________________________________________________________
 163 UShort_t
 164 AliForwardUtil::ParseCenterOfMassEnergy(UShort_t sys, Float_t beam)
 165 {
 166   //
 167   // Parse the center of mass energy given as a float and return known
 168   // values as a unsigned integer
 169   //
 170   // Parameters:
 171   //    sys   Collision system (needed for AA)
 172   //    beam  Center of mass energy * total charge
 173   //
 174   // Return:
 175   //    Center of mass energy per nucleon
 176   //
 177   Float_t energy = beam;
 178   // Below no longer needed apparently
 179   // if (sys == AliForwardUtil::kPbPb) energy = energy / 208 * 82;
 180   if (sys == AliForwardUtil::kPPb)
 181     energy = CenterOfMassEnergy(beam, 82, 208, 1, 1);
 182   else if (sys == AliForwardUtil::kPbPb)
 183     energy = CenterOfMassEnergy(beam, 82, 208, 82, 208);
 184   if (TMath::Abs(energy - 900.)   < 10)  return 900;
 185   if (TMath::Abs(energy - 2400.)  < 10)  return 2400;
 186   if (TMath::Abs(energy - 2760.)  < 20)  return 2760;
 187   if (TMath::Abs(energy - 4400.)  < 10)  return 4400;
 188   if (TMath::Abs(energy - 5022.)  < 10)  return 5023;
 189   if (TMath::Abs(energy - 5500.)  < 40)  return 5500;
 190   if (TMath::Abs(energy - 7000.)  < 10)  return 7000;
 191   if (TMath::Abs(energy - 8000.)  < 10)  return 8000;
 192   if (TMath::Abs(energy - 10000.) < 10)  return 10000;
 193   if (TMath::Abs(energy - 14000.) < 10)  return 14000;
 194   return 0;
 195 }
 196 //____________________________________________________________________
 197 const char*
 198 AliForwardUtil::CenterOfMassEnergyString(UShort_t cms)
 199 {
 200   //
 201   // Get a string representation of the center of mass energy per nuclean
 202   //
 203   // Parameters:
 204   //    cms  Center of mass energy per nucleon
 205   //
 206   // Return:
 207   //    String representation of the center of mass energy per nuclean
 208   //
 209   return Form("%04dGeV", cms);
 210 }
 211 //____________________________________________________________________
 212 Short_t
 213 AliForwardUtil::ParseMagneticField(Float_t v)
 214 {
 215   //
 216   // Parse the magnetic field (in kG) as given by a floating point number
 217   //
 218   // Parameters:
 219   //    field  Magnetic field in kG
 220   //
 221   // Return:
 222   //    Short integer value of magnetic field in kG
 223   //
 224   if (TMath::Abs(v - 5.) < 1 ) return +5;
 225   if (TMath::Abs(v + 5.) < 1 ) return -5;
 226   if (TMath::Abs(v) < 1)       return 0;
 227   return 999;
 228 }
 229 //____________________________________________________________________
 230 const char*
 231 AliForwardUtil::MagneticFieldString(Short_t f)
 232 {
 233   //
 234   // Get a string representation of the magnetic field
 235   //
 236   // Parameters:
 237   //    field Magnetic field in kG
 238   //
 239   // Return:
 240   //    String representation of the magnetic field
 241   //
 242   return Form("%01dkG", f);
 243 }
 244 //_____________________________________________________________________
 245 AliAODEvent* AliForwardUtil::GetAODEvent(AliAnalysisTaskSE* task)
 246 {
 247   // Check if AOD is the output event
 248   if (!task) ::Fatal("GetAODEvent", "Null task given, cannot do that");
 249
 250   AliAODEvent* ret = task->AODEvent();
 251   if (ret) return ret;
 252
 253   // Check if AOD is the input event
 254   ret = dynamic_cast<AliAODEvent*>(task->InputEvent());
 255   if (!ret) ::Warning("GetAODEvent", "No AOD event found");
 256
 257   return ret;
 258 }
 259 //_____________________________________________________________________
 260 UShort_t AliForwardUtil::CheckForAOD()
 261 {
 262   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 263   if (dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler())) {
 264     ::Info("CheckForAOD", "Found AOD Input handler");
 265     return 1;
 266   }
 267   if (dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler())) {
 268     ::Info("CheckForAOD", "Found AOD Output handler");
 269     return 2;
 270   }
 271
 272   ::Warning("CheckForAOD",
 273             "Neither and input nor output AOD handler is specified");
 274   return 0;
 275 }
 276 //_____________________________________________________________________
 277 Bool_t AliForwardUtil::CheckForTask(const char* clsOrName, Bool_t cls)
 278 {
 279   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 280   if (!cls) {
 281     AliAnalysisTask* t = am->GetTask(clsOrName);
 282     if (!t) {
 283       ::Warning("CheckForTask", "Task %s not found in manager", clsOrName);
 284       return false;
 285     }
 286     ::Info("CheckForTask", "Found task %s", clsOrName);
 287     return true;
 288   }
 289   TClass* dep = gROOT->GetClass(clsOrName);
 290   if (!dep) {
 291     ::Warning("CheckForTask", "Unknown class %s for needed task", clsOrName);
 292     return false;
 293   }
 294   TIter next(am->GetTasks());
 295   TObject* o = 0;
 296   while ((o = next())) {
 297     if (o->IsA()->InheritsFrom(dep)) {
 298       ::Info("CheckForTask", "Found task of class %s: %s",
 299              clsOrName, o->GetName());
 300       return true;
 301     }
 302   }
 303   ::Warning("CheckForTask", "No task of class %s was found", clsOrName);
 304   return false;
 305 }
 306
 307 //_____________________________________________________________________
 308 TObject* AliForwardUtil::MakeParameter(const Char_t* name, UShort_t value)
 309 {
 310   TParameter<int>* ret = new TParameter<int>(name, value);
 311   ret->SetMergeMode('f');
 312   ret->SetUniqueID(value);
 313   return ret;
 314 }
 315 //_____________________________________________________________________
 316 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Int_t value)
 317 {
 318   TParameter<int>* ret = new TParameter<int>(name, value);
 319   ret->SetMergeMode('f');
 320   ret->SetUniqueID(value);
 321   return ret;
 322 }
 323 //_____________________________________________________________________
 324 TObject* AliForwardUtil::MakeParameter(const Char_t* name, ULong_t value)
 325 {
 326   TParameter<Long_t>* ret = new TParameter<Long_t>(name, value);
 327   ret->SetMergeMode('f');
 328   ret->SetUniqueID(value);
 329   return ret;
 330 }
 331 //_____________________________________________________________________
 332 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Double_t value)
 333 {
 334   TParameter<double>* ret = new TParameter<double>(name, value);
 335   // Float_t v = value;
 336   // UInt_t* tmp = reinterpret_cast<UInt_t*>(&v);
 337   ret->SetMergeMode('f');
 338   // ret->SetUniqueID(*tmp);
 339   return ret;
 340 }
 341 //_____________________________________________________________________
 342 TObject* AliForwardUtil::MakeParameter(const Char_t* name, Bool_t value)
 343 {
 344   TParameter<bool>* ret = new TParameter<bool>(name, value);
 345   ret->SetMergeMode('f');
 346   ret->SetUniqueID(value);
 347   return ret;
 348 }
 349
 350 //_____________________________________________________________________
 351 void AliForwardUtil::GetParameter(TObject* o, UShort_t& value)
 352 {
 353   if (!o) return;
 354   TParameter<int>* p = static_cast<TParameter<int>*>(o);
 355   if (p->TestBit(BIT(17)))
 356     value = p->GetVal();
 357   else
 358     value = o->GetUniqueID();
 359 }
 360 //_____________________________________________________________________
 361 void AliForwardUtil::GetParameter(TObject* o, Int_t& value)
 362 {
 363   if (!o) return;
 364   TParameter<int>* p = static_cast<TParameter<int>*>(o);
 365   if (p->TestBit(BIT(17)))
 366     value = p->GetVal();
 367   else
 368     value = o->GetUniqueID();
 369 }
 370 //_____________________________________________________________________
 371 void AliForwardUtil::GetParameter(TObject* o, ULong_t& value)
 372 {
 373   if (!o) return;
 374   TParameter<Long_t>* p = static_cast<TParameter<Long_t>*>(o);
 375   if (p->TestBit(BIT(17)))
 376     value = p->GetVal();
 377   else
 378     value = o->GetUniqueID();
 379 }
 380 //_____________________________________________________________________
 381 void AliForwardUtil::GetParameter(TObject* o, Double_t& value)
 382 {
 383   if (!o) return;
 384   TParameter<double>* p = static_cast<TParameter<double>*>(o);
 385   if (p->TestBit(BIT(17)))
 386     value = p->GetVal(); // o->GetUniqueID();
 387   else {
 388     UInt_t  i = o->GetUniqueID();
 389     Float_t v = *reinterpret_cast<Float_t*>(&i);
 390     value = v;
 391   }
 392 }
 393 //_____________________________________________________________________
 394 void AliForwardUtil::GetParameter(TObject* o, Bool_t& value)
 395 {
 396   if (!o) return;
 397   TParameter<bool>* p = static_cast<TParameter<bool>*>(o);
 398   if (p->TestBit(BIT(17)))
 399     value = p->GetVal(); // o->GetUniqueID();
 400   else
 401     value = o->GetUniqueID();
 402 }
 403
 404 #if 0
 405 //_____________________________________________________________________
 406 Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
 407 {
 408   // Get max R of ring
 409   //
 410   // Optimized version that has a cache
 411   static TArrayD inner;
 412   static TArrayD outer;
 413   if (inner.GetSize() <= 0 || outer.GetSize() <= 0) {
 414     const Double_t minR[] = {  4.5213, 15.4 };
 415     const Double_t maxR[] = { 17.2,    28.0 };
 416     const Int_t    nStr[] = { 512,     256  };
 417     for (Int_t q = 0; q < 2; q++) {
 418       TArrayD& a = (q == 0 ? inner : outer);
 419       a.Set(nStr[q]);
 420
 421       for (Int_t it = 0; it < nStr[q]; it++) {
 422         Double_t   rad     = maxR[q] - minR[q];
 423         Double_t   segment = rad / nStr[q];
 424         Double_t   r       = minR[q] + segment*strip;
 425         a[it]              = r;
 426       }
 427     }
 428   }
 429   if (ring == 'I' || ring == 'i') return inner.At(strip);
 430   return outer.At(strip);
 431 }
 432 #else
 433 //_____________________________________________________________________
 434 Double_t AliForwardUtil::GetStripR(Char_t ring, UShort_t strip)
 435 {
 436   // Get max R of ring
 437   //
 438   // New implementation has only one branch
 439   const Double_t minR[] = {  4.5213, 15.4 };
 440   const Double_t maxR[] = { 17.2,    28.0 };
 441   const Int_t    nStr[] = { 512,     256  };
 442
 443   Int_t      q       = (ring == 'I' || ring == 'i') ? 0 : 1;
 444   Double_t   rad     = maxR[q] - minR[q];
 445   Double_t   segment = rad / nStr[q];
 446   Double_t   r       = minR[q] + segment*strip;
 447
 448   return r;
 449 }
 450 #endif
 451
 452 #if 1
 453 //_____________________________________________________________________
 454 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring,
 455                                          UShort_t sec, UShort_t strip,
 456                                          Double_t zvtx)
 457 {
 458   // Calculate eta from strip with vertex (redundant with
 459   // AliESDFMD::Eta but support displaced vertices)
 460   //
 461   // Slightly more optimized version that uses less branching
 462
 463   // Get R of the strip
 464   Double_t   r         = GetStripR(ring, strip);
 465   Int_t      hybrid    = sec / 2;
 466   Int_t      q        = (ring == 'I' || ring == 'i') ? 0 : 1;
 467
 468   const Double_t zs[][2] = { { 320.266, -999999 },
 469                             {  83.666,  74.966 },
 470                             { -63.066, -74.966 } };
 471   if (det > 3 || zs[det-1][q] == -999999) return -999999;
 472
 473   Double_t z = zs[det-1][q];
 474   if ((hybrid % 2) == 0) z -= .5;
 475
 476   Double_t   theta = TMath::ATan2(r,z-zvtx);
 477   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 478
 479   return eta;
 480 }
 481 #else
 482 //_____________________________________________________________________
 483 Double_t AliForwardUtil::GetEtaFromStrip(UShort_t det, Char_t ring,
 484                                          UShort_t sec, UShort_t strip,
 485                                          Double_t zvtx)
 486 {
 487   // Calculate eta from strip with vertex (redundant with
 488   // AliESDFMD::Eta but support displaced vertices)
 489
 490   //Get max R of ring
 491   Double_t   r         = GetStripR(ring, strip);
 492   Int_t      hybrid    = sec / 2;
 493   Bool_t     inner     = (ring == 'I' || ring == 'i');
 494   Double_t   z         = 0;
 495
 496
 497   switch (det) {
 498   case 1: z = 320.266;                     break;
 499   case 2: z = (inner ?  83.666 :  74.966); break;
 500   case 3: z = (inner ? -63.066 : -74.966); break;
 501   default: return -999999;
 502   }
 503   if ((hybrid % 2) == 0) z -= .5;
 504
 505   Double_t   theta = TMath::ATan2(r,z-zvtx);
 506   Double_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 507
 508   return eta;
 509 }
 510 #endif
 511
 512 //_____________________________________________________________________
 513 Double_t AliForwardUtil::GetPhiFromStrip(Char_t ring, UShort_t strip,
 514                                          Double_t phi,
 515                                          Double_t xvtx, Double_t yvtx)
 516 {
 517   // Calculate eta from strip with vertex (redundant with
 518   // AliESDFMD::Eta but support displaced vertices)
 519
 520   // Unknown x,y -> no change
 521   if (yvtx > 999 || xvtx > 999) return phi;
 522
 523   //Get max R of ring
 524   Double_t r   = GetStripR(ring, strip);
 525   Double_t amp = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx) / r;
 526   Double_t pha = (TMath::Abs(yvtx) < 1e-12  ? 0 : TMath::ATan2(xvtx, yvtx));
 527   Double_t cha = amp * TMath::Cos(phi+pha);
 528   phi += cha;
 529   if (phi < 0)              phi += TMath::TwoPi();
 530   if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
 531   return phi;
 532 }
 533
 534 //====================================================================
 535 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
 536 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
 537 namespace {
 538   //
 539   // The shift of the most probable value for the ROOT function TMath::Landau
 540   //
 541   const Double_t  mpshift  = -0.22278298;
 542   //
 543   // Integration normalisation
 544   //
 545   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
 546
 547   //
 548   // Utility function to use in TF1 defintition
 549   //
 550   Double_t landauGaus1(Double_t* xp, Double_t* pp)
 551   {
 552     Double_t x        = xp[0];
 553     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
 554     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
 555     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
 556     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
 557     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
 558
 559     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
 560   }
 561
 562   Double_t landauGausComposite(Double_t* xp, Double_t* pp)
 563   {
 564     Double_t x           = xp[0];
 565     Double_t cP          = pp[AliForwardUtil::ELossFitter::kC];
 566     Double_t deltaP      = pp[AliForwardUtil::ELossFitter::kDelta];
 567     Double_t xiP         = pp[AliForwardUtil::ELossFitter::kXi];
 568     Double_t sigmaP      = pp[AliForwardUtil::ELossFitter::kSigma];
 569     Double_t cS          = pp[AliForwardUtil::ELossFitter::kSigma+1];
 570     Double_t deltaS      = pp[AliForwardUtil::ELossFitter::kSigma+2];
 571     Double_t xiS         = pp[AliForwardUtil::ELossFitter::kSigma+3];
 572     Double_t sigmaS      = pp[AliForwardUtil::ELossFitter::kSigma+4];
 573
 574     return (cP * AliForwardUtil::LandauGaus(x,deltaP,xiP,sigmaP,0) +
 575             cS * AliForwardUtil::LandauGaus(x,deltaS,xiS,sigmaS,0));
 576   }
 577
 578   //
 579   // Utility function to use in TF1 defintition
 580   //
 581   Double_t landauGausN(Double_t* xp, Double_t* pp)
 582   {
 583     Double_t  x        = xp[0];
 584     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 585     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 586     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 587     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 588     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 589     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 590     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
 591
 592     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
 593                                                   n, a);
 594   }
 595   //
 596   // Utility function to use in TF1 defintition
 597   //
 598   Double_t landauGausI(Double_t* xp, Double_t* pp)
 599   {
 600     Double_t x         = xp[0];
 601     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 602     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 603     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 604     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 605     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 606     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 607
 608     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 609   }
 610
 611
 612 }
 613 //____________________________________________________________________
 614 Double_t
 615 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
 616 {
 617   //
 618   // Calculate the shifted Landau
 619   // @f[
 620   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 621   // @f]
 622   //
 623   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 624   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 625   // @f$\Delta=0,\xi=1@f$.
 626   //
 627   // Parameters:
 628   //    x      Where to evaluate @f$ f'_{L}@f$
 629   //    delta  Most probable value
 630   //    xi     The 'width' of the distribution
 631   //
 632   // Return:
 633   //    @f$ f'_{L}(x;\Delta,\xi) @f$
 634   //
 635   return TMath::Landau(x, delta - xi * mpshift, xi);
 636 }
 637 //____________________________________________________________________
 638 Double_t
 639 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
 640                            Double_t sigma, Double_t sigmaN)
 641 {
 642   //
 643   // Calculate the value of a Landau convolved with a Gaussian
 644   //
 645   // @f[
 646   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 647   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 648   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 649   // @f]
 650   //
 651   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 652   // energy loss, @f$ \xi@f$ the width of the Landau, and
 653   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 654   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 655   // noise in the detector.
 656   //
 657   // Note that this function uses the constants fgConvolutionSteps and
 658   // fgConvolutionNSigma
 659   //
 660   // References:
 661   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 662   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 663   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 664   //
 665   // Parameters:
 666   //    x         where to evaluate @f$ f@f$
 667   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 668   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 669   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 670   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 671   //
 672   // Return:
 673   //    @f$ f@f$ evaluated at @f$ x@f$.
 674   //
 675   Double_t deltap = delta - xi * mpshift;
 676   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
 677   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
 678   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
 679   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
 680   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
 681   Double_t sum    = 0;
 682
 683   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
 684     Double_t x1 = xlow  + (i - .5) * step;
 685     Double_t x2 = xhigh - (i - .5) * step;
 686
 687     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
 688     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
 689   }
 690   return step * sum * invSq2pi / sigma1;
 691 }
 692
 693 //____________________________________________________________________
 694 Double_t
 695 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 696                             Double_t sigma, Double_t sigmaN, Int_t i)
 697 {
 698   //
 699   // Evaluate
 700   // @f[
 701   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 702   // @f]
 703   // corresponding to @f$ i@f$ particles i.e., with the substitutions
 704   // @f{eqnarray*}{
 705   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
 706   //    \xi       \rightarrow \xi_i       &=& i \xi
 707   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
 708   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 709   // @f}
 710   //
 711   // Parameters:
 712   //    x        Where to evaluate
 713   //    delta    @f$ \Delta@f$
 714   //    xi       @f$ \xi@f$
 715   //    sigma    @f$ \sigma@f$
 716   //    sigma_n  @f$ \sigma_n@f$
 717   //    i        @f$ i @f$
 718   //
 719   // Return:
 720   //    @f$ f_i @f$ evaluated
 721   //
 722   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
 723   Double_t xiI    =  i * xi;
 724   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
 725   if (sigmaI < 1e-10) {
 726     // Fall back to landau
 727     return AliForwardUtil::Landau(x, deltaI, xiI);
 728   }
 729   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
 730 }
 731
 732 //____________________________________________________________________
 733 Double_t
 734 AliForwardUtil::IdLandauGausdPar(Double_t x,
 735                                  UShort_t par,   Double_t dPar,
 736                                  Double_t delta, Double_t xi,
 737                                  Double_t sigma, Double_t sigmaN,
 738                                  Int_t    i)
 739 {
 740   //
 741   // Numerically evaluate
 742   // @f[
 743   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 744   // @f]
 745   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 746   // of the parameters is given by
 747   //
 748   // - 0: @f$\Delta@f$
 749   // - 1: @f$\xi@f$
 750   // - 2: @f$\sigma@f$
 751   // - 3: @f$\sigma_n@f$
 752   //
 753   // This is the partial derivative with respect to the parameter of
 754   // the response function corresponding to @f$ i@f$ particles i.e.,
 755   // with the substitutions
 756   // @f[
 757   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
 758   //    \xi       \rightarrow \xi_i       = i \xi
 759   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
 760   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 761   // @f]
 762   //
 763   // Parameters:
 764   //    x        Where to evaluate
 765   //    ipar     Parameter number
 766   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 767   //    delta    @f$ \Delta@f$
 768   //    xi       @f$ \xi@f$
 769   //    sigma    @f$ \sigma@f$
 770   //    sigma_n  @f$ \sigma_n@f$
 771   //    i        @f$ i@f$
 772   //
 773   // Return:
 774   //    @f$ f_i@f$ evaluated
 775   //
 776   if (dPar == 0) return 0;
 777   Double_t dp      = dPar;
 778   Double_t d2      = dPar / 2;
 779   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
 780   Double_t xiI     =  i * xi;
 781   Double_t si      =  TMath::Sqrt(Double_t(i));
 782   Double_t sigmaI  =  si*sigma;
 783   Double_t y1      = 0;
 784   Double_t y2      = 0;
 785   Double_t y3      = 0;
 786   Double_t y4      = 0;
 787   switch (par) {
 788   case 0:
 789     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
 790     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
 791     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
 792     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
 793     break;
 794   case 1:
 795     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
 796     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
 797     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
 798     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
 799     break;
 800   case 2:
 801     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
 802     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
 803     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
 804     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
 805     break;
 806   case 3:
 807     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
 808     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
 809     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
 810     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
 811     break;
 812   default:
 813     return 0;
 814   }
 815
 816   Double_t d0  = y1 - y4;
 817   Double_t d1  = 2 * (y2 - y3);
 818
 819   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
 820
 821   return g;
 822 }
 823
 824 //____________________________________________________________________
 825 Double_t
 826 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 827                             Double_t sigma, Double_t sigmaN, Int_t n,
 828                             const Double_t* a)
 829 {
 830   //
 831   // Evaluate
 832   // @f[
 833   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 834   // @f]
 835   //
 836   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 837   // Landau with a Gaussian (see LandauGaus).  Note that
 838   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 839   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 840   //
 841   // References:
 842   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 843   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 844   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 845   //
 846   // Parameters:
 847   //    x        Where to evaluate @f$ f_N@f$
 848   //    delta    @f$ \Delta_1@f$
 849   //    xi       @f$ \xi_1@f$
 850   //    sigma    @f$ \sigma_1@f$
 851   //    sigma_n  @f$ \sigma_n@f$
 852   //    n        @f$ N@f$ in the sum above.
 853   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 854   //                 @f$ i > 1@f$
 855   //
 856   // Return:
 857   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$
 858   //
 859   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
 860   for (Int_t i = 2; i <= n; i++)
 861     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 862   return result;
 863 }
 864 namespace {
 865   const Int_t kColors[] = { kRed+1,
 866                             kPink+3,
 867                             kMagenta+2,
 868                             kViolet+2,
 869                             kBlue+1,
 870                             kAzure+3,
 871                             kCyan+1,
 872                             kTeal+2,
 873                             kGreen+2,
 874                             kSpring+3,
 875                             kYellow+2,
 876                             kOrange+2 };
 877 }
 878
 879 //____________________________________________________________________
 880 TF1*
 881 AliForwardUtil::MakeNLandauGaus(Double_t  c,
 882                                 Double_t  delta, Double_t xi,
 883                                 Double_t  sigma, Double_t sigmaN, Int_t n,
 884                                 const Double_t* a,
 885                                 Double_t  xmin, Double_t xmax)
 886 {
 887   //
 888   // Generate a TF1 object of @f$ f_N@f$
 889   //
 890   // Parameters:
 891   //    c         Constant
 892   //    delta     @f$ \Delta@f$
 893   //    xi            @f$ \xi_1@f$
 894   //    sigma     @f$ \sigma_1@f$
 895   //    sigma_n   @f$ \sigma_n@f$
 896   //    n             @f$ N@f$ - how many particles to sum to
 897   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
 898   //                  @f$ i > 1@f$
 899   //    xmin      Least value of range
 900   //    xmax      Largest value of range
 901   //
 902   // Return:
 903   //    Newly allocated TF1 object
 904   //
 905   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
 906   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
 907   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
 908   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
 909   landaun->SetLineWidth(2);
 910   landaun->SetNpx(500);
 911   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
 912
 913   // Set the initial parameters from the seed fit
 914   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 915   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 916   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 917   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 918   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
 919   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);
 920
 921   // Set the range and name of the scale parameters
 922   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
 923     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
 924     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
 925   }
 926   return landaun;
 927 }
 928 //____________________________________________________________________
 929 TF1*
 930 AliForwardUtil::MakeILandauGaus(Double_t  c,
 931                                 Double_t  delta, Double_t xi,
 932                                 Double_t  sigma, Double_t sigmaN, Int_t i,
 933                                 Double_t  xmin, Double_t xmax)
 934 {
 935   //
 936   // Generate a TF1 object of @f$ f_I@f$
 937   //
 938   // Parameters:
 939   //    c        Constant
 940   //    delta    @f$ \Delta@f$
 941   //    xi       @f$ \xi_1@f$
 942   //    sigma    @f$ \sigma_1@f$
 943   //    sigma_n  @f$ \sigma_n@f$
 944   //    i            @f$ i@f$ - the number of particles
 945   //    xmin     Least value of range
 946   //    xmax     Largest value of range
 947   //
 948   // Return:
 949   //    Newly allocated TF1 object
 950   //
 951   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
 952   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
 953   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
 954   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
 955   landaui->SetLineWidth(1);
 956   landaui->SetNpx(500);
 957   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
 958
 959   // Set the initial parameters from the seed fit
 960   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
 961   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
 962   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
 963   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
 964   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
 965   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);
 966
 967   return landaui;
 968 }
 969
 970 //====================================================================
 971 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
 972                                          Double_t maxRange,
 973                                          UShort_t minusBins)
 974   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
 975     fFitResults(0), fFunctions(0), fDebug(false)
 976 {
 977   //
 978   // Constructor
 979   //
 980   // Parameters:
 981   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
 982   //    maxRange   Maximum range to fit to
 983   //    minusBins  The number of bins below maximum to use
 984   //
 985   fFitResults.SetOwner();
 986   fFunctions.SetOwner();
 987 }
 988 //____________________________________________________________________
 989 AliForwardUtil::ELossFitter::~ELossFitter()
 990 {
 991   //
 992   // Destructor
 993   //
 994   //
 995   fFitResults.Delete();
 996   fFunctions.Delete();
 997 }
 998 //____________________________________________________________________
 999 void
1000 AliForwardUtil::ELossFitter::Clear()
1001 {
1002   //
1003   // Clear internal arrays
1004   //
1005   //
1006   fFitResults.Clear();
1007   fFunctions.Clear();
1008 }
1009 //____________________________________________________________________
1010 TF1*
1011 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
1012 {
1013   //
1014   // Fit a 1-particle signal to the passed energy loss distribution
1015   //
1016   // Note that this function clears the internal arrays first
1017   //
1018   // Parameters:
1019   //    dist    Data to fit the function to
1020   //    sigman If larger than zero, the initial guess of the
1021   //               detector induced noise. If zero or less, then this
1022   //               parameter is ignored in the fit (fixed at 0)
1023   //
1024   // Return:
1025   //    The function fitted to the data
1026   //
1027
1028   // Clear the cache
1029   Clear();
1030
1031   // Find the fit range
1032   // Find the fit range
1033   Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
1034   Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
1035                                 dist->GetNbinsX());
1036   dist->GetXaxis()->SetRange(cutBin, maxBin);
1037   // dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
1038
1039   // Get the bin with maximum
1040   Int_t    peakBin = dist->GetMaximumBin();
1041   Double_t peakE   = dist->GetBinLowEdge(peakBin);
1042
1043   // Get the low edge
1044   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
1045   Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
1046   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
1047   Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
1048
1049   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1050   Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
1051   Double_t intg  = dist->Integral(minEb, maxEb);
1052   if (intg <= 0) {
1053     ::Warning("Fit1Particle",
1054               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
1055               dist->GetName(), minE, maxE, minEb, maxEb, intg);
1056     return 0;
1057   }
1058
1059   // Restore the range
1060   dist->GetXaxis()->SetRange(1, maxBin);
1061
1062   // Define the function to fit
1063   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxE,kSigmaN+1);
1064
1065   // Set initial guesses, parameter names, and limits
1066   landau1->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
1067   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
1068   landau1->SetNpx(500);
1069   if (peakE >= minE && peakE <= fMaxRange) {
1070     // printf("Fit1: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
1071     landau1->SetParLimits(kDelta, minE, fMaxRange);
1072   }
1073   if (peakE/10 >= 0 && peakE <= 0.1) {
1074     // printf("Fit1: Set par limits on xi: %f, %f\n", 0., 0.1);
1075     landau1->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
1076   }
1077   if (peakE/5 >= 0 && peakE/5 <= 0.1) {
1078     // printf("Fit1: Set par limits on sigma: %f, %f\n", 0., 0.1);
1079     landau1->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
1080   }
1081   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
1082   else {
1083     // printf("Fit1: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
1084     landau1->SetParLimits(kSigmaN, 0, fMaxRange);
1085   }
1086
1087   // Do the fit, getting the result object
1088   if (fDebug)
1089     ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxE);
1090   TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxE);
1091   // landau1->SetRange(minE, fMaxRange);
1092   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
1093   fFunctions.AddAtAndExpand(landau1, 0);
1094
1095   return landau1;
1096 }
1097 //____________________________________________________________________
1098 TF1*
1099 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
1100                                           Double_t sigman)
1101 {
1102   //
1103   // Fit a N-particle signal to the passed energy loss distribution
1104   //
1105   // If there's no 1-particle fit present, it does that first
1106   //
1107   // Parameters:
1108   //    dist   Data to fit the function to
1109   //    n      Number of particle signals to fit
1110   //    sigman If larger than zero, the initial guess of the
1111   //               detector induced noise. If zero or less, then this
1112   //               parameter is ignored in the fit (fixed at 0)
1113   //
1114   // Return:
1115   //    The function fitted to the data
1116   //
1117
1118   // Get the seed fit result
1119   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
1120   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
1121   if (!r || !f) {
1122     f = Fit1Particle(dist, sigman);
1123     r = static_cast<TFitResult*>(fFitResults.At(0));
1124     if (!r || !f) {
1125       ::Warning("FitNLandau", "No first shot at landau fit");
1126       return 0;
1127     }
1128   }
1129
1130   // Get some parameters from seed fit
1131   Double_t delta1  = r->Parameter(kDelta);
1132   Double_t xi1     = r->Parameter(kXi);
1133   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
1134   Double_t minE    = f->GetXmin();
1135
1136   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1137   Int_t    maxEb = dist->GetXaxis()->FindBin(maxEi);
1138   Double_t intg  = dist->Integral(minEb, maxEb);
1139   if (intg <= 0) {
1140     ::Warning("FitNParticle",
1141               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
1142               dist->GetName(), minE, maxEi, minEb, maxEb, intg);
1143     return 0;
1144   }
1145
1146   // Array of weights
1147   TArrayD a(n-1);
1148   for (UShort_t i = 2; i <= n; i++)
1149     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
1150   // Make the fit function
1151   TF1* landaun = MakeNLandauGaus(r->Parameter(kC),
1152                                  r->Parameter(kDelta),
1153                                  r->Parameter(kXi),
1154                                  r->Parameter(kSigma),
1155                                  r->Parameter(kSigmaN),
1156                                  n, a.fArray, minE, maxEi);
1157   if (minE      <= r->Parameter(kDelta) &&
1158       fMaxRange >= r->Parameter(kDelta)) {
1159     // Protect against warning from ParameterSettings
1160     // printf("FitN: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
1161     landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
1162   }
1163   if (r->Parameter(kXi) >= 0 && r->Parameter(kXi) <= 0.1) {
1164     // printf("FitN: Set par limits on xi: %f, %f\n", 0., 0.1);
1165     landaun->SetParLimits(kXi,     0.00, 0.1);  // was fMaxRange - too wide
1166   }
1167   if (r->Parameter(kSigma) >= 1e-5 && r->Parameter(kSigma) <= 0.1) {
1168     // printf("FitN: Set par limits on sigma: %f, %f\n", 1e-5, 0.1);
1169     landaun->SetParLimits(kSigma,  1e-5, 0.1);  // was fMaxRange - too wide
1170   }
1171   // Check if we're using the noise sigma
1172   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
1173   else {
1174     // printf("FitN: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
1175     landaun->SetParLimits(kSigmaN, 0, fMaxRange);
1176   }
1177
1178   // Set the range and name of the scale parameters
1179   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
1180     if (a[i-2] >= 0 && a[i-2] <= 1) {
1181       // printf("FitN: Set par limits on a_%d: %f, %f\n", i, 0., 1.);
1182       landaun->SetParLimits(kA+i-2, 0,1);
1183     }
1184   }
1185
1186   // Do the fit
1187   if (fDebug)
1188     ::Info("FitNParticle", "Fitting in the range %f,%f (%d)", minE, maxEi, n);
1189   TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
1190
1191   // landaun->SetRange(minE, fMaxRange);
1192   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
1193   fFunctions.AddAtAndExpand(landaun, n-1);
1194
1195   return landaun;
1196 }
1197 //____________________________________________________________________
1198 TF1*
1199 AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
1200 {
1201   //
1202   // Fit a composite particle signal to the passed energy loss
1203   // distribution
1204   //
1205   // Parameters:
1206   //    dist    Data to fit the function to
1207   //    sigman If larger than zero, the initial guess of the
1208   //               detector induced noise. If zero or less, then this
1209   //               parameter is ignored in the fit (fixed at 0)
1210   //
1211   // Return:
1212   //    The function fitted to the data
1213   //
1214
1215   // Find the fit range
1216   Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
1217   Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
1218                                 dist->GetNbinsX());
1219   dist->GetXaxis()->SetRange(cutBin, maxBin);
1220
1221   // Get the bin with maximum
1222   Int_t    peakBin = dist->GetMaximumBin();
1223   Double_t peakE   = dist->GetBinLowEdge(peakBin);
1224
1225   // Get the low edge
1226   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
1227   Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
1228   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
1229   Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
1230
1231   // Get the range in bins and the integral of that range
1232   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1233   Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
1234   Double_t intg  = dist->Integral(minEb, maxEb);
1235   if (intg <= 0) {
1236     ::Warning("Fit1Particle",
1237               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
1238               dist->GetName(), minE, maxE, minEb, maxEb, intg);
1239     return 0;
1240   }
1241
1242   // Restore the range
1243   dist->GetXaxis()->SetRange(1, maxBin);
1244
1245   // Define the function to fit
1246   TF1* seed = new TF1("landauSeed", landauGaus1, minE,maxE,kSigmaN+1);
1247
1248   // Set initial guesses, parameter names, and limits
1249   seed->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
1250   seed->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
1251   seed->SetNpx(500);
1252   seed->SetParLimits(kDelta, minE, fMaxRange);
1253   seed->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
1254   seed->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
1255   if (sigman <= 0)  seed->FixParameter(kSigmaN, 0);
1256   else              seed->SetParLimits(kSigmaN, 0, fMaxRange);
1257
1258   // Do the fit, getting the result object
1259   if (fDebug)
1260     ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
1261   /* TFitResultPtr r = */ dist->Fit(seed, "RNQS", "", minE, maxE);
1262
1263   maxE = dist->GetXaxis()->GetXmax();
1264   TF1* comp = new TF1("composite", landauGausComposite,
1265                       minE, maxE, kSigma+1+4);
1266   comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
1267                     "C#prime", "#Delta_{p}#prime", "#xi#prime", "#sigma#prim");
1268   comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight
1269                       seed->GetParameter(kDelta),    // 1 Primary Delta
1270                       seed->GetParameter(kDelta)/10, // 2 primary Xi
1271                       seed->GetParameter(kDelta)/5,  // 3 primary sigma
1272                       1.20 * seed->GetParameter(kC), // 5 Secondary weight
1273                       seed->GetParameter(kDelta),    // 6 secondary Delta
1274                       seed->GetParameter(kXi),       // 7 secondary Xi
1275                       seed->GetParameter(kSigma));   // 8 secondary sigma
1276
1277   // comp->SetParLimits(kC,       minE, fMaxRange); // C
1278   comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
1279   comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi
1280   comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
1281   // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
1282   comp->SetParLimits(kSigma+2,    minE/10, fMaxRange); // Delta
1283   comp->SetParLimits(kSigma+3,    0.00,    fMaxRange); // Xi
1284   comp->SetParLimits(kSigma+4,    1e-6,    fMaxRange); // Sigma
1285   comp->SetLineColor(kRed+1);
1286   comp->SetLineWidth(3);
1287
1288   // Do the fit, getting the result object
1289   if (fDebug)
1290     ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
1291   /* TFitResultPtr r = */ dist->Fit(comp, "RNQS", "", minE, maxE);
1292
1293 #if 0
1294   TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
1295   part1->SetLineColor(kGreen+1);
1296   part1->SetLineWidth(4);
1297   part1->SetRange(minE, maxE);
1298   part1->SetParameters(comp->GetParameter(0), // C
1299                        comp->GetParameter(1), // Delta
1300                        comp->GetParameter(2), // Xi
1301                        comp->GetParameter(3), // sigma
1302                        0);
1303   part1->Save(minE,maxE,0,0,0,0);
1304   dist->GetListOfFunctions()->Add(part1);
1305
1306   TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
1307   part2->SetLineColor(kBlue+1);
1308   part2->SetLineWidth(4);
1309   part2->SetRange(minE, maxE);
1310   part2->SetParameters(comp->GetParameter(4), // C
1311                        comp->GetParameter(5), // Delta
1312                        comp->GetParameter(6), // Xi
1313                        comp->GetParameter(7), // sigma
1314                        0);
1315   part2->Save(minE,maxE,0,0,0,0);
1316   dist->GetListOfFunctions()->Add(part2);
1317 #endif
1318   return comp;
1319 }
1320
1321 //====================================================================
1322 AliForwardUtil::Histos::~Histos()
1323 {
1324   //
1325   // Destructor
1326   //
1327 }
1328
1329 //____________________________________________________________________
1330 void
1331 AliForwardUtil::Histos::Delete(Option_t* opt)
1332 {
1333   if (fFMD1i) delete fFMD1i;
1334   if (fFMD2i) delete fFMD2i;
1335   if (fFMD2o) delete fFMD2o;
1336   if (fFMD3i) delete fFMD3i;
1337   if (fFMD3o) delete fFMD3o;
1338   fFMD1i = 0;
1339   fFMD2i = 0;
1340   fFMD2o = 0;
1341   fFMD3i = 0;
1342   fFMD3o = 0;
1343   TObject::Delete(opt);
1344 }
1345
1346 //____________________________________________________________________
1347 TH2D*
1348 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, const TAxis& etaAxis)
1349 {
1350   //
1351   // Make a histogram
1352   //
1353   // Parameters:
1354   //    d        Detector
1355   //    r        Ring
1356   //    etaAxis  Eta axis to use
1357   //
1358   // Return:
1359   //    Newly allocated histogram
1360   //
1361   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
1362   TH2D* hist = 0;
1363   if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1364     hist = new TH2D(Form("FMD%d%c_cache", d, r),
1365                     Form("FMD%d%c cache", d, r),
1366                     etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(),
1367                     ns, 0, TMath::TwoPi());
1368   else
1369     hist = new TH2D(Form("FMD%d%c_cache", d, r),
1370                     Form("FMD%d%c cache", d, r),
1371                     etaAxis.GetNbins(), etaAxis.GetXmin(),
1372                     etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
1373   hist->SetXTitle("#eta");
1374   hist->SetYTitle("#phi [radians]");
1375   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
1376   hist->Sumw2();
1377   hist->SetDirectory(0);
1378
1379   return hist;
1380 }
1381 //____________________________________________________________________
1382 void
1383 AliForwardUtil::Histos::RebinEta(TH2D* hist, const TAxis& etaAxis)
1384 {
1385   TAxis* xAxis = hist->GetXaxis();
1386   if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1387     xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
1388   else
1389     xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
1390   hist->Rebuild();
1391 }
1392
1393
1394 //____________________________________________________________________
1395 void
1396 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
1397 {
1398   //
1399   // Initialize the object
1400   //
1401   // Parameters:
1402   //    etaAxis Eta axis to use
1403   //
1404   fFMD1i = Make(1, 'I', etaAxis);
1405   fFMD2i = Make(2, 'I', etaAxis);
1406   fFMD2o = Make(2, 'O', etaAxis);
1407   fFMD3i = Make(3, 'I', etaAxis);
1408   fFMD3o = Make(3, 'O', etaAxis);
1409 }
1410 //____________________________________________________________________
1411 void
1412 AliForwardUtil::Histos::ReInit(const TAxis& etaAxis)
1413 {
1414   //
1415   // Initialize the object
1416   //
1417   // Parameters:
1418   //    etaAxis Eta axis to use
1419   //
1420   RebinEta(fFMD1i, etaAxis);
1421   RebinEta(fFMD2i, etaAxis);
1422   RebinEta(fFMD2o, etaAxis);
1423   RebinEta(fFMD3i, etaAxis);
1424   RebinEta(fFMD3o, etaAxis);
1425 }
1426
1427 //____________________________________________________________________
1428 void
1429 AliForwardUtil::Histos::Clear(Option_t* option)
1430 {
1431   //
1432   // Clear data
1433   //
1434   // Parameters:
1435   //    option Not used
1436   //
1437   if (fFMD1i) fFMD1i->Reset(option);
1438   if (fFMD2i) fFMD2i->Reset(option);
1439   if (fFMD2o) fFMD2o->Reset(option);
1440   if (fFMD3i) fFMD3i->Reset(option);
1441   if (fFMD3o) fFMD3o->Reset(option);
1442 }
1443
1444 //____________________________________________________________________
1445 TH2D*
1446 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
1447 {
1448   //
1449   // Get the histogram for a particular detector,ring
1450   //
1451   // Parameters:
1452   //    d Detector
1453   //    r Ring
1454   //
1455   // Return:
1456   //    Histogram for detector,ring or nul
1457   //
1458   switch (d) {
1459   case 1: return fFMD1i;
1460   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1461   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1462   }
1463   return 0;
1464 }
1465 //====================================================================
1466 TList*
1467 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
1468 {
1469   //
1470   // Define the outout list in @a d
1471   //
1472   // Parameters:
1473   //    d Where to put the output list
1474   //
1475   // Return:
1476   //    Newly allocated TList object or null
1477   //
1478   if (!d) return 0;
1479   TList* list = new TList;
1480   list->SetOwner();
1481   list->SetName(fName.Data());
1482   d->Add(list);
1483   return list;
1484 }
1485 //____________________________________________________________________
1486 TList*
1487 AliForwardUtil::RingHistos::GetOutputList(const TList* d) const
1488 {
1489   //
1490   // Get our output list from the container @a d
1491   //
1492   // Parameters:
1493   //    d where to get the output list from
1494   //
1495   // Return:
1496   //    The found TList or null
1497   //
1498   if (!d) return 0;
1499   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1500   return list;
1501 }
1502
1503 //____________________________________________________________________
1504 TH1*
1505 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1506 {
1507   //
1508   // Find a specific histogram in the source list @a d
1509   //
1510   // Parameters:
1511   //    d     (top)-container
1512   //    name  Name of histogram
1513   //
1514   // Return:
1515   //    Found histogram or null
1516   //
1517   return static_cast<TH1*>(d->FindObject(name));
1518 }
1519
1520 //====================================================================
1521 AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl,
1522                                        const char* format, ...)
1523   : fMsg("")
1524 {
1525   if (lvl < msgLvl) return;
1526   va_list ap;
1527   va_start(ap, format);
1528   Format(fMsg, format, ap);
1529   va_end(ap);
1530   Output(+1, fMsg);
1531 }
1532 //____________________________________________________________________
1533 AliForwardUtil::DebugGuard::~DebugGuard()
1534 {
1535   if (fMsg.IsNull()) return;
1536   Output(-1, fMsg);
1537 }
1538 //____________________________________________________________________
1539 void
1540 AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl,
1541                                     const char* format, ...)
1542 {
1543   if (lvl < msgLvl) return;
1544   TString msg;
1545   va_list ap;
1546   va_start(ap, format);
1547   Format(msg, format, ap);
1548   va_end(ap);
1549   Output(0, msg);
1550 }
1551
1552 //____________________________________________________________________
1553 void
1554 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1555 {
1556   static char buf[512];
1557   Int_t n = gROOT->GetDirLevel() + 2;
1558   for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1559   vsnprintf(&(buf[n]), 511-n, format, ap);
1560   buf[511] = '\0';
1561   out = buf;
1562 }
1563 //____________________________________________________________________
1564 void
1565 AliForwardUtil::DebugGuard::Output(int in, TString& msg)
1566 {
1567   msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
1568   AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1569   if      (in > 0) gROOT->IncreaseDirLevel();
1570   else if (in < 0) gROOT->DecreaseDirLevel();
1571 }
1572
1573
1574
1575 //
1576 // EOF
1577 //