PWGLF/FORWARD/analysis2/AliFMDEnergyFitter.cxx

   1 //
   2 // Class to do the energy correction of FMD ESD data
   3 //
   4 #include "AliFMDEnergyFitter.h"
   5 #include "AliForwardCorrectionManager.h"
   6 #include <AliESDFMD.h>
   7 #include <TAxis.h>
   8 #include <TList.h>
   9 #include <TH1.h>
  10 #include <TF1.h>
  11 #include <TMath.h>
  12 #include <AliLog.h>
  13 #include <TClonesArray.h>
  14 #include <TFitResult.h>
  15 #include <THStack.h>
  16 #include <TROOT.h>
  17 #include <iostream>
  18 #include <iomanip>
  19
  20 ClassImp(AliFMDEnergyFitter)
  21 #if 0
  22 ; // This is for Emacs
  23 #endif
  24 namespace {
  25   const char* fgkEDistFormat = "%s_etabin%03d";
  26 }
  27
  28
  29 //____________________________________________________________________
  30 AliFMDEnergyFitter::AliFMDEnergyFitter()
  31   : TNamed(),
  32     fRingHistos(),
  33     fLowCut(0.4),
  34     fNParticles(5),
  35     fMinEntries(1000),
  36     fFitRangeBinWidth(4),
  37     fDoFits(false),
  38     fDoMakeObject(false),
  39     fEtaAxis(),
  40     fCentralityAxis(),
  41     fMaxE(10),
  42     fNEbins(300),
  43     fUseIncreasingBins(true),
  44     fMaxRelParError(.25),
  45     fMaxChi2PerNDF(20),
  46     fMinWeight(1e-7),
  47     fDebug(0)
  48 {
  49   //
  50   // Default Constructor - do not use
  51   //
  52   DGUARD(fDebug, 0, "Default CTOR of AliFMDEnergyFitter");
  53 }
  54
  55 //____________________________________________________________________
  56 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
  57   : TNamed("fmdEnergyFitter", title),
  58     fRingHistos(),
  59     fLowCut(0.4),
  60     fNParticles(5),
  61     fMinEntries(1000),
  62     fFitRangeBinWidth(4),
  63     fDoFits(false),
  64     fDoMakeObject(false),
  65     fEtaAxis(0,0,0),
  66     fCentralityAxis(0,0,0),
  67     fMaxE(10),
  68     fNEbins(300),
  69     fUseIncreasingBins(true),
  70     fMaxRelParError(.25),
  71     fMaxChi2PerNDF(20),
  72     fMinWeight(1e-7),
  73     fDebug(3)
  74 {
  75   //
  76   // Constructor
  77   //
  78   // Parameters:
  79   //    title Title of object  - not significant
  80   //
  81   DGUARD(fDebug, 0, "Named CTOR of AliFMDEnergyFitter: %s", title);
  82   fEtaAxis.SetName("etaAxis");
  83   fEtaAxis.SetTitle("#eta");
  84   fCentralityAxis.SetName("centralityAxis");
  85   fCentralityAxis.SetTitle("Centrality [%]");
  86   fRingHistos.Add(new RingHistos(1, 'I'));
  87   fRingHistos.Add(new RingHistos(2, 'I'));
  88   fRingHistos.Add(new RingHistos(2, 'O'));
  89   fRingHistos.Add(new RingHistos(3, 'I'));
  90   fRingHistos.Add(new RingHistos(3, 'O'));
  91 }
  92
  93 //____________________________________________________________________
  94 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
  95   : TNamed(o),
  96     fRingHistos(),
  97     fLowCut(o.fLowCut),
  98     fNParticles(o.fNParticles),
  99     fMinEntries(o.fMinEntries),
 100     fFitRangeBinWidth(o.fFitRangeBinWidth),
 101     fDoFits(o.fDoFits),
 102     fDoMakeObject(o.fDoMakeObject),
 103     fEtaAxis(o.fEtaAxis),
 104     fCentralityAxis(o.fCentralityAxis),
 105     fMaxE(o.fMaxE),
 106     fNEbins(o.fNEbins),
 107     fUseIncreasingBins(o.fUseIncreasingBins),
 108     fMaxRelParError(o.fMaxRelParError),
 109     fMaxChi2PerNDF(o.fMaxChi2PerNDF),
 110     fMinWeight(o.fMinWeight),
 111     fDebug(o.fDebug)
 112 {
 113   //
 114   // Copy constructor
 115   //
 116   // Parameters:
 117   //    o Object to copy from
 118   //
 119   DGUARD(fDebug, 0, "Copy CTOR of AliFMDEnergyFitter");
 120   TIter    next(&o.fRingHistos);
 121   TObject* obj = 0;
 122   while ((obj = next())) fRingHistos.Add(obj);
 123 }
 124
 125 //____________________________________________________________________
 126 AliFMDEnergyFitter::~AliFMDEnergyFitter()
 127 {
 128   //
 129   // Destructor
 130   //
 131   fRingHistos.Delete();
 132 }
 133
 134 //____________________________________________________________________
 135 AliFMDEnergyFitter&
 136 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
 137 {
 138   //
 139   // Assignment operator
 140   //
 141   // Parameters:
 142   //    o Object to assign from
 143   //
 144   // Return:
 145   //    Reference to this
 146   //
 147   DGUARD(fDebug, 3, "Assignment of AliFMDEnergyFitter");
 148   if (&o == this) return *this;
 149   TNamed::operator=(o);
 150
 151   fLowCut        = o.fLowCut;
 152   fNParticles       = o.fNParticles;
 153   fMinEntries    = o.fMinEntries;
 154   fFitRangeBinWidth= o.fFitRangeBinWidth;
 155   fDoFits        = o.fDoFits;
 156   fDoMakeObject  = o.fDoMakeObject;
 157   fEtaAxis.Set(o.fEtaAxis.GetNbins(),
 158                o.fEtaAxis.GetXmin(),
 159                o.fEtaAxis.GetXmax());
 160   if (o.fCentralityAxis.GetXbins()) {
 161     const TArrayD* bins = o.fCentralityAxis.GetXbins();
 162     fCentralityAxis.Set(bins->GetSize()-1,bins->GetArray());
 163   }
 164   else
 165     fCentralityAxis.Set(o.fCentralityAxis.GetNbins(),
 166                         o.fCentralityAxis.GetXmin(),
 167                         o.fCentralityAxis.GetXmax());
 168   fDebug         = o.fDebug;
 169   fMaxRelParError= o.fMaxRelParError;
 170   fMaxChi2PerNDF = o.fMaxChi2PerNDF;
 171   fMinWeight     = o.fMinWeight;
 172
 173   fRingHistos.Delete();
 174   TIter    next(&o.fRingHistos);
 175   TObject* obj = 0;
 176   while ((obj = next())) fRingHistos.Add(obj);
 177
 178   return *this;
 179 }
 180
 181 //____________________________________________________________________
 182 AliFMDEnergyFitter::RingHistos*
 183 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
 184 {
 185   //
 186   // Get the ring histogram container
 187   //
 188   // Parameters:
 189   //    d Detector
 190   //    r Ring
 191   //
 192   // Return:
 193   //    Ring histogram container
 194   //
 195   Int_t idx = -1;
 196   switch (d) {
 197   case 1: idx = 0; break;
 198   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 199   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 200   }
 201   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 202
 203   return static_cast<RingHistos*>(fRingHistos.At(idx));
 204 }
 205
 206 //____________________________________________________________________
 207 void
 208 AliFMDEnergyFitter::Init(const TAxis& eAxis)
 209 {
 210   //
 211   // Initialise the task
 212   //
 213   // Parameters:
 214   //    etaAxis The eta axis to use.  Note, that if the eta axis
 215   // has already been set (using SetEtaAxis), then this parameter will be
 216   // ignored
 217   //
 218   DGUARD(fDebug, 1, "Initialize of AliFMDEnergyFitter");
 219   if (fEtaAxis.GetNbins() == 0 ||
 220       TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7)
 221     SetEtaAxis(eAxis);
 222   if (fCentralityAxis.GetNbins() == 0) {
 223     UShort_t n = 12;
 224     Double_t bins[] = {  0.,  5., 10., 15., 20., 30.,
 225                          40., 50., 60., 70., 80., 100. };
 226     SetCentralityAxis(n, bins);
 227   }
 228   TIter    next(&fRingHistos);
 229   RingHistos* o = 0;
 230   while ((o = static_cast<RingHistos*>(next())))
 231     o->Init(fEtaAxis, fCentralityAxis, fMaxE, fNEbins, fUseIncreasingBins);
 232 }
 233 //____________________________________________________________________
 234 void
 235 AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
 236 {
 237   //
 238   // Set the eta axis to use.  This will force the code to use this
 239   // eta axis definition - irrespective of whatever axis is passed to
 240   // the Init member function.  Therefore, this member function can be
 241   // used to force another eta axis than one found in the correction
 242   // objects.
 243   //
 244   // Parameters:
 245   //    etaAxis Eta axis to use
 246   //
 247   SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
 248 }
 249 //____________________________________________________________________
 250 void
 251 AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
 252 {
 253   //
 254   // Set the eta axis to use.  This will force the code to use this
 255   // eta axis definition - irrespective of whatever axis is passed to
 256   // the Init member function.  Therefore, this member function can be
 257   // used to force another eta axis than one found in the correction
 258   // objects.
 259   //
 260   // Parameters:
 261   //    nBins  Number of bins
 262   //    etaMin Minimum of the eta axis
 263   //    etaMax Maximum of the eta axis
 264   //
 265   fEtaAxis.Set(nBins,etaMin,etaMax);
 266 }
 267 //____________________________________________________________________
 268 void
 269 AliFMDEnergyFitter::SetCentralityAxis(UShort_t n, Double_t* bins)
 270 {
 271   fCentralityAxis.Set(n-1, bins);
 272 }
 273
 274
 275 //____________________________________________________________________
 276 Bool_t
 277 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
 278                                Double_t         cent,
 279                                Bool_t           empty)
 280 {
 281   //
 282   // Fitter the input AliESDFMD object
 283   //
 284   // Parameters:
 285   //    input     Input
 286   //    cent      Centrality
 287   //    empty     Whether the event is 'empty'
 288   //
 289   // Return:
 290   //    True on success, false otherwise
 291   //
 292   DGUARD(fDebug, 1, "Accumulate statistics in AliFMDEnergyFitter");
 293   Int_t icent = fCentralityAxis.FindBin(cent);
 294   if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0;
 295
 296   for(UShort_t d = 1; d <= 3; d++) {
 297     Int_t nRings = (d == 1 ? 1 : 2);
 298     for (UShort_t q = 0; q < nRings; q++) {
 299       Char_t      r    = (q == 0 ? 'I' : 'O');
 300       UShort_t    nsec = (q == 0 ?  20 :  40);
 301       UShort_t    nstr = (q == 0 ? 512 : 256);
 302
 303       RingHistos* histos = GetRingHistos(d, r);
 304
 305       for(UShort_t s = 0; s < nsec;  s++) {
 306         for(UShort_t t = 0; t < nstr; t++) {
 307           Float_t mult = input.Multiplicity(d,r,s,t);
 308
 309           // Keep dead-channel information.
 310           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0)
 311             continue;
 312
 313           // Get the pseudo-rapidity
 314           Double_t eta1 = input.Eta(d,r,s,t);
 315           Int_t ieta = fEtaAxis.FindBin(eta1);
 316           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue;
 317
 318           histos->Fill(empty, ieta-1, icent, mult);
 319
 320         } // for strip
 321       } // for sector
 322     } // for ring
 323   } // for detector
 324
 325   return kTRUE;
 326 }
 327
 328 //____________________________________________________________________
 329 void
 330 AliFMDEnergyFitter::Fit(const TList* dir)
 331 {
 332   //
 333   // Scale the histograms to the total number of events
 334   //
 335   // Parameters:
 336   //    dir Where the histograms are
 337   //
 338   DGUARD(fDebug, 1, "Fit distributions in AliFMDEnergyFitter");
 339   if (!fDoFits) return;
 340
 341   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 342   if (!d) return;
 343
 344   // +1          for chi^2
 345   // +3          for 1 landau
 346   // +1          for N
 347   // +fNParticles-1 for weights
 348   Int_t nStack = kN+fNParticles;
 349   THStack* stack[nStack];
 350   stack[0]         = new THStack("chi2",   "#chi^{2}/#nu");
 351   stack[kC     +1] = new THStack("c",      "Constant");
 352   stack[kDelta +1] = new THStack("delta",  "#Delta_{p}");
 353   stack[kXi    +1] = new THStack("xi",     "#xi");
 354   stack[kSigma +1] = new THStack("sigma",  "#sigma");
 355   stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}");
 356   stack[kN     +1] = new THStack("n",      "# Particles");
 357   for (Int_t i = 2; i <= fNParticles; i++)
 358     stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i));
 359   for (Int_t i = 0; i < nStack; i++)
 360     d->Add(stack[i]);
 361
 362   TIter    next(&fRingHistos);
 363   RingHistos* o = 0;
 364   while ((o = static_cast<RingHistos*>(next()))) {
 365     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
 366                           fMinEntries, fFitRangeBinWidth,
 367                           fMaxRelParError, fMaxChi2PerNDF);
 368     if (!l) continue;
 369     for (Int_t i = 0; i < l->GetEntriesFast(); i++) {
 370       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
 371     }
 372   }
 373
 374   if (!fDoMakeObject) return;
 375
 376   MakeCorrectionsObject(d);
 377 }
 378
 379 //____________________________________________________________________
 380 void
 381 AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
 382 {
 383   //
 384   // Generate the corrections object
 385   //
 386   // Parameters:
 387   //    dir List to analyse
 388   //
 389   DGUARD(fDebug, 1, "Make the correction objec in AliFMDEnergyFitter");
 390   AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
 391
 392   AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
 393   obj->SetEtaAxis(fEtaAxis);
 394   obj->SetLowCut(fLowCut);
 395
 396   TIter    next(&fRingHistos);
 397   RingHistos* o = 0;
 398   while ((o = static_cast<RingHistos*>(next()))) {
 399     o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError,
 400                     fMaxChi2PerNDF, fMinWeight);
 401   }
 402
 403   TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
 404   TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
 405   AliInfo(Form("Object %s created in output - should be extracted and copied "
 406                "to %s", oName.Data(), fileName.Data()));
 407   d->Add(new TNamed("filename", fileName.Data()));
 408   d->Add(obj, oName.Data());
 409 }
 410
 411 //____________________________________________________________________
 412 void
 413 AliFMDEnergyFitter::DefineOutput(TList* dir)
 414 {
 415   //
 416   // Define the output histograms.  These are put in a sub list of the
 417   // passed list.   The histograms are merged before the parent task calls
 418   // AliAnalysisTaskSE::Terminate
 419   //
 420   // Parameters:
 421   //    dir Directory to add to
 422   //
 423   DGUARD(fDebug, 1, "Define output in AliFMDEnergyFitter");
 424   TList* d = new TList;
 425   d->SetName(GetName());
 426   dir->Add(d);
 427   d->Add(&fEtaAxis);
 428   TIter    next(&fRingHistos);
 429   RingHistos* o = 0;
 430   while ((o = static_cast<RingHistos*>(next()))) {
 431     o->Output(d);
 432   }
 433 }
 434 //____________________________________________________________________
 435 void
 436 AliFMDEnergyFitter::SetDebug(Int_t dbg)
 437 {
 438   //
 439   // Set the debug level.  The higher the value the more output
 440   //
 441   // Parameters:
 442   //    dbg Debug level
 443   //
 444   fDebug = dbg;
 445   TIter    next(&fRingHistos);
 446   RingHistos* o = 0;
 447   while ((o = static_cast<RingHistos*>(next())))
 448     o->fDebug = dbg;
 449 }
 450
 451 //____________________________________________________________________
 452 void
 453 AliFMDEnergyFitter::Print(Option_t*) const
 454 {
 455   //
 456   // Print information
 457   //
 458   // Parameters:
 459   //    option Not used
 460   //
 461   char ind[gROOT->GetDirLevel()+1];
 462   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 463   ind[gROOT->GetDirLevel()] = '\0';
 464   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 465             << ind << " Low cut:                " << fLowCut << " E/E_mip\n"
 466             << ind << " Max(particles):         " << fNParticles << '\n'
 467             << ind << " Min(entries):           " << fMinEntries << '\n'
 468             << ind << " Fit range:              "
 469             << fFitRangeBinWidth << " bins\n"
 470             << ind << " Make fits:              "
 471             << (fDoFits ? "yes\n" : "no\n")
 472             << ind << " Make object:            "
 473             << (fDoMakeObject ? "yes\n" : "no\n")
 474             << ind << " Max E/E_mip:            " << fMaxE << '\n'
 475             << ind << " N bins:                 " << fNEbins << '\n'
 476             << ind << " Increasing bins:        "
 477             << (fUseIncreasingBins ?"yes\n" : "no\n")
 478             << ind << " max(delta p/p):         " << fMaxRelParError << '\n'
 479             << ind << " max(chi^2/nu):          " << fMaxChi2PerNDF << '\n'
 480             << ind << " min(a_i):               " << fMinWeight << std::endl;
 481 }
 482
 483 //====================================================================
 484 AliFMDEnergyFitter::RingHistos::RingHistos()
 485   : AliForwardUtil::RingHistos(),
 486     fEDist(0),
 487     fEmpty(0),
 488     fEtaEDists(),
 489     fList(0),
 490     fFits("AliFMDCorrELossFit::ELossFit"),
 491     fDebug(0)
 492 {
 493   //
 494   // Default CTOR
 495   //
 496   DGUARD(fDebug, 0, "Default CTOR AliFMDEnergyFitter::RingHistos");
 497 }
 498
 499 //____________________________________________________________________
 500 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
 501   : AliForwardUtil::RingHistos(d,r),
 502     fEDist(0),
 503     fEmpty(0),
 504     fEtaEDists(),
 505     fList(0),
 506     fFits("AliFMDCorrELossFit::ELossFit"),
 507     fDebug(0)
 508 {
 509   //
 510   // Constructor
 511   //
 512   // Parameters:
 513   //    d detector
 514   //    r ring
 515   //
 516   fEtaEDists.SetName("EDists");
 517   DGUARD(fDebug, 0, "Named CTOR AliFMDEnergyFitter::RingHistos: FMD%d%c",
 518          d, r);
 519 }
 520 //____________________________________________________________________
 521 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
 522   : AliForwardUtil::RingHistos(o),
 523     fEDist(o.fEDist),
 524     fEmpty(o.fEmpty),
 525     fEtaEDists(),
 526     fList(0),
 527     fFits("AliFMDCorrELossFit::ELossFit"),
 528     fDebug(0)
 529 {
 530   //
 531   // Copy constructor
 532   //
 533   // Parameters:
 534   //    o Object to copy from
 535   //
 536   DGUARD(fDebug, 0, "Copy CTOR AliFMDEnergyFitter::RingHistos");
 537   fFits.Clear();
 538   TIter next(&o.fEtaEDists);
 539   TObject* obj = 0;
 540   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 541   if (o.fList) {
 542     fList = new TList;
 543     fList->SetName(fName);
 544     TIter next2(o.fList);
 545     while ((obj = next2())) fList->Add(obj->Clone());
 546   }
 547 }
 548
 549 //____________________________________________________________________
 550 AliFMDEnergyFitter::RingHistos&
 551 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
 552 {
 553   //
 554   // Assignment operator
 555   //
 556   // Parameters:
 557   //    o Object to assign from
 558   //
 559   // Return:
 560   //    Reference to this
 561   //
 562   DGUARD(fDebug, 3, "Assignment of AliFMDEnergyFitter::RingHistos");
 563   if (&o == this) return *this;
 564   AliForwardUtil::RingHistos::operator=(o);
 565
 566   if (fEDist) delete fEDist;
 567   if (fEmpty) delete fEmpty;
 568   fEtaEDists.Delete();
 569   if (fList) fList->Delete();
 570   fFits.Clear();
 571
 572   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
 573   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
 574
 575   TIter next(&o.fEtaEDists);
 576   TObject* obj = 0;
 577   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 578
 579   if (o.fList) {
 580     fList = new TList;
 581     fList->SetName(fName);
 582     TIter next2(o.fList);
 583     while ((obj = next2())) fList->Add(obj->Clone());
 584   }
 585
 586   return *this;
 587 }
 588 //____________________________________________________________________
 589 AliFMDEnergyFitter::RingHistos::~RingHistos()
 590 {
 591   //
 592   // Destructor
 593   //
 594   DGUARD(fDebug, 3, "DTOR of AliFMDEnergyFitter::RingHistos");
 595   if (fEDist) delete fEDist;
 596   fEtaEDists.Delete();
 597 }
 598
 599 //____________________________________________________________________
 600 void
 601 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta,
 602                                      Int_t /* icent */,  Double_t mult)
 603 {
 604   //
 605   // Fill histogram
 606   //
 607   // Parameters:
 608   //    empty  True if event is empty
 609   //    ieta   Eta bin (0-based)
 610   //    icent  Centrality bin (1-based)
 611   //    mult   Signal
 612   //
 613   DGUARD(fDebug, 2, "Filling in AliFMDEnergyFitter::RingHistos");
 614   if (empty) {
 615     fEmpty->Fill(mult);
 616     return;
 617   }
 618   fEDist->Fill(mult);
 619
 620   // Here, we should first look up in a centrality array, and then in
 621   // that array look up the eta bin - or perhaps we should do it the
 622   // other way around?
 623   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
 624
 625   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
 626   if (!h) return;
 627
 628   h->Fill(mult);
 629 }
 630
 631 //__________________________________________________________________
 632 TArrayD
 633 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
 634                                                    Double_t max) const
 635 {
 636   //
 637   // Make an axis with increasing bins
 638   //
 639   // Parameters:
 640   //    n    Number of bins
 641   //    min  Minimum
 642   //    max  Maximum
 643   //
 644   // Return:
 645   //    An axis with quadratically increasing bin size
 646   //
 647
 648   TArrayD bins(n+1);
 649   Double_t dx = (max-min) / n;
 650   bins[0]     = min;
 651   Int_t    i  = 1;
 652   for (i = 1; i < n+1; i++) {
 653     Double_t dI   = float(i)/n;
 654     Double_t next = bins[i-1] + dx + dI * dI * dx;
 655     bins[i]       = next;
 656     if (next > max) break;
 657   }
 658   bins.Set(i+1);
 659   return bins;
 660 }
 661
 662 //____________________________________________________________________
 663 void
 664 AliFMDEnergyFitter::RingHistos::Make(Int_t    ieta,
 665                                      Double_t emin,
 666                                      Double_t emax,
 667                                      Double_t deMax,
 668                                      Int_t    nDeBins,
 669                                      Bool_t   incr)
 670 {
 671   //
 672   // Make E/E_mip histogram
 673   //
 674   // Parameters:
 675   //    ieta    Eta bin
 676   //    eMin    Least signal
 677   //    eMax    Largest signal
 678   //    deMax   Maximum energy loss
 679   //    nDeBins Number energy loss bins
 680   //    incr    Whether to make bins of increasing size
 681   //
 682   TH1D* h = 0;
 683   TString name  = Form(fgkEDistFormat, fName.Data(), ieta);
 684   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
 685                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
 686   if (!incr)
 687     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
 688   else {
 689     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
 690     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
 691   }
 692
 693   h->SetDirectory(0);
 694   h->SetXTitle("#DeltaE/#DeltaE_{mip}");
 695   h->SetFillColor(Color());
 696   h->SetMarkerColor(Color());
 697   h->SetLineColor(Color());
 698   h->SetFillStyle(3001);
 699   h->Sumw2();
 700
 701   fEtaEDists.AddAt(h, ieta-1);
 702 }
 703 //____________________________________________________________________
 704 TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
 705                                               const char* title,
 706                                               const TAxis& eta) const
 707 {
 708   //
 709   // Make a parameter histogram
 710   //
 711   // Parameters:
 712   //    name   Name of histogram.
 713   //    title  Title of histogram.
 714   //    eta    Eta axis
 715   //
 716   // Return:
 717   //
 718   //
 719   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 720                      Form("%s for %s", title, fName.Data()),
 721                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
 722   h->SetXTitle("#eta");
 723   h->SetYTitle(title);
 724   h->SetDirectory(0);
 725   h->SetFillColor(Color());
 726   h->SetMarkerColor(Color());
 727   h->SetLineColor(Color());
 728   h->SetFillStyle(3001);
 729   h->Sumw2();
 730
 731   return h;
 732 }
 733 //____________________________________________________________________
 734 TH1D*
 735 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
 736                                           const char* title,
 737                                           const TAxis& eta,
 738                                           Int_t low,
 739                                           Int_t high,
 740                                           Double_t val,
 741                                           Double_t err) const
 742 {
 743   //
 744   // Make a histogram that contains the results of the fit over the full ring
 745   //
 746   // Parameters:
 747   //    name  Name
 748   //    title Title
 749   //    eta   Eta axis
 750   //    low   Least bin
 751   //    high  Largest bin
 752   //    val   Value of parameter
 753   //    err   Error on parameter
 754   //
 755   // Return:
 756   //    The newly allocated histogram
 757   //
 758   Double_t xlow  = eta.GetBinLowEdge(low);
 759   Double_t xhigh = eta.GetBinUpEdge(high);
 760   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 761                      Form("%s for %s", title, fName.Data()),
 762                      1, xlow, xhigh);
 763   h->SetBinContent(1, val);
 764   h->SetBinError(1, err);
 765   h->SetXTitle("#eta");
 766   h->SetYTitle(title);
 767   h->SetDirectory(0);
 768   h->SetFillColor(0);
 769   h->SetMarkerColor(Color());
 770   h->SetLineColor(Color());
 771   h->SetFillStyle(0);
 772   h->SetLineStyle(2);
 773   h->SetLineWidth(2);
 774
 775   return h;
 776 }
 777
 778 //____________________________________________________________________
 779 void
 780 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
 781                                      const TAxis& /* cAxis */,
 782                                      Double_t maxDE, Int_t nDEbins,
 783                                      Bool_t useIncrBin)
 784 {
 785   //
 786   // Initialise object
 787   //
 788   // Parameters:
 789   //    eAxis      Eta axis
 790   //    maxDE      Max energy loss to consider
 791   //    nDEbins    Number of bins
 792   //    useIncrBin Whether to use an increasing bin size
 793   //
 794   DGUARD(fDebug, 2, "Initialize in AliFMDEnergyFitter::RingHistos");
 795   fEDist = new TH1D(Form("%s_edist", fName.Data()),
 796                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()),
 797                     200, 0, 6);
 798   fEmpty = new TH1D(Form("%s_empty", fName.Data()),
 799                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)",
 800                          fName.Data()), 200, 0, 6);
 801   fList->Add(fEDist);
 802   fList->Add(fEmpty);
 803   // Here, we should make an array of centrality bins ...
 804   // In that case, the structure will be
 805   //
 806   //    -+- Ring1 -+- Centrality1 -+- Eta1
 807   //     |         |               +- Eta2
 808   //     ...       ...             ...
 809   //     |         |               +- EtaM
 810   //     |         +- Centrality2 -+- Eta1
 811   //     ...       ...             ...
 812   //     |         |               +- EtaM
 813   //     ...       ...
 814   //     |         +- CentralityN -+- Eta1
 815   //     ...       ...             ...
 816   //     |                         +- EtaM
 817   //    -+- Ring2 -+- Centrality1 -+- Eta1
 818   //     |         |               +- Eta2
 819   //     ...       ...             ...
 820   //     |         |               +- EtaM
 821   //     |         +- Centrality2 -+- Eta1
 822   //     ...       ...             ...
 823   //     |         |               +- EtaM
 824   //     ...       ...
 825   //     |         +- CentralityN -+- Eta1
 826   //     ...       ...             ...
 827   //     |                         +- EtaM
 828   //     ...       ...             ...
 829   //    -+- RingO -+- Centrality1 -+- Eta1
 830   //               |               +- Eta2
 831   //               ...             ...
 832   //               |               +- EtaM
 833   //               +- Centrality2 -+- Eta1
 834   //               ...             ...
 835   //               |               +- EtaM
 836   //               ...
 837   //               +- CentralityN -+- Eta1
 838   //               ...             ...
 839   //                               +- EtaM
 840   //
 841   // fEtaEDists.Expand(eAxis.GetNbins());
 842   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) {
 843     Double_t min = eAxis.GetBinLowEdge(i);
 844     Double_t max = eAxis.GetBinUpEdge(i);
 845     // Or may we should do it here?  In that case we'd have the
 846     // following structure:
 847     //     Ring1 -+- Eta1 -+- Centrality1
 848     //            |        +- Centrality2
 849     //            ...      ...
 850     //            |        +- CentralityN
 851     //            +- Eta2 -+- Centrality1
 852     //            |        +- Centrality2
 853     //            ...      ...
 854     //            |        +- CentralityN
 855     //            ...
 856     //            +- EtaM -+- Centrality1
 857     //                     +- Centrality2
 858     //                     ...
 859     //                     +- CentralityN
 860     Make(i, min, max, maxDE, nDEbins, useIncrBin);
 861   }
 862   fList->Add(&fEtaEDists);
 863   // fEtaEDists.ls();
 864 }
 865 //____________________________________________________________________
 866 TObjArray*
 867 AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
 868                                     const TAxis&     eta,
 869                                     Double_t         lowCut,
 870                                     UShort_t         nParticles,
 871                                     UShort_t         minEntries,
 872                                     UShort_t         minusBins,
 873                                     Double_t         relErrorCut,
 874                                     Double_t         chi2nuCut) const
 875 {
 876   //
 877   // Fit each histogram to up to @a nParticles particle responses.
 878   //
 879   // Parameters:
 880   //    dir         Output list
 881   //    eta         Eta axis
 882   //    lowCut      Lower cut
 883   //    nParticles  Max number of convolved landaus to fit
 884   //    minEntries  Minimum number of entries
 885   //    minusBins   Number of bins from peak to subtract to
 886   //                    get the fit range
 887   //    relErrorCut Cut applied to relative error of parameter.
 888   //                    Note, for multi-particle weights, the cut
 889   //                    is loosend by a factor of 2
 890   //    chi2nuCut   Cut on @f$ \chi^2/\nu@f$ -
 891   //                    the reduced @f$\chi^2@f$
 892   //
 893   DGUARD(fDebug, 2, "Fit in AliFMDEnergyFitter::RingHistos");
 894
 895   // Get our ring list from the output container
 896   TList* l = GetOutputList(dir);
 897   if (!l) return 0;
 898
 899   // Get the energy distributions from the output container
 900   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
 901   if (!dists) {
 902     AliWarning(Form("Didn't find %s_EtaEDists in %s",
 903                     fName.Data(), l->GetName()));
 904     l->ls();
 905     return 0;
 906   }
 907
 908   // Container of the fit results histograms
 909   TObjArray* pars  = new TObjArray(3+nParticles+1);
 910   pars->SetName("FitResults");
 911   l->Add(pars);
 912
 913   // Result objects stored in separate list on the output
 914   TH1* hChi2   = 0;
 915   TH1* hC      = 0;
 916   TH1* hDelta  = 0;
 917   TH1* hXi     = 0;
 918   TH1* hSigma  = 0;
 919   TH1* hSigmaN = 0;
 920   TH1* hN      = 0;
 921   TH1* hA[nParticles-1];
 922   pars->Add(hChi2   = MakePar("chi2",   "#chi^{2}/#nu", eta));
 923   pars->Add(hC      = MakePar("c",      "Constant",     eta));
 924   pars->Add(hDelta  = MakePar("delta",  "#Delta_{p}",   eta));
 925   pars->Add(hXi     = MakePar("xi",     "#xi",          eta));
 926   pars->Add(hSigma  = MakePar("sigma",  "#sigma",       eta));
 927   pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}",   eta));
 928   pars->Add(hN      = MakePar("n",      "N", eta));
 929   for (UShort_t i = 1; i < nParticles; i++)
 930     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
 931
 932
 933   Int_t nDists = dists->GetEntries();
 934   Int_t low    = nDists;
 935   Int_t high   = 0;
 936   Int_t nEmpty = 0;
 937   Int_t nLow   = 0;
 938   Int_t nFitted= 0;
 939   for (Int_t i = 0; i < nDists; i++) {
 940     TH1D* dist = static_cast<TH1D*>(dists->At(i));
 941     // Ignore empty histograms altoghether
 942     if (!dist || dist->GetEntries() <= 0) {
 943       nEmpty++;
 944       continue;
 945     }
 946
 947     // Scale to the bin-width
 948     dist->Scale(1., "width");
 949
 950     // Normalize peak to 1
 951     Double_t max = dist->GetMaximum();
 952     if (max <= 0) continue;
 953     dist->Scale(1/max);
 954
 955     // Check that we have enough entries
 956     Int_t nEntries = Int_t(dist->GetEntries());
 957     if (nEntries <= minEntries) {
 958       AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
 959                       i, dist->GetName(), nEntries, minEntries));
 960       nLow++;
 961       continue;
 962     }
 963
 964     // Now fit
 965     TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
 966                        relErrorCut,chi2nuCut);
 967     if (!res) continue;
 968     nFitted++;
 969     // dist->GetListOfFunctions()->Add(res);
 970
 971     // Store eta limits
 972     low   = TMath::Min(low,i+1);
 973     high  = TMath::Max(high,i+1);
 974
 975     // Get the reduced chi square
 976     Double_t chi2 = res->GetChisquare();
 977     Int_t    ndf  = res->GetNDF();
 978
 979     // Store results of best fit in output histograms
 980     res->SetLineWidth(4);
 981     hChi2   ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
 982     hC      ->SetBinContent(i+1, res->GetParameter(kC));
 983     hDelta  ->SetBinContent(i+1, res->GetParameter(kDelta));
 984     hXi     ->SetBinContent(i+1, res->GetParameter(kXi));
 985     hSigma  ->SetBinContent(i+1, res->GetParameter(kSigma));
 986     hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));
 987     hN      ->SetBinContent(i+1, res->GetParameter(kN));
 988
 989     hC     ->SetBinError(i+1, res->GetParError(kC));
 990     hDelta ->SetBinError(i+1, res->GetParError(kDelta));
 991     hXi    ->SetBinError(i+1, res->GetParError(kXi));
 992     hSigma ->SetBinError(i+1, res->GetParError(kSigma));
 993     hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
 994     hN     ->SetBinError(i+1, res->GetParError(kN));
 995
 996     for (UShort_t j = 0; j < nParticles-1; j++) {
 997       hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
 998       hA[j]->SetBinError(i+1, res->GetParError(kA+j));
 999     }
1000   }
1001   printf("%s: Out of %d histograms, %d where empty, %d had too little data,"
1002          "leaving %d to be fitted, of which %d succeeded\n",
1003          GetName(), nDists, nEmpty, nLow, nDists-nEmpty-nLow, nFitted);
1004
1005   TH1* status = new TH1I("status", "Status of Fits", 5, 0, 5);
1006   status->GetXaxis()->SetBinLabel(1, "Total");
1007   status->GetXaxis()->SetBinLabel(2, "Empty");
1008   status->GetXaxis()->SetBinLabel(3, Form("<%d", minEntries));
1009   status->GetXaxis()->SetBinLabel(4, "Candidates");
1010   status->GetXaxis()->SetBinLabel(5, "Fitted");
1011   status->SetXTitle("Status");
1012   status->SetYTitle("# of #Delta distributions");
1013   status->SetBinContent(1, nDists);
1014   status->SetBinContent(2, nEmpty);
1015   status->SetBinContent(3, nLow);
1016   status->SetBinContent(4, nDists-nLow-nEmpty);
1017   status->SetBinContent(5, nFitted);
1018   status->SetFillColor(Color());
1019   status->SetFillStyle(3001);
1020   status->SetLineColor(Color());
1021   status->SetDirectory(0);
1022   status->SetStats(0);
1023   pars->Add(status);
1024
1025   // Fit the full-ring histogram
1026   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
1027   if (total && total->GetEntries() >= minEntries) {
1028
1029     // Scale to the bin-width
1030     total->Scale(1., "width");
1031
1032     // Normalize peak to 1
1033     Double_t max = total->GetMaximum();
1034     if (max > 0) total->Scale(1/max);
1035
1036     TF1* res = FitHist(total,lowCut,nParticles,minusBins,
1037                        relErrorCut,chi2nuCut);
1038     if (res) {
1039       // Make histograms for the result of this fit
1040       Double_t chi2 = res->GetChisquare();
1041       Int_t    ndf  = res->GetNDF();
1042       pars->Add(MakeTotal("t_chi2",   "#chi^{2}/#nu", eta, low, high,
1043                           ndf > 0 ? chi2/ndf : 0, 0));
1044       pars->Add(MakeTotal("t_c",      "Constant",     eta, low, high,
1045                           res->GetParameter(kC),
1046                           res->GetParError(kC)));
1047       pars->Add(MakeTotal("t_delta",  "#Delta_{p}",   eta, low, high,
1048                           res->GetParameter(kDelta),
1049                           res->GetParError(kDelta)));
1050       pars->Add(MakeTotal("t_xi",     "#xi",          eta, low, high,
1051                           res->GetParameter(kXi),
1052                           res->GetParError(kXi)));
1053       pars->Add(MakeTotal("t_sigma",  "#sigma",       eta, low, high,
1054                           res->GetParameter(kSigma),
1055                           res->GetParError(kSigma)));
1056       pars->Add(MakeTotal("t_sigman", "#sigma_{n}",   eta, low, high,
1057                           res->GetParameter(kSigmaN),
1058                           res->GetParError(kSigmaN)));
1059       pars->Add(MakeTotal("t_n",    "N",              eta, low, high,
1060                           res->GetParameter(kN),0));
1061       for (UShort_t j = 0; j < nParticles-1; j++)
1062         pars->Add(MakeTotal(Form("t_a%d",j+2),
1063                             Form("a_{%d}",j+2), eta, low, high,
1064                             res->GetParameter(kA+j),
1065                             res->GetParError(kA+j)));
1066     }
1067   }
1068
1069   // Clean up list of histogram.  Histograms with no entries or
1070   // no functions are deleted.  We have to do this using the TObjLink
1071   // objects stored in the list since ROOT cannot guaranty the validity
1072   // of iterators when removing from a list - tsck.  Should just implement
1073   // TIter::Remove().
1074   TObjLink* lnk = dists->FirstLink();
1075   while (lnk) {
1076     TH1* h = static_cast<TH1*>(lnk->GetObject());
1077     bool remove = false;
1078     if (h->GetEntries() <= 0) {
1079       // AliWarning(Form("No entries in %s - removing", h->GetName()));
1080       remove = true;
1081     }
1082     else if (h->GetListOfFunctions()->GetEntries() <= 0) {
1083       // AliWarning(Form("No fuctions associated with %s - removing",
1084       //            h->GetName()));
1085       remove = true;
1086     }
1087     if (remove) {
1088       TObjLink* keep = lnk->Next();
1089       dists->Remove(lnk);
1090       lnk = keep;
1091       continue;
1092     }
1093     lnk = lnk->Next();
1094   }
1095
1096   return pars;
1097 }
1098
1099 //____________________________________________________________________
1100 TF1*
1101 AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
1102                                         Double_t lowCut,
1103                                         UShort_t nParticles,
1104                                         UShort_t minusBins,
1105                                         Double_t relErrorCut,
1106                                         Double_t chi2nuCut) const
1107 {
1108   //
1109   // Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
1110   // maximum bin content in the range @f$ [E_{min},\infty]@f$ is
1111   // found.  Then the fit range is set to the bin range
1112   // @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
1113   // particle signal is fitted to that.  The parameters of that fit
1114   // is then used as seeds for a fit of the @f$ N@f$ particle response
1115   // to the data in the range
1116   // @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
1117   //
1118   // Parameters:
1119   //    dist        Histogram to fit
1120   //    lowCut      Lower cut @f$ E_{min}@f$ on signal
1121   //    nParticles  Max number @f$ N@f$ of convolved landaus to fit
1122   //    minusBins   Number of bins @f$ \Delta b@f$ from peak to
1123   //                    subtract to get the fit range
1124   //    relErrorCut Cut applied to relative error of parameter.
1125   //                    Note, for multi-particle weights, the cut
1126   //                    is loosend by a factor of 2
1127   //    chi2nuCut   Cut on @f$ \chi^2/\nu@f$ -
1128   //                    the reduced @f$\chi^2@f$
1129   //
1130   // Return:
1131   //    The best fit function
1132   //
1133   DGUARD(fDebug, 3, "Fit histogram in AliFMDEnergyFitter::RingHistos: %s",
1134          dist->GetName());
1135   Double_t maxRange = 10;
1136
1137   // Create a fitter object
1138   AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins);
1139   f.Clear();
1140
1141
1142   // If we are only asked to fit a single particle, return this fit,
1143   // no matter what.
1144   if (nParticles == 1) {
1145     TF1* r = f.Fit1Particle(dist, 0);
1146     if (!r) return 0;
1147     TF1* ret = new TF1(*r);
1148     dist->GetListOfFunctions()->Add(ret);
1149     return ret;
1150   }
1151
1152   // Fit from 2 upto n particles
1153   for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
1154
1155
1156   // Now, we need to select the best fit
1157   Int_t nFits = f.GetFitResults().GetEntriesFast();
1158   TF1*  good[nFits];
1159   for (Int_t i = nFits-1; i >= 0; i--) {
1160     good[i] = 0;
1161     TF1* ff = static_cast<TF1*>(f.GetFunctions().At(i));
1162     // ff->SetLineColor(Color());
1163     ff->SetRange(0, maxRange);
1164     dist->GetListOfFunctions()->Add(new TF1(*ff));
1165     if (CheckResult(static_cast<TFitResult*>(f.GetFitResults().At(i)),
1166                     relErrorCut, chi2nuCut)) {
1167       good[i] = ff;
1168       ff->SetLineWidth(2);
1169       // f.fFitResults.At(i)->Print("V");
1170     }
1171   }
1172   // If all else fails, use the 1 particle fit
1173   TF1* ret = static_cast<TF1*>(f.GetFunctions().At(0));
1174
1175   // Find the fit with the most valid particles
1176   for (Int_t i = nFits-1; i > 0; i--) {
1177     if (!good[i]) continue;
1178     if (fDebug > 1) {
1179       AliInfo(Form("Choosing fit with n=%d", i+1));
1180       f.GetFitResults().At(i)->Print();
1181     }
1182     ret = good[i];
1183     break;
1184   }
1185   // Give a warning if we're using fall-back
1186   if (ret == f.GetFunctions().At(0)) {
1187     AliWarning("Choosing fall-back 1 particle fit");
1188   }
1189   // Copy our result and return (the functions are owned by the fitter)
1190   TF1* fret = new TF1(*ret);
1191   return fret;
1192 }
1193
1194 //____________________________________________________________________
1195 Bool_t
1196 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
1197                                             Double_t relErrorCut,
1198                                             Double_t chi2nuCut) const
1199 {
1200   //
1201   // Check the result of the fit. Returns true if
1202   // - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
1203   // - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note,
1204   //   for multi-particle fits, this requirement is relaxed by a
1205   //   factor of 2
1206   // - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$
1207   //   particle response
1208   //
1209   // Parameters:
1210   //    r           Result to check
1211   //    relErrorCut Cut @f$ \delta_e@f$ applied to relative error
1212   //                    of parameter.
1213   //    chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$
1214   //
1215   // Return:
1216   //    true if fit is good.
1217   //
1218   if (fDebug > 10) r->Print();
1219   TString  n    = r->GetName();
1220   n.ReplaceAll("TFitResult-", "");
1221   Double_t chi2 = r->Chi2();
1222   Int_t    ndf  = r->Ndf();
1223   // Double_t prob = r.Prob();
1224   Bool_t ret = kTRUE;
1225
1226   // Check that the reduced chi square isn't larger than 5
1227   if (ndf <= 0 || chi2 / ndf > chi2nuCut) {
1228     if (fDebug > 2) {
1229       AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f",
1230                       n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
1231                       chi2nuCut)); }
1232     ret = kFALSE;
1233   }
1234
1235   // Check each parameter
1236   UShort_t nPar = r->NPar();
1237   for (UShort_t i = 0; i < nPar; i++) {
1238     if (i == kN) continue;  // Skip the number parameter
1239
1240     // Get value and error and check value
1241     Double_t v  = r->Parameter(i);
1242     Double_t e  = r->ParError(i);
1243     if (v == 0) continue;
1244
1245     // Calculate the relative error and check it
1246     Double_t re  = e / v;
1247     Double_t cut = relErrorCut * (i < kN ? 1 : 2);
1248     if (re > cut) {
1249       if (fDebug > 2) {
1250         AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
1251                         n.Data(), r->ParName(i).c_str(),
1252                         r->ParName(i).c_str(), e, v,
1253                         100*(v == 0 ? 0 : e/v),
1254                         100*(cut))); }
1255       ret = kFALSE;
1256     }
1257   }
1258
1259   // Check if we have scale parameters
1260   if (nPar > kN) {
1261
1262     // Check that the last particle has a significant contribution
1263     Double_t lastScale = r->Parameter(nPar-1);
1264     if (lastScale <= 1e-7) {
1265       if (fDebug) {
1266         AliWarning(Form("%s: %s=%9.6f<1e-7",
1267                         n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
1268       ret = kFALSE;
1269     }
1270   }
1271   return ret;
1272 }
1273
1274
1275 //__________________________________________________________________
1276 void
1277 AliFMDEnergyFitter::RingHistos::FindBestFits(const TList*        d,
1278                                              AliFMDCorrELossFit& obj,
1279                                              const TAxis&        eta,
1280                                              Double_t            relErrorCut,
1281                                              Double_t            chi2nuCut,
1282                                              Double_t            minWeightCut)
1283 {
1284   //
1285   // Find the best fits
1286   //
1287   // Parameters:
1288   //    d              Parent list
1289   //    obj            Object to add fits to
1290   //    eta            Eta axis
1291   //    relErrorCut    Cut applied to relative error of parameter.
1292   //                       Note, for multi-particle weights, the cut
1293   //                       is loosend by a factor of 2
1294   //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
1295   //                       the reduced @f$\chi^2@f$
1296   //    minWeightCut   Least valid @f$ a_i@f$
1297   //
1298
1299   // Get our ring list from the output container
1300   TList* l = GetOutputList(d);
1301   if (!l) return;
1302
1303   // Get the energy distributions from the output container
1304   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
1305   if (!dists) {
1306     AliWarning(Form("Didn't find %s_EtaEDists in %s",
1307                     fName.Data(), l->GetName()));
1308     l->ls();
1309     return;
1310   }
1311   Int_t nBin = eta.GetNbins();
1312
1313   for (Int_t b = 1; b <= nBin; b++) {
1314     TString n(Form(fgkEDistFormat, fName.Data(), b));
1315     TH1D*   dist = static_cast<TH1D*>(dists->FindObject(n));
1316     // Ignore empty histograms altoghether
1317     if (!dist || dist->GetEntries() <= 0) continue;
1318
1319     AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
1320                                                      relErrorCut,
1321                                                      chi2nuCut,
1322                                                      minWeightCut);
1323     best->fDet  = fDet;
1324     best->fRing = fRing;
1325     best->fBin  = b; //
1326     // Double_t eta = fAxis->GetBinCenter(b);
1327     obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
1328   }
1329 }
1330
1331 //__________________________________________________________________
1332 AliFMDCorrELossFit::ELossFit*
1333 AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist,
1334                                             Double_t relErrorCut,
1335                                             Double_t chi2nuCut,
1336                                             Double_t minWeightCut)
1337 {
1338   //
1339   // Find the best fit
1340   //
1341   // Parameters:
1342   //    dist           Histogram
1343   //    relErrorCut    Cut applied to relative error of parameter.
1344   //                       Note, for multi-particle weights, the cut
1345   //                       is loosend by a factor of 2
1346   //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
1347   //                       the reduced @f$\chi^2@f$
1348   //    minWeightCut   Least valid @f$ a_i@f$
1349   //
1350   // Return:
1351   //    Best fit
1352   //
1353   TList* funcs = dist->GetListOfFunctions();
1354   TIter  next(funcs);
1355   TF1*   func = 0;
1356   fFits.Clear();
1357   Int_t  i = 0;
1358   // Info("FindBestFit", "%s", dist->GetName());
1359   while ((func = static_cast<TF1*>(next()))) {
1360     AliFMDCorrELossFit::ELossFit* fit =
1361       new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
1362     fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
1363   }
1364   fFits.Sort(false);
1365   // fFits.Print();
1366   return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
1367 }
1368
1369
1370 //____________________________________________________________________
1371 void
1372 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
1373 {
1374   //
1375   // Define outputs
1376   //
1377   // Parameters:
1378   //    dir
1379   //
1380   DGUARD(fDebug, 2, "Define output in AliFMDEnergyFitter::RingHistos");
1381   fList = DefineOutputList(dir);
1382 }
1383
1384 //____________________________________________________________________
1385 //
1386 // EOF
1387 //