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