]> 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
Renamed
[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(),
33 fLowCut(0.3),
c389303e 34 fNParticles(3),
f8715167 35 fMinEntries(100),
c389303e 36 fFitRangeBinWidth(4),
f8715167 37 fDoFits(false),
0bd4b00f 38 fDoMakeObject(false),
f8715167 39 fEtaAxis(),
4b9857f3 40 fMaxE(10),
41 fNEbins(300),
42 fUseIncreasingBins(true),
c389303e 43 fMaxRelParError(.25),
44 fMaxChi2PerNDF(20),
0bd4b00f 45 fMinWeight(1e-7),
f8715167 46 fDebug(0)
47{}
48
49//____________________________________________________________________
50AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
51 : TNamed("fmdEnergyFitter", title),
52 fRingHistos(),
53 fLowCut(0.3),
c389303e 54 fNParticles(3),
f8715167 55 fMinEntries(100),
c389303e 56 fFitRangeBinWidth(4),
f8715167 57 fDoFits(false),
0bd4b00f 58 fDoMakeObject(false),
4b9857f3 59 fEtaAxis(0,0,0),
60 fMaxE(10),
61 fNEbins(300),
62 fUseIncreasingBins(true),
c389303e 63 fMaxRelParError(.25),
0bd4b00f 64 fMaxChi2PerNDF(20),
65 fMinWeight(1e-7),
f8715167 66 fDebug(3)
67{
68 fEtaAxis.SetName("etaAxis");
69 fEtaAxis.SetTitle("#eta");
70 fRingHistos.Add(new RingHistos(1, 'I'));
71 fRingHistos.Add(new RingHistos(2, 'I'));
72 fRingHistos.Add(new RingHistos(2, 'O'));
73 fRingHistos.Add(new RingHistos(3, 'I'));
74 fRingHistos.Add(new RingHistos(3, 'O'));
75}
76
77//____________________________________________________________________
78AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
79 : TNamed(o),
80 fRingHistos(),
81 fLowCut(o.fLowCut),
c389303e 82 fNParticles(o.fNParticles),
f8715167 83 fMinEntries(o.fMinEntries),
c389303e 84 fFitRangeBinWidth(o.fFitRangeBinWidth),
f8715167 85 fDoFits(o.fDoFits),
0bd4b00f 86 fDoMakeObject(o.fDoMakeObject),
f8715167 87 fEtaAxis(o.fEtaAxis),
4b9857f3 88 fMaxE(o.fMaxE),
89 fNEbins(o.fNEbins),
90 fUseIncreasingBins(o.fUseIncreasingBins),
c389303e 91 fMaxRelParError(o.fMaxRelParError),
92 fMaxChi2PerNDF(o.fMaxChi2PerNDF),
0bd4b00f 93 fMinWeight(o.fMinWeight),
f8715167 94 fDebug(o.fDebug)
95{
96 TIter next(&o.fRingHistos);
97 TObject* obj = 0;
98 while ((obj = next())) fRingHistos.Add(obj);
99}
100
101//____________________________________________________________________
102AliFMDEnergyFitter::~AliFMDEnergyFitter()
103{
104 fRingHistos.Delete();
105}
106
107//____________________________________________________________________
108AliFMDEnergyFitter&
109AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
110{
111 TNamed::operator=(o);
112
113 fLowCut = o.fLowCut;
c389303e 114 fNParticles = o.fNParticles;
f8715167 115 fMinEntries = o.fMinEntries;
c389303e 116 fFitRangeBinWidth= o.fFitRangeBinWidth;
f8715167 117 fDoFits = o.fDoFits;
0bd4b00f 118 fDoMakeObject = o.fDoMakeObject;
f8715167 119 fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
120 fDebug = o.fDebug;
0bd4b00f 121 fMaxRelParError= o.fMaxRelParError;
122 fMaxChi2PerNDF = o.fMaxChi2PerNDF;
123 fMinWeight = o.fMinWeight;
f8715167 124
125 fRingHistos.Delete();
126 TIter next(&o.fRingHistos);
127 TObject* obj = 0;
128 while ((obj = next())) fRingHistos.Add(obj);
129
130 return *this;
131}
132
133//____________________________________________________________________
134AliFMDEnergyFitter::RingHistos*
135AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
136{
137 Int_t idx = -1;
138 switch (d) {
139 case 1: idx = 0; break;
140 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
141 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
142 }
143 if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
144
145 return static_cast<RingHistos*>(fRingHistos.At(idx));
146}
147
148//____________________________________________________________________
149void
150AliFMDEnergyFitter::Init(const TAxis& eAxis)
151{
4b9857f3 152 if (fEtaAxis.GetNbins() == 0 ||
7c1a1f1d 153 TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7)
4b9857f3 154 SetEtaAxis(eAxis);
f8715167 155 TIter next(&fRingHistos);
156 RingHistos* o = 0;
157 while ((o = static_cast<RingHistos*>(next())))
4b9857f3 158 o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
f8715167 159}
4b9857f3 160//____________________________________________________________________
161void
162AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
163{
164 SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
165}
166//____________________________________________________________________
167void
168AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
169{
170 fEtaAxis.Set(nBins,etaMin,etaMax);
171}
172
f8715167 173//____________________________________________________________________
174Bool_t
175AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
176 Bool_t empty)
177{
178 for(UShort_t d = 1; d <= 3; d++) {
179 Int_t nRings = (d == 1 ? 1 : 2);
180 for (UShort_t q = 0; q < nRings; q++) {
181 Char_t r = (q == 0 ? 'I' : 'O');
182 UShort_t nsec = (q == 0 ? 20 : 40);
183 UShort_t nstr = (q == 0 ? 512 : 256);
184
185 RingHistos* histos = GetRingHistos(d, r);
186
187 for(UShort_t s = 0; s < nsec; s++) {
188 for(UShort_t t = 0; t < nstr; t++) {
189 Float_t mult = input.Multiplicity(d,r,s,t);
190
191 // Keep dead-channel information.
192 if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0)
193 continue;
194
195 // Get the pseudo-rapidity
196 Double_t eta1 = input.Eta(d,r,s,t);
197 Int_t ieta = fEtaAxis.FindBin(eta1);
198 if (ieta < 1 || ieta > fEtaAxis.GetNbins()) continue;
199
200 histos->Fill(empty, ieta-1, mult);
201
202 } // for strip
203 } // for sector
204 } // for ring
205 } // for detector
206
207 return kTRUE;
208}
209
210//____________________________________________________________________
211void
7c1a1f1d 212AliFMDEnergyFitter::Fit(const TList* dir)
f8715167 213{
214 if (!fDoFits) return;
215
216 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
217 if (!d) return;
218
219 // +1 for chi^2
220 // +3 for 1 landau
221 // +1 for N
c389303e 222 // +fNParticles-1 for weights
223 Int_t nStack = kN+fNParticles;
f8715167 224 THStack* stack[nStack];
c389303e 225 stack[0] = new THStack("chi2", "#chi^{2}/#nu");
226 stack[kC +1] = new THStack("c", "Constant");
227 stack[kDelta +1] = new THStack("delta", "#Delta_{p}");
228 stack[kXi +1] = new THStack("xi", "#xi");
229 stack[kSigma +1] = new THStack("sigma", "#sigma");
230 stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}");
231 stack[kN +1] = new THStack("n", "# Particles");
232 for (Int_t i = 2; i <= fNParticles; i++)
233 stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i));
f8715167 234 for (Int_t i = 0; i < nStack; i++)
235 d->Add(stack[i]);
236
237 TIter next(&fRingHistos);
238 RingHistos* o = 0;
239 while ((o = static_cast<RingHistos*>(next()))) {
c389303e 240 TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
241 fMinEntries, fFitRangeBinWidth,
242 fMaxRelParError, fMaxChi2PerNDF);
f8715167 243 if (!l) continue;
244 for (Int_t i = 0; i < l->GetEntriesFast(); i++) {
245 stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
246 }
247 }
0bd4b00f 248
249 if (!fDoMakeObject) return;
250
251 MakeCorrectionsObject(d);
252}
253
254//____________________________________________________________________
255void
256AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
257{
258 AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
259
260 AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
261 obj->SetEtaAxis(fEtaAxis);
262 obj->SetLowCut(fLowCut);
263
264 TIter next(&fRingHistos);
265 RingHistos* o = 0;
266 while ((o = static_cast<RingHistos*>(next()))) {
267 o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError,
268 fMaxChi2PerNDF, fMinWeight);
269 }
270
271 TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
272 TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
273 AliInfo(Form("Object %s created in output - should be extracted and copied "
274 "to %s", oName.Data(), fileName.Data()));
275 d->Add(obj, oName.Data());
f8715167 276}
277
278//____________________________________________________________________
279void
280AliFMDEnergyFitter::DefineOutput(TList* dir)
281{
282 TList* d = new TList;
283 d->SetName(GetName());
284 dir->Add(d);
285 d->Add(&fEtaAxis);
286 TIter next(&fRingHistos);
287 RingHistos* o = 0;
288 while ((o = static_cast<RingHistos*>(next()))) {
289 o->Output(d);
290 }
291}
292//____________________________________________________________________
293void
294AliFMDEnergyFitter::SetDebug(Int_t dbg)
295{
296 fDebug = dbg;
297 TIter next(&fRingHistos);
298 RingHistos* o = 0;
299 while ((o = static_cast<RingHistos*>(next())))
300 o->fDebug = dbg;
301}
0bd4b00f 302
303//____________________________________________________________________
304void
305AliFMDEnergyFitter::Print(Option_t*) const
306{
307 char ind[gROOT->GetDirLevel()+1];
308 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
309 ind[gROOT->GetDirLevel()] = '\0';
310 std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n'
311 << ind << " Low cut: " << fLowCut << " E/E_mip\n"
312 << ind << " Max(particles): " << fNParticles << '\n'
313 << ind << " Min(entries): " << fMinEntries << '\n'
314 << ind << " Fit range: "
315 << fFitRangeBinWidth << " bins\n"
316 << ind << " Make fits: "
317 << (fDoFits ? "yes\n" : "no\n")
318 << ind << " Make object: "
319 << (fDoMakeObject ? "yes\n" : "no\n")
320 << ind << " Max E/E_mip: " << fMaxE << '\n'
321 << ind << " N bins: " << fNEbins << '\n'
322 << ind << " Increasing bins: "
323 << (fUseIncreasingBins ?"yes\n" : "no\n")
324 << ind << " max(delta p/p): " << fMaxRelParError << '\n'
325 << ind << " max(chi^2/nu): " << fMaxChi2PerNDF << '\n'
326 << ind << " min(a_i): " << fMinWeight << std::endl;
327}
f8715167 328
329//====================================================================
330AliFMDEnergyFitter::RingHistos::RingHistos()
331 : AliForwardUtil::RingHistos(),
332 fEDist(0),
333 fEmpty(0),
334 fEtaEDists(),
335 fList(0),
0bd4b00f 336 fFits("AliFMDCorrELossFit::ELossFit"),
f8715167 337 fDebug(0)
338{}
339
340//____________________________________________________________________
341AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
342 : AliForwardUtil::RingHistos(d,r),
343 fEDist(0),
344 fEmpty(0),
345 fEtaEDists(),
346 fList(0),
0bd4b00f 347 fFits("AliFMDCorrELossFit::ELossFit"),
f8715167 348 fDebug(0)
349{
350 fEtaEDists.SetName("EDists");
351}
352//____________________________________________________________________
353AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
354 : AliForwardUtil::RingHistos(o),
355 fEDist(o.fEDist),
356 fEmpty(o.fEmpty),
357 fEtaEDists(),
358 fList(0),
0bd4b00f 359 fFits("AliFMDCorrELossFit::ELossFit"),
f8715167 360 fDebug(0)
361{
0bd4b00f 362 fFits.Clear();
f8715167 363 TIter next(&o.fEtaEDists);
364 TObject* obj = 0;
365 while ((obj = next())) fEtaEDists.Add(obj->Clone());
366 if (o.fList) {
367 fList = new TList;
368 fList->SetName(fName);
369 TIter next2(o.fList);
370 while ((obj = next2())) fList->Add(obj->Clone());
371 }
372}
373
374//____________________________________________________________________
375AliFMDEnergyFitter::RingHistos&
376AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
377{
378 AliForwardUtil::RingHistos::operator=(o);
379
380 if (fEDist) delete fEDist;
381 if (fEmpty) delete fEmpty;
382 fEtaEDists.Delete();
383 if (fList) fList->Delete();
0bd4b00f 384 fFits.Clear();
f8715167 385
386 fEDist = static_cast<TH1D*>(o.fEDist->Clone());
387 fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
388
389 TIter next(&o.fEtaEDists);
390 TObject* obj = 0;
391 while ((obj = next())) fEtaEDists.Add(obj->Clone());
392
393 if (o.fList) {
394 fList = new TList;
395 fList->SetName(fName);
396 TIter next2(o.fList);
397 while ((obj = next2())) fList->Add(obj->Clone());
398 }
399
400 return *this;
401}
402//____________________________________________________________________
403AliFMDEnergyFitter::RingHistos::~RingHistos()
404{
405 if (fEDist) delete fEDist;
406 fEtaEDists.Delete();
407}
408
409//____________________________________________________________________
410void
411AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
412{
413 if (empty) {
414 fEmpty->Fill(mult);
415 return;
416 }
417 fEDist->Fill(mult);
418
419 if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
420
421 TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
422 if (!h) return;
423
424 h->Fill(mult);
425}
426
427//__________________________________________________________________
428TArrayD
66b34429 429AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
430 Double_t max) const
f8715167 431{
432 // Service function to define a logarithmic axis.
433 // Parameters:
434 // n Number of bins
435 // min Minimum of axis
436 // max Maximum of axis
437 TArrayD bins(n+1);
438 Double_t dx = (max-min) / n;
439 bins[0] = min;
440 Int_t i = 1;
441 for (i = 1; i < n+1; i++) {
442 Double_t dI = float(i)/n;
66b34429 443 Double_t next = bins[i-1] + dx + dI * dI * dx;
f8715167 444 bins[i] = next;
445 if (next > max) break;
446 }
447 bins.Set(i+1);
448 return bins;
449}
450
451//____________________________________________________________________
452void
c389303e 453AliFMDEnergyFitter::RingHistos::Make(Int_t ieta,
454 Double_t emin,
455 Double_t emax,
456 Double_t deMax,
457 Int_t nDeBins,
458 Bool_t incr)
f8715167 459{
460 TH1D* h = 0;
0bd4b00f 461 TString name = Form(fgkEDistFormat, fName.Data(), ieta);
f8715167 462 TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
463 "(#eta bin %d)", fName.Data(), emin, emax, ieta);
464 if (!incr)
465 h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
466 else {
467 TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
468 h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
469 }
470
471 h->SetDirectory(0);
472 h->SetXTitle("#DeltaE/#DeltaE_{mip}");
473 h->SetFillColor(Color());
474 h->SetMarkerColor(Color());
475 h->SetLineColor(Color());
476 h->SetFillStyle(3001);
477 h->Sumw2();
478
479 fEtaEDists.AddAt(h, ieta-1);
480}
481//____________________________________________________________________
c389303e 482TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
483 const char* title,
484 const TAxis& eta) const
f8715167 485{
486 TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
487 Form("%s for %s", title, fName.Data()),
488 eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
489 h->SetXTitle("#eta");
490 h->SetYTitle(title);
491 h->SetDirectory(0);
492 h->SetFillColor(Color());
493 h->SetMarkerColor(Color());
494 h->SetLineColor(Color());
495 h->SetFillStyle(3001);
496 h->Sumw2();
497
498 return h;
499}
500//____________________________________________________________________
501TH1D*
502AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
503 const char* title,
504 const TAxis& eta,
505 Int_t low,
506 Int_t high,
507 Double_t val,
508 Double_t err) const
509{
510 Double_t xlow = eta.GetBinLowEdge(low);
511 Double_t xhigh = eta.GetBinUpEdge(high);
512 TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name),
513 Form("%s for %s", title, fName.Data()),
514 1, xlow, xhigh);
515 h->SetBinContent(1, val);
516 h->SetBinError(1, err);
517 h->SetXTitle("#eta");
518 h->SetYTitle(title);
519 h->SetDirectory(0);
520 h->SetFillColor(0);
521 h->SetMarkerColor(Color());
522 h->SetLineColor(Color());
523 h->SetFillStyle(0);
524 h->SetLineStyle(2);
525 h->SetLineWidth(2);
526
527 return h;
528}
529
530//____________________________________________________________________
531void
532AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
533 Double_t maxDE, Int_t nDEbins,
534 Bool_t useIncrBin)
535{
536 fEDist = new TH1D(Form("%s_edist", fName.Data()),
537 Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()),
538 200, 0, 6);
539 fEmpty = new TH1D(Form("%s_empty", fName.Data()),
540 Form("#DeltaE/#DeltaE_{mip} for %s (empty events)",
541 fName.Data()), 200, 0, 6);
542 fList->Add(fEDist);
543 fList->Add(fEmpty);
544 // fEtaEDists.Expand(eAxis.GetNbins());
545 for (Int_t i = 1; i <= eAxis.GetNbins(); i++) {
546 Double_t min = eAxis.GetBinLowEdge(i);
547 Double_t max = eAxis.GetBinUpEdge(i);
548 Make(i, min, max, maxDE, nDEbins, useIncrBin);
549 }
550 fList->Add(&fEtaEDists);
551}
552//____________________________________________________________________
553TObjArray*
0bd4b00f 554AliFMDEnergyFitter::RingHistos::Fit(TList* dir,
555 const TAxis& eta,
556 Double_t lowCut,
557 UShort_t nParticles,
558 UShort_t minEntries,
559 UShort_t minusBins,
560 Double_t relErrorCut,
561 Double_t chi2nuCut) const
f8715167 562{
c389303e 563 // Get our ring list from the output container
f8715167 564 TList* l = GetOutputList(dir);
565 if (!l) return 0;
566
c389303e 567 // Get the energy distributions from the output container
f8715167 568 TList* dists = static_cast<TList*>(l->FindObject("EDists"));
569 if (!dists) {
570 AliWarning(Form("Didn't find %s_EtaEDists in %s",
571 fName.Data(), l->GetName()));
572 l->ls();
573 return 0;
574 }
575
c389303e 576 // Container of the fit results histograms
577 TObjArray* pars = new TObjArray(3+nParticles+1);
f8715167 578 pars->SetName("FitResults");
579 l->Add(pars);
580
c389303e 581 // Result objects stored in separate list on the output
582 TH1* hChi2 = 0;
583 TH1* hC = 0;
584 TH1* hDelta = 0;
585 TH1* hXi = 0;
586 TH1* hSigma = 0;
587 TH1* hSigmaN = 0;
588 TH1* hN = 0;
589 TH1* hA[nParticles-1];
590 pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta));
591 pars->Add(hC = MakePar("c", "Constant", eta));
592 pars->Add(hDelta = MakePar("delta", "#Delta_{p}", eta));
593 pars->Add(hXi = MakePar("xi", "#xi", eta));
594 pars->Add(hSigma = MakePar("sigma", "#sigma", eta));
595 pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}", eta));
596 pars->Add(hN = MakePar("n", "N", eta));
597 for (UShort_t i = 1; i < nParticles; i++)
f8715167 598 pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
599
600
601 Int_t nDists = dists->GetEntries();
602 Int_t low = nDists;
603 Int_t high = 0;
604 for (Int_t i = 0; i < nDists; i++) {
605 TH1D* dist = static_cast<TH1D*>(dists->At(i));
c389303e 606 // Ignore empty histograms altoghether
607 if (!dist || dist->GetEntries() <= 0) continue;
f8715167 608
c389303e 609 // Scale to the bin-width
610 dist->Scale(1., "width");
611
612 // Normalize peak to 1
613 Double_t max = dist->GetMaximum();
0bd4b00f 614 if (max <= 0) continue;
c389303e 615 dist->Scale(1/max);
616
617 // Check that we have enough entries
0bd4b00f 618 if (dist->GetEntries() <= minEntries) {
619 AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
620 i, dist->GetName(), Int_t(dist->GetEntries()),
621 minEntries));
622 continue;
623 }
f8715167 624
c389303e 625 // Now fit
626 TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
627 relErrorCut,chi2nuCut);
f8715167 628 if (!res) continue;
c389303e 629 // dist->GetListOfFunctions()->Add(res);
0bd4b00f 630
c389303e 631 // Store eta limits
f8715167 632 low = TMath::Min(low,i+1);
633 high = TMath::Max(high,i+1);
634
c389303e 635 // Get the reduced chi square
f8715167 636 Double_t chi2 = res->GetChisquare();
637 Int_t ndf = res->GetNDF();
c389303e 638
639 // Store results of best fit in output histograms
640 res->SetLineWidth(4);
641 hChi2 ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
642 hC ->SetBinContent(i+1, res->GetParameter(kC));
643 hDelta ->SetBinContent(i+1, res->GetParameter(kDelta));
644 hXi ->SetBinContent(i+1, res->GetParameter(kXi));
645 hSigma ->SetBinContent(i+1, res->GetParameter(kSigma));
646 hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));
647 hN ->SetBinContent(i+1, res->GetParameter(kN));
648
649 hC ->SetBinError(i+1, res->GetParError(kC));
650 hDelta ->SetBinError(i+1, res->GetParError(kDelta));
651 hXi ->SetBinError(i+1, res->GetParError(kXi));
652 hSigma ->SetBinError(i+1, res->GetParError(kSigma));
653 hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
654 hN ->SetBinError(i+1, res->GetParError(kN));
655
656 for (UShort_t j = 0; j < nParticles-1; j++) {
657 hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
658 hA[j]->SetBinError(i+1, res->GetParError(kA+j));
f8715167 659 }
660 }
661
c389303e 662 // Fit the full-ring histogram
f8715167 663 TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
4b9857f3 664 if (total && total->GetEntries() >= minEntries) {
0bd4b00f 665
666 // Scale to the bin-width
667 total->Scale(1., "width");
668
669 // Normalize peak to 1
670 Double_t max = total->GetMaximum();
671 if (max > 0) total->Scale(1/max);
672
c389303e 673 TF1* res = FitHist(total,lowCut,nParticles,minusBins,
674 relErrorCut,chi2nuCut);
f8715167 675 if (res) {
c389303e 676 // Make histograms for the result of this fit
f8715167 677 Double_t chi2 = res->GetChisquare();
678 Int_t ndf = res->GetNDF();
c389303e 679 pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
f8715167 680 ndf > 0 ? chi2/ndf : 0, 0));
c389303e 681 pars->Add(MakeTotal("t_c", "Constant", eta, low, high,
682 res->GetParameter(kC),
683 res->GetParError(kC)));
684 pars->Add(MakeTotal("t_delta", "#Delta_{p}", eta, low, high,
685 res->GetParameter(kDelta),
686 res->GetParError(kDelta)));
687 pars->Add(MakeTotal("t_xi", "#xi", eta, low, high,
688 res->GetParameter(kXi),
689 res->GetParError(kXi)));
690 pars->Add(MakeTotal("t_sigma", "#sigma", eta, low, high,
691 res->GetParameter(kSigma),
692 res->GetParError(kSigma)));
693 pars->Add(MakeTotal("t_sigman", "#sigma_{n}", eta, low, high,
694 res->GetParameter(kSigmaN),
695 res->GetParError(kSigmaN)));
696 pars->Add(MakeTotal("t_n", "N", eta, low, high,
697 res->GetParameter(kN),0));
698 for (UShort_t j = 0; j < nParticles-1; j++)
699 pars->Add(MakeTotal(Form("t_a%d",j+2),
700 Form("a_{%d}",j+2), eta, low, high,
701 res->GetParameter(kA+j),
702 res->GetParError(kA+j)));
f8715167 703 }
704 }
705
c389303e 706 // Clean up list of histogram. Histograms with no entries or
707 // no functions are deleted. We have to do this using the TObjLink
708 // objects stored in the list since ROOT cannot guaranty the validity
709 // of iterators when removing from a list - tsck. Should just implement
710 // TIter::Remove().
f8715167 711 TObjLink* lnk = dists->FirstLink();
712 while (lnk) {
713 TH1* h = static_cast<TH1*>(lnk->GetObject());
0bd4b00f 714 bool remove = false;
715 if (h->GetEntries() <= 0) {
716 // AliWarning(Form("No entries in %s - removing", h->GetName()));
717 remove = true;
718 }
719 else if (h->GetListOfFunctions()->GetEntries() <= 0) {
720 // AliWarning(Form("No fuctions associated with %s - removing",
721 // h->GetName()));
722 remove = true;
723 }
724 if (remove) {
f8715167 725 TObjLink* keep = lnk->Next();
726 dists->Remove(lnk);
727 lnk = keep;
728 continue;
729 }
730 lnk = lnk->Next();
731 }
732
733 return pars;
734}
735
736//____________________________________________________________________
737TF1*
4b9857f3 738AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,
739 Double_t lowCut,
c389303e 740 UShort_t nParticles,
741 UShort_t minusBins,
742 Double_t relErrorCut,
743 Double_t chi2nuCut) const
f8715167 744{
745 Double_t maxRange = 10;
746
c389303e 747 // Create a fitter object
4b9857f3 748 AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins);
749 f.Clear();
f8715167 750
c389303e 751
4b9857f3 752 // If we are only asked to fit a single particle, return this fit,
753 // no matter what.
c389303e 754 if (nParticles == 1) {
4b9857f3 755 TF1* r = f.Fit1Particle(dist, 0);
756 if (!r) return 0;
757 return new TF1(*r);
f8715167 758 }
f8715167 759
4b9857f3 760 // Fit from 2 upto n particles
c389303e 761 for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
f8715167 762
f8715167 763
4b9857f3 764 // Now, we need to select the best fit
765 Int_t nFits = f.fFitResults.GetEntriesFast();
766 TF1* good[nFits];
767 for (Int_t i = nFits-1; i >= 0; i--) {
c389303e 768 good[i] = 0;
769 TF1* ff = static_cast<TF1*>(f.fFunctions.At(i));
770 // ff->SetLineColor(Color());
771 ff->SetRange(0, maxRange);
772 dist->GetListOfFunctions()->Add(new TF1(*ff));
773 if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i)),
774 relErrorCut, chi2nuCut)) {
775 good[i] = ff;
776 ff->SetLineWidth(2);
777 // f.fFitResults.At(i)->Print("V");
778 }
4b9857f3 779 }
780 // If all else fails, use the 1 particle fit
781 TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
c389303e 782
783 // Find the fit with the most valid particles
4b9857f3 784 for (Int_t i = nFits-1; i > 0; i--) {
785 if (!good[i]) continue;
c389303e 786 if (fDebug > 1) {
787 AliInfo(Form("Choosing fit with n=%d", i+1));
788 f.fFitResults.At(i)->Print();
789 }
4b9857f3 790 ret = good[i];
f8715167 791 break;
792 }
c389303e 793 // Give a warning if we're using fall-back
0bd4b00f 794 if (ret == f.fFunctions.At(0)) {
c389303e 795 AliWarning("Choosing fall-back 1 particle fit");
0bd4b00f 796 }
c389303e 797 // Copy our result and return (the functions are owned by the fitter)
798 TF1* fret = new TF1(*ret);
799 return fret;
f8715167 800}
801
802//____________________________________________________________________
803Bool_t
c389303e 804AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
805 Double_t relErrorCut,
806 Double_t chi2nuCut) const
f8715167 807{
c389303e 808 if (fDebug > 10) r->Print();
809 TString n = r->GetName();
810 n.ReplaceAll("TFitResult-", "");
4b9857f3 811 Double_t chi2 = r->Chi2();
812 Int_t ndf = r->Ndf();
f8715167 813 // Double_t prob = r.Prob();
c389303e 814 Bool_t ret = kTRUE;
815
816 // Check that the reduced chi square isn't larger than 5
817 if (ndf <= 0 || chi2 / ndf > chi2nuCut) {
0bd4b00f 818 if (fDebug > 2) {
c389303e 819 AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f",
820 n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
0bd4b00f 821 chi2nuCut)); }
c389303e 822 ret = kFALSE;
f8715167 823 }
824
c389303e 825 // Check each parameter
4b9857f3 826 UShort_t nPar = r->NPar();
f8715167 827 for (UShort_t i = 0; i < nPar; i++) {
c389303e 828 if (i == kN) continue; // Skip the number parameter
f8715167 829
c389303e 830 // Get value and error and check value
831 Double_t v = r->Parameter(i);
832 Double_t e = r->ParError(i);
f8715167 833 if (v == 0) continue;
c389303e 834
835 // Calculate the relative error and check it
836 Double_t re = e / v;
837 Double_t cut = relErrorCut * (i < kN ? 1 : 2);
838 if (re > cut) {
0bd4b00f 839 if (fDebug > 2) {
c389303e 840 AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
841 n.Data(), r->ParName(i).c_str(),
842 r->ParName(i).c_str(), e, v,
843 100*(v == 0 ? 0 : e/v),
0bd4b00f 844 100*(cut))); }
c389303e 845 ret = kFALSE;
f8715167 846 }
847 }
c389303e 848
849 // Check if we have scale parameters
850 if (nPar > kN) {
851
852 // Check that the last particle has a significant contribution
4b9857f3 853 Double_t lastScale = r->Parameter(nPar-1);
854 if (lastScale <= 1e-7) {
0bd4b00f 855 if (fDebug) {
c389303e 856 AliWarning(Form("%s: %s=%9.6f<1e-7",
0bd4b00f 857 n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
c389303e 858 ret = kFALSE;
f8715167 859 }
860 }
c389303e 861 return ret;
f8715167 862}
863
864
0bd4b00f 865//__________________________________________________________________
866void
867AliFMDEnergyFitter::RingHistos::FindBestFits(TList* d,
868 AliFMDCorrELossFit& obj,
869 const TAxis& eta,
870 Double_t relErrorCut,
871 Double_t chi2nuCut,
872 Double_t minWeightCut)
873{
874 // Get our ring list from the output container
875 TList* l = GetOutputList(d);
876 if (!l) return;
877
878 // Get the energy distributions from the output container
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;
885 }
886 Int_t nBin = eta.GetNbins();
887
888 for (Int_t b = 1; b <= nBin; b++) {
889 TString n(Form(fgkEDistFormat, fName.Data(), b));
890 TH1D* dist = static_cast<TH1D*>(dists->FindObject(n));
891 // Ignore empty histograms altoghether
892 if (!dist || dist->GetEntries() <= 0) continue;
893
894 AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
895 relErrorCut,
896 chi2nuCut,
897 minWeightCut);
898 best->fDet = fDet;
899 best->fRing = fRing;
900 best->fBin = b; //
901 // Double_t eta = fAxis->GetBinCenter(b);
902 obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
903 }
904}
905
906//__________________________________________________________________
907AliFMDCorrELossFit::ELossFit*
908AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist,
909 Double_t relErrorCut,
910 Double_t chi2nuCut,
911 Double_t minWeightCut)
912{
913 TList* funcs = dist->GetListOfFunctions();
914 TIter next(funcs);
915 TF1* func = 0;
916 fFits.Clear();
917 Int_t i = 0;
918 // Info("FindBestFit", "%s", dist->GetName());
919 while ((func = static_cast<TF1*>(next()))) {
920 AliFMDCorrELossFit::ELossFit* fit =
921 new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
922 fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
923 }
924 fFits.Sort(false);
925 // fFits.Print();
926 return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
927}
928
929
f8715167 930//____________________________________________________________________
931void
932AliFMDEnergyFitter::RingHistos::Output(TList* dir)
933{
934 fList = DefineOutputList(dir);
935}
936
937//____________________________________________________________________
938//
939// EOF
940//