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