BayesianAnalysisToolkit: BCIntegrate.cxx Source File

Go to the documentation of this file.
00001 /*
00002  * Copyright (C) 2008-2011, Daniel Kollar and Kevin Kroeninger.
00003  * All rights reserved.
00004  *
00005  * For the licensing terms see doc/COPYING.
00006  */
00007 
00008 // ---------------------------------------------------------
00009 
00010 #include "config.h"
00011 
00012 #include "BCIntegrate.h"
00013 #include "BCLog.h"
00014 #include "BCMath.h"
00015 
00016 #include <TH1D.h>
00017 #include <TH2D.h>
00018 #include <TMinuit.h>
00019 #include <TRandom3.h>
00020 #include <TTree.h>
00021 
00022 #ifdef HAVE_CUBA_H
00023 #include "cuba.h"
00024 #endif
00025 
00026 // ---------------------------------------------------------
00027 
00028 class BCIntegrate * global_this;
00029 
00030 // ---------------------------------------------------------
00031 BCIntegrate::BCIntegrate()
00032    : BCEngineMCMC()
00033    , fNvar(0)
00034    , fNbins(100)
00035    , fNSamplesPer2DBin(100)
00036    , fMarkovChainStepSize(0.1)
00037    , fMarkovChainAutoN(true)
00038    , fDataPointLowerBoundaries(0)
00039    , fDataPointUpperBoundaries(0)
00040    , fFillErrorBand(false)
00041    , fFitFunctionIndexX(-1)
00042    , fFitFunctionIndexY(-1)
00043    , fErrorBandContinuous(true)
00044    , fErrorBandNbinsX(100)
00045    , fErrorBandNbinsY(500)
00046    , fMinuit(0)
00047    , fFlagIgnorePrevOptimization(false)
00048    , fFlagWriteMarkovChain(false)
00049    , fMarkovChainTree(0)
00050    , fSAT0(100)
00051    , fSATmin(0.1)
00052    , fTreeSA(0)
00053    , fFlagWriteSAToFile(false)
00054 
00055    , fNiterPerDimension(100)
00056 #ifdef HAVE_CUBA_H
00057    , fIntegrationMethod(BCIntegrate::kIntCuba)
00058 #else
00059    , fIntegrationMethod(BCIntegrate::kIntMonteCarlo)
00060 #endif
00061    , fMarginalizationMethod(BCIntegrate::kMargMetropolis)
00062    , fOptimizationMethod(BCIntegrate::kOptMinuit)
00063    , fOptimizationMethodMode(BCIntegrate::kOptMinuit)
00064    , fSASchedule(BCIntegrate::kSACauchy)
00065    , fNIterationsMax(1000000)
00066    , fNIterations(0)
00067    , fRelativePrecision(1e-3)
00068    , fAbsolutePrecision(1e-12)
00069    , fCubaIntegrationMethod(BCIntegrate::kCubaVegas)
00070    , fCubaMinEval(0)
00071    , fCubaMaxEval(2000000)
00072    , fCubaVerbosity(0)
00073    , fCubaVegasNStart(25000)
00074    , fCubaVegasNIncrease(25000)
00075    , fCubaSuaveNNew(10000)
00076    , fCubaSuaveFlatness(50)
00077    , fError(-999.)
00078 {
00079    fMinuitArglist[0] = 20000;
00080    fMinuitArglist[1] = 0.01;
00081 }
00082 
00083 // ---------------------------------------------------------
00084 BCIntegrate::BCIntegrate(BCParameterSet * par)
00085    : BCEngineMCMC()
00086    , fNvar(0)
00087    , fNbins(100)
00088    , fNSamplesPer2DBin(100)
00089    , fMarkovChainStepSize(0.1)
00090    , fMarkovChainAutoN(true)
00091    , fDataPointLowerBoundaries(0)
00092    , fDataPointUpperBoundaries(0)
00093    , fFillErrorBand(false)
00094    , fFitFunctionIndexX(-1)
00095    , fFitFunctionIndexY(-1)
00096    , fErrorBandContinuous(true)
00097    , fErrorBandNbinsX(100)
00098    , fErrorBandNbinsY(500)
00099    , fMinuit(0)
00100    , fFlagIgnorePrevOptimization(false)
00101    , fFlagWriteMarkovChain(false)
00102    , fMarkovChainTree(0)
00103    , fSAT0(100)
00104    , fSATmin(0.1)
00105    , fTreeSA(0)
00106    , fFlagWriteSAToFile(false)
00107 
00108    , fNiterPerDimension(100)
00109 #ifdef HAVE_CUBA_H
00110    , fIntegrationMethod(BCIntegrate::kIntCuba)
00111 #else
00112    , fIntegrationMethod(BCIntegrate::kIntMonteCarlo)
00113 #endif
00114    , fMarginalizationMethod(BCIntegrate::kMargMetropolis)
00115    , fOptimizationMethod(BCIntegrate::kOptMinuit)
00116    , fOptimizationMethodMode(BCIntegrate::kOptMinuit)
00117    , fSASchedule(BCIntegrate::kSACauchy)
00118    , fNIterationsMax(1000000)
00119    , fNIterations(0)
00120    , fRelativePrecision(1e-3)
00121    , fAbsolutePrecision(1e-12)
00122    , fCubaIntegrationMethod(BCIntegrate::kCubaVegas)
00123    , fCubaMinEval(0)
00124    , fCubaMaxEval(2000000)
00125    , fCubaVerbosity(0)
00126    , fCubaVegasNStart(25000)
00127    , fCubaVegasNIncrease(25000)
00128    , fCubaSuaveNNew(10000)
00129    , fCubaSuaveFlatness(50)
00130    , fError(-999.)
00131 {
00132    SetParameters(par);
00133 
00134    fMinuitArglist[0] = 20000;
00135    fMinuitArglist[1] = 0.01;
00136 }
00137 
00138 // ---------------------------------------------------------
00139 BCIntegrate::BCIntegrate(const BCIntegrate & bcintegrate) : BCEngineMCMC(bcintegrate)
00140  {
00141     fNvar                     = bcintegrate.fNvar;
00142     fNbins                    = bcintegrate.fNbins;
00143     fNSamplesPer2DBin         = bcintegrate.fNSamplesPer2DBin;
00144     fMarkovChainStepSize      = bcintegrate.fMarkovChainStepSize;
00145     fMarkovChainNIterations   = bcintegrate.fMarkovChainNIterations;
00146     fMarkovChainAutoN         = bcintegrate.fMarkovChainAutoN;
00147     if (bcintegrate.fDataPointLowerBoundaries)
00148        fDataPointLowerBoundaries = new BCDataPoint(*bcintegrate.fDataPointLowerBoundaries);
00149     else
00150        fDataPointLowerBoundaries = 0;
00151     if (bcintegrate.fDataPointUpperBoundaries)
00152        fDataPointUpperBoundaries = new BCDataPoint(*bcintegrate.fDataPointUpperBoundaries);
00153     else
00154        fDataPointUpperBoundaries = 0;
00155     fDataFixedValues          = bcintegrate.fDataFixedValues;
00156     fBestFitParameters        = bcintegrate.fBestFitParameters;
00157     fBestFitParameterErrors   = bcintegrate.fBestFitParameterErrors;
00158     fBestFitParametersMarginalized = bcintegrate.fBestFitParametersMarginalized;
00159     for (int i = 0; i < int(bcintegrate.fHProb1D.size()); ++i) {
00160        if (bcintegrate.fHProb1D.at(i))
00161           fHProb1D.push_back(new TH1D(*(bcintegrate.fHProb1D.at(i))));
00162        else
00163           fHProb1D.push_back(0);
00164     }
00165     for (int i = 0; i < int(bcintegrate.fHProb2D.size()); ++i) {
00166        if (bcintegrate.fHProb2D.at(i))
00167           fHProb2D.push_back(new TH2D(*(fHProb2D.at(i))));
00168        else
00169           fHProb2D.push_back(0);
00170     }
00171     fFillErrorBand            = bcintegrate.fFillErrorBand;
00172     fFitFunctionIndexX        = bcintegrate.fFitFunctionIndexX;
00173     fFitFunctionIndexY        = bcintegrate.fFitFunctionIndexY;
00174     fErrorBandX               = bcintegrate.fErrorBandX;
00175     if (bcintegrate.fErrorBandXY)
00176        fErrorBandXY = new TH2D(*(bcintegrate.fErrorBandXY));
00177     else
00178        fErrorBandXY = 0;
00179     fErrorBandNbinsX          = bcintegrate.fErrorBandNbinsX;
00180     fErrorBandNbinsY          = bcintegrate.fErrorBandNbinsY;
00181     fMinuit                   = new TMinuit();
00182     // debugKK
00183     //    *fMinuit = *(bcintegrate.fMinuit);
00184     fMinuitArglist[0]         = bcintegrate.fMinuitArglist[0];
00185     fMinuitArglist[1]         = bcintegrate.fMinuitArglist[1];
00186     fMinuitErrorFlag          = bcintegrate.fMinuitErrorFlag;
00187     fFlagIgnorePrevOptimization = bcintegrate.fFlagIgnorePrevOptimization;
00188     fFlagWriteMarkovChain     = bcintegrate.fFlagWriteMarkovChain;
00189     fMarkovChainTree          = bcintegrate.fMarkovChainTree;
00190     fMCMCIteration            = bcintegrate.fMCMCIteration;
00191     fSAT0                     = bcintegrate.fSAT0;
00192     fSATmin                   = bcintegrate.fSATmin; 
00193     // debugKK
00194     fTreeSA = 0;
00195     fFlagWriteSAToFile        = bcintegrate.fFlagWriteSAToFile;
00196     fSANIterations            = bcintegrate.fSANIterations;
00197     fSATemperature            = bcintegrate.fSATemperature;
00198     fSALogProb                = bcintegrate.fSALogProb;
00199     fSAx                      = bcintegrate.fSAx;
00200     if (bcintegrate.fx)
00201        fx = new BCParameterSet(*(bcintegrate.fx));
00202     else
00203        fx = 0;
00204     fMin                      = new double[fNvar]; 
00205     fMax                      = new double[fNvar];
00206     fVarlist                  = new int[fNvar];
00207     fMin                      = bcintegrate.fMin;
00208     fMax                      = bcintegrate.fMax;
00209     fVarlist                  = bcintegrate.fVarlist;
00210     fNiterPerDimension        = bcintegrate.fNiterPerDimension;
00211     fIntegrationMethod        = bcintegrate.fIntegrationMethod;
00212     fMarginalizationMethod    = bcintegrate.fMarginalizationMethod;
00213     fOptimizationMethod       = bcintegrate.fOptimizationMethod;
00214     fOptimizationMethodMode   = bcintegrate.fOptimizationMethodMode;
00215     fSASchedule               = bcintegrate.fSASchedule;
00216     fNIterationsMax           = bcintegrate.fNIterationsMax;
00217     fNIterations              = bcintegrate.fNIterations;
00218     fRelativePrecision        = bcintegrate.fRelativePrecision;
00219     fAbsolutePrecision        = bcintegrate.fAbsolutePrecision;
00220     fCubaIntegrationMethod    = bcintegrate.fCubaIntegrationMethod;
00221     fCubaMinEval              = bcintegrate.fCubaMinEval;
00222     fCubaMaxEval              = bcintegrate.fCubaMaxEval;
00223     fCubaVerbosity            = bcintegrate.fCubaVerbosity;
00224     fCubaVegasNStart          = bcintegrate.fCubaVegasNStart;
00225     fCubaVegasNIncrease       = bcintegrate.fCubaVegasNIncrease;
00226     fCubaSuaveNNew            = bcintegrate.fCubaSuaveNNew;
00227     fCubaSuaveFlatness        = bcintegrate.fCubaSuaveFlatness;
00228     fError                    = bcintegrate.fError;
00229     fNmetro                   = bcintegrate.fNmetro;
00230     fNacceptedMCMC            = bcintegrate.fNacceptedMCMC;
00231     fXmetro0                  = bcintegrate.fXmetro0;
00232     fXmetro1                  = bcintegrate.fXmetro1;
00233     fMarkovChainValue         = bcintegrate.fMarkovChainValue;
00234 }
00235 
00236 // ---------------------------------------------------------
00237 BCIntegrate & BCIntegrate::operator = (const BCIntegrate & bcintegrate)
00238 {
00239     BCEngineMCMC::operator=(bcintegrate);
00240 
00241     fNvar                     = bcintegrate.fNvar;
00242     fNbins                    = bcintegrate.fNbins;
00243     fNSamplesPer2DBin         = bcintegrate.fNSamplesPer2DBin;
00244     fMarkovChainStepSize      = bcintegrate.fMarkovChainStepSize;
00245     fMarkovChainNIterations   = bcintegrate.fMarkovChainNIterations;
00246     fMarkovChainAutoN         = bcintegrate.fMarkovChainAutoN;
00247     if (bcintegrate.fDataPointLowerBoundaries)
00248        fDataPointLowerBoundaries = new BCDataPoint(*bcintegrate.fDataPointLowerBoundaries);
00249     else
00250        fDataPointLowerBoundaries = 0;
00251     if (bcintegrate.fDataPointUpperBoundaries)
00252        fDataPointUpperBoundaries = new BCDataPoint(*bcintegrate.fDataPointUpperBoundaries);
00253     else
00254        fDataPointUpperBoundaries = 0;
00255     fDataFixedValues          = bcintegrate.fDataFixedValues;
00256     fBestFitParameters        = bcintegrate.fBestFitParameters;
00257     fBestFitParameterErrors   = bcintegrate.fBestFitParameterErrors;
00258     fBestFitParametersMarginalized = bcintegrate.fBestFitParametersMarginalized;
00259     for (int i = 0; i < int(bcintegrate.fHProb1D.size()); ++i) {
00260        if (bcintegrate.fHProb1D.at(i))
00261           fHProb1D.push_back(new TH1D(*(bcintegrate.fHProb1D.at(i))));
00262        else
00263           fHProb1D.push_back(0);
00264     }
00265     for (int i = 0; i < int(bcintegrate.fHProb2D.size()); ++i) {
00266        if (bcintegrate.fHProb2D.at(i))
00267           fHProb2D.push_back(new TH2D(*(fHProb2D.at(i))));
00268        else
00269           fHProb2D.push_back(0);
00270     }
00271     fFillErrorBand            = bcintegrate.fFillErrorBand;
00272     fFitFunctionIndexX        = bcintegrate.fFitFunctionIndexX;
00273     fFitFunctionIndexY        = bcintegrate.fFitFunctionIndexY;
00274     fErrorBandX               = bcintegrate.fErrorBandX;
00275     if (bcintegrate.fErrorBandXY)
00276        fErrorBandXY = new TH2D(*(bcintegrate.fErrorBandXY));
00277     else
00278        fErrorBandXY = 0;
00279     fErrorBandNbinsX          = bcintegrate.fErrorBandNbinsX;
00280     fErrorBandNbinsY          = bcintegrate.fErrorBandNbinsY;
00281     fMinuit                   = new TMinuit();
00282     // debugKK
00283     //    *fMinuit = *(bcintegrate.fMinuit);
00284     fMinuitArglist[0]         = bcintegrate.fMinuitArglist[0];
00285     fMinuitArglist[1]         = bcintegrate.fMinuitArglist[1];
00286     fMinuitErrorFlag          = bcintegrate.fMinuitErrorFlag;
00287     fFlagIgnorePrevOptimization = bcintegrate.fFlagIgnorePrevOptimization;
00288     fFlagWriteMarkovChain     = bcintegrate.fFlagWriteMarkovChain;
00289     fMarkovChainTree          = bcintegrate.fMarkovChainTree;
00290     fMCMCIteration            = bcintegrate.fMCMCIteration;
00291     fSAT0                     = bcintegrate.fSAT0;
00292     fSATmin                   = bcintegrate.fSATmin;
00293     // debugKK
00294     fTreeSA = 0;
00295     fFlagWriteSAToFile        = bcintegrate.fFlagWriteSAToFile;
00296     fSANIterations            = bcintegrate.fSANIterations;
00297     fSATemperature            = bcintegrate.fSATemperature;
00298     fSALogProb                = bcintegrate.fSALogProb;
00299     fSAx                      = bcintegrate.fSAx;
00300     if (bcintegrate.fx)
00301        fx = new BCParameterSet(*(bcintegrate.fx));
00302     else
00303        fx = 0;
00304     fMin                      = new double[fNvar];
00305     fMax                      = new double[fNvar];
00306     fVarlist                  = new int[fNvar];
00307     fNiterPerDimension        = bcintegrate.fNiterPerDimension;
00308     fIntegrationMethod        = bcintegrate.fIntegrationMethod;
00309     fMarginalizationMethod    = bcintegrate.fMarginalizationMethod;
00310     fOptimizationMethod       = bcintegrate.fOptimizationMethod;
00311     fOptimizationMethodMode   = bcintegrate.fOptimizationMethodMode;
00312     fSASchedule               = bcintegrate.fSASchedule;
00313 
00314     fNIterationsMax           = bcintegrate.fNIterationsMax;
00315     fNIterations              = bcintegrate.fNIterations;
00316     fRelativePrecision        = bcintegrate.fRelativePrecision;
00317     fAbsolutePrecision        = bcintegrate.fAbsolutePrecision;
00318     fCubaIntegrationMethod    = bcintegrate.fCubaIntegrationMethod;
00319     fCubaMinEval              = bcintegrate.fCubaMinEval;
00320     fCubaMaxEval              = bcintegrate.fCubaMaxEval;
00321     fCubaVerbosity            = bcintegrate.fCubaVerbosity;
00322     fCubaVegasNStart          = bcintegrate.fCubaVegasNStart;
00323     fCubaVegasNIncrease       = bcintegrate.fCubaVegasNIncrease;
00324     fCubaSuaveNNew            = bcintegrate.fCubaSuaveNNew;
00325     fCubaSuaveFlatness        = bcintegrate.fCubaSuaveFlatness;
00326     fError                    = bcintegrate.fError;
00327     fNmetro                   = bcintegrate.fNmetro;
00328     fNacceptedMCMC            = bcintegrate.fNacceptedMCMC;
00329     fXmetro0                  = bcintegrate.fXmetro0;
00330     fXmetro1                  = bcintegrate.fXmetro1;
00331     fMarkovChainValue         = bcintegrate.fMarkovChainValue;
00332 
00333    // return this
00334    return *this;
00335 }
00336 
00337 // ---------------------------------------------------------
00338 BCIntegrate::~BCIntegrate()
00339 {
00340    DeleteVarList();
00341 
00342    fx=0;
00343 
00344    if (fMinuit)
00345       delete fMinuit;
00346 
00347    int n1 = fHProb1D.size();
00348    if(n1>0) {
00349       for (int i=0;i<n1;i++)
00350          delete fHProb1D.at(i);
00351    }
00352 
00353    int n2 = fHProb2D.size();
00354    if(n2>0) {
00355       for (int i=0;i<n2;i++)
00356          delete fHProb2D.at(i);
00357    }
00358 }
00359 
00360 // ---------------------------------------------------------
00361 void BCIntegrate::SetParameters(BCParameterSet * par)
00362 {
00363    DeleteVarList();
00364 
00365    fx = par;
00366    fNvar = fx->size();
00367    fMin = new double[fNvar];
00368    fMax = new double[fNvar];
00369    fVarlist = new int[fNvar];
00370 
00371    ResetVarlist(1);
00372 
00373    for(int i=0;i<fNvar;i++) {
00374       fMin[i]=(fx->at(i))->GetLowerLimit();
00375       fMax[i]=(fx->at(i))->GetUpperLimit();
00376    }
00377 
00378    fXmetro0.clear();
00379    for(int i=0;i<fNvar;i++)
00380       fXmetro0.push_back((fMin[i]+fMax[i])/2.0);
00381 
00382    fXmetro1.clear();
00383    fXmetro1.assign(fNvar,0.);
00384 
00385    fMCMCBoundaryMin.clear();
00386    fMCMCBoundaryMax.clear();
00387 
00388    for(int i=0;i<fNvar;i++) {
00389       fMCMCBoundaryMin.push_back(fMin[i]);
00390       fMCMCBoundaryMax.push_back(fMax[i]);
00391       fMCMCFlagsFillHistograms.push_back(true);
00392    }
00393 
00394    for (int i = int(fMCMCH1NBins.size()); i<fNvar; ++i)
00395       fMCMCH1NBins.push_back(100);
00396 
00397    fMCMCNParameters = fNvar;
00398 }
00399 
00400 // ---------------------------------------------------------
00401 void BCIntegrate::SetMarkovChainInitialPosition(std::vector<double> position)
00402 {
00403    int n = position.size();
00404 
00405    fXmetro0.clear();
00406 
00407    for (int i = 0; i < n; ++i)
00408       fXmetro0.push_back(position.at(i));
00409 }
00410 
00411 // ---------------------------------------------------------
00412 void BCIntegrate::DeleteVarList()
00413 {
00414    if(fNvar) {
00415       delete[] fVarlist;
00416       fVarlist=0;
00417 
00418       delete[] fMin;
00419       fMin=0;
00420 
00421       delete[] fMax;
00422       fMax=0;
00423 
00424       fx=0;
00425       fNvar=0;
00426    }
00427 }
00428 
00429 // ---------------------------------------------------------
00430 void BCIntegrate::SetNbins(int nbins, int index)
00431 {
00432    if (fNvar == 0)
00433       return;
00434 
00435    // check if index is in range
00436    if (index >= fNvar) {
00437       BCLog::OutWarning("BCIntegrate::SetNbins : Index out of range.");
00438       return;
00439    }
00440    // set for all parameters at once
00441    else if (index < 0) {
00442       for (int i = 0; i < fNvar; ++i)
00443          SetNbins(nbins, i);
00444       return;
00445    }
00446 
00447    // sanity check for nbins
00448    if (nbins <= 0)
00449       nbins = 10;
00450 
00451    fMCMCH1NBins[index] = nbins;
00452 
00453    return;
00454 
00455 //    if(n<2) {
00456 //       BCLog::OutWarning("BCIntegrate::SetNbins. Number of bins less than 2 makes no sense. Setting to 2.");
00457 //       n=2;
00458 //    }
00459 //    MCMCSetH1NBins(n, -1);
00460 
00461    //   fNbins=n;
00462 
00463    //   fMCMCH1NBins = n;
00464    //   fMCMCH2NBinsX = n;
00465    //   fMCMCH2NBinsY = n;
00466 }
00467 
00468 // ---------------------------------------------------------
00469 // void BCIntegrate::SetNbinsX(int n)
00470 // {
00471 //    if(n<2) {
00472 //       BCLog::OutWarning("BCIntegrate::SetNbins. Number of bins less than 2 makes no sense. Setting to 2.");
00473 //       n=2;
00474 //    }
00475 //    fMCMCH2NBinsX = n;
00476 // }
00477 
00478 // ---------------------------------------------------------
00479 // void BCIntegrate::SetNbinsY(int n)
00480 // {
00481 //    if(n<2) {
00482 //       BCLog::OutWarning("BCIntegrate::SetNbins. Number of bins less than 2 makes no sense. Setting to 2.");
00483 //       n=2;
00484 //    }
00485 //    fNbins=n;
00486 
00487 //    fMCMCH2NBinsY = n;
00488 // }
00489 
00490 // ---------------------------------------------------------
00491 void BCIntegrate::SetVarList(int * varlist)
00492 {
00493    for(int i=0;i<fNvar;i++)
00494       fVarlist[i]=varlist[i];
00495 }
00496 
00497 // ---------------------------------------------------------
00498 void BCIntegrate::ResetVarlist(int v)
00499 {
00500    for(int i=0;i<fNvar;i++)
00501       fVarlist[i]=v;
00502 }
00503 
00504 // ---------------------------------------------------------
00505 double BCIntegrate::Eval(std::vector <double> x)
00506 {
00507    BCLog::OutWarning( "BCIntegrate::Eval. No function. Function needs to be overloaded.");
00508    return 0;
00509 }
00510 
00511 // ---------------------------------------------------------
00512 double BCIntegrate::LogEval(std::vector <double> x)
00513 {
00514    // this method should better also be overloaded
00515    return log(Eval(x));
00516 }
00517 
00518 // ---------------------------------------------------------
00519 double BCIntegrate::EvalSampling(std::vector <double> x)
00520 {
00521    BCLog::OutWarning( "BCIntegrate::EvalSampling. No function. Function needs to be overloaded.");
00522    return 0;
00523 }
00524 
00525 // ---------------------------------------------------------
00526 double BCIntegrate::LogEvalSampling(std::vector <double> x)
00527 {
00528    return log(EvalSampling(x));
00529 }
00530 
00531 // ---------------------------------------------------------
00532 double BCIntegrate::EvalPrint(std::vector <double> x)
00533 {
00534    double val=Eval(x);
00535    BCLog::OutDebug(Form("BCIntegrate::EvalPrint. Value: %g.", val));
00536 
00537    return val;
00538 }
00539 
00540 // ---------------------------------------------------------
00541 void BCIntegrate::SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method)
00542 {
00543 #ifdef HAVE_CUBA_H
00544    fIntegrationMethod = method;
00545 #else
00546    BCLog::OutWarning("!!! This version of BAT is compiled without Cuba.");
00547    BCLog::OutWarning("    Monte Carlo Sampled Mean integration method will be used.");
00548    BCLog::OutWarning("    To be able to use Cuba integration, install Cuba and recompile BAT.");
00549 #endif
00550 }
00551 
00552 // ---------------------------------------------------------
00553 int BCIntegrate::IntegrateResetResults()
00554 {
00555    fBestFitParameters.clear();
00556    fBestFitParameterErrors.clear();
00557    fBestFitParametersMarginalized.clear();
00558 
00559    // no errors
00560    return 1;
00561 }
00562 
00563 // ---------------------------------------------------------
00564 double BCIntegrate::Integrate()
00565 {
00566    std::vector <double> parameter;
00567    parameter.assign(fNvar, 0.0);
00568 
00569    BCLog::OutSummary(
00570          Form("Running numerical integration using %s (%s)",
00571          DumpIntegrationMethod().c_str(),
00572          DumpCubaIntegrationMethod().c_str()));
00573 
00574    switch(fIntegrationMethod)
00575    {
00576       case BCIntegrate::kIntMonteCarlo:
00577          return IntegralMC(parameter);
00578 
00579       case BCIntegrate::kIntMetropolis:
00580          return IntegralMetro(parameter);
00581 
00582       case BCIntegrate::kIntImportance:
00583          return IntegralImportance(parameter);
00584 
00585       case BCIntegrate::kIntCuba:
00586 #ifdef HAVE_CUBA_H
00587          return CubaIntegrate();
00588 #else
00589          BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
00590          BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
00591 #endif
00592       default:
00593          BCLog::OutError(
00594             Form("BCIntegrate::Integrate : Invalid integration method: %d", fIntegrationMethod));
00595    }
00596 
00597    return 0;
00598 }
00599 
00600 // ---------------------------------------------------------
00601 double BCIntegrate::IntegralMC(std::vector <double> x, int * varlist)
00602 {
00603    SetVarList(varlist);
00604    return IntegralMC(x);
00605 }
00606 
00607 // ---------------------------------------------------------
00608 double BCIntegrate::IntegralMC(std::vector <double> x)
00609 {
00610    // count the variables to integrate over
00611    int NvarNow=0;
00612 
00613    for(int i = 0; i < fNvar; i++)
00614       if(fVarlist[i])
00615          NvarNow++;
00616 
00617    // print to log
00618    BCLog::OutDebug(Form(" --> Running MC integation over %i dimensions.", NvarNow));
00619    BCLog::OutDebug(Form(" --> Maximum number of iterations: %i", GetNIterationsMax()));
00620    BCLog::OutDebug(Form(" --> Aimed relative precision:     %e", GetRelativePrecision()));
00621 
00622    // calculate D (the integrated volume)
00623    double D = 1.0;
00624    for(int i = 0; i < fNvar; i++)
00625       if (fVarlist[i])
00626          D = D * (fMax[i] - fMin[i]);
00627 
00628    // reset variables
00629    double pmax = 0.0;
00630    double sumW  = 0.0;
00631    double sumW2 = 0.0;
00632    double integral = 0.0;
00633    double variance = 0.0;
00634    double relprecision = 1.0;
00635 
00636    std::vector <double> randx;
00637    randx.assign(fNvar, 0.0);
00638 
00639    // reset number of iterations
00640    fNIterations = 0;
00641 
00642    // iterate while precision is not reached and the number of iterations is lower than maximum number of iterations
00643    while ((fRelativePrecision < relprecision && fNIterationsMax > fNIterations) || fNIterations < 10) {
00644       // increase number of iterations
00645       fNIterations++;
00646 
00647       // get random numbers
00648       GetRandomVector(randx);
00649 
00650       // scale random numbers
00651       for(int i = 0; i < fNvar; i++) {
00652          if(fVarlist[i])
00653             randx[i]=fMin[i]+randx[i]*(fMax[i]-fMin[i]);
00654          else
00655             randx[i]=x[i];
00656       }
00657 
00658       // evaluate function at sampled point
00659       double value = Eval(randx);
00660 
00661       // add value to sum and sum of squares
00662       sumW  += value;
00663       sumW2 += value * value;
00664 
00665       // search for maximum probability
00666       if (value > pmax) {
00667          // set new maximum value
00668          pmax = value;
00669 
00670          // delete old best fit parameter values
00671          fBestFitParameters.clear();
00672 
00673          // write best fit parameters
00674          for(int i = 0; i < fNvar; i++)
00675             fBestFitParameters.push_back(randx.at(i));
00676       }
00677 
00678       // calculate integral and variance
00679       integral = D * sumW / fNIterations;
00680 
00681       if (fNIterations%10000 == 0) {
00682          variance = (1.0 / double(fNIterations)) * (D * D * sumW2 / double(fNIterations) - integral * integral);
00683          double error = sqrt(variance);
00684          relprecision = error / integral;
00685 
00686          BCLog::OutDebug(Form("BCIntegrate::IntegralMC. Iteration %i, integral: %e +- %e.", fNIterations, integral, error));
00687       }
00688    }
00689 
00690    fError = variance / fNIterations;
00691 
00692    // print to log
00693    BCLog::OutDebug(Form(" --> Result of integration:        %e +- %e   in %i iterations.", integral, sqrt(variance), fNIterations));
00694    BCLog::OutDebug(Form(" --> Obtained relative precision:  %e. ", sqrt(variance) / integral));
00695 
00696    return integral;
00697 }
00698 
00699 
00700 // ---------------------------------------------------------
00701 double BCIntegrate::IntegralMetro(std::vector <double> x)
00702 {
00703    // print debug information
00704    BCLog::OutDebug(Form("BCIntegrate::IntegralMetro. Integate over %i dimensions.", fNvar));
00705 
00706    // get total number of iterations
00707    double Niter = pow(fNiterPerDimension, fNvar);
00708 
00709    // print if more than 100,000 iterations
00710    if(Niter>1e5)
00711       BCLog::OutDebug(Form(" --> Total number of iterations in Metropolis: %d.", (int)Niter));
00712 
00713    // reset sum
00714    double sumI = 0;
00715 
00716    // prepare Metropolis algorithm
00717    std::vector <double> randx;
00718    randx.assign(fNvar,0.);
00719    InitMetro();
00720 
00721    // prepare maximum value
00722    double pmax = 0.0;
00723 
00724    // loop over iterations
00725    for(int i = 0; i <= Niter; i++) {
00726       // get random point from sampling distribution
00727       GetRandomPointSamplingMetro(randx);
00728 
00729       // calculate probability at random point
00730       double val_f = Eval(randx);
00731 
00732       // calculate sampling distributions at that point
00733       double val_g = EvalSampling(randx);
00734 
00735       // add ratio to sum
00736       if (val_g > 0)
00737          sumI += val_f / val_g;
00738 
00739       // search for maximum probability
00740       if (val_f > pmax) {
00741          // set new maximum value
00742          pmax = val_f;
00743 
00744          // delete old best fit parameter values
00745          fBestFitParameters.clear();
00746 
00747          // write best fit parameters
00748          for(int i = 0; i < fNvar; i++)
00749             fBestFitParameters.push_back(randx.at(i));
00750       }
00751 
00752       // write intermediate results
00753       if((i+1)%100000 == 0)
00754          BCLog::OutDebug(Form("BCIntegrate::IntegralMetro. Iteration %d, integral: %g.", i+1, sumI/double(i+1)));
00755    }
00756 
00757    // calculate integral
00758    double result=sumI/Niter;
00759 
00760    // print debug information
00761    BCLog::OutDebug(Form(" --> Integral is %g in %g iterations.", result, Niter));
00762 
00763    return result;
00764 }
00765 
00766 // ---------------------------------------------------------
00767 double BCIntegrate::IntegralImportance(std::vector <double> x)
00768 {
00769    // print debug information
00770    BCLog::OutDebug(Form("BCIntegrate::IntegralImportance. Integate over %i dimensions.", fNvar));
00771 
00772    // get total number of iterations
00773    double Niter = pow(fNiterPerDimension, fNvar);
00774 
00775    // print if more than 100,000 iterations
00776    if(Niter>1e5)
00777       BCLog::OutDebug(Form("BCIntegrate::IntegralImportance. Total number of iterations: %d.", (int)Niter));
00778 
00779    // reset sum
00780    double sumI = 0;
00781 
00782    std::vector <double> randx;
00783    randx.assign(fNvar,0.);
00784 
00785    // prepare maximum value
00786    double pmax = 0.0;
00787 
00788    // loop over iterations
00789    for(int i = 0; i <= Niter; i++) {
00790       // get random point from sampling distribution
00791       GetRandomPointImportance(randx);
00792 
00793       // calculate probability at random point
00794       double val_f = Eval(randx);
00795 
00796       // calculate sampling distributions at that point
00797       double val_g = EvalSampling(randx);
00798 
00799       // add ratio to sum
00800       if (val_g > 0)
00801          sumI += val_f / val_g;
00802 
00803       // search for maximum probability
00804       if (val_f > pmax) {
00805          // set new maximum value
00806          pmax = val_f;
00807 
00808          // delete old best fit parameter values
00809          fBestFitParameters.clear();
00810 
00811          // write best fit parameters
00812          for(int i = 0; i < fNvar; i++)
00813             fBestFitParameters.push_back(randx.at(i));
00814       }
00815 
00816       // write intermediate results
00817       if((i+1)%100000 == 0)
00818          BCLog::OutDebug(Form("BCIntegrate::IntegralImportance : Iteration %d, integral: %g.", i+1, sumI/double(i+1)));
00819    }
00820 
00821    // calculate integral
00822    double result=sumI/Niter;
00823 
00824    // print debug information
00825    BCLog::OutDebug(Form("BCIntegrate::IntegralImportance : Integral %g in %i iterations.", result, (int)Niter));
00826 
00827    return result;
00828 }
00829 
00830 // ---------------------------------------------------------
00831 TH1D* BCIntegrate::Marginalize(BCParameter * parameter)
00832 {
00833    BCLog::OutDebug(Form(" --> Marginalizing model wrt. parameter %s using %s.", parameter->GetName().data(), DumpMarginalizationMethod().c_str()));
00834 
00835    switch(fMarginalizationMethod)
00836    {
00837       case BCIntegrate::kMargMonteCarlo:
00838          return MarginalizeByIntegrate(parameter);
00839 
00840       case BCIntegrate::kMargMetropolis:
00841          return MarginalizeByMetro(parameter);
00842 
00843       default:
00844          BCLog::OutError(
00845             Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod));
00846 
00847    }
00848 
00849    return 0;
00850 }
00851 
00852 // ---------------------------------------------------------
00853 TH2D * BCIntegrate::Marginalize(BCParameter * parameter1, BCParameter * parameter2)
00854 {
00855    switch(fMarginalizationMethod)
00856    {
00857       case BCIntegrate::kMargMonteCarlo:
00858          return MarginalizeByIntegrate(parameter1,parameter2);
00859 
00860       case BCIntegrate::kMargMetropolis:
00861          return MarginalizeByMetro(parameter1,parameter2);
00862 
00863       default:
00864          BCLog::OutError(
00865             Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod));
00866    }
00867 
00868    return 0;
00869 }
00870 
00871 // ---------------------------------------------------------
00872 TH1D* BCIntegrate::MarginalizeByIntegrate(BCParameter * parameter)
00873 {
00874    // set parameter to marginalize
00875    ResetVarlist(1);
00876    int index = parameter->GetIndex();
00877    UnsetVar(index);
00878 
00879    // define histogram
00880    double hmin = parameter->GetLowerLimit();
00881    double hmax = parameter->GetUpperLimit();
00882    double hdx  = (hmax - hmin) / double(fNbins);
00883    TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax);
00884 
00885    // integrate
00886    std::vector <double> randx;
00887    randx.assign(fNvar, 0.);
00888 
00889    for(int i=0;i<=fNbins;i++) {
00890       randx[index] = hmin + (double)i * hdx;
00891 
00892       double val = IntegralMC(randx);
00893       hist->Fill(randx[index], val);
00894    }
00895 
00896    // normalize
00897    hist->Scale( 1./hist->Integral("width") );
00898 
00899    return hist;
00900 }
00901 
00902 // ---------------------------------------------------------
00903 TH2D * BCIntegrate::MarginalizeByIntegrate(BCParameter * parameter1, BCParameter * parameter2)
00904 {
00905    // set parameter to marginalize
00906    ResetVarlist(1);
00907    int index1 = parameter1->GetIndex();
00908    UnsetVar(index1);
00909    int index2 = parameter2->GetIndex();
00910    UnsetVar(index2);
00911 
00912    // define histogram
00913    double hmin1 = parameter1->GetLowerLimit();
00914    double hmax1 = parameter1->GetUpperLimit();
00915    double hdx1  = (hmax1 - hmin1) / double(fNbins);
00916 
00917    double hmin2 = parameter2->GetLowerLimit();
00918    double hmax2 = parameter2->GetUpperLimit();
00919    double hdx2  = (hmax2 - hmin2) / double(fNbins);
00920 
00921    TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1->GetName().data(), parameter2->GetName().data()),"",
00922          fNbins, hmin1, hmax1,
00923          fNbins, hmin2, hmax2);
00924 
00925    // integrate
00926    std::vector <double> randx;
00927    randx.assign(fNvar, 0.0);
00928 
00929    for(int i=0;i<=fNbins;i++) {
00930       randx[index1] = hmin1 + (double)i * hdx1;
00931       for(int j=0;j<=fNbins;j++) {
00932          randx[index2] = hmin2 + (double)j * hdx2;
00933 
00934          double val = IntegralMC(randx);
00935          hist->Fill(randx[index1],randx[index2], val);
00936       }
00937    }
00938 
00939    // normalize
00940    hist->Scale(1.0/hist->Integral("width"));
00941 
00942    return hist;
00943 }
00944 
00945 // ---------------------------------------------------------
00946 TH1D * BCIntegrate::MarginalizeByMetro(BCParameter * parameter)
00947 {
00948    int niter = fMarkovChainNIterations;
00949 
00950    if (fMarkovChainAutoN == true)
00951       niter = fNbins*fNbins*fNSamplesPer2DBin*fNvar;
00952 
00953    BCLog::OutDetail(Form(" --> Number of samples in Metropolis marginalization: %d.", niter));
00954 
00955    // set parameter to marginalize
00956    int index = parameter->GetIndex();
00957 
00958    // define histogram
00959    double hmin = parameter->GetLowerLimit();
00960    double hmax = parameter->GetUpperLimit();
00961    TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax);
00962 
00963    // prepare Metro
00964    std::vector <double> randx;
00965    randx.assign(fNvar, 0.0);
00966    InitMetro();
00967 
00968    for(int i=0;i<=niter;i++) {
00969       GetRandomPointMetro(randx);
00970       hist->Fill(randx[index]);
00971    }
00972 
00973    // normalize
00974    hist->Scale(1.0/hist->Integral("width"));
00975 
00976    return hist;
00977 }
00978 
00979 // ---------------------------------------------------------
00980 TH2D * BCIntegrate::MarginalizeByMetro(BCParameter * parameter1, BCParameter * parameter2)
00981 {
00982    int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar;
00983 
00984    // set parameter to marginalize
00985    int index1 = parameter1->GetIndex();
00986    int index2 = parameter2->GetIndex();
00987 
00988    // define histogram
00989    double hmin1 = parameter1->GetLowerLimit();
00990    double hmax1 = parameter1->GetUpperLimit();
00991 
00992    double hmin2 = parameter2->GetLowerLimit();
00993    double hmax2 = parameter2->GetUpperLimit();
00994 
00995    TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1->GetName().data(), parameter2->GetName().data()),"",
00996          fNbins, hmin1, hmax1,
00997          fNbins, hmin2, hmax2);
00998 
00999    // prepare Metro
01000    std::vector <double> randx;
01001    randx.assign(fNvar, 0.0);
01002    InitMetro();
01003 
01004    for(int i=0;i<=niter;i++) {
01005       GetRandomPointMetro(randx);
01006       hist->Fill(randx[index1],randx[index2]);
01007    }
01008 
01009    // normalize
01010    hist->Scale(1.0/hist->Integral("width"));
01011 
01012    return hist;
01013 }
01014 
01015 // ---------------------------------------------------------
01016 int BCIntegrate::MarginalizeAllByMetro(const char * name="")
01017 {
01018    int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar;
01019 
01020    BCLog::OutDetail(Form(" --> Number of samples in Metropolis marginalization: %d.", niter));
01021 
01022    // define 1D histograms
01023    for(int i=0;i<fNvar;i++)
01024    {
01025       double hmin1 = fx->at(i)->GetLowerLimit();
01026       double hmax1 = fx->at(i)->GetUpperLimit();
01027 
01028       TH1D * h1 = new TH1D(
01029             TString::Format("h%s_%s", name, fx->at(i)->GetName().data()),"",
01030             fNbins, hmin1, hmax1);
01031 
01032       fHProb1D.push_back(h1);
01033    }
01034 
01035    // define 2D histograms
01036    for(int i=0;i<fNvar-1;i++)
01037       for(int j=i+1;j<fNvar;j++) {
01038          double hmin1 = fx->at(i)->GetLowerLimit();
01039          double hmax1 = fx->at(i)->GetUpperLimit();
01040 
01041          double hmin2 = fx->at(j)->GetLowerLimit();
01042          double hmax2 = fx->at(j)->GetUpperLimit();
01043 
01044          TH2D * h2 = new TH2D(
01045             TString::Format("h%s_%s_%s", name, fx->at(i)->GetName().data(), fx->at(j)->GetName().data()),"",
01046             fNbins, hmin1, hmax1,
01047             fNbins, hmin2, hmax2);
01048 
01049          fHProb2D.push_back(h2);
01050       }
01051 
01052    // get number of 2d distributions
01053    int nh2d = fHProb2D.size();
01054 
01055    BCLog::OutDetail(Form(" --> Marginalizing %d 1D distributions and %d 2D distributions.", fNvar, nh2d));
01056 
01057    // prepare function fitting
01058    double dx = 0.;
01059    double dy = 0.;
01060 
01061    if (fFitFunctionIndexX >= 0) {
01062       dx = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX) -
01063             fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX))
01064             / double(fErrorBandNbinsX);
01065 
01066       dx = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY) -
01067             fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY))
01068             / double(fErrorBandNbinsY);
01069 
01070       fErrorBandXY = new TH2D(
01071             TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "",
01072             fErrorBandNbinsX,
01073             fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX) - 0.5 * dx,
01074             fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX) + 0.5 * dx,
01075             fErrorBandNbinsY,
01076             fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY) - 0.5 * dy,
01077             fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY) + 0.5 * dy);
01078       fErrorBandXY->SetStats(kFALSE);
01079 
01080       for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
01081          for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
01082             fErrorBandXY->SetBinContent(ix, iy, 0.0);
01083    }
01084 
01085    // prepare Metro
01086    std::vector <double> randx;
01087 
01088    randx.assign(fNvar, 0.0);
01089    InitMetro();
01090 
01091    // reset counter
01092    fNacceptedMCMC = 0;
01093 
01094    // run Metro and fill histograms
01095    for(int i=0;i<=niter;i++) {
01096       GetRandomPointMetro(randx);
01097 
01098       // save this point to the markov chain in the ROOT file
01099       if (fFlagWriteMarkovChain) {
01100          fMCMCIteration = i;
01101          fMarkovChainTree->Fill();
01102       }
01103 
01104       for(int j=0;j<fNvar;j++)
01105          fHProb1D[j]->Fill( randx[j] );
01106 
01107       int ih=0;
01108       for(int j=0;j<fNvar-1;j++)
01109          for(int k=j+1;k<fNvar;k++) {
01110             fHProb2D[ih]->Fill(randx[j],randx[k]);
01111             ih++;
01112          }
01113 
01114       if((i+1)%100000==0)
01115          BCLog::OutDebug(Form("BCIntegrate::MarginalizeAllByMetro. %d samples done.",i+1));
01116 
01117       // function fitting
01118 
01119       if (fFitFunctionIndexX >= 0) {
01120          // loop over all possible x values ...
01121          if (fErrorBandContinuous) {
01122             double x = 0;
01123 
01124             for (int ix = 0; ix < fErrorBandNbinsX; ix++) {
01125                // calculate x
01126                x = fErrorBandXY->GetXaxis()->GetBinCenter(ix + 1);
01127 
01128                // calculate y
01129                std::vector <double> xvec;
01130                xvec.push_back(x);
01131                double y = FitFunction(xvec, randx);
01132                xvec.clear();
01133 
01134                // fill histogram
01135                fErrorBandXY->Fill(x, y);
01136             }
01137          }
01138 
01139          // ... or evaluate at the data point x-values
01140          else {
01141             int ndatapoints = int(fErrorBandX.size());
01142             double x = 0;
01143 
01144             for (int ix = 0; ix < ndatapoints; ++ix) {
01145                // calculate x
01146                x = fErrorBandX.at(ix);
01147 
01148                // calculate y
01149                std::vector <double> xvec;
01150                xvec.push_back(x);
01151                double y = FitFunction(xvec, randx);
01152                xvec.clear();
01153 
01154                // fill histogram
01155                fErrorBandXY->Fill(x, y);
01156             }
01157          }
01158       }
01159    }
01160 
01161    // normalize histograms
01162    for(int i=0;i<fNvar;i++)
01163       fHProb1D[i]->Scale( 1./fHProb1D[i]->Integral("width") );
01164 
01165    for (int i=0;i<nh2d;i++)
01166       fHProb2D[i]->Scale( 1./fHProb2D[i]->Integral("width") );
01167 
01168    if (fFitFunctionIndexX >= 0)
01169       fErrorBandXY->Scale(1.0/fErrorBandXY->Integral() * fErrorBandXY->GetNbinsX());
01170 
01171    if (fFitFunctionIndexX >= 0) {
01172       for (int ix = 1; ix <= fErrorBandNbinsX; ix++) {
01173          double sum = 0;
01174 
01175          for (int iy = 1; iy <= fErrorBandNbinsY; iy++)
01176             sum += fErrorBandXY->GetBinContent(ix, iy);
01177 
01178          for (int iy = 1; iy <= fErrorBandNbinsY; iy++) {
01179             double newvalue = fErrorBandXY->GetBinContent(ix, iy) / sum;
01180             fErrorBandXY->SetBinContent(ix, iy, newvalue);
01181          }
01182       }
01183    }
01184 
01185    BCLog::OutDebug(Form("BCIntegrate::MarginalizeAllByMetro done with %i trials and %i accepted trials. Efficiency is %f",fNmetro, fNacceptedMCMC, double(fNacceptedMCMC)/double(fNmetro)));
01186 
01187    return fNvar+nh2d;
01188 }
01189 
01190 // ---------------------------------------------------------
01191 TH1D * BCIntegrate::GetH1D(int parIndex)
01192 {
01193    if(fHProb1D.size()==0) {
01194       BCLog::OutWarning("BCModel::GetH1D. MarginalizeAll() has to be run prior to this to fill the distributions.");
01195       return 0;
01196    }
01197 
01198    if(parIndex<0 || parIndex>fNvar) {
01199       BCLog::OutWarning(Form("BCIntegrate::GetH1D. Parameter index %d is invalid.",parIndex));
01200       return 0;
01201    }
01202 
01203    return fHProb1D[parIndex];
01204 }
01205 
01206 // ---------------------------------------------------------
01207 int BCIntegrate::GetH2DIndex(int parIndex1, int parIndex2)
01208 {
01209    if(parIndex1>fNvar-1 || parIndex1<0) {
01210       BCLog::OutWarning(Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex1));
01211       return -1;
01212    }
01213 
01214    if(parIndex2>fNvar-1 || parIndex2<0) {
01215       BCLog::OutWarning(Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex2));
01216       return -1;
01217    }
01218 
01219    if(parIndex1==parIndex2) {
01220       BCLog::OutWarning(Form("BCIntegrate::GetH2DIndex. Parameters have equal indeces: %d , %d", parIndex1, parIndex2));
01221       return -1;
01222    }
01223 
01224    if(parIndex1>parIndex2) {
01225       BCLog::OutWarning("BCIntegrate::GetH2DIndex. First parameters must be smaller than second (sorry).");
01226       return -1;
01227    }
01228 
01229    int index=0;
01230    for(int i=0;i<fNvar-1;i++)
01231       for(int j=i+1;j<fNvar;j++) {
01232          if(i==parIndex1 && j==parIndex2)
01233             return index;
01234          index++;
01235       }
01236 
01237    BCLog::OutWarning("BCIntegrate::GetH2DIndex. Invalid index combination.");
01238 
01239    return -1;
01240 }
01241 
01242 // ---------------------------------------------------------
01243 TH2D * BCIntegrate::GetH2D(int parIndex1, int parIndex2)
01244 {
01245    if(fHProb2D.size()==0) {
01246       BCLog::OutWarning("BCModel::GetH2D. MarginalizeAll() has to be run prior to this to fill the distributions.");
01247       return 0;
01248    }
01249 
01250    int hindex = GetH2DIndex(parIndex1, parIndex2);
01251    if(hindex==-1)
01252       return 0;
01253 
01254    if(hindex>(int)fHProb2D.size()-1) {
01255       BCLog::OutWarning("BCIntegrate::GetH2D. Got invalid index from GetH2DIndex(). Something went wrong.");
01256       return 0;
01257    }
01258 
01259    return fHProb2D[hindex];
01260 }
01261 
01262 // ---------------------------------------------------------
01263 double BCIntegrate::GetRandomPoint(std::vector <double> &x)
01264 {
01265    GetRandomVector(x);
01266 
01267    for(int i=0;i<fNvar;i++)
01268       x[i]=fMin[i]+x[i]*(fMax[i]-fMin[i]);
01269 
01270    return Eval(x);
01271 }
01272 
01273 // ---------------------------------------------------------
01274 double BCIntegrate::GetRandomPointImportance(std::vector <double> &x)
01275 {
01276    double p = 1.1;
01277    double g = 1.0;
01278 
01279    while (p>g) {
01280       GetRandomVector(x);
01281 
01282       for(int i=0;i<fNvar;i++)
01283          x[i] = fMin[i] + x[i] * (fMax[i]-fMin[i]);
01284 
01285       p = fRandom->Rndm();
01286 
01287       g = EvalSampling(x);
01288    }
01289 
01290    return Eval(x);
01291 }
01292 
01293 // ---------------------------------------------------------
01294 void BCIntegrate::InitMetro()
01295 {
01296    fNmetro=0;
01297 
01298    if (fXmetro0.size() <= 0) {
01299       // start in the center of the phase space
01300       for(int i=0;i<fNvar;i++)
01301          fXmetro0.push_back((fMin[i]+fMax[i])/2.0);
01302    }
01303 
01304    fMarkovChainValue =  LogEval(fXmetro0);
01305 
01306    // run metropolis for a few times and dump the result... (to forget the initial position)
01307    std::vector <double> x;
01308    x.assign(fNvar,0.);
01309 
01310    for(int i=0;i<1000;i++)
01311       GetRandomPointMetro(x);
01312 
01313    fNmetro = 0;
01314 }
01315 
01316 // ---------------------------------------------------------
01317 void BCIntegrate::GetRandomPointMetro(std::vector <double> &x)
01318 {
01319    // get new point
01320    GetRandomVectorMetro(fXmetro1);
01321 
01322    // scale the point to the allowed region and stepsize
01323    int in=1;
01324    for(int i=0;i<fNvar;i++) {
01325       fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]);
01326 
01327       // check whether the generated point is inside the allowed region
01328       if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] )
01329          in=0;
01330    }
01331 
01332    // calculate the log probabilities and compare old and new point
01333 
01334    double p0 = fMarkovChainValue; // old point
01335    double p1 = 0; // new point
01336    int accept=0;
01337 
01338    if(in) {
01339       p1 = LogEval(fXmetro1);
01340 
01341       if(p1>=p0)
01342          accept=1;
01343       else {
01344          double r=log(fRandom->Rndm());
01345          if(r<p1-p0)
01346             accept=1;
01347       }
01348    }
01349 
01350    // fill the return point after the decision
01351    if(accept) {
01352       // increase counter
01353       fNacceptedMCMC++;
01354 
01355       for(int i=0;i<fNvar;i++) {
01356          fXmetro0[i]=fXmetro1[i];
01357          x[i]=fXmetro1[i];
01358          fMarkovChainValue = p1;
01359       }
01360    }
01361    else
01362       for(int i=0;i<fNvar;i++) {
01363          x[i]=fXmetro0[i];
01364          fXmetro1[i] = x[i];
01365          fMarkovChainValue = p0;
01366       }
01367 
01368    fNmetro++;
01369 }
01370 
01371 // ---------------------------------------------------------
01372 void BCIntegrate::GetRandomPointSamplingMetro(std::vector <double> &x)
01373 {
01374    // get new point
01375    GetRandomVectorMetro(fXmetro1);
01376 
01377    // scale the point to the allowed region and stepsize
01378    int in=1;
01379    for(int i=0;i<fNvar;i++) {
01380       fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]);
01381 
01382       // check whether the generated point is inside the allowed region
01383       if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] )
01384          in=0;
01385    }
01386 
01387    // calculate the log probabilities and compare old and new point
01388    double p0 = LogEvalSampling(fXmetro0); // old point
01389    double p1=0; // new point
01390    if(in)
01391       p1= LogEvalSampling(fXmetro1);
01392 
01393    // compare log probabilities
01394    int accept=0;
01395    if(in) {
01396       if(p1>=p0)
01397          accept=1;
01398       else {
01399          double r=log(fRandom->Rndm());
01400          if(r<p1-p0)
01401             accept=1;
01402       }
01403    }
01404 
01405    // fill the return point after the decision
01406    if(accept)
01407       for(int i=0;i<fNvar;i++) {
01408          fXmetro0[i]=fXmetro1[i];
01409          x[i]=fXmetro1[i];
01410       }
01411    else
01412       for(int i=0;i<fNvar;i++)
01413          x[i]=fXmetro0[i];
01414 
01415    fNmetro++;
01416 }
01417 
01418 // ---------------------------------------------------------
01419 void BCIntegrate::GetRandomVector(std::vector <double> &x)
01420 {
01421    double * randx = new double[fNvar];
01422 
01423    fRandom->RndmArray(fNvar, randx);
01424 
01425    for(int i=0;i<fNvar;i++)
01426       x[i] = randx[i];
01427 
01428    delete[] randx;
01429    randx = 0;
01430 }
01431 
01432 // ---------------------------------------------------------
01433 void BCIntegrate::GetRandomVectorMetro(std::vector <double> &x)
01434 {
01435    double * randx = new double[fNvar];
01436 
01437    fRandom->RndmArray(fNvar, randx);
01438 
01439    for(int i=0;i<fNvar;i++)
01440       x[i] = 2.0 * (randx[i] - 0.5) * fMarkovChainStepSize;
01441 
01442    delete[] randx;
01443    randx = 0;
01444 }
01445 
01446 // ---------------------------------------------------------
01447 TMinuit * BCIntegrate::GetMinuit()
01448 {
01449    if (!fMinuit)
01450       fMinuit = new TMinuit();
01451 
01452    return fMinuit;
01453 }
01454 
01455 // ---------------------------------------------------------
01456 
01457 void BCIntegrate::FindModeMinuit(std::vector<double> start, int printlevel)
01458 {
01459    bool have_start = true;
01460 
01461    // check start values
01462    if (int(start.size()) != fNvar)
01463       have_start = false;
01464 
01465    // set global this
01466    global_this = this;
01467 
01468    // define minuit
01469    if (fMinuit)
01470       delete fMinuit;
01471    fMinuit = new TMinuit(fNvar);
01472 
01473    // set print level
01474    fMinuit->SetPrintLevel(printlevel);
01475 
01476    // set function
01477    fMinuit->SetFCN(&BCIntegrate::FCNLikelihood);
01478 
01479    // set UP for likelihood
01480    fMinuit->SetErrorDef(0.5);
01481 
01482    // set parameters
01483    int flag;
01484    for (int i = 0; i < fNvar; i++) {
01485       double starting_point = (fMin[i]+fMax[i])/2.;
01486       if(have_start)
01487          starting_point = start[i];
01488          fMinuit->mnparm(i,
01489                          fx->at(i)->GetName().data(),
01490                          starting_point,
01491                          (fMax[i]-fMin[i])/100.,
01492                          fMin[i],
01493                          fMax[i],
01494                          flag);
01495    }
01496 
01497    // do mcmc minimization
01498    //   fMinuit->mnseek();
01499 
01500    // do minimization
01501    fMinuit->mnexcm("MIGRAD", fMinuitArglist, 2, flag);
01502 
01503    // improve search for local minimum
01504    //   fMinuit->mnimpr();
01505 
01506    // copy flag
01507    fMinuitErrorFlag = flag;
01508 
01509    // check if mode found by minuit is better than previous estimation
01510    double probmax = 0;
01511    bool valid = false;
01512 
01513    if ( int(fBestFitParameters.size()) == fNvar) {
01514       valid = true;
01515       for (int i = 0; i < fNvar; ++i)
01516          if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i])
01517             valid= false;
01518 
01519       if (valid)
01520          probmax = Eval( fBestFitParameters );
01521    }
01522 
01523    std::vector<double> tempvec;
01524    for (int i = 0; i < fNvar; i++) {
01525       double par;
01526       double parerr;
01527       fMinuit->GetParameter(i, par, parerr);
01528       tempvec.push_back(par);
01529    }
01530    double probmaxminuit = Eval( tempvec );
01531 
01532    // set best fit parameters
01533    if ( (probmaxminuit > probmax && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization) {
01534       fBestFitParameters.clear();
01535       fBestFitParameterErrors.clear();
01536 
01537       for (int i = 0; i < fNvar; i++) {
01538          double par;
01539          double parerr;
01540          fMinuit->GetParameter(i, par, parerr);
01541          fBestFitParameters.push_back(par);
01542          fBestFitParameterErrors.push_back(parerr);
01543       }
01544 
01545       // set optimization method used to find the mode
01546       fOptimizationMethodMode = BCIntegrate::kOptMinuit;
01547    }
01548 
01549    // delete minuit
01550    delete fMinuit;
01551    fMinuit = 0;
01552 
01553    return;
01554 }
01555 
01556 // ---------------------------------------------------------
01557 void BCIntegrate::InitializeSATree()
01558 {
01559   if (fTreeSA)
01560     delete fTreeSA;
01561   fTreeSA = new TTree("SATree", "SATree");
01562 
01563    fTreeSA->Branch("Iteration",      &fSANIterations,   "iteration/I");
01564    fTreeSA->Branch("NParameters",    &fNvar,            "parameters/I");
01565    fTreeSA->Branch("Temperature",    &fSATemperature,   "temperature/D");
01566    fTreeSA->Branch("LogProbability", &fSALogProb,       "log(probability)/D");
01567 
01568    for (int i = 0; i < fNvar; ++i)
01569       fTreeSA->Branch(TString::Format("Parameter%i", i), &fSAx[i], TString::Format("parameter %i/D", i));
01570 }
01571 
01572 // ---------------------------------------------------------
01573 void BCIntegrate::FindModeSA(std::vector<double> start)
01574 {
01575    // note: if f(x) is the function to be minimized, then
01576    // f(x) := - LogEval(parameters)
01577 
01578    bool have_start = true;
01579    std::vector<double> x, y, best_fit; // vectors for current state, new proposed state and best fit up to now
01580    double fval_x, fval_y, fval_best_fit; // function values at points x, y and best_fit (we save them rather than to re-calculate them every time)
01581    int t = 1; // time iterator
01582 
01583    // check start values
01584    if (int(start.size()) != fNvar)
01585       have_start = false;
01586 
01587    // if no starting point is given, set to center of parameter space
01588    if ( !have_start ) {
01589       start.clear();
01590       for (int i = 0; i < fNvar; i++)
01591          start.push_back((fMin[i]+fMax[i])/2.);
01592    }
01593 
01594    // set current state and best fit to starting point
01595    x.clear();
01596    best_fit.clear();
01597    for (int i = 0; i < fNvar; i++) {
01598       x.push_back(start[i]);
01599       best_fit.push_back(start[i]);
01600    }
01601    // calculate function value at starting point
01602    fval_x = fval_best_fit = LogEval(x);
01603 
01604    // run while still "hot enough"
01605    while ( SATemperature(t) > fSATmin ) {
01606       // generate new state
01607       y = GetProposalPointSA(x, t);
01608 
01609       // check if the proposed point is inside the phase space
01610       // if not, reject it
01611       bool is_in_ranges = true;
01612 
01613       for (int i = 0; i < fNvar; i++)
01614          if (y[i] > fMax[i] || y[i] < fMin[i])
01615             is_in_ranges = false;
01616 
01617       if ( !is_in_ranges )
01618          ; // do nothing...
01619       else {
01620          // calculate function value at new point
01621          fval_y = LogEval(y);
01622 
01623          // is it better than the last one?
01624          // if so, update state and chef if it is the new best fit...
01625          if (fval_y >= fval_x) {
01626             x.clear();
01627             for (int i = 0; i < fNvar; i++)
01628                x.push_back(y[i]);
01629 
01630             fval_x = fval_y;
01631 
01632             if (fval_y > fval_best_fit) {
01633                best_fit.clear();
01634                for (int i = 0; i < fNvar; i++)
01635                   best_fit.push_back(y[i]);
01636 
01637                fval_best_fit = fval_y;
01638             }
01639          }
01640          // ...else, only accept new state w/ certain probability
01641          else {
01642             if (fRandom->Rndm() <= exp( (fval_y - fval_x) / SATemperature(t) )) {
01643                x.clear();
01644                for (int i = 0; i < fNvar; i++)
01645                   x.push_back(y[i]);
01646 
01647                fval_x = fval_y;
01648             }
01649          }
01650       }
01651 
01652       // update tree variables
01653       fSANIterations = t;
01654       fSATemperature = SATemperature(t);
01655       fSALogProb = fval_x;
01656       fSAx.clear();
01657       for (int i = 0; i < fNvar; ++i)
01658          fSAx.push_back(x[i]);
01659 
01660       // fill tree
01661       if (fFlagWriteSAToFile)
01662          fTreeSA->Fill();
01663 
01664       // increate t
01665       t++;
01666    }
01667 
01668    // check if mode found by minuit is better than previous estimation
01669    double probmax = 0;
01670    bool valid = false;
01671 
01672    if ( int(fBestFitParameters.size()) == fNvar) {
01673       valid = true;
01674       for (int i = 0; i < fNvar; ++i)
01675          if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i])
01676             valid= false;
01677 
01678       if (valid)
01679          probmax = Eval( fBestFitParameters );
01680    }
01681 
01682    double probmaxsa = Eval( best_fit);
01683 
01684    if ( (probmaxsa > probmax  && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization) {
01685       // set best fit parameters
01686       fBestFitParameters.clear();
01687       fBestFitParameterErrors.clear();
01688 
01689       for (int i = 0; i < fNvar; i++) {
01690          fBestFitParameters.push_back(best_fit[i]);
01691          fBestFitParameterErrors.push_back(-1.); // error undefined
01692       }
01693 
01694       // set optimization moethod used to find the mode
01695       fOptimizationMethodMode = BCIntegrate::kOptSA;
01696    }
01697 
01698    return;
01699 }
01700 
01701 // ---------------------------------------------------------
01702 
01703 double BCIntegrate::SATemperature(double t)
01704 {
01705    // do we have Cauchy (default), Boltzmann or custom annealing schedule?
01706    if (fSASchedule == BCIntegrate::kSABoltzmann)
01707       return SATemperatureBoltzmann(t);
01708    else if (fSASchedule == BCIntegrate::kSACauchy)
01709       return SATemperatureCauchy(t);
01710    else
01711       return SATemperatureCustom(t);
01712 }
01713 
01714 // ---------------------------------------------------------
01715 
01716 double BCIntegrate::SATemperatureBoltzmann(double t)
01717 {
01718    return fSAT0 / log((double)(t + 1));
01719 }
01720 
01721 // ---------------------------------------------------------
01722 
01723 double BCIntegrate::SATemperatureCauchy(double t)
01724 {
01725    return fSAT0 / (double)t;
01726 }
01727 
01728 // ---------------------------------------------------------
01729 
01730 double BCIntegrate::SATemperatureCustom(double t)
01731 {
01732    BCLog::OutError("BCIntegrate::SATemperatureCustom : No custom temperature schedule defined");
01733    return 0.;
01734 }
01735 
01736 // ---------------------------------------------------------
01737 
01738 std::vector<double> BCIntegrate::GetProposalPointSA(std::vector<double> x, int t)
01739 {
01740    // do we have Cauchy (default), Boltzmann or custom annealing schedule?
01741    if (fSASchedule == BCIntegrate::kSABoltzmann)
01742       return GetProposalPointSABoltzmann(x, t);
01743    else if (fSASchedule == BCIntegrate::kSACauchy)
01744       return GetProposalPointSACauchy(x, t);
01745    else
01746       return GetProposalPointSACustom(x, t);
01747 }
01748 
01749 // ---------------------------------------------------------
01750 
01751 std::vector<double> BCIntegrate::GetProposalPointSABoltzmann(std::vector<double> x, int t)
01752 {
01753    std::vector<double> y;
01754    y.clear();
01755    double new_val, norm;
01756 
01757    for (int i = 0; i < fNvar; i++) {
01758       norm = (fMax[i] - fMin[i]) * SATemperature(t) / 2.;
01759       new_val = x[i] + norm * fRandom->Gaus();
01760       y.push_back(new_val);
01761    }
01762 
01763    return y;
01764 }
01765 
01766 // ---------------------------------------------------------
01767 
01768 std::vector<double> BCIntegrate::GetProposalPointSACauchy(std::vector<double> x, int t)
01769 {
01770    std::vector<double> y;
01771    y.clear();
01772 
01773    if (fNvar == 1) {
01774       double cauchy, new_val, norm;
01775 
01776       norm = (fMax[0] - fMin[0]) * SATemperature(t) / 2.;
01777       cauchy = tan(3.14159 * (fRandom->Rndm() - 0.5));
01778       new_val = x[0] + norm * cauchy;
01779       y.push_back(new_val);
01780    }
01781    else {
01782       // use sampling to get radial n-dim Cauchy distribution
01783 
01784       // first generate a random point uniformly distributed on a
01785       // fNvar-dimensional hypersphere
01786       y = SAHelperGetRandomPointOnHypersphere();
01787 
01788       // scale the vector by a random factor determined by the radial
01789       // part of the fNvar-dimensional Cauchy distribution
01790       double radial = SATemperature(t) * SAHelperGetRadialCauchy();
01791 
01792       // scale y by radial part and the size of dimension i in phase space
01793       // afterwards, move by x
01794       for (int i = 0; i < fNvar; i++)
01795          y[i] = (fMax[i] - fMin[i]) * y[i] * radial / 2. + x[i];
01796    }
01797 
01798    return y;
01799 }
01800 
01801 // ---------------------------------------------------------
01802 
01803 std::vector<double> BCIntegrate::GetProposalPointSACustom(std::vector<double> x, int t)
01804 {
01805    BCLog::OutError("BCIntegrate::GetProposalPointSACustom : No custom proposal function defined");
01806    return std::vector<double>(fNvar);
01807 }
01808 
01809 // ---------------------------------------------------------
01810 
01811 std::vector<double> BCIntegrate::SAHelperGetRandomPointOnHypersphere()
01812 {
01813    std::vector<double> rand_point(fNvar);
01814 
01815    // This method can only be called with fNvar >= 2 since the 1-dim case
01816    // is already hard wired into the Cauchy annealing proposal function.
01817    // To speed things up, hard-code fast method for 2 and dimensions.
01818    // The algorithm for 2D can be found at
01819    // http://mathworld.wolfram.com/CirclePointPicking.html
01820    // For 3D just using ROOT's algorithm.
01821    if (fNvar == 2) {
01822       double x1, x2, s;
01823       do {
01824          x1 = fRandom->Rndm() * 2. - 1.;
01825          x2 = fRandom->Rndm() * 2. - 1.;
01826          s = x1*x1 + x2*x2;
01827       }
01828       while (s >= 1);
01829 
01830       rand_point[0] = (x1*x1 - x2*x2) / s;
01831       rand_point[1] = (2.*x1*x2) / s;
01832    }
01833    else if (fNvar == 3) {
01834       fRandom->Sphere(rand_point[0], rand_point[1], rand_point[2], 1.0);
01835    }
01836    else {
01837       double s = 0.,
01838          gauss_num;
01839 
01840       for (int i = 0; i < fNvar; i++) {
01841          gauss_num = fRandom->Gaus();
01842          rand_point[i] = gauss_num;
01843          s += gauss_num * gauss_num;
01844       }
01845       s = sqrt(s);
01846 
01847       for (int i = 0; i < fNvar; i++)
01848          rand_point[i] = rand_point[i] / s;
01849    }
01850 
01851    return rand_point;
01852 }
01853 
01854 // ---------------------------------------------------------
01855 
01856 double BCIntegrate::SAHelperGetRadialCauchy()
01857 {
01858    // theta is sampled from a rather complicated distribution,
01859    // so first we create a lookup table with 10000 random numbers
01860    // once and then, each time we need a new random number,
01861    // we just look it up in the table.
01862    double theta;
01863 
01864    // static vectors for theta-sampling-map
01865    static double map_u[10001];
01866    static double map_theta[10001];
01867    static bool initialized = false;
01868    static int map_dimension = 0;
01869 
01870    // is the lookup-table already initialized? if not, do it!
01871    if (!initialized || map_dimension != fNvar) {
01872       double init_theta;
01873       double beta = SAHelperSinusToNIntegral(fNvar - 1, 1.57079632679);
01874 
01875       for (int i = 0; i <= 10000; i++) {
01876          init_theta = 3.14159265 * (double)i / 5000.;
01877          map_theta[i] = init_theta;
01878 
01879          map_u[i] = SAHelperSinusToNIntegral(fNvar - 1, init_theta) / beta;
01880       }
01881 
01882       map_dimension = fNvar;
01883       initialized = true;
01884    } // initializing is done.
01885 
01886    // generate uniform random number for sampling
01887    double u = fRandom->Rndm();
01888 
01889    // Find the two elements just greater than and less than u
01890    // using a binary search (O(log(N))).
01891    int lo = 0;
01892    int up = 10000;
01893    int mid;
01894 
01895    while (up != lo) {
01896       mid = ((up - lo + 1) / 2) + lo;
01897 
01898       if (u >= map_u[mid])
01899          lo = mid;
01900       else
01901          up = mid - 1;
01902    }
01903    up++;
01904 
01905    // perform linear interpolation:
01906    theta = map_theta[lo] + (u - map_u[lo]) / (map_u[up] - map_u[lo]) * (map_theta[up] - map_theta[lo]);
01907 
01908    return tan(theta);
01909 }
01910 
01911 // ---------------------------------------------------------
01912 
01913 double BCIntegrate::SAHelperSinusToNIntegral(int dim, double theta)
01914 {
01915    if (dim < 1)
01916       return theta;
01917    else if (dim == 1)
01918       return (1. - cos(theta));
01919    else if (dim == 2)
01920       return 0.5 * (theta - sin(theta) * cos(theta));
01921    else if (dim == 3)
01922       return (2. - sin(theta) * sin(theta) * cos(theta) - 2. * cos(theta)) / 3.;
01923    else
01924       return - pow(sin(theta), (double)(dim - 1)) * cos(theta) / (double)dim
01925          + (double)(dim - 1) / (double)dim
01926          * SAHelperSinusToNIntegral(dim - 2, theta);
01927 }
01928 
01929 // ---------------------------------------------------------
01930 
01931 void BCIntegrate::SetMode(std::vector <double> mode)
01932 {
01933    if((int)mode.size() == fNvar) {
01934       fBestFitParameters.clear();
01935       for (int i = 0; i < fNvar; i++)
01936          fBestFitParameters.push_back(mode[i]);
01937    }
01938 }
01939 
01940 // ---------------------------------------------------------
01941 
01942 void BCIntegrate::FCNLikelihood(int &npar, double * grad, double &fval, double * par, int flag)
01943 {
01944    // copy parameters
01945    std::vector <double> parameters;
01946 
01947    int n = global_this->GetNvar();
01948 
01949    for (int i = 0; i < n; i++)
01950       parameters.push_back(par[i]);
01951 
01952    fval = - global_this->LogEval(parameters);
01953 
01954    // delete parameters
01955    parameters.clear();
01956 }
01957 
01958 // ---------------------------------------------------------
01959 
01960 //void BCIntegrate::FindModeMCMC(int flag_run)
01961 void BCIntegrate::FindModeMCMC()
01962 {
01963    // call PreRun
01964 //   if (flag_run == 0)
01965    if (!fMCMCFlagPreRun)
01966       MCMCMetropolisPreRun();
01967 
01968    // find global maximum
01969    //   double probmax = (MCMCGetMaximumLogProb()).at(0);
01970 
01971    double probmax = 0;
01972 
01973    if ( int(fBestFitParameters.size()) == fNvar) {
01974       probmax = Eval( fBestFitParameters );
01975 //      fBestFitParameters = MCMCGetMaximumPoint(0);
01976    }
01977 
01978    // loop over all chains and find the maximum point
01979    for (int i = 0; i < fMCMCNChains; ++i) {
01980       double prob = exp( (MCMCGetMaximumLogProb()).at(i));
01981 
01982       // copy the point into the vector
01983       if ( (prob >= probmax && !fFlagIgnorePrevOptimization) || fFlagIgnorePrevOptimization) {
01984          probmax = prob;
01985 
01986          fBestFitParameters.clear();
01987          fBestFitParameterErrors.clear();
01988          fBestFitParameters = MCMCGetMaximumPoint(i);
01989 
01990          for (int j = 0; j < fNvar; ++j)
01991             fBestFitParameterErrors.push_back(-1.); // error undefined
01992 
01993          fOptimizationMethodMode = BCIntegrate::kOptMetropolis;
01994       }
01995    }
01996 }
01997 
01998 // ---------------------------------------------------------
01999 void BCIntegrate::SetCubaIntegrationMethod(BCIntegrate::BCCubaMethod type)
02000 {
02001 #ifdef HAVE_CUBA_H
02002    switch(type) {
02003       case BCIntegrate::kCubaVegas:
02004       case BCIntegrate::kCubaSuave:
02005          fCubaIntegrationMethod = type;
02006          return;
02007       case BCIntegrate::kCubaDivonne:
02008       case BCIntegrate::kCubaCuhre:
02009          BCLog::OutError(TString::Format(
02010             "BAT does not yet support global setting of Cuba integration method to %s. "
02011             "To use this method use explicit call to CubaIntegrate() with arguments.",
02012             DumpCubaIntegrationMethod(type).c_str()));
02013          return;
02014       default:
02015          BCLog::OutError(TString::Format("Integration method of type %d is not defined for Cuba",type));
02016          return;
02017    }
02018 #else
02019    BCLog::OutWarning("!!! This version of BAT is compiled without Cuba.");
02020    BCLog::OutWarning("    Setting Cuba integration method will have no effect.");
02021 #endif
02022 }
02023 
02024 // ---------------------------------------------------------
02025 int BCIntegrate::CubaIntegrand(const int *ndim, const double xx[],
02026                                               const int *ncomp, double ff[], void *userdata)
02027 {
02028 #ifdef HAVE_CUBA_H
02029    // scale variables
02030    double jacobian = 1.0;
02031 
02032    std::vector<double> scaled_parameters;
02033 
02034    for (int i = 0; i < *ndim; i++) {
02035       double range = global_this->fx->at(i)->GetUpperLimit() -  global_this->fx->at(i)->GetLowerLimit();
02036 
02037       // multiply range to jacobian
02038       jacobian *= range;
02039 
02040       // get the scaled parameter value
02041       scaled_parameters.push_back(global_this->fx->at(i)->GetLowerLimit() + xx[i] * range);
02042    }
02043 
02044    // call function to integrate
02045    ff[0] = global_this->Eval(scaled_parameters);
02046 
02047    // multiply jacobian
02048    ff[0] *= jacobian;
02049 
02050    // multiply fudge factor
02051    ff[0] *= 1e99;
02052 
02053    // remove parameter vector
02054    scaled_parameters.clear();
02055 #else
02056    BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
02057    BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
02058 #endif
02059 
02060     return 0;
02061 }
02062 
02063 // ---------------------------------------------------------
02064 double BCIntegrate::CubaIntegrate()
02065 {
02066 #ifdef HAVE_CUBA_H
02067    std::vector<double> parameters_double;
02068    std::vector<double>    parameters_int;
02069 
02070    parameters_double.push_back(fRelativePrecision);
02071    parameters_double.push_back(fAbsolutePrecision);
02072 
02073    parameters_int.push_back(fCubaVerbosity);
02074    parameters_int.push_back(fCubaMinEval);
02075    parameters_int.push_back(fCubaMaxEval);
02076 
02077    switch (fCubaIntegrationMethod) {
02078       case BCIntegrate::kCubaSuave:
02079          parameters_int.push_back(fCubaSuaveNNew);
02080          parameters_double.push_back(fCubaSuaveFlatness);
02081          break;
02082       case BCIntegrate::kCubaDivonne:
02083          break;
02084       case BCIntegrate::kCubaCuhre:
02085          break;
02086       default: // if unknown method run Vegas. Shouldn't ever happen anyway
02087       case BCIntegrate::kCubaVegas:
02088          parameters_int.push_back(fCubaVegasNStart);
02089          parameters_int.push_back(fCubaVegasNIncrease);
02090    }
02091 
02092    // print to log
02093    BCLog::OutDebug( Form(" --> Running Cuba/%s integation over %i dimensions.", DumpCubaIntegrationMethod().c_str(), fNvar));
02094    BCLog::OutDebug( Form(" --> Maximum number of iterations: %i", fCubaMaxEval));
02095    BCLog::OutDebug( Form(" --> Aimed relative precision:     %e", fRelativePrecision));
02096 
02097    return CubaIntegrate(fCubaIntegrationMethod, parameters_double, parameters_int);
02098 #else
02099    BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
02100    BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
02101    return 0.;
02102 #endif
02103 }
02104 
02105 // ---------------------------------------------------------
02106 double BCIntegrate::CubaIntegrate(BCIntegrate::BCCubaMethod method, std::vector<double> parameters_double, std::vector<double> parameters_int)
02107 {
02108 #ifdef HAVE_CUBA_H
02109    const int NDIM      = int(fx ->size());
02110    const int NCOMP     = 1;
02111    const int USERDATA  = 0;
02112    const int SEED      = 0;
02113    const int NBATCH    = 1000;
02114    const int GRIDNO    = -1;
02115    const char*STATEFILE = "";
02116 
02117    const double EPSREL = parameters_double[0];
02118    const double EPSABS = parameters_double[1];
02119    const int VERBOSE   = int(parameters_int[0]);
02120    const int MINEVAL   = int(parameters_int[1]);
02121    const int MAXEVAL   = int(parameters_int[2]);
02122 
02123    int neval;
02124    int fail;
02125    int nregions;
02126    double integral[NCOMP];
02127    double error[NCOMP];
02128    double prob[NCOMP];
02129 
02130    global_this = this;
02131 
02132    integrand_t an_integrand = &BCIntegrate::CubaIntegrand;
02133 
02134    if (method == 0) {
02135       // set VEGAS specific parameters
02136       const int NSTART    = int(parameters_int[3]);
02137       const int NINCREASE = int(parameters_int[4]);
02138 
02139       // call VEGAS integration method
02140          Vegas(NDIM, NCOMP, an_integrand, USERDATA,
02141                   EPSREL, EPSABS, VERBOSE, SEED,
02142                   MINEVAL, MAXEVAL, NSTART, NINCREASE, NBATCH,
02143                   GRIDNO, STATEFILE,
02144                   &neval, &fail, integral, error, prob);
02145 
02146          // interface for Cuba version 1.5
02147          /*
02148             Vegas(NDIM, NCOMP, an_integrand,
02149             EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
02150             NSTART, NINCREASE,
02151             &neval, &fail, integral, error, prob);
02152          */
02153    }
02154    else if (method == 1) {
02155       // set SUAVE specific parameters
02156        //      const int LAST     = int(parameters_int[3]);
02157       const int NNEW        = int(parameters_int[3]);
02158       const double FLATNESS = parameters_double[2];
02159 
02160       // call SUAVE integration method
02161       Suave(NDIM, NCOMP, an_integrand,
02162                   USERDATA, EPSREL, EPSABS, VERBOSE, SEED,
02163                   MINEVAL, MAXEVAL,
02164                   NNEW, FLATNESS,
02165                   &nregions, &neval, &fail, integral, error, prob);
02166 
02167          // interface for Cuba version 1.5
02168          /*
02169       Suave(NDIM, NCOMP, an_integrand,
02170          EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
02171          NNEW, FLATNESS,
02172          &nregions, &neval, &fail, integral, error, prob);
02173          */
02174    }
02175    else if (method == 2) {
02176       // set DIVONNE specific parameters
02177       const int KEY1         = int(parameters_int[3]);
02178       const int KEY2         = int(parameters_int[4]);
02179       const int KEY3         = int(parameters_int[5]);
02180       const int MAXPASS      = int(parameters_int[6]);
02181       const int BORDER       = int(parameters_int[7]);
02182       const int MAXCHISQ     = int(parameters_int[8]);
02183       const int MINDEVIATION = int(parameters_int[9]);
02184       const int NGIVEN       = int(parameters_int[10]);
02185       const int LDXGIVEN     = int(parameters_int[11]);
02186       const int NEXTRA       = int(parameters_int[12]);
02187          const int FLAGS    = 0;
02188 
02189          // call DIVONNE integration method
02190          Divonne(NDIM, NCOMP, an_integrand, USERDATA,
02191                      EPSREL, EPSABS, FLAGS, SEED, MINEVAL, MAXEVAL,
02192                      KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
02193                      NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
02194                      &nregions, &neval, &fail, integral, error, prob);
02195 
02196 
02197       // interface for Cuba version 1.5
02198          /*
02199             Divonne(NDIM, NCOMP, an_integrand,
02200             EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
02201             KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
02202             NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
02203             &nregions, &neval, &fail, integral, error, prob);
02204          */
02205    }
02206    else if (method == 3) {
02207       // set CUHRE specific parameters
02208        //      const int LAST = int(parameters_int[3]);
02209       const int KEY  = int(parameters_int[4]);
02210          const int FLAGS    = 0;
02211 
02212       // call CUHRE integration method
02213          Cuhre(NDIM, NCOMP, an_integrand, USERDATA,
02214                   EPSREL, EPSABS, FLAGS, MINEVAL, MAXEVAL, KEY,
02215                   &nregions, &neval, &fail, integral, error, prob);
02216 
02217          // interface for Cuba version 1.5
02218          /*
02219             Cuhre(NDIM, NCOMP, an_integrand,
02220             EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, KEY,
02221             &nregions, &neval, &fail, integral, error, prob);
02222          */
02223    }
02224    else {
02225       std::cout << " Integration method not available. " << std::endl;
02226       integral[0] = -1e99;
02227    }
02228 
02229    if (fail != 0) {
02230       std::cout << " Warning, integral did not converge with the given set of parameters. "<< std::endl;
02231       std::cout << " neval    = " << neval       << std::endl;
02232       std::cout << " fail     = " << fail        << std::endl;
02233       std::cout << " integral = " << integral[0] << std::endl;
02234       std::cout << " error    = " << error[0]    << std::endl;
02235       std::cout << " prob     = " << prob[0]     << std::endl;
02236    }
02237 
02238    return integral[0] / 1e99;
02239 #else
02240    BCLog::OutError("!!! This version of BAT is compiled without Cuba.");
02241    BCLog::OutError("    Use other integration methods or install Cuba and recompile BAT.");
02242    return 0.;
02243 #endif
02244 }
02245 
02246 // ---------------------------------------------------------
02247 void BCIntegrate::MCMCIterationInterface()
02248 {
02249    // what's within this method will be executed
02250    // for every iteration of the MCMC
02251 
02252    // fill error band
02253    MCMCFillErrorBand();
02254 
02255    // do user defined stuff
02256    MCMCUserIterationInterface();
02257 }
02258 
02259 // ---------------------------------------------------------
02260 void BCIntegrate::MCMCFillErrorBand()
02261 {
02262    if (!fFillErrorBand)
02263       return;
02264 
02265    // function fitting
02266    if (fFitFunctionIndexX < 0)
02267       return;
02268 
02269    // loop over all possible x values ...
02270    if (fErrorBandContinuous) {
02271       double x = 0;
02272       for (int ix = 0; ix < fErrorBandNbinsX; ix++) {
02273          // calculate x
02274          x = fErrorBandXY->GetXaxis()->GetBinCenter(ix + 1);
02275 
02276          // calculate y
02277          std::vector <double> xvec;
02278          xvec.push_back(x);
02279 
02280          // loop over all chains
02281          for (int ichain = 0; ichain < MCMCGetNChains(); ++ichain) {
02282             // calculate y
02283             double y = FitFunction(xvec, MCMCGetx(ichain));
02284 
02285             // fill histogram
02286             fErrorBandXY->Fill(x, y);
02287          }
02288 
02289          xvec.clear();
02290       }
02291    }
02292    // ... or evaluate at the data point x-values
02293    else {
02294       int ndatapoints = int(fErrorBandX.size());
02295       double x = 0;
02296 
02297       for (int ix = 0; ix < ndatapoints; ++ix) {
02298          // calculate x
02299          x = fErrorBandX.at(ix);
02300 
02301          // calculate y
02302          std::vector <double> xvec;
02303          xvec.push_back(x);
02304 
02305          // loop over all chains
02306          for (int ichain = 0; ichain < MCMCGetNChains(); ++ichain) {
02307             // calculate y
02308             double y = FitFunction(xvec, MCMCGetx(ichain));
02309 
02310             // fill histogram
02311             fErrorBandXY->Fill(x, y);
02312          }
02313 
02314          xvec.clear();
02315       }
02316    }
02317 }
02318 
02319 // ---------------------------------------------------------
02320 void BCIntegrate::SAInitialize()
02321 {
02322    fSAx.clear();
02323    fSAx.assign(fNvar, 0.0);
02324 }
02325 
02326 // ---------------------------------------------------------
02327 std::string BCIntegrate::DumpIntegrationMethod(BCIntegrate::BCIntegrationMethod type)
02328 {
02329    switch(type) {
02330       case BCIntegrate::kIntMonteCarlo:
02331          return "Sampled Mean Monte Carlo";
02332       case BCIntegrate::kIntImportance:
02333          return "Importance Sampling";
02334       case BCIntegrate::kIntMetropolis:
02335          return "Metropolis";
02336       case BCIntegrate::kIntCuba:
02337          return "Cuba";
02338       default:
02339          return "Undefined";
02340    }
02341 }
02342 
02343 // ---------------------------------------------------------
02344 std::string BCIntegrate::DumpMarginalizationMethod(BCIntegrate::BCMarginalizationMethod type)
02345 {
02346    switch(type) {
02347       case BCIntegrate::kMargMonteCarlo:
02348          return "Monte Carlo Integration";
02349       case BCIntegrate::kMargMetropolis:
02350          return "Metropolis MCMC";
02351       default:
02352          return "Undefined";
02353    }
02354 }
02355 
02356 // ---------------------------------------------------------
02357 std::string BCIntegrate::DumpOptimizationMethod(BCIntegrate::BCOptimizationMethod type)
02358 {
02359    switch(type) {
02360       case BCIntegrate::kOptSA:
02361          return "Simulated Annealing";
02362       case BCIntegrate::kOptMetropolis:
02363          return "Metropolis MCMC";
02364       case BCIntegrate::kOptMinuit:
02365          return "Minuit";
02366       default:
02367          return "Undefined";
02368    }
02369 }
02370 
02371 // ---------------------------------------------------------
02372 std::string BCIntegrate::DumpCubaIntegrationMethod(BCIntegrate::BCCubaMethod type)
02373 {
02374    switch(type) {
02375       case BCIntegrate::kCubaVegas:
02376          return "Vegas";
02377       case BCIntegrate::kCubaSuave:
02378          return "Suave";
02379       case BCIntegrate::kCubaDivonne:
02380          return "Divonne";
02381       case BCIntegrate::kCubaCuhre:
02382          return "Cuhre";
02383       default:
02384          return "Undefined";
02385    }
02386 }
02387 
02388 // ---------------------------------------------------------
02389