StandardProfileInspectorDemo.C

Go to the documentation of this file.
00001 // Standard demo of the ProfileInspector class
00002 /*
00003 StandardProfileInspectorDemo
00004 
00005 Author: Kyle Cranmer
00006 date: Dec. 2010
00007 
00008 This is a standard demo that can be used with any ROOT file 
00009 prepared in the standard way.  You specify:
00010  - name for input ROOT file
00011  - name of workspace inside ROOT file that holds model and data
00012  - name of ModelConfig that specifies details for calculator tools
00013  - name of dataset 
00014 
00015 With default parameters the macro will attempt to run the 
00016 standard hist2workspace example and read the ROOT file
00017 that it produces.
00018 
00019 The actual heart of the demo is only about 10 lines long.
00020 
00021 The ProfileInspector plots the conditional maximum likelihood estimate
00022 of each nuisance parameter in the model vs. the parameter of interest.
00023 (aka. profiled value of nuisance parameter vs. parameter of interest)
00024 (aka. best fit nuisance parameter with p.o.i fixed vs. parameter of interest)
00025 
00026 */
00027 
00028 #include "TFile.h"
00029 #include "TROOT.h"
00030 #include "TCanvas.h"
00031 #include "TList.h"
00032 #include "RooWorkspace.h"
00033 #include "RooAbsData.h"
00034 
00035 #include "RooStats/ModelConfig.h"
00036 #include "RooStats/ProfileInspector.h"
00037 
00038 using namespace RooFit;
00039 using namespace RooStats;
00040 
00041 void StandardProfileInspectorDemo(const char* infile = "",
00042                       const char* workspaceName = "proto",
00043                       const char* modelConfigName = "ModelConfig",
00044                       const char* dataName = "expData"){
00045 
00046   /////////////////////////////////////////////////////////////
00047   // First part is just to access a user-defined file 
00048   // or create the standard example file if it doesn't exist
00049   ////////////////////////////////////////////////////////////
00050   const char* filename = "";
00051   if (!strcmp(infile,""))
00052     filename = "results/example_channel1_GammaExample_model.root";
00053   else
00054     filename = infile;
00055   // Check if example input file exists
00056   TFile *file = TFile::Open(filename);
00057 
00058   // if input file was specified byt not found, quit
00059   if(!file && strcmp(infile,"")){
00060     cout <<"file not found" << endl;
00061     return;
00062   } 
00063 
00064   // if default file not found, try to create it
00065   if(!file ){
00066     // Normally this would be run on the command line
00067     cout <<"will run standard hist2workspace example"<<endl;
00068     gROOT->ProcessLine(".! prepareHistFactory .");
00069     gROOT->ProcessLine(".! hist2workspace config/example.xml");
00070     cout <<"\n\n---------------------"<<endl;
00071     cout <<"Done creating example input"<<endl;
00072     cout <<"---------------------\n\n"<<endl;
00073   }
00074 
00075   // now try to access the file again
00076   file = TFile::Open(filename);
00077   if(!file){
00078     // if it is still not there, then we can't continue
00079     cout << "Not able to run hist2workspace to create example input" <<endl;
00080     return;
00081   }
00082 
00083   
00084   /////////////////////////////////////////////////////////////
00085   // Tutorial starts here
00086   ////////////////////////////////////////////////////////////
00087 
00088   // get the workspace out of the file
00089   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
00090   if(!w){
00091     cout <<"workspace not found" << endl;
00092     return;
00093   }
00094 
00095   // get the modelConfig out of the file
00096   ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);
00097 
00098   // get the modelConfig out of the file
00099   RooAbsData* data = w->data(dataName);
00100 
00101   // make sure ingredients are found
00102   if(!data || !mc){
00103     w->Print();
00104     cout << "data or ModelConfig was not found" <<endl;
00105     return;
00106   }
00107 
00108   //////////////////////////////////////////////
00109   // now use the profile inspector
00110   ProfileInspector p;
00111   TList* list = p.GetListOfProfilePlots(*data,mc);
00112   
00113    // now make plots
00114   TCanvas* c1 = new TCanvas("c1","ProfileInspectorDemo",800,200);
00115   c1->Divide(list->GetSize());
00116   for(int i=0; i<list->GetSize(); ++i){
00117     c1->cd(i+1);
00118     list->At(i)->Draw("al");
00119   }
00120   
00121 }

Generated on Tue Jul 5 15:45:10 2011 for ROOT_528-00b_version by  doxygen 1.5.1