[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliLandauGausFitter.h

#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:
Commit	Line	Data
59eb7e1a	1	#ifndef ALILANDAUGAUSFITTER_H
	2	#define ALILANDAUGAUSFITTER_H
	3	/**
	4	* @file AliLandauGausFitter.h
	5	* @author Christian Holm Christensen <cholm@nbi.dk>
	6	* @date Tue Mar 11 08:56:03 2014
	7	*
	8	* @brief Declaration and implementation of fitter of Landau-Gauss
	9	* distributions to energy loss spectra.
	10	*
	11	* @ingroup pwglf_forward
	12	*/
	13	#include <AliLandauGaus.h>
	14	#include <TH1.h>
	15	#include <TF1.h>
	16	#include <TMath.h>
	17	#include <TObjArray.h>
	18	#include <TString.h>
	19	#include <TArray.h>
	20	#include <TFitResult.h>
	21	#include <TError.h>
	22
	23	/**
	24	* Fit Landau-Gauss distributions to energy loss distributions
	25	*
	26	* @see AliLandauGaus
	27	*
	28	* @ingroup pwglf_forward
	29	*/
	30	class AliLandauGausFitter
	31	{
	32	public:
	33	/**
	34	* Enumeration of parameters
	35	*/
	36	enum {
	37	/** Index of pre-constant @f$ C@f$ */
	38	kC = AliLandauGaus::kC,
	39	/** Index of most probable value @f$ \Delta_p@f$ */
	40	kDelta = AliLandauGaus::kDelta,
	41	/** Index of Landau width @f$ \xi@f$ */
	42	kXi = AliLandauGaus::kXi,
	43	/** Index of Gaussian width @f$ \sigma@f$ */
	44	kSigma = AliLandauGaus::kSigma,
	45	/** Index of Gaussian additional width @f$ \sigma_n@f$ */
	46	kSigmaN = AliLandauGaus::kSigmaN,
	47	/** Index of Number of particles @f$ N@f$ */
	48	kN = AliLandauGaus::kN,
	49	/** Base index of particle strengths @f$ a_i@f$ for
	50	@f$i=2,\ldots,N@f$ */
	51	kA = AliLandauGaus::kA
	52	};
	53	/**
	54	* Get the fit options to use
	55	*
	56	* @return Fit options used through-out
	57	*/
	58	static const char* GetFitOptions() { return "RNS"; }
	59	/**
	60	* Constructor
	61	*
	62	* @param lowCut Lower cut of spectrum - data below this cuts is ignored
	63	* @param maxRange Maximum range to fit to
	64	* @param minusBins The number of bins below maximum to use
65	*/
66	AliLandauGausFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins)
67	: fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins),
68	fFitResults(0), fFunctions(0), fDebug(false)
69	{
70	fFitResults.SetOwner();
71	fFunctions.SetOwner();
72	}
73	/**
74	* Destructor
75	*
76	*/
77	virtual ~AliLandauGausFitter()
78	{
79	fFitResults.Delete();
80	fFunctions.Delete();
81	}
82	/**
83	* Enable/disable debugging output.
84	*
85	* @param debug If true, enable debugging output
86	*/
87	void SetDebug(Bool_t debug=true) { fDebug = debug; }
88	/**
89	* Clear internal arrays
90	*
91	*/
92	void Clear()
93	{
94	fFitResults.Clear();
95	fFunctions.Clear();
96	}
97	/**
98	* Fit a 1-particle signal to the passed energy loss distribution
99	*
100	* Note that this function clears the internal arrays first
101	*
102	* @param dist Data to fit the function to
103	* @param sigman If larger than zero, the initial guess of the
104	* detector induced noise. If zero or less, then this
105	* parameter is ignored in the fit (fixed at 0)
106	*
107	* @return The function fitted to the data
108	*/
109	TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
110	/**
111	* Fit a N-particle signal to the passed energy loss distribution
112	*
113	* If there's no 1-particle fit present, it does that first
114	*
115	* @param dist Data to fit the function to
116	* @param n Number of particle signals to fit
117	* @param sigman If larger than zero, the initial guess of the
118	* detector induced noise. If zero or less, then this
119	* parameter is ignored in the fit (fixed at 0)
120	*
121	* @return The function fitted to the data
122	*/
123	TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
124	/**
125	* Fit a composite distribution of energy loss from both primaries
126	* and secondaries
127	*
128	* @param dist Distribution
129	* @param sigman If larger than zero, the initial guess of the
130	* detector included noise. If zero or less this
131	* parameter is fixed to 0.
132	*
133	* @return Function fitted to the data
134	*/
135	TF1* FitComposite(TH1* dist, Double_t sigman);
136	/**
137	* Get Lower cut on data
138	*
139	* @return Lower cut on data
140	*/
141	Double_t GetLowCut() const { return fLowCut; }
142	/**
143	* Get Maximum range to fit
144	*
145	* @return Maximum range to fit
146	*/
147	Double_t GetMaxRange() const { return fMaxRange; }
148	/**
149	* Get Number of bins from maximum to fit 1st peak
150	*
151	* @return Number of bins from maximum to fit 1st peak
152	*/
153	UShort_t GetMinusBins() const { return fMinusBins; }
154	/**
155	* Get Array of fit results
156	*
157	* @return Array of fit results
158	*/
159	const TObjArray& GetFitResults() const { return fFitResults; }
160	/**
161	* Get Array of fit results
162	*
163	* @return Array of fit results
164	*/
165	TObjArray& GetFitResults() { return fFitResults; }
166	/**
167	* Get Array of functions
168	*
169	* @return Array of functions
170	*/
171	const TObjArray& GetFunctions() const { return fFunctions; }
172	/**
173	* Get Array of functions
174	*
175	* @return Array of functions
176	*/
177	TObjArray& GetFunctions() { return fFunctions; }
178	protected:
179	/**
180	* Set parameter limits on the function @a f. The limits are only
181	* set if @f$ low \le test \le high@f$
182	*
183	* @param f Function object
184	* @param iPar Parameter number
185	* @param test Initial guess
186	* @param low Low limit
187	* @param high high limit
188	*/
189	void SetParLimits(TF1* f, Int_t iPar, Double_t test,
190	Double_t low, Double_t high)
191	{
192	if (test < low \|\| test > high) {
193	::Warning("","Initial value of %s=%f not in [%f,%f]",
194	f->GetParName(iPar), test, low, high);
195	return;
196	}
197	if (fDebug)
198	::Info(/"SetParLimits"/"", "Set par limits on %-12s=%f: [%f,%f]",
199	f->GetParName(iPar), test, low, high);
200	f->SetParLimits(iPar, low, high);
201	}
202	const Double_t fLowCut; // Lower cut on data
203	const Double_t fMaxRange; // Maximum range to fit
204	const UShort_t fMinusBins; // Number of bins from maximum to fit 1st peak
205	TObjArray fFitResults; // Array of fit results
206	TObjArray fFunctions; // Array of functions
207	Bool_t fDebug; // Debug flag
208	};
209
210
211	//____________________________________________________________________
212	inline TF1*
213	AliLandauGausFitter::Fit1Particle(TH1* dist, Double_t sigman)
214	{
215	// Clear the cache
216	Clear();
217
218	// Find the fit range
219	Int_t cutBin = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
220	Int_t maxBin = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
221	dist->GetNbinsX());
222	dist->GetXaxis()->SetRange(cutBin, maxBin);
223	// dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
224
225	// Get the bin with maximum
226	Int_t peakBin = dist->GetMaximumBin();
227	Double_t peakE = dist->GetBinLowEdge(peakBin);
228	Double_t rmsE = dist->GetRMS();
229
230	// Get the low edge
231	// dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
232	Int_t minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
233	Double_t minE = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
234	Double_t maxE = dist->GetBinCenter(peakBin+2*fMinusBins);
235
236	Int_t minEb = dist->GetXaxis()->FindBin(minE);
237	Int_t maxEb = dist->GetXaxis()->FindBin(maxE);
238	Double_t intg = dist->Integral(minEb, maxEb);
239	if (intg <= 0) {
240	::Warning("Fit1Particle",
241	"Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
242	dist->GetName(), minE, maxE, minEb, maxEb, intg);
243	return 0;
244	}
245
246	// Restore the range
247	dist->GetXaxis()->SetRange(1, maxBin);
248
249	// Define the function to fit
250	TF1* f = AliLandauGaus::MakeF1(intg,peakE,peakE/10,peakE/5,sigman,minE,maxE);
251	SetParLimits(f, kDelta, peakE, minE, fMaxRange);
252	SetParLimits(f, kXi, peakE, 0, rmsE); // 0.1
253	SetParLimits(f, kSigma, peakE/5, 1e-5, rmsE); // 0.1
254	if (sigman <= 0)
255	f->FixParameter(kSigmaN, 0);
256	else
257	SetParLimits(f, kSigmaN, peakE, 0, rmsE);
258
259
260	TString opts(Form("%s%s", GetFitOptions(), fDebug ? "" : "Q"));
261	// Do the fit, getting the result object
262	if (fDebug)
263	::Info(/"Fit1Particle"/"", "Fitting in the range %f,%f", minE, maxE);
264	TFitResultPtr r = dist->Fit(f, opts, "", minE, maxE);
265	if (!r.Get()) {
266	::Warning("Fit1Particle",
267	"No fit returned when processing %s in the range [%f,%f] "
268	"options %s", dist->GetName(), minE, maxE, opts.Data());
269	return 0;
270	}
271	// f->SetRange(minE, fMaxRange);
272	fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
273	fFunctions.AddAtAndExpand(f, 0);
274
275	return f;
276	}
277	//____________________________________________________________________
278	inline TF1*
279	AliLandauGausFitter::FitNParticle(TH1* dist, UShort_t n, Double_t sigman)
280	{
281	// Get the seed fit result
282	TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
283	TF1* f = static_cast<TF1*>(fFunctions.At(0));
284	if (!r \|\| !f) {
285	f = Fit1Particle(dist, sigman);
286	r = static_cast<TFitResult*>(fFitResults.At(0));
287	if (!r \|\| !f) {
288	::Warning("FitNLandau", "No first shot at landau fit");
289	return 0;
290	}
291	}
292
293	// Get some parameters from seed fit
294	Double_t delta1 = r->Parameter(kDelta);
295	Double_t xi1 = r->Parameter(kXi);
296	Double_t maxEi = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
297	Double_t minE = f->GetXmin();
298
299	Int_t minEb = dist->GetXaxis()->FindBin(minE);
300	Int_t maxEb = dist->GetXaxis()->FindBin(maxEi);
301	Double_t rmsE = dist->GetRMS();
302	Double_t intg = dist->Integral(minEb, maxEb);
303	if (intg <= 0) {
304	::Warning("FitNParticle",
305	"Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
306	dist->GetName(), minE, maxEi, minEb, maxEb, intg);
307	return 0;
308	}
309
310	// Array of weights
311	TArrayD a(n-1);
312	for (UShort_t i = 2; i <= n; i++)
313	a.fArray[i-2] = (n == 2 ? 0.05 : 0.000001);
314	// Make the fit function
315	f = AliLandauGaus::MakeFn(r->Parameter(kC),
316	r->Parameter(kDelta),
317	r->Parameter(kXi),
318	r->Parameter(kSigma),
319	r->Parameter(kSigmaN),
320	n, a.fArray, minE, maxEi);
321	SetParLimits(f,kDelta, f->GetParameter(kDelta), minE, fMaxRange);
322	SetParLimits(f,kXi, f->GetParameter(kXi), 0, 2*rmsE); //0.1
323	SetParLimits(f,kSigma, f->GetParameter(kSigma), 1e-5, 2*rmsE); //0.1
324	if (sigman <= 0) f->FixParameter(kSigmaN, 0);
325	else
326	SetParLimits(f,kSigmaN, f->GetParameter(kSigmaN), 0, rmsE);
327
328	// Set the range and name of the scale parameters
329	for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit
330	SetParLimits(f,kA+i-2, a[i-2], 0, 1);
331	}
332
333	// Do the fit
334	TString opts(Form("%s%s", GetFitOptions(), fDebug ? "" : "Q"));
335	if (fDebug)
336	::Info(/"FitNParticle"/"",
337	"Fitting in the range %f,%f (%d)", minE, maxEi, n);
338	TFitResultPtr tr = dist->Fit(f, opts, "", minE, maxEi);
339
340	// f->SetRange(minE, fMaxRange);
341	fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
342	fFunctions.AddAtAndExpand(f, n-1);
343
344	return f;
345	}
346	//____________________________________________________________________
347	inline TF1*
348	AliLandauGausFitter::FitComposite(TH1* dist, Double_t sigman)
349	{
350	// Find the fit range
351	Int_t cutBin = TMath::Max(dist->GetXaxis()->FindBin(fLowCut),3);
352	Int_t maxBin = TMath::Min(dist->GetXaxis()->FindBin(fMaxRange),
353	dist->GetNbinsX());
354	dist->GetXaxis()->SetRange(cutBin, maxBin);
355
356	// Get the bin with maximum
357	Int_t peakBin = dist->GetMaximumBin();
358	Double_t peakE = dist->GetBinLowEdge(peakBin);
359
360	// Get the low edge
361	// dist->GetXaxis()->SetRangeUser(fLowCut, peakE);
362	Int_t minBin = peakBin - fMinusBins; // dist->GetMinimumBin();
363	Double_t minE = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
364	Double_t maxE = dist->GetBinCenter(peakBin+2*fMinusBins);
365
366	// Get the range in bins and the integral of that range
367	Int_t minEb = dist->GetXaxis()->FindBin(minE);
368	Int_t maxEb = dist->GetXaxis()->FindBin(maxE);
369	Double_t intg = dist->Integral(minEb, maxEb);
370	if (intg <= 0) {
371	::Warning("Fit1Particle",
372	"Integral of %s between [%f,%f] [%03d,%03d] = %f < 0",
373	dist->GetName(), minE, maxE, minEb, maxEb, intg);
374	return 0;
375	}
376
377	// Restore the range
378	dist->GetXaxis()->SetRange(1, maxBin);
379
380	// Define the function to fit
381	TF1* seed = AliLandauGaus::MakeF1(1,peakE,peakE/10,peakE/5,sigman,minE,maxE);
382
383	// Set initial guesses, parameter names, and limits
384	seed->SetParLimits(kDelta, minE, fMaxRange);
385	seed->SetParLimits(kXi, 0.00, 0.1);
386	seed->SetParLimits(kSigma, 1e-5, 0.1);
387	if (sigman <= 0) seed->FixParameter(kSigmaN, 0);
388	else seed->SetParLimits(kSigmaN, 0, fMaxRange);
389
390	// Do the fit, getting the result object
391	if (fDebug)
392	::Info(/"FitComposite"/"", "Fitting seed in the range %f,%f", minE, maxE);
393	/* TFitResultPtr r = */ dist->Fit(seed, GetFitOptions(), "", minE, maxE);
394
395	maxE = dist->GetXaxis()->GetXmax();
396	TF1* comp =
397	AliLandauGaus::MakeComposite(0.8 * seed->GetParameter(kC),
398	seed->GetParameter(kDelta),
399	seed->GetParameter(kDelta)/10,
400	seed->GetParameter(kDelta)/5,
401	1.20 * seed->GetParameter(kC),
402	seed->GetParameter(kXi),
403	minE, maxE);
404	// comp->SetParLimits(kC, minE, fMaxRange); // C
405	comp->SetParLimits(kDelta, minE, fMaxRange); // Delta
406	comp->SetParLimits(kXi, 0.00, fMaxRange); // Xi
407	comp->SetParLimits(kSigma, 1e-5, fMaxRange); // Sigma
408	// comp->SetParLimits(kSigma+1, minE, fMaxRange); // C
409	comp->SetParLimits(kSigma+2, 0.00, fMaxRange); // Xi'
410	comp->SetLineColor(kRed+1);
411	comp->SetLineWidth(3);
412
413	// Do the fit, getting the result object
414	TString opts(Form("%s%s", GetFitOptions(), fDebug ? "" : "Q"));
415	if (fDebug)
416	::Info("FitComposite", "Fitting composite in the range %f,%f", minE, maxE);
417	/* TFitResultPtr r = */ dist->Fit(comp, opts, "", minE, maxE);
418
419	#if 0
420	// This is to store the two components with the output
421	TF1* part1 = static_cast<TF1*>(seed->Clone("part1"));
422	part1->SetLineColor(kGreen+1);
423	part1->SetLineWidth(4);
424	part1->SetRange(minE, maxE);
425	part1->SetParameters(comp->GetParameter(0), // C
426	comp->GetParameter(1), // Delta
427	comp->GetParameter(2), // Xi
428	comp->GetParameter(3), // sigma
429	0);
430	part1->Save(minE,maxE,0,0,0,0);
431	dist->GetListOfFunctions()->Add(part1);
432
433	TF1* part2 = static_cast<TF1*>(seed->Clone("part2"));
434	part2->SetLineColor(kBlue+1);
435	part2->SetLineWidth(4);
436	part2->SetRange(minE, maxE);
437	part2->SetParameters(comp->GetParameter(4), // C
438	comp->GetParameter(1), // Delta
439	comp->GetParameter(5), // Xi
440	comp->GetParameter(3), // sigma
441	0);
442	part2->Save(minE,maxE,0,0,0,0);
443	dist->GetListOfFunctions()->Add(part2);
444	#endif
445	return comp;
446	}
447
448	#endif
449	// Local Variables:
450	// mode: C++
451	// End: