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