AnalyticalIntegrals.cxx

Go to the documentation of this file.

// @(#)root/hist:$Id$
// Authors: L. Moneta, A. Flandi 2015
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2015  ROOT  Team, CERN/PH-SFT                        *
 *                                                                    *
 *                                                                    *
 **********************************************************************/
//
//  AnalyticalIntegrals.cxx
//
//
//  Created by Aurélie Flandi on 09.09.14.
//
//
 
#include "AnalyticalIntegrals.h"
 
#include "TROOT.h"
#include "TF1.h"
#include "TFormula.h"
#include "TMath.h"
#include "Math/DistFuncMathCore.h" //for cdf
 
#include <cstdio>
 
Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b)
{
 
   Double_t xmin = a;
   Double_t xmax = b;
   Int_t    num  = f->GetNumber();
   Double_t *p   = f->GetParameters();
   Double_t result = 0.;
 
   TFormula * formula = f->GetFormula();
   if (!formula) {
      Error("TF1::AnalyticalIntegral","Invalid formula number - return a NaN");
      return TMath::QuietNaN();
   }
 
   else if (num == 200) { // expo: exp(p0 + p1*x)
      const double p0 = p[0];
      const double p1 = p[1];
 
      if (p1 == 0) {
         // Limit p1 -> 0: integral of constant exp(p0)
         result = std::exp(p0) * (xmax - xmin);
      } else {
         const double ea = p0 + p1 * xmin;
         const double eb = p0 + p1 * xmax;
 
         // Use expm1 for improved numerical stability when eb ~ ea
         result = std::exp(ea) * std::expm1(eb - ea) / p1;
      }
   }
 
   else if (num == 100) // gaus: [0]*exp(-0.5*((x-[1])/[2])^2))
   {
      double amp   = p[0];
      double mean  = p[1];
      double sigma = p[2];
      if (sigma <= 0)
         return TMath::QuietNaN();
      if (formula->TestBit(TFormula::kNormalized))
         result = amp * (ROOT::Math::gaussian_cdf(xmax, sigma, mean) - ROOT::Math::gaussian_cdf(xmin, sigma, mean));
      else
         result = amp * sqrt(2 * TMath::Pi()) * sigma *
                  (ROOT::Math::gaussian_cdf(xmax, sigma, mean) - ROOT::Math::gaussian_cdf(xmin, sigma, mean)); //
   } else if (num == 400) // landau: root::math::landau(x,mpv=0,sigma=1,bool norm=false)
   {
 
      double amp   = p[0];
      double mean  = p[1];
      double sigma = p[2];
      if (sigma <= 0)
         return TMath::QuietNaN();
      //printf("computing integral for landau in [%f,%f] for m=%f s = %f \n",xmin,xmax,mean,sigma);
      if (formula->TestBit(TFormula::kNormalized) )
         result = amp*(ROOT::Math::landau_cdf(xmax,sigma,mean) - ROOT::Math::landau_cdf(xmin,sigma,mean));
      else
         result = amp*sigma*(ROOT::Math::landau_cdf(xmax,sigma,mean) - ROOT::Math::landau_cdf(xmin,sigma,mean));
   } else if (num == 500) // crystal ball
   {
      double amp   = p[0];
      double mean  = p[1];
      double sigma = p[2];
      double alpha = p[3];
      double n     = p[4];
      if (sigma <= 0 || n <= 0)
         return TMath::QuietNaN();
      //printf("computing integral for CB in [%f,%f] for m=%f s = %f alpha = %f n = %f\n",xmin,xmax,mean,sigma,alpha,n);
      if (alpha > 0)
         result = amp*( ROOT::Math::crystalball_integral(xmin,alpha,n,sigma,mean) -  ROOT::Math::crystalball_integral(xmax,alpha,n,sigma,mean) );
      else {
         result = amp*( ROOT::Math::crystalball_integral(xmax,alpha,n,sigma,mean) -  ROOT::Math::crystalball_integral(xmin,alpha,n,sigma,mean) );
      }
   }
 
   else if (num >= 300 && num < 400) // polN
   {
      Int_t n = num - 300;
      for (int i=0;i<n+1;i++)
      {
         result += p[i]/(i+1)*(std::pow(xmax,i+1)-std::pow(xmin,i+1));
      }
   } else
      result = TMath::QuietNaN();
 
   return result;
}