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