X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG2%2FFORWARD%2Fanalysis2%2FAliFMDEnergyFitter.cxx;h=8c7eb80e004ca4aa33c8106c21492eecc8d33776;hb=f4a1b795965b59e87ae7c5b8b62485b592600a52;hp=219d925207681a518000cd9cb221874c735c5a8b;hpb=c389303ef9eccf425db64297088ff00c303389bf;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx b/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx index 219d9252076..8c7eb80e004 100644 --- a/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx @@ -2,61 +2,85 @@ // 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 +namespace { + const char* fgkEDistFormat = "%s_etabin%03d"; +} //____________________________________________________________________ AliFMDEnergyFitter::AliFMDEnergyFitter() : TNamed(), fRingHistos(), - fLowCut(0.3), - fNParticles(3), - fMinEntries(100), + fLowCut(0.4), + fNParticles(5), + fMinEntries(1000), fFitRangeBinWidth(4), fDoFits(false), + fDoMakeObject(false), fEtaAxis(), + fCentralityAxis(), fMaxE(10), fNEbins(300), fUseIncreasingBins(true), fMaxRelParError(.25), fMaxChi2PerNDF(20), + fMinWeight(1e-7), fDebug(0) -{} +{ + // + // Default Constructor - do not use + // +} //____________________________________________________________________ AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title) : TNamed("fmdEnergyFitter", title), fRingHistos(), - fLowCut(0.3), - fNParticles(3), - fMinEntries(100), + fLowCut(0.4), + fNParticles(5), + fMinEntries(1000), fFitRangeBinWidth(4), fDoFits(false), + fDoMakeObject(false), fEtaAxis(0,0,0), + fCentralityAxis(0,0,0), fMaxE(10), fNEbins(300), fUseIncreasingBins(true), fMaxRelParError(.25), - fMaxChi2PerNDF(20), + fMaxChi2PerNDF(20), + fMinWeight(1e-7), fDebug(3) { + // + // Constructor + // + // Parameters: + // title Title of object - not significant + // fEtaAxis.SetName("etaAxis"); fEtaAxis.SetTitle("#eta"); + fCentralityAxis.SetName("centralityAxis"); + fCentralityAxis.SetTitle("Centrality [%]"); fRingHistos.Add(new RingHistos(1, 'I')); fRingHistos.Add(new RingHistos(2, 'I')); fRingHistos.Add(new RingHistos(2, 'O')); @@ -73,14 +97,23 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o) fMinEntries(o.fMinEntries), fFitRangeBinWidth(o.fFitRangeBinWidth), fDoFits(o.fDoFits), + fDoMakeObject(o.fDoMakeObject), fEtaAxis(o.fEtaAxis), + fCentralityAxis(o.fCentralityAxis), fMaxE(o.fMaxE), fNEbins(o.fNEbins), fUseIncreasingBins(o.fUseIncreasingBins), fMaxRelParError(o.fMaxRelParError), fMaxChi2PerNDF(o.fMaxChi2PerNDF), + fMinWeight(o.fMinWeight), fDebug(o.fDebug) { + // + // Copy constructor + // + // Parameters: + // o Object to copy from + // TIter next(&o.fRingHistos); TObject* obj = 0; while ((obj = next())) fRingHistos.Add(obj); @@ -89,6 +122,9 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o) //____________________________________________________________________ AliFMDEnergyFitter::~AliFMDEnergyFitter() { + // + // Destructor + // fRingHistos.Delete(); } @@ -96,6 +132,15 @@ AliFMDEnergyFitter::~AliFMDEnergyFitter() AliFMDEnergyFitter& AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o) { + // + // Assignment operator + // + // Parameters: + // o Object to assign from + // + // Return: + // Reference to this + // TNamed::operator=(o); fLowCut = o.fLowCut; @@ -103,8 +148,22 @@ AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o) fMinEntries = o.fMinEntries; fFitRangeBinWidth= o.fFitRangeBinWidth; fDoFits = o.fDoFits; - fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()); + fDoMakeObject = o.fDoMakeObject; + fEtaAxis.Set(o.fEtaAxis.GetNbins(), + o.fEtaAxis.GetXmin(), + o.fEtaAxis.GetXmax()); + if (o.fCentralityAxis.GetXbins()) { + const TArrayD* bins = o.fCentralityAxis.GetXbins(); + fCentralityAxis.Set(bins->GetSize()-1,bins->GetArray()); + } + else + fCentralityAxis.Set(o.fCentralityAxis.GetNbins(), + o.fCentralityAxis.GetXmin(), + o.fCentralityAxis.GetXmax()); fDebug = o.fDebug; + fMaxRelParError= o.fMaxRelParError; + fMaxChi2PerNDF = o.fMaxChi2PerNDF; + fMinWeight = o.fMinWeight; fRingHistos.Delete(); TIter next(&o.fRingHistos); @@ -118,6 +177,16 @@ AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o) AliFMDEnergyFitter::RingHistos* AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const { + // + // Get the ring histogram container + // + // Parameters: + // d Detector + // r Ring + // + // Return: + // Ring histogram container + // Int_t idx = -1; switch (d) { case 1: idx = 0; break; @@ -133,32 +202,90 @@ AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const void AliFMDEnergyFitter::Init(const TAxis& eAxis) { + // + // Initialise the task + // + // Parameters: + // etaAxis The eta axis to use. Note, that if the eta axis + // has already been set (using SetEtaAxis), then this parameter will be + // ignored + // if (fEtaAxis.GetNbins() == 0 || - fEtaAxis.GetXmin() == fEtaAxis.GetXmax()) + TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7) SetEtaAxis(eAxis); + if (fCentralityAxis.GetNbins() == 0) { + UShort_t n = 12; + Double_t bins[] = { 0., 5., 10., 15., 20., 30., + 40., 50., 60., 70., 80., 100. }; + SetCentralityAxis(n, bins); + } TIter next(&fRingHistos); RingHistos* o = 0; while ((o = static_cast(next()))) - o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins); + o->Init(fEtaAxis, fCentralityAxis, fMaxE, fNEbins, fUseIncreasingBins); } //____________________________________________________________________ void AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis) { + // + // Set the eta axis to use. This will force the code to use this + // eta axis definition - irrespective of whatever axis is passed to + // the Init member function. Therefore, this member function can be + // used to force another eta axis than one found in the correction + // objects. + // + // Parameters: + // etaAxis Eta axis to use + // SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax()); } //____________________________________________________________________ void AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax) { + // + // Set the eta axis to use. This will force the code to use this + // eta axis definition - irrespective of whatever axis is passed to + // the Init member function. Therefore, this member function can be + // used to force another eta axis than one found in the correction + // objects. + // + // Parameters: + // nBins Number of bins + // etaMin Minimum of the eta axis + // etaMax Maximum of the eta axis + // fEtaAxis.Set(nBins,etaMin,etaMax); } +//____________________________________________________________________ +void +AliFMDEnergyFitter::SetCentralityAxis(UShort_t n, Double_t* bins) +{ + fCentralityAxis.Set(n-1, bins); +} + //____________________________________________________________________ Bool_t -AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, +AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, + Double_t cent, Bool_t empty) { + // + // Fitter the input AliESDFMD object + // + // Parameters: + // input Input + // cent Centrality + // empty Whether the event is 'empty' + // + // Return: + // True on success, false otherwise + // + Int_t icent = fCentralityAxis.FindBin(cent); + if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0; + for(UShort_t d = 1; d <= 3; d++) { Int_t nRings = (d == 1 ? 1 : 2); for (UShort_t q = 0; q < nRings; q++) { @@ -181,7 +308,7 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, Int_t ieta = fEtaAxis.FindBin(eta1); if (ieta < 1 || ieta > fEtaAxis.GetNbins()) continue; - histos->Fill(empty, ieta-1, mult); + histos->Fill(empty, ieta-1, icent, mult); } // for strip } // for sector @@ -193,8 +320,14 @@ AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, //____________________________________________________________________ void -AliFMDEnergyFitter::Fit(TList* dir) +AliFMDEnergyFitter::Fit(const TList* dir) { + // + // Scale the histograms to the total number of events + // + // Parameters: + // dir Where the histograms are + // if (!fDoFits) return; TList* d = static_cast(dir->FindObject(GetName())); @@ -229,12 +362,55 @@ AliFMDEnergyFitter::Fit(TList* dir) stack[i % nStack]->Add(static_cast(l->At(i))); } } + + if (!fDoMakeObject) return; + + MakeCorrectionsObject(d); +} + +//____________________________________________________________________ +void +AliFMDEnergyFitter::MakeCorrectionsObject(TList* d) +{ + // + // Generate the corrections object + // + // Parameters: + // dir List to analyse + // + 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(new TNamed("filename", fileName.Data())); + d->Add(obj, oName.Data()); } //____________________________________________________________________ void AliFMDEnergyFitter::DefineOutput(TList* dir) { + // + // Define the output histograms. These are put in a sub list of the + // passed list. The histograms are merged before the parent task calls + // AliAnalysisTaskSE::Terminate + // + // Parameters: + // dir Directory to add to + // TList* d = new TList; d->SetName(GetName()); dir->Add(d); @@ -249,12 +425,50 @@ AliFMDEnergyFitter::DefineOutput(TList* dir) void AliFMDEnergyFitter::SetDebug(Int_t dbg) { + // + // Set the debug level. The higher the value the more output + // + // Parameters: + // dbg Debug level + // fDebug = dbg; TIter next(&fRingHistos); RingHistos* o = 0; while ((o = static_cast(next()))) o->fDebug = dbg; } + +//____________________________________________________________________ +void +AliFMDEnergyFitter::Print(Option_t*) const +{ + // + // Print information + // + // Parameters: + // option Not used + // + char ind[gROOT->GetDirLevel()+1]; + for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; + ind[gROOT->GetDirLevel()] = '\0'; + std::cout << ind << ClassName() << ": " << 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() @@ -263,8 +477,13 @@ AliFMDEnergyFitter::RingHistos::RingHistos() fEmpty(0), fEtaEDists(), fList(0), + fFits("AliFMDCorrELossFit::ELossFit"), fDebug(0) -{} +{ + // + // Default CTOR + // +} //____________________________________________________________________ AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r) @@ -273,8 +492,16 @@ AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r) fEmpty(0), fEtaEDists(), fList(0), + fFits("AliFMDCorrELossFit::ELossFit"), fDebug(0) { + // + // Constructor + // + // Parameters: + // d detector + // r ring + // fEtaEDists.SetName("EDists"); } //____________________________________________________________________ @@ -284,8 +511,16 @@ AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o) fEmpty(o.fEmpty), fEtaEDists(), fList(0), + fFits("AliFMDCorrELossFit::ELossFit"), fDebug(0) { + // + // Copy constructor + // + // Parameters: + // o Object to copy from + // + fFits.Clear(); TIter next(&o.fEtaEDists); TObject* obj = 0; while ((obj = next())) fEtaEDists.Add(obj->Clone()); @@ -301,12 +536,22 @@ AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o) AliFMDEnergyFitter::RingHistos& AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o) { + // + // Assignment operator + // + // Parameters: + // o Object to assign from + // + // Return: + // Reference to this + // AliForwardUtil::RingHistos::operator=(o); if (fEDist) delete fEDist; if (fEmpty) delete fEmpty; fEtaEDists.Delete(); if (fList) fList->Delete(); + fFits.Clear(); fEDist = static_cast(o.fEDist->Clone()); fEmpty = static_cast(o.fEmpty->Clone()); @@ -327,20 +572,36 @@ AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o) //____________________________________________________________________ AliFMDEnergyFitter::RingHistos::~RingHistos() { + // + // Destructor + // if (fEDist) delete fEDist; fEtaEDists.Delete(); } //____________________________________________________________________ void -AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult) +AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, + Int_t /* icent */, Double_t mult) { + // + // Fill histogram + // + // Parameters: + // empty True if event is empty + // ieta Eta bin (0-based) + // icent Centrality bin (1-based) + // mult Signal + // if (empty) { fEmpty->Fill(mult); return; } fEDist->Fill(mult); + // Here, we should first look up in a centrality array, and then in + // that array look up the eta bin - or perhaps we should do it the + // other way around? if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return; TH1D* h = static_cast(fEtaEDists.At(ieta)); @@ -354,11 +615,18 @@ TArrayD AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const { - // Service function to define a logarithmic axis. - // Parameters: - // n Number of bins - // min Minimum of axis - // max Maximum of axis + // + // Make an axis with increasing bins + // + // Parameters: + // n Number of bins + // min Minimum + // max Maximum + // + // Return: + // An axis with quadratically increasing bin size + // + TArrayD bins(n+1); Double_t dx = (max-min) / n; bins[0] = min; @@ -382,8 +650,19 @@ AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Int_t nDeBins, Bool_t incr) { + // + // Make E/E_mip histogram + // + // Parameters: + // ieta Eta bin + // eMin Least signal + // eMax Largest signal + // deMax Maximum energy loss + // nDeBins Number energy loss bins + // incr Whether to make bins of increasing size + // 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) @@ -408,6 +687,17 @@ TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name, const char* title, const TAxis& eta) const { + // + // Make a parameter histogram + // + // Parameters: + // name Name of histogram. + // title Title of histogram. + // eta Eta axis + // + // Return: + // + // TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), Form("%s for %s", title, fName.Data()), eta.GetNbins(), eta.GetXmin(), eta.GetXmax()); @@ -432,6 +722,21 @@ AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, Double_t val, Double_t err) const { + // + // Make a histogram that contains the results of the fit over the full ring + // + // Parameters: + // name Name + // title Title + // eta Eta axis + // low Least bin + // high Largest bin + // val Value of parameter + // err Error on parameter + // + // Return: + // The newly allocated histogram + // Double_t xlow = eta.GetBinLowEdge(low); Double_t xhigh = eta.GetBinUpEdge(high); TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), @@ -455,9 +760,19 @@ AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, //____________________________________________________________________ void AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, + const TAxis& /* cAxis */, Double_t maxDE, Int_t nDEbins, Bool_t useIncrBin) { + // + // Initialise object + // + // Parameters: + // eAxis Eta axis + // maxDE Max energy loss to consider + // nDEbins Number of bins + // useIncrBin Whether to use an increasing bin size + // fEDist = new TH1D(Form("%s_edist", fName.Data()), Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 200, 0, 6); @@ -466,24 +781,97 @@ AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, fName.Data()), 200, 0, 6); fList->Add(fEDist); fList->Add(fEmpty); + // Here, we should make an array of centrality bins ... + // In that case, the structure will be + // + // -+- Ring1 -+- Centrality1 -+- Eta1 + // | | +- Eta2 + // ... ... ... + // | | +- EtaM + // | +- Centrality2 -+- Eta1 + // ... ... ... + // | | +- EtaM + // ... ... + // | +- CentralityN -+- Eta1 + // ... ... ... + // | +- EtaM + // -+- Ring2 -+- Centrality1 -+- Eta1 + // | | +- Eta2 + // ... ... ... + // | | +- EtaM + // | +- Centrality2 -+- Eta1 + // ... ... ... + // | | +- EtaM + // ... ... + // | +- CentralityN -+- Eta1 + // ... ... ... + // | +- EtaM + // ... ... ... + // -+- RingO -+- Centrality1 -+- Eta1 + // | +- Eta2 + // ... ... + // | +- EtaM + // +- Centrality2 -+- Eta1 + // ... ... + // | +- EtaM + // ... + // +- CentralityN -+- Eta1 + // ... ... + // +- EtaM + // // fEtaEDists.Expand(eAxis.GetNbins()); for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { Double_t min = eAxis.GetBinLowEdge(i); Double_t max = eAxis.GetBinUpEdge(i); + // Or may we should do it here? In that case we'd have the + // following structure: + // Ring1 -+- Eta1 -+- Centrality1 + // | +- Centrality2 + // ... ... + // | +- CentralityN + // +- Eta2 -+- Centrality1 + // | +- Centrality2 + // ... ... + // | +- CentralityN + // ... + // +- EtaM -+- Centrality1 + // +- Centrality2 + // ... + // +- CentralityN Make(i, min, max, maxDE, nDEbins, useIncrBin); } fList->Add(&fEtaEDists); + // fEtaEDists.ls(); } //____________________________________________________________________ TObjArray* -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 +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 { + // + // Fit each histogram to up to @a nParticles particle responses. + // + // Parameters: + // dir Output list + // eta Eta axis + // lowCut Lower cut + // nParticles Max number of convolved landaus to fit + // minEntries Minimum number of entries + // minusBins Number of bins from peak to subtract to + // get the fit range + // relErrorCut Cut applied to relative error of parameter. + // Note, for multi-particle weights, the cut + // is loosend by a factor of 2 + // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - + // the reduced @f$\chi^2@f$ + // + // Get our ring list from the output container TList* l = GetOutputList(dir); if (!l) return 0; @@ -535,17 +923,23 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, // 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) continue; + 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; + } // 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); @@ -580,6 +974,14 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, // Fit the full-ring histogram TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data())); 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) { @@ -621,8 +1023,17 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta, 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; @@ -643,6 +1054,31 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist, Double_t relErrorCut, Double_t chi2nuCut) const { + // + // Fit a signal histogram. First, the bin @f$ b_{min}@f$ with + // maximum bin content in the range @f$ [E_{min},\infty]@f$ is + // found. Then the fit range is set to the bin range + // @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1 + // particle signal is fitted to that. The parameters of that fit + // is then used as seeds for a fit of the @f$ N@f$ particle response + // to the data in the range + // @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$ + // + // Parameters: + // dist Histogram to fit + // lowCut Lower cut @f$ E_{min}@f$ on signal + // nParticles Max number @f$ N@f$ of convolved landaus to fit + // minusBins Number of bins @f$ \Delta b@f$ from peak to + // subtract to get the fit range + // relErrorCut Cut applied to relative error of parameter. + // Note, for multi-particle weights, the cut + // is loosend by a factor of 2 + // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - + // the reduced @f$\chi^2@f$ + // + // Return: + // The best fit function + // Double_t maxRange = 10; // Create a fitter object @@ -655,7 +1091,9 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist, if (nParticles == 1) { TF1* r = f.Fit1Particle(dist, 0); if (!r) return 0; - return new TF1(*r); + TF1* ret = new TF1(*r); + dist->GetListOfFunctions()->Add(ret); + return ret; } // Fit from 2 upto n particles @@ -663,15 +1101,15 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist, // Now, we need to select the best fit - Int_t nFits = f.fFitResults.GetEntriesFast(); + Int_t nFits = f.GetFitResults().GetEntriesFast(); TF1* good[nFits]; for (Int_t i = nFits-1; i >= 0; i--) { good[i] = 0; - TF1* ff = static_cast(f.fFunctions.At(i)); + TF1* ff = static_cast(f.GetFunctions().At(i)); // ff->SetLineColor(Color()); ff->SetRange(0, maxRange); dist->GetListOfFunctions()->Add(new TF1(*ff)); - if (CheckResult(static_cast(f.fFitResults.At(i)), + if (CheckResult(static_cast(f.GetFitResults().At(i)), relErrorCut, chi2nuCut)) { good[i] = ff; ff->SetLineWidth(2); @@ -679,22 +1117,22 @@ AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist, } } // If all else fails, use the 1 particle fit - TF1* ret = static_cast(f.fFunctions.At(0)); + TF1* ret = static_cast(f.GetFunctions().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(); + f.GetFitResults().At(i)->Print(); } ret = good[i]; break; } // Give a warning if we're using fall-back - if (ret == f.fFunctions.At(0)) + if (ret == f.GetFunctions().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; @@ -706,6 +1144,24 @@ AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, Double_t relErrorCut, Double_t chi2nuCut) const { + // + // Check the result of the fit. Returns true if + // - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$ + // - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters. Note, + // for multi-particle fits, this requirement is relaxed by a + // factor of 2 + // - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$ + // particle response + // + // Parameters: + // r Result to check + // relErrorCut Cut @f$ \delta_e@f$ applied to relative error + // of parameter. + // chi2nuCut Cut @f$ \max{\chi^2/\nu}@f$ + // + // Return: + // true if fit is good. + // if (fDebug > 10) r->Print(); TString n = r->GetName(); n.ReplaceAll("TFitResult-", ""); @@ -716,10 +1172,10 @@ AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, // Check that the reduced chi square isn't larger than 5 if (ndf <= 0 || chi2 / ndf > chi2nuCut) { - if (fDebug > 2) + 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)); + chi2nuCut)); } ret = kFALSE; } @@ -737,12 +1193,12 @@ AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, Double_t re = e / v; Double_t cut = relErrorCut * (i < kN ? 1 : 2); if (re > cut) { - if (fDebug > 2) + 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))); + 100*(cut))); } ret = kFALSE; } } @@ -753,9 +1209,9 @@ AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, // Check that the last particle has a significant contribution Double_t lastScale = r->Parameter(nPar-1); if (lastScale <= 1e-7) { - if (fDebug) + if (fDebug) { AliWarning(Form("%s: %s=%9.6f<1e-7", - n.Data(), r->ParName(nPar-1).c_str(), lastScale)); + n.Data(), r->ParName(nPar-1).c_str(), lastScale)); } ret = kFALSE; } } @@ -763,10 +1219,111 @@ AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r, } +//__________________________________________________________________ +void +AliFMDEnergyFitter::RingHistos::FindBestFits(const TList* d, + AliFMDCorrELossFit& obj, + const TAxis& eta, + Double_t relErrorCut, + Double_t chi2nuCut, + Double_t minWeightCut) +{ + // + // Find the best fits + // + // Parameters: + // d Parent list + // obj Object to add fits to + // eta Eta axis + // relErrorCut Cut applied to relative error of parameter. + // Note, for multi-particle weights, the cut + // is loosend by a factor of 2 + // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - + // the reduced @f$\chi^2@f$ + // minWeightCut Least valid @f$ a_i@f$ + // + + // 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(const TH1* dist, + Double_t relErrorCut, + Double_t chi2nuCut, + Double_t minWeightCut) +{ + // + // Find the best fit + // + // Parameters: + // dist Histogram + // relErrorCut Cut applied to relative error of parameter. + // Note, for multi-particle weights, the cut + // is loosend by a factor of 2 + // chi2nuCut Cut on @f$ \chi^2/\nu@f$ - + // the reduced @f$\chi^2@f$ + // minWeightCut Least valid @f$ a_i@f$ + // + // Return: + // Best fit + // + 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)); +} + + //____________________________________________________________________ void AliFMDEnergyFitter::RingHistos::Output(TList* dir) { + // + // Define outputs + // + // Parameters: + // dir + // fList = DefineOutputList(dir); }