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 "RNQS"
  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 Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
 691 Double_t AliForwardUtil::fgConvolutionNSigma = 5;
 692 namespace {
 693   //
 694   // The shift of the most probable value for the ROOT function TMath::Landau
 695   //
 696   const Double_t  mpshift  = -0.22278298;
 697   //
 698   // Integration normalisation
 699   //
 700   const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
 701
 702   //
 703   // Utility function to use in TF1 defintition
 704   //
 705   Double_t landauGaus1(Double_t* xp, Double_t* pp)
 706   {
 707     Double_t x        = xp[0];
 708     Double_t constant = pp[AliForwardUtil::ELossFitter::kC];
 709     Double_t delta    = pp[AliForwardUtil::ELossFitter::kDelta];
 710     Double_t xi       = pp[AliForwardUtil::ELossFitter::kXi];
 711     Double_t sigma    = pp[AliForwardUtil::ELossFitter::kSigma];
 712     Double_t sigmaN   = pp[AliForwardUtil::ELossFitter::kSigmaN];
 713
 714     return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigmaN);
 715   }
 716
 717   Double_t landauGausComposite(Double_t* xp, Double_t* pp)
 718   {
 719     Double_t x           = xp[0];
 720     Double_t cP          = pp[AliForwardUtil::ELossFitter::kC];
 721     Double_t deltaP      = pp[AliForwardUtil::ELossFitter::kDelta];
 722     Double_t xiP         = pp[AliForwardUtil::ELossFitter::kXi];
 723     Double_t sigmaP      = pp[AliForwardUtil::ELossFitter::kSigma];
 724     Double_t cS          = pp[AliForwardUtil::ELossFitter::kSigma+1];
 725     Double_t deltaS      = deltaP; // pp[AliForwardUtil::ELossFitter::kSigma+2];
 726     Double_t xiS         = pp[AliForwardUtil::ELossFitter::kSigma+2/*3*/];
 727     Double_t sigmaS      = sigmaP; // pp[AliForwardUtil::ELossFitter::kSigma+4];
 728
 729     return (cP * AliForwardUtil::LandauGaus(x,deltaP,xiP,sigmaP,0) +
 730             cS * AliForwardUtil::LandauGaus(x,deltaS,xiS,sigmaS,0));
 731   }
 732
 733   //
 734   // Utility function to use in TF1 defintition
 735   //
 736   Double_t landauGausN(Double_t* xp, Double_t* pp)
 737   {
 738     Double_t  x        = xp[0];
 739     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 740     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 741     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 742     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 743     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 744     Int_t     n        = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 745     Double_t* a        = &(pp[AliForwardUtil::ELossFitter::kA]);
 746
 747     return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigmaN,
 748                                                   n, a);
 749   }
 750   //
 751   // Utility function to use in TF1 defintition
 752   //
 753   Double_t landauGausI(Double_t* xp, Double_t* pp)
 754   {
 755     Double_t x         = xp[0];
 756     Double_t constant  = pp[AliForwardUtil::ELossFitter::kC];
 757     Double_t delta     = pp[AliForwardUtil::ELossFitter::kDelta];
 758     Double_t xi        = pp[AliForwardUtil::ELossFitter::kXi];
 759     Double_t sigma     = pp[AliForwardUtil::ELossFitter::kSigma];
 760     Double_t sigmaN    = pp[AliForwardUtil::ELossFitter::kSigmaN];
 761     Int_t    i         = Int_t(pp[AliForwardUtil::ELossFitter::kN]);
 762
 763     return constant * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
 764   }
 765
 766
 767 }
 768 //____________________________________________________________________
 769 Double_t
 770 AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
 771 {
 772   //
 773   // Calculate the shifted Landau
 774   // @f[
 775   //    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
 776   // @f]
 777   //
 778   // where @f$ f_{L}@f$ is the ROOT implementation of the Landau
 779   // distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
 780   // @f$\Delta=0,\xi=1@f$.
 781   //
 782   // Parameters:
 783   //    x      Where to evaluate @f$ f'_{L}@f$
 784   //    delta  Most probable value
 785   //    xi     The 'width' of the distribution
 786   //
 787   // Return:
 788   //    @f$ f'_{L}(x;\Delta,\xi) @f$
 789   //
 790   return TMath::Landau(x, delta - xi * mpshift, xi);
 791 }
 792 //____________________________________________________________________
 793 Double_t
 794 AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
 795                            Double_t sigma, Double_t sigmaN)
 796 {
 797   //
 798   // Calculate the value of a Landau convolved with a Gaussian
 799   //
 800   // @f[
 801   // f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
 802   //    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
 803   //    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
 804   // @f]
 805   //
 806   // where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
 807   // energy loss, @f$ \xi@f$ the width of the Landau, and
 808   // @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
 809   // variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling
 810   // noise in the detector.
 811   //
 812   // Note that this function uses the constants fgConvolutionSteps and
 813   // fgConvolutionNSigma
 814   //
 815   // References:
 816   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 817   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 818   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
 819   //
 820   // Parameters:
 821   //    x         where to evaluate @f$ f@f$
 822   //    delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 823   //    xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
 824   //    sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 825   //    sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
 826   //
 827   // Return:
 828   //    @f$ f@f$ evaluated at @f$ x@f$.
 829   //
 830   Double_t deltap = delta - xi * mpshift;
 831   Double_t sigma2 = sigmaN*sigmaN + sigma*sigma;
 832   Double_t sigma1 = sigmaN == 0 ? sigma : TMath::Sqrt(sigma2);
 833   Double_t xlow   = x - fgConvolutionNSigma * sigma1;
 834   Double_t xhigh  = x + fgConvolutionNSigma * sigma1;
 835   Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
 836   Double_t sum    = 0;
 837
 838   for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) {
 839     Double_t x1 = xlow  + (i - .5) * step;
 840     Double_t x2 = xhigh - (i - .5) * step;
 841
 842     sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
 843     sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
 844   }
 845   return step * sum * invSq2pi / sigma1;
 846 }
 847
 848 //____________________________________________________________________
 849 Double_t
 850 AliForwardUtil::ILandauGaus(Double_t x, Double_t delta, Double_t xi,
 851                             Double_t sigma, Double_t sigmaN, Int_t i)
 852 {
 853   //
 854   // Evaluate
 855   // @f[
 856   //    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
 857   // @f]
 858   // corresponding to @f$ i@f$ particles i.e., with the substitutions
 859   // @f{eqnarray*}{
 860   //    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))
 861   //    \xi       \rightarrow \xi_i       &=& i \xi
 862   //    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma
 863   //    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
 864   // @f}
 865   //
 866   // Parameters:
 867   //    x        Where to evaluate
 868   //    delta    @f$ \Delta@f$
 869   //    xi       @f$ \xi@f$
 870   //    sigma    @f$ \sigma@f$
 871   //    sigma_n  @f$ \sigma_n@f$
 872   //    i        @f$ i @f$
 873   //
 874   // Return:
 875   //    @f$ f_i @f$ evaluated
 876   //
 877   Double_t deltaI =  (i == 1 ? delta : i * (delta + xi * TMath::Log(i)));
 878   Double_t xiI    =  i * xi;
 879   Double_t sigmaI =  (i == 1 ? sigma : TMath::Sqrt(Double_t(i))*sigma);
 880   if (sigmaI < 1e-10) {
 881     // Fall back to landau
 882     return AliForwardUtil::Landau(x, deltaI, xiI);
 883   }
 884   return AliForwardUtil::LandauGaus(x, deltaI, xiI, sigmaI, sigmaN);
 885 }
 886
 887 //____________________________________________________________________
 888 Double_t
 889 AliForwardUtil::IdLandauGausdPar(Double_t x,
 890                                  UShort_t par,   Double_t dPar,
 891                                  Double_t delta, Double_t xi,
 892                                  Double_t sigma, Double_t sigmaN,
 893                                  Int_t    i)
 894 {
 895   //
 896   // Numerically evaluate
 897   // @f[
 898   //    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
 899   // @f]
 900   // where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping
 901   // of the parameters is given by
 902   //
 903   // - 0: @f$\Delta@f$
 904   // - 1: @f$\xi@f$
 905   // - 2: @f$\sigma@f$
 906   // - 3: @f$\sigma_n@f$
 907   //
 908   // This is the partial derivative with respect to the parameter of
 909   // the response function corresponding to @f$ i@f$ particles i.e.,
 910   // with the substitutions
 911   // @f[
 912   //    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))
 913   //    \xi       \rightarrow \xi_i       = i \xi
 914   //    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma
 915   //    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
 916   // @f]
 917   //
 918   // Parameters:
 919   //    x        Where to evaluate
 920   //    ipar     Parameter number
 921   //    dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
 922   //    delta    @f$ \Delta@f$
 923   //    xi       @f$ \xi@f$
 924   //    sigma    @f$ \sigma@f$
 925   //    sigma_n  @f$ \sigma_n@f$
 926   //    i        @f$ i@f$
 927   //
 928   // Return:
 929   //    @f$ f_i@f$ evaluated
 930   //
 931   if (dPar == 0) return 0;
 932   Double_t dp      = dPar;
 933   Double_t d2      = dPar / 2;
 934   Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
 935   Double_t xiI     =  i * xi;
 936   Double_t si      =  TMath::Sqrt(Double_t(i));
 937   Double_t sigmaI  =  si*sigma;
 938   Double_t y1      = 0;
 939   Double_t y2      = 0;
 940   Double_t y3      = 0;
 941   Double_t y4      = 0;
 942   switch (par) {
 943   case 0:
 944     y1 = ILandauGaus(x, deltaI+i*dp, xiI, sigmaI, sigmaN, i);
 945     y2 = ILandauGaus(x, deltaI+i*d2, xiI, sigmaI, sigmaN, i);
 946     y3 = ILandauGaus(x, deltaI-i*d2, xiI, sigmaI, sigmaN, i);
 947     y4 = ILandauGaus(x, deltaI-i*dp, xiI, sigmaI, sigmaN, i);
 948     break;
 949   case 1:
 950     y1 = ILandauGaus(x, deltaI, xiI+i*dp, sigmaI, sigmaN, i);
 951     y2 = ILandauGaus(x, deltaI, xiI+i*d2, sigmaI, sigmaN, i);
 952     y3 = ILandauGaus(x, deltaI, xiI-i*d2, sigmaI, sigmaN, i);
 953     y4 = ILandauGaus(x, deltaI, xiI-i*dp, sigmaI, sigmaN, i);
 954     break;
 955   case 2:
 956     y1 = ILandauGaus(x, deltaI, xiI, sigmaI+si*dp, sigmaN, i);
 957     y2 = ILandauGaus(x, deltaI, xiI, sigmaI+si*d2, sigmaN, i);
 958     y3 = ILandauGaus(x, deltaI, xiI, sigmaI-si*d2, sigmaN, i);
 959     y4 = ILandauGaus(x, deltaI, xiI, sigmaI-si*dp, sigmaN, i);
 960     break;
 961   case 3:
 962     y1 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+dp, i);
 963     y2 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN+d2, i);
 964     y3 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-d2, i);
 965     y4 = ILandauGaus(x, deltaI, xiI, sigmaI, sigmaN-dp, i);
 966     break;
 967   default:
 968     return 0;
 969   }
 970
 971   Double_t d0  = y1 - y4;
 972   Double_t d1  = 2 * (y2 - y3);
 973
 974   Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
 975
 976   return g;
 977 }
 978
 979 //____________________________________________________________________
 980 Double_t
 981 AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi,
 982                             Double_t sigma, Double_t sigmaN, Int_t n,
 983                             const Double_t* a)
 984 {
 985   //
 986   // Evaluate
 987   // @f[
 988   //   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
 989   // @f]
 990   //
 991   // where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
 992   // Landau with a Gaussian (see LandauGaus).  Note that
 993   // @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$,
 994   // @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$.
 995   //
 996   // References:
 997   //  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
 998   //  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
 999   //  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
1000   //
1001   // Parameters:
1002   //    x        Where to evaluate @f$ f_N@f$
1003   //    delta    @f$ \Delta_1@f$
1004   //    xi       @f$ \xi_1@f$
1005   //    sigma    @f$ \sigma_1@f$
1006   //    sigma_n  @f$ \sigma_n@f$
1007   //    n        @f$ N@f$ in the sum above.
1008   //    a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
1009   //                 @f$ i > 1@f$
1010   //
1011   // Return:
1012   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$
1013   //
1014   Double_t result = ILandauGaus(x, delta, xi, sigma, sigmaN, 1);
1015   for (Int_t i = 2; i <= n; i++)
1016     result += a[i-2] * AliForwardUtil::ILandauGaus(x,delta,xi,sigma,sigmaN,i);
1017   return result;
1018 }
1019 namespace {
1020   const Int_t kColors[] = { kRed+1,
1021                             kPink+3,
1022                             kMagenta+2,
1023                             kViolet+2,
1024                             kBlue+1,
1025                             kAzure+3,
1026                             kCyan+1,
1027                             kTeal+2,
1028                             kGreen+2,
1029                             kSpring+3,
1030                             kYellow+2,
1031                             kOrange+2 };
1032 }
1033
1034 //____________________________________________________________________
1035 TF1*
1036 AliForwardUtil::MakeNLandauGaus(Double_t  c,
1037                                 Double_t  delta, Double_t xi,
1038                                 Double_t  sigma, Double_t sigmaN, Int_t n,
1039                                 const Double_t* a,
1040                                 Double_t  xmin, Double_t xmax)
1041 {
1042   //
1043   // Generate a TF1 object of @f$ f_N@f$
1044   //
1045   // Parameters:
1046   //    c         Constant
1047   //    delta     @f$ \Delta@f$
1048   //    xi            @f$ \xi_1@f$
1049   //    sigma     @f$ \sigma_1@f$
1050   //    sigma_n   @f$ \sigma_n@f$
1051   //    n             @f$ N@f$ - how many particles to sum to
1052   //    a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for
1053   //                  @f$ i > 1@f$
1054   //    xmin      Least value of range
1055   //    xmax      Largest value of range
1056   //
1057   // Return:
1058   //    Newly allocated TF1 object
1059   //
1060   Int_t npar       = AliForwardUtil::ELossFitter::kN+n;
1061   TF1* landaun     = new TF1(Form("nlandau%d", n), &landauGausN,xmin,xmax,npar);
1062   // landaun->SetLineStyle(((n-2) % 10)+2); // start at dashed
1063   landaun->SetLineColor(kColors[((n-1) % 12)]); // start at red
1064   landaun->SetLineWidth(2);
1065   landaun->SetNpx(500);
1066   landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
1067
1068   // Set the initial parameters from the seed fit
1069   landaun->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
1070   landaun->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
1071   landaun->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
1072   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
1073   landaun->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
1074   landaun->FixParameter(AliForwardUtil::ELossFitter::kN,      n);
1075
1076   // Set the range and name of the scale parameters
1077   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
1078     landaun->SetParameter(AliForwardUtil::ELossFitter::kA+i-2, a[i-2]);
1079     landaun->SetParName(AliForwardUtil::ELossFitter::kA+i-2, Form("a_{%d}", i));
1080   }
1081   return landaun;
1082 }
1083 //____________________________________________________________________
1084 TF1*
1085 AliForwardUtil::MakeILandauGaus(Double_t  c,
1086                                 Double_t  delta, Double_t xi,
1087                                 Double_t  sigma, Double_t sigmaN, Int_t i,
1088                                 Double_t  xmin, Double_t xmax)
1089 {
1090   //
1091   // Generate a TF1 object of @f$ f_I@f$
1092   //
1093   // Parameters:
1094   //    c        Constant
1095   //    delta    @f$ \Delta@f$
1096   //    xi       @f$ \xi_1@f$
1097   //    sigma    @f$ \sigma_1@f$
1098   //    sigma_n  @f$ \sigma_n@f$
1099   //    i            @f$ i@f$ - the number of particles
1100   //    xmin     Least value of range
1101   //    xmax     Largest value of range
1102   //
1103   // Return:
1104   //    Newly allocated TF1 object
1105   //
1106   Int_t npar       = AliForwardUtil::ELossFitter::kN+1;
1107   TF1* landaui     = new TF1(Form("ilandau%d", i), &landauGausI,xmin,xmax,npar);
1108   // landaui->SetLineStyle(((i-2) % 10)+2); // start at dashed
1109   landaui->SetLineColor(kColors[((i-1) % 12)]); // start at red
1110   landaui->SetLineWidth(1);
1111   landaui->SetNpx(500);
1112   landaui->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "i");
1113
1114   // Set the initial parameters from the seed fit
1115   landaui->SetParameter(AliForwardUtil::ELossFitter::kC,      c);
1116   landaui->SetParameter(AliForwardUtil::ELossFitter::kDelta,  delta);
1117   landaui->SetParameter(AliForwardUtil::ELossFitter::kXi,     xi);
1118   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigma,  sigma);
1119   landaui->SetParameter(AliForwardUtil::ELossFitter::kSigmaN, sigmaN);
1120   landaui->FixParameter(AliForwardUtil::ELossFitter::kN,      i);
1121
1122   return landaui;
1123 }
1124
1125 //====================================================================
1126 AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut,
1127                                          Double_t maxRange,
1128                                          UShort_t minusBins)
1129   : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
1130     fFitResults(0), fFunctions(0), fDebug(false)
1131 {
1132   //
1133   // Constructor
1134   //
1135   // Parameters:
1136   //    lowCut     Lower cut of spectrum - data below this cuts is ignored
1137   //    maxRange   Maximum range to fit to
1138   //    minusBins  The number of bins below maximum to use
1139   //
1140   fFitResults.SetOwner();
1141   fFunctions.SetOwner();
1142 }
1143 //____________________________________________________________________
1144 AliForwardUtil::ELossFitter::~ELossFitter()
1145 {
1146   //
1147   // Destructor
1148   //
1149   //
1150   fFitResults.Delete();
1151   fFunctions.Delete();
1152 }
1153 //____________________________________________________________________
1154 void
1155 AliForwardUtil::ELossFitter::Clear()
1156 {
1157   //
1158   // Clear internal arrays
1159   //
1160   //
1161   fFitResults.Clear();
1162   fFunctions.Clear();
1163 }
1164 //____________________________________________________________________
1165 TF1*
1166 AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
1167 {
1168   //
1169   // Fit a 1-particle signal to the passed energy loss distribution
1170   //
1171   // Note that this function clears the internal arrays first
1172   //
1173   // Parameters:
1174   //    dist    Data to fit the function to
1175   //    sigman If larger than zero, the initial guess of the
1176   //               detector induced noise. If zero or less, then this
1177   //               parameter is ignored in the fit (fixed at 0)
1178   //
1179   // Return:
1180   //    The function fitted to the data
1181   //
1182
1183   // Clear the cache
1184   Clear();
1185
1186   // Find the fit range
1187   // Find the fit range
1188   Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
1189   Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
1190                                 dist->GetNbinsX());
1191   dist->GetXaxis()->SetRange(cutBin, maxBin);
1192   // dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
1193
1194   // Get the bin with maximum
1195   Int_t    peakBin = dist->GetMaximumBin();
1196   Double_t peakE   = dist->GetBinLowEdge(peakBin);
1197
1198   // Get the low edge
1199   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
1200   Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
1201   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
1202   Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
1203
1204   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1205   Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
1206   Double_t intg  = dist->Integral(minEb, maxEb);
1207   if (intg <= 0) {
1208     ::Warning("Fit1Particle",
1209               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
1210               dist->GetName(), minE, maxE, minEb, maxEb, intg);
1211     return 0;
1212   }
1213
1214   // Restore the range
1215   dist->GetXaxis()->SetRange(1, maxBin);
1216
1217   // Define the function to fit
1218   TF1* landau1 = new TF1("landau1", landauGaus1, minE,maxE,kSigmaN+1);
1219
1220   // Set initial guesses, parameter names, and limits
1221   landau1->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
1222   landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
1223   landau1->SetNpx(500);
1224   if (peakE >= minE && peakE <= fMaxRange) {
1225     // printf("Fit1: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
1226     landau1->SetParLimits(kDelta, minE, fMaxRange);
1227   }
1228   if (peakE/10 >= 0 && peakE <= 0.1) {
1229     // printf("Fit1: Set par limits on xi: %f, %f\n", 0., 0.1);
1230     landau1->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
1231   }
1232   if (peakE/5 >= 0 && peakE/5 <= 0.1) {
1233     // printf("Fit1: Set par limits on sigma: %f, %f\n", 0., 0.1);
1234     landau1->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
1235   }
1236   if (sigman <= 0)  landau1->FixParameter(kSigmaN, 0);
1237   else {
1238     // printf("Fit1: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
1239     landau1->SetParLimits(kSigmaN, 0, fMaxRange);
1240   }
1241
1242   // Do the fit, getting the result object
1243   if (fDebug)
1244     ::Info("Fit1Particle", "Fitting in the range %f,%f", minE, maxE);
1245   TFitResultPtr r = dist->Fit(landau1, FIT_OPTIONS, "", minE, maxE);
1246   if (!r.Get()) {
1247     ::Warning("Fit1Particle",
1248               "No fit returned when processing %s in the range [%f,%f] "
1249               "options %s", dist->GetName(), minE, maxE, FIT_OPTIONS);
1250     return 0;
1251   }
1252   // landau1->SetRange(minE, fMaxRange);
1253   fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
1254   fFunctions.AddAtAndExpand(landau1, 0);
1255
1256   return landau1;
1257 }
1258 //____________________________________________________________________
1259 TF1*
1260 AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n,
1261                                           Double_t sigman)
1262 {
1263   //
1264   // Fit a N-particle signal to the passed energy loss distribution
1265   //
1266   // If there's no 1-particle fit present, it does that first
1267   //
1268   // Parameters:
1269   //    dist   Data to fit the function to
1270   //    n      Number of particle signals to fit
1271   //    sigman If larger than zero, the initial guess of the
1272   //               detector induced noise. If zero or less, then this
1273   //               parameter is ignored in the fit (fixed at 0)
1274   //
1275   // Return:
1276   //    The function fitted to the data
1277   //
1278
1279   // Get the seed fit result
1280   TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
1281   TF1*        f = static_cast<TF1*>(fFunctions.At(0));
1282   if (!r || !f) {
1283     f = Fit1Particle(dist, sigman);
1284     r = static_cast<TFitResult*>(fFitResults.At(0));
1285     if (!r || !f) {
1286       ::Warning("FitNLandau", "No first shot at landau fit");
1287       return 0;
1288     }
1289   }
1290
1291   // Get some parameters from seed fit
1292   Double_t delta1  = r->Parameter(kDelta);
1293   Double_t xi1     = r->Parameter(kXi);
1294   Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
1295   Double_t minE    = f->GetXmin();
1296
1297   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1298   Int_t    maxEb = dist->GetXaxis()->FindBin(maxEi);
1299   Double_t intg  = dist->Integral(minEb, maxEb);
1300   if (intg <= 0) {
1301     ::Warning("FitNParticle",
1302               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
1303               dist->GetName(), minE, maxEi, minEb, maxEb, intg);
1304     return 0;
1305   }
1306
1307   // Array of weights
1308   TArrayD a(n-1);
1309   for (UShort_t i = 2; i <= n; i++)
1310     a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
1311   // Make the fit function
1312   TF1* landaun = MakeNLandauGaus(r->Parameter(kC),
1313                                  r->Parameter(kDelta),
1314                                  r->Parameter(kXi),
1315                                  r->Parameter(kSigma),
1316                                  r->Parameter(kSigmaN),
1317                                  n, a.fArray, minE, maxEi);
1318   if (minE      <= r->Parameter(kDelta) &&
1319       fMaxRange >= r->Parameter(kDelta)) {
1320     // Protect against warning from ParameterSettings
1321     // printf("FitN: Set par limits on Delta: %f, %f\n", minE, fMaxRange);
1322     landaun->SetParLimits(kDelta,  minE, fMaxRange);       // Delta
1323   }
1324   if (r->Parameter(kXi) >= 0 && r->Parameter(kXi) <= 0.1) {
1325     // printf("FitN: Set par limits on xi: %f, %f\n", 0., 0.1);
1326     landaun->SetParLimits(kXi,     0.00, 0.1);  // was fMaxRange - too wide
1327   }
1328   if (r->Parameter(kSigma) >= 1e-5 && r->Parameter(kSigma) <= 0.1) {
1329     // printf("FitN: Set par limits on sigma: %f, %f\n", 1e-5, 0.1);
1330     landaun->SetParLimits(kSigma,  1e-5, 0.1);  // was fMaxRange - too wide
1331   }
1332   // Check if we're using the noise sigma
1333   if (sigman <= 0)  landaun->FixParameter(kSigmaN, 0);
1334   else {
1335     // printf("FitN: Set par limits on sigmaN: %f, %f\n", 0., fMaxRange);
1336     landaun->SetParLimits(kSigmaN, 0, fMaxRange);
1337   }
1338
1339   // Set the range and name of the scale parameters
1340   for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
1341     if (a[i-2] >= 0 && a[i-2] <= 1) {
1342       // printf("FitN: Set par limits on a_%d: %f, %f\n", i, 0., 1.);
1343       landaun->SetParLimits(kA+i-2, 0,1);
1344     }
1345   }
1346
1347   // Do the fit
1348   if (fDebug)
1349     ::Info("FitNParticle", "Fitting in the range %f,%f (%d)", minE, maxEi, n);
1350   TFitResultPtr tr = dist->Fit(landaun, FIT_OPTIONS, "", minE, maxEi);
1351
1352   // landaun->SetRange(minE, fMaxRange);
1353   fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
1354   fFunctions.AddAtAndExpand(landaun, n-1);
1355
1356   return landaun;
1357 }
1358 //____________________________________________________________________
1359 TF1*
1360 AliForwardUtil::ELossFitter::FitComposite(TH1* dist, Double_t sigman)
1361 {
1362   //
1363   // Fit a composite particle signal to the passed energy loss
1364   // distribution
1365   //
1366   // Parameters:
1367   //    dist    Data to fit the function to
1368   //    sigman If larger than zero, the initial guess of the
1369   //               detector induced noise. If zero or less, then this
1370   //               parameter is ignored in the fit (fixed at 0)
1371   //
1372   // Return:
1373   //    The function fitted to the data
1374   //
1375
1376   // Find the fit range
1377   Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
1378   Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
1379                                 dist->GetNbinsX());
1380   dist->GetXaxis()->SetRange(cutBin, maxBin);
1381
1382   // Get the bin with maximum
1383   Int_t    peakBin = dist->GetMaximumBin();
1384   Double_t peakE   = dist->GetBinLowEdge(peakBin);
1385
1386   // Get the low edge
1387   // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
1388   Int_t    minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
1389   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
1390   Double_t maxE   = dist->GetBinCenter(peakBin+2*fMinusBins);
1391
1392   // Get the range in bins and the integral of that range
1393   Int_t    minEb = dist->GetXaxis()->FindBin(minE);
1394   Int_t    maxEb = dist->GetXaxis()->FindBin(maxE);
1395   Double_t intg  = dist->Integral(minEb, maxEb);
1396   if (intg <= 0) {
1397     ::Warning("Fit1Particle",
1398               "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
1399               dist->GetName(), minE, maxE, minEb, maxEb, intg);
1400     return 0;
1401   }
1402
1403   // Restore the range
1404   dist->GetXaxis()->SetRange(1, maxBin);
1405
1406   // Define the function to fit
1407   TF1* seed = new TF1("landauSeed", landauGaus1, minE,maxE,kSigmaN+1);
1408
1409   // Set initial guesses, parameter names, and limits
1410   seed->SetParameters(1,peakE,peakE/10,peakE/5,sigman);
1411   seed->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
1412   seed->SetNpx(500);
1413   seed->SetParLimits(kDelta, minE, fMaxRange);
1414   seed->SetParLimits(kXi,    0.00, 0.1); // Was fMaxRange - too wide
1415   seed->SetParLimits(kSigma, 1e-5, 0.1); // Was fMaxRange - too wide
1416   if (sigman <= 0)  seed->FixParameter(kSigmaN, 0);
1417   else              seed->SetParLimits(kSigmaN, 0, fMaxRange);
1418
1419   // Do the fit, getting the result object
1420   if (fDebug)
1421     ::Info("FitComposite", "Fitting seed in the range %f,%f", minE, maxE);
1422   /* TFitResultPtr r = */ dist->Fit(seed, FIT_OPTIONS, "", minE, maxE);
1423
1424   maxE = dist->GetXaxis()->GetXmax();
1425 #if 1
1426   TF1* comp = new TF1("composite", landauGausComposite,
1427                       minE, maxE, kSigma+1+2);
1428   comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
1429                     "C#prime", "#xi#prime");
1430   comp->SetParameters(0.8 * seed->GetParameter(kC),  // 0 Primary weight
1431                       seed->GetParameter(kDelta),    // 1 Primary Delta
1432                       seed->GetParameter(kDelta)/10, // 2 primary Xi
1433                       seed->GetParameter(kDelta)/5,  // 3 primary sigma
1434                       1.20 * seed->GetParameter(kC), // 5 Secondary weight
1435                       seed->GetParameter(kXi));      // 7 secondary Xi
1436   // comp->SetParLimits(kC,       minE, fMaxRange); // C
1437   comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
1438   comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi
1439   comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
1440   // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
1441   comp->SetParLimits(kSigma+2,    0.00, fMaxRange); // Xi'
1442 #else
1443   TF1* comp = new TF1("composite", landauGausComposite,
1444                       minE, maxE, kSigma+1+4);
1445   comp->SetParNames("C",       "#Delta_{p}",       "#xi",       "#sigma",
1446                     "C#prime", "#Delta_{p}#prime", "#xi#prime", "#sigma#prim");
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(kDelta),    // 6 secondary Delta
1453                       seed->GetParameter(kXi),       // 7 secondary Xi
1454                       seed->GetParameter(kSigma));   // 8 secondary sigma
1455   // comp->SetParLimits(kC,       minE, fMaxRange); // C
1456   comp->SetParLimits(kDelta,      minE, fMaxRange); // Delta
1457   comp->SetParLimits(kXi,         0.00, fMaxRange); // Xi
1458   comp->SetParLimits(kSigma,      1e-5, fMaxRange); // Sigma
1459   // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
1460   comp->SetParLimits(kSigma+2,    minE/10, fMaxRange); // Delta
1461   comp->SetParLimits(kSigma+3,    0.00,    fMaxRange); // Xi
1462   comp->SetParLimits(kSigma+4,    1e-6,    fMaxRange); // Sigma
1463 #endif
1464   comp->SetLineColor(kRed+1);
1465   comp->SetLineWidth(3);
1466
1467   // Do the fit, getting the result object
1468   if (fDebug)
1469     ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
1470   /* TFitResultPtr r = */ dist->Fit(comp, FIT_OPTIONS, "", minE, maxE);
1471
1472 #if 0
1473   TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
1474   part1->SetLineColor(kGreen+1);
1475   part1->SetLineWidth(4);
1476   part1->SetRange(minE, maxE);
1477   part1->SetParameters(comp->GetParameter(0), // C
1478                        comp->GetParameter(1), // Delta
1479                        comp->GetParameter(2), // Xi
1480                        comp->GetParameter(3), // sigma
1481                        0);
1482   part1->Save(minE,maxE,0,0,0,0);
1483   dist->GetListOfFunctions()->Add(part1);
1484
1485   TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
1486   part2->SetLineColor(kBlue+1);
1487   part2->SetLineWidth(4);
1488   part2->SetRange(minE, maxE);
1489   part2->SetParameters(comp->GetParameter(4), // C
1490                        comp->GetParameter(5), // Delta
1491                        comp->GetParameter(6), // Xi
1492                        comp->GetParameter(7), // sigma
1493                        0);
1494   part2->Save(minE,maxE,0,0,0,0);
1495   dist->GetListOfFunctions()->Add(part2);
1496 #endif
1497   return comp;
1498 }
1499
1500 //====================================================================
1501 AliForwardUtil::Histos::~Histos()
1502 {
1503   //
1504   // Destructor
1505   //
1506 }
1507
1508 //____________________________________________________________________
1509 void
1510 AliForwardUtil::Histos::Delete(Option_t* opt)
1511 {
1512   if (fFMD1i) delete fFMD1i;
1513   if (fFMD2i) delete fFMD2i;
1514   if (fFMD2o) delete fFMD2o;
1515   if (fFMD3i) delete fFMD3i;
1516   if (fFMD3o) delete fFMD3o;
1517   fFMD1i = 0;
1518   fFMD2i = 0;
1519   fFMD2o = 0;
1520   fFMD3i = 0;
1521   fFMD3o = 0;
1522   TObject::Delete(opt);
1523 }
1524
1525 //____________________________________________________________________
1526 TH2D*
1527 AliForwardUtil::Histos::Make(UShort_t d, Char_t r, const TAxis& etaAxis)
1528 {
1529   //
1530   // Make a histogram
1531   //
1532   // Parameters:
1533   //    d        Detector
1534   //    r        Ring
1535   //    etaAxis  Eta axis to use
1536   //
1537   // Return:
1538   //    Newly allocated histogram
1539   //
1540   Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
1541   TH2D* hist = 0;
1542   if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1543     hist = new TH2D(Form("FMD%d%c_cache", d, r),
1544                     Form("FMD%d%c cache", d, r),
1545                     etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray(),
1546                     ns, 0, TMath::TwoPi());
1547   else
1548     hist = new TH2D(Form("FMD%d%c_cache", d, r),
1549                     Form("FMD%d%c cache", d, r),
1550                     etaAxis.GetNbins(), etaAxis.GetXmin(),
1551                     etaAxis.GetXmax(), ns, 0, TMath::TwoPi());
1552   hist->SetXTitle("#eta");
1553   hist->SetYTitle("#phi [radians]");
1554   hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
1555   hist->Sumw2();
1556   hist->SetDirectory(0);
1557
1558   return hist;
1559 }
1560 //____________________________________________________________________
1561 void
1562 AliForwardUtil::Histos::RebinEta(TH2D* hist, const TAxis& etaAxis)
1563 {
1564   TAxis* xAxis = hist->GetXaxis();
1565   if (etaAxis.GetXbins() && etaAxis.GetXbins()->GetArray())
1566     xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXbins()->GetArray());
1567   else
1568     xAxis->Set(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax());
1569   hist->Rebuild();
1570 }
1571
1572
1573 //____________________________________________________________________
1574 void
1575 AliForwardUtil::Histos::Init(const TAxis& etaAxis)
1576 {
1577   //
1578   // Initialize the object
1579   //
1580   // Parameters:
1581   //    etaAxis Eta axis to use
1582   //
1583   fFMD1i = Make(1, 'I', etaAxis);
1584   fFMD2i = Make(2, 'I', etaAxis);
1585   fFMD2o = Make(2, 'O', etaAxis);
1586   fFMD3i = Make(3, 'I', etaAxis);
1587   fFMD3o = Make(3, 'O', etaAxis);
1588 }
1589 //____________________________________________________________________
1590 void
1591 AliForwardUtil::Histos::ReInit(const TAxis& etaAxis)
1592 {
1593   //
1594   // Initialize the object
1595   //
1596   // Parameters:
1597   //    etaAxis Eta axis to use
1598   //
1599   if (!fFMD1i) fFMD1i = Make(1, 'i', etaAxis); else RebinEta(fFMD1i, etaAxis);
1600   if (!fFMD2i) fFMD2i = Make(2, 'i', etaAxis); else RebinEta(fFMD2i, etaAxis);
1601   if (!fFMD2o) fFMD2o = Make(2, 'o', etaAxis); else RebinEta(fFMD2o, etaAxis);
1602   if (!fFMD3i) fFMD3i = Make(3, 'i', etaAxis); else RebinEta(fFMD3i, etaAxis);
1603   if (!fFMD3o) fFMD3o = Make(3, 'o', etaAxis); else RebinEta(fFMD3o, etaAxis);
1604 }
1605
1606 //____________________________________________________________________
1607 void
1608 AliForwardUtil::Histos::Clear(Option_t* option)
1609 {
1610   //
1611   // Clear data
1612   //
1613   // Parameters:
1614   //    option Not used
1615   //
1616   if (fFMD1i) { fFMD1i->Reset(option); fFMD1i->ResetBit(kSkipRing); }
1617   if (fFMD2i) { fFMD2i->Reset(option); fFMD2i->ResetBit(kSkipRing); }
1618   if (fFMD2o) { fFMD2o->Reset(option); fFMD2o->ResetBit(kSkipRing); }
1619   if (fFMD3i) { fFMD3i->Reset(option); fFMD3i->ResetBit(kSkipRing); }
1620   if (fFMD3o) { fFMD3o->Reset(option); fFMD3o->ResetBit(kSkipRing); }
1621 }
1622
1623 //____________________________________________________________________
1624 TH2D*
1625 AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
1626 {
1627   //
1628   // Get the histogram for a particular detector,ring
1629   //
1630   // Parameters:
1631   //    d Detector
1632   //    r Ring
1633   //
1634   // Return:
1635   //    Histogram for detector,ring or nul
1636   //
1637   switch (d) {
1638   case 1: return fFMD1i;
1639   case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
1640   case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
1641   }
1642   return 0;
1643 }
1644 //====================================================================
1645 TList*
1646 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
1647 {
1648   //
1649   // Define the outout list in @a d
1650   //
1651   // Parameters:
1652   //    d Where to put the output list
1653   //
1654   // Return:
1655   //    Newly allocated TList object or null
1656   //
1657   if (!d) return 0;
1658   TList* list = new TList;
1659   list->SetOwner();
1660   list->SetName(fName.Data());
1661   d->Add(list);
1662   return list;
1663 }
1664 //____________________________________________________________________
1665 TList*
1666 AliForwardUtil::RingHistos::GetOutputList(const TList* d) const
1667 {
1668   //
1669   // Get our output list from the container @a d
1670   //
1671   // Parameters:
1672   //    d where to get the output list from
1673   //
1674   // Return:
1675   //    The found TList or null
1676   //
1677   if (!d) return 0;
1678   TList* list = static_cast<TList*>(d->FindObject(fName.Data()));
1679   return list;
1680 }
1681
1682 //____________________________________________________________________
1683 TH1*
1684 AliForwardUtil::RingHistos::GetOutputHist(const TList* d, const char* name) const
1685 {
1686   //
1687   // Find a specific histogram in the source list @a d
1688   //
1689   // Parameters:
1690   //    d     (top)-container
1691   //    name  Name of histogram
1692   //
1693   // Return:
1694   //    Found histogram or null
1695   //
1696   return static_cast<TH1*>(d->FindObject(name));
1697 }
1698
1699 //====================================================================
1700 AliForwardUtil::DebugGuard::DebugGuard(Int_t lvl, Int_t msgLvl,
1701                                        const char* format, ...)
1702   : fMsg("")
1703 {
1704   if (lvl < msgLvl) return;
1705   va_list ap;
1706   va_start(ap, format);
1707   Format(fMsg, format, ap);
1708   va_end(ap);
1709   Output(+1, fMsg);
1710 }
1711 //____________________________________________________________________
1712 AliForwardUtil::DebugGuard::~DebugGuard()
1713 {
1714   if (fMsg.IsNull()) return;
1715   Output(-1, fMsg);
1716 }
1717 //____________________________________________________________________
1718 void
1719 AliForwardUtil::DebugGuard::Message(Int_t lvl, Int_t msgLvl,
1720                                     const char* format, ...)
1721 {
1722   if (lvl < msgLvl) return;
1723   TString msg;
1724   va_list ap;
1725   va_start(ap, format);
1726   Format(msg, format, ap);
1727   va_end(ap);
1728   Output(0, msg);
1729 }
1730
1731 //____________________________________________________________________
1732 void
1733 AliForwardUtil::DebugGuard::Format(TString& out, const char* format, va_list ap)
1734 {
1735   static char buf[512];
1736   Int_t n = gROOT->GetDirLevel() + 2;
1737   for (Int_t i = 0; i < n; i++) buf[i] = ' ';
1738   vsnprintf(&(buf[n]), 511-n, format, ap);
1739   buf[511] = '\0';
1740   out = buf;
1741 }
1742 //____________________________________________________________________
1743 void
1744 AliForwardUtil::DebugGuard::Output(int in, TString& msg)
1745 {
1746   msg[0] = (in > 0 ? '>' :  in < 0 ? '<' : '=');
1747   AliLog::Message(AliLog::kInfo, msg, 0, 0, "PWGLF/forward", 0, 0);
1748   if      (in > 0) gROOT->IncreaseDirLevel();
1749   else if (in < 0) gROOT->DecreaseDirLevel();
1750 }
1751
1752
1753
1754 //
1755 // EOF
1756 //