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 #define N_A(N)  (4+N-2)
  23 #define N2_A(N) (4+(N-2)*3)
  24 #define N2_D(N) (4+(N-2)*3+1)
  25 #define N2_X(N) (4+(N-2)*3+2)
  26
  27 //____________________________________________________________________
  28 namespace {
  29   Double_t
  30   NLandau(Double_t* xp, Double_t* pp)
  31   {
  32     Double_t  e        = xp[0];
  33     Double_t  constant = pp[0];
  34     Double_t  mpv      = pp[1];
  35     Double_t  fwhm     = pp[2];
  36     Int_t     n        = Int_t(pp[3]);
  37     Double_t  result   = 0;
  38     for (Int_t i = 1; i <= n; i++) {
  39       Double_t mpvi  =  i*(mpv+fwhm*TMath::Log(i));
  40       Double_t fwhmi =  i*fwhm;
  41       Double_t ai    =  (i == 1 ? 1 : pp[N_A(i)]);
  42       result         += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE);
  43     }
  44     result *= constant;
  45     return result;
  46   }
  47
  48   Double_t
  49   NLandau2(Double_t* xp, Double_t* pp)
  50   {
  51     Double_t  e        = xp[0];
  52     Double_t  constant = pp[0];
  53     Double_t  mpv      = pp[1];
  54     Double_t  fwhm     = pp[2];
  55     Int_t     n        = Int_t(pp[3]);
  56     Double_t  result   = TMath::Landau(e,mpv,fwhm,kTRUE);
  57     for (Int_t i = 2; i <= n; i++) {
  58       Double_t ai    =  pp[N2_A(i)];
  59       Double_t mpvi  =  pp[N2_D(i)];
  60       Double_t fwhmi =  pp[N2_X(i)];
  61       result         += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE);
  62     }
  63     result *= constant;
  64     return result;
  65   }
  66
  67   Double_t
  68   TripleLandau(Double_t *x, Double_t *par)
  69   {
  70     Double_t energy   = x[0];
  71     Double_t constant = par[0];
  72     Double_t mpv      = par[1];
  73     Double_t fwhm     = par[2];
  74     Double_t alpha    = par[3];
  75     Double_t beta     = par[4];
  76     Double_t mpv2     = 2*mpv+2*fwhm*TMath::Log(2);
  77     Double_t mpv3     = 3*mpv+3*fwhm*TMath::Log(3);
  78
  79     Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+
  80                              alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE)+
  81                              beta  * TMath::Landau(energy,mpv3,3*fwhm,kTRUE));
  82
  83     return f;
  84   }
  85   Double_t
  86   DoubleLandau(Double_t *x, Double_t *par)
  87   {
  88     Double_t energy   = x[0];
  89     Double_t constant = par[0];
  90     Double_t mpv      = par[1];
  91     Double_t fwhm     = par[2];
  92     Double_t alpha    = par[3];
  93     Double_t mpv2     = 2*mpv+2*fwhm*TMath::Log(2);
  94
  95     Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+
  96                              alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE));
  97
  98     return f;
  99   }
 100   Double_t
 101   SingleLandau(Double_t *x, Double_t *par)
 102   {
 103     Double_t energy   = x[0];
 104     Double_t constant = par[0];
 105     Double_t mpv      = par[1];
 106     Double_t fwhm     = par[2];
 107
 108     Double_t f = constant * TMath::Landau(energy,mpv,fwhm,kTRUE);
 109
 110     return f;
 111   }
 112 }
 113
 114 //____________________________________________________________________
 115 AliFMDEnergyFitter::AliFMDEnergyFitter()
 116   : TNamed(),
 117     fRingHistos(),
 118     fLowCut(0.3),
 119     fNLandau(3),
 120     fMinEntries(100),
 121     fBinsToSubtract(4),
 122     fDoFits(false),
 123     fEtaAxis(),
 124     fDebug(0)
 125 {}
 126
 127 //____________________________________________________________________
 128 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
 129   : TNamed("fmdEnergyFitter", title),
 130     fRingHistos(),
 131     fLowCut(0.3),
 132     fNLandau(3),
 133     fMinEntries(100),
 134     fBinsToSubtract(4),
 135     fDoFits(false),
 136     fEtaAxis(100,-4,6),
 137     fDebug(3)
 138 {
 139   fEtaAxis.SetName("etaAxis");
 140   fEtaAxis.SetTitle("#eta");
 141   fRingHistos.Add(new RingHistos(1, 'I'));
 142   fRingHistos.Add(new RingHistos(2, 'I'));
 143   fRingHistos.Add(new RingHistos(2, 'O'));
 144   fRingHistos.Add(new RingHistos(3, 'I'));
 145   fRingHistos.Add(new RingHistos(3, 'O'));
 146 }
 147
 148 //____________________________________________________________________
 149 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
 150   : TNamed(o),
 151     fRingHistos(),
 152     fLowCut(o.fLowCut),
 153     fNLandau(o.fNLandau),
 154     fMinEntries(o.fMinEntries),
 155     fBinsToSubtract(o.fBinsToSubtract),
 156     fDoFits(o.fDoFits),
 157     fEtaAxis(o.fEtaAxis),
 158     fDebug(o.fDebug)
 159 {
 160   TIter    next(&o.fRingHistos);
 161   TObject* obj = 0;
 162   while ((obj = next())) fRingHistos.Add(obj);
 163 }
 164
 165 //____________________________________________________________________
 166 AliFMDEnergyFitter::~AliFMDEnergyFitter()
 167 {
 168   fRingHistos.Delete();
 169 }
 170
 171 //____________________________________________________________________
 172 AliFMDEnergyFitter&
 173 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
 174 {
 175   TNamed::operator=(o);
 176
 177   fLowCut        = o.fLowCut;
 178   fNLandau       = o.fNLandau;
 179   fMinEntries    = o.fMinEntries;
 180   fBinsToSubtract= o.fBinsToSubtract;
 181   fDoFits        = o.fDoFits;
 182   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
 183   fDebug         = o.fDebug;
 184
 185   fRingHistos.Delete();
 186   TIter    next(&o.fRingHistos);
 187   TObject* obj = 0;
 188   while ((obj = next())) fRingHistos.Add(obj);
 189
 190   return *this;
 191 }
 192
 193 //____________________________________________________________________
 194 AliFMDEnergyFitter::RingHistos*
 195 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
 196 {
 197   Int_t idx = -1;
 198   switch (d) {
 199   case 1: idx = 0; break;
 200   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 201   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 202   }
 203   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 204
 205   return static_cast<RingHistos*>(fRingHistos.At(idx));
 206 }
 207
 208 //____________________________________________________________________
 209 void
 210 AliFMDEnergyFitter::Init(const TAxis& eAxis)
 211 {
 212   fEtaAxis.Set(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
 213   TIter    next(&fRingHistos);
 214   RingHistos* o = 0;
 215   while ((o = static_cast<RingHistos*>(next())))
 216     o->Init(eAxis);
 217 }
 218 //____________________________________________________________________
 219 Bool_t
 220 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
 221                                Bool_t           empty)
 222 {
 223   for(UShort_t d = 1; d <= 3; d++) {
 224     Int_t nRings = (d == 1 ? 1 : 2);
 225     for (UShort_t q = 0; q < nRings; q++) {
 226       Char_t      r    = (q == 0 ? 'I' : 'O');
 227       UShort_t    nsec = (q == 0 ?  20 :  40);
 228       UShort_t    nstr = (q == 0 ? 512 : 256);
 229
 230       RingHistos* histos = GetRingHistos(d, r);
 231
 232       for(UShort_t s = 0; s < nsec;  s++) {
 233         for(UShort_t t = 0; t < nstr; t++) {
 234           Float_t mult = input.Multiplicity(d,r,s,t);
 235
 236           // Keep dead-channel information.
 237           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0)
 238             continue;
 239
 240           // Get the pseudo-rapidity
 241           Double_t eta1 = input.Eta(d,r,s,t);
 242           Int_t ieta = fEtaAxis.FindBin(eta1);
 243           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue;
 244
 245           histos->Fill(empty, ieta-1, mult);
 246
 247         } // for strip
 248       } // for sector
 249     } // for ring
 250   } // for detector
 251
 252   return kTRUE;
 253 }
 254
 255 //____________________________________________________________________
 256 void
 257 AliFMDEnergyFitter::Fit(TList* dir)
 258 {
 259   if (!fDoFits) return;
 260
 261   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 262   if (!d) return;
 263
 264   // +1          for chi^2
 265   // +3          for 1 landau
 266   // +1          for N
 267   // +fNLandau-1 for weights
 268   Int_t nStack = 1+3+1+fNLandau-1;
 269   THStack* stack[nStack];
 270   stack[0] = new THStack("chi2", "#chi^{2}/#nu");
 271   stack[1] = new THStack("c",    "constant");
 272   stack[2] = new THStack("mpv",  "#Delta_{p}");
 273   stack[3] = new THStack("w",    "w");
 274   stack[4] = new THStack("n",    "# of Landaus");
 275   for (Int_t i = 2; i <= fNLandau; i++)
 276     stack[i-1+4] = new THStack(Form("a%d", i),
 277                                Form("Weight of %d signal", i));
 278   for (Int_t i = 0; i < nStack; i++)
 279     d->Add(stack[i]);
 280
 281   TIter    next(&fRingHistos);
 282   RingHistos* o = 0;
 283   while ((o = static_cast<RingHistos*>(next()))) {
 284     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNLandau,
 285                           fMinEntries, fBinsToSubtract);
 286     if (!l) continue;
 287     for (Int_t i = 0; i < l->GetEntriesFast(); i++) {
 288       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
 289     }
 290   }
 291 }
 292
 293 //____________________________________________________________________
 294 void
 295 AliFMDEnergyFitter::DefineOutput(TList* dir)
 296 {
 297   TList* d = new TList;
 298   d->SetName(GetName());
 299   dir->Add(d);
 300   d->Add(&fEtaAxis);
 301   TIter    next(&fRingHistos);
 302   RingHistos* o = 0;
 303   while ((o = static_cast<RingHistos*>(next()))) {
 304     o->Output(d);
 305   }
 306 }
 307 //____________________________________________________________________
 308 void
 309 AliFMDEnergyFitter::SetDebug(Int_t dbg)
 310 {
 311   fDebug = dbg;
 312   TIter    next(&fRingHistos);
 313   RingHistos* o = 0;
 314   while ((o = static_cast<RingHistos*>(next())))
 315     o->fDebug = dbg;
 316 }
 317
 318 //====================================================================
 319 AliFMDEnergyFitter::RingHistos::RingHistos()
 320   : AliForwardUtil::RingHistos(),
 321     fEDist(0),
 322     fEmpty(0),
 323     fEtaEDists(),
 324     fList(0),
 325     fDebug(0)
 326 {}
 327
 328 //____________________________________________________________________
 329 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
 330   : AliForwardUtil::RingHistos(d,r),
 331     fEDist(0),
 332     fEmpty(0),
 333     fEtaEDists(),
 334     fList(0),
 335     fDebug(0)
 336 {
 337   fEtaEDists.SetName("EDists");
 338 }
 339 //____________________________________________________________________
 340 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
 341   : AliForwardUtil::RingHistos(o),
 342     fEDist(o.fEDist),
 343     fEmpty(o.fEmpty),
 344     fEtaEDists(),
 345     fList(0),
 346     fDebug(0)
 347 {
 348   TIter next(&o.fEtaEDists);
 349   TObject* obj = 0;
 350   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 351   if (o.fList) {
 352     fList = new TList;
 353     fList->SetName(fName);
 354     TIter next2(o.fList);
 355     while ((obj = next2())) fList->Add(obj->Clone());
 356   }
 357 }
 358
 359 //____________________________________________________________________
 360 AliFMDEnergyFitter::RingHistos&
 361 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
 362 {
 363   AliForwardUtil::RingHistos::operator=(o);
 364
 365   if (fEDist) delete fEDist;
 366   if (fEmpty) delete fEmpty;
 367   fEtaEDists.Delete();
 368   if (fList) fList->Delete();
 369
 370   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
 371   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
 372
 373   TIter next(&o.fEtaEDists);
 374   TObject* obj = 0;
 375   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 376
 377   if (o.fList) {
 378     fList = new TList;
 379     fList->SetName(fName);
 380     TIter next2(o.fList);
 381     while ((obj = next2())) fList->Add(obj->Clone());
 382   }
 383
 384   return *this;
 385 }
 386 //____________________________________________________________________
 387 AliFMDEnergyFitter::RingHistos::~RingHistos()
 388 {
 389   if (fEDist) delete fEDist;
 390   fEtaEDists.Delete();
 391 }
 392
 393 //____________________________________________________________________
 394 void
 395 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
 396 {
 397   if (empty) {
 398     fEmpty->Fill(mult);
 399     return;
 400   }
 401   fEDist->Fill(mult);
 402
 403   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
 404
 405   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
 406   if (!h) return;
 407
 408   h->Fill(mult);
 409 }
 410
 411 //__________________________________________________________________
 412 TArrayD
 413 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
 414                                                    Double_t max) const
 415 {
 416   // Service function to define a logarithmic axis.
 417   // Parameters:
 418   //   n    Number of bins
 419   //   min  Minimum of axis
 420   //   max  Maximum of axis
 421   TArrayD bins(n+1);
 422   Double_t dx = (max-min) / n;
 423   bins[0]     = min;
 424   Int_t    i  = 1;
 425   for (i = 1; i < n+1; i++) {
 426     Double_t dI   = float(i)/n;
 427     Double_t next = bins[i-1] + dx + dI * dI * dx;
 428     bins[i]       = next;
 429     if (next > max) break;
 430   }
 431   bins.Set(i+1);
 432   return bins;
 433 }
 434
 435 //____________________________________________________________________
 436 void
 437 AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax,
 438                                      Double_t deMax, Int_t nDeBins,
 439                                      Bool_t incr)
 440 {
 441   TH1D* h = 0;
 442   TString name  = Form("%s_etabin%03d", fName.Data(), ieta);
 443   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
 444                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
 445   if (!incr)
 446     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
 447   else {
 448     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
 449     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
 450   }
 451
 452   h->SetDirectory(0);
 453   h->SetXTitle("#DeltaE/#DeltaE_{mip}");
 454   h->SetFillColor(Color());
 455   h->SetMarkerColor(Color());
 456   h->SetLineColor(Color());
 457   h->SetFillStyle(3001);
 458   h->Sumw2();
 459
 460   fEtaEDists.AddAt(h, ieta-1);
 461 }
 462 //____________________________________________________________________
 463 TH1D*
 464 AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
 465                                         const char* title,
 466                                         const TAxis& eta) const
 467 {
 468   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 469                      Form("%s for %s", title, fName.Data()),
 470                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
 471   h->SetXTitle("#eta");
 472   h->SetYTitle(title);
 473   h->SetDirectory(0);
 474   h->SetFillColor(Color());
 475   h->SetMarkerColor(Color());
 476   h->SetLineColor(Color());
 477   h->SetFillStyle(3001);
 478   h->Sumw2();
 479
 480   return h;
 481 }
 482 //____________________________________________________________________
 483 TH1D*
 484 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
 485                                           const char* title,
 486                                           const TAxis& eta,
 487                                           Int_t low,
 488                                           Int_t high,
 489                                           Double_t val,
 490                                           Double_t err) const
 491 {
 492   Double_t xlow  = eta.GetBinLowEdge(low);
 493   Double_t xhigh = eta.GetBinUpEdge(high);
 494   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 495                      Form("%s for %s", title, fName.Data()),
 496                      1, xlow, xhigh);
 497   h->SetBinContent(1, val);
 498   h->SetBinError(1, err);
 499   h->SetXTitle("#eta");
 500   h->SetYTitle(title);
 501   h->SetDirectory(0);
 502   h->SetFillColor(0);
 503   h->SetMarkerColor(Color());
 504   h->SetLineColor(Color());
 505   h->SetFillStyle(0);
 506   h->SetLineStyle(2);
 507   h->SetLineWidth(2);
 508
 509   return h;
 510 }
 511
 512 //____________________________________________________________________
 513 void
 514 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
 515                                      Double_t maxDE, Int_t nDEbins,
 516                                      Bool_t useIncrBin)
 517 {
 518   fEDist = new TH1D(Form("%s_edist", fName.Data()),
 519                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()),
 520                     200, 0, 6);
 521   fEmpty = new TH1D(Form("%s_empty", fName.Data()),
 522                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)",
 523                          fName.Data()), 200, 0, 6);
 524   fList->Add(fEDist);
 525   fList->Add(fEmpty);
 526   // fEtaEDists.Expand(eAxis.GetNbins());
 527   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) {
 528     Double_t min = eAxis.GetBinLowEdge(i);
 529     Double_t max = eAxis.GetBinUpEdge(i);
 530     Make(i, min, max, maxDE, nDEbins, useIncrBin);
 531   }
 532   fList->Add(&fEtaEDists);
 533 }
 534 //____________________________________________________________________
 535 TObjArray*
 536 AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
 537                                     Double_t lowCut, UShort_t nLandau,
 538                                     UShort_t minEntries,
 539                                     UShort_t minusBins) const
 540 {
 541   TList* l = GetOutputList(dir);
 542   if (!l) return 0;
 543
 544   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
 545   if (!dists) {
 546     AliWarning(Form("Didn't find %s_EtaEDists in %s",
 547                     fName.Data(), l->GetName()));
 548     l->ls();
 549     return 0;
 550   }
 551
 552   TObjArray* pars  = new TObjArray(3+nLandau+1);
 553   pars->SetName("FitResults");
 554   l->Add(pars);
 555
 556   TH1* hChi2 = 0;
 557   TH1* hC    = 0;
 558   TH1* hMpv  = 0;
 559   TH1* hW    = 0;
 560   TH1* hS    = 0;
 561   TH1* hN    = 0;
 562   TH1* hA[nLandau-1];
 563   pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta));
 564   pars->Add(hC    = MakePar("c",    "Constant", eta));
 565   pars->Add(hMpv  = MakePar("mpv",  "#Delta_{p}", eta));
 566   pars->Add(hW    = MakePar("w",    "#xi", eta));
 567   pars->Add(hS    = MakePar("s",    "#sigma", eta));
 568   pars->Add(hN    = MakePar("n",    "N", eta));
 569   for (UShort_t i = 1; i < nLandau; i++)
 570     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
 571
 572
 573   Int_t nDists = dists->GetEntries();
 574   Int_t low    = nDists;
 575   Int_t high   = 0;
 576   for (Int_t i = 0; i < nDists; i++) {
 577     TH1D* dist = static_cast<TH1D*>(dists->At(i));
 578
 579     if (!dist) continue;
 580
 581     TF1* res = FitHist(dist,lowCut,nLandau,minEntries,minusBins);
 582     if (!res) continue;
 583
 584     low   = TMath::Min(low,i+1);
 585     high  = TMath::Max(high,i+1);
 586
 587     Double_t chi2 = res->GetChisquare();
 588     Int_t    ndf  = res->GetNDF();
 589     hChi2->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
 590     hC  ->SetBinContent(i+1, res->GetParameter(0));
 591     hMpv->SetBinContent(i+1, res->GetParameter(1));
 592     hW  ->SetBinContent(i+1, res->GetParameter(2));
 593     hN  ->SetBinContent(i+1, res->GetParameter(3));
 594
 595     hC  ->SetBinError(i+1, res->GetParError(1));
 596     hMpv->SetBinError(i+1, res->GetParError(2));
 597     hW  ->SetBinError(i+1, res->GetParError(2));
 598     // hN  ->SetBinError(i, res->GetParError(3));
 599
 600     for (UShort_t j = 0; j < nLandau-1; j++) {
 601       hA[j]->SetBinContent(i+1, res->GetParameter(4+j));
 602       hA[j]->SetBinError(i+1, res->GetParError(4+j));
 603     }
 604   }
 605
 606   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
 607   if (total) {
 608     TF1* res = FitHist(total,lowCut,nLandau,minEntries,minusBins);
 609     if (res) {
 610       Double_t chi2 = res->GetChisquare();
 611       Int_t    ndf  = res->GetNDF();
 612       pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
 613                           ndf > 0 ? chi2/ndf : 0, 0));
 614       pars->Add(MakeTotal("t_c",    "Constant",     eta, low, high,
 615                           res->GetParameter(0),res->GetParError(0)));
 616       pars->Add(MakeTotal("t_mpv",  "#Delta_{p}",   eta, low, high,
 617                           res->GetParameter(1),res->GetParError(1)));
 618       pars->Add(MakeTotal("t_w",    "#xi",          eta, low, high,
 619                           res->GetParameter(2),res->GetParError(2)));
 620       pars->Add(MakeTotal("t_n",    "N",            eta, low, high,
 621                           res->GetParameter(3),0));
 622       for (UShort_t j = 1; j < nLandau; j++)
 623         pars->Add(MakeTotal(Form("a%d_t",j+1),
 624                             Form("a_{%d}",j+1), eta, low, high,
 625                             res->GetParameter(3+j), res->GetParError(3+j)));
 626     }
 627   }
 628
 629   TObjLink* lnk = dists->FirstLink();
 630   while (lnk) {
 631     TH1* h = static_cast<TH1*>(lnk->GetObject());
 632     if (h->GetEntries() <= 0 ||
 633         h->GetListOfFunctions()->GetEntries() <= 0) {
 634       TObjLink* keep = lnk->Next();
 635       dists->Remove(lnk);
 636       lnk = keep;
 637       continue;
 638     }
 639     lnk = lnk->Next();
 640   }
 641
 642   return pars;
 643 }
 644
 645 //____________________________________________________________________
 646 TF1*
 647 AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,Double_t lowCut,
 648                                         UShort_t nLandau,
 649                                         UShort_t minEntries,
 650                                         UShort_t minusBins) const
 651 {
 652   Double_t maxRange = 10;
 653
 654   if (dist->GetEntries() <= minEntries) return 0;
 655
 656   // Find the fit range
 657   dist->GetXaxis()->SetRangeUser(lowCut, maxRange);
 658
 659   // Normalize peak to 1
 660   Double_t max = dist->GetMaximum();
 661   dist->Scale(1/max);
 662
 663   // Get the bin with maximum
 664   Int_t    maxBin = dist->GetMaximumBin();
 665   Double_t maxE   = dist->GetBinLowEdge(maxBin);
 666
 667   // Get the low edge
 668   dist->GetXaxis()->SetRangeUser(lowCut, maxE);
 669   Int_t    minBin = maxBin - minusBins; // dist->GetMinimumBin();
 670   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),lowCut);
 671   Double_t maxEE  = dist->GetBinCenter(maxBin+2*minusBins);
 672
 673   // Restore the range
 674   dist->GetXaxis()->SetRangeUser(0, maxRange);
 675
 676   // First do a single landau fit
 677   TF1*          landau1 = new TF1("landau1", "landau", minE, maxEE);
 678   // landau1->SetParameters(1,1,1,1);
 679   landau1->SetParNames("C","#Delta_{p}","#xi");
 680   landau1->SetParLimits(1,minE,maxRange);
 681   landau1->SetParLimits(2,0,maxRange);
 682   landau1->SetLineColor(Color());
 683   // Tight fit around peak - make sure we get that right.
 684   TFitResultPtr r = dist->Fit(landau1, "RQNS", "", minE, maxEE);
 685
 686   return FitMore2(dist, nLandau, *r, landau1, minE, maxRange);
 687 }
 688
 689 //____________________________________________________________________
 690 TF1*
 691 AliFMDEnergyFitter::RingHistos::FitMore(TH1*        dist,
 692                                         UShort_t    nLandau,
 693                                         TFitResult& r,
 694                                         TF1*        landau1,
 695                                         Double_t    minE,
 696                                         Double_t    maxRange) const
 697 {
 698   static TClonesArray res("TFitResult");
 699   static TObjArray funcs;
 700   res.Clear();
 701   funcs.SetOwner();
 702   funcs.Clear();
 703   Int_t nRes = 0;
 704
 705   //r->Print();
 706   new(res[nRes++]) TFitResult(r);
 707   funcs.AddAtAndExpand(landau1, 0);
 708
 709   // Now try to fit
 710   for (UShort_t n = 2; n <= nLandau; n++) {
 711     TFitResult* rr = static_cast<TFitResult*>(res.At(nRes-1));
 712     if (!rr) {
 713       AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n));
 714       return 0;
 715     }
 716     // Create the function object
 717     Double_t mpvi  = rr->Parameter(1);
 718     Double_t wi    = rr->Parameter(2);
 719     Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi;
 720     TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,3+n);
 721     landaui->SetLineStyle((n % 10)+1);
 722     landaui->SetLineWidth(2);
 723     landaui->SetLineColor(Color());
 724     landaui->SetParameter(0, rr->Parameter(0));
 725     landaui->SetParameter(1, rr->Parameter(1));
 726     landaui->SetParameter(2, rr->Parameter(2));
 727     landaui->SetParLimits(1, minE, maxRange);
 728     landaui->SetParLimits(2,0,maxRange);
 729     landaui->FixParameter(3, n);
 730     landaui->SetParNames("C","#Delta_{p}","#xi", "N");
 731     for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit
 732       landaui->SetParameter(N_A(i), rr->Parameter(N_A(i)));
 733       landaui->SetParLimits(N_A(i), 0,1);
 734       landaui->SetParName(i, Form("a_{%d}", i));
 735     }
 736     landaui->SetParameter(N_A(n), (n == 2 ? 0.05 : rr->Parameter(N_A(n-1))/5));
 737     landaui->SetParLimits(N_A(n), 0, 1);
 738     landaui->SetParName(N_A(n), Form("a_{%d}", n));
 739     // landaui->Print();
 740
 741     TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi);
 742     // tr->Print();
 743     if (CheckResult(*tr)) {
 744       new(res[nRes++]) TFitResult(*tr);
 745       funcs.AddAtAndExpand(landaui,n-1);
 746       continue;
 747     }
 748     // Stop on bad fit
 749     break;
 750   }
 751   TF1* ret = 0;
 752   if (funcs.At(nRes-1)) ret = static_cast<TF1*>(funcs.At(nRes-1)->Clone());
 753   dist->GetListOfFunctions()->Add(ret);
 754
 755   res.Clear();
 756   funcs.Clear();
 757
 758   return ret;
 759 }
 760
 761 //____________________________________________________________________
 762 TF1*
 763 AliFMDEnergyFitter::RingHistos::FitMore2(TH1*        dist,
 764                                          UShort_t    nLandau,
 765                                          TFitResult& r,
 766                                          TF1*        landau1,
 767                                          Double_t    minE,
 768                                          Double_t    maxRange) const
 769 {
 770   static TClonesArray res("TFitResult");
 771   static TObjArray funcs;
 772   res.Clear();
 773   funcs.SetOwner();
 774   funcs.Clear();
 775   Int_t nRes = 0;
 776
 777   //r->Print();
 778   new(res[nRes++]) TFitResult(r);
 779   funcs.AddAtAndExpand(landau1, 0);
 780
 781   // Now try to fit
 782   for (UShort_t n = 2; n <= nLandau; n++) {
 783     TFitResult* rr = static_cast<TFitResult*>(res.At(nRes-1));
 784     if (!rr) {
 785       AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n));
 786       return 0;
 787     }
 788     // Create the function object
 789     Double_t mpvi  = rr->Parameter(1);
 790     Double_t wi    = rr->Parameter(2);
 791     Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi;
 792     TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,
 793                            n*3+1);
 794     landaui->SetLineStyle((n % 10)+1);
 795     landaui->SetLineWidth(2);
 796     landaui->SetLineColor(Color());
 797     landaui->SetParameter(0, rr->Parameter(0));
 798     landaui->SetParameter(1, rr->Parameter(1));
 799     landaui->SetParameter(2, rr->Parameter(2));
 800     landaui->SetParLimits(1, minE, maxRange);
 801     landaui->SetParLimits(2,0,maxRange);
 802     landaui->FixParameter(3, n);
 803     landaui->SetParNames("C","#Delta_{p}","#xi", "N");
 804     for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit
 805       landaui->SetParameter(N2_A(i), rr->Parameter(N2_A(i)));
 806       landaui->SetParameter(N2_D(i), rr->Parameter(N2_D(i)));
 807       landaui->SetParameter(N2_X(i), rr->Parameter(N2_X(i)));
 808       landaui->SetParLimits(N2_A(i), 0,1);
 809       landaui->SetParLimits(N2_D(i), minE,maxEi);
 810       landaui->SetParLimits(N2_X(i), 0,maxRange);
 811       landaui->SetParName(N2_A(i), Form("a_{%d}", i));
 812       landaui->SetParName(N2_D(i), Form("#Delta_{p,%d}", i));
 813       landaui->SetParName(N2_X(i), Form("#xi_{%d}", i));
 814     }
 815     landaui->SetParameter(N2_A(n), n == 2 ? 0.05 : rr->Parameter(N2_A(n-1))/5);
 816     landaui->SetParameter(N2_D(n), mpvi);
 817     landaui->SetParameter(N2_X(n), wi);
 818     landaui->SetParLimits(N2_A(n), 0,1);
 819     landaui->SetParLimits(N2_D(n), minE,maxEi);
 820     landaui->SetParLimits(N2_X(n), 0,maxRange);
 821     landaui->SetParName(N2_A(n), Form("a_{%d}", n));
 822     landaui->SetParName(N2_D(n), Form("#Delta_{p,%d}", n));
 823     landaui->SetParName(N2_X(n), Form("#xi_{%d}", n));
 824     if (fDebug > 2) landaui->Print();
 825
 826     TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi);
 827     if (fDebug > 2)  tr->Print();
 828     if (CheckResult(*tr)) {
 829       new(res[nRes++]) TFitResult(*tr);
 830       funcs.AddAtAndExpand(landaui,n-1);
 831       continue;
 832     }
 833     // Stop on bad fit
 834     break;
 835   }
 836   TF1* ret = 0;
 837   if (funcs.At(nRes-1)) ret = static_cast<TF1*>(funcs.At(nRes-1)->Clone());
 838   dist->GetListOfFunctions()->Add(ret);
 839
 840   res.Clear();
 841   funcs.Clear();
 842
 843   return ret;
 844 }
 845
 846 //____________________________________________________________________
 847 Bool_t
 848 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult& r) const
 849 {
 850   Double_t chi2 = r.Chi2();
 851   Int_t    ndf  = r.Ndf();
 852   // Double_t prob = r.Prob();
 853   if (ndf <= 0 || chi2 / ndf > 5) {
 854     if (fDebug > 2)
 855       AliWarning(Form("Fit %s has chi^2/ndf=%f/%d=%f",
 856                       r.GetName(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf)));
 857     return kFALSE;
 858   }
 859
 860   UShort_t nPar = r.NPar();
 861   for (UShort_t i = 0; i < nPar; i++) {
 862     if (i == 3) continue;
 863
 864     Double_t v = r.Parameter(i);
 865     Double_t e = r.ParError(i);
 866     if (v == 0) continue;
 867     if (v == 0 || e / v > 0.2) {
 868       if (fDebug > 2)
 869         AliWarning(Form("Fit %s has Delta %s/%s=%f/%f=%f%%",
 870                         r.GetName(), r.ParName(i).c_str(),
 871                         r.ParName(i).c_str(), e, v, 100*(v == 0 ? 0 : e/v)));
 872       return kFALSE;
 873     }
 874   }
 875   if (nPar > 4) {
 876     if (r.Parameter(nPar-1) <= 1e-10) {
 877       if (fDebug)
 878         AliWarning(Form("Last scale %s is too small %f<1e-10",
 879                         r.ParName(nPar-1).c_str(), r.Parameter(nPar-1)));
 880       return kFALSE;
 881     }
 882   }
 883   return kTRUE;
 884 }
 885
 886
 887 //____________________________________________________________________
 888 void
 889 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
 890 {
 891   fList = DefineOutputList(dir);
 892 }
 893
 894 //____________________________________________________________________
 895 //
 896 // EOF
 897 //