TestSignificance.C

// StandardFrequentistDiscovery

/*
 StandardFrequentistDiscovery

 Author: Sven Kreiss, Kyle Cranmer
 date: May 2012

 This is a standard demo that can be used with any ROOT file
 prepared in the standard way.  You specify:
 - name for input ROOT file
 - name of workspace inside ROOT file that holds model and data
 - name of ModelConfig that specifies details for calculator tools
 - name of dataset

 With default parameters the macro will attempt to run the
 standard hist2workspace example and read the ROOT file
 that it produces.
 */

#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/RatioOfProfiledLikelihoodsTestStat.h"
#include "RooStats/SimpleLikelihoodRatioTestStat.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/AsymptoticCalculator.h"

#include "RooStats/FrequentistCalculator.h"

#include <vector>

using namespace RooFit;
using namespace RooStats;
Double_t median(TH1 *h1) {
  //compute the median for 1-d histogram h1
  Int_t nbins = h1->GetXaxis()->GetNbins();
  Double_t *x = new Double_t[nbins];
  Double_t *y = new Double_t[nbins];
  for (Int_t i=0;i<nbins;i++) {
    x[i] = h1->GetXaxis()->GetBinCenter(i+1);
    y[i] = h1->GetBinContent(i+1);
  }
  Double_t median = TMath::Median(nbins,x,y);
  delete [] x;
  delete [] y;
  return median;
}

double TestSignificance ( TString infile = "hh_cp_12_combination_125.root",
			  const char* workspaceName = "workspace_hh_cp_12_combination_125",
			  const char* modelConfigNameSB = "ModelConfig",
			  const char* dataName = "asimovData",
			  int toys = 500,
			  double poiValueForCPodd = 0.0,
			  double poiValueForCPeven = 1.0 )  {

  // The workspace contains the model for s+b. The b model is "autogenerated"
  // by copying s+b and setting the one parameter of interest to zero.
  // To keep the script simple, multiple parameters of interest or different
  // functional forms of the b model are not supported.

  // for now, assume there is only one parameter of interest, and these are
  // its values:

  /////////////////////////////////////////////////////////////
  // First part is just to access a user-defined file
  // or create the standard example file if it doesn't exist
  ////////////////////////////////////////////////////////////
  TString filename = "";
  if (!strcmp(infile, ""))
    {
      cout << "Input file not specified!!!" << endl;
      exit(1);
    }
  else filename = infile;

  std::cout<<"Working with file: "<<filename<<std::endl;

  // Check if example input file exists
  TFile *file = TFile::Open(filename);

  // if input file was specified but not found, quit
  if (!file && strcmp(infile, "")) {
    cout << "File not found!!!" << endl;
    exit(1);
  }


  /////////////////////////////////////////////////////////////
  // Tutorial starts here
  ////////////////////////////////////////////////////////////

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

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

  // get the modelConfig out of the file
  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigNameSB);

  // get the data out of the file
  RooAbsData* data = w->data(dataName);

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

  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  firstPOI->setVal(poiValueForCPeven);
  mc->SetSnapshot(*mc->GetParametersOfInterest());
  RooArgSet altParams(*mc->GetSnapshot());

  // create null model
  ModelConfig *mcNull = mc->Clone("ModelConfigNull");
  firstPOI->setVal(poiValueForCPodd);
  mcNull->SetSnapshot(*(RooArgSet*)mcNull->GetParametersOfInterest()->snapshot());
  RooArgSet nullParams(*mcNull->GetSnapshot());

  // make asimov data for alt
  mc->LoadSnapshot();
  cout << "Variables for AsimovData:" << endl;
  w->allVars().Print("V");
  cout << endl;
  data = AsymptoticCalculator::GenerateAsimovData( *mc->GetPdf(), *mc->GetObservables() );

  // make asimov data for null
  mcNull->LoadSnapshot();
  dataNull = AsymptoticCalculator::GenerateAsimovData( *mcNull->GetPdf(), *mcNull->GetObservables() );

  // ----------------------------------------------------
  // Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
  // to use simultaneously with ToyMCSampler
  ProfileLikelihoodTestStat* plts =  new ProfileLikelihoodTestStat(*mc->GetPdf());
  //   plts->SetOneSidedDiscovery(true);

  SimpleLikelihoodRatioTestStat * slrts =  new SimpleLikelihoodRatioTestStat( *mcNull->GetPdf(),  *mc->GetPdf());
  slrts->SetAltParameters( *mc->GetSnapshot() );
  slrts->SetNullParameters( *mcNull->GetSnapshot() );
  RooArgSet nullParams(*mcNull->GetSnapshot());
  if (mcNull->GetNuisanceParameters()) nullParams.add(*mcNull->GetNuisanceParameters());
  slrts->SetNullParameters(nullParams);

  RooArgSet altParams(*mc->GetSnapshot());
  if (mc->GetNuisanceParameters()) altParams.add(*mc->GetNuisanceParameters());
  slrts->SetAltParameters(altParams);

  RatioOfProfiledLikelihoodsTestStat * ropl = new RatioOfProfiledLikelihoodsTestStat( *mcNull->GetPdf(), *mc->GetPdf(), mc->GetSnapshot() );
  //   plts->SetSubtractMLE( false );


  // ----------------------------------------------------
  // configure the ToyMCImportanceSampler with two test statistics
  ToyMCSampler toymcs(*ropl, 50);

  // Since this tool needs to throw toy MC the PDF needs to be
  // extended or the tool needs to know how many entries in a dataset
  // per pseudo experiment.
  // In the 'number counting form' where the entries in the dataset
  // are counts, and not values of discriminating variables, the
  // datasets typically only have one entry and the PDF is not
  // extended.
  if (!mc->GetPdf()->canBeExtended()) {
    if (data->numEntries() == 1) {
      toymcs.SetNEventsPerToy(1);
    } else cout << "Not sure what to do about this model" << endl;
  }

  // We can use PROOF to speed things along in parallel
  ProofConfig pc(*w, 1e2, "workers=5", false);
  toymcs.SetProofConfig(&pc);    // enable proof

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

  // // instantiate the calculator
  // HybridCalculator hybridCalc(*data, *mc, *mcNull, &toymcs);
  // hybridCalc.SetToys( toys, toys ); // null toys, alt toys


  // Run the calculator and print result
  HypoTestResult* CalcResult = freqCalc.GetHypoTest();
  CalcResult->GetNullDistribution()->SetTitle( "Hypothesis 1" );
  CalcResult->GetAltDistribution()->SetTitle( "Hypotehsis 2" );
  CalcResult->Print();
  double pvalue = CalcResult->NullPValue();

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

  // plot ===================
  // Get alt / null distribution
  int Nbin = 100; double xmin = -4.9999; double xmax = 5.0001;
  SamplingDistribution *alt = CalcResult->GetAltDistribution();
  SamplingDistribution *null = CalcResult->GetNullDistribution();
  TH1F *althist = new TH1F("ModelConfig","",Nbin,xmin,xmax);
  TH1F *nullhist = new TH1F("ModelConfigNull","",Nbin,xmin,xmax);
  althist->Sumw2();
  nullhist->Sumw2();
  vector<Double_t> altdistr = alt->GetSamplingDistribution();
  vector<Double_t> nulldistr = null->GetSamplingDistribution();
  for (int i=0;i<altdistr.size();i++) althist->Fill(altdistr.at(i));
  for (int i=0;i<nulldistr.size();i++) nullhist->Fill(nulldistr.at(i));
  double sum = 1.0;
  sum = althist->Integral();  althist->Scale(1./sum);
  sum = nullhist->Integral(); nullhist->Scale(1./sum);
  althist->SetLineWidth(2);
  nullhist->SetLineWidth(2);
  althist->SetLineColor(kBlue);
  nullhist->SetLineColor(kRed);

  TCanvas* c1 = new TCanvas();
  double ymin = 0.01;
  double ymax = TMath::Max(althist->GetMaximum(),nullhist->GetMaximum()) * 1.5;
  TH1 *frame = (TH1*) gPad->DrawFrame(xmin,ymin,xmax,ymax);
  frame->SetXTitle("H(0^{+})/H(0^{-})");
  frame->SetYTitle("Entries");
  frame->Draw();

  althist->Draw("hist sames");
  nullhist->Draw("hist sames");

  // Set shade distribution
  double altMedian = median(althist); //CalcResult->GetTestStatisticData();
  double nullMedian = median(nullhist); //CalcResult->GetTestStatisticData();

  TH1F *altshade = (TH1F*)althist->Clone();
  TH1F *nullshade = (TH1F*)nullhist->Clone();
  TString altshadename = althist->GetName(); altshadename += "_shaded";
  TString nullshadename = nullhist->GetName(); nullshadename += "_shaded";
  altshade->SetName(altshadename);
  altshade->SetFillStyle(3004);
  altshade->SetFillColor(kBlue);
  altshade->SetLineWidth(0);
  nullshade->SetName(nullshadename);
  nullshade->SetFillStyle(3005);
  nullshade->SetFillColor(kRed);
  nullshade->SetLineWidth(0);
  altshade->Draw("hist sames");
  nullshade->Draw("hist sames");

  int altbin = althist->GetXaxis()->FindBin(altMedian);
  int nullbin = nullhist->GetXaxis()->FindBin(nullMedian);
  for (int i=0; i<altshade->GetNbinsX(); ++i) {
    if (i > altbin)  altshade->SetBinContent(i,0);
    if (i < altbin)  nullshade->SetBinContent(i,0);
  }

  // Set Asimov Line
  TLine *altline = new TLine(altMedian, ymin, altMedian, althist->GetBinContent(altbin));
  altline->SetLineWidth(3);
  altline->SetLineColor(kBlack);
  altline->Draw();

  TLine *nullline = new TLine(nullMedian, ymin, nullMedian, nullhist->GetBinContent(nullbin));
  nullline->SetLineWidth(2);
  nullline->SetLineStyle(9);
  nullline->SetLineColor(kRed);
  nullline->Draw();

  TLegend *leg = new TLegend(0.6,0.7,0.9,0.9);
  leg->AddEntry(ModelConfig,"J^{P} = 0^{+} (H^{0})","L");
  leg->AddEntry(ModelConfigNull,"J^{P} = 0^{-} (A^{0})","L");
  leg->AddEntry(altline,"Asimov Data","L");
  leg->SetTextSize(0.035);
  leg->SetLineColor(0);
  leg->SetFillColor(0);
  leg->Draw();

  ATLAS_WorkInProgress(0.19,0.87,0.035);
  myText(0.19,0.8,0.035,"#int L dt = 20.3fb^{-1} #sqrt{s} = 8TeV");
  myText(0.19,0.73,0.035,"H#rightarrow#tau#tau");

}