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