#include <VerticalInterpHistPdf.h>

Public Member Functions
virtual TObject *	clone (const char *newname) const
const RooArgList &	coefList () const
Double_t	evaluate () const
const RooArgList &	funcList () const
Bool_t	selfNormalized () const
	VerticalInterpHistPdf (const VerticalInterpHistPdf &other, const char *name=0)
	VerticalInterpHistPdf (const char name, const char title, const RooRealVar &x, const RooArgList &funcList, const RooArgList &coefList, Double_t smoothRegion=1., Int_t smoothAlgo=1)
	VerticalInterpHistPdf ()
const RooRealProxy &	x () const
virtual	~VerticalInterpHistPdf ()
Protected Attributes
TH1 **	_cacheSingle
	not to be serialized
int *	_cacheSingleGood
	not to be serialized
TH1 *	_cacheTotal
TIterator *	_coefIter
	Iterator over FUNC list.
RooListProxy	_coefList
TIterator *	_funcIter
RooListProxy	_funcList
SimpleCacheSentry	_sentry
	Iterator over coefficient list.
Int_t	_smoothAlgo
Double_t	_smoothRegion
RooRealProxy	_x
Private Member Functions
double	smoothStepFunc (double x) const
void	syncComponent (int which) const
ClassDef(VerticalInterpHistPdf, 1) protected void	syncTotal () const
	not to be serialized
Friends
class	FastVerticalInterpHistPdf

Detailed Description

Definition at line 20 of file VerticalInterpHistPdf.h.

Constructor & Destructor Documentation

VerticalInterpHistPdf::VerticalInterpHistPdf ( )

Referenced by clone().

VerticalInterpHistPdf::VerticalInterpHistPdf	(	const char *	name,
		const char *	title,
		const RooRealVar &	x,
		const RooArgList &	funcList,
		const RooArgList &	coefList,
		Double_t	smoothRegion = `1.`,
		Int_t	smoothAlgo = `1`
	)

Definition at line 37 of file VerticalInterpHistPdf.cc.

References _coefIter, _coefList, _funcIter, and _funcList.

                                                                                                                                                                                                          :
  RooAbsPdf(name,title),
  _x("x","Independent variable",this,const_cast<RooRealVar&>(x)),
  _funcList("funcList","List of pdfs",this),
  _coefList("coefList","List of coefficients",this), // we should get shapeDirty when coefficients change
  _smoothRegion(smoothRegion),
  _smoothAlgo(smoothAlgo),
  _cacheTotal(0),
  _cacheSingle(0),
  _cacheSingleGood(0) 
{ 

  if (inFuncList.getSize()!=2*inCoefList.getSize()+1) {
    coutE(InputArguments) << "VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() 
                          << ") number of pdfs and coefficients inconsistent, must have Nfunc=1+2*Ncoef" << endl ;
    assert(0);
  }

  TIterator* funcIter = inFuncList.createIterator() ;
  RooAbsArg* func;
  while((func = (RooAbsArg*)funcIter->Next())) {
    RooAbsPdf *pdf = dynamic_cast<RooAbsPdf*>(func);
    if (!pdf) {
      coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") function  " << func->GetName() << " is not of type RooAbsPdf" << endl;
      assert(0);
    }
    RooArgSet *params = pdf->getParameters(RooArgSet(x));
    if (params->getSize() > 0) {
      coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") pdf  " << func->GetName() << " has some parameters." << endl;
      assert(0);
    }
    delete params;
    _funcList.add(*func) ;
  }
  delete funcIter;

  TIterator* coefIter = inCoefList.createIterator() ;
  RooAbsArg* coef;
  while((coef = (RooAbsArg*)coefIter->Next())) {
    if (!dynamic_cast<RooAbsReal*>(coef)) {
      coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal" << endl;
      assert(0);
    }
    _coefList.add(*coef) ;    
  }
  delete coefIter;

  _funcIter  = _funcList.createIterator() ;
  _coefIter = _coefList.createIterator() ;

}

VerticalInterpHistPdf::VerticalInterpHistPdf	(	const VerticalInterpHistPdf &	other,
		const char *	name = `0`
	)

Definition at line 93 of file VerticalInterpHistPdf.cc.

References _coefIter, _coefList, _funcIter, and _funcList.

                                                                                                 :
  RooAbsPdf(other,name),
  _x("x",this,other._x),
  _funcList("funcList",this,other._funcList),
  _coefList("coefList",this,other._coefList),
  _smoothRegion(other._smoothRegion),
  _smoothAlgo(other._smoothAlgo),
  _cacheTotal(0),
  _cacheSingle(0),
  _cacheSingleGood(0) 
{
  // Copy constructor

  _funcIter  = _funcList.createIterator() ;
  _coefIter = _coefList.createIterator() ;
}

VerticalInterpHistPdf::~VerticalInterpHistPdf ( ) [virtual]

Definition at line 113 of file VerticalInterpHistPdf.cc.

References _cacheSingle, _cacheSingleGood, _cacheTotal, _coefIter, _funcIter, _funcList, and i.

{
  // Destructor
  delete _funcIter ;
  delete _coefIter ;
  if (_cacheTotal) {
      delete _cacheTotal;
      for (int i = 0; i < _funcList.getSize(); ++i) delete _cacheSingle[i]; 
      delete [] _cacheSingle;
      delete [] _cacheSingleGood;
  }
}

Member Function Documentation

virtual TObject* VerticalInterpHistPdf::clone ( const char * newname ) const [inline, virtual]

Definition at line 26 of file VerticalInterpHistPdf.h.

References VerticalInterpHistPdf().

{ return new VerticalInterpHistPdf(*this,newname) ; }

const RooArgList& VerticalInterpHistPdf::coefList ( ) const [inline]

Definition at line 34 of file VerticalInterpHistPdf.h.

References _coefList.

{ return _coefList ; }

Double_t VerticalInterpHistPdf::evaluate ( ) const

Definition at line 127 of file VerticalInterpHistPdf.cc.

References _cacheTotal, _coefList, _sentry, _x, gather_cfg::cout, SimpleCacheSentry::good(), i, and syncTotal().

{
  if (_cacheTotal == 0) setupCaches();
#if 0
  printf("Evaluate called for x %f, cache status %d: ", _x.arg().getVal(), _sentry.good());
  int ndim = _coefList.getSize();
  for (int i = 0; i < ndim; ++i) { 
    RooAbsReal *rar = dynamic_cast<RooAbsReal *>(_coefList.at(i));
    printf("%s = %+6.4f  ", rar->GetName(), rar->getVal());
  }
  std::cout << std::endl;
#endif

  if (!_sentry.good()) syncTotal();
  int nbin = _cacheTotal->GetNbinsX(), ibin = _cacheTotal->FindBin(_x);
  if (ibin < 1) ibin = 1;
  else if (ibin > nbin) ibin = nbin;
  return _cacheTotal->GetBinContent(ibin);
}

const RooArgList& VerticalInterpHistPdf::funcList ( ) const [inline]

Definition at line 33 of file VerticalInterpHistPdf.h.

References _funcList.

{ return _funcList ; }

Bool_t VerticalInterpHistPdf::selfNormalized ( ) const [inline]

Definition at line 29 of file VerticalInterpHistPdf.h.

{ return kTRUE; }

double VerticalInterpHistPdf::smoothStepFunc ( double x ) const [inline, private]

Definition at line 63 of file VerticalInterpHistPdf.h.

References _smoothRegion.

Referenced by syncTotal().

                                               { 
    if (fabs(x) >= _smoothRegion) return x > 0 ? +1 : -1;
    double xnorm = x/_smoothRegion, xnorm2 = xnorm*xnorm;
    return 0.125 * xnorm * (xnorm2 * (3.*xnorm2 - 10.) + 15);
  }

void VerticalInterpHistPdf::syncComponent ( int which ) const [private]

Definition at line 148 of file VerticalInterpHistPdf.cc.

References _cacheSingle, _cacheSingleGood, _funcList, _smoothAlgo, _x, b, i, create_public_lumi_plots::log, and detailsBasic3DVector::y.

Referenced by syncTotal().

                                                     {
    RooAbsPdf *pdfi = dynamic_cast<RooAbsPdf *>(_funcList.at(i));
    if (_cacheSingle[i] != 0) delete _cacheSingle[i];
    _cacheSingle[i] = pdfi->createHistogram("",dynamic_cast<const RooRealVar &>(_x.arg())); 
    _cacheSingle[i]->SetDirectory(0);
    if (_cacheSingle[i]->Integral("width")) { _cacheSingle[i]->Scale(1.0/_cacheSingle[i]->Integral("width")); }
    if (i > 0) {
        for (int b = 1, nb = _cacheSingle[i]->GetNbinsX(); b <= nb; ++b) {
            double y  = _cacheSingle[i]->GetBinContent(b);
            double y0 = _cacheSingle[0]->GetBinContent(b);
            if (_smoothAlgo < 0) {
                if (y > 0 && y0 > 0) {
                    double logk = log(y/y0);
                    // odd numbers correspond to up variations, even numbers to down variations,
                    // and down variations need -log(kappa) instead of log(kappa)
                    _cacheSingle[i]->SetBinContent(b, logk);
                } else {
                    _cacheSingle[i]->SetBinContent(b, 0);
                }
            } else {
                _cacheSingle[i]->SetBinContent(b, y - y0);
            }
        }
    }
    _cacheSingleGood[i] = true;
}

void VerticalInterpHistPdf::syncTotal ( ) const [private]

not to be serialized

Definition at line 175 of file VerticalInterpHistPdf.cc.

References _cacheSingle, _cacheSingleGood, _cacheTotal, _coefIter, _coefList, _sentry, _smoothAlgo, alpha, b, alignCSCRings::e, create_public_lumi_plots::exp, i, SimpleCacheSentry::reset(), smoothStepFunc(), syncComponent(), and x().

Referenced by evaluate().

                                            {
    int ndim = _coefList.getSize();
    for (int i = 0; i < 2*ndim+1; ++i) {
        if (!_cacheSingleGood[i]) syncComponent(i);
    }
    for (int b = 1, nb = _cacheTotal->GetNbinsX(); b <= nb; ++b) {
        double val = _cacheSingle[0]->GetBinContent(b);
        _coefIter->Reset();
        for (int i = 0; i < ndim; ++i) {
            double dhi = _cacheSingle[2*i+1]->GetBinContent(b);
            double dlo = _cacheSingle[2*i+2]->GetBinContent(b);
            double x = (dynamic_cast<RooAbsReal *>(_coefIter->Next()))->getVal();
            double alpha = x * 0.5 * ((dhi-dlo) + (dhi+dlo)*smoothStepFunc(x));
            // alpha(0) = 0
            // alpha(+1) = dhi 
            // alpha(-1) = dlo
            // alpha(x >= +1) = |x|*dhi
            // alpha(x <= -1) = |x|*dlo
            // alpha is continuous and has continuous first and second derivative, as smoothStepFunc has them
            if (_smoothAlgo < 0) {
                val *= exp(alpha);
            } else {
                val += alpha; 
            }
        }    
        if (val <= 0) val = 1e-9;
        _cacheTotal->SetBinContent(b, val);
    }
    double norm = _cacheTotal->Integral("width");
    if (norm > 0) _cacheTotal->Scale(1.0/norm);
    _sentry.reset();
}