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