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