rs401d_FeldmanCousins.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roostats
/// \notebook
/// Neutrino Oscillation Example from Feldman & Cousins
///
/// This tutorial shows a more complex example using the FeldmanCousins utility
/// to create a confidence interval for a toy neutrino oscillation experiment.
/// The example attempts to faithfully reproduce the toy example described in Feldman & Cousins'
/// original paper, Phys.Rev.D57:3873-3889,1998.
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Kyle Cranmer

#include "RooGlobalFunc.h"
#include "RooStats/ConfInterval.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/MCMCCalculator.h"
#include "RooStats/UniformProposal.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/MCMCIntervalPlot.h"
#include "RooStats/MCMCInterval.h"

#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"

#include "TROOT.h"
#include "RooPolynomial.h"
#include "RooRandom.h"

#include "RooNLLVar.h"
#include "RooProfileLL.h"

#include "RooPlot.h"

#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TTree.h"
#include "TMarker.h"
#include "TStopwatch.h"

#include <iostream>

// PDF class created for this macro
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "../tutorials/roostats/NuMuToNuE_Oscillation.h"
#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx" // so that it can be executed directly
#else
#include "../tutorials/roostats/NuMuToNuE_Oscillation.cxx+" // so that it can be executed directly
#endif

// use this order for safety on library loading
using namespace RooFit;
using namespace RooStats;

void rs401d_FeldmanCousins(bool doFeldmanCousins = false, bool doMCMC = true)
{

   // to time the macro
   TStopwatch t;
   t.Start();

   // Taken from Feldman & Cousins paper, Phys.Rev.D57:3873-3889,1998.
   // e-Print: physics/9711021 (see page 13.)
   //
   // Quantum mechanics dictates that the probability of such a transformation is given by the formula
   // $P (\nu\mu \rightarrow \nu e ) = sin^2 (2\theta) sin^2 (1.27 \Delta m^2 L /E )$
   // where P is the probability for a $\nu\mu$ to transform into a $\nu e$ , L is the distance in km between
   // the creation of the neutrino from meson decay and its interaction in the detector, E is the
   // neutrino energy in GeV, and $\Delta m^2 = |m^2 - m^2 |$ in $(eV/c^2 )^2$ .
   //
   // To demonstrate how this works in practice, and how it compares to alternative approaches
   // that have been used, we consider a toy model of a typical neutrino oscillation experiment.
   // The toy model is defined by the following parameters: Mesons are assumed to decay to
   // neutrinos uniformly in a region 600 m to 1000 m from the detector. The expected background
   // from conventional $\nu e$ interactions and misidentified $\nu\mu$ interactions is assumed to be 100
   // events in each of 5 energy bins which span the region from 10 to 60 GeV. We assume that
   // the $\nu\mu$ flux is such that if $P (\nu\mu \rightarrow \nu e ) = 0.01$ averaged over any bin, then that bin
   // would
   // have an expected additional contribution of 100 events due to $\nu\mu \rightarrow \nu e$ oscillations.

   // Make signal model model
   RooRealVar E("E", "", 15, 10, 60, "GeV");
   RooRealVar L("L", "", .800, .600, 1.0, "km"); // need these units in formula
   RooRealVar deltaMSq("deltaMSq", "#Delta m^{2}", 40, 1, 300, "eV/c^{2}");
   RooRealVar sinSq2theta("sinSq2theta", "sin^{2}(2#theta)", .006, .0, .02);
   // RooRealVar deltaMSq("deltaMSq","#Delta m^{2}",40,20,70,"eV/c^{2}");
   //  RooRealVar sinSq2theta("sinSq2theta","sin^{2}(2#theta)", .006,.001,.01);
   // PDF for oscillation only describes deltaMSq dependence, sinSq2theta goes into sigNorm
   // 1) The code for this PDF was created by issuing these commands
   //    root [0] RooClassFactory x
   //    root [1] x.makePdf("NuMuToNuE_Oscillation","L,E,deltaMSq","","pow(sin(1.27*deltaMSq*L/E),2)")
   NuMuToNuE_Oscillation PnmuTone("PnmuTone", "P(#nu_{#mu} #rightarrow #nu_{e}", L, E, deltaMSq);

   // only E is observable, so create the signal model by integrating out L
   RooAbsPdf *sigModel = PnmuTone.createProjection(L);

   // create  $ \int dE' dL' P(E',L' | \Delta m^2)$.
   // Given RooFit will renormalize the PDF in the range of the observables,
   // the average probability to oscillate in the experiment's acceptance
   // needs to be incorporated into the extended term in the likelihood.
   // Do this by creating a RooAbsReal representing the integral and divide by
   // the area in the E-L plane.
   // The integral should be over "primed" observables, so we need
   // an independent copy of PnmuTone not to interfere with the original.

   // Independent copy for Integral
   RooRealVar EPrime("EPrime", "", 15, 10, 60, "GeV");
   RooRealVar LPrime("LPrime", "", .800, .600, 1.0, "km"); // need these units in formula
   NuMuToNuE_Oscillation PnmuTonePrime("PnmuTonePrime", "P(#nu_{#mu} #rightarrow #nu_{e}", LPrime, EPrime, deltaMSq);
   RooAbsReal *intProbToOscInExp = PnmuTonePrime.createIntegral(RooArgSet(EPrime, LPrime));

   // Getting the flux is a bit tricky.  It is more clear to include a cross section term that is not
   // explicitly referred to in the text, eg.
   // number events in bin = flux * cross-section for nu_e interaction in E bin * average prob nu_mu osc. to nu_e in bin
   // let maxEventsInBin = flux * cross-section for nu_e interaction in E bin
   // maxEventsInBin * 1% chance per bin =  100 events / bin
   // therefore maxEventsInBin = 10,000.
   // for 5 bins, this means maxEventsTot = 50,000
   RooConstVar maxEventsTot("maxEventsTot", "maximum number of sinal events", 50000);
   RooConstVar inverseArea("inverseArea", "1/(#Delta E #Delta L)",
                           1. / (EPrime.getMax() - EPrime.getMin()) / (LPrime.getMax() - LPrime.getMin()));

   // $sigNorm = maxEventsTot \cdot \int dE dL \frac{P_{oscillate\ in\ experiment}}{Area} \cdot {sin}^2(2\theta)$
   RooProduct sigNorm("sigNorm", "", RooArgSet(maxEventsTot, *intProbToOscInExp, inverseArea, sinSq2theta));
   // bkg = 5 bins * 100 events / bin
   RooConstVar bkgNorm("bkgNorm", "normalization for background", 500);

   // flat background (0th order polynomial, so no arguments for coefficients)
   RooPolynomial bkgEShape("bkgEShape", "flat bkg shape", E);

   // total model
   RooAddPdf model("model", "", RooArgList(*sigModel, bkgEShape), RooArgList(sigNorm, bkgNorm));

   // for debugging, check model tree
   //  model.printCompactTree();
   //  model.graphVizTree("model.dot");

   // turn off some messages
   RooMsgService::instance().setStreamStatus(0, kFALSE);
   RooMsgService::instance().setStreamStatus(1, kFALSE);
   RooMsgService::instance().setStreamStatus(2, kFALSE);

   // --------------------------------------
   // n events in data to data, simply sum of sig+bkg
   Int_t nEventsData = bkgNorm.getVal() + sigNorm.getVal();
   cout << "generate toy data with nEvents = " << nEventsData << endl;
   // adjust random seed to get a toy dataset similar to one in paper.
   // Found by trial and error (3 trials, so not very "fine tuned")
   RooRandom::randomGenerator()->SetSeed(3);
   // create a toy dataset
   RooDataSet *data = model.generate(RooArgSet(E), nEventsData);

   // --------------------------------------
   // make some plots
   TCanvas *dataCanvas = new TCanvas("dataCanvas");
   dataCanvas->Divide(2, 2);

   // plot the PDF
   dataCanvas->cd(1);
   TH1 *hh = PnmuTone.createHistogram("hh", E, Binning(40), YVar(L, Binning(40)), Scaling(kFALSE));
   hh->SetLineColor(kBlue);
   hh->SetTitle("True Signal Model");
   hh->Draw("surf");

   // plot the data with the best fit
   dataCanvas->cd(2);
   RooPlot *Eframe = E.frame();
   data->plotOn(Eframe);
   model.fitTo(*data, Extended());
   model.plotOn(Eframe);
   model.plotOn(Eframe, Components(*sigModel), LineColor(kRed));
   model.plotOn(Eframe, Components(bkgEShape), LineColor(kGreen));
   model.plotOn(Eframe);
   Eframe->SetTitle("toy data with best fit model (and sig+bkg components)");
   Eframe->Draw();

   // plot the likelihood function
   dataCanvas->cd(3);
   RooNLLVar nll("nll", "nll", model, *data, Extended());
   RooProfileLL pll("pll", "", nll, RooArgSet(deltaMSq, sinSq2theta));
   //  TH1* hhh = nll.createHistogram("hhh",sinSq2theta,Binning(40),YVar(deltaMSq,Binning(40))) ;
   TH1 *hhh = pll.createHistogram("hhh", sinSq2theta, Binning(40), YVar(deltaMSq, Binning(40)), Scaling(kFALSE));
   hhh->SetLineColor(kBlue);
   hhh->SetTitle("Likelihood Function");
   hhh->Draw("surf");

   dataCanvas->Update();

   // --------------------------------------------------------------
   // show use of Feldman-Cousins utility in RooStats
   // set the distribution creator, which encodes the test statistic
   RooArgSet parameters(deltaMSq, sinSq2theta);
   RooWorkspace *w = new RooWorkspace();

   ModelConfig modelConfig;
   modelConfig.SetWorkspace(*w);
   modelConfig.SetPdf(model);
   modelConfig.SetParametersOfInterest(parameters);

   RooStats::FeldmanCousins fc(*data, modelConfig);
   fc.SetTestSize(.1); // set size of test
   fc.UseAdaptiveSampling(true);
   fc.SetNBins(10); // number of points to test per parameter

   // use the Feldman-Cousins tool
   ConfInterval *interval = 0;
   if (doFeldmanCousins)
      interval = fc.GetInterval();

   // ---------------------------------------------------------
   // show use of ProfileLikeihoodCalculator utility in RooStats
   RooStats::ProfileLikelihoodCalculator plc(*data, modelConfig);
   plc.SetTestSize(.1);

   ConfInterval *plcInterval = plc.GetInterval();

   // --------------------------------------------
   // show use of MCMCCalculator utility in RooStats
   MCMCInterval *mcInt = NULL;

   if (doMCMC) {
      // turn some messages back on
      RooMsgService::instance().setStreamStatus(0, kTRUE);
      RooMsgService::instance().setStreamStatus(1, kTRUE);

      TStopwatch mcmcWatch;
      mcmcWatch.Start();

      RooArgList axisList(deltaMSq, sinSq2theta);
      MCMCCalculator mc(*data, modelConfig);
      mc.SetNumIters(5000);
      mc.SetNumBurnInSteps(100);
      mc.SetUseKeys(true);
      mc.SetTestSize(.1);
      mc.SetAxes(axisList); // set which is x and y axis in posterior histogram
      // mc.SetNumBins(50);
      mcInt = (MCMCInterval *)mc.GetInterval();

      mcmcWatch.Stop();
      mcmcWatch.Print();
   }
   // -------------------------------
   // make plot of resulting interval

   dataCanvas->cd(4);

   // first plot a small dot for every point tested
   if (doFeldmanCousins) {
      RooDataHist *parameterScan = (RooDataHist *)fc.GetPointsToScan();
      TH2F *hist = (TH2F *)parameterScan->createHistogram("sinSq2theta:deltaMSq", 30, 30);
      //  hist->Draw();
      TH2F *forContour = (TH2F *)hist->Clone();

      // now loop through the points and put a marker if it's in the interval
      RooArgSet *tmpPoint;
      // loop over points to test
      for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
         // get a parameter point from the list of points to test.
         tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");

         if (interval) {
            if (interval->IsInInterval(*tmpPoint)) {
               forContour->SetBinContent(
                  hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 1);
            } else {
               forContour->SetBinContent(
                  hist->FindBin(tmpPoint->getRealValue("sinSq2theta"), tmpPoint->getRealValue("deltaMSq")), 0);
            }
         }

         delete tmpPoint;
      }

      if (interval) {
         Double_t level = 0.5;
         forContour->SetContour(1, &level);
         forContour->SetLineWidth(2);
         forContour->SetLineColor(kRed);
         forContour->Draw("cont2,same");
      }
   }

   MCMCIntervalPlot *mcPlot = NULL;
   if (mcInt) {
      cout << "MCMC actual confidence level: " << mcInt->GetActualConfidenceLevel() << endl;
      mcPlot = new MCMCIntervalPlot(*mcInt);
      mcPlot->SetLineColor(kMagenta);
      mcPlot->Draw();
   }
   dataCanvas->Update();

   LikelihoodIntervalPlot plotInt((LikelihoodInterval *)plcInterval);
   plotInt.SetTitle("90% Confidence Intervals");
   if (mcInt)
      plotInt.Draw("same");
   else
      plotInt.Draw();
   dataCanvas->Update();

   /// print timing info
   t.Stop();
   t.Print();
}