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