// @(#)root/tmva $Id: TMVA_MethodBase.cxx,v 1.2 2006/05/09 08:37:06 brun Exp $
// Author: Andreas Hoecker, Helge Voss, Kai Voss

/**********************************************************************************
 * Project: TMVA - a Root-integrated toolkit for multivariate data analysis       *
 * Package: TMVA                                                                  *
 * Class  : TMVA_MethodBase                                                       *
 *                                                                                *
 * Description:                                                                   *
 *      Implementation (see header for description)                               *
 *                                                                                *
 * Authors (alphabetical):                                                        *
 *      Andreas Hoecker <Andreas.Hocker@cern.ch> - CERN, Switzerland              *
 *      Xavier Prudent  <prudent@lapp.in2p3.fr>  - LAPP, France                   *
 *      Helge Voss      <Helge.Voss@cern.ch>     - MPI-KP Heidelberg, Germany     *
 *      Kai Voss        <Kai.Voss@cern.ch>       - U. of Victoria, Canada         *
 *                                                                                *
 * Copyright (c) 2005:                                                            *
 *      CERN, Switzerland,                                                        *
 *      U. of Victoria, Canada,                                                   *
 *      MPI-KP Heidelberg, Germany,                                               *
 *      LAPP, Annecy, France                                                      *
 *                                                                                *
 * Redistribution and use in source and binary forms, with or without             *
 * modification, are permitted according to the terms listed in LICENSE           *
 * (http://mva.sourceforge.net/license.txt)                                       *
 *                                                                                *
 * File and Version Information:                                                  *
 * $Id: TMVA_MethodBase.cxx,v 1.2 2006/05/09 08:37:06 brun Exp $
 **********************************************************************************/

//_______________________________________________________________________
//
// Virtual base class for all MVA method
//
//_______________________________________________________________________

#include "string"
#include "TMVA_MethodBase.h"
#include "TMVA_Event.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TMVA_Timer.h"
#include <stdlib.h>
#include "TMVA_Tools.h"
#include "TMVA_RootFinder.h"
#include "TMVA_PDF.h"
#include "TObjString.h"
#include "TQObject.h"
#include "TSpline.h"
#include "TMatrix.h"
#include "TMath.h"

#define DEBUG_TMVA_MethodBase          kFALSE
#define TMVA_MethodBase_MaxIterations_ 200
#define Use_Splines_for_Eff_          kTRUE
#define Thats_Big__                   1.0e30

#define NBIN_HIST_PLOT    100
#define NBIN_HIST_HIGH    10000

ClassImp(TMVA_MethodBase)

//_______________________________________________________________________
TMVA_MethodBase::TMVA_MethodBase( TString jobName,
				  vector<TString>* theVariables,
				  TTree*  theTree,
				  TString theOption,
				  TDirectory*  theBaseDir) :
  fJobName      ( jobName ),
  fTrainingTree ( theTree ),
  fInputVars    ( theVariables ),
  fOptions      ( theOption ),
  fBaseDir      ( theBaseDir ),
  fWeightFile   ( "" )
{
// default constructur
  this->Init();
  // parse option string and search for verbose
  // after that, remove the verbose option to not interfere with method-specific options
  TList*  list = TMVA_Tools::ParseFormatLine( fOptions );
  TString opt;
  for (Int_t i=0; i<list->GetSize(); i++) {
    TString s = ((TObjString*)list->At(i))->GetString();
    s.ToUpper();
    if (s == "V") {
      fVerbose = kTRUE;
      if (i == list->GetSize()-1) opt.Chop();
    }
    else {
      opt += (TString)((TObjString*)list->At(i))->GetString();
      if (i < list->GetSize()-1) opt += ":";
    }
  }
  fOptions = opt;

  // default extension for weight files
  fFileExtension = "weights";
  fFileDir       = "weights";
  gSystem->MakeDirectory( fFileDir );


  // init the normalization vectors
  InitNorm( fTrainingTree );
}

//_______________________________________________________________________
TMVA_MethodBase::TMVA_MethodBase( vector<TString> *theVariables,
				  TString weightFile,
				  TDirectory*  theBaseDir)
  : fJobName      ( "" ),
    fTrainingTree ( NULL ),
    fInputVars    ( theVariables ),
    fOptions      ( "" ),
    fBaseDir      ( theBaseDir ),
    fWeightFile   ( weightFile )
{
// constructor used for Testing + Application of the MVA, only (no training), using given WeightFiles
  this->Init();
  fJobName       = "";   //not used
}

//_______________________________________________________________________
 void TMVA_MethodBase::Init(){
  fVerbose       = kFALSE;
  fIsOK          = kTRUE;
  fNvar = fInputVars->size();
  fXminNorm      = 0;
  fXmaxNorm      = 0;
  fMeanS         = -1; // it is nice to have them "initialized". Every method
  fMeanB         = -1; // but "MethodCuts" sets them later
  fRmsS          = -1;
  fRmsB          = -1;

  fNbins         = NBIN_HIST_PLOT;
  fNbinsH        = NBIN_HIST_HIGH;

  fHistS_plotbin = NULL;
  fHistB_plotbin = NULL;
  fHistS_highbin = NULL;
  fHistB_highbin = NULL;
  fEffS          = NULL;
  fEffB          = NULL;
  fEffBvsS       = NULL;
  fRejBvsS       = NULL;
  fHistBhatS     = NULL;
  fHistBhatB     = NULL;
  fHistMuS       = NULL;
  fHistMuB       = NULL;
  fTestvarPrefix = "MVA_";

  // init variable bounds
  fXminNorm = new vector<Double_t>( fNvar );
  fXmaxNorm = new vector<Double_t>( fNvar );
  for (Int_t ivar=0; ivar<fNvar; ivar++) {
    (*fXminNorm)[ivar] = +Thats_Big__;
    (*fXmaxNorm)[ivar] = -Thats_Big__;
  }

  // define "this" pointer
  ResetThisBase();
}

//_______________________________________________________________________
 TMVA_MethodBase::~TMVA_MethodBase( void )
{
/// default destructur
  if (Verbose()) cout << "--- TMVA_MethodCuts: Destructor called " << endl;

   if (NULL != fXminNorm) delete fXminNorm;
   if (NULL != fXmaxNorm) delete fXmaxNorm;
}

//_______________________________________________________________________
 void TMVA_MethodBase::InitNorm( TTree* theTree )
{
  // if trainingsTree exists, fill min/max vector
  if (NULL != theTree) {
    for (Int_t ivar=0; ivar<fNvar; ivar++) {
      this->SetXminNorm( ivar, theTree->GetMinimum( (*fInputVars)[ivar] ) );
      this->SetXmaxNorm( ivar, theTree->GetMaximum( (*fInputVars)[ivar] ) );
    }
  }
  else {
    cout << "--- " << GetName()
	 << ":InitNorm Error: tree has zero pointer ==> abort" << endl;
    exit(1);
  }
  if (Verbose()) {
    cout << "--- " << GetName() << " <verbose>: set minNorm/maxNorm to: " << endl;
    cout << setprecision(3);
    for (Int_t ivar=0; ivar<fNvar; ivar++)
      cout << "    " << (*fInputVars)[ivar]
	   << "\t: [" << GetXminNorm( ivar ) << "\t, " << GetXmaxNorm( ivar ) << "\t] " << endl;
    cout << setprecision(5); // reset to better value
  }
}

//_______________________________________________________________________
 void TMVA_MethodBase::SetWeightFileName( void )
{
  fWeightFile =  fFileDir + "/" +fJobName + "_" + fMethodName + "." + fFileExtension;
}

//_______________________________________________________________________
 void TMVA_MethodBase::SetWeightFileName( TString theWeightFile)
{
  fWeightFile = theWeightFile;
}

//_______________________________________________________________________
 TString TMVA_MethodBase::GetWeightFileName( void )
{
  if (fWeightFile == "") this->SetWeightFileName();
  return fWeightFile;
}

//_______________________________________________________________________
 Bool_t TMVA_MethodBase::CheckSanity( TTree* theTree )
{
  // if no tree is given, use the trainingTree
  TTree* tree = (0 != theTree) ? theTree : fTrainingTree;

  // the input variables must exist in the tree
  vector<TString>::iterator itrVar    = fInputVars->begin();
  vector<TString>::iterator itrVarEnd = fInputVars->end();
  Bool_t found = kTRUE;
  for (; itrVar != itrVarEnd; itrVar++)
    if (0 == tree->FindBranch( *itrVar )) found = kFALSE;

  return found;
}

//_______________________________________________________________________
 void TMVA_MethodBase::AppendToMethodName( TString methodNameSuffix )
{
  fMethodName += "_";
  fTestvar += "_";
  fMethodName += methodNameSuffix;
  fTestvar += methodNameSuffix;
}

//_______________________________________________________________________
 void TMVA_MethodBase::SetWeightFileDir( TString fileDir )
{
  fFileDir = fileDir;
  gSystem->MakeDirectory( fFileDir );
}

// ---------------------------------------------------------------------------------------
// ----- methods related to renormalization of variables ---------------------------------
// ---------------------------------------------------------------------------------------

//_______________________________________________________________________
 Double_t TMVA_MethodBase::Norm( TString var, Double_t x ) const
{
  return TMVA_Tools::NormVariable( x, GetXminNorm( var ), GetXmaxNorm( var ) );
}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::Norm( Int_t ivar, Double_t x ) const
{
  return TMVA_Tools::NormVariable( x, GetXminNorm( ivar ), GetXmaxNorm( ivar ) );
}

//_______________________________________________________________________
 void TMVA_MethodBase::UpdateNorm( Int_t ivar, Double_t x )
{
  if (x < GetXminNorm( ivar )) SetXminNorm( ivar, x );
  if (x > GetXmaxNorm( ivar )) SetXmaxNorm( ivar, x );
}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetXminNorm( TString var ) const
{
  for (Int_t ivar=0; ivar<fNvar; ivar++)
    if (var == (*fInputVars)[ivar]) return (*fXminNorm)[ivar];

  cout << "--- " << GetName() << ": Error in ::GetXminNorm: variable not found ==> abort "
       << var << endl;
  exit(1);
}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetXmaxNorm( TString var ) const
{
  for (Int_t ivar=0; ivar<fNvar; ivar++)
    if (var == (*fInputVars)[ivar]) return (*fXmaxNorm)[ivar];

  cout << "--- " << GetName() << ": Error in ::GetXmaxNorm: variable not found ==> abort "
       << var << endl;
  exit(1);
}

//_______________________________________________________________________
 void TMVA_MethodBase::SetXminNorm( TString var, Double_t x )
{
  for (Int_t ivar=0; ivar<fNvar; ivar++) {
    if (var == (*fInputVars)[ivar]) {
      (*fXminNorm)[ivar] = x;
      return;
    }
  }

  cout << "--- " << GetName() << ": Error in ::SetXminNorm: variable not found ==> abort "
       << var << endl;
  exit(1);
}

//_______________________________________________________________________
 void TMVA_MethodBase::SetXmaxNorm( TString var, Double_t x )
{
  for (Int_t ivar=0; ivar<fNvar; ivar++) {
    if (var == (*fInputVars)[ivar]) {
      (*fXmaxNorm)[ivar] = x;
      return;
    }
  }

  cout << "--- " << GetName() << ": Error in ::SetXmaxNorm: variable not found ==> abort "
       << var << endl;
  exit(1);
}
// ---------------------------------------------------------------------------------------

//_______________________________________________________________________
 void TMVA_MethodBase::TestInit(TTree* theTestTree)
{

  //  fTestTree       = theTestTree;
  fHistS_plotbin  = fHistB_plotbin = 0;
  fHistS_highbin  = fHistB_highbin = 0;
  fEffS           = fEffB = fEffBvsS = fRejBvsS = 0;
  fGraphS         = fGraphB = 0;
  fEffSatB        = 0;
  fSeparation     = 0;
  fCutOrientation = kPositive;
  fSplS           = fSplB = 0;
  fSplRefS        = fSplRefB = 0;


  // sanity checks: tree must exist, and theVar must be in tree
  if (0 == theTestTree ||
      ( 0 == theTestTree->FindBranch( fTestvar ) && !(GetMethodName().Contains("Cuts")))){
    cout<<"--- "<< GetName() << ": Error in TestInit: test variable "<<fTestvar
	<<" not found in tree"<<endl;

    fIsOK = kFALSE;
  }

  // now call the TestInitLocal for possible individual initialisation
  // of each method
  this->TestInitLocal(theTestTree);

}

//_______________________________________________________________________
 void TMVA_MethodBase::PrepareEvaluationTree( TTree* testTree )
{
  // sanity checks
  if (0 == testTree) {
    cout << "--- " << GetName()
	 << ": PrepareEvaluationTree Error: testTree has zero pointer ==> exit(1)"
	 << endl;
    exit(1);
  }

  // checks that all variables in input vector indeed exist in the testTree
  if (!CheckSanity( testTree )) {
    cout << "--- " << GetName()
	 << ": PrepareEvaluationTree Error: sanity check failed" << endl;
    exit(1);
  }

  // read the coefficients
  this->ReadWeightsFromFile();

  // fill a new branch into the testTree with the MVA-value of the method
  Double_t myMVA;
  TBranch *newBranch = testTree->Branch( fTestvar, &myMVA, fTestvar + "/D" );

  // use timer
  TMVA_Timer timer( testTree->GetEntries(), GetName(), kTRUE );

  for (Int_t ievt=0; ievt<testTree->GetEntries(); ievt++) {
    if ((Int_t)ievt%100 == 0) timer.DrawProgressBar( ievt );
    TMVA_Event *e = new TMVA_Event( testTree, ievt, fInputVars );
    myMVA = this->GetMvaValue( e );
    newBranch->Fill();
    delete e;
  }
  cout << "--- " << GetName() << ": elapsed time for evaluation of "
       << testTree->GetEntries() <<  " events: "
       << timer.GetElapsedTime() << "       " << endl;
}

//_______________________________________________________________________
 void TMVA_MethodBase::Test( TTree *theTestTree )
{
  // basic statistics operations are made in base class
  // note: cannot directly modify private class members
  Double_t meanS, meanB, rmsS, rmsB, xmin, xmax;
  TMVA_Tools::ComputeStat( theTestTree, fTestvar, meanS, meanB, rmsS, rmsB, xmin, xmax );

  // choose reasonable histogram ranges, by removing outliers
  Double_t nrms = 4;
  xmin = TMath::Max( TMath::Min(meanS - nrms*rmsS, meanB - nrms*rmsB ), xmin );
  xmax = TMath::Min( TMath::Max(meanS + nrms*rmsS, meanB + nrms*rmsB ), xmax );

  fMeanS = meanS; fMeanB = meanB;
  fRmsS  = rmsS;  fRmsB  = rmsB;
  fXmin  = xmin;  fXmax  = xmax;

  // determine cut orientation
  fCutOrientation = (fMeanS > fMeanB) ? kPositive : kNegative;

  // fill 2 types of histograms for the various analyses
  // this one is for actual plotting
  fHistS_plotbin = TMVA_Tools::projNormTH1F( theTestTree, fTestvar,
					     fTestvar + "_S",
					     fNbins, fXmin, fXmax, "type == 1" );
  fHistB_plotbin = TMVA_Tools::projNormTH1F( theTestTree, fTestvar,
					     fTestvar + "_B",
					     fNbins, fXmin, fXmax, "type == 0" );

  // need histograms with even more bins for efficiency calculation and integration
  fHistS_highbin = TMVA_Tools::projNormTH1F( theTestTree, fTestvar,
					     fTestvar + "_S_high",
					     fNbinsH, fXmin, fXmax, "type == 1" );
  fHistB_highbin = TMVA_Tools::projNormTH1F( theTestTree, fTestvar,
					     fTestvar + "_B_high",
					     fNbinsH, fXmin, fXmax, "type == 0" );

  // create PDFs from histograms, using default splines, and no additional smoothing
  fSplS = new TMVA_PDF( fHistS_plotbin, TMVA_PDF::kSpline2, 0 );
  fSplB = new TMVA_PDF( fHistB_plotbin, TMVA_PDF::kSpline2, 0  );

}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetEfficiency( TString theString, TTree *theTree )
{
  // parse input string for required background efficiency
  TList*  list  = TMVA_Tools::ParseFormatLine( theString );
  // sanity check

  if (list->GetSize() != 2) {
    cout << "--- " << GetName() << ": Error in::GetEfficiency: wrong number of arguments"
	 << " in string: " << theString
	 << " | required format, e.g., Efficiency:0.05" << endl;
    return -1;
  }
  //that will be the value of the efficiency retured (does not affect
  //the efficiency-vs-bkg plot which is done anyway.
  Float_t effBref  = atof( ((TObjString*)list->At(1))->GetString() );

  if (DEBUG_TMVA_MethodBase)
    cout << "--- " << GetName() << "::GetEfficiency : compute eff(S) at eff(B) = " << effBref << endl;

  // sanity check
  if (fHistS_highbin->GetNbinsX() != fHistB_highbin->GetNbinsX() ||
      fHistS_plotbin->GetNbinsX() != fHistB_plotbin->GetNbinsX()) {
    cout << "--- " << GetName()
	 << "WARNING: in GetEfficiency() binning mismatch between signal and background histos"<<endl;
    fIsOK = kFALSE;
    return -1.0;
  }

  // create histogram

  // first, get efficiency histograms for signal and background
  Double_t xmin = fHistS_highbin->GetXaxis()->GetXmin();
  Double_t xmax = fHistS_highbin->GetXaxis()->GetXmax();

  // first round ? --> create histograms
  Bool_t firstPass = kFALSE;
  if (NULL == fEffS && NULL == fEffB) firstPass = kTRUE;

  if (firstPass) {

    fEffS = new TH1F( fTestvar + "_effS", fTestvar + " (signal)",     fNbinsH, xmin, xmax );
    fEffB = new TH1F( fTestvar + "_effB", fTestvar + " (background)", fNbinsH, xmin, xmax );

    // sign if cut
    Int_t sign = (fCutOrientation == kPositive) ? +1 : -1;

    // this method is unbinned
    for (Int_t ievt=0; ievt<theTree->GetEntries(); ievt++) {

      TH1* theHist = 0;
      if ((Int_t)TMVA_Tools::GetValue( theTree, ievt, "type" ) == 1) { // this is signal
	theHist = fEffS;
      }
      else { // this is background
	theHist = fEffB;
      }

      Double_t theVal = TMVA_Tools::GetValue( theTree, ievt, fTestvar );
      for (Int_t bin=1; bin<=fNbinsH; bin++)
	if (sign*theVal > sign*theHist->GetBinCenter( bin )) theHist->AddBinContent( bin );
    }

    // renormalize to maximum
    fEffS->Scale( 1.0/(fEffS->GetMaximum() > 0 ? fEffS->GetMaximum() : 1) );
    fEffB->Scale( 1.0/(fEffB->GetMaximum() > 0 ? fEffB->GetMaximum() : 1) );

    // now create efficiency curve: background versus signal
    fEffBvsS = new TH1F( fTestvar + "_effBvsS", fTestvar + "", fNbins, 0, 1 );
    fRejBvsS = new TH1F( fTestvar + "_rejBvsS", fTestvar + "", fNbins, 0, 1 );
    // use root finder
    // spline background efficiency plot
    // note that there is a bin shift when going from a TH1F object to a TGraph :-(
    if (Use_Splines_for_Eff_) {
      fGraphS   = new TGraph( fEffS );
      fGraphB   = new TGraph( fEffB );
      fSplRefS  = new TMVA_TSpline1( "spline2_signal",     fGraphS );
      fSplRefB  = new TMVA_TSpline1( "spline2_background", fGraphB );

      // verify spline sanity
      if (Verbose())
	cout << "--- " << GetName()
	     << "::GetEfficiency <verbose>: verify signal and background eff. splines" << endl;
      TMVA_Tools::CheckSplines( fEffS, fSplRefS );
      TMVA_Tools::CheckSplines( fEffB, fSplRefB );
    }

    // make the background-vs-signal efficiency plot

    // create root finder
    // reset static "this" pointer before calling external function
    ResetThisBase();
    TMVA_RootFinder rootFinder( &IGetEffForRoot, fXmin, fXmax );

    Double_t effB = 0;
    for (Int_t bini=1; bini<=fNbins; bini++) {

      // find cut value corresponding to a given signal efficiency
      Double_t effS = fEffBvsS->GetBinCenter( bini );

      Double_t cut  = rootFinder.Root( effS );

      // retrieve background efficiency for given cut
      if (Use_Splines_for_Eff_)
	effB = fSplRefB->Eval( cut );
      else
	effB = fEffB->GetBinContent( fEffB->FindBin( cut ) );

      // and fill histograms
      fEffBvsS->SetBinContent( bini, effB     );
      fRejBvsS->SetBinContent( bini, 1.0-effB );
    }

    // create splines for histogram
    fGrapheffBvsS = new TGraph( fEffBvsS );
    fSpleffBvsS   = new TMVA_TSpline1( "effBvsS", fGrapheffBvsS );
  }

  // must exist...
  if (NULL == fSpleffBvsS) return 0.0;

  // now find signal efficiency that corresponds to required background efficiency
  Double_t effS, effB, effS_ = 0, effB_ = 0;
  Int_t    nbins_ = 1000;
  for (Int_t bini=1; bini<=nbins_; bini++) {

    // get corresponding signal and background efficiencies
    effS = (bini - 0.5)/Float_t(nbins_);
    effB = fSpleffBvsS->Eval( effS );

    // find signal efficiency that corresponds to required background efficiency
    if ((effB - effBref)*(effB_ - effBref) < 0) break;
    effS_ = effS;
    effB_ = effB;
  }

  return fEffSatB = 0.5*(effS + effS_);//the mean between bin above and bin below

}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetSignificance( void )
{
  // compute significance of mean difference
  // significance = |<S> - <B>|/Sqrt(RMS_S2 + RMS_B2)
  Double_t rms = sqrt(pow(fRmsS,2) + pow(fRmsB,2));

  return fSignificance = (rms > 0) ? fabs(fMeanS - fMeanB)/rms : 0;
}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetSeparation( void )
{
  // compute "separation" defined as
  // <s2> = (1/2) Int_-oo..+oo { (S(x)2 - B(x)2)/(S(x) + B(x)) dx }
  fSeparation = 0;

  Int_t nstep  = 1000;
  Double_t intBin = (fXmax - fXmin)/nstep;
  for (Int_t bin=0; bin<nstep; bin++) {
    Double_t x = (bin + 0.5)*intBin + fXmin;
    Double_t s = fSplS->GetVal( x );
    Double_t b = fSplB->GetVal( x );
    // separation
    if (s + b > 0) fSeparation += 0.5*pow(s - b,2)/(s + b);
  }
  fSeparation *= intBin;

  return fSeparation;
}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetmuTransform( TTree *theTree )
{
  //---------------------------------------------------------------------------------------
  // Authors     : Francois Le Diberder and Muriel Pivk
  // Reference   : Muriel Pivk,
  //               "Etude de la violation de CP dans la désintégration
  //                B0 -> h+ h- (h = pi, K) auprès du détecteur BaBar à SLAC",
  //               PhD thesis at Universite de Paris VI-VII, LPNHE (IN2P3/CNRS), Paris, 2003
  //               http://tel.ccsd.cnrs.fr/documents/archives0/00/00/29/91/index_fr.html
  //
  // Definitions : Bhat = PDFbackground(x)/(PDFbackground(x) + PDFsignal(x))
  //               mu   = mu(b) = Int_0B Bhat[b'] db'
  //---------------------------------------------------------------------------------------

  // create Bhat distribution function
  Int_t nbin  = 70;
  fHistBhatS = new TH1F( fTestvar + "_BhatS", fTestvar + ": Bhat (S)", nbin, 0.0, 1.0 );
  fHistBhatB = new TH1F( fTestvar + "_BhatB", fTestvar + ": Bhat (B)", nbin, 0.0, 1.0 );

  fHistBhatS->Sumw2();
  fHistBhatB->Sumw2();

  vector<Double_t>* aBhatB = new vector<Double_t>;
  vector<Double_t>* aBhatS = new vector<Double_t>;
  Int_t ievt;
  for (ievt=0; ievt<theTree->GetEntries(); ievt++) {
    Double_t x    = TMVA_Tools::GetValue( theTree, ievt, fTestvar );
    Double_t s    = fSplS->GetVal( x );
    Double_t b    = fSplB->GetVal( x );
    Double_t aBhat = 0;
    if (b + s > 0) aBhat = b/(b + s);

    if ((Int_t)TMVA_Tools::GetValue( theTree, ievt, "type" ) == 1) { // this is signal
      aBhatS->push_back ( aBhat );
      fHistBhatS->Fill( aBhat );
    }
    else {
      aBhatB->push_back ( aBhat );
      fHistBhatB->Fill( aBhat );
    }
  }

  // normalize histograms
  fHistBhatS->Scale( 1.0/((fHistBhatS->GetEntries() > 0 ? fHistBhatS->GetEntries() : 1) / nbin) );
  fHistBhatB->Scale( 1.0/((fHistBhatB->GetEntries() > 0 ? fHistBhatB->GetEntries() : 1) / nbin) );

  TMVA_PDF* yB = new TMVA_PDF( fHistBhatB, TMVA_PDF::kSpline2, 100 );

  Int_t nevtS = aBhatS->size();
  Int_t nevtB = aBhatB->size();

  // get the mu-transform
  Int_t nbinMu = 50;
  fHistMuS = new TH1F( fTestvar + "_muTransform_S",
			fTestvar + ": mu-Transform (S)", nbinMu, 0.0, 1.0 );
  fHistMuB = new TH1F( fTestvar + "_muTransform_B",
			fTestvar + ": mu-Transform (B)", nbinMu, 0.0, 1.0 );

  // signal
  for (ievt=0; ievt<nevtS; ievt++) {
    Double_t w = yB->GetVal( (*aBhatS)[ievt] );
    if (w > 0) fHistMuS->Fill( 1.0 - (*aBhatS)[ievt], 1.0/w );
  }

  // background (must be flat)
  for (ievt=0; ievt<nevtB; ievt++) {
    Double_t w = yB->GetVal( (*aBhatB)[ievt] );
    if (w > 0) fHistMuB->Fill( 1.0 - (*aBhatB)[ievt], 1.0/w );
  }

  // normalize mu-transforms
  TMVA_Tools::NormHist( fHistMuS );
  TMVA_Tools::NormHist( fHistMuB );

  // determine the mu-transform value, which is defined as
  // the average of the signal mu-transform Int_[0,1] { S(mu) dmu }
  // this average is 0.5 for background, by definition
  TMVA_PDF* thePdf = new TMVA_PDF( fHistMuS, TMVA_PDF::kSpline2 );
  Double_t intS = 0;
  Int_t    nstp = 10000;
  for (Int_t istp=0; istp<nstp; istp++) {
    Double_t x = (istp + 0.5)/Double_t(nstp);
    intS += x*thePdf->GetVal( x );
  }
  intS /= Double_t(nstp);

  delete yB;
  delete thePdf;
  delete aBhatB;
  delete aBhatS;

  return intS; // return average mu-transform for signal
}

//_______________________________________________________________________
 void TMVA_MethodBase::WriteHistosToFile( TDirectory* targetDir )
{
  targetDir->cd();
  if (0 != fHistS_plotbin) fHistS_plotbin->Write();
  if (0 != fHistB_plotbin) fHistB_plotbin->Write();
  if (0 != fHistS_highbin) fHistS_highbin->Write();
  if (0 != fHistB_highbin) fHistB_highbin->Write();
  if (0 != fEffS         ) fEffS->Write();
  if (0 != fEffB         ) fEffB->Write();
  if (0 != fEffBvsS      ) fEffBvsS->Write();
  if (0 != fRejBvsS      ) fRejBvsS->Write();
  if (0 != fHistBhatS    ) fHistBhatS->Write();
  if (0 != fHistBhatB    ) fHistBhatB->Write();
  if (0 != fHistMuS      ) fHistMuS->Write();
  if (0 != fHistMuB      ) fHistMuB->Write();
}

// ----------------------- r o o t   f i n d i n g ----------------------------

TMVA_MethodBase* TMVA_MethodBase::fgThisBase = NULL;

//_______________________________________________________________________
 Double_t TMVA_MethodBase::IGetEffForRoot( Double_t theCut )
{
  return TMVA_MethodBase::GetThisBase()->GetEffForRoot( theCut );
}

//_______________________________________________________________________
 Double_t TMVA_MethodBase::GetEffForRoot( Double_t theCut )
{
  Double_t retval;

  // retrieve the class object
  if (Use_Splines_for_Eff_)
    retval = fSplRefS->Eval( theCut );
  else
    retval = fEffS->GetBinContent( fEffS->FindBin( theCut ) );

  // caution: here we take some "forbidden" action to hide a problem:
  // in some cases, in particular for likelihood, the binned efficiency distributions
  // do not equal 1, at xmin, and 0 at xmax; of course, in principle we have the
  // unbinned information available in the trees, but the unbinned minimization is
  // too slow, and we don't need to do a precision measurement here. Hence, we force
  // this property.
  Double_t eps = 1.0e-5;
  if      (theCut-fXmin < eps) retval = (GetCutOrientation() == kPositive) ? 1.0 : 0.0;
  else if (fXmax-theCut < eps) retval = (GetCutOrientation() == kPositive) ? 0.0 : 1.0;

  return retval;
};



ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.