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