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