--- /dev/null
+#ifndef ALILANDAUGAUSFITTER_H
+#define ALILANDAUGAUSFITTER_H
+/**
+ * @file AliLandauGausFitter.h
+ * @author Christian Holm Christensen <cholm@nbi.dk>
+ * @date Tue Mar 11 08:56:03 2014
+ *
+ * @brief Declaration and implementation of fitter of Landau-Gauss
+ * distributions to energy loss spectra.
+ *
+ * @ingroup pwglf_forward
+ */
+#include <AliLandauGaus.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TObjArray.h>
+#include <TString.h>
+#include <TArray.h>
+#include <TFitResult.h>
+#include <TError.h>
+
+/**
+ * Fit Landau-Gauss distributions to energy loss distributions
+ *
+ * @see AliLandauGaus
+ *
+ * @ingroup pwglf_forward
+ */
+class AliLandauGausFitter
+{
+public:
+ /**
+ * Enumeration of parameters
+ */
+ enum {
+ /** Index of pre-constant @f$ C@f$ */
+ kC = AliLandauGaus::kC,
+ /** Index of most probable value @f$ \Delta_p@f$ */
+ kDelta = AliLandauGaus::kDelta,
+ /** Index of Landau width @f$ \xi@f$ */
+ kXi = AliLandauGaus::kXi,
+ /** Index of Gaussian width @f$ \sigma@f$ */
+ kSigma = AliLandauGaus::kSigma,
+ /** Index of Gaussian additional width @f$ \sigma_n@f$ */
+ kSigmaN = AliLandauGaus::kSigmaN,
+ /** Index of Number of particles @f$ N@f$ */
+ kN = AliLandauGaus::kN,
+ /** Base index of particle strengths @f$ a_i@f$ for
+ @f$i=2,\ldots,N@f$ */
+ kA = AliLandauGaus::kA
+ };
+ /**
+ * Get the fit options to use
+ *
+ * @return Fit options used through-out
+ */
+ static const char* GetFitOptions() { return "RNS"; }
+ /**
+ * Constructor
+ *
+ * @param lowCut Lower cut of spectrum - data below this cuts is ignored
+ * @param maxRange Maximum range to fit to
+ * @param minusBins The number of bins below maximum to use
+ */
+ AliLandauGausFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins)
+ : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
+ fFitResults(0), fFunctions(0), fDebug(false)
+ {
+ fFitResults.SetOwner();
+ fFunctions.SetOwner();
+ }
+ /**
+ * Destructor
+ *
+ */
+ virtual ~AliLandauGausFitter()
+ {
+ fFitResults.Delete();
+ fFunctions.Delete();
+ }
+ /**
+ * Enable/disable debugging output.
+ *
+ * @param debug If true, enable debugging output
+ */
+ void SetDebug(Bool_t debug=true) { fDebug = debug; }
+ /**
+ * Clear internal arrays
+ *
+ */
+ void Clear()
+ {
+ fFitResults.Clear();
+ fFunctions.Clear();
+ }
+ /**
+ * Fit a 1-particle signal to the passed energy loss distribution
+ *
+ * Note that this function clears the internal arrays first
+ *
+ * @param dist Data to fit the function to
+ * @param sigman If larger than zero, the initial guess of the
+ * detector induced noise. If zero or less, then this
+ * parameter is ignored in the fit (fixed at 0)
+ *
+ * @return The function fitted to the data
+ */
+ TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
+ /**
+ * Fit a N-particle signal to the passed energy loss distribution
+ *
+ * If there's no 1-particle fit present, it does that first
+ *
+ * @param dist Data to fit the function to
+ * @param n Number of particle signals to fit
+ * @param sigman If larger than zero, the initial guess of the
+ * detector induced noise. If zero or less, then this
+ * parameter is ignored in the fit (fixed at 0)
+ *
+ * @return The function fitted to the data
+ */
+ TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
+ /**
+ * Fit a composite distribution of energy loss from both primaries
+ * and secondaries
+ *
+ * @param dist Distribution
+ * @param sigman If larger than zero, the initial guess of the
+ * detector included noise. If zero or less this
+ * parameter is fixed to 0.
+ *
+ * @return Function fitted to the data
+ */
+ TF1* FitComposite(TH1* dist, Double_t sigman);
+ /**
+ * Get Lower cut on data
+ *
+ * @return Lower cut on data
+ */
+ Double_t GetLowCut() const { return fLowCut; }
+ /**
+ * Get Maximum range to fit
+ *
+ * @return Maximum range to fit
+ */
+ Double_t GetMaxRange() const { return fMaxRange; }
+ /**
+ * Get Number of bins from maximum to fit 1st peak
+ *
+ * @return Number of bins from maximum to fit 1st peak
+ */
+ UShort_t GetMinusBins() const { return fMinusBins; }
+ /**
+ * Get Array of fit results
+ *
+ * @return Array of fit results
+ */
+ const TObjArray& GetFitResults() const { return fFitResults; }
+ /**
+ * Get Array of fit results
+ *
+ * @return Array of fit results
+ */
+ TObjArray& GetFitResults() { return fFitResults; }
+ /**
+ * Get Array of functions
+ *
+ * @return Array of functions
+ */
+ const TObjArray& GetFunctions() const { return fFunctions; }
+ /**
+ * Get Array of functions
+ *
+ * @return Array of functions
+ */
+ TObjArray& GetFunctions() { return fFunctions; }
+protected:
+ /**
+ * Set parameter limits on the function @a f. The limits are only
+ * set if @f$ low \le test \le high@f$
+ *
+ * @param f Function object
+ * @param iPar Parameter number
+ * @param test Initial guess
+ * @param low Low limit
+ * @param high high limit
+ */
+ void SetParLimits(TF1* f, Int_t iPar, Double_t test,
+ Double_t low, Double_t high)
+ {
+ if (test < low || test > high) {
+ ::Warning("","Initial value of %s=%f not in [%f,%f]",
+ f->GetParName(iPar), test, low, high);
+ return;
+ }
+ if (fDebug)
+ ::Info(/*"SetParLimits"*/"", "Set par limits on %-12s=%f: [%f,%f]",
+ f->GetParName(iPar), test, low, high);
+ f->SetParLimits(iPar, low, high);
+ }
+ const Double_t fLowCut; // Lower cut on data
+ const Double_t fMaxRange; // Maximum range to fit
+ const UShort_t fMinusBins; // Number of bins from maximum to fit 1st peak
+ TObjArray fFitResults; // Array of fit results
+ TObjArray fFunctions; // Array of functions
+ Bool_t fDebug; // Debug flag
+};
+
+
+//____________________________________________________________________
+inline TF1*
+AliLandauGausFitter::Fit1Particle(TH1* dist, Double_t sigman)
+{
+ // Clear the cache
+ Clear();
+
+ // Find the fit range
+ Int_t cutBin = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
+ Int_t maxBin = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
+ dist->GetNbinsX());
+ dist->GetXaxis()->SetRange(cutBin, maxBin);
+ // dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
+
+ // Get the bin with maximum
+ Int_t peakBin = dist->GetMaximumBin();
+ Double_t peakE = dist->GetBinLowEdge(peakBin);
+ Double_t rmsE = dist->GetRMS();
+
+ // Get the low edge
+ // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+ Int_t minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
+ Double_t minE = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
+ Double_t maxE = dist->GetBinCenter(peakBin+2*fMinusBins);
+
+ Int_t minEb = dist->GetXaxis()->FindBin(minE);
+ Int_t maxEb = dist->GetXaxis()->FindBin(maxE);
+ Double_t intg = dist->Integral(minEb, maxEb);
+ if (intg <= 0) {
+ ::Warning("Fit1Particle",
+ "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
+ dist->GetName(), minE, maxE, minEb, maxEb, intg);
+ return 0;
+ }
+
+ // Restore the range
+ dist->GetXaxis()->SetRange(1, maxBin);
+
+ // Define the function to fit
+ TF1* f = AliLandauGaus::MakeF1(intg,peakE,peakE/10,peakE/5,sigman,minE,maxE);
+ SetParLimits(f, kDelta, peakE, minE, fMaxRange);
+ SetParLimits(f, kXi, peakE, 0, rmsE); // 0.1
+ SetParLimits(f, kSigma, peakE/5, 1e-5, rmsE); // 0.1
+ if (sigman <= 0)
+ f->FixParameter(kSigmaN, 0);
+ else
+ SetParLimits(f, kSigmaN, peakE, 0, rmsE);
+
+
+ TString opts(Form("%s%s", GetFitOptions(), fDebug ? "" : "Q"));
+ // Do the fit, getting the result object
+ if (fDebug)
+ ::Info(/*"Fit1Particle"*/"", "Fitting in the range %f,%f", minE, maxE);
+ TFitResultPtr r = dist->Fit(f, opts, "", minE, maxE);
+ if (!r.Get()) {
+ ::Warning("Fit1Particle",
+ "No fit returned when processing %s in the range [%f,%f] "
+ "options %s", dist->GetName(), minE, maxE, opts.Data());
+ return 0;
+ }
+ // f->SetRange(minE, fMaxRange);
+ fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
+ fFunctions.AddAtAndExpand(f, 0);
+
+ return f;
+}
+//____________________________________________________________________
+inline TF1*
+AliLandauGausFitter::FitNParticle(TH1* dist, UShort_t n, Double_t sigman)
+{
+ // Get the seed fit result
+ TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
+ TF1* f = static_cast<TF1*>(fFunctions.At(0));
+ if (!r || !f) {
+ f = Fit1Particle(dist, sigman);
+ r = static_cast<TFitResult*>(fFitResults.At(0));
+ if (!r || !f) {
+ ::Warning("FitNLandau", "No first shot at landau fit");
+ return 0;
+ }
+ }
+
+ // Get some parameters from seed fit
+ Double_t delta1 = r->Parameter(kDelta);
+ Double_t xi1 = r->Parameter(kXi);
+ Double_t maxEi = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
+ Double_t minE = f->GetXmin();
+
+ Int_t minEb = dist->GetXaxis()->FindBin(minE);
+ Int_t maxEb = dist->GetXaxis()->FindBin(maxEi);
+ Double_t rmsE = dist->GetRMS();
+ Double_t intg = dist->Integral(minEb, maxEb);
+ if (intg <= 0) {
+ ::Warning("FitNParticle",
+ "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
+ dist->GetName(), minE, maxEi, minEb, maxEb, intg);
+ return 0;
+ }
+
+ // Array of weights
+ TArrayD a(n-1);
+ for (UShort_t i = 2; i <= n; i++)
+ a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
+ // Make the fit function
+ f = AliLandauGaus::MakeFn(r->Parameter(kC),
+ r->Parameter(kDelta),
+ r->Parameter(kXi),
+ r->Parameter(kSigma),
+ r->Parameter(kSigmaN),
+ n, a.fArray, minE, maxEi);
+ SetParLimits(f,kDelta, f->GetParameter(kDelta), minE, fMaxRange);
+ SetParLimits(f,kXi, f->GetParameter(kXi), 0, 2*rmsE); //0.1
+ SetParLimits(f,kSigma, f->GetParameter(kSigma), 1e-5, 2*rmsE); //0.1
+ if (sigman <= 0) f->FixParameter(kSigmaN, 0);
+ else
+ SetParLimits(f,kSigmaN, f->GetParameter(kSigmaN), 0, rmsE);
+
+ // Set the range and name of the scale parameters
+ for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
+ SetParLimits(f,kA+i-2, a[i-2], 0, 1);
+ }
+
+ // Do the fit
+ TString opts(Form("%s%s", GetFitOptions(), fDebug ? "" : "Q"));
+ if (fDebug)
+ ::Info(/*"FitNParticle"*/"",
+ "Fitting in the range %f,%f (%d)", minE, maxEi, n);
+ TFitResultPtr tr = dist->Fit(f, opts, "", minE, maxEi);
+
+ // f->SetRange(minE, fMaxRange);
+ fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
+ fFunctions.AddAtAndExpand(f, n-1);
+
+ return f;
+}
+//____________________________________________________________________
+inline TF1*
+AliLandauGausFitter::FitComposite(TH1* dist, Double_t sigman)
+{
+ // Find the fit range
+ Int_t cutBin = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
+ Int_t maxBin = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
+ dist->GetNbinsX());
+ dist->GetXaxis()->SetRange(cutBin, maxBin);
+
+ // Get the bin with maximum
+ Int_t peakBin = dist->GetMaximumBin();
+ Double_t peakE = dist->GetBinLowEdge(peakBin);
+
+ // Get the low edge
+ // dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
+ Int_t minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
+ Double_t minE = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
+ Double_t maxE = dist->GetBinCenter(peakBin+2*fMinusBins);
+
+ // Get the range in bins and the integral of that range
+ Int_t minEb = dist->GetXaxis()->FindBin(minE);
+ Int_t maxEb = dist->GetXaxis()->FindBin(maxE);
+ Double_t intg = dist->Integral(minEb, maxEb);
+ if (intg <= 0) {
+ ::Warning("Fit1Particle",
+ "Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
+ dist->GetName(), minE, maxE, minEb, maxEb, intg);
+ return 0;
+ }
+
+ // Restore the range
+ dist->GetXaxis()->SetRange(1, maxBin);
+
+ // Define the function to fit
+ TF1* seed = AliLandauGaus::MakeF1(1,peakE,peakE/10,peakE/5,sigman,minE,maxE);
+
+ // Set initial guesses, parameter names, and limits
+ seed->SetParLimits(kDelta, minE, fMaxRange);
+ seed->SetParLimits(kXi, 0.00, 0.1);
+ seed->SetParLimits(kSigma, 1e-5, 0.1);
+ if (sigman <= 0) seed->FixParameter(kSigmaN, 0);
+ else seed->SetParLimits(kSigmaN, 0, fMaxRange);
+
+ // Do the fit, getting the result object
+ if (fDebug)
+ ::Info(/*"FitComposite"*/"", "Fitting seed in the range %f,%f", minE, maxE);
+ /* TFitResultPtr r = */ dist->Fit(seed, GetFitOptions(), "", minE, maxE);
+
+ maxE = dist->GetXaxis()->GetXmax();
+ TF1* comp =
+ AliLandauGaus::MakeComposite(0.8 * seed->GetParameter(kC),
+ seed->GetParameter(kDelta),
+ seed->GetParameter(kDelta)/10,
+ seed->GetParameter(kDelta)/5,
+ 1.20 * seed->GetParameter(kC),
+ seed->GetParameter(kXi),
+ minE, maxE);
+ // comp->SetParLimits(kC, minE, fMaxRange); // C
+ comp->SetParLimits(kDelta, minE, fMaxRange); // Delta
+ comp->SetParLimits(kXi, 0.00, fMaxRange); // Xi
+ comp->SetParLimits(kSigma, 1e-5, fMaxRange); // Sigma
+ // comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
+ comp->SetParLimits(kSigma+2, 0.00, fMaxRange); // Xi'
+ comp->SetLineColor(kRed+1);
+ comp->SetLineWidth(3);
+
+ // Do the fit, getting the result object
+ TString opts(Form("%s%s", GetFitOptions(), fDebug ? "" : "Q"));
+ if (fDebug)
+ ::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
+ /* TFitResultPtr r = */ dist->Fit(comp, opts, "", minE, maxE);
+
+#if 0
+ // This is to store the two components with the output
+ TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
+ part1->SetLineColor(kGreen+1);
+ part1->SetLineWidth(4);
+ part1->SetRange(minE, maxE);
+ part1->SetParameters(comp->GetParameter(0), // C
+ comp->GetParameter(1), // Delta
+ comp->GetParameter(2), // Xi
+ comp->GetParameter(3), // sigma
+ 0);
+ part1->Save(minE,maxE,0,0,0,0);
+ dist->GetListOfFunctions()->Add(part1);
+
+ TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
+ part2->SetLineColor(kBlue+1);
+ part2->SetLineWidth(4);
+ part2->SetRange(minE, maxE);
+ part2->SetParameters(comp->GetParameter(4), // C
+ comp->GetParameter(1), // Delta
+ comp->GetParameter(5), // Xi
+ comp->GetParameter(3), // sigma
+ 0);
+ part2->Save(minE,maxE,0,0,0,0);
+ dist->GetListOfFunctions()->Add(part2);
+#endif
+ return comp;
+}
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff