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