X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PWG2%2FFORWARD%2Fanalysis2%2FAliFMDEnergyFitter.cxx;h=76575473ca8fc07970b510add1284273e584c505;hp=34464d80aecc8863a4675862edfb2e05909be8fb;hb=0bd4b00f0c31630362b92dc62b75eb8252459127;hpb=66b34429bc513dde71251ad76062880776bbdff7 diff --git a/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx b/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx index 34464d80aec..76575473ca8 100644 --- a/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx @@ -2,125 +2,47 @@ // Class to do the energy correction of FMD ESD data // #include "AliFMDEnergyFitter.h" +#include "AliForwardCorrectionManager.h" #include #include #include #include #include #include -#include "AliFMDAnaParameters.h" #include #include #include #include +#include +#include +#include ClassImp(AliFMDEnergyFitter) #if 0 ; // This is for Emacs #endif - -#define N_A(N) (4+N-2) -#define N2_A(N) (4+(N-2)*3) -#define N2_D(N) (4+(N-2)*3+1) -#define N2_X(N) (4+(N-2)*3+2) - -//____________________________________________________________________ namespace { - Double_t - NLandau(Double_t* xp, Double_t* pp) - { - Double_t e = xp[0]; - Double_t constant = pp[0]; - Double_t mpv = pp[1]; - Double_t fwhm = pp[2]; - Int_t n = Int_t(pp[3]); - Double_t result = 0; - for (Int_t i = 1; i <= n; i++) { - Double_t mpvi = i*(mpv+fwhm*TMath::Log(i)); - Double_t fwhmi = i*fwhm; - Double_t ai = (i == 1 ? 1 : pp[N_A(i)]); - result += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE); - } - result *= constant; - return result; - } - - Double_t - NLandau2(Double_t* xp, Double_t* pp) - { - Double_t e = xp[0]; - Double_t constant = pp[0]; - Double_t mpv = pp[1]; - Double_t fwhm = pp[2]; - Int_t n = Int_t(pp[3]); - Double_t result = TMath::Landau(e,mpv,fwhm,kTRUE); - for (Int_t i = 2; i <= n; i++) { - Double_t ai = pp[N2_A(i)]; - Double_t mpvi = pp[N2_D(i)]; - Double_t fwhmi = pp[N2_X(i)]; - result += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE); - } - result *= constant; - return result; - } - - Double_t - TripleLandau(Double_t *x, Double_t *par) - { - Double_t energy = x[0]; - Double_t constant = par[0]; - Double_t mpv = par[1]; - Double_t fwhm = par[2]; - Double_t alpha = par[3]; - Double_t beta = par[4]; - Double_t mpv2 = 2*mpv+2*fwhm*TMath::Log(2); - Double_t mpv3 = 3*mpv+3*fwhm*TMath::Log(3); - - Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+ - alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE)+ - beta * TMath::Landau(energy,mpv3,3*fwhm,kTRUE)); - - return f; - } - Double_t - DoubleLandau(Double_t *x, Double_t *par) - { - Double_t energy = x[0]; - Double_t constant = par[0]; - Double_t mpv = par[1]; - Double_t fwhm = par[2]; - Double_t alpha = par[3]; - Double_t mpv2 = 2*mpv+2*fwhm*TMath::Log(2); - - Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+ - alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE)); - - return f; - } - Double_t - SingleLandau(Double_t *x, Double_t *par) - { - Double_t energy = x[0]; - Double_t constant = par[0]; - Double_t mpv = par[1]; - Double_t fwhm = par[2]; - - Double_t f = constant * TMath::Landau(energy,mpv,fwhm,kTRUE); - - return f; - } + const char* fgkEDistFormat = "%s_etabin%03d"; } + //____________________________________________________________________ AliFMDEnergyFitter::AliFMDEnergyFitter() : TNamed(), fRingHistos(), fLowCut(0.3), - fNLandau(3), + fNParticles(3), fMinEntries(100), - fBinsToSubtract(4), + fFitRangeBinWidth(4), fDoFits(false), + fDoMakeObject(false), fEtaAxis(), + fMaxE(10), + fNEbins(300), + fUseIncreasingBins(true), + fMaxRelParError(.25), + fMaxChi2PerNDF(20), + fMinWeight(1e-7), fDebug(0) {} @@ -129,11 +51,18 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title) : TNamed("fmdEnergyFitter", title), fRingHistos(), fLowCut(0.3), - fNLandau(3), + fNParticles(3), fMinEntries(100), - fBinsToSubtract(4), + fFitRangeBinWidth(4), fDoFits(false), - fEtaAxis(100,-4,6), + fDoMakeObject(false), + fEtaAxis(0,0,0), + fMaxE(10), + fNEbins(300), + fUseIncreasingBins(true), + fMaxRelParError(.25), + fMaxChi2PerNDF(20), + fMinWeight(1e-7), fDebug(3) { fEtaAxis.SetName("etaAxis"); @@ -150,11 +79,18 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o) : TNamed(o), fRingHistos(), fLowCut(o.fLowCut), - fNLandau(o.fNLandau), + fNParticles(o.fNParticles), fMinEntries(o.fMinEntries), - fBinsToSubtract(o.fBinsToSubtract), + fFitRangeBinWidth(o.fFitRangeBinWidth), fDoFits(o.fDoFits), + fDoMakeObject(o.fDoMakeObject), fEtaAxis(o.fEtaAxis), + fMaxE(o.fMaxE), + fNEbins(o.fNEbins), + fUseIncreasingBins(o.fUseIncreasingBins), + fMaxRelParError(o.fMaxRelParError), + fMaxChi2PerNDF(o.fMaxChi2PerNDF), + fMinWeight(o.fMinWeight), fDebug(o.fDebug) { TIter next(&o.fRingHistos); @@ -175,12 +111,16 @@ AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o) TNamed::operator=(o); fLowCut = o.fLowCut; - fNLandau = o.fNLandau; + fNParticles = o.fNParticles; fMinEntries = o.fMinEntries; - fBinsToSubtract= o.fBinsToSubtract; + fFitRangeBinWidth= o.fFitRangeBinWidth; fDoFits = o.fDoFits; + fDoMakeObject = o.fDoMakeObject; fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()); fDebug = o.fDebug; + fMaxRelParError= o.fMaxRelParError; + fMaxChi2PerNDF = o.fMaxChi2PerNDF; + fMinWeight = o.fMinWeight; fRingHistos.Delete(); TIter next(&o.fRingHistos); @@ -209,12 +149,27 @@ AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const void AliFMDEnergyFitter::Init(const TAxis& eAxis) { - fEtaAxis.Set(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax()); + if (fEtaAxis.GetNbins() == 0 || + fEtaAxis.GetXmin() == fEtaAxis.GetXmax()) + SetEtaAxis(eAxis); TIter next(&fRingHistos); RingHistos* o = 0; while ((o = static_cast(next()))) - o->Init(eAxis); + o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins); } +//____________________________________________________________________ +void +AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis) +{ + SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax()); +} +//____________________________________________________________________ +void +AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax) +{ + fEtaAxis.Set(nBins,etaMin,etaMax); +} + //____________________________________________________________________ Bool_t AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, @@ -264,30 +219,60 @@ AliFMDEnergyFitter::Fit(TList* dir) // +1 for chi^2 // +3 for 1 landau // +1 for N - // +fNLandau-1 for weights - Int_t nStack = 1+3+1+fNLandau-1; + // +fNParticles-1 for weights + Int_t nStack = kN+fNParticles; THStack* stack[nStack]; - stack[0] = new THStack("chi2", "#chi^{2}/#nu"); - stack[1] = new THStack("c", "constant"); - stack[2] = new THStack("mpv", "#Delta_{p}"); - stack[3] = new THStack("w", "w"); - stack[4] = new THStack("n", "# of Landaus"); - for (Int_t i = 2; i <= fNLandau; i++) - stack[i-1+4] = new THStack(Form("a%d", i), - Form("Weight of %d signal", i)); + stack[0] = new THStack("chi2", "#chi^{2}/#nu"); + stack[kC +1] = new THStack("c", "Constant"); + stack[kDelta +1] = new THStack("delta", "#Delta_{p}"); + stack[kXi +1] = new THStack("xi", "#xi"); + stack[kSigma +1] = new THStack("sigma", "#sigma"); + stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}"); + stack[kN +1] = new THStack("n", "# Particles"); + for (Int_t i = 2; i <= fNParticles; i++) + stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i)); for (Int_t i = 0; i < nStack; i++) d->Add(stack[i]); TIter next(&fRingHistos); RingHistos* o = 0; while ((o = static_cast(next()))) { - TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNLandau, - fMinEntries, fBinsToSubtract); + TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles, + fMinEntries, fFitRangeBinWidth, + fMaxRelParError, fMaxChi2PerNDF); if (!l) continue; for (Int_t i = 0; i < l->GetEntriesFast(); i++) { stack[i % nStack]->Add(static_cast(l->At(i))); } } + + if (!fDoMakeObject) return; + + MakeCorrectionsObject(d); +} + +//____________________________________________________________________ +void +AliFMDEnergyFitter::MakeCorrectionsObject(TList* d) +{ + AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance(); + + AliFMDCorrELossFit* obj = new AliFMDCorrELossFit; + obj->SetEtaAxis(fEtaAxis); + obj->SetLowCut(fLowCut); + + TIter next(&fRingHistos); + RingHistos* o = 0; + while ((o = static_cast(next()))) { + o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError, + fMaxChi2PerNDF, fMinWeight); + } + + TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits)); + TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits)); + AliInfo(Form("Object %s created in output - should be extracted and copied " + "to %s", oName.Data(), fileName.Data())); + d->Add(obj, oName.Data()); } //____________________________________________________________________ @@ -314,6 +299,32 @@ AliFMDEnergyFitter::SetDebug(Int_t dbg) while ((o = static_cast(next()))) o->fDebug = dbg; } + +//____________________________________________________________________ +void +AliFMDEnergyFitter::Print(Option_t*) const +{ + char ind[gROOT->GetDirLevel()+1]; + for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; + ind[gROOT->GetDirLevel()] = '\0'; + std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n' + << ind << " Low cut: " << fLowCut << " E/E_mip\n" + << ind << " Max(particles): " << fNParticles << '\n' + << ind << " Min(entries): " << fMinEntries << '\n' + << ind << " Fit range: " + << fFitRangeBinWidth << " bins\n" + << ind << " Make fits: " + << (fDoFits ? "yes\n" : "no\n") + << ind << " Make object: " + << (fDoMakeObject ? "yes\n" : "no\n") + << ind << " Max E/E_mip: " << fMaxE << '\n' + << ind << " N bins: " << fNEbins << '\n' + << ind << " Increasing bins: " + << (fUseIncreasingBins ?"yes\n" : "no\n") + << ind << " max(delta p/p): " << fMaxRelParError << '\n' + << ind << " max(chi^2/nu): " << fMaxChi2PerNDF << '\n' + << ind << " min(a_i): " << fMinWeight << std::endl; +} //==================================================================== AliFMDEnergyFitter::RingHistos::RingHistos() @@ -322,6 +333,7 @@ AliFMDEnergyFitter::RingHistos::RingHistos() fEmpty(0), fEtaEDists(), fList(0), + fFits("AliFMDCorrELossFit::ELossFit"), fDebug(0) {} @@ -332,6 +344,7 @@ AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r) fEmpty(0), fEtaEDists(), fList(0), + fFits("AliFMDCorrELossFit::ELossFit"), fDebug(0) { fEtaEDists.SetName("EDists"); @@ -343,8 +356,10 @@ AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o) fEmpty(o.fEmpty), fEtaEDists(), fList(0), + fFits("AliFMDCorrELossFit::ELossFit"), fDebug(0) { + fFits.Clear(); TIter next(&o.fEtaEDists); TObject* obj = 0; while ((obj = next())) fEtaEDists.Add(obj->Clone()); @@ -366,6 +381,7 @@ AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o) if (fEmpty) delete fEmpty; fEtaEDists.Delete(); if (fList) fList->Delete(); + fFits.Clear(); fEDist = static_cast(o.fEDist->Clone()); fEmpty = static_cast(o.fEmpty->Clone()); @@ -434,12 +450,15 @@ AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, //____________________________________________________________________ void -AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax, - Double_t deMax, Int_t nDeBins, - Bool_t incr) +AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, + Double_t emin, + Double_t emax, + Double_t deMax, + Int_t nDeBins, + Bool_t incr) { TH1D* h = 0; - TString name = Form("%s_etabin%03d", fName.Data(), ieta); + TString name = Form(fgkEDistFormat, fName.Data(), ieta); TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f " "(#eta bin %d)", fName.Data(), emin, emax, ieta); if (!incr) @@ -460,10 +479,9 @@ AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax, fEtaEDists.AddAt(h, ieta-1); } //____________________________________________________________________ -TH1D* -AliFMDEnergyFitter::RingHistos::MakePar(const char* name, - const char* title, - const TAxis& eta) const +TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name, + const char* title, + const TAxis& eta) const { TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), Form("%s for %s", title, fName.Data()), @@ -533,14 +551,20 @@ AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, } //____________________________________________________________________ TObjArray* -AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, - Double_t lowCut, UShort_t nLandau, - UShort_t minEntries, - UShort_t minusBins) const +AliFMDEnergyFitter::RingHistos::Fit(TList* dir, + const TAxis& eta, + Double_t lowCut, + UShort_t nParticles, + UShort_t minEntries, + UShort_t minusBins, + Double_t relErrorCut, + Double_t chi2nuCut) const { + // Get our ring list from the output container TList* l = GetOutputList(dir); if (!l) return 0; + // Get the energy distributions from the output container TList* dists = static_cast(l->FindObject("EDists")); if (!dists) { AliWarning(Form("Didn't find %s_EtaEDists in %s", @@ -549,24 +573,28 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, return 0; } - TObjArray* pars = new TObjArray(3+nLandau+1); + // Container of the fit results histograms + TObjArray* pars = new TObjArray(3+nParticles+1); pars->SetName("FitResults"); l->Add(pars); - TH1* hChi2 = 0; - TH1* hC = 0; - TH1* hMpv = 0; - TH1* hW = 0; - TH1* hS = 0; - TH1* hN = 0; - TH1* hA[nLandau-1]; - pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta)); - pars->Add(hC = MakePar("c", "Constant", eta)); - pars->Add(hMpv = MakePar("mpv", "#Delta_{p}", eta)); - pars->Add(hW = MakePar("w", "#xi", eta)); - pars->Add(hS = MakePar("s", "#sigma", eta)); - pars->Add(hN = MakePar("n", "N", eta)); - for (UShort_t i = 1; i < nLandau; i++) + // Result objects stored in separate list on the output + TH1* hChi2 = 0; + TH1* hC = 0; + TH1* hDelta = 0; + TH1* hXi = 0; + TH1* hSigma = 0; + TH1* hSigmaN = 0; + TH1* hN = 0; + TH1* hA[nParticles-1]; + pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta)); + pars->Add(hC = MakePar("c", "Constant", eta)); + pars->Add(hDelta = MakePar("delta", "#Delta_{p}", eta)); + pars->Add(hXi = MakePar("xi", "#xi", eta)); + pars->Add(hSigma = MakePar("sigma", "#sigma", eta)); + pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}", eta)); + pars->Add(hN = MakePar("n", "N", eta)); + for (UShort_t i = 1; i < nParticles; i++) pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta)); @@ -575,62 +603,125 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, Int_t high = 0; for (Int_t i = 0; i < nDists; i++) { TH1D* dist = static_cast(dists->At(i)); + // Ignore empty histograms altoghether + if (!dist || dist->GetEntries() <= 0) continue; - if (!dist) continue; + // Scale to the bin-width + dist->Scale(1., "width"); + + // Normalize peak to 1 + Double_t max = dist->GetMaximum(); + if (max <= 0) continue; + dist->Scale(1/max); + + // Check that we have enough entries + if (dist->GetEntries() <= minEntries) { + AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)", + i, dist->GetName(), Int_t(dist->GetEntries()), + minEntries)); + continue; + } - TF1* res = FitHist(dist,lowCut,nLandau,minEntries,minusBins); + // Now fit + TF1* res = FitHist(dist,lowCut,nParticles,minusBins, + relErrorCut,chi2nuCut); if (!res) continue; - + // dist->GetListOfFunctions()->Add(res); + + // Store eta limits low = TMath::Min(low,i+1); high = TMath::Max(high,i+1); + // Get the reduced chi square Double_t chi2 = res->GetChisquare(); Int_t ndf = res->GetNDF(); - hChi2->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0); - hC ->SetBinContent(i+1, res->GetParameter(0)); - hMpv->SetBinContent(i+1, res->GetParameter(1)); - hW ->SetBinContent(i+1, res->GetParameter(2)); - hN ->SetBinContent(i+1, res->GetParameter(3)); - - hC ->SetBinError(i+1, res->GetParError(1)); - hMpv->SetBinError(i+1, res->GetParError(2)); - hW ->SetBinError(i+1, res->GetParError(2)); - // hN ->SetBinError(i, res->GetParError(3)); - - for (UShort_t j = 0; j < nLandau-1; j++) { - hA[j]->SetBinContent(i+1, res->GetParameter(4+j)); - hA[j]->SetBinError(i+1, res->GetParError(4+j)); + + // Store results of best fit in output histograms + res->SetLineWidth(4); + hChi2 ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0); + hC ->SetBinContent(i+1, res->GetParameter(kC)); + hDelta ->SetBinContent(i+1, res->GetParameter(kDelta)); + hXi ->SetBinContent(i+1, res->GetParameter(kXi)); + hSigma ->SetBinContent(i+1, res->GetParameter(kSigma)); + hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN)); + hN ->SetBinContent(i+1, res->GetParameter(kN)); + + hC ->SetBinError(i+1, res->GetParError(kC)); + hDelta ->SetBinError(i+1, res->GetParError(kDelta)); + hXi ->SetBinError(i+1, res->GetParError(kXi)); + hSigma ->SetBinError(i+1, res->GetParError(kSigma)); + hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN)); + hN ->SetBinError(i+1, res->GetParError(kN)); + + for (UShort_t j = 0; j < nParticles-1; j++) { + hA[j]->SetBinContent(i+1, res->GetParameter(kA+j)); + hA[j]->SetBinError(i+1, res->GetParError(kA+j)); } } + // Fit the full-ring histogram TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data())); - if (total) { - TF1* res = FitHist(total,lowCut,nLandau,minEntries,minusBins); + if (total && total->GetEntries() >= minEntries) { + + // Scale to the bin-width + total->Scale(1., "width"); + + // Normalize peak to 1 + Double_t max = total->GetMaximum(); + if (max > 0) total->Scale(1/max); + + TF1* res = FitHist(total,lowCut,nParticles,minusBins, + relErrorCut,chi2nuCut); if (res) { + // Make histograms for the result of this fit Double_t chi2 = res->GetChisquare(); Int_t ndf = res->GetNDF(); - pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high, + pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high, ndf > 0 ? chi2/ndf : 0, 0)); - pars->Add(MakeTotal("t_c", "Constant", eta, low, high, - res->GetParameter(0),res->GetParError(0))); - pars->Add(MakeTotal("t_mpv", "#Delta_{p}", eta, low, high, - res->GetParameter(1),res->GetParError(1))); - pars->Add(MakeTotal("t_w", "#xi", eta, low, high, - res->GetParameter(2),res->GetParError(2))); - pars->Add(MakeTotal("t_n", "N", eta, low, high, - res->GetParameter(3),0)); - for (UShort_t j = 1; j < nLandau; j++) - pars->Add(MakeTotal(Form("a%d_t",j+1), - Form("a_{%d}",j+1), eta, low, high, - res->GetParameter(3+j), res->GetParError(3+j))); + pars->Add(MakeTotal("t_c", "Constant", eta, low, high, + res->GetParameter(kC), + res->GetParError(kC))); + pars->Add(MakeTotal("t_delta", "#Delta_{p}", eta, low, high, + res->GetParameter(kDelta), + res->GetParError(kDelta))); + pars->Add(MakeTotal("t_xi", "#xi", eta, low, high, + res->GetParameter(kXi), + res->GetParError(kXi))); + pars->Add(MakeTotal("t_sigma", "#sigma", eta, low, high, + res->GetParameter(kSigma), + res->GetParError(kSigma))); + pars->Add(MakeTotal("t_sigman", "#sigma_{n}", eta, low, high, + res->GetParameter(kSigmaN), + res->GetParError(kSigmaN))); + pars->Add(MakeTotal("t_n", "N", eta, low, high, + res->GetParameter(kN),0)); + for (UShort_t j = 0; j < nParticles-1; j++) + pars->Add(MakeTotal(Form("t_a%d",j+2), + Form("a_{%d}",j+2), eta, low, high, + res->GetParameter(kA+j), + res->GetParError(kA+j))); } } + // Clean up list of histogram. Histograms with no entries or + // no functions are deleted. We have to do this using the TObjLink + // objects stored in the list since ROOT cannot guaranty the validity + // of iterators when removing from a list - tsck. Should just implement + // TIter::Remove(). TObjLink* lnk = dists->FirstLink(); while (lnk) { TH1* h = static_cast(lnk->GetObject()); - if (h->GetEntries() <= 0 || - h->GetListOfFunctions()->GetEntries() <= 0) { + bool remove = false; + if (h->GetEntries() <= 0) { + // AliWarning(Form("No entries in %s - removing", h->GetName())); + remove = true; + } + else if (h->GetListOfFunctions()->GetEntries() <= 0) { + // AliWarning(Form("No fuctions associated with %s - removing", + // h->GetName())); + remove = true; + } + if (remove) { TObjLink* keep = lnk->Next(); dists->Remove(lnk); lnk = keep; @@ -644,243 +735,195 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, //____________________________________________________________________ TF1* -AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,Double_t lowCut, - UShort_t nLandau, - UShort_t minEntries, - UShort_t minusBins) const +AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist, + Double_t lowCut, + UShort_t nParticles, + UShort_t minusBins, + Double_t relErrorCut, + Double_t chi2nuCut) const { Double_t maxRange = 10; - if (dist->GetEntries() <= minEntries) return 0; - - // Find the fit range - dist->GetXaxis()->SetRangeUser(lowCut, maxRange); - - // Normalize peak to 1 - Double_t max = dist->GetMaximum(); - dist->Scale(1/max); + // Create a fitter object + AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins); + f.Clear(); - // Get the bin with maximum - Int_t maxBin = dist->GetMaximumBin(); - Double_t maxE = dist->GetBinLowEdge(maxBin); - // Get the low edge - dist->GetXaxis()->SetRangeUser(lowCut, maxE); - Int_t minBin = maxBin - minusBins; // dist->GetMinimumBin(); - Double_t minE = TMath::Max(dist->GetBinCenter(minBin),lowCut); - Double_t maxEE = dist->GetBinCenter(maxBin+2*minusBins); - - // Restore the range - dist->GetXaxis()->SetRangeUser(0, maxRange); - - // First do a single landau fit - TF1* landau1 = new TF1("landau1", "landau", minE, maxEE); - // landau1->SetParameters(1,1,1,1); - landau1->SetParNames("C","#Delta_{p}","#xi"); - landau1->SetParLimits(1,minE,maxRange); - landau1->SetParLimits(2,0,maxRange); - landau1->SetLineColor(Color()); - // Tight fit around peak - make sure we get that right. - TFitResultPtr r = dist->Fit(landau1, "RQNS", "", minE, maxEE); - - return FitMore2(dist, nLandau, *r, landau1, minE, maxRange); -} - -//____________________________________________________________________ -TF1* -AliFMDEnergyFitter::RingHistos::FitMore(TH1* dist, - UShort_t nLandau, - TFitResult& r, - TF1* landau1, - Double_t minE, - Double_t maxRange) const -{ - static TClonesArray res("TFitResult"); - static TObjArray funcs; - res.Clear(); - funcs.SetOwner(); - funcs.Clear(); - Int_t nRes = 0; - - //r->Print(); - new(res[nRes++]) TFitResult(r); - funcs.AddAtAndExpand(landau1, 0); - - // Now try to fit - for (UShort_t n = 2; n <= nLandau; n++) { - TFitResult* rr = static_cast(res.At(nRes-1)); - if (!rr) { - AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n)); - return 0; - } - // Create the function object - Double_t mpvi = rr->Parameter(1); - Double_t wi = rr->Parameter(2); - Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi; - TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,3+n); - landaui->SetLineStyle((n % 10)+1); - landaui->SetLineWidth(2); - landaui->SetLineColor(Color()); - landaui->SetParameter(0, rr->Parameter(0)); - landaui->SetParameter(1, rr->Parameter(1)); - landaui->SetParameter(2, rr->Parameter(2)); - landaui->SetParLimits(1, minE, maxRange); - landaui->SetParLimits(2,0,maxRange); - landaui->FixParameter(3, n); - landaui->SetParNames("C","#Delta_{p}","#xi", "N"); - for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit - landaui->SetParameter(N_A(i), rr->Parameter(N_A(i))); - landaui->SetParLimits(N_A(i), 0,1); - landaui->SetParName(i, Form("a_{%d}", i)); - } - landaui->SetParameter(N_A(n), (n == 2 ? 0.05 : rr->Parameter(N_A(n-1))/5)); - landaui->SetParLimits(N_A(n), 0, 1); - landaui->SetParName(N_A(n), Form("a_{%d}", n)); - // landaui->Print(); - - TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi); - // tr->Print(); - if (CheckResult(*tr)) { - new(res[nRes++]) TFitResult(*tr); - funcs.AddAtAndExpand(landaui,n-1); - continue; - } - // Stop on bad fit - break; + // If we are only asked to fit a single particle, return this fit, + // no matter what. + if (nParticles == 1) { + TF1* r = f.Fit1Particle(dist, 0); + if (!r) return 0; + return new TF1(*r); } - TF1* ret = 0; - if (funcs.At(nRes-1)) ret = static_cast(funcs.At(nRes-1)->Clone()); - dist->GetListOfFunctions()->Add(ret); - res.Clear(); - funcs.Clear(); - - return ret; -} - -//____________________________________________________________________ -TF1* -AliFMDEnergyFitter::RingHistos::FitMore2(TH1* dist, - UShort_t nLandau, - TFitResult& r, - TF1* landau1, - Double_t minE, - Double_t maxRange) const -{ - static TClonesArray res("TFitResult"); - static TObjArray funcs; - res.Clear(); - funcs.SetOwner(); - funcs.Clear(); - Int_t nRes = 0; - - //r->Print(); - new(res[nRes++]) TFitResult(r); - funcs.AddAtAndExpand(landau1, 0); - - // Now try to fit - for (UShort_t n = 2; n <= nLandau; n++) { - TFitResult* rr = static_cast(res.At(nRes-1)); - if (!rr) { - AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n)); - return 0; + // Fit from 2 upto n particles + for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0); + + + // Now, we need to select the best fit + Int_t nFits = f.fFitResults.GetEntriesFast(); + TF1* good[nFits]; + for (Int_t i = nFits-1; i >= 0; i--) { + good[i] = 0; + TF1* ff = static_cast(f.fFunctions.At(i)); + // ff->SetLineColor(Color()); + ff->SetRange(0, maxRange); + dist->GetListOfFunctions()->Add(new TF1(*ff)); + if (CheckResult(static_cast(f.fFitResults.At(i)), + relErrorCut, chi2nuCut)) { + good[i] = ff; + ff->SetLineWidth(2); + // f.fFitResults.At(i)->Print("V"); } - // Create the function object - Double_t mpvi = rr->Parameter(1); - Double_t wi = rr->Parameter(2); - Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi; - TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi, - n*3+1); - landaui->SetLineStyle((n % 10)+1); - landaui->SetLineWidth(2); - landaui->SetLineColor(Color()); - landaui->SetParameter(0, rr->Parameter(0)); - landaui->SetParameter(1, rr->Parameter(1)); - landaui->SetParameter(2, rr->Parameter(2)); - landaui->SetParLimits(1, minE, maxRange); - landaui->SetParLimits(2,0,maxRange); - landaui->FixParameter(3, n); - landaui->SetParNames("C","#Delta_{p}","#xi", "N"); - for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit - landaui->SetParameter(N2_A(i), rr->Parameter(N2_A(i))); - landaui->SetParameter(N2_D(i), rr->Parameter(N2_D(i))); - landaui->SetParameter(N2_X(i), rr->Parameter(N2_X(i))); - landaui->SetParLimits(N2_A(i), 0,1); - landaui->SetParLimits(N2_D(i), minE,maxEi); - landaui->SetParLimits(N2_X(i), 0,maxRange); - landaui->SetParName(N2_A(i), Form("a_{%d}", i)); - landaui->SetParName(N2_D(i), Form("#Delta_{p,%d}", i)); - landaui->SetParName(N2_X(i), Form("#xi_{%d}", i)); - } - landaui->SetParameter(N2_A(n), n == 2 ? 0.05 : rr->Parameter(N2_A(n-1))/5); - landaui->SetParameter(N2_D(n), mpvi); - landaui->SetParameter(N2_X(n), wi); - landaui->SetParLimits(N2_A(n), 0,1); - landaui->SetParLimits(N2_D(n), minE,maxEi); - landaui->SetParLimits(N2_X(n), 0,maxRange); - landaui->SetParName(N2_A(n), Form("a_{%d}", n)); - landaui->SetParName(N2_D(n), Form("#Delta_{p,%d}", n)); - landaui->SetParName(N2_X(n), Form("#xi_{%d}", n)); - if (fDebug > 2) landaui->Print(); - - TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi); - if (fDebug > 2) tr->Print(); - if (CheckResult(*tr)) { - new(res[nRes++]) TFitResult(*tr); - funcs.AddAtAndExpand(landaui,n-1); - continue; + } + // If all else fails, use the 1 particle fit + TF1* ret = static_cast(f.fFunctions.At(0)); + + // Find the fit with the most valid particles + for (Int_t i = nFits-1; i > 0; i--) { + if (!good[i]) continue; + if (fDebug > 1) { + AliInfo(Form("Choosing fit with n=%d", i+1)); + f.fFitResults.At(i)->Print(); } - // Stop on bad fit + ret = good[i]; break; } - TF1* ret = 0; - if (funcs.At(nRes-1)) ret = static_cast(funcs.At(nRes-1)->Clone()); - dist->GetListOfFunctions()->Add(ret); - - res.Clear(); - funcs.Clear(); - - return ret; + // Give a warning if we're using fall-back + if (ret == f.fFunctions.At(0)) { + AliWarning("Choosing fall-back 1 particle fit"); + } + // Copy our result and return (the functions are owned by the fitter) + TF1* fret = new TF1(*ret); + return fret; } //____________________________________________________________________ Bool_t -AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult& r) const +AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, + Double_t relErrorCut, + Double_t chi2nuCut) const { - Double_t chi2 = r.Chi2(); - Int_t ndf = r.Ndf(); + if (fDebug > 10) r->Print(); + TString n = r->GetName(); + n.ReplaceAll("TFitResult-", ""); + Double_t chi2 = r->Chi2(); + Int_t ndf = r->Ndf(); // Double_t prob = r.Prob(); - if (ndf <= 0 || chi2 / ndf > 5) { - if (fDebug > 2) - AliWarning(Form("Fit %s has chi^2/ndf=%f/%d=%f", - r.GetName(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf))); - return kFALSE; + Bool_t ret = kTRUE; + + // Check that the reduced chi square isn't larger than 5 + if (ndf <= 0 || chi2 / ndf > chi2nuCut) { + if (fDebug > 2) { + AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f", + n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf), + chi2nuCut)); } + ret = kFALSE; } - UShort_t nPar = r.NPar(); + // Check each parameter + UShort_t nPar = r->NPar(); for (UShort_t i = 0; i < nPar; i++) { - if (i == 3) continue; + if (i == kN) continue; // Skip the number parameter - Double_t v = r.Parameter(i); - Double_t e = r.ParError(i); + // Get value and error and check value + Double_t v = r->Parameter(i); + Double_t e = r->ParError(i); if (v == 0) continue; - if (v == 0 || e / v > 0.2) { - if (fDebug > 2) - AliWarning(Form("Fit %s has Delta %s/%s=%f/%f=%f%%", - r.GetName(), r.ParName(i).c_str(), - r.ParName(i).c_str(), e, v, 100*(v == 0 ? 0 : e/v))); - return kFALSE; + + // Calculate the relative error and check it + Double_t re = e / v; + Double_t cut = relErrorCut * (i < kN ? 1 : 2); + if (re > cut) { + if (fDebug > 2) { + AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%", + n.Data(), r->ParName(i).c_str(), + r->ParName(i).c_str(), e, v, + 100*(v == 0 ? 0 : e/v), + 100*(cut))); } + ret = kFALSE; } } - if (nPar > 4) { - if (r.Parameter(nPar-1) <= 1e-10) { - if (fDebug) - AliWarning(Form("Last scale %s is too small %f<1e-10", - r.ParName(nPar-1).c_str(), r.Parameter(nPar-1))); - return kFALSE; + + // Check if we have scale parameters + if (nPar > kN) { + + // Check that the last particle has a significant contribution + Double_t lastScale = r->Parameter(nPar-1); + if (lastScale <= 1e-7) { + if (fDebug) { + AliWarning(Form("%s: %s=%9.6f<1e-7", + n.Data(), r->ParName(nPar-1).c_str(), lastScale)); } + ret = kFALSE; } } - return kTRUE; + return ret; +} + + +//__________________________________________________________________ +void +AliFMDEnergyFitter::RingHistos::FindBestFits(TList* d, + AliFMDCorrELossFit& obj, + const TAxis& eta, + Double_t relErrorCut, + Double_t chi2nuCut, + Double_t minWeightCut) +{ + // Get our ring list from the output container + TList* l = GetOutputList(d); + if (!l) return; + + // Get the energy distributions from the output container + TList* dists = static_cast(l->FindObject("EDists")); + if (!dists) { + AliWarning(Form("Didn't find %s_EtaEDists in %s", + fName.Data(), l->GetName())); + l->ls(); + return; + } + Int_t nBin = eta.GetNbins(); + + for (Int_t b = 1; b <= nBin; b++) { + TString n(Form(fgkEDistFormat, fName.Data(), b)); + TH1D* dist = static_cast(dists->FindObject(n)); + // Ignore empty histograms altoghether + if (!dist || dist->GetEntries() <= 0) continue; + + AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist, + relErrorCut, + chi2nuCut, + minWeightCut); + best->fDet = fDet; + best->fRing = fRing; + best->fBin = b; // + // Double_t eta = fAxis->GetBinCenter(b); + obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best)); + } +} + +//__________________________________________________________________ +AliFMDCorrELossFit::ELossFit* +AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist, + Double_t relErrorCut, + Double_t chi2nuCut, + Double_t minWeightCut) +{ + TList* funcs = dist->GetListOfFunctions(); + TIter next(funcs); + TF1* func = 0; + fFits.Clear(); + Int_t i = 0; + // Info("FindBestFit", "%s", dist->GetName()); + while ((func = static_cast(next()))) { + AliFMDCorrELossFit::ELossFit* fit = + new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func); + fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut); + } + fFits.Sort(false); + // fFits.Print(); + return static_cast(fFits.At(0)); }