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   // Default Constructor - do not use
  50   //
  51 }
  52
  53 //____________________________________________________________________
  54 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
  55   : TNamed("fmdEnergyFitter", title),
  56     fRingHistos(),
  57     fLowCut(0.3),
  58     fNParticles(3),
  59     fMinEntries(100),
  60     fFitRangeBinWidth(4),
  61     fDoFits(false),
  62     fDoMakeObject(false),
  63     fEtaAxis(0,0,0),
  64     fMaxE(10),
  65     fNEbins(300),
  66     fUseIncreasingBins(true),
  67     fMaxRelParError(.25),
  68     fMaxChi2PerNDF(20),
  69     fMinWeight(1e-7),
  70     fDebug(3)
  71 {
  72   //
  73   // Constructor
  74   //
  75   // Parameters:
  76   //    title Title of object  - not significant
  77   //
  78   fEtaAxis.SetName("etaAxis");
  79   fEtaAxis.SetTitle("#eta");
  80   fRingHistos.Add(new RingHistos(1, 'I'));
  81   fRingHistos.Add(new RingHistos(2, 'I'));
  82   fRingHistos.Add(new RingHistos(2, 'O'));
  83   fRingHistos.Add(new RingHistos(3, 'I'));
  84   fRingHistos.Add(new RingHistos(3, 'O'));
  85 }
  86
  87 //____________________________________________________________________
  88 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
  89   : TNamed(o),
  90     fRingHistos(),
  91     fLowCut(o.fLowCut),
  92     fNParticles(o.fNParticles),
  93     fMinEntries(o.fMinEntries),
  94     fFitRangeBinWidth(o.fFitRangeBinWidth),
  95     fDoFits(o.fDoFits),
  96     fDoMakeObject(o.fDoMakeObject),
  97     fEtaAxis(o.fEtaAxis),
  98     fMaxE(o.fMaxE),
  99     fNEbins(o.fNEbins),
 100     fUseIncreasingBins(o.fUseIncreasingBins),
 101     fMaxRelParError(o.fMaxRelParError),
 102     fMaxChi2PerNDF(o.fMaxChi2PerNDF),
 103     fMinWeight(o.fMinWeight),
 104     fDebug(o.fDebug)
 105 {
 106   //
 107   // Copy constructor
 108   //
 109   // Parameters:
 110   //    o Object to copy from
 111   //
 112   TIter    next(&o.fRingHistos);
 113   TObject* obj = 0;
 114   while ((obj = next())) fRingHistos.Add(obj);
 115 }
 116
 117 //____________________________________________________________________
 118 AliFMDEnergyFitter::~AliFMDEnergyFitter()
 119 {
 120   //
 121   // Destructor
 122   //
 123   fRingHistos.Delete();
 124 }
 125
 126 //____________________________________________________________________
 127 AliFMDEnergyFitter&
 128 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
 129 {
 130   //
 131   // Assignment operator
 132   //
 133   // Parameters:
 134   //    o Object to assign from
 135   //
 136   // Return:
 137   //    Reference to this
 138   //
 139   TNamed::operator=(o);
 140
 141   fLowCut        = o.fLowCut;
 142   fNParticles       = o.fNParticles;
 143   fMinEntries    = o.fMinEntries;
 144   fFitRangeBinWidth= o.fFitRangeBinWidth;
 145   fDoFits        = o.fDoFits;
 146   fDoMakeObject  = o.fDoMakeObject;
 147   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
 148   fDebug         = o.fDebug;
 149   fMaxRelParError= o.fMaxRelParError;
 150   fMaxChi2PerNDF = o.fMaxChi2PerNDF;
 151   fMinWeight     = o.fMinWeight;
 152
 153   fRingHistos.Delete();
 154   TIter    next(&o.fRingHistos);
 155   TObject* obj = 0;
 156   while ((obj = next())) fRingHistos.Add(obj);
 157
 158   return *this;
 159 }
 160
 161 //____________________________________________________________________
 162 AliFMDEnergyFitter::RingHistos*
 163 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
 164 {
 165   //
 166   // Get the ring histogram container
 167   //
 168   // Parameters:
 169   //    d Detector
 170   //    r Ring
 171   //
 172   // Return:
 173   //    Ring histogram container
 174   //
 175   Int_t idx = -1;
 176   switch (d) {
 177   case 1: idx = 0; break;
 178   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 179   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 180   }
 181   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 182
 183   return static_cast<RingHistos*>(fRingHistos.At(idx));
 184 }
 185
 186 //____________________________________________________________________
 187 void
 188 AliFMDEnergyFitter::Init(const TAxis& eAxis)
 189 {
 190   //
 191   // Initialise the task
 192   //
 193   // Parameters:
 194   //    etaAxis The eta axis to use.  Note, that if the eta axis
 195   // has already been set (using SetEtaAxis), then this parameter will be
 196   // ignored
 197   //
 198   if (fEtaAxis.GetNbins() == 0 ||
 199       TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7)
 200     SetEtaAxis(eAxis);
 201   TIter    next(&fRingHistos);
 202   RingHistos* o = 0;
 203   while ((o = static_cast<RingHistos*>(next())))
 204     o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
 205 }
 206 //____________________________________________________________________
 207 void
 208 AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
 209 {
 210   //
 211   // Set the eta axis to use.  This will force the code to use this
 212   // eta axis definition - irrespective of whatever axis is passed to
 213   // the Init member function.  Therefore, this member function can be
 214   // used to force another eta axis than one found in the correction
 215   // objects.
 216   //
 217   // Parameters:
 218   //    etaAxis Eta axis to use
 219   //
 220   SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
 221 }
 222 //____________________________________________________________________
 223 void
 224 AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
 225 {
 226   //
 227   // Set the eta axis to use.  This will force the code to use this
 228   // eta axis definition - irrespective of whatever axis is passed to
 229   // the Init member function.  Therefore, this member function can be
 230   // used to force another eta axis than one found in the correction
 231   // objects.
 232   //
 233   // Parameters:
 234   //    nBins  Number of bins
 235   //    etaMin Minimum of the eta axis
 236   //    etaMax Maximum of the eta axis
 237   //
 238   fEtaAxis.Set(nBins,etaMin,etaMax);
 239 }
 240
 241 //____________________________________________________________________
 242 Bool_t
 243 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
 244                                Bool_t           empty)
 245 {
 246   //
 247   // Fitter the input AliESDFMD object
 248   //
 249   // Parameters:
 250   //    input     Input
 251   //    empty     Whether the event is 'empty'
 252   //
 253   // Return:
 254   //    True on success, false otherwise
 255   //
 256   for(UShort_t d = 1; d <= 3; d++) {
 257     Int_t nRings = (d == 1 ? 1 : 2);
 258     for (UShort_t q = 0; q < nRings; q++) {
 259       Char_t      r    = (q == 0 ? 'I' : 'O');
 260       UShort_t    nsec = (q == 0 ?  20 :  40);
 261       UShort_t    nstr = (q == 0 ? 512 : 256);
 262
 263       RingHistos* histos = GetRingHistos(d, r);
 264
 265       for(UShort_t s = 0; s < nsec;  s++) {
 266         for(UShort_t t = 0; t < nstr; t++) {
 267           Float_t mult = input.Multiplicity(d,r,s,t);
 268
 269           // Keep dead-channel information.
 270           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0)
 271             continue;
 272
 273           // Get the pseudo-rapidity
 274           Double_t eta1 = input.Eta(d,r,s,t);
 275           Int_t ieta = fEtaAxis.FindBin(eta1);
 276           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue;
 277
 278           histos->Fill(empty, ieta-1, mult);
 279
 280         } // for strip
 281       } // for sector
 282     } // for ring
 283   } // for detector
 284
 285   return kTRUE;
 286 }
 287
 288 //____________________________________________________________________
 289 void
 290 AliFMDEnergyFitter::Fit(const TList* dir)
 291 {
 292   //
 293   // Scale the histograms to the total number of events
 294   //
 295   // Parameters:
 296   //    dir Where the histograms are
 297   //
 298   if (!fDoFits) return;
 299
 300   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 301   if (!d) return;
 302
 303   // +1          for chi^2
 304   // +3          for 1 landau
 305   // +1          for N
 306   // +fNParticles-1 for weights
 307   Int_t nStack = kN+fNParticles;
 308   THStack* stack[nStack];
 309   stack[0]         = new THStack("chi2",   "#chi^{2}/#nu");
 310   stack[kC     +1] = new THStack("c",      "Constant");
 311   stack[kDelta +1] = new THStack("delta",  "#Delta_{p}");
 312   stack[kXi    +1] = new THStack("xi",     "#xi");
 313   stack[kSigma +1] = new THStack("sigma",  "#sigma");
 314   stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}");
 315   stack[kN     +1] = new THStack("n",      "# Particles");
 316   for (Int_t i = 2; i <= fNParticles; i++)
 317     stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i));
 318   for (Int_t i = 0; i < nStack; i++)
 319     d->Add(stack[i]);
 320
 321   TIter    next(&fRingHistos);
 322   RingHistos* o = 0;
 323   while ((o = static_cast<RingHistos*>(next()))) {
 324     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
 325                           fMinEntries, fFitRangeBinWidth,
 326                           fMaxRelParError, fMaxChi2PerNDF);
 327     if (!l) continue;
 328     for (Int_t i = 0; i < l->GetEntriesFast(); i++) {
 329       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
 330     }
 331   }
 332
 333   if (!fDoMakeObject) return;
 334
 335   MakeCorrectionsObject(d);
 336 }
 337
 338 //____________________________________________________________________
 339 void
 340 AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
 341 {
 342   //
 343   // Generate the corrections object
 344   //
 345   // Parameters:
 346   //    dir List to analyse
 347   //
 348   AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
 349
 350   AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
 351   obj->SetEtaAxis(fEtaAxis);
 352   obj->SetLowCut(fLowCut);
 353
 354   TIter    next(&fRingHistos);
 355   RingHistos* o = 0;
 356   while ((o = static_cast<RingHistos*>(next()))) {
 357     o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError,
 358                     fMaxChi2PerNDF, fMinWeight);
 359   }
 360
 361   TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
 362   TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
 363   AliInfo(Form("Object %s created in output - should be extracted and copied "
 364                "to %s", oName.Data(), fileName.Data()));
 365   d->Add(obj, oName.Data());
 366 }
 367
 368 //____________________________________________________________________
 369 void
 370 AliFMDEnergyFitter::DefineOutput(TList* dir)
 371 {
 372   //
 373   // Define the output histograms.  These are put in a sub list of the
 374   // passed list.   The histograms are merged before the parent task calls
 375   // AliAnalysisTaskSE::Terminate
 376   //
 377   // Parameters:
 378   //    dir Directory to add to
 379   //
 380   TList* d = new TList;
 381   d->SetName(GetName());
 382   dir->Add(d);
 383   d->Add(&fEtaAxis);
 384   TIter    next(&fRingHistos);
 385   RingHistos* o = 0;
 386   while ((o = static_cast<RingHistos*>(next()))) {
 387     o->Output(d);
 388   }
 389 }
 390 //____________________________________________________________________
 391 void
 392 AliFMDEnergyFitter::SetDebug(Int_t dbg)
 393 {
 394   //
 395   // Set the debug level.  The higher the value the more output
 396   //
 397   // Parameters:
 398   //    dbg Debug level
 399   //
 400   fDebug = dbg;
 401   TIter    next(&fRingHistos);
 402   RingHistos* o = 0;
 403   while ((o = static_cast<RingHistos*>(next())))
 404     o->fDebug = dbg;
 405 }
 406
 407 //____________________________________________________________________
 408 void
 409 AliFMDEnergyFitter::Print(Option_t*) const
 410 {
 411   //
 412   // Print information
 413   //
 414   // Parameters:
 415   //    option Not used
 416   //
 417   char ind[gROOT->GetDirLevel()+1];
 418   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 419   ind[gROOT->GetDirLevel()] = '\0';
 420   std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n'
 421             << ind << " Low cut:                " << fLowCut << " E/E_mip\n"
 422             << ind << " Max(particles):         " << fNParticles << '\n'
 423             << ind << " Min(entries):           " << fMinEntries << '\n'
 424             << ind << " Fit range:              "
 425             << fFitRangeBinWidth << " bins\n"
 426             << ind << " Make fits:              "
 427             << (fDoFits ? "yes\n" : "no\n")
 428             << ind << " Make object:            "
 429             << (fDoMakeObject ? "yes\n" : "no\n")
 430             << ind << " Max E/E_mip:            " << fMaxE << '\n'
 431             << ind << " N bins:                 " << fNEbins << '\n'
 432             << ind << " Increasing bins:        "
 433             << (fUseIncreasingBins ?"yes\n" : "no\n")
 434             << ind << " max(delta p/p):         " << fMaxRelParError << '\n'
 435             << ind << " max(chi^2/nu):          " << fMaxChi2PerNDF << '\n'
 436             << ind << " min(a_i):               " << fMinWeight << std::endl;
 437 }
 438
 439 //====================================================================
 440 AliFMDEnergyFitter::RingHistos::RingHistos()
 441   : AliForwardUtil::RingHistos(),
 442     fEDist(0),
 443     fEmpty(0),
 444     fEtaEDists(),
 445     fList(0),
 446     fFits("AliFMDCorrELossFit::ELossFit"),
 447     fDebug(0)
 448 {
 449   //
 450   // Default CTOR
 451   //
 452 }
 453
 454 //____________________________________________________________________
 455 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
 456   : AliForwardUtil::RingHistos(d,r),
 457     fEDist(0),
 458     fEmpty(0),
 459     fEtaEDists(),
 460     fList(0),
 461     fFits("AliFMDCorrELossFit::ELossFit"),
 462     fDebug(0)
 463 {
 464   //
 465   // Constructor
 466   //
 467   // Parameters:
 468   //    d detector
 469   //    r ring
 470   //
 471   fEtaEDists.SetName("EDists");
 472 }
 473 //____________________________________________________________________
 474 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
 475   : AliForwardUtil::RingHistos(o),
 476     fEDist(o.fEDist),
 477     fEmpty(o.fEmpty),
 478     fEtaEDists(),
 479     fList(0),
 480     fFits("AliFMDCorrELossFit::ELossFit"),
 481     fDebug(0)
 482 {
 483   //
 484   // Copy constructor
 485   //
 486   // Parameters:
 487   //    o Object to copy from
 488   //
 489   fFits.Clear();
 490   TIter next(&o.fEtaEDists);
 491   TObject* obj = 0;
 492   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 493   if (o.fList) {
 494     fList = new TList;
 495     fList->SetName(fName);
 496     TIter next2(o.fList);
 497     while ((obj = next2())) fList->Add(obj->Clone());
 498   }
 499 }
 500
 501 //____________________________________________________________________
 502 AliFMDEnergyFitter::RingHistos&
 503 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
 504 {
 505   //
 506   // Assignment operator
 507   //
 508   // Parameters:
 509   //    o Object to assign from
 510   //
 511   // Return:
 512   //    Reference to this
 513   //
 514   AliForwardUtil::RingHistos::operator=(o);
 515
 516   if (fEDist) delete fEDist;
 517   if (fEmpty) delete fEmpty;
 518   fEtaEDists.Delete();
 519   if (fList) fList->Delete();
 520   fFits.Clear();
 521
 522   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
 523   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
 524
 525   TIter next(&o.fEtaEDists);
 526   TObject* obj = 0;
 527   while ((obj = next())) fEtaEDists.Add(obj->Clone());
 528
 529   if (o.fList) {
 530     fList = new TList;
 531     fList->SetName(fName);
 532     TIter next2(o.fList);
 533     while ((obj = next2())) fList->Add(obj->Clone());
 534   }
 535
 536   return *this;
 537 }
 538 //____________________________________________________________________
 539 AliFMDEnergyFitter::RingHistos::~RingHistos()
 540 {
 541   //
 542   // Destructor
 543   //
 544   if (fEDist) delete fEDist;
 545   fEtaEDists.Delete();
 546 }
 547
 548 //____________________________________________________________________
 549 void
 550 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
 551 {
 552   //
 553   // Fill histogram
 554   //
 555   // Parameters:
 556   //    empty  True if event is empty
 557   //    ieta   Eta bin
 558   //    mult   Signal
 559   //
 560   if (empty) {
 561     fEmpty->Fill(mult);
 562     return;
 563   }
 564   fEDist->Fill(mult);
 565
 566   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
 567
 568   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
 569   if (!h) return;
 570
 571   h->Fill(mult);
 572 }
 573
 574 //__________________________________________________________________
 575 TArrayD
 576 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
 577                                                    Double_t max) const
 578 {
 579   //
 580   // Make an axis with increasing bins
 581   //
 582   // Parameters:
 583   //    n    Number of bins
 584   //    min  Minimum
 585   //    max  Maximum
 586   //
 587   // Return:
 588   //    An axis with quadratically increasing bin size
 589   //
 590
 591   TArrayD bins(n+1);
 592   Double_t dx = (max-min) / n;
 593   bins[0]     = min;
 594   Int_t    i  = 1;
 595   for (i = 1; i < n+1; i++) {
 596     Double_t dI   = float(i)/n;
 597     Double_t next = bins[i-1] + dx + dI * dI * dx;
 598     bins[i]       = next;
 599     if (next > max) break;
 600   }
 601   bins.Set(i+1);
 602   return bins;
 603 }
 604
 605 //____________________________________________________________________
 606 void
 607 AliFMDEnergyFitter::RingHistos::Make(Int_t    ieta,
 608                                      Double_t emin,
 609                                      Double_t emax,
 610                                      Double_t deMax,
 611                                      Int_t    nDeBins,
 612                                      Bool_t   incr)
 613 {
 614   //
 615   // Make E/E_mip histogram
 616   //
 617   // Parameters:
 618   //    ieta    Eta bin
 619   //    eMin    Least signal
 620   //    eMax    Largest signal
 621   //    deMax   Maximum energy loss
 622   //    nDeBins Number energy loss bins
 623   //    incr    Whether to make bins of increasing size
 624   //
 625   TH1D* h = 0;
 626   TString name  = Form(fgkEDistFormat, fName.Data(), ieta);
 627   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
 628                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
 629   if (!incr)
 630     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
 631   else {
 632     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
 633     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
 634   }
 635
 636   h->SetDirectory(0);
 637   h->SetXTitle("#DeltaE/#DeltaE_{mip}");
 638   h->SetFillColor(Color());
 639   h->SetMarkerColor(Color());
 640   h->SetLineColor(Color());
 641   h->SetFillStyle(3001);
 642   h->Sumw2();
 643
 644   fEtaEDists.AddAt(h, ieta-1);
 645 }
 646 //____________________________________________________________________
 647 TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
 648                                               const char* title,
 649                                               const TAxis& eta) const
 650 {
 651   //
 652   // Make a parameter histogram
 653   //
 654   // Parameters:
 655   //    name   Name of histogram.
 656   //    title  Title of histogram.
 657   //    eta    Eta axis
 658   //
 659   // Return:
 660   //
 661   //
 662   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 663                      Form("%s for %s", title, fName.Data()),
 664                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
 665   h->SetXTitle("#eta");
 666   h->SetYTitle(title);
 667   h->SetDirectory(0);
 668   h->SetFillColor(Color());
 669   h->SetMarkerColor(Color());
 670   h->SetLineColor(Color());
 671   h->SetFillStyle(3001);
 672   h->Sumw2();
 673
 674   return h;
 675 }
 676 //____________________________________________________________________
 677 TH1D*
 678 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
 679                                           const char* title,
 680                                           const TAxis& eta,
 681                                           Int_t low,
 682                                           Int_t high,
 683                                           Double_t val,
 684                                           Double_t err) const
 685 {
 686   //
 687   // Make a histogram that contains the results of the fit over the full ring
 688   //
 689   // Parameters:
 690   //    name  Name
 691   //    title Title
 692   //    eta   Eta axis
 693   //    low   Least bin
 694   //    high  Largest bin
 695   //    val   Value of parameter
 696   //    err   Error on parameter
 697   //
 698   // Return:
 699   //    The newly allocated histogram
 700   //
 701   Double_t xlow  = eta.GetBinLowEdge(low);
 702   Double_t xhigh = eta.GetBinUpEdge(high);
 703   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
 704                      Form("%s for %s", title, fName.Data()),
 705                      1, xlow, xhigh);
 706   h->SetBinContent(1, val);
 707   h->SetBinError(1, err);
 708   h->SetXTitle("#eta");
 709   h->SetYTitle(title);
 710   h->SetDirectory(0);
 711   h->SetFillColor(0);
 712   h->SetMarkerColor(Color());
 713   h->SetLineColor(Color());
 714   h->SetFillStyle(0);
 715   h->SetLineStyle(2);
 716   h->SetLineWidth(2);
 717
 718   return h;
 719 }
 720
 721 //____________________________________________________________________
 722 void
 723 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
 724                                      Double_t maxDE, Int_t nDEbins,
 725                                      Bool_t useIncrBin)
 726 {
 727   //
 728   // Initialise object
 729   //
 730   // Parameters:
 731   //    eAxis      Eta axis
 732   //    maxDE      Max energy loss to consider
 733   //    nDEbins    Number of bins
 734   //    useIncrBin Whether to use an increasing bin size
 735   //
 736   fEDist = new TH1D(Form("%s_edist", fName.Data()),
 737                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()),
 738                     200, 0, 6);
 739   fEmpty = new TH1D(Form("%s_empty", fName.Data()),
 740                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)",
 741                          fName.Data()), 200, 0, 6);
 742   fList->Add(fEDist);
 743   fList->Add(fEmpty);
 744   // fEtaEDists.Expand(eAxis.GetNbins());
 745   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) {
 746     Double_t min = eAxis.GetBinLowEdge(i);
 747     Double_t max = eAxis.GetBinUpEdge(i);
 748     Make(i, min, max, maxDE, nDEbins, useIncrBin);
 749   }
 750   fList->Add(&fEtaEDists);
 751 }
 752 //____________________________________________________________________
 753 TObjArray*
 754 AliFMDEnergyFitter::RingHistos::Fit(TList*           dir,
 755                                     const TAxis&     eta,
 756                                     Double_t         lowCut,
 757                                     UShort_t         nParticles,
 758                                     UShort_t         minEntries,
 759                                     UShort_t         minusBins,
 760                                     Double_t         relErrorCut,
 761                                     Double_t         chi2nuCut) const
 762 {
 763   //
 764   // Fit each histogram to up to @a nParticles particle responses.
 765   //
 766   // Parameters:
 767   //    dir         Output list
 768   //    eta         Eta axis
 769   //    lowCut      Lower cut
 770   //    nParticles  Max number of convolved landaus to fit
 771   //    minEntries  Minimum number of entries
 772   //    minusBins   Number of bins from peak to subtract to
 773   //                    get the fit range
 774   //    relErrorCut Cut applied to relative error of parameter.
 775   //                    Note, for multi-particle weights, the cut
 776   //                    is loosend by a factor of 2
 777   //    chi2nuCut   Cut on @f$ \chi^2/\nu@f$ -
 778   //                    the reduced @f$\chi^2@f$
 779   //
 780
 781   // Get our ring list from the output container
 782   TList* l = GetOutputList(dir);
 783   if (!l) return 0;
 784
 785   // Get the energy distributions from the output container
 786   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
 787   if (!dists) {
 788     AliWarning(Form("Didn't find %s_EtaEDists in %s",
 789                     fName.Data(), l->GetName()));
 790     l->ls();
 791     return 0;
 792   }
 793
 794   // Container of the fit results histograms
 795   TObjArray* pars  = new TObjArray(3+nParticles+1);
 796   pars->SetName("FitResults");
 797   l->Add(pars);
 798
 799   // Result objects stored in separate list on the output
 800   TH1* hChi2   = 0;
 801   TH1* hC      = 0;
 802   TH1* hDelta  = 0;
 803   TH1* hXi     = 0;
 804   TH1* hSigma  = 0;
 805   TH1* hSigmaN = 0;
 806   TH1* hN      = 0;
 807   TH1* hA[nParticles-1];
 808   pars->Add(hChi2   = MakePar("chi2",   "#chi^{2}/#nu", eta));
 809   pars->Add(hC      = MakePar("c",      "Constant",     eta));
 810   pars->Add(hDelta  = MakePar("delta",  "#Delta_{p}",   eta));
 811   pars->Add(hXi     = MakePar("xi",     "#xi",          eta));
 812   pars->Add(hSigma  = MakePar("sigma",  "#sigma",       eta));
 813   pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}",   eta));
 814   pars->Add(hN      = MakePar("n",      "N", eta));
 815   for (UShort_t i = 1; i < nParticles; i++)
 816     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
 817
 818
 819   Int_t nDists = dists->GetEntries();
 820   Int_t low    = nDists;
 821   Int_t high   = 0;
 822   for (Int_t i = 0; i < nDists; i++) {
 823     TH1D* dist = static_cast<TH1D*>(dists->At(i));
 824     // Ignore empty histograms altoghether
 825     if (!dist || dist->GetEntries() <= 0) continue;
 826
 827     // Scale to the bin-width
 828     dist->Scale(1., "width");
 829
 830     // Normalize peak to 1
 831     Double_t max = dist->GetMaximum();
 832     if (max <= 0) continue;
 833     dist->Scale(1/max);
 834
 835     // Check that we have enough entries
 836     if (dist->GetEntries() <= minEntries) {
 837       AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
 838                       i, dist->GetName(), Int_t(dist->GetEntries()),
 839                       minEntries));
 840       continue;
 841     }
 842
 843     // Now fit
 844     TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
 845                        relErrorCut,chi2nuCut);
 846     if (!res) continue;
 847     // dist->GetListOfFunctions()->Add(res);
 848
 849     // Store eta limits
 850     low   = TMath::Min(low,i+1);
 851     high  = TMath::Max(high,i+1);
 852
 853     // Get the reduced chi square
 854     Double_t chi2 = res->GetChisquare();
 855     Int_t    ndf  = res->GetNDF();
 856
 857     // Store results of best fit in output histograms
 858     res->SetLineWidth(4);
 859     hChi2   ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
 860     hC      ->SetBinContent(i+1, res->GetParameter(kC));
 861     hDelta  ->SetBinContent(i+1, res->GetParameter(kDelta));
 862     hXi     ->SetBinContent(i+1, res->GetParameter(kXi));
 863     hSigma  ->SetBinContent(i+1, res->GetParameter(kSigma));
 864     hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));
 865     hN      ->SetBinContent(i+1, res->GetParameter(kN));
 866
 867     hC     ->SetBinError(i+1, res->GetParError(kC));
 868     hDelta ->SetBinError(i+1, res->GetParError(kDelta));
 869     hXi    ->SetBinError(i+1, res->GetParError(kXi));
 870     hSigma ->SetBinError(i+1, res->GetParError(kSigma));
 871     hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
 872     hN     ->SetBinError(i+1, res->GetParError(kN));
 873
 874     for (UShort_t j = 0; j < nParticles-1; j++) {
 875       hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
 876       hA[j]->SetBinError(i+1, res->GetParError(kA+j));
 877     }
 878   }
 879
 880   // Fit the full-ring histogram
 881   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
 882   if (total && total->GetEntries() >= minEntries) {
 883
 884     // Scale to the bin-width
 885     total->Scale(1., "width");
 886
 887     // Normalize peak to 1
 888     Double_t max = total->GetMaximum();
 889     if (max > 0) total->Scale(1/max);
 890
 891     TF1* res = FitHist(total,lowCut,nParticles,minusBins,
 892                        relErrorCut,chi2nuCut);
 893     if (res) {
 894       // Make histograms for the result of this fit
 895       Double_t chi2 = res->GetChisquare();
 896       Int_t    ndf  = res->GetNDF();
 897       pars->Add(MakeTotal("t_chi2",   "#chi^{2}/#nu", eta, low, high,
 898                           ndf > 0 ? chi2/ndf : 0, 0));
 899       pars->Add(MakeTotal("t_c",      "Constant",     eta, low, high,
 900                           res->GetParameter(kC),
 901                           res->GetParError(kC)));
 902       pars->Add(MakeTotal("t_delta",  "#Delta_{p}",   eta, low, high,
 903                           res->GetParameter(kDelta),
 904                           res->GetParError(kDelta)));
 905       pars->Add(MakeTotal("t_xi",     "#xi",          eta, low, high,
 906                           res->GetParameter(kXi),
 907                           res->GetParError(kXi)));
 908       pars->Add(MakeTotal("t_sigma",  "#sigma",       eta, low, high,
 909                           res->GetParameter(kSigma),
 910                           res->GetParError(kSigma)));
 911       pars->Add(MakeTotal("t_sigman", "#sigma_{n}",   eta, low, high,
 912                           res->GetParameter(kSigmaN),
 913                           res->GetParError(kSigmaN)));
 914       pars->Add(MakeTotal("t_n",    "N",              eta, low, high,
 915                           res->GetParameter(kN),0));
 916       for (UShort_t j = 0; j < nParticles-1; j++)
 917         pars->Add(MakeTotal(Form("t_a%d",j+2),
 918                             Form("a_{%d}",j+2), eta, low, high,
 919                             res->GetParameter(kA+j),
 920                             res->GetParError(kA+j)));
 921     }
 922   }
 923
 924   // Clean up list of histogram.  Histograms with no entries or
 925   // no functions are deleted.  We have to do this using the TObjLink
 926   // objects stored in the list since ROOT cannot guaranty the validity
 927   // of iterators when removing from a list - tsck.  Should just implement
 928   // TIter::Remove().
 929   TObjLink* lnk = dists->FirstLink();
 930   while (lnk) {
 931     TH1* h = static_cast<TH1*>(lnk->GetObject());
 932     bool remove = false;
 933     if (h->GetEntries() <= 0) {
 934       // AliWarning(Form("No entries in %s - removing", h->GetName()));
 935       remove = true;
 936     }
 937     else if (h->GetListOfFunctions()->GetEntries() <= 0) {
 938       // AliWarning(Form("No fuctions associated with %s - removing",
 939       //            h->GetName()));
 940       remove = true;
 941     }
 942     if (remove) {
 943       TObjLink* keep = lnk->Next();
 944       dists->Remove(lnk);
 945       lnk = keep;
 946       continue;
 947     }
 948     lnk = lnk->Next();
 949   }
 950
 951   return pars;
 952 }
 953
 954 //____________________________________________________________________
 955 TF1*
 956 AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
 957                                         Double_t lowCut,
 958                                         UShort_t nParticles,
 959                                         UShort_t minusBins,
 960                                         Double_t relErrorCut,
 961                                         Double_t chi2nuCut) const
 962 {
 963   //
 964   // Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
 965   // maximum bin content in the range @f$ [E_{min},\infty]@f$ is
 966   // found.  Then the fit range is set to the bin range
 967   // @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
 968   // particle signal is fitted to that.  The parameters of that fit
 969   // is then used as seeds for a fit of the @f$ N@f$ particle response
 970   // to the data in the range
 971   // @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
 972   //
 973   // Parameters:
 974   //    dist        Histogram to fit
 975   //    lowCut      Lower cut @f$ E_{min}@f$ on signal
 976   //    nParticles  Max number @f$ N@f$ of convolved landaus to fit
 977   //    minusBins   Number of bins @f$ \Delta b@f$ from peak to
 978   //                    subtract to get the fit range
 979   //    relErrorCut Cut applied to relative error of parameter.
 980   //                    Note, for multi-particle weights, the cut
 981   //                    is loosend by a factor of 2
 982   //    chi2nuCut   Cut on @f$ \chi^2/\nu@f$ -
 983   //                    the reduced @f$\chi^2@f$
 984   //
 985   // Return:
 986   //    The best fit function
 987   //
 988   Double_t maxRange = 10;
 989
 990   // Create a fitter object
 991   AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins);
 992   f.Clear();
 993
 994
 995   // If we are only asked to fit a single particle, return this fit,
 996   // no matter what.
 997   if (nParticles == 1) {
 998     TF1* r = f.Fit1Particle(dist, 0);
 999     if (!r) return 0;
1000     return new TF1(*r);
1001   }
1002
1003   // Fit from 2 upto n particles
1004   for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
1005
1006
1007   // Now, we need to select the best fit
1008   Int_t nFits = f.fFitResults.GetEntriesFast();
1009   TF1*  good[nFits];
1010   for (Int_t i = nFits-1; i >= 0; i--) {
1011     good[i] = 0;
1012     TF1* ff = static_cast<TF1*>(f.fFunctions.At(i));
1013     // ff->SetLineColor(Color());
1014     ff->SetRange(0, maxRange);
1015     dist->GetListOfFunctions()->Add(new TF1(*ff));
1016     if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i)),
1017                     relErrorCut, chi2nuCut)) {
1018       good[i] = ff;
1019       ff->SetLineWidth(2);
1020       // f.fFitResults.At(i)->Print("V");
1021     }
1022   }
1023   // If all else fails, use the 1 particle fit
1024   TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
1025
1026   // Find the fit with the most valid particles
1027   for (Int_t i = nFits-1; i > 0; i--) {
1028     if (!good[i]) continue;
1029     if (fDebug > 1) {
1030       AliInfo(Form("Choosing fit with n=%d", i+1));
1031       f.fFitResults.At(i)->Print();
1032     }
1033     ret = good[i];
1034     break;
1035   }
1036   // Give a warning if we're using fall-back
1037   if (ret == f.fFunctions.At(0)) {
1038     AliWarning("Choosing fall-back 1 particle fit");
1039   }
1040   // Copy our result and return (the functions are owned by the fitter)
1041   TF1* fret = new TF1(*ret);
1042   return fret;
1043 }
1044
1045 //____________________________________________________________________
1046 Bool_t
1047 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
1048                                             Double_t relErrorCut,
1049                                             Double_t chi2nuCut) const
1050 {
1051   //
1052   // Check the result of the fit. Returns true if
1053   // - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
1054   // - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note,
1055   //   for multi-particle fits, this requirement is relaxed by a
1056   //   factor of 2
1057   // - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$
1058   //   particle response
1059   //
1060   // Parameters:
1061   //    r           Result to check
1062   //    relErrorCut Cut @f$ \delta_e@f$ applied to relative error
1063   //                    of parameter.
1064   //    chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$
1065   //
1066   // Return:
1067   //    true if fit is good.
1068   //
1069   if (fDebug > 10) r->Print();
1070   TString  n    = r->GetName();
1071   n.ReplaceAll("TFitResult-", "");
1072   Double_t chi2 = r->Chi2();
1073   Int_t    ndf  = r->Ndf();
1074   // Double_t prob = r.Prob();
1075   Bool_t ret = kTRUE;
1076
1077   // Check that the reduced chi square isn't larger than 5
1078   if (ndf <= 0 || chi2 / ndf > chi2nuCut) {
1079     if (fDebug > 2) {
1080       AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f",
1081                       n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
1082                       chi2nuCut)); }
1083     ret = kFALSE;
1084   }
1085
1086   // Check each parameter
1087   UShort_t nPar = r->NPar();
1088   for (UShort_t i = 0; i < nPar; i++) {
1089     if (i == kN) continue;  // Skip the number parameter
1090
1091     // Get value and error and check value
1092     Double_t v  = r->Parameter(i);
1093     Double_t e  = r->ParError(i);
1094     if (v == 0) continue;
1095
1096     // Calculate the relative error and check it
1097     Double_t re  = e / v;
1098     Double_t cut = relErrorCut * (i < kN ? 1 : 2);
1099     if (re > cut) {
1100       if (fDebug > 2) {
1101         AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
1102                         n.Data(), r->ParName(i).c_str(),
1103                         r->ParName(i).c_str(), e, v,
1104                         100*(v == 0 ? 0 : e/v),
1105                         100*(cut))); }
1106       ret = kFALSE;
1107     }
1108   }
1109
1110   // Check if we have scale parameters
1111   if (nPar > kN) {
1112
1113     // Check that the last particle has a significant contribution
1114     Double_t lastScale = r->Parameter(nPar-1);
1115     if (lastScale <= 1e-7) {
1116       if (fDebug) {
1117         AliWarning(Form("%s: %s=%9.6f<1e-7",
1118                         n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
1119       ret = kFALSE;
1120     }
1121   }
1122   return ret;
1123 }
1124
1125
1126 //__________________________________________________________________
1127 void
1128 AliFMDEnergyFitter::RingHistos::FindBestFits(TList*              d,
1129                                              AliFMDCorrELossFit& obj,
1130                                              const TAxis&        eta,
1131                                              Double_t            relErrorCut,
1132                                              Double_t            chi2nuCut,
1133                                              Double_t            minWeightCut)
1134 {
1135   //
1136   // Find the best fits
1137   //
1138   // Parameters:
1139   //    d              Parent list
1140   //    obj            Object to add fits to
1141   //    eta            Eta axis
1142   //    relErrorCut    Cut applied to relative error of parameter.
1143   //                       Note, for multi-particle weights, the cut
1144   //                       is loosend by a factor of 2
1145   //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
1146   //                       the reduced @f$\chi^2@f$
1147   //    minWeightCut   Least valid @f$ a_i@f$
1148   //
1149
1150   // Get our ring list from the output container
1151   TList* l = GetOutputList(d);
1152   if (!l) return;
1153
1154   // Get the energy distributions from the output container
1155   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
1156   if (!dists) {
1157     AliWarning(Form("Didn't find %s_EtaEDists in %s",
1158                     fName.Data(), l->GetName()));
1159     l->ls();
1160     return;
1161   }
1162   Int_t nBin = eta.GetNbins();
1163
1164   for (Int_t b = 1; b <= nBin; b++) {
1165     TString n(Form(fgkEDistFormat, fName.Data(), b));
1166     TH1D*   dist = static_cast<TH1D*>(dists->FindObject(n));
1167     // Ignore empty histograms altoghether
1168     if (!dist || dist->GetEntries() <= 0) continue;
1169
1170     AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
1171                                                      relErrorCut,
1172                                                      chi2nuCut,
1173                                                      minWeightCut);
1174     best->fDet  = fDet;
1175     best->fRing = fRing;
1176     best->fBin  = b; //
1177     // Double_t eta = fAxis->GetBinCenter(b);
1178     obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
1179   }
1180 }
1181
1182 //__________________________________________________________________
1183 AliFMDCorrELossFit::ELossFit*
1184 AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist,
1185                                             Double_t relErrorCut,
1186                                             Double_t chi2nuCut,
1187                                             Double_t minWeightCut)
1188 {
1189   //
1190   // Find the best fit
1191   //
1192   // Parameters:
1193   //    dist           Histogram
1194   //    relErrorCut    Cut applied to relative error of parameter.
1195   //                       Note, for multi-particle weights, the cut
1196   //                       is loosend by a factor of 2
1197   //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
1198   //                       the reduced @f$\chi^2@f$
1199   //    minWeightCut   Least valid @f$ a_i@f$
1200   //
1201   // Return:
1202   //    Best fit
1203   //
1204   TList* funcs = dist->GetListOfFunctions();
1205   TIter  next(funcs);
1206   TF1*   func = 0;
1207   fFits.Clear();
1208   Int_t  i = 0;
1209   // Info("FindBestFit", "%s", dist->GetName());
1210   while ((func = static_cast<TF1*>(next()))) {
1211     AliFMDCorrELossFit::ELossFit* fit =
1212       new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
1213     fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
1214   }
1215   fFits.Sort(false);
1216   // fFits.Print();
1217   return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
1218 }
1219
1220
1221 //____________________________________________________________________
1222 void
1223 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
1224 {
1225   //
1226   // Define outputs
1227   //
1228   // Parameters:
1229   //    dir
1230   //
1231   fList = DefineOutputList(dir);
1232 }
1233
1234 //____________________________________________________________________
1235 //
1236 // EOF
1237 //