PWGLF/FORWARD/analysis2/AliFMDEnergyFitter.cxx

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