//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Storage for Fopi PID information from TPC/CDC dE/dx
//      Contains a (normalized) RooAbsPdf 
//      Aware of its momentum range
//
//
// Environment:
//      Software developed for GEM-TPC detector in FOPI.
//
// Author List:
//      Felix Boehmer    TUM            (original author)
//
//
//-----------------------------------------------------------


#include "FopiPidProb.h"

#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooArgSet.h"
#include "TF1.h"


ClassImp(FopiPidProb)

FopiPidProb::FopiPidProb() 
:
  fPdf(NULL),
  fPartType(TString("")),
  fDetType(TString("")),
  fLow(0.),
  fHigh(0.)
{;}

FopiPidProb::FopiPidProb(const FopiPidProb& other) 
:
  fPartType(other.fPartType),
  fDetType(other.fDetType),
  fLow(other.fLow),
  fHigh(other.fHigh),
  fVarName(other.fVarName)
{
  RooAbsPdf* theothers = (RooAbsPdf*) other.getProbDens();
  fPdf = (RooAbsPdf*) theothers->clone("meh");
}


FopiPidProb::FopiPidProb(double momLow, double momHigh, 
			 const TString& partType,
			 const TString& detType,
			 const RooAbsPdf* pdf,
			 const TString& varname) 
:
  fPartType(partType),
  fDetType(detType),
  fLow(momLow),
  fHigh(momHigh),
  fVarName(varname)
{
  fPdf = (RooAbsPdf*) pdf->clone("meh");
}

 

FopiPidProb::~FopiPidProb() 
{
  delete fPdf;
  //delete fVar;
}

void
FopiPidProb::setProb(const RooAbsPdf* pdf) {
  if(fPdf!=NULL) 
    delete fPdf;
  fPdf=(RooAbsPdf*) pdf->clone("meh");
}

void
FopiPidProb::getMomRange(double& low, double& high) const {
  low = fLow;
  high = fHigh;
}


TF1*
FopiPidProb::getTF1() const {
  //TODO: somehow find the name of the var intelligently 
  //OR build in exception handling
  //or simply give the name to the constructor
  RooRealVar* e = (RooRealVar*) fPdf->getVariables()->find(fVarName);
  TF1* tf = fPdf->asTF(RooArgList(*e), RooArgList(), RooArgSet());
  tf->SetNpx(1000);
  return tf;
}

double
FopiPidProb::getLikelihoodValue(double min, double max) const {
  RooRealVar* x = (RooRealVar*) fPdf->getVariables()->find(fVarName);
  x->setRange("intRange",min,max);
  RooAbsReal* integral = fPdf->createIntegral(*x,RooArgSet(*x),"intRange");
  double res = integral->getVal();
  delete integral;
  //delete x;
  return res;
}