PWGLF/FORWARD/analysis2/AliFMDEnergyFitter.cxx

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