Profile likelihood confidence limits

#include "TROOT.h"
#include "TSystem.h"
#include "TRolke.h"
#include "Riostream.h"
   
   
void Rolke()
{
//////////////////////////////////////////////////
// Routine computes the profile likelihood confidence
// limits for  7 different model assumptions
// on systematic/statistical uncertainties
//
// You must load libPhysics before executing this script.
//
// Author : Jan Conrad (CERN) <Jan.Conrad@cern.ch>
//////////////////////////////////////////////////



/////////////////////////////////////////////////////////////
// Model 1 assumes:
//
// Poisson uncertainty in the background estimate
// Binomial uncertainty in the efficiency estimate
//
// y = 10      events observed in the background region
// x = 5       events in the signal region
// tau = 2.5   ratio between size of signal/background region
// m = 100     MC events have been produced  (signal)
// z = 50      MC events have been observed (signal)
// alpha = 0.9 Confidence Level
//////////////////////////////////////////////////////////////

    //gSystem->Load("libPhysics");
 Double_t bm = 0.0;
 Double_t tau = 2.5;
 Int_t mid = 1;
 Int_t m = 100;
 Int_t z = 50;
 Int_t y = 10;
 Int_t x = 5;
 // Initialize parameters not used.
 Double_t e = 0.0;
 Double_t em = 0.0;
 Double_t sde=0.0;
 Double_t sdb=0.0;
 Double_t b = 0.0;


 TRolke g;
 
 g.SetCL(0.90);
 
 Double_t ul = g.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
 Double_t ll = g.GetLowerLimit();
 
 cout << "Assuming Model 1" << endl; 
 cout <<  "the Profile Likelihood interval is :" << endl;
 cout << "[" << ll << "," << ul << "]" << endl;

 
/////////////////////////////////////////////////////////////
// Model 2 assumes:
//
// Poisson uncertainty in the background estimate
// Gaussian  uncertainty in the efficiency estimate
//
// y = 3        events observed in the background region
// x = 10       events in the signal region
// tau = 2.5    ratio between size of signal/background region
// em = 0.9     measured efficiency
// sde = 0.05   standard deviation of efficiency
// alpha =0.95  Confidence L evel
//////////////////////////////////////////////////////////////


 tau = 2.5;
 mid = 2;
 y = 3;
 x = 10;
 em=0.9;
 sde=0.05;

 g.SetCL(0.95);

 ul = g.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
 ll = g.GetLowerLimit();
 
  cout << "Assuming MODEL 2" << endl; 
  cout <<  "the Profile Likelihood interval is :" << endl;
  cout << "[" << ll << "," << ul << "]" << endl;
  

/////////////////////////////////////////////////////////////
// Model 3 assumes:
//
// Gaussian uncertainty in the background estimate
// Gaussian  uncertainty in the efficiency estimate
//
// bm = 5       expected background
// x = 10       events in the signal region
// sdb = 0.5    standard deviation in background estimate
// em = 0.9     measured efficiency
// sde = 0.05   standard deviation of efficiency
// alpha =0.99 Confidence Level
//////////////////////////////////////////////////////////////



 mid = 3;
 bm = 5.0;
 x = 10;
 em = 0.9;
 sde=0.05;
 sdb=0.5;

 g.SetCL(0.99);


 ul = g.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
 ll = g.GetLowerLimit();
 
cout << "Assuming Model 3" << endl; 
cout <<  "the Profile Likelihood interval is :" << endl;
cout << "[" << ll << "," << ul << "]" << endl;





/////////////////////////////////////////////////////////////
// Model 4 assumes:
//
// Poisson uncertainty in the background estimate
// known efficiency
//
// y = 7       events observed in the background region
// x = 1       events in the signal region
// tau = 5     ratio between size of signal/background region
//
// alpha =0.68  Confidence L evel
//////////////////////////////////////////////////////////////


 tau = 5;
 mid = 4;
 y = 7;
 x = 1;
 e = 0.25;


 g.SetCL(0.68);

 ul = g.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
 ll = g.GetLowerLimit();
 
  cout << "Assuming Model 4" << endl; 
    cout <<  "the Profile Likelihood interval is :" << endl;
  cout << "[" << ll << "," << ul << "]" << endl;





////////////////////////////////////////////////////////
// Model 5 assumes:
//
// Gaussian uncertainty in the background estimate
// Known efficiency
//
// bm = 0           measured background expectation
// x = 1            events in the signal region
// e = 0.65
// sdb = 1.         standard deviation of background estimate
// alpha =0.799999  Confidence Level
///////////////////////////////////////////////////////


 mid = 5;
 bm = 0.0;
 x = 1;
 e = 0.65;
 sdb=1.0;

 g.SetCL(0.80);

 ul = g.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
 ll = g.GetLowerLimit();
 
  cout << "Assuming Model 5" << endl; 
  cout <<  "the Profile Likelihood interval is :" << endl;
  cout << "[" << ll << "," << ul << "]" << endl;
  

////////////////////////////////////////////////////////
// Model 6 assumes:
//
// Known background 
// Binomial uncertainty in the efficiency estimate
//
// b = 10       known background
// x = 25       events in the signal region
// z = 500      Number of observed signal MC events
// m = 750      Number of produced MC signal events
// alpha =0.9   Confidence L evel
///////////////////////////////////////////////////////

 y = 1;
 mid = 6;
 m = 750;
 z = 500;
 x = 25;
 b = 10.0;

 TRolke p; 
 p.SetCL(0.90);
 ul = p.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
 Double_t newll = p.GetLowerLimit();
  
  cout << "Assuming Model 6" << endl; 
  cout <<  "the Profile Likelihood interval is :" << endl;
  cout << "[" << newll << "," << ul << "]" << endl;
  


////////////////////////////////////////////////////////
// Model 7 assumes:
//
// Known Background
// Gaussian  uncertainty in the efficiency estimate
//
// x = 15       events in the signal region
// em = 0.77    measured efficiency
// sde = 0.15   standard deviation of efficiency estimate
// b = 10       known background
// alpha =0.95  Confidence L evel
///////////////////////////////////////////////////////


 mid = 7;
 x = 15;
 em = 0.77;
 sde=0.15;
 b = 10.0;

 g.SetCL(0.95);
 
  ul = g.CalculateInterval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
  ll = g.GetLowerLimit();
  
  cout << "Assuming Model 7 " << endl; 
  cout <<  "the Profile Likelihood interval is :" << endl;
  cout << "[" << ll << "," << ul << "]" << endl;
}