%%cpp -d
#include "TFile.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TStopwatch.h"

#include "RooWorkspace.h"
#include "RooAbsData.h"
#include "RooRandom.h"
#include "RooRealSumPdf.h"
#include "RooNumIntConfig.h"

#include "RooStats/ModelConfig.h"
#include "RooStats/ToyMCImportanceSampler.h"
#include "RooStats/HypoTestResult.h"
#include "RooStats/HypoTestPlot.h"
#include "RooStats/SamplingDistribution.h"
#include "RooStats/ProfileLikelihoodTestStat.h"
#include "RooStats/SimpleLikelihoodRatioTestStat.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"

#include "RooStats/FrequentistCalculator.h"
#include "TSystem.h"

#include <vector>

using namespace RooFit;
using namespace RooStats;

const char *infile = "";
const char *workspaceName = "channel1";
const char *modelConfigNameSB = "ModelConfig";
const char *dataName = "obsData";
int toys = 1000;
double poiValueForBackground = 0.0;
double poiValueForSignal = 1.0;

const char *filename = "";
if (!strcmp(infile, "")) {
   filename = "results/example_channel1_GammaExample_model.root";
   bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
   // if file does not exists generate with histfactory
   if (!fileExist) {
      // Normally this would be run on the command line
      cout << "will run standard hist2workspace example" << endl;
      gROOT->ProcessLine(".! prepareHistFactory .");
      gROOT->ProcessLine(".! hist2workspace config/example.xml");
      cout << "\n\n---------------------" << endl;
      cout << "Done creating example input" << endl;
      cout << "---------------------\n\n" << endl;
   }

} else
   filename = infile;

TFile *file = TFile::Open(filename);

if (!file) {
   cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
   return -1;
}

TStopwatch *mn_t = new TStopwatch;
mn_t->Start();

RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
   cout << "workspace not found" << endl;
   return -1.0;
}

ModelConfig *mc = (ModelConfig *)w->obj(modelConfigNameSB);

RooAbsData *data = w->data(dataName);

if (!data || !mc) {
   w->Print();
   cout << "data or ModelConfig was not found" << endl;
   return -1.0;
}

RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
firstPOI->setVal(poiValueForSignal);
mc->SetSnapshot(*mc->GetParametersOfInterest());

ModelConfig *mcNull = mc->Clone("ModelConfigNull");
firstPOI->setVal(poiValueForBackground);
mcNull->SetSnapshot(*(RooArgSet *)mcNull->GetParametersOfInterest()->snapshot());

ProfileLikelihoodTestStat *plts = new ProfileLikelihoodTestStat(*mc->GetPdf());
plts->SetOneSidedDiscovery(true);
plts->SetVarName("q_{0}/2");

ToyMCSampler toymcs(*plts, 50);

if (!mc->GetPdf()->canBeExtended()) {
   if (data->numEntries() == 1) {
      toymcs.SetNEventsPerToy(1);
   } else
      cout << "Not sure what to do about this model" << endl;
}

ProofConfig pc(*w, 2, "", false);

FrequentistCalculator freqCalc(*data, *mc, *mcNull, &toymcs);
freqCalc.SetToys(toys, toys); // null toys, alt toys

HypoTestResult *freqCalcResult = freqCalc.GetHypoTest();
freqCalcResult->GetNullDistribution()->SetTitle("b only");
freqCalcResult->GetAltDistribution()->SetTitle("s+b");
freqCalcResult->Print();
double pvalue = freqCalcResult->NullPValue();

mn_t->Stop();
cout << "total CPU time: " << mn_t->CpuTime() << endl;
cout << "total real time: " << mn_t->RealTime() << endl;

TCanvas *c1 = new TCanvas();
HypoTestPlot *plot = new HypoTestPlot(*freqCalcResult, 100, -0.49, 9.51);
plot->SetLogYaxis(true);

int nPOI = 1;
TF1 *f = new TF1("f", TString::Format("1*ROOT::Math::chisquared_pdf(2*x,%d,0)", nPOI), 0, 20);
f->SetLineColor(kBlack);
f->SetLineStyle(7);
plot->AddTF1(f, TString::Format("#chi^{2}(2x,%d)", nPOI));

plot->Draw();
c1->SaveAs("standard_discovery_output.pdf");

return pvalue;

gROOT->GetListOfCanvases()->Draw()