PWGLF/FORWARD/analysis2/AliLandauGausFitter.h

   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: