Bayesian Analysis Toolkit: BCModelOutput.cxx Source File

00001 /*    
00002  * Copyright (C) 2008, Daniel Kollar and Kevin Kroeninger.    
00003  * All rights reserved.    
00004  *    
00005  * For the licensing terms see doc/COPYING.    
00006  */    
00007   
00008 // ---------------------------------------------------------   
00009  
00010 #include "BCModelOutput.h" 
00011 
00012 #include "TDirectory.h"
00013 
00014 // --------------------------------------------------------- 
00015 
00016 BCModelOutput::BCModelOutput() 
00017 {
00018 
00019    fIndex = 0; 
00020 
00021    fOutputFile = 0; 
00022    fMarkovChainTree = 0; 
00023    fAnalysisTree = 0; 
00024    fModel = 0; 
00025    fOutputFile = 0; 
00026 
00027 }
00028 
00029 // --------------------------------------------------------- 
00030 
00031 BCModelOutput::BCModelOutput(BCModel * model, const char * filename) 
00032 {
00033 
00034    // remeber current directory 
00035 
00036    TDirectory * dir = gDirectory; 
00037 
00038    fFilename = filename; 
00039 
00040    // create a new file 
00041 
00042    fOutputFile = new TFile(fFilename, "RECREATE"); 
00043 
00044    BCModelOutput(); 
00045 
00046    fModel = model;
00047 
00048    // initialize trees 
00049 
00050    this -> InitializeAnalysisTree(); 
00051 
00052    this -> InitializeMarkovChainTrees(); 
00053 
00054    // change again to the old directory 
00055 
00056    gDirectory = dir; 
00057 
00058 }
00059 
00060 // --------------------------------------------------------- 
00061 
00062 BCModelOutput::~BCModelOutput() 
00063 {
00064 
00065    if (fOutputFile)
00066       {
00067          fOutputFile -> Close(); 
00068 
00069          delete fOutputFile; 
00070       } 
00071 
00072    if (fAnalysisTree)
00073       delete fAnalysisTree; 
00074 
00075    if (fMarkovChainTree)
00076       delete fMarkovChainTree; 
00077 
00078    if (fOutputFile)
00079       delete fOutputFile; 
00080 
00081    fModel = 0; 
00082 
00083 }
00084 
00085 // --------------------------------------------------------- 
00086 
00087 BCModelOutput::BCModelOutput(const BCModelOutput & modeloutput)
00088 {
00089 
00090    modeloutput.Copy(*this); 
00091 
00092 }
00093 
00094 // --------------------------------------------------------- 
00095 
00096 BCModelOutput & BCModelOutput::operator = (const BCModelOutput & modeloutput)
00097 {
00098 
00099    if (this != &modeloutput) 
00100       modeloutput.Copy(* this); 
00101 
00102    return * this; 
00103 
00104 }
00105 
00106 // --------------------------------------------------------- 
00107 
00108 void BCModelOutput::SetFile(const char * filename) 
00109 {
00110 
00111    // delete the old file 
00112 
00113    if (fOutputFile)
00114       {
00115          fOutputFile -> Close(); 
00116 
00117          delete fOutputFile; 
00118       }
00119 
00120    // remember current directory 
00121 
00122    TDirectory * dir = gDirectory; 
00123 
00124    // create a new file 
00125 
00126    fOutputFile = new TFile(fFilename, "RECREATE");    
00127 
00128    // change back to the old directory 
00129 
00130    gDirectory = dir; 
00131 
00132    // initialize tree 
00133 
00134    this -> InitializeAnalysisTree(); 
00135 
00136 }
00137 
00138 // --------------------------------------------------------- 
00139 
00140 void BCModelOutput::WriteMarkovChain(bool flag)
00141 {
00142 
00143    if (fModel)
00144       fModel -> WriteMarkovChain(flag); 
00145 
00146 }
00147 
00148 // --------------------------------------------------------- 
00149    
00150 void BCModelOutput::FillAnalysisTree()
00151 {
00152 
00153    // get output values from model 
00154 
00155    fNParameters = fModel -> GetNParameters(); 
00156    fProbability_apriori   = fModel -> GetModelAPrioriProbability(); 
00157    fProbability_aposteriori = fModel -> GetModelAPosterioriProbability(); 
00158    
00159    // loop over parameters 
00160 
00161    int nparameters = fModel -> GetNParameters(); 
00162 
00163    for (int i = 0; i < nparameters; ++i)
00164       {
00165          BCParameter * parameter = fModel -> GetParameter(i); 
00166 
00167          if (fModel -> GetBestFitParameters().size() > 0)
00168             fMode_global[i] = fModel -> GetBestFitParameters().at(i); 
00169 
00170          if (fModel -> GetMarginalized(parameter -> GetName().data()))
00171             {
00172                fMode_marginalized[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetMode(); 
00173                fMean_marginalized[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetMean(); 
00174                fMedian_marginalized[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetMedian(); 
00175                fQuantile_05[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetQuantile(0.05); 
00176                fQuantile_10[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetQuantile(0.10); 
00177                fQuantile_16[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetQuantile(0.16); 
00178                fQuantile_84[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetQuantile(0.84); 
00179                fQuantile_90[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetQuantile(0.90); 
00180                fQuantile_95[i] = fModel -> GetMarginalized(parameter -> GetName().data()) -> GetQuantile(0.95); 
00181             }
00182       }
00183 
00184    // fill tree 
00185 
00186    fAnalysisTree -> Fill(); 
00187 
00188    // increase index 
00189 
00190    fIndex++; 
00191 
00192 }
00193 
00194 // ---------------------------------------------------------
00195 
00196 void BCModelOutput::WriteMarginalizedDistributions()
00197 {
00198 
00199    // remember current directory
00200    TDirectory * dir = gDirectory;
00201 
00202    // change to file
00203    fOutputFile -> cd();
00204 
00205    int nparameters = fModel -> GetNParameters();
00206    for (int i = 0; i < nparameters; ++i)
00207       fModel -> GetMarginalized(fModel -> GetParameter(i)) -> GetHistogram() -> Write();
00208 
00209    if (nparameters > 1)
00210    {
00211       for (int i = 0; i < nparameters - 1; ++i)
00212          for (int j = i + 1; j < nparameters; ++j)
00213             fModel -> GetMarginalized(fModel -> GetParameter(i),
00214                   fModel -> GetParameter(j)) -> GetHistogram() -> Write();
00215    }
00216 
00217    // return to old directory
00218    gDirectory = dir;
00219 
00220 }
00221 
00222 // ---------------------------------------------------------
00223    
00224 void BCModelOutput::Close()
00225 {
00226 
00227    // remember current directory 
00228 
00229    TDirectory * dir = gDirectory;
00230 
00231    // change to file 
00232 
00233    fOutputFile -> cd(); 
00234 
00235    // write analysis tree to file 
00236 
00237    if (fAnalysisTree -> GetEntries() > 0) 
00238       fAnalysisTree -> Write(); 
00239 
00240    // write markov chain tree to file 
00241 
00242    for (int i = 0; i < fModel -> MCMCGetNChains(); ++i)
00243       if (fMarkovChainTrees[i] -> GetEntries() > 0) 
00244          fMarkovChainTrees[i] -> Write(); 
00245 
00246    // write control plots to file 
00247 
00248    // if (fModel -> MCMCGetH1RValue())
00249    //   fModel -> MCMCGetH1RValue() -> Write(); 
00250 
00251    // if (fModel -> MCMCGetH1Efficiency())
00252    //   fModel -> MCMCGetH1Efficiency() -> Write(); 
00253 
00254    // close file 
00255 
00256    fOutputFile -> Close(); 
00257 
00258    // return to old directory 
00259 
00260    gDirectory = dir; 
00261 
00262 }
00263 
00264 // --------------------------------------------------------- 
00265    
00266 void BCModelOutput::InitializeAnalysisTree() 
00267 {
00268 
00269    // create new tree 
00270 
00271    fAnalysisTree = new TTree("AnalysisTree", "AnalysisTree"); 
00272 
00273    // set branch addresses 
00274 
00275    fAnalysisTree -> Branch("fIndex",                   &fIndex,                   "index/I"); 
00276    fAnalysisTree -> Branch("fNParameters",             &fNParameters,             "parameters/I"); 
00277    fAnalysisTree -> Branch("fProbability_apriori" ,    &fProbability_apriori,     "apriori probability/D"); 
00278    fAnalysisTree -> Branch("fProbability_aposteriori", &fProbability_aposteriori, "aposteriori probability/D"); 
00279    fAnalysisTree -> Branch("fMode_global",              fMode_global,             "mode (global) [parameters]/D"); 
00280    fAnalysisTree -> Branch("fMode_marginalized",        fMode_marginalized,       "mode (marginalized) [parameters]/D"); 
00281    fAnalysisTree -> Branch("fMean_marginalized",        fMean_marginalized,       "mean (marginalized)[parameters]/D"); 
00282    fAnalysisTree -> Branch("fMedian_marginalized",      fMedian_marginalized,     "median (marginalized)[parameters]/D"); 
00283    fAnalysisTree -> Branch("fQuantile_05" ,             fQuantile_05,             "quantile 5% [parameters]/D"); 
00284    fAnalysisTree -> Branch("fQuantile_10" ,             fQuantile_10,             "quantile 10% [parameters]/D"); 
00285    fAnalysisTree -> Branch("fQuantile_16" ,             fQuantile_16,             "quantile 16% [parameters]/D"); 
00286    fAnalysisTree -> Branch("fQuantile_84" ,             fQuantile_84,             "quantile 84% [parameters]/D"); 
00287    fAnalysisTree -> Branch("fQuantile_90" ,             fQuantile_90,             "quantile 90% [parameters]/D"); 
00288    fAnalysisTree -> Branch("fQuantile_95" ,             fQuantile_95,             "quantile 95% [parameters]/D"); 
00289 
00290 }
00291 
00292 // --------------------------------------------------------- 
00293    
00294 void BCModelOutput::InitializeMarkovChainTrees() 
00295 {
00296 
00297    // create new tree 
00298 
00299    fMarkovChainTrees.clear(); 
00300 
00301    for (int i = 0; i < fModel -> MCMCGetNChains(); ++i)
00302       {
00303          TTree * tree = new TTree(Form("MarkovChainTree_%i", i), "MarkovChainTree"); 
00304          fMarkovChainTrees.push_back(tree); 
00305       }
00306 
00307 //    // connect pointer to parameter vectors 
00308 
00309    fParameters    = fModel -> MCMCGetP2x(); 
00310    fIteration     = fModel -> MCMCGetP2NIterations(); 
00311    fLogLikelihood = fModel -> MCMCGetP2LogProbx(); 
00312 
00313    fNParameters = fModel -> MCMCGetNParameters(); 
00314 
00315    for (int i = 0; i < fModel -> MCMCGetNChains(); ++i)
00316       {
00317          fMarkovChainTrees[i] -> Branch("fIteration",      &(*fIteration)[i],   "iteration/I"); 
00318          fMarkovChainTrees[i] -> Branch("fNParameters",    &fNParameters,       "parameters/I"); 
00319          fMarkovChainTrees[i] -> Branch("fLogLikelihood",  &(*fLogLikelihood)[i], "log (likelihood)/D"); 
00320       }
00321 
00322    // loop over all parameters 
00323    
00324    for (int i = 0; i < fModel -> MCMCGetNChains(); ++i)
00325       for (int j = 0; j <  fModel -> MCMCGetNParameters(); ++j)
00326       {
00327          fMarkovChainTrees[i] -> Branch(Form("fParameter%i", j), 
00328                                                        &(*fParameters)[i * fModel -> MCMCGetNParameters() + j], 
00329                                                        Form("parameter %i/D", j)); 
00330       }
00331 
00332    fModel -> MCMCSetMarkovChainTrees(fMarkovChainTrees); 
00333 
00334 }
00335 
00336 // --------------------------------------------------------- 
00337 
00338 void BCModelOutput::Copy(BCModelOutput & modeloutput) const 
00339 {
00340 
00341    // don't copy the content 
00342 
00343    modeloutput.fModel           = this -> fModel; 
00344    modeloutput.fAnalysisTree    = this -> fAnalysisTree; 
00345    modeloutput.fMarkovChainTree = this -> fMarkovChainTree; 
00346 
00347 }
00348 
00349 // ---------------------------------------------------------