/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitModels                                                     *
 *    File: $Id: RooChebychev.cc,v 1.11 2005/06/20 15:51:06 wverkerke Exp $
 * Authors:                                                                  *
 *   GR, Gerhard Raven,   UC San Diego, Gerhard.Raven@slac.stanford.edu
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

// -- CLASS DESCRIPTION [PDF] --

#include "RooFit.h"

#include "Riostream.h"
#include "Riostream.h"
#include <math.h>

#include "RooChebychev.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooArgList.h"

ClassImp(RooChebychev)
;

 RooChebychev::RooChebychev()
{
}

 RooChebychev::RooChebychev(const char* name, const char* title, 
                           RooAbsReal& x, const RooArgList& coefList): 
  RooAbsPdf(name, title),
  _x("x", "Dependent", this, x),
  _coefList("coefList","List of coefficients",this)
{
  // Constructor
  TIterator* coefIter = coefList.createIterator() ;
  RooAbsArg* coef ;
  while((coef = (RooAbsArg*)coefIter->Next())) {
    if (!dynamic_cast<RooAbsReal*>(coef)) {
      cout << "RooChebychev::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName() 
	   << " is not of type RooAbsReal" << endl ;
      assert(0) ;
    }
    _coefList.add(*coef) ;
  }
  delete coefIter ;
}


 RooChebychev::RooChebychev(const RooChebychev& other, const char* name) :
  RooAbsPdf(other, name), 
  _x("x", this, other._x), 
  _coefList("coefList",this,other._coefList)
{
}

inline static double p0(double ,double a) {  return a; }
inline static double p1(double t,double a,double b) {  return a*t+b; }
inline static double p2(double t,double a,double b,double c) {  return p1(t,p1(t,a,b),c); }
inline static double p3(double t,double a,double b,double c,double d) {  return p2(t,p1(t,a,b),c,d); }
inline static double p4(double t,double a,double b,double c,double d,double e) {  return p3(t,p1(t,a,b),c,d,e); }

 Double_t RooChebychev::evaluate() const 
{

  Double_t xmin = _x.min(); Double_t xmax = _x.max();
  Double_t x(-1+2*(_x-xmin)/(xmax-xmin));
  Double_t x2(x*x);
  Double_t sum(0) ;
  switch (_coefList.getSize()) {
  case  7: sum+=((RooAbsReal&)_coefList[6]).getVal()*x*p3(x2,64,-112,56,-7);
  case  6: sum+=((RooAbsReal&)_coefList[5]).getVal()*p3(x2,32,-48,18,-1);
  case  5: sum+=((RooAbsReal&)_coefList[4]).getVal()*x*p2(x2,16,-20,5);
  case  4: sum+=((RooAbsReal&)_coefList[3]).getVal()*p2(x2,8,-8,1);
  case  3: sum+=((RooAbsReal&)_coefList[2]).getVal()*x*p1(x2,4,-3);
  case  2: sum+=((RooAbsReal&)_coefList[1]).getVal()*p1(x2,2,-1);
  case  1: sum+=((RooAbsReal&)_coefList[0]).getVal()*x;
  case  0: sum+=1;
  }
  return sum;
}

 Int_t RooChebychev::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const 
{
  if (matchArgs(allVars, analVars, _x)) return 1;
  return 0;
}

 Double_t RooChebychev::analyticalIntegral(Int_t code, const char* rangeName) const 
{
  assert(code==1) ;
  // WVE check that this works all right for sub ranges
  Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName);
  Double_t norm(0) ;
  switch(_coefList.getSize()) {
  case  7: case  6: norm+=((RooAbsReal&)_coefList[5]).getVal()*(-1 + 18./3. - 48./5. + 32./7.);
  case  5: case  4: norm+=((RooAbsReal&)_coefList[3]).getVal()*( 1 -  8./3. +  8./5.);
  case  3: case  2: norm+=((RooAbsReal&)_coefList[1]).getVal()*(-1 +  2./3.);
  case  1: case  0: norm+= 1;
  }
  norm *= xmax-xmin;
  return norm;
}