BayesianAnalysisToolkit: BCIntegrate.cxx Source File

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