Bayesian Analysis Toolkit: BCIntegrate.cxx Source File

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