PWGLF/FORWARD/analysis2/AliForwardUtil.cxx

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