AliPhysics  ebc57ae (ebc57ae)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliLandauGausFitter.h
Go to the documentation of this file.
1 #ifndef ALILANDAUGAUSFITTER_H
2 #define ALILANDAUGAUSFITTER_H
3 
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 
31 {
32 public:
36  enum {
52  };
58  static const char* GetFitOptions() { return "RNS"; }
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  }
78  {
79  fFitResults.Delete();
80  fFunctions.Delete();
81  }
87  void SetDebug(Bool_t debug=true) { fDebug = debug; }
92  void Clear()
93  {
94  fFitResults.Clear();
95  fFunctions.Clear();
96  }
109  TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
123  TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
135  TF1* FitComposite(TH1* dist, Double_t sigman);
141  Double_t GetLowCut() const { return fLowCut; }
147  Double_t GetMaxRange() const { return fMaxRange; }
153  UShort_t GetMinusBins() const { return fMinusBins; }
159  const TObjArray& GetFitResults() const { return fFitResults; }
171  const TObjArray& GetFunctions() const { return fFunctions; }
178 protected:
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*
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, 2*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*
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*
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:
double Double_t
Definition: External.C:58
const int debug
Definition: scanAll.C:15
void SetParLimits(TF1 *f, Int_t iPar, Double_t test, Double_t low, Double_t high)
static TF1 * MakeFn(Double_t c, Double_t delta, Double_t xi, Double_t sigma, Double_t sigma_n, Int_t n, const Double_t *a, Double_t xmin, Double_t xmax)
Double_t GetMaxRange() const
static TF1 * MakeF1(Double_t c, Double_t delta, Double_t xi, Double_t sigma, Double_t sigma_n, Double_t xmin, Double_t xmax)
const TObjArray & GetFitResults() const
Double_t GetLowCut() const
Declaration and implementation of Landau-Gauss distributions.
static TF1 * MakeComposite(Double_t c1, Double_t delta, Double_t xi1, Double_t sigma, Double_t c2, Double_t xi2, Double_t xmin, Double_t xmax)
const TObjArray & GetFunctions() const
int Int_t
Definition: External.C:63
TF1 * FitComposite(TH1 *dist, Double_t sigman)
static const char * GetFitOptions()
TF1 * Fit1Particle(TH1 *dist, Double_t sigman=-1)
AliLandauGausFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins)
unsigned short UShort_t
Definition: External.C:28
TF1 * FitNParticle(TH1 *dist, UShort_t n, Double_t sigman=-1)
void test(int runnumber=195345)
bool Bool_t
Definition: External.C:53
UShort_t GetMinusBins() const
void SetDebug(Bool_t debug=true)
Definition: External.C:196