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