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