2 // Class to do the energy correction of FMD ESD data
4 #include "AliFMDEnergyFitter.h"
5 #include "AliForwardUtil.h"
14 #include <TClonesArray.h>
15 #include <TFitResult.h>
21 ClassImp(AliFMDEnergyFitter)
23 ; // This is for Emacs
26 const char* fgkEDistFormat = "%s_etabin%03d";
30 //____________________________________________________________________
31 AliFMDEnergyFitter::AliFMDEnergyFitter()
44 fUseIncreasingBins(true),
45 fMaxRelParError(AliFMDCorrELossFit::ELossFit::fgMaxRelError),
46 fMaxChi2PerNDF(AliFMDCorrELossFit::ELossFit::fgMaxChi2nu),
47 fMinWeight(AliFMDCorrELossFit::ELossFit::fgLeastWeight),
49 fResidualMethod(kNoResiduals),
51 fRegularizationCut(3e6)
54 // Default Constructor - do not use
56 // DGUARD(fDebug, 1, "Default CTOR of AliFMDEnergyFitter");
57 fRingHistos.SetOwner();
60 //____________________________________________________________________
61 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
62 : TNamed("fmdEnergyFitter", title),
71 fCentralityAxis(0,0,0),
74 fUseIncreasingBins(true),
75 fMaxRelParError(AliFMDCorrELossFit::ELossFit::fgMaxRelError),
76 fMaxChi2PerNDF(AliFMDCorrELossFit::ELossFit::fgMaxChi2nu),
77 fMinWeight(AliFMDCorrELossFit::ELossFit::fgLeastWeight),
79 fResidualMethod(kNoResiduals),
81 fRegularizationCut(3e6)
87 // title Title of object - not significant
89 // DGUARD(fDebug, 1, "Named CTOR of AliFMDEnergyFitter: %s", title);
90 fEtaAxis.SetName("etaAxis");
91 fEtaAxis.SetTitle("#eta");
92 fCentralityAxis.SetName("centralityAxis");
93 fCentralityAxis.SetTitle("Centrality [%]");
96 //____________________________________________________________________
97 AliFMDEnergyFitter::~AliFMDEnergyFitter()
102 DGUARD(fDebug, 1, "DTOR of AliFMDEnergyFitter");
103 // fRingHistos.Delete();
106 //____________________________________________________________________
107 AliFMDEnergyFitter::RingHistos*
108 AliFMDEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const
110 return new RingHistos(d,r);
113 //____________________________________________________________________
114 AliFMDEnergyFitter::RingHistos*
115 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
118 // Get the ring histogram container
125 // Ring histogram container
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;
133 if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
135 return static_cast<RingHistos*>(fRingHistos.At(idx));
138 //____________________________________________________________________
140 AliFMDEnergyFitter::Init()
142 // Create the ring histograms.
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);
153 while ((o = static_cast<RingHistos*>(next()))) {
158 //____________________________________________________________________
160 AliFMDEnergyFitter::CreateOutputObjects(TList* dir)
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.
168 // dir Directory to add to
170 DGUARD(fDebug, 1, "Define output in AliFMDEnergyFitter");
171 TList* d = new TList;
172 d->SetName(GetName());
176 // Store eta axis as a histogram, since that can be merged!
178 if (fEtaAxis.GetXbins()->GetArray())
179 hEta = new TH1I(fEtaAxis.GetName(), fEtaAxis.GetTitle(),
180 fEtaAxis.GetNbins(), fEtaAxis.GetXbins()->GetArray());
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);
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));
203 if (fRingHistos.GetEntries() <= 0) {
204 AliFatal("No ring histograms where defined - giving up!");
207 TIter next(&fRingHistos);
209 while ((o = static_cast<RingHistos*>(next()))) {
211 o->CreateOutputObjects(d);
215 //____________________________________________________________________
217 AliFMDEnergyFitter::SetupForData(const TAxis& eAxis)
220 // Initialise the task - called at first event
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
227 DGUARD(fDebug, 1, "Initialize of AliFMDEnergyFitter");
228 if (fEtaAxis.GetNbins() == 0 ||
229 TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7)
231 if (fCentralityAxis.GetNbins() == 0) {
233 Double_t bins[] = { 0., 5., 10., 15., 20., 30.,
234 40., 50., 60., 70., 80., 100. };
235 SetCentralityAxis(n, bins);
237 TIter next(&fRingHistos);
239 while ((o = static_cast<RingHistos*>(next())))
240 o->SetupForData(fEtaAxis, fCentralityAxis, fMaxE,
241 fNEbins, fUseIncreasingBins);
243 //____________________________________________________________________
245 AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
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
255 // etaAxis Eta axis to use
257 SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
259 //____________________________________________________________________
261 AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
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
271 // nBins Number of bins
272 // etaMin Minimum of the eta axis
273 // etaMax Maximum of the eta axis
275 fEtaAxis.Set(nBins,etaMin,etaMax);
277 //____________________________________________________________________
279 AliFMDEnergyFitter::SetCentralityAxis(UShort_t n, Double_t* bins)
281 fCentralityAxis.Set(n-1, bins);
285 //____________________________________________________________________
287 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
292 // Fitter the input AliESDFMD object
297 // empty Whether the event is 'empty'
300 // True on success, false otherwise
301 DGUARD(fDebug, 5, "Accumulate statistics in AliFMDEnergyFitter - cholm");
302 Int_t icent = fCentralityAxis.FindBin(cent);
303 if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0;
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);
313 RingHistos* histos = GetRingHistos(d, r);
314 if (!histos) continue;
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);
320 // Keep dead-channel information.
321 if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0)
324 // Get the pseudo-rapidity
325 Double_t eta1 = input.Eta(d,r,s,t);
326 // Int_t ieta = fEtaAxis.FindBin(eta1);
327 // if (ieta < 1 || ieta > fEtaAxis.GetNbins()) continue;
329 // histos->Fill(empty, ieta-1, icent, mult);
330 histos->Fill(empty, eta1, icent, mult);
337 DMSG(fDebug, 3, "Found a total of %d signals for c=%f, and %sempty event",
338 nFills, cent, (empty ? "" : "non-"));
342 //____________________________________________________________________
344 AliFMDEnergyFitter::Fit(const TList* dir)
347 // Scale the histograms to the total number of events
350 // dir Where the histograms are
352 DGUARD(fDebug, 1, "Fit distributions in AliFMDEnergyFitter");
353 if (!fDoFits) return;
355 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
357 AliWarningF("No list named %s found in %s", GetName(), dir->GetName());
364 // +fNParticles-1 for weights
365 Int_t nStack = kN+fNParticles;
366 THStack* stack[nStack];
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));
376 for (Int_t i = 0; i < nStack; i++)
379 TIter next(&fRingHistos);
381 while ((o = static_cast<RingHistos*>(next()))) {
382 if (CheckSkip(o->fDet, o->fRing, fSkips)) {
383 AliWarningF("Skipping FMD%d%c for fitting", o->fDet, o->fRing);
387 TObjArray* l = o->Fit(d, fLowCut, fNParticles,
388 fMinEntries, fFitRangeBinWidth,
389 fMaxRelParError, fMaxChi2PerNDF,
390 fMinWeight, fRegularizationCut,
393 for (Int_t i = 0; i < l->GetEntriesFast()-1; i++) { // Last is status
394 stack[i % nStack]->Add(static_cast<TH1*>(l->At(i)));
398 if (!fDoMakeObject) return;
400 MakeCorrectionsObject(d);
403 //____________________________________________________________________
405 AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
408 // Generate the corrections object
411 // dir List to analyse
413 DGUARD(fDebug, 1, "Make the correction objec in AliFMDEnergyFitter");
415 AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
416 obj->SetEtaAxis(fEtaAxis);
417 obj->SetLowCut(fLowCut);
419 TIter next(&fRingHistos);
421 while ((o = static_cast<RingHistos*>(next()))) {
422 if (CheckSkip(o->fDet, o->fRing, fSkips)) {
423 AliWarningF("Skipping FMD%d%c for correction object", o->fDet, o->fRing);
427 o->FindBestFits(d, *obj, fEtaAxis);
429 d->Add(obj, "elossFits");
432 //____________________________________________________________________
434 AliFMDEnergyFitter::SetDebug(Int_t dbg)
437 // Set the debug level. The higher the value the more output
443 TIter next(&fRingHistos);
445 while ((o = static_cast<RingHistos*>(next())))
449 //____________________________________________________________________
451 template <typename T>
452 void GetParam(Bool_t& ret, const TCollection* col,
453 const TString& name, T& var)
455 TObject* o = col->FindObject(name);
456 if (o) AliForwardUtil::GetParameter(o,var);
461 //____________________________________________________________________
463 AliFMDEnergyFitter::ReadParameters(const TCollection* col)
465 // Read parameters of this object from a collection
468 // col Collection to read parameters from
471 // true on success, false otherwise
473 if (!col) return false;
475 TAxis* axis = static_cast<TAxis*>(col->FindObject("etaAxis"));
476 if (!axis) ret = false;
478 if (axis->GetXbins()->GetArray())
479 fEtaAxis.Set(axis->GetNbins(), axis->GetXbins()->GetArray());
481 fEtaAxis.Set(axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
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);
496 GetParam(dummy,col,"regCut", fRegularizationCut);
501 //____________________________________________________________________
503 AliFMDEnergyFitter::CheckSkip(UShort_t d, Char_t r, UShort_t skips)
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);
509 return (t & skips) == t;
512 //____________________________________________________________________
514 AliFMDEnergyFitter::Print(Option_t*) const
522 char ind[gROOT->GetDirLevel()+1];
523 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
524 ind[gROOT->GetDirLevel()] = '\0';
525 std::cout << ind << ClassName() << ": " << GetName() << '\n'
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'
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;
549 std::cout << std::endl;
552 //====================================================================
553 AliFMDEnergyFitter::RingHistos::RingHistos()
554 : AliForwardUtil::RingHistos(),
561 fFits("AliFMDCorrELossFit::ELossFit", 200),
567 DGUARD(fDebug, 3, "Default CTOR AliFMDEnergyFitter::RingHistos");
571 //____________________________________________________________________
572 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
573 : AliForwardUtil::RingHistos(d,r),
580 fFits("AliFMDCorrELossFit::ELossFit", 200),
590 DGUARD(fDebug, 3, "Named CTOR AliFMDEnergyFitter::RingHistos: FMD%d%c",
594 //____________________________________________________________________
595 AliFMDEnergyFitter::RingHistos::~RingHistos()
600 DGUARD(fDebug, 3, "DTOR of AliFMDEnergyFitter::RingHistos");
601 // if (fEDist) delete fEDist;
602 // fEtaEDists.Delete();
604 //__________________________________________________________________
606 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
610 // Make an axis with increasing bins
618 // An axis with quadratically increasing bin size
622 Double_t dx = (max-min) / n;
625 for (i = 1; i < n+1; i++) {
626 Double_t dI = float(i)/n;
627 Double_t next = bins[i-1] + dx + dI * dI * dx;
629 if (next > max) break;
635 //____________________________________________________________________
637 AliFMDEnergyFitter::RingHistos::Make(const char* name,
645 // Make E/E_mip histogram
648 // deMax Maximum energy loss
649 // nDeBins Number energy loss bins
650 // incr Whether to make bins of increasing size
655 TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
656 mAxis.Set(deAxis.GetSize()-1, deAxis.GetArray());
659 mAxis.Set(nDeBins, 0, deMax);
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());
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());
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());
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());
692 h->SetYTitle("#sum#DeltaE/#DeltaE_{mip}");
693 h->SetXTitle("#eta");
694 h->SetFillColor(Color());
695 h->SetMarkerColor(Color());
696 h->SetLineColor(Color());
697 h->SetFillStyle(3001);
698 h->SetMarkerStyle(20);
703 //____________________________________________________________________
705 AliFMDEnergyFitter::RingHistos::CreateOutputObjects(TList* dir)
713 DGUARD(fDebug, 2, "Define output in AliFMDEnergyFitter::RingHistos");
714 fList = DefineOutputList(dir);
716 // fEtaEDists = new TList;
717 // fEtaEDists->SetOwner();
718 // fEtaEDists->SetName("EDists");
720 // fList->Add(fEtaEDists);
722 //____________________________________________________________________
724 AliFMDEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis,
725 const TAxis& /* cAxis */,
735 // maxDE Max energy loss to consider
736 // nDEbins Number of bins
737 // useIncrBin Whether to use an increasing bin size
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()),
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);
750 fEDist->SetDirectory(0);
752 fEmpty = static_cast<TH1D*>(fEDist->Clone(Form("%s_empty", fName.Data())));
753 fEmpty->SetTitle(Form("#sum#DeltaE/#DeltaE_{mip} for %s (empty events)",
755 fEmpty->SetDirectory(0);
759 fHist = Make("eloss", "#sum#Delta/#Delta_{mip}", eAxis,
760 maxDE, nDEbins, useIncrBin);
766 //____________________________________________________________________
768 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty,
777 // empty True if event is empty
778 // ieta Eta bin (0-based)
779 // icent Centrality bin (1-based)
782 DGUARD(fDebug, 10, "Filling in AliFMDEnergyFitter::RingHistos");
789 // if (!fEtaEDists) {
791 Warning("Fill", "No list of E dists defined");
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
798 // if (ieta < 0 || ieta >= fEtaEDists->GetEntries()) return;
800 // TH1D* h = static_cast<TH1D*>(fEtaEDists->At(ieta));
806 //____________________________________________________________________
808 AliFMDEnergyFitter::RingHistos::MakePar(const char* name,
810 const TAxis& eta) const
813 // Make a parameter histogram
816 // name Name of histogram.
817 // title Title of histogram.
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");
829 h->SetFillColor(Color());
830 h->SetMarkerColor(Color());
831 h->SetLineColor(Color());
832 h->SetFillStyle(3001);
837 //____________________________________________________________________
839 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name,
848 // Make a histogram that contains the results of the fit over the full ring
856 // val Value of parameter
857 // err Error on parameter
860 // The newly allocated histogram
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()),
867 h->SetBinContent(1, val);
868 h->SetBinError(1, err);
869 h->SetXTitle("#eta");
873 h->SetMarkerColor(Color());
874 h->SetLineColor(Color());
882 //____________________________________________________________________
884 AliFMDEnergyFitter::RingHistos::Fit(TList* dir,
889 Double_t relErrorCut,
893 EResidualMethod residuals) const
895 return FitSlices(dir, "eloss", lowCut, nParticles, minEntries,
896 minusBins, relErrorCut, chi2nuCut, minWeight, regCut,
897 residuals, true, &fBest);
900 //____________________________________________________________________
902 AliFMDEnergyFitter::RingHistos::FitSlices(TList* dir,
908 Double_t relErrorCut,
912 EResidualMethod residuals,
914 TObjArray* best) const
917 // Fit each histogram to up to @a nParticles particle responses.
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
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$
933 DGUARD(fDebug, 2, "Fit in AliFMDEnergyFitter::RingHistos");
935 // Get our ring list from the output container
936 TList* l = GetOutputList(dir);
939 // Get the energy distributions from the output container
940 // TList* dists = static_cast<TList*>(l->FindObject("EDists"));
942 // AliWarning(Form("Didn't find EtaEDists (%s) in %s",
943 // fName.Data(), l->GetName()));
947 TH2* h = static_cast<TH2*>(l->FindObject(name));
949 AliWarningF("Didn't find 2D histogram '%s' in %s", name, l->GetName());
953 const TAxis& eta = *(h->GetXaxis());
955 // Create an output list for the fitted distributions
956 TList* out = new TList;
958 out->SetName(Form("%sDists", name));
961 // Optional container for residuals
963 if (residuals != kNoResiduals) {
965 resi->SetName(Form("%sResiduals", name));
970 // Container of the fit results histograms
971 TObjArray* pars = new TObjArray(3+nParticles+1);
972 pars->SetName(Form("%sResults", name));
975 // Result objects stored in separate list on the output
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++)
992 pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
995 Int_t nDists = h->GetNbinsX(); // dists->GetEntries();
1002 best->Expand(nDists+1);
1004 best->SetOwner(false);
1006 for (Int_t i = 0; i < nDists; i++) {
1007 // TH1D* dist = static_cast<TH1D*>(dists->At(i));
1008 // Ignore empty histograms altoghether
1010 TH1D* dist = h->ProjectionY(Form(fgkEDistFormat,GetName(),b),b,b,"e");
1012 // If we got the null pointer, return 0
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)));
1024 UShort_t status1 = 0;
1025 ELossFit_t* res = FitHist(dist,
1038 case 1: nEmpty++; break;
1039 case 2: nLow++; break;
1045 // Now count as fitted, store as best fits, and add distribution
1046 // to the dedicated output list
1048 res->fBin = b; // We only have the bin information here
1049 if (best) best->AddAt(res, b);
1051 // dist->GetListOfFunctions()->Add(res);
1053 // If asked to calculate residuals, do so, and store result on the
1054 // dedicated output list
1055 if (residuals != kNoResiduals && resi)
1056 CalculateResiduals(residuals, lowCut, dist, res, resi);
1059 low = TMath::Min(low,b);
1060 high = TMath::Max(high,b);
1062 // Get the reduced chi square
1063 Double_t chi2 = res->fChi2; // GetChisquare();
1064 Int_t ndf = res->fNu; // GetNDF();
1066 // Store results of best fit in output histograms
1067 // res->SetLineWidth(4);
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);
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));
1083 for (UShort_t j = 0; j < nParticles-1; j++) {
1084 hA[j]->SetBinContent(b, res->fA[j]);
1085 hA[j]->SetBinError(b, res->fEA[j]);
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);
1092 // Fit the full-ring histogram
1093 TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
1095 UShort_t statusT = 0;
1096 ELossFit_t* resT = FitHist(total,
1108 // Make histograms for the result of this fit
1109 Double_t chi2 = resT->GetChi2();
1110 Int_t ndf = resT->GetNu();
1111 pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
1112 ndf > 0 ? chi2/ndf : 0, 0));
1113 pars->Add(MakeTotal("t_c", "Constant", eta, low, high,
1114 resT->GetC(), resT->GetEC()));
1115 pars->Add(MakeTotal("t_delta", "#Delta_{p}", eta, low, high,
1116 resT->GetDelta(), resT->GetEDelta()));
1117 pars->Add(MakeTotal("t_xi", "#xi", eta, low, high,
1118 resT->GetXi(), resT->GetEXi()));
1119 pars->Add(MakeTotal("t_sigma", "#sigma", eta, low, high,
1120 resT->GetSigma(), resT->GetESigma()));
1121 pars->Add(MakeTotal("t_sigman", "#sigma_{n}", eta, low, high,
1122 resT->GetSigmaN(), resT->GetESigmaN()));
1123 pars->Add(MakeTotal("t_n", "N", eta, low, high,
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,
1128 resT->GetA(j), resT->GetEA(j)));
1132 TH1* status = new TH1I(Form("%sStatus",name), "Status of Fits", 5, 0, 5);
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);
1155 //____________________________________________________________________
1156 AliFMDEnergyFitter::RingHistos::ELossFit_t*
1157 AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,
1159 UShort_t nParticles,
1160 UShort_t minEntries,
1162 Double_t relErrorCut,
1167 UShort_t& status) const
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$
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$
1192 // The best fit function
1194 DGUARD(fDebug, 2, "Fit histogram in AliFMDEnergyFitter::RingHistos: %s",
1196 Double_t maxRange = 10;
1199 if (dist->GetEntries() <= 0) {
1200 status = 1; // `empty'
1204 // Scale to the bin-width
1205 dist->Scale(1., "width");
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),
1211 dist->GetXaxis()->SetRange(cutBin, maxBin);
1213 // Get the bin with maximum
1214 Int_t peakBin = dist->GetMaximumBin();
1216 // Normalize peak to 1
1217 // Double_t max = dist->GetMaximum();
1218 Double_t max = dist->GetBinContent(peakBin); // Maximum();
1220 status = 1; // `empty'
1223 if (scaleToPeak) dist->Scale(1/max);
1224 DMSG(fDebug,5,"max(%s) -> %f", dist->GetName(), max);
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));
1235 // Create a fitter object
1236 AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins);
1238 f.SetDebug(fDebug > 2);
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);
1250 // If we are only asked to fit a single particle, return this fit,
1252 if (nParticles == 1) {
1253 TF1* r = f.Fit1Particle(dist, 0);
1255 status = 3; // No-fit
1258 TF1* ff = new TF1(*r);
1259 dist->GetListOfFunctions()->Add(ff);
1260 ELossFit_t* ret = new ELossFit_t(0, *ff);
1261 ret->CalculateQuality(chi2nuCut, relErrorCut, minWeight);
1266 // Fit from 2 upto n particles
1267 for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
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));
1278 // Here, we use the real quality assesor instead of the old
1279 // `CheckResult' to ensure consitency in all output.
1280 ELossFit_t* ret = FindBestFit(dist, relErrorCut, chi2nuCut, minWeight);
1281 if (!ret) status = 3;
1285 //__________________________________________________________________
1286 AliFMDEnergyFitter::RingHistos::ELossFit_t*
1287 AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist,
1288 Double_t relErrorCut,
1290 Double_t minWeightCut) const
1293 // Find the best fit
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$
1307 TList* funcs = dist->GetListOfFunctions();
1311 fFits.Clear(); // This is only ever used here
1313 if (fDebug) printf("Find best fit for %s ... ", dist->GetName());
1314 if (fDebug > 2) printf("\n");
1316 // Loop over all functions stored in distribution,
1317 // and calculate the quality
1318 while ((func = static_cast<TF1*>(next()))) {
1319 ELossFit_t* fit = new(fFits[i++]) ELossFit_t(0,*func);
1323 fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
1325 Printf("%10s: %3d (chi^2/nu: %6.3f)",
1326 func->GetName(), fit->fQuality,
1327 (fit->fNu > 0 ? fit->fChi2 / fit->fNu : 999));
1330 // Sort all the found fit objects in increasing quality
1332 if (fDebug > 1) fFits.Print("s");
1334 // Get the top-most fit
1335 ELossFit_t* ret = static_cast<ELossFit_t*>(fFits.At(i-1));
1337 AliWarningF("No fit found for %s", GetName());
1340 if (ret && fDebug > 0) {
1341 if (fDebug > 1) printf(" %d: ", i-1);
1344 // We have to make a copy here, because other wise the clones array
1345 // will overwrite the address
1346 return new ELossFit_t(*ret);
1349 //____________________________________________________________________
1351 AliFMDEnergyFitter::RingHistos::CalculateResiduals(EResidualMethod mode,
1355 TCollection* out) const
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);
1365 // Set title on Y axis
1367 case kResidualDifference:
1368 resi->SetYTitle("h_{i}-f(#Delta_{i}) #pm #delta_{i}");
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}");
1376 resi->SetYTitle("Unknown");
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);
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;
1399 Double_t h = dist->GetBinContent(i);
1400 Double_t e = dist->GetBinError(i);
1403 if (h > 0 && e > 0) {
1404 Double_t f = fit->GetC() * fit->Evaluate(x);
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;
1415 resi->SetBinContent(i, r);
1416 resi->SetBinError(i, er);
1420 //__________________________________________________________________
1422 AliFMDEnergyFitter::RingHistos::FindBestFits(const TList* d,
1423 AliFMDCorrELossFit& obj,
1427 // Find the best fits
1431 // obj Object to add fits to
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$
1441 // Get our ring list from the output container
1442 TList* l = GetOutputList(d);
1445 // Get the energy distributions from the output container
1446 TList* dists = static_cast<TList*>(l->FindObject("elossDists"));
1448 AliWarningF("Didn't find elossDists in %s", l->GetName());
1452 Int_t nBin = eta.GetNbins();
1453 if (fBest.GetEntriesFast() <= 0) {
1454 AliWarningF("No fits found for %s", GetName());
1458 for (Int_t b = 1; b <= nBin; b++) {
1459 ELossFit_t* best = static_cast<ELossFit_t*>(fBest.At(b));
1461 // AliErrorF("No best fit found @ %d for %s", b, GetName());
1464 // FindBestFit(b, dist, relErrorCut, chi2nuCut, minWeightCut);
1466 best->fRing = fRing;
1469 printf("Bin # %3d: ", b);
1472 // Double_t eta = fAxis->GetBinCenter(b);
1473 obj.SetFit(fDet, fRing, b, best);
1474 // new ELossFit_t(*best));
1480 //____________________________________________________________________