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