]> 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 (Ruben)
[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()));
f7cfc454 398 d->Add(new TNamed("filename", fileName.Data()));
0bd4b00f 399 d->Add(obj, oName.Data());
f8715167 400}
401
402//____________________________________________________________________
403void
404AliFMDEnergyFitter::DefineOutput(TList* dir)
405{
7984e5f7 406 //
407 // Define the output histograms. These are put in a sub list of the
408 // passed list. The histograms are merged before the parent task calls
409 // AliAnalysisTaskSE::Terminate
410 //
411 // Parameters:
412 // dir Directory to add to
413 //
f8715167 414 TList* d = new TList;
415 d->SetName(GetName());
416 dir->Add(d);
417 d->Add(&fEtaAxis);
418 TIter next(&fRingHistos);
419 RingHistos* o = 0;
420 while ((o = static_cast<RingHistos*>(next()))) {
421 o->Output(d);
422 }
423}
424//____________________________________________________________________
425void
426AliFMDEnergyFitter::SetDebug(Int_t dbg)
427{
7984e5f7 428 //
429 // Set the debug level. The higher the value the more output
430 //
431 // Parameters:
432 // dbg Debug level
433 //
f8715167 434 fDebug = dbg;
435 TIter next(&fRingHistos);
436 RingHistos* o = 0;
437 while ((o = static_cast<RingHistos*>(next())))
438 o->fDebug = dbg;
439}
0bd4b00f 440
441//____________________________________________________________________
442void
443AliFMDEnergyFitter::Print(Option_t*) const
444{
7984e5f7 445 //
446 // Print information
447 //
448 // Parameters:
449 // option Not used
450 //
0bd4b00f 451 char ind[gROOT->GetDirLevel()+1];
452 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
453 ind[gROOT->GetDirLevel()] = '\0';
1f480471 454 std::cout << ind << ClassName() << ": " << GetName() << '\n'
0bd4b00f 455 << ind << " Low cut: " << fLowCut << " E/E_mip\n"
456 << ind << " Max(particles): " << fNParticles << '\n'
457 << ind << " Min(entries): " << fMinEntries << '\n'
458 << ind << " Fit range: "
459 << fFitRangeBinWidth << " bins\n"
460 << ind << " Make fits: "
461 << (fDoFits ? "yes\n" : "no\n")
462 << ind << " Make object: "
463 << (fDoMakeObject ? "yes\n" : "no\n")
464 << ind << " Max E/E_mip: " << fMaxE << '\n'
465 << ind << " N bins: " << fNEbins << '\n'
466 << ind << " Increasing bins: "
467 << (fUseIncreasingBins ?"yes\n" : "no\n")
468 << ind << " max(delta p/p): " << fMaxRelParError << '\n'
469 << ind << " max(chi^2/nu): " << fMaxChi2PerNDF << '\n'
470 << ind << " min(a_i): " << fMinWeight << std::endl;
471}
f8715167 472
473//====================================================================
474AliFMDEnergyFitter::RingHistos::RingHistos()
475 : AliForwardUtil::RingHistos(),
476 fEDist(0),
477 fEmpty(0),
478 fEtaEDists(),
479 fList(0),
0bd4b00f 480 fFits("AliFMDCorrELossFit::ELossFit"),
f8715167 481 fDebug(0)
7984e5f7 482{
483 //
484 // Default CTOR
485 //
486}
f8715167 487
488//____________________________________________________________________
489AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
490 : AliForwardUtil::RingHistos(d,r),
491 fEDist(0),
492 fEmpty(0),
493 fEtaEDists(),
494 fList(0),
0bd4b00f 495 fFits("AliFMDCorrELossFit::ELossFit"),
f8715167 496 fDebug(0)
497{
7984e5f7 498 //
499 // Constructor
500 //
501 // Parameters:
502 // d detector
503 // r ring
504 //
f8715167 505 fEtaEDists.SetName("EDists");
506}
507//____________________________________________________________________
508AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
509 : AliForwardUtil::RingHistos(o),
510 fEDist(o.fEDist),
511 fEmpty(o.fEmpty),
512 fEtaEDists(),
513 fList(0),
0bd4b00f 514 fFits("AliFMDCorrELossFit::ELossFit"),
f8715167 515 fDebug(0)
516{
7984e5f7 517 //
518 // Copy constructor
519 //
520 // Parameters:
521 // o Object to copy from
522 //
0bd4b00f 523 fFits.Clear();
f8715167 524 TIter next(&o.fEtaEDists);
525 TObject* obj = 0;
526 while ((obj = next())) fEtaEDists.Add(obj->Clone());
527 if (o.fList) {
528 fList = new TList;
529 fList->SetName(fName);
530 TIter next2(o.fList);
531 while ((obj = next2())) fList->Add(obj->Clone());
532 }
533}
534
535//____________________________________________________________________
536AliFMDEnergyFitter::RingHistos&
537AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
538{
7984e5f7 539 //
540 // Assignment operator
541 //
542 // Parameters:
543 // o Object to assign from
544 //
545 // Return:
546 // Reference to this
547 //
f8715167 548 AliForwardUtil::RingHistos::operator=(o);
549
550 if (fEDist) delete fEDist;
551 if (fEmpty) delete fEmpty;
552 fEtaEDists.Delete();
553 if (fList) fList->Delete();
0bd4b00f 554 fFits.Clear();
f8715167 555
556 fEDist = static_cast<TH1D*>(o.fEDist->Clone());
557 fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
558
559 TIter next(&o.fEtaEDists);
560 TObject* obj = 0;
561 while ((obj = next())) fEtaEDists.Add(obj->Clone());
562
563 if (o.fList) {
564 fList = new TList;
565 fList->SetName(fName);
566 TIter next2(o.fList);
567 while ((obj = next2())) fList->Add(obj->Clone());
568 }
569
570 return *this;
571}
572//____________________________________________________________________
573AliFMDEnergyFitter::RingHistos::~RingHistos()
574{
7984e5f7 575 //
576 // Destructor
577 //
f8715167 578 if (fEDist) delete fEDist;
579 fEtaEDists.Delete();
580}
581
582//____________________________________________________________________
583void
5e4d905e 584AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta,
585 Int_t /* icent */, Double_t mult)
f8715167 586{
7984e5f7 587 //
588 // Fill histogram
589 //
590 // Parameters:
591 // empty True if event is empty
5e4d905e 592 // ieta Eta bin (0-based)
593 // icent Centrality bin (1-based)
7984e5f7 594 // mult Signal
595 //
f8715167 596 if (empty) {
597 fEmpty->Fill(mult);
598 return;
599 }
600 fEDist->Fill(mult);
601
5e4d905e 602 // Here, we should first look up in a centrality array, and then in
603 // that array look up the eta bin - or perhaps we should do it the
604 // other way around?
f8715167 605 if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
606
607 TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
608 if (!h) return;
609
610 h->Fill(mult);
611}
612
613//__________________________________________________________________
614TArrayD
66b34429 615AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
616 Double_t max) const
f8715167 617{
7984e5f7 618 //
619 // Make an axis with increasing bins
620 //
621 // Parameters:
622 // n Number of bins
623 // min Minimum
624 // max Maximum
625 //
626 // Return:
627 // An axis with quadratically increasing bin size
628 //
629
f8715167 630 TArrayD bins(n+1);
631 Double_t dx = (max-min) / n;
632 bins[0] = min;
633 Int_t i = 1;
634 for (i = 1; i < n+1; i++) {
635 Double_t dI = float(i)/n;
66b34429 636 Double_t next = bins[i-1] + dx + dI * dI * dx;
f8715167 637 bins[i] = next;
638 if (next > max) break;
639 }
640 bins.Set(i+1);
641 return bins;
642}
643
644//____________________________________________________________________
645void
c389303e 646AliFMDEnergyFitter::RingHistos::Make(Int_t ieta,
647 Double_t emin,
648 Double_t emax,
649 Double_t deMax,
650 Int_t nDeBins,
651 Bool_t incr)
f8715167 652{
7984e5f7 653 //
654 // Make E/E_mip histogram
655 //
656 // Parameters:
657 // ieta Eta bin
658 // eMin Least signal
659 // eMax Largest signal
660 // deMax Maximum energy loss
661 // nDeBins Number energy loss bins
662 // incr Whether to make bins of increasing size
663 //
f8715167 664 TH1D* h = 0;
0bd4b00f 665 TString name = Form(fgkEDistFormat, fName.Data(), ieta);
f8715167 666 TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
667 "(#eta bin %d)", fName.Data(), emin, emax, ieta);
668 if (!incr)
669 h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
670 else {
671 TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
672 h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
673 }
674
675 h->SetDirectory(0);
676 h->SetXTitle("#DeltaE/#DeltaE_{mip}");
677 h->SetFillColor(Color());
678 h->SetMarkerColor(Color());
679 h->SetLineColor(Color());
680 h->SetFillStyle(3001);
681 h->Sumw2();
682
683 fEtaEDists.AddAt(h, ieta-1);
684}
685//____________________________________________________________________
c389303e 686TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
687 const char* title,
688 const TAxis& eta) const
f8715167 689{
7984e5f7 690 //
691 // Make a parameter histogram
692 //
693 // Parameters:
694 // name Name of histogram.
695 // title Title of histogram.
696 // eta Eta axis
697 //
698 // Return:
699 //
700 //
f8715167 701 TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
702 Form("%s for %s", title, fName.Data()),
703 eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
704 h->SetXTitle("#eta");
705 h->SetYTitle(title);
706 h->SetDirectory(0);
707 h->SetFillColor(Color());
708 h->SetMarkerColor(Color());
709 h->SetLineColor(Color());
710 h->SetFillStyle(3001);
711 h->Sumw2();
712
713 return h;
714}
715//____________________________________________________________________
716TH1D*
717AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
718 const char* title,
719 const TAxis& eta,
720 Int_t low,
721 Int_t high,
722 Double_t val,
723 Double_t err) const
724{
7984e5f7 725 //
726 // Make a histogram that contains the results of the fit over the full ring
727 //
728 // Parameters:
729 // name Name
730 // title Title
731 // eta Eta axis
732 // low Least bin
733 // high Largest bin
734 // val Value of parameter
735 // err Error on parameter
736 //
737 // Return:
738 // The newly allocated histogram
739 //
f8715167 740 Double_t xlow = eta.GetBinLowEdge(low);
741 Double_t xhigh = eta.GetBinUpEdge(high);
742 TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
743 Form("%s for %s", title, fName.Data()),
744 1, xlow, xhigh);
745 h->SetBinContent(1, val);
746 h->SetBinError(1, err);
747 h->SetXTitle("#eta");
748 h->SetYTitle(title);
749 h->SetDirectory(0);
750 h->SetFillColor(0);
751 h->SetMarkerColor(Color());
752 h->SetLineColor(Color());
753 h->SetFillStyle(0);
754 h->SetLineStyle(2);
755 h->SetLineWidth(2);
756
757 return h;
758}
759
760//____________________________________________________________________
761void
762AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
5e4d905e 763 const TAxis& /* cAxis */,
f8715167 764 Double_t maxDE, Int_t nDEbins,
765 Bool_t useIncrBin)
766{
7984e5f7 767 //
768 // Initialise object
769 //
770 // Parameters:
771 // eAxis Eta axis
772 // maxDE Max energy loss to consider
773 // nDEbins Number of bins
774 // useIncrBin Whether to use an increasing bin size
775 //
f8715167 776 fEDist = new TH1D(Form("%s_edist", fName.Data()),
777 Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()),
778 200, 0, 6);
779 fEmpty = new TH1D(Form("%s_empty", fName.Data()),
780 Form("#DeltaE/#DeltaE_{mip} for %s (empty events)",
781 fName.Data()), 200, 0, 6);
782 fList->Add(fEDist);
783 fList->Add(fEmpty);
5e4d905e 784 // Here, we should make an array of centrality bins ...
785 // In that case, the structure will be
786 //
787 // -+- Ring1 -+- Centrality1 -+- Eta1
788 // | | +- Eta2
789 // ... ... ...
790 // | | +- EtaM
791 // | +- Centrality2 -+- Eta1
792 // ... ... ...
793 // | | +- EtaM
794 // ... ...
795 // | +- CentralityN -+- Eta1
796 // ... ... ...
797 // | +- EtaM
798 // -+- Ring2 -+- Centrality1 -+- Eta1
799 // | | +- Eta2
800 // ... ... ...
801 // | | +- EtaM
802 // | +- Centrality2 -+- Eta1
803 // ... ... ...
804 // | | +- EtaM
805 // ... ...
806 // | +- CentralityN -+- Eta1
807 // ... ... ...
808 // | +- EtaM
809 // ... ... ...
810 // -+- RingO -+- Centrality1 -+- Eta1
811 // | +- Eta2
812 // ... ...
813 // | +- EtaM
814 // +- Centrality2 -+- Eta1
815 // ... ...
816 // | +- EtaM
817 // ...
818 // +- CentralityN -+- Eta1
819 // ... ...
820 // +- EtaM
821 //
f8715167 822 // fEtaEDists.Expand(eAxis.GetNbins());
823 for (Int_t i = 1; i <= eAxis.GetNbins(); i++) {
824 Double_t min = eAxis.GetBinLowEdge(i);
825 Double_t max = eAxis.GetBinUpEdge(i);
5e4d905e 826 // Or may we should do it here? In that case we'd have the
827 // following structure:
828 // Ring1 -+- Eta1 -+- Centrality1
829 // | +- Centrality2
830 // ... ...
831 // | +- CentralityN
832 // +- Eta2 -+- Centrality1
833 // | +- Centrality2
834 // ... ...
835 // | +- CentralityN
836 // ...
837 // +- EtaM -+- Centrality1
838 // +- Centrality2
839 // ...
840 // +- CentralityN
f8715167 841 Make(i, min, max, maxDE, nDEbins, useIncrBin);
842 }
843 fList->Add(&fEtaEDists);
844}
845//____________________________________________________________________
846TObjArray*
0bd4b00f 847AliFMDEnergyFitter::RingHistos::Fit(TList* dir,
848 const TAxis& eta,
849 Double_t lowCut,
850 UShort_t nParticles,
851 UShort_t minEntries,
852 UShort_t minusBins,
853 Double_t relErrorCut,
854 Double_t chi2nuCut) const
f8715167 855{
7984e5f7 856 //
857 // Fit each histogram to up to @a nParticles particle responses.
858 //
859 // Parameters:
860 // dir Output list
861 // eta Eta axis
862 // lowCut Lower cut
863 // nParticles Max number of convolved landaus to fit
864 // minEntries Minimum number of entries
865 // minusBins Number of bins from peak to subtract to
866 // get the fit range
867 // relErrorCut Cut applied to relative error of parameter.
868 // Note, for multi-particle weights, the cut
869 // is loosend by a factor of 2
870 // chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
871 // the reduced @f$\chi^2@f$
872 //
873
c389303e 874 // Get our ring list from the output container
f8715167 875 TList* l = GetOutputList(dir);
876 if (!l) return 0;
877
c389303e 878 // Get the energy distributions from the output container
f8715167 879 TList* dists = static_cast<TList*>(l->FindObject("EDists"));
880 if (!dists) {
881 AliWarning(Form("Didn't find %s_EtaEDists in %s",
882 fName.Data(), l->GetName()));
883 l->ls();
884 return 0;
885 }
886
c389303e 887 // Container of the fit results histograms
888 TObjArray* pars = new TObjArray(3+nParticles+1);
f8715167 889 pars->SetName("FitResults");
890 l->Add(pars);
891
c389303e 892 // Result objects stored in separate list on the output
893 TH1* hChi2 = 0;
894 TH1* hC = 0;
895 TH1* hDelta = 0;
896 TH1* hXi = 0;
897 TH1* hSigma = 0;
898 TH1* hSigmaN = 0;
899 TH1* hN = 0;
900 TH1* hA[nParticles-1];
901 pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta));
902 pars->Add(hC = MakePar("c", "Constant", eta));
903 pars->Add(hDelta = MakePar("delta", "#Delta_{p}", eta));
904 pars->Add(hXi = MakePar("xi", "#xi", eta));
905 pars->Add(hSigma = MakePar("sigma", "#sigma", eta));
906 pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}", eta));
907 pars->Add(hN = MakePar("n", "N", eta));
908 for (UShort_t i = 1; i < nParticles; i++)
f8715167 909 pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
910
911
912 Int_t nDists = dists->GetEntries();
913 Int_t low = nDists;
914 Int_t high = 0;
915 for (Int_t i = 0; i < nDists; i++) {
916 TH1D* dist = static_cast<TH1D*>(dists->At(i));
c389303e 917 // Ignore empty histograms altoghether
918 if (!dist || dist->GetEntries() <= 0) continue;
f8715167 919
c389303e 920 // Scale to the bin-width
921 dist->Scale(1., "width");
922
923 // Normalize peak to 1
924 Double_t max = dist->GetMaximum();
0bd4b00f 925 if (max <= 0) continue;
c389303e 926 dist->Scale(1/max);
927
928 // Check that we have enough entries
0bd4b00f 929 if (dist->GetEntries() <= minEntries) {
930 AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
931 i, dist->GetName(), Int_t(dist->GetEntries()),
932 minEntries));
933 continue;
934 }
f8715167 935
c389303e 936 // Now fit
937 TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
938 relErrorCut,chi2nuCut);
f8715167 939 if (!res) continue;
c389303e 940 // dist->GetListOfFunctions()->Add(res);
0bd4b00f 941
c389303e 942 // Store eta limits
f8715167 943 low = TMath::Min(low,i+1);
944 high = TMath::Max(high,i+1);
945
c389303e 946 // Get the reduced chi square
f8715167 947 Double_t chi2 = res->GetChisquare();
948 Int_t ndf = res->GetNDF();
c389303e 949
950 // Store results of best fit in output histograms
951 res->SetLineWidth(4);
952 hChi2 ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
953 hC ->SetBinContent(i+1, res->GetParameter(kC));
954 hDelta ->SetBinContent(i+1, res->GetParameter(kDelta));
955 hXi ->SetBinContent(i+1, res->GetParameter(kXi));
956 hSigma ->SetBinContent(i+1, res->GetParameter(kSigma));
957 hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));
958 hN ->SetBinContent(i+1, res->GetParameter(kN));
959
960 hC ->SetBinError(i+1, res->GetParError(kC));
961 hDelta ->SetBinError(i+1, res->GetParError(kDelta));
962 hXi ->SetBinError(i+1, res->GetParError(kXi));
963 hSigma ->SetBinError(i+1, res->GetParError(kSigma));
964 hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
965 hN ->SetBinError(i+1, res->GetParError(kN));
966
967 for (UShort_t j = 0; j < nParticles-1; j++) {
968 hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
969 hA[j]->SetBinError(i+1, res->GetParError(kA+j));
f8715167 970 }
971 }
972
c389303e 973 // Fit the full-ring histogram
f8715167 974 TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
4b9857f3 975 if (total && total->GetEntries() >= minEntries) {
0bd4b00f 976
977 // Scale to the bin-width
978 total->Scale(1., "width");
979
980 // Normalize peak to 1
981 Double_t max = total->GetMaximum();
982 if (max > 0) total->Scale(1/max);
983
c389303e 984 TF1* res = FitHist(total,lowCut,nParticles,minusBins,
985 relErrorCut,chi2nuCut);
f8715167 986 if (res) {
c389303e 987 // Make histograms for the result of this fit
f8715167 988 Double_t chi2 = res->GetChisquare();
989 Int_t ndf = res->GetNDF();
c389303e 990 pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
f8715167 991 ndf > 0 ? chi2/ndf : 0, 0));
c389303e 992 pars->Add(MakeTotal("t_c", "Constant", eta, low, high,
993 res->GetParameter(kC),
994 res->GetParError(kC)));
995 pars->Add(MakeTotal("t_delta", "#Delta_{p}", eta, low, high,
996 res->GetParameter(kDelta),
997 res->GetParError(kDelta)));
998 pars->Add(MakeTotal("t_xi", "#xi", eta, low, high,
999 res->GetParameter(kXi),
1000 res->GetParError(kXi)));
1001 pars->Add(MakeTotal("t_sigma", "#sigma", eta, low, high,
1002 res->GetParameter(kSigma),
1003 res->GetParError(kSigma)));
1004 pars->Add(MakeTotal("t_sigman", "#sigma_{n}", eta, low, high,
1005 res->GetParameter(kSigmaN),
1006 res->GetParError(kSigmaN)));
1007 pars->Add(MakeTotal("t_n", "N", eta, low, high,
1008 res->GetParameter(kN),0));
1009 for (UShort_t j = 0; j < nParticles-1; j++)
1010 pars->Add(MakeTotal(Form("t_a%d",j+2),
1011 Form("a_{%d}",j+2), eta, low, high,
1012 res->GetParameter(kA+j),
1013 res->GetParError(kA+j)));
f8715167 1014 }
1015 }
1016
c389303e 1017 // Clean up list of histogram. Histograms with no entries or
1018 // no functions are deleted. We have to do this using the TObjLink
1019 // objects stored in the list since ROOT cannot guaranty the validity
1020 // of iterators when removing from a list - tsck. Should just implement
1021 // TIter::Remove().
f8715167 1022 TObjLink* lnk = dists->FirstLink();
1023 while (lnk) {
1024 TH1* h = static_cast<TH1*>(lnk->GetObject());
0bd4b00f 1025 bool remove = false;
1026 if (h->GetEntries() <= 0) {
1027 // AliWarning(Form("No entries in %s - removing", h->GetName()));
1028 remove = true;
1029 }
1030 else if (h->GetListOfFunctions()->GetEntries() <= 0) {
1031 // AliWarning(Form("No fuctions associated with %s - removing",
1032 // h->GetName()));
1033 remove = true;
1034 }
1035 if (remove) {
f8715167 1036 TObjLink* keep = lnk->Next();
1037 dists->Remove(lnk);
1038 lnk = keep;
1039 continue;
1040 }
1041 lnk = lnk->Next();
1042 }
1043
1044 return pars;
1045}
1046
1047//____________________________________________________________________
1048TF1*
4b9857f3 1049AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,
1050 Double_t lowCut,
c389303e 1051 UShort_t nParticles,
1052 UShort_t minusBins,
1053 Double_t relErrorCut,
1054 Double_t chi2nuCut) const
f8715167 1055{
7984e5f7 1056 //
1057 // Fit a signal histogram. First, the bin @f$ b_{min}@f$ with
1058 // maximum bin content in the range @f$ [E_{min},\infty]@f$ is
1059 // found. Then the fit range is set to the bin range
1060 // @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
1061 // particle signal is fitted to that. The parameters of that fit
1062 // is then used as seeds for a fit of the @f$ N@f$ particle response
1063 // to the data in the range
1064 // @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
1065 //
1066 // Parameters:
1067 // dist Histogram to fit
1068 // lowCut Lower cut @f$ E_{min}@f$ on signal
1069 // nParticles Max number @f$ N@f$ of convolved landaus to fit
1070 // minusBins Number of bins @f$ \Delta b@f$ from peak to
1071 // subtract to get the fit range
1072 // relErrorCut Cut applied to relative error of parameter.
1073 // Note, for multi-particle weights, the cut
1074 // is loosend by a factor of 2
1075 // chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
1076 // the reduced @f$\chi^2@f$
1077 //
1078 // Return:
1079 // The best fit function
1080 //
f8715167 1081 Double_t maxRange = 10;
1082
c389303e 1083 // Create a fitter object
4b9857f3 1084 AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins);
1085 f.Clear();
f8715167 1086
c389303e 1087
4b9857f3 1088 // If we are only asked to fit a single particle, return this fit,
1089 // no matter what.
c389303e 1090 if (nParticles == 1) {
4b9857f3 1091 TF1* r = f.Fit1Particle(dist, 0);
1092 if (!r) return 0;
1093 return new TF1(*r);
f8715167 1094 }
f8715167 1095
4b9857f3 1096 // Fit from 2 upto n particles
c389303e 1097 for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
f8715167 1098
f8715167 1099
4b9857f3 1100 // Now, we need to select the best fit
fb3430ac 1101 Int_t nFits = f.GetFitResults().GetEntriesFast();
4b9857f3 1102 TF1* good[nFits];
1103 for (Int_t i = nFits-1; i >= 0; i--) {
c389303e 1104 good[i] = 0;
fb3430ac 1105 TF1* ff = static_cast<TF1*>(f.GetFunctions().At(i));
c389303e 1106 // ff->SetLineColor(Color());
1107 ff->SetRange(0, maxRange);
1108 dist->GetListOfFunctions()->Add(new TF1(*ff));
fb3430ac 1109 if (CheckResult(static_cast<TFitResult*>(f.GetFitResults().At(i)),
c389303e 1110 relErrorCut, chi2nuCut)) {
1111 good[i] = ff;
1112 ff->SetLineWidth(2);
1113 // f.fFitResults.At(i)->Print("V");
1114 }
4b9857f3 1115 }
1116 // If all else fails, use the 1 particle fit
fb3430ac 1117 TF1* ret = static_cast<TF1*>(f.GetFunctions().At(0));
c389303e 1118
1119 // Find the fit with the most valid particles
4b9857f3 1120 for (Int_t i = nFits-1; i > 0; i--) {
1121 if (!good[i]) continue;
c389303e 1122 if (fDebug > 1) {
1123 AliInfo(Form("Choosing fit with n=%d", i+1));
fb3430ac 1124 f.GetFitResults().At(i)->Print();
c389303e 1125 }
4b9857f3 1126 ret = good[i];
f8715167 1127 break;
1128 }
c389303e 1129 // Give a warning if we're using fall-back
fb3430ac 1130 if (ret == f.GetFunctions().At(0)) {
c389303e 1131 AliWarning("Choosing fall-back 1 particle fit");
0bd4b00f 1132 }
c389303e 1133 // Copy our result and return (the functions are owned by the fitter)
1134 TF1* fret = new TF1(*ret);
1135 return fret;
f8715167 1136}
1137
1138//____________________________________________________________________
1139Bool_t
c389303e 1140AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
1141 Double_t relErrorCut,
1142 Double_t chi2nuCut) const
f8715167 1143{
7984e5f7 1144 //
1145 // Check the result of the fit. Returns true if
1146 // - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
1147 // - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters. Note,
1148 // for multi-particle fits, this requirement is relaxed by a
1149 // factor of 2
1150 // - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$
1151 // particle response
1152 //
1153 // Parameters:
1154 // r Result to check
1155 // relErrorCut Cut @f$ \delta_e@f$ applied to relative error
1156 // of parameter.
1157 // chi2nuCut Cut @f$ \max{\chi^2/\nu}@f$
1158 //
1159 // Return:
1160 // true if fit is good.
1161 //
c389303e 1162 if (fDebug > 10) r->Print();
1163 TString n = r->GetName();
1164 n.ReplaceAll("TFitResult-", "");
4b9857f3 1165 Double_t chi2 = r->Chi2();
1166 Int_t ndf = r->Ndf();
f8715167 1167 // Double_t prob = r.Prob();
c389303e 1168 Bool_t ret = kTRUE;
1169
1170 // Check that the reduced chi square isn't larger than 5
1171 if (ndf <= 0 || chi2 / ndf > chi2nuCut) {
0bd4b00f 1172 if (fDebug > 2) {
c389303e 1173 AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f",
1174 n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
0bd4b00f 1175 chi2nuCut)); }
c389303e 1176 ret = kFALSE;
f8715167 1177 }
1178
c389303e 1179 // Check each parameter
4b9857f3 1180 UShort_t nPar = r->NPar();
f8715167 1181 for (UShort_t i = 0; i < nPar; i++) {
c389303e 1182 if (i == kN) continue; // Skip the number parameter
f8715167 1183
c389303e 1184 // Get value and error and check value
1185 Double_t v = r->Parameter(i);
1186 Double_t e = r->ParError(i);
f8715167 1187 if (v == 0) continue;
c389303e 1188
1189 // Calculate the relative error and check it
1190 Double_t re = e / v;
1191 Double_t cut = relErrorCut * (i < kN ? 1 : 2);
1192 if (re > cut) {
0bd4b00f 1193 if (fDebug > 2) {
c389303e 1194 AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
1195 n.Data(), r->ParName(i).c_str(),
1196 r->ParName(i).c_str(), e, v,
1197 100*(v == 0 ? 0 : e/v),
0bd4b00f 1198 100*(cut))); }
c389303e 1199 ret = kFALSE;
f8715167 1200 }
1201 }
c389303e 1202
1203 // Check if we have scale parameters
1204 if (nPar > kN) {
1205
1206 // Check that the last particle has a significant contribution
4b9857f3 1207 Double_t lastScale = r->Parameter(nPar-1);
1208 if (lastScale <= 1e-7) {
0bd4b00f 1209 if (fDebug) {
c389303e 1210 AliWarning(Form("%s: %s=%9.6f<1e-7",
0bd4b00f 1211 n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
c389303e 1212 ret = kFALSE;
f8715167 1213 }
1214 }
c389303e 1215 return ret;
f8715167 1216}
1217
1218
0bd4b00f 1219//__________________________________________________________________
1220void
fb3430ac 1221AliFMDEnergyFitter::RingHistos::FindBestFits(const TList* d,
0bd4b00f 1222 AliFMDCorrELossFit& obj,
1223 const TAxis& eta,
1224 Double_t relErrorCut,
1225 Double_t chi2nuCut,
1226 Double_t minWeightCut)
1227{
7984e5f7 1228 //
1229 // Find the best fits
1230 //
1231 // Parameters:
1232 // d Parent list
1233 // obj Object to add fits to
1234 // eta Eta axis
1235 // relErrorCut Cut applied to relative error of parameter.
1236 // Note, for multi-particle weights, the cut
1237 // is loosend by a factor of 2
1238 // chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
1239 // the reduced @f$\chi^2@f$
1240 // minWeightCut Least valid @f$ a_i@f$
1241 //
1242
0bd4b00f 1243 // Get our ring list from the output container
1244 TList* l = GetOutputList(d);
1245 if (!l) return;
1246
1247 // Get the energy distributions from the output container
1248 TList* dists = static_cast<TList*>(l->FindObject("EDists"));
1249 if (!dists) {
1250 AliWarning(Form("Didn't find %s_EtaEDists in %s",
1251 fName.Data(), l->GetName()));
1252 l->ls();
1253 return;
1254 }
1255 Int_t nBin = eta.GetNbins();
1256
1257 for (Int_t b = 1; b <= nBin; b++) {
1258 TString n(Form(fgkEDistFormat, fName.Data(), b));
1259 TH1D* dist = static_cast<TH1D*>(dists->FindObject(n));
1260 // Ignore empty histograms altoghether
1261 if (!dist || dist->GetEntries() <= 0) continue;
1262
1263 AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
1264 relErrorCut,
1265 chi2nuCut,
1266 minWeightCut);
1267 best->fDet = fDet;
1268 best->fRing = fRing;
1269 best->fBin = b; //
1270 // Double_t eta = fAxis->GetBinCenter(b);
1271 obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
1272 }
1273}
1274
1275//__________________________________________________________________
1276AliFMDCorrELossFit::ELossFit*
fb3430ac 1277AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist,
0bd4b00f 1278 Double_t relErrorCut,
1279 Double_t chi2nuCut,
1280 Double_t minWeightCut)
1281{
7984e5f7 1282 //
1283 // Find the best fit
1284 //
1285 // Parameters:
1286 // dist Histogram
1287 // relErrorCut Cut applied to relative error of parameter.
1288 // Note, for multi-particle weights, the cut
1289 // is loosend by a factor of 2
1290 // chi2nuCut Cut on @f$ \chi^2/\nu@f$ -
1291 // the reduced @f$\chi^2@f$
1292 // minWeightCut Least valid @f$ a_i@f$
1293 //
1294 // Return:
1295 // Best fit
1296 //
0bd4b00f 1297 TList* funcs = dist->GetListOfFunctions();
1298 TIter next(funcs);
1299 TF1* func = 0;
1300 fFits.Clear();
1301 Int_t i = 0;
1302 // Info("FindBestFit", "%s", dist->GetName());
1303 while ((func = static_cast<TF1*>(next()))) {
1304 AliFMDCorrELossFit::ELossFit* fit =
1305 new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
1306 fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
1307 }
1308 fFits.Sort(false);
1309 // fFits.Print();
1310 return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
1311}
1312
1313
f8715167 1314//____________________________________________________________________
1315void
1316AliFMDEnergyFitter::RingHistos::Output(TList* dir)
1317{
7984e5f7 1318 //
1319 // Define outputs
1320 //
1321 // Parameters:
1322 // dir
1323 //
f8715167 1324 fList = DefineOutputList(dir);
1325}
1326
1327//____________________________________________________________________
1328//
1329// EOF
1330//