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