Bayesian Analysis Toolkit: BCTemplateEnsembleTest.cxx Source File

Go to the documentation of this file.
00001 #include <TH1D.h>
00002 #include <TFile.h>
00003 #include <TTree.h>
00004 #include <TROOT.h>
00005 #include <TRandom3.h>
00006 
00007 #include <math.h>
00008 #include <vector>
00009 #include <iostream>
00010 
00011 #include <BAT/BCLog.h>
00012 #include <BAT/BCH1D.h>
00013 
00014 #include "BCTemplateFitter.h"
00015 #include "BCTemplateEnsembleTest.h"
00016 
00017 using namespace std;
00018 
00019 // ---------------------------------------------------------
00020 BCTemplateEnsembleTest::BCTemplateEnsembleTest()
00021    : fTemplateFitter(0)
00022    , fFile(0)
00023    , fTree(0)
00024    , fEnsembleCounter(0)
00025    , fEnsembleExpectation(0)
00026    , fNEnsembles(0)
00027    , fFlagMCMC(false)
00028 {
00029    fRandom = new TRandom3(0);
00030 }
00031 
00032 // ---------------------------------------------------------
00033 BCTemplateEnsembleTest::~BCTemplateEnsembleTest()
00034 {
00035    if (fRandom)
00036       delete fRandom;
00037 }
00038 
00039 // ---------------------------------------------------------
00040 int BCTemplateEnsembleTest::SetEnsembleTemplate(TH1D hist)
00041 {
00042    // calculate integral
00043    double integral = hist.Integral();
00044 
00045    // check if integral is ok
00046    if (integral <= 0) {
00047       std::cout << "Template not valid. Integral is lower or equal to 0." << std::endl;
00048       return 0;
00049    }
00050 
00051    // set template
00052    fEnsembleTemplate = hist;
00053 
00054    // scale template
00055    fEnsembleTemplate.Scale(1.0/integral);
00056 
00057    // no error
00058    return 1;
00059 };
00060 
00061 // ---------------------------------------------------------
00062 int BCTemplateEnsembleTest::PerformEnsembleTest()
00063 {
00064    // set log level to nothing
00065    BCLog::LogLevel ll = BCLog::GetLogLevelScreen();
00066    BCLog::SetLogLevel(BCLog::nothing);
00067 
00068    // Prepare tree
00069    PrepareTree();
00070 
00071    // loop over ensembles
00072    for(int j = 0; j < fNEnsembles; j++){
00073 
00074       // print status
00075       if ((j+1) % 100 == 0 && j > 0)
00076          cout << "Fraction of ensembles analyzed: " << double(j+1) / double(fNEnsembles) * 100 << "%" << std::endl;
00077 
00078       // create new ensemble
00079       TH1D * ensemble = BuildEnsemble();
00080 
00081       // set ensemble as new data set
00082       fTemplateFitter->SetData(*ensemble);
00083 
00084       // find mode
00085       fTemplateFitter->FindMode();
00086 
00087       // perform MCMC
00088       if(fFlagMCMC) {
00089          fTemplateFitter->MarginalizeAll();
00090 
00091          // get number of parameters
00092          int npar = fTemplateFitter->GetNParameters();
00093 
00094          // loop over parameters and set tree variables
00095          for (int i = 0; i < npar; ++i) {
00096             BCH1D * hist = fTemplateFitter->GetMarginalized(fTemplateFitter->GetParameter(i));
00097             fOutParModeMarg[i]       = hist->GetMode();
00098             fOutParMedianMarg[i]     = hist->GetMedian();
00099             fOutParMeanMarg[i]       = hist->GetMean();
00100             fOutParRMSMarg[i]        = hist->GetRMS();
00101             fOutParErrorUpMarg[i]    = hist->GetQuantile(0.84)-hist->GetMode();
00102             fOutParErrorDownMarg[i]  = hist->GetMode()-hist->GetQuantile(0.16);
00103             fOutParQuantile5Marg[i]  = hist->GetQuantile(0.05);
00104             fOutParQuantile10Marg[i] = hist->GetQuantile(0.10);
00105             fOutParQuantile90Marg[i] = hist->GetQuantile(0.90);
00106             fOutParQuantile95Marg[i] = hist->GetQuantile(0.95);
00107          }
00108       }
00109 
00110       if (fFlagMCMC) {
00111          int nratios = fTemplateFitter->GetNRatios();
00112          for (int i = 0; i < nratios; ++i) {
00113             TH1D histtemp = fTemplateFitter->GetHistRatio1D(i);
00114             BCH1D * hist = new BCH1D( &histtemp );
00115             fOutRatioModeMarg[i]       = hist->GetMode();
00116             fOutRatioMedianMarg[i]     = hist->GetMedian();
00117             fOutRatioMeanMarg[i]       = hist->GetMean();
00118             fOutRatioRMSMarg[i]        = hist->GetRMS();
00119             fOutRatioErrorUpMarg[i]    = hist->GetQuantile(0.84)-hist->GetMode();
00120             fOutRatioErrorDownMarg[i]  = hist->GetMode()-hist->GetQuantile(0.16);
00121             fOutRatioQuantile5Marg[i]  = hist->GetQuantile(0.05);
00122             fOutRatioQuantile10Marg[i] = hist->GetQuantile(0.10);
00123             fOutRatioQuantile90Marg[i] = hist->GetQuantile(0.90);
00124             fOutRatioQuantile95Marg[i] = hist->GetQuantile(0.95);
00125          }
00126       }
00127 
00128       // set tree variables
00129       fOutParModeGlobal      = fTemplateFitter->GetBestFitParameters();
00130       fOutParErrorUpGlobal   = fTemplateFitter->GetBestFitParameterErrors();
00131       fOutParErrorDownGlobal = fTemplateFitter->GetBestFitParameterErrors();
00132       fOutChi2               = fTemplateFitter->CalculateChi2();
00133       fOutNDF                = fTemplateFitter->GetNDF();
00134       fOutChi2Prob           = fTemplateFitter->CalculateChi2Prob();
00135       fOutKSProb             = fTemplateFitter->CalculateKSProb();
00136       fOutPValue             = fTemplateFitter->CalculatePValue();
00137       fOutNEvents            = int(fTemplateFitter->GetData().Integral());
00138 
00139       // fill the tree
00140       fTree->Fill();
00141    }
00142 
00143    // reset log level
00144    BCLog::SetLogLevel(ll);
00145 
00146    // no error
00147    return 1;
00148 }
00149 
00150 //---------------------------------------------------------------------------------------------------------
00151 TH1D* BCTemplateEnsembleTest::BuildEnsemble()
00152 {
00153    // get histogram parameters
00154   int nbins   = fEnsembleTemplate.GetNbinsX();
00155   double xmin = fEnsembleTemplate.GetXaxis()->GetXmin();
00156   double xmax = fEnsembleTemplate.GetXaxis()->GetXmax();
00157 
00158    // create new ensemble
00159   TH1D* ensemble = new TH1D("", "", nbins, xmin, xmax);
00160 
00161    // increase ensemble counter
00162    fEnsembleCounter++;
00163 
00164    // loop over bins and fill them
00165   for(int i = 1; i <= nbins; ++i){
00166     double p      = fEnsembleTemplate.GetBinContent(i);
00167     double lambda = p * fEnsembleExpectation;
00168       double n      = gRandom -> Poisson(lambda);
00169 
00170       // set the bin content
00171     ensemble->SetBinContent(i, n);
00172   }
00173 
00174    // return the ensemble histogram
00175   return ensemble;
00176 }
00177 
00178 //---------------------------------------------------------------------------------------------------------
00179 int BCTemplateEnsembleTest::Write(const char * filename)
00180 {
00181    // open file
00182    fFile = new TFile(filename, "RECREATE");
00183 
00184    // write tree
00185    fTree->Write();
00186 
00187    // close file
00188    fFile->Close();
00189 
00190    // free memory
00191    delete fFile;
00192 
00193    // no error
00194    return 1;
00195 }
00196 
00197 //---------------------------------------------------------------------------------------------------------
00198 int BCTemplateEnsembleTest::PrepareTree()
00199 {
00200    // delete old tree if necessary
00201    if (fTree)
00202       delete fTree;
00203 
00204    // create new tree
00205    fTree = new TTree("fTree", "fTree");
00206 
00207    // get number of parameters and ratios
00208    int npar = fTemplateFitter->GetNParameters();
00209    int nratios = fTemplateFitter->GetNRatios();
00210 
00211    // initialize variables
00212    fOutParModeGlobal.assign(npar, 0);
00213    fOutParErrorUpGlobal.assign(npar, 0);
00214    fOutParErrorDownGlobal.assign(npar, 0);
00215    fOutParModeMarg.assign(npar, 0);
00216    fOutParMeanMarg.assign(npar, 0);
00217    fOutParMedianMarg.assign(npar, 0);
00218    fOutParRMSMarg.assign(npar, 0);
00219    fOutParErrorUpMarg.assign(npar, 0);
00220    fOutParErrorDownMarg.assign(npar, 0);
00221    fOutParQuantile5Marg.assign(npar, 0);
00222    fOutParQuantile10Marg.assign(npar, 0);
00223    fOutParQuantile90Marg.assign(npar, 0);
00224    fOutParQuantile95Marg.assign(npar, 0);
00225    fOutRatioModeMarg.assign(nratios, 0);
00226    fOutRatioMeanMarg.assign(nratios, 0);
00227    fOutRatioMedianMarg.assign(nratios, 0);
00228    fOutRatioRMSMarg.assign(nratios, 0);
00229    fOutRatioErrorUpMarg.assign(nratios, 0);
00230    fOutRatioErrorDownMarg.assign(nratios, 0);
00231    fOutRatioQuantile5Marg.assign(nratios, 0);
00232    fOutRatioQuantile10Marg.assign(nratios, 0);
00233    fOutRatioQuantile90Marg.assign(nratios, 0);
00234    fOutRatioQuantile95Marg.assign(nratios, 0);
00235 
00236    fTree->Branch("chi2",     &fOutChi2,     "chi2/D");
00237    fTree->Branch("ndf",      &fOutNDF,      "ndf/I");
00238    fTree->Branch("chi2prob", &fOutChi2Prob, "chi2 prob probability/D");
00239    fTree->Branch("KSprob",   &fOutKSProb,   "KS probability/D");
00240    fTree->Branch("pvalue",   &fOutPValue,   "p-value/D");
00241    fTree->Branch("nevents",  &fOutNEvents,  "n events/I");
00242 
00243    for (int i = 0; i < npar; ++i) {
00244       // add branches
00245       fTree->Branch(Form("par_global_mode_par_%i", i),       &fOutParModeGlobal[i],      Form("par_global_Mode_par_%i/D", i));
00246       fTree->Branch(Form("par_global_error_up_par_%i", i),   &fOutParErrorUpGlobal[i],   Form("par_global_error_up_par_%i/D", i));
00247       fTree->Branch(Form("par_global_error_down_par_%i", i), &fOutParErrorDownGlobal[i], Form("par_global_error_down_par_%i/D", i));
00248 
00249       if(fFlagMCMC) {
00250          fTree->Branch(Form("par_marg_mode_par_%i", i),       &fOutParModeMarg[i],       Form("par_marg_mode_par_%i/D", i));
00251          fTree->Branch(Form("par_marg_mean_par_%i", i),       &fOutParMeanMarg[i],       Form("par_marg_mean_par_%i/D", i));
00252          fTree->Branch(Form("par_marg_median_par_%i", i),     &fOutParMedianMarg[i],     Form("par_marg_median_par_%i/D", i));
00253          fTree->Branch(Form("par_marg_rms_par_%i", i),        &fOutParRMSMarg[i],        Form("par_marg_rms_par_%i/D", i));
00254          fTree->Branch(Form("par_marg_error_up_par_%i", i),   &fOutParErrorUpMarg[i],    Form("par_marg_ErrorUp_par_%i/D", i));
00255          fTree->Branch(Form("par_marg_error_down_par_%i", i), &fOutParErrorDownMarg[i],  Form("par_marg_error_down_par_%i/D", i));
00256          fTree->Branch(Form("par_marg_quantile5_par_%i", i),  &fOutParQuantile5Marg[i],  Form("par_marg_Quantile5_par_%i/D", i));
00257          fTree->Branch(Form("par_marg_quantile10_par_%i", i), &fOutParQuantile10Marg[i], Form("par_marg_Quantile10_par_%i/D", i));
00258          fTree->Branch(Form("par_marg_quantile90_par_%i", i), &fOutParQuantile90Marg[i], Form("par_marg_Quantile90_par_%i/D", i));
00259          fTree->Branch(Form("par_marg_quantile95_par_%i", i), &fOutParQuantile95Marg[i], Form("par_marg_Quantile95_par_%i/D", i));
00260       }
00261    }
00262 
00263    if (fFlagMCMC) {
00264       for (int i = 0; i < nratios; ++i) {
00265          fTree->Branch(Form("ratio_marg_mode_ratio_%i", i),       &fOutRatioModeMarg[i],       Form("ratio_marg_mode_ratio_%i/D", i));
00266          fTree->Branch(Form("ratio_marg_mean_ratio_%i", i),       &fOutRatioMeanMarg[i],       Form("ratio_marg_mean_ratio_%i/D", i));
00267          fTree->Branch(Form("ratio_marg_median_ratio_%i", i),     &fOutRatioMedianMarg[i],     Form("ratio_marg_median_ratio_%i/D", i));
00268          fTree->Branch(Form("ratio_marg_rms_ratio_%i", i),        &fOutRatioRMSMarg[i],        Form("ratio_marg_rms_ratio_%i/D", i));
00269          fTree->Branch(Form("ratio_marg_error_up_ratio_%i", i),   &fOutRatioErrorUpMarg[i],    Form("ratio_marg_ErrorUp_ratio_%i/D", i));
00270          fTree->Branch(Form("ratio_marg_error_down_ratio_%i", i), &fOutRatioErrorDownMarg[i],  Form("ratio_marg_error_down_ratio_%i/D", i));
00271          fTree->Branch(Form("ratio_marg_quantile5_ratio_%i", i),  &fOutRatioQuantile5Marg[i],  Form("ratio_marg_Quantile5_ratio_%i/D", i));
00272          fTree->Branch(Form("ratio_marg_quantile10_ratio_%i", i), &fOutRatioQuantile10Marg[i], Form("ratio_marg_Quantile10_ratio_%i/D", i));
00273          fTree->Branch(Form("ratio_marg_quantile90_ratio_%i", i), &fOutRatioQuantile90Marg[i], Form("ratio_marg_Quantile90_ratio_%i/D", i));
00274          fTree->Branch(Form("ratio_marg_quantile95_ratio_%i", i), &fOutRatioQuantile95Marg[i], Form("ratio_marg_Quantile95_ratio_%i/D", i));
00275       }
00276    }
00277 
00278    // no error
00279    return 1;
00280 }
00281 
00282 //---------------------------------------------------------------------------------------------------------