PWG2/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 "AliFMDAnaParameters.h"
  12 #include <AliLog.h>
  13 #include <TClonesArray.h>
  14 #include <TFitResult.h>
  15 #include <THStack.h>
  16
  17 ClassImp(AliFMDEnergyFitter)
  18 #if 0
  19 ; // This is for Emacs
  20 #endif
  21
  22
  23 //____________________________________________________________________
  24 AliFMDEnergyFitter::AliFMDEnergyFitter()
  25   : TNamed(),
  26     fRingHistos(),
  27     fLowCut(0.3),
  28     fNParticles(3),
  29     fMinEntries(100),
  30     fFitRangeBinWidth(4),
  31     fDoFits(false),
  32     fEtaAxis(),
  33     fMaxE(10),
  34     fNEbins(300),
  35     fUseIncreasingBins(true),
  36     fMaxRelParError(.25),
  37     fMaxChi2PerNDF(20),
  38     fDebug(0)
  39 {}
  40
  41 //____________________________________________________________________
  42 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
  43   : TNamed("fmdEnergyFitter", title),
  44     fRingHistos(),
  45     fLowCut(0.3),
  46     fNParticles(3),
  47     fMinEntries(100),
  48     fFitRangeBinWidth(4),
  49     fDoFits(false),
  50     fEtaAxis(0,0,0),
  51     fMaxE(10),
  52     fNEbins(300),
  53     fUseIncreasingBins(true),
  54     fMaxRelParError(.25),
  55     fMaxChi2PerNDF(20),
  56     fDebug(3)
  57 {
  58   fEtaAxis.SetName("etaAxis");
  59   fEtaAxis.SetTitle("#eta");
  60   fRingHistos.Add(new RingHistos(1, 'I'));
  61   fRingHistos.Add(new RingHistos(2, 'I'));
  62   fRingHistos.Add(new RingHistos(2, 'O'));
  63   fRingHistos.Add(new RingHistos(3, 'I'));
  64   fRingHistos.Add(new RingHistos(3, 'O'));
  65 }
  66
  67 //____________________________________________________________________
  68 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
  69   : TNamed(o),
  70     fRingHistos(),
  71     fLowCut(o.fLowCut),
  72     fNParticles(o.fNParticles),
  73     fMinEntries(o.fMinEntries),
  74     fFitRangeBinWidth(o.fFitRangeBinWidth),
  75     fDoFits(o.fDoFits),
  76     fEtaAxis(o.fEtaAxis),
  77     fMaxE(o.fMaxE),
  78     fNEbins(o.fNEbins),
  79     fUseIncreasingBins(o.fUseIncreasingBins),
  80     fMaxRelParError(o.fMaxRelParError),
  81     fMaxChi2PerNDF(o.fMaxChi2PerNDF),
  82     fDebug(o.fDebug)
  83 {
  84   TIter    next(&o.fRingHistos);
  85   TObject* obj = 0;
  86   while ((obj = next())) fRingHistos.Add(obj);
  87 }
  88
  89 //____________________________________________________________________
  90 AliFMDEnergyFitter::~AliFMDEnergyFitter()
  91 {
  92   fRingHistos.Delete();
  93 }
  94
  95 //____________________________________________________________________
  96 AliFMDEnergyFitter&
  97 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
  98 {
  99   TNamed::operator=(o);
 100
 101   fLowCut        = o.fLowCut;
 102   fNParticles       = o.fNParticles;
 103   fMinEntries    = o.fMinEntries;
 104   fFitRangeBinWidth= o.fFitRangeBinWidth;
 105   fDoFits        = o.fDoFits;
 106   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
 107   fDebug         = o.fDebug;
 108
 109   fRingHistos.Delete();
 110   TIter    next(&o.fRingHistos);
 111   TObject* obj = 0;
 112   while ((obj = next())) fRingHistos.Add(obj);
 113
 114   return *this;
 115 }
 116
 117 //____________________________________________________________________
 118 AliFMDEnergyFitter::RingHistos*
 119 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
 120 {
 121   Int_t idx = -1;
 122   switch (d) {
 123   case 1: idx = 0; break;
 124   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 125   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 126   }
 127   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 128
 129   return static_cast<RingHistos*>(fRingHistos.At(idx));
 130 }
 131
 132 //____________________________________________________________________
 133 void
 134 AliFMDEnergyFitter::Init(const TAxis& eAxis)
 135 {
 136   if (fEtaAxis.GetNbins() == 0 ||
 137       fEtaAxis.GetXmin() == fEtaAxis.GetXmax())
 138     SetEtaAxis(eAxis);
 139   TIter    next(&fRingHistos);
 140   RingHistos* o = 0;
 141   while ((o = static_cast<RingHistos*>(next())))
 142     o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
 143 }
 144 //____________________________________________________________________
 145 void
 146 AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
 147 {
 148   SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
 149 }
 150 //____________________________________________________________________
 151 void
 152 AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
 153 {
 154   fEtaAxis.Set(nBins,etaMin,etaMax);
 155 }
 156
 157 //____________________________________________________________________
 158 Bool_t
 159 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
 160                                Bool_t           empty)
 161 {
 162   for(UShort_t d = 1; d <= 3; d++) {
 163     Int_t nRings = (d == 1 ? 1 : 2);
 164     for (UShort_t q = 0; q < nRings; q++) {
 165       Char_t      r    = (q == 0 ? 'I' : 'O');
 166       UShort_t    nsec = (q == 0 ?  20 :  40);
 167       UShort_t    nstr = (q == 0 ? 512 : 256);
 168
 169       RingHistos* histos = GetRingHistos(d, r);
 170
 171       for(UShort_t s = 0; s < nsec;  s++) {
 172         for(UShort_t t = 0; t < nstr; t++) {
 173           Float_t mult = input.Multiplicity(d,r,s,t);
 174
 175           // Keep dead-channel information.
 176           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0)
 177             continue;
 178
 179           // Get the pseudo-rapidity
 180           Double_t eta1 = input.Eta(d,r,s,t);
 181           Int_t ieta = fEtaAxis.FindBin(eta1);
 182           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue;
 183
 184           histos->Fill(empty, ieta-1, mult);
 185
 186         } // for strip
 187       } // for sector
 188     } // for ring
 189   } // for detector
 190
 191   return kTRUE;
 192 }
 193
 194 //____________________________________________________________________
 195 void
 196 AliFMDEnergyFitter::Fit(TList* dir)
 197 {
 198   if (!fDoFits) return;
 199
 200   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 201   if (!d) return;
 202
 203   // +1          for chi^2
 204   // +3          for 1 landau
 205   // +1          for N
 206   // +fNParticles-1 for weights
 207   Int_t nStack = kN+fNParticles;
 208   THStack* stack[nStack];
 209   stack[0]         = new THStack("chi2",   "#chi^{2}/#nu");
 210   stack[kC     +1] = new THStack("c",      "Constant");
 211   stack[kDelta +1] = new THStack("delta",  "#Delta_{p}");
 212   stack[kXi    +1] = new THStack("xi",     "#xi");
 213   stack[kSigma +1] = new THStack("sigma",  "#sigma");
 214   stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}");
 215   stack[kN     +1] = new THStack("n",      "# Particles");
 216   for (Int_t i = 2; i <= fNParticles; i++)
 217     stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i));
 218   for (Int_t i = 0; i < nStack; i++)
 219     d->Add(stack[i]);
 220
 221   TIter    next(&fRingHistos);
 222   RingHistos* o = 0;
 223   while ((o = static_cast<RingHistos*>(next()))) {
 224     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
 225                           fMinEntries, fFitRangeBinWidth,
 226                           fMaxRelParError, fMaxChi2PerNDF);
 227     if (!l) continue;
 228     for (Int_t i = 0; i < l->GetEntriesFast(); i++) {
 229       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
 230     }
 231   }
 232 }
 233
 234 //____________________________________________________________________
 235 void
 236 AliFMDEnergyFitter::DefineOutput(TList* dir)
 237 {
 238   TList* d = new TList;
 239   d->SetName(GetName());
 240   dir->Add(d);
 241   d->Add(&fEtaAxis);
 242   TIter    next(&fRingHistos);
 243   RingHistos* o = 0;
 244   while ((o = static_cast<RingHistos*>(next()))) {
 245     o->Output(d);
 246   }
 247 }
 248 //____________________________________________________________________
 249 void
 250 AliFMDEnergyFitter::SetDebug(Int_t dbg)
 251 {
 252   fDebug = dbg;
 253   TIter    next(&fRingHistos);
 254   RingHistos* o = 0;
 255   while ((o = static_cast<RingHistos*>(next())))
 256     o->fDebug = dbg;
 257 }
 258
 259 //====================================================================
 260 AliFMDEnergyFitter::RingHistos::RingHistos()
 261   : AliForwardUtil::RingHistos(),
 262     fEDist(0),
 263     fEmpty(0),
 264     fEtaEDists(),
 265     fList(0),
 266     fDebug(0)
 267 {}
 268
 269 //____________________________________________________________________
 270 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
 271   : AliForwardUtil::RingHistos(d,r),
 272     fEDist(0),
 273     fEmpty(0),
 274     fEtaEDists(),
 275     fList(0),
 276     fDebug(0)
 277 {
 278   fEtaEDists.SetName("EDists");
 279 }
 280 //____________________________________________________________________
 281 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
 282   : AliForwardUtil::RingHistos(o),
 283     fEDist(o.fEDist),
 284     fEmpty(o.fEmpty),
 285     fEtaEDists(),
 286     fList(0),
 287     fDebug(0)
 288 {
 289   TIter next(&o.fEtaEDists);
 290   TObject* obj = 0;
 291   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 292   if (o.fList) {
 293     fList = new TList;
 294     fList->SetName(fName);
 295     TIter next2(o.fList);
 296     while ((obj = next2())) fList->Add(obj->Clone());
 297   }
 298 }
 299
 300 //____________________________________________________________________
 301 AliFMDEnergyFitter::RingHistos&
 302 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
 303 {
 304   AliForwardUtil::RingHistos::operator=(o);
 305
 306   if (fEDist) delete fEDist;
 307   if (fEmpty) delete fEmpty;
 308   fEtaEDists.Delete();
 309   if (fList) fList->Delete();
 310
 311   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
 312   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
 313
 314   TIter next(&o.fEtaEDists);
 315   TObject* obj = 0;
 316   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 317
 318   if (o.fList) {
 319     fList = new TList;
 320     fList->SetName(fName);
 321     TIter next2(o.fList);
 322     while ((obj = next2())) fList->Add(obj->Clone());
 323   }
 324
 325   return *this;
 326 }
 327 //____________________________________________________________________
 328 AliFMDEnergyFitter::RingHistos::~RingHistos()
 329 {
 330   if (fEDist) delete fEDist;
 331   fEtaEDists.Delete();
 332 }
 333
 334 //____________________________________________________________________
 335 void
 336 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
 337 {
 338   if (empty) {
 339     fEmpty->Fill(mult);
 340     return;
 341   }
 342   fEDist->Fill(mult);
 343
 344   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
 345
 346   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
 347   if (!h) return;
 348
 349   h->Fill(mult);
 350 }
 351
 352 //__________________________________________________________________
 353 TArrayD
 354 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
 355                                                    Double_t max) const
 356 {
 357   // Service function to define a logarithmic axis.
 358   // Parameters:
 359   //   n    Number of bins
 360   //   min  Minimum of axis
 361   //   max  Maximum of axis
 362   TArrayD bins(n+1);
 363   Double_t dx = (max-min) / n;
 364   bins[0]     = min;
 365   Int_t    i  = 1;
 366   for (i = 1; i < n+1; i++) {
 367     Double_t dI   = float(i)/n;
 368     Double_t next = bins[i-1] + dx + dI * dI * dx;
 369     bins[i]       = next;
 370     if (next > max) break;
 371   }
 372   bins.Set(i+1);
 373   return bins;
 374 }
 375
 376 //____________________________________________________________________
 377 void
 378 AliFMDEnergyFitter::RingHistos::Make(Int_t    ieta,
 379                                      Double_t emin,
 380                                      Double_t emax,
 381                                      Double_t deMax,
 382                                      Int_t    nDeBins,
 383                                      Bool_t   incr)
 384 {
 385   TH1D* h = 0;
 386   TString name  = Form("%s_etabin%03d", fName.Data(), ieta);
 387   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
 388                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
 389   if (!incr)
 390     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
 391   else {
 392     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
 393     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
 394   }
 395
 396   h->SetDirectory(0);
 397   h->SetXTitle("#DeltaE/#DeltaE_{mip}");
 398   h->SetFillColor(Color());
 399   h->SetMarkerColor(Color());
 400   h->SetLineColor(Color());
 401   h->SetFillStyle(3001);
 402   h->Sumw2();
 403
 404   fEtaEDists.AddAt(h, ieta-1);
 405 }
 406 //____________________________________________________________________
 407 TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
 408                                               const char* title,
 409                                               const TAxis& eta) const
 410 {
 411   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 412                      Form("%s for %s", title, fName.Data()),
 413                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
 414   h->SetXTitle("#eta");
 415   h->SetYTitle(title);
 416   h->SetDirectory(0);
 417   h->SetFillColor(Color());
 418   h->SetMarkerColor(Color());
 419   h->SetLineColor(Color());
 420   h->SetFillStyle(3001);
 421   h->Sumw2();
 422
 423   return h;
 424 }
 425 //____________________________________________________________________
 426 TH1D*
 427 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
 428                                           const char* title,
 429                                           const TAxis& eta,
 430                                           Int_t low,
 431                                           Int_t high,
 432                                           Double_t val,
 433                                           Double_t err) const
 434 {
 435   Double_t xlow  = eta.GetBinLowEdge(low);
 436   Double_t xhigh = eta.GetBinUpEdge(high);
 437   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 438                      Form("%s for %s", title, fName.Data()),
 439                      1, xlow, xhigh);
 440   h->SetBinContent(1, val);
 441   h->SetBinError(1, err);
 442   h->SetXTitle("#eta");
 443   h->SetYTitle(title);
 444   h->SetDirectory(0);
 445   h->SetFillColor(0);
 446   h->SetMarkerColor(Color());
 447   h->SetLineColor(Color());
 448   h->SetFillStyle(0);
 449   h->SetLineStyle(2);
 450   h->SetLineWidth(2);
 451
 452   return h;
 453 }
 454
 455 //____________________________________________________________________
 456 void
 457 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
 458                                      Double_t maxDE, Int_t nDEbins,
 459                                      Bool_t useIncrBin)
 460 {
 461   fEDist = new TH1D(Form("%s_edist", fName.Data()),
 462                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()),
 463                     200, 0, 6);
 464   fEmpty = new TH1D(Form("%s_empty", fName.Data()),
 465                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)",
 466                          fName.Data()), 200, 0, 6);
 467   fList->Add(fEDist);
 468   fList->Add(fEmpty);
 469   // fEtaEDists.Expand(eAxis.GetNbins());
 470   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) {
 471     Double_t min = eAxis.GetBinLowEdge(i);
 472     Double_t max = eAxis.GetBinUpEdge(i);
 473     Make(i, min, max, maxDE, nDEbins, useIncrBin);
 474   }
 475   fList->Add(&fEtaEDists);
 476 }
 477 //____________________________________________________________________
 478 TObjArray*
 479 AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
 480                                     Double_t lowCut,
 481                                     UShort_t nParticles,
 482                                     UShort_t minEntries,
 483                                     UShort_t minusBins,
 484                                     Double_t relErrorCut,
 485                                     Double_t chi2nuCut) const
 486 {
 487   // Get our ring list from the output container
 488   TList* l = GetOutputList(dir);
 489   if (!l) return 0;
 490
 491   // Get the energy distributions from the output container
 492   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
 493   if (!dists) {
 494     AliWarning(Form("Didn't find %s_EtaEDists in %s",
 495                     fName.Data(), l->GetName()));
 496     l->ls();
 497     return 0;
 498   }
 499
 500   // Container of the fit results histograms
 501   TObjArray* pars  = new TObjArray(3+nParticles+1);
 502   pars->SetName("FitResults");
 503   l->Add(pars);
 504
 505   // Result objects stored in separate list on the output
 506   TH1* hChi2   = 0;
 507   TH1* hC      = 0;
 508   TH1* hDelta  = 0;
 509   TH1* hXi     = 0;
 510   TH1* hSigma  = 0;
 511   TH1* hSigmaN = 0;
 512   TH1* hN      = 0;
 513   TH1* hA[nParticles-1];
 514   pars->Add(hChi2   = MakePar("chi2",   "#chi^{2}/#nu", eta));
 515   pars->Add(hC      = MakePar("c",      "Constant",     eta));
 516   pars->Add(hDelta  = MakePar("delta",  "#Delta_{p}",   eta));
 517   pars->Add(hXi     = MakePar("xi",     "#xi",          eta));
 518   pars->Add(hSigma  = MakePar("sigma",  "#sigma",       eta));
 519   pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}",   eta));
 520   pars->Add(hN      = MakePar("n",      "N", eta));
 521   for (UShort_t i = 1; i < nParticles; i++)
 522     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
 523
 524
 525   Int_t nDists = dists->GetEntries();
 526   Int_t low    = nDists;
 527   Int_t high   = 0;
 528   for (Int_t i = 0; i < nDists; i++) {
 529     TH1D* dist = static_cast<TH1D*>(dists->At(i));
 530     // Ignore empty histograms altoghether
 531     if (!dist || dist->GetEntries() <= 0) continue;
 532
 533     // Scale to the bin-width
 534     dist->Scale(1., "width");
 535
 536     // Normalize peak to 1
 537     Double_t max = dist->GetMaximum();
 538     dist->Scale(1/max);
 539
 540     // Check that we have enough entries
 541     if (dist->GetEntries() <= minEntries) continue;
 542
 543     // Now fit
 544     TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
 545                        relErrorCut,chi2nuCut);
 546     if (!res) continue;
 547     // dist->GetListOfFunctions()->Add(res);
 548
 549     // Store eta limits
 550     low   = TMath::Min(low,i+1);
 551     high  = TMath::Max(high,i+1);
 552
 553     // Get the reduced chi square
 554     Double_t chi2 = res->GetChisquare();
 555     Int_t    ndf  = res->GetNDF();
 556
 557     // Store results of best fit in output histograms
 558     res->SetLineWidth(4);
 559     hChi2   ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
 560     hC      ->SetBinContent(i+1, res->GetParameter(kC));
 561     hDelta  ->SetBinContent(i+1, res->GetParameter(kDelta));
 562     hXi     ->SetBinContent(i+1, res->GetParameter(kXi));
 563     hSigma  ->SetBinContent(i+1, res->GetParameter(kSigma));
 564     hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));
 565     hN      ->SetBinContent(i+1, res->GetParameter(kN));
 566
 567     hC     ->SetBinError(i+1, res->GetParError(kC));
 568     hDelta ->SetBinError(i+1, res->GetParError(kDelta));
 569     hXi    ->SetBinError(i+1, res->GetParError(kXi));
 570     hSigma ->SetBinError(i+1, res->GetParError(kSigma));
 571     hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
 572     hN     ->SetBinError(i+1, res->GetParError(kN));
 573
 574     for (UShort_t j = 0; j < nParticles-1; j++) {
 575       hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
 576       hA[j]->SetBinError(i+1, res->GetParError(kA+j));
 577     }
 578   }
 579
 580   // Fit the full-ring histogram
 581   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
 582   if (total && total->GetEntries() >= minEntries) {
 583     TF1* res = FitHist(total,lowCut,nParticles,minusBins,
 584                        relErrorCut,chi2nuCut);
 585     if (res) {
 586       // Make histograms for the result of this fit
 587       Double_t chi2 = res->GetChisquare();
 588       Int_t    ndf  = res->GetNDF();
 589       pars->Add(MakeTotal("t_chi2",   "#chi^{2}/#nu", eta, low, high,
 590                           ndf > 0 ? chi2/ndf : 0, 0));
 591       pars->Add(MakeTotal("t_c",      "Constant",     eta, low, high,
 592                           res->GetParameter(kC),
 593                           res->GetParError(kC)));
 594       pars->Add(MakeTotal("t_delta",  "#Delta_{p}",   eta, low, high,
 595                           res->GetParameter(kDelta),
 596                           res->GetParError(kDelta)));
 597       pars->Add(MakeTotal("t_xi",     "#xi",          eta, low, high,
 598                           res->GetParameter(kXi),
 599                           res->GetParError(kXi)));
 600       pars->Add(MakeTotal("t_sigma",  "#sigma",       eta, low, high,
 601                           res->GetParameter(kSigma),
 602                           res->GetParError(kSigma)));
 603       pars->Add(MakeTotal("t_sigman", "#sigma_{n}",   eta, low, high,
 604                           res->GetParameter(kSigmaN),
 605                           res->GetParError(kSigmaN)));
 606       pars->Add(MakeTotal("t_n",    "N",              eta, low, high,
 607                           res->GetParameter(kN),0));
 608       for (UShort_t j = 0; j < nParticles-1; j++)
 609         pars->Add(MakeTotal(Form("t_a%d",j+2),
 610                             Form("a_{%d}",j+2), eta, low, high,
 611                             res->GetParameter(kA+j),
 612                             res->GetParError(kA+j)));
 613     }
 614   }
 615
 616   // Clean up list of histogram.  Histograms with no entries or
 617   // no functions are deleted.  We have to do this using the TObjLink
 618   // objects stored in the list since ROOT cannot guaranty the validity
 619   // of iterators when removing from a list - tsck.  Should just implement
 620   // TIter::Remove().
 621   TObjLink* lnk = dists->FirstLink();
 622   while (lnk) {
 623     TH1* h = static_cast<TH1*>(lnk->GetObject());
 624     if (h->GetEntries() <= 0 ||
 625         h->GetListOfFunctions()->GetEntries() <= 0) {
 626       TObjLink* keep = lnk->Next();
 627       dists->Remove(lnk);
 628       lnk = keep;
 629       continue;
 630     }
 631     lnk = lnk->Next();
 632   }
 633
 634   return pars;
 635 }
 636
 637 //____________________________________________________________________
 638 TF1*
 639 AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
 640                                         Double_t lowCut,
 641                                         UShort_t nParticles,
 642                                         UShort_t minusBins,
 643                                         Double_t relErrorCut,
 644                                         Double_t chi2nuCut) const
 645 {
 646   Double_t maxRange = 10;
 647
 648   // Create a fitter object
 649   AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins);
 650   f.Clear();
 651
 652
 653   // If we are only asked to fit a single particle, return this fit,
 654   // no matter what.
 655   if (nParticles == 1) {
 656     TF1* r = f.Fit1Particle(dist, 0);
 657     if (!r) return 0;
 658     return new TF1(*r);
 659   }
 660
 661   // Fit from 2 upto n particles
 662   for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
 663
 664
 665   // Now, we need to select the best fit
 666   Int_t nFits = f.fFitResults.GetEntriesFast();
 667   TF1*  good[nFits];
 668   for (Int_t i = nFits-1; i >= 0; i--) {
 669     good[i] = 0;
 670     TF1* ff = static_cast<TF1*>(f.fFunctions.At(i));
 671     // ff->SetLineColor(Color());
 672     ff->SetRange(0, maxRange);
 673     dist->GetListOfFunctions()->Add(new TF1(*ff));
 674     if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i)),
 675                     relErrorCut, chi2nuCut)) {
 676       good[i] = ff;
 677       ff->SetLineWidth(2);
 678       // f.fFitResults.At(i)->Print("V");
 679     }
 680   }
 681   // If all else fails, use the 1 particle fit
 682   TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
 683
 684   // Find the fit with the most valid particles
 685   for (Int_t i = nFits-1; i > 0; i--) {
 686     if (!good[i]) continue;
 687     if (fDebug > 1) {
 688       AliInfo(Form("Choosing fit with n=%d", i+1));
 689       f.fFitResults.At(i)->Print();
 690     }
 691     ret = good[i];
 692     break;
 693   }
 694   // Give a warning if we're using fall-back
 695   if (ret == f.fFunctions.At(0))
 696     AliWarning("Choosing fall-back 1 particle fit");
 697
 698   // Copy our result and return (the functions are owned by the fitter)
 699   TF1* fret = new TF1(*ret);
 700   return fret;
 701 }
 702
 703 //____________________________________________________________________
 704 Bool_t
 705 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
 706                                             Double_t relErrorCut,
 707                                             Double_t chi2nuCut) const
 708 {
 709   if (fDebug > 10) r->Print();
 710   TString  n    = r->GetName();
 711   n.ReplaceAll("TFitResult-", "");
 712   Double_t chi2 = r->Chi2();
 713   Int_t    ndf  = r->Ndf();
 714   // Double_t prob = r.Prob();
 715   Bool_t ret = kTRUE;
 716
 717   // Check that the reduced chi square isn't larger than 5
 718   if (ndf <= 0 || chi2 / ndf > chi2nuCut) {
 719     if (fDebug > 2)
 720       AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f",
 721                       n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
 722                       chi2nuCut));
 723     ret = kFALSE;
 724   }
 725
 726   // Check each parameter
 727   UShort_t nPar = r->NPar();
 728   for (UShort_t i = 0; i < nPar; i++) {
 729     if (i == kN) continue;  // Skip the number parameter
 730
 731     // Get value and error and check value
 732     Double_t v  = r->Parameter(i);
 733     Double_t e  = r->ParError(i);
 734     if (v == 0) continue;
 735
 736     // Calculate the relative error and check it
 737     Double_t re  = e / v;
 738     Double_t cut = relErrorCut * (i < kN ? 1 : 2);
 739     if (re > cut) {
 740       if (fDebug > 2)
 741         AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
 742                         n.Data(), r->ParName(i).c_str(),
 743                         r->ParName(i).c_str(), e, v,
 744                         100*(v == 0 ? 0 : e/v),
 745                         100*(cut)));
 746       ret = kFALSE;
 747     }
 748   }
 749
 750   // Check if we have scale parameters
 751   if (nPar > kN) {
 752
 753     // Check that the last particle has a significant contribution
 754     Double_t lastScale = r->Parameter(nPar-1);
 755     if (lastScale <= 1e-7) {
 756       if (fDebug)
 757         AliWarning(Form("%s: %s=%9.6f<1e-7",
 758                         n.Data(), r->ParName(nPar-1).c_str(), lastScale));
 759       ret = kFALSE;
 760     }
 761   }
 762   return ret;
 763 }
 764
 765
 766 //____________________________________________________________________
 767 void
 768 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
 769 {
 770   fList = DefineOutputList(dir);
 771 }
 772
 773 //____________________________________________________________________
 774 //
 775 // EOF
 776 //