```// @(#)root/roostats:\$Id\$
// Author: Kyle Cranmer   January 2009

/*************************************************************************
*                                                                       *
* For the licensing terms see \$ROOTSYS/LICENSE.                         *
* For the list of contributors see \$ROOTSYS/README/CREDITS.             *
*************************************************************************/

//_________________________________________________
/*
BEGIN_HTML
<p>
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration
of the more general NeymanConstruction.  It is a concrete implementation of the IntervalCalculator interface that, which uses the NeymanConstruction in a particular way.  As the name suggests, it returns a ConfidenceInterval.  In particular, it produces a RooStats::PointSetInterval, which is a concrete implementation of the ConfInterval interface.
</p>
<p>
The Neyman Construction is not a uniquely defined statistical technique, it requires that one specify an ordering rule
or ordering principle, which is usually incoded by choosing a specific test statistic and limits of integration
(corresponding to upper/lower/central limits).  As a result, this class must be configured with the corresponding
information before it can produce an interval.
</p>
<p>In the case of the Feldman-Cousins approach, the ordering principle is the likelihood ratio -- motivated
by the Neyman-Pearson lemma.  When nuisance parameters are involved, the profile likelihood ratio is the natural generalization.  One may either choose to perform the construction over the full space of the nuisance parameters, or restrict the nusiance parameters to their conditional MLE (eg. profiled values).
</p>
END_HTML
*/
//

#ifndef RooStats_FeldmanCousins
#include "RooStats/FeldmanCousins.h"
#endif

#ifndef RooStats_RooStatsUtils
#include "RooStats/RooStatsUtils.h"
#endif

#ifndef RooStats_PointSetInterval
#include "RooStats/PointSetInterval.h"
#endif

#include "RooStats/ModelConfig.h"

#include "RooStats/SamplingDistribution.h"

#include "RooStats/ProfileLikelihoodTestStat.h"
#include "RooStats/NeymanConstruction.h"
#include "RooStats/RooStatsUtils.h"

#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGlobalFunc.h"
#include "RooMsgService.h"
#include "TFile.h"
#include "TTree.h"

ClassImp(RooStats::FeldmanCousins) ;

using namespace RooFit;
using namespace RooStats;
using namespace std;

/*
//_______________________________________________________
FeldmanCousins::FeldmanCousins() :
//  fModel(NULL),
fData(0),
fTestStatSampler(0),
fPointsToTest(0),
fNbins(10),
fFluctuateData(true),
fDoProfileConstruction(true),
fSaveBeltToFile(false),
fCreateBelt(false)
{
// default constructor
//   fWS = new RooWorkspace("FeldmanCousinsWS");
//   fOwnsWorkspace = true;
//   fDataName = "";
//   fPdfName = "";
}
*/

//_______________________________________________________
FeldmanCousins::FeldmanCousins(RooAbsData& data, ModelConfig& model) :
fSize(0.05),
fModel(model),
fData(data),
fTestStatSampler(0),
fPointsToTest(0),
fPOIToTest(0),
fConfBelt(0),
fNbins(10),
fFluctuateData(true),
fDoProfileConstruction(true),
fSaveBeltToFile(false),
fCreateBelt(false)
{
// standard constructor
}

//_______________________________________________________
FeldmanCousins::~FeldmanCousins() {
// destructor
//if(fOwnsWorkspace && fWS) delete fWS;
if(fPointsToTest) delete fPointsToTest;
if(fPOIToTest) delete fPOIToTest;
if(fTestStatSampler) delete fTestStatSampler;
}

//_______________________________________________________
void FeldmanCousins::SetModel(const ModelConfig & model) {
// set the model
fModel = model;
}

//_______________________________________________________
TestStatSampler*  FeldmanCousins::GetTestStatSampler() const{
if(!fTestStatSampler)
this->CreateTestStatSampler();
return fTestStatSampler;
}

//_______________________________________________________
void FeldmanCousins::CreateTestStatSampler() const{
// specify the Test Statistic and create a ToyMC test statistic sampler

// use the profile likelihood ratio as the test statistic
ProfileLikelihoodTestStat* testStatistic = new ProfileLikelihoodTestStat(*fModel.GetPdf());

// create the ToyMC test statistic sampler
fTestStatSampler->SetParametersForTestStat(*fModel.GetParametersOfInterest() );
if(fModel.GetObservables())
fTestStatSampler->SetObservables(*fModel.GetObservables());
fTestStatSampler->SetPdf(*fModel.GetPdf());

ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point = " << (int) (fAdditionalNToysFactor*50./fSize) << endl;
} else{
ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point: adaptive" << endl;
}
if(fFluctuateData){
ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will fluctuate about  expectation" << endl;
} else{
ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will not fluctuate, will always be " << fData.numEntries() << endl;
fTestStatSampler->SetNEventsPerToy(fData.numEntries());
}
}

//_______________________________________________________
void FeldmanCousins::CreateParameterPoints() const{
// specify the parameter points to perform the construction.
// allow ability to profile on some nuisance paramters

// get ingredients
RooAbsPdf* pdf   = fModel.GetPdf();
if (!pdf ){
ooccoutE(&fModel,Generation) << "FeldmanCousins: ModelConfig has no PDF" << endl;
return;
}

// get list of all paramters
RooArgSet* parameters = new RooArgSet(*fModel.GetParametersOfInterest());
if(fModel.GetNuisanceParameters())

if( fModel.GetNuisanceParameters() && ! fModel.GetParametersOfInterest()->equals(*parameters) && fDoProfileConstruction) {
// if parameters include nuisance parameters, do profile construction
ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has nuisance parameters, will do profile construction" << endl;

// set nbins for the POI
TIter it2 = fModel.GetParametersOfInterest()->createIterator();
RooRealVar *myarg2;
while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
myarg2->setBins(fNbins);
}

// get dataset for POI scan
//     RooDataHist* parameterScan = NULL;
RooAbsData* parameterScan = NULL;
if(fPOIToTest)
parameterScan = fPOIToTest;
else
parameterScan = new RooDataHist("parameterScan", "", *fModel.GetParametersOfInterest());

ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
// make profile construction
RooFit::MsgLevel previous  = RooMsgService::instance().globalKillBelow();
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL) ;
RooAbsReal* nll = pdf->createNLL(fData,RooFit::CloneData(false));
RooAbsReal* profile = nll->createProfile(*fModel.GetParametersOfInterest());

RooDataSet* profileConstructionPoints = new RooDataSet("profileConstruction",
"profileConstruction",
*parameters);

for(Int_t i=0; i<parameterScan->numEntries(); ++i){
// here's where we figure out the profiled value of nuisance parameters
*parameters = *parameterScan->get(i);
profile->getVal();
}
RooMsgService::instance().setGlobalKillBelow(previous) ;
delete profile;
delete nll;
if(!fPOIToTest) delete parameterScan;

// done
fPointsToTest = profileConstructionPoints;

} else{
// Do full construction
ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has no nuisance parameters" << endl;

TIter it = parameters->createIterator();
RooRealVar *myarg;
while ((myarg = dynamic_cast<RooRealVar*>(it.Next()))) {
myarg->setBins(fNbins);
}

RooDataHist* parameterScan = new RooDataHist("parameterScan", "", *parameters);
ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;

fPointsToTest = parameterScan;
}

delete parameters;

}

//_______________________________________________________
PointSetInterval* FeldmanCousins::GetInterval() const {
// Main interface to get a RooStats::ConfInterval.
// It constructs a RooStats::PointSetInterval.

// local variables
//  RooAbsData* data = fData; //fWS->data(fDataName);

// fill in implied variables given data
fModel.GuessObsAndNuisance(fData);

// create the test statistic sampler (private data member fTestStatSampler)
if(!fTestStatSampler)
this->CreateTestStatSampler();

fTestStatSampler->SetObservables(*fModel.GetObservables());

if(!fFluctuateData)
fTestStatSampler->SetNEventsPerToy(fData.numEntries());

// create paramter points to perform construction (private data member fPointsToTest)
this->CreateParameterPoints();

// Create a Neyman Construction
NeymanConstruction nc(fData,fModel);
// configure it
nc.SetTestStatSampler(*fTestStatSampler);
nc.SetTestSize(fSize); // set size of test
//  nc.SetParameters( fModel.GetParametersOfInterest);
nc.SetParameterPointsToTest( *fPointsToTest );
nc.SetLeftSideTailFraction(0.); // part of definition of Feldman-Cousins
nc.SetData(fData);
nc.SaveBeltToFile(fSaveBeltToFile);
nc.CreateConfBelt(fCreateBelt);
fConfBelt = nc.GetConfidenceBelt();
// use it
return nc.GetInterval();
}
```
FeldmanCousins.cxx:1
FeldmanCousins.cxx:2
FeldmanCousins.cxx:3
FeldmanCousins.cxx:4
FeldmanCousins.cxx:5
FeldmanCousins.cxx:6
FeldmanCousins.cxx:7
FeldmanCousins.cxx:8
FeldmanCousins.cxx:9
FeldmanCousins.cxx:10
FeldmanCousins.cxx:11
FeldmanCousins.cxx:12
FeldmanCousins.cxx:13
FeldmanCousins.cxx:14
FeldmanCousins.cxx:15
FeldmanCousins.cxx:16
FeldmanCousins.cxx:17
FeldmanCousins.cxx:18
FeldmanCousins.cxx:19
FeldmanCousins.cxx:20
FeldmanCousins.cxx:21
FeldmanCousins.cxx:22
FeldmanCousins.cxx:23
FeldmanCousins.cxx:24
FeldmanCousins.cxx:25
FeldmanCousins.cxx:26
FeldmanCousins.cxx:27
FeldmanCousins.cxx:28
FeldmanCousins.cxx:29
FeldmanCousins.cxx:30
FeldmanCousins.cxx:31
FeldmanCousins.cxx:32
FeldmanCousins.cxx:33
FeldmanCousins.cxx:34
FeldmanCousins.cxx:35
FeldmanCousins.cxx:36
FeldmanCousins.cxx:37
FeldmanCousins.cxx:38
FeldmanCousins.cxx:39
FeldmanCousins.cxx:40
FeldmanCousins.cxx:41
FeldmanCousins.cxx:42
FeldmanCousins.cxx:43
FeldmanCousins.cxx:44
FeldmanCousins.cxx:45
FeldmanCousins.cxx:46
FeldmanCousins.cxx:47
FeldmanCousins.cxx:48
FeldmanCousins.cxx:49
FeldmanCousins.cxx:50
FeldmanCousins.cxx:51
FeldmanCousins.cxx:52
FeldmanCousins.cxx:53
FeldmanCousins.cxx:54
FeldmanCousins.cxx:55
FeldmanCousins.cxx:56
FeldmanCousins.cxx:57
FeldmanCousins.cxx:58
FeldmanCousins.cxx:59
FeldmanCousins.cxx:60
FeldmanCousins.cxx:61
FeldmanCousins.cxx:62
FeldmanCousins.cxx:63
FeldmanCousins.cxx:64
FeldmanCousins.cxx:65
FeldmanCousins.cxx:66
FeldmanCousins.cxx:67
FeldmanCousins.cxx:68
FeldmanCousins.cxx:69
FeldmanCousins.cxx:70
FeldmanCousins.cxx:71
FeldmanCousins.cxx:72
FeldmanCousins.cxx:73
FeldmanCousins.cxx:74
FeldmanCousins.cxx:75
FeldmanCousins.cxx:76
FeldmanCousins.cxx:77
FeldmanCousins.cxx:78
FeldmanCousins.cxx:79
FeldmanCousins.cxx:80
FeldmanCousins.cxx:81
FeldmanCousins.cxx:82
FeldmanCousins.cxx:83
FeldmanCousins.cxx:84
FeldmanCousins.cxx:85
FeldmanCousins.cxx:86
FeldmanCousins.cxx:87
FeldmanCousins.cxx:88
FeldmanCousins.cxx:89
FeldmanCousins.cxx:90
FeldmanCousins.cxx:91
FeldmanCousins.cxx:92
FeldmanCousins.cxx:93
FeldmanCousins.cxx:94
FeldmanCousins.cxx:95
FeldmanCousins.cxx:96
FeldmanCousins.cxx:97
FeldmanCousins.cxx:98
FeldmanCousins.cxx:99
FeldmanCousins.cxx:100
FeldmanCousins.cxx:101
FeldmanCousins.cxx:102
FeldmanCousins.cxx:103
FeldmanCousins.cxx:104
FeldmanCousins.cxx:105
FeldmanCousins.cxx:106
FeldmanCousins.cxx:107
FeldmanCousins.cxx:108
FeldmanCousins.cxx:109
FeldmanCousins.cxx:110
FeldmanCousins.cxx:111
FeldmanCousins.cxx:112
FeldmanCousins.cxx:113
FeldmanCousins.cxx:114
FeldmanCousins.cxx:115
FeldmanCousins.cxx:116
FeldmanCousins.cxx:117
FeldmanCousins.cxx:118
FeldmanCousins.cxx:119
FeldmanCousins.cxx:120
FeldmanCousins.cxx:121
FeldmanCousins.cxx:122
FeldmanCousins.cxx:123
FeldmanCousins.cxx:124
FeldmanCousins.cxx:125
FeldmanCousins.cxx:126
FeldmanCousins.cxx:127
FeldmanCousins.cxx:128
FeldmanCousins.cxx:129
FeldmanCousins.cxx:130
FeldmanCousins.cxx:131
FeldmanCousins.cxx:132
FeldmanCousins.cxx:133
FeldmanCousins.cxx:134
FeldmanCousins.cxx:135
FeldmanCousins.cxx:136
FeldmanCousins.cxx:137
FeldmanCousins.cxx:138
FeldmanCousins.cxx:139
FeldmanCousins.cxx:140
FeldmanCousins.cxx:141
FeldmanCousins.cxx:142
FeldmanCousins.cxx:143
FeldmanCousins.cxx:144
FeldmanCousins.cxx:145
FeldmanCousins.cxx:146
FeldmanCousins.cxx:147
FeldmanCousins.cxx:148
FeldmanCousins.cxx:149
FeldmanCousins.cxx:150
FeldmanCousins.cxx:151
FeldmanCousins.cxx:152
FeldmanCousins.cxx:153
FeldmanCousins.cxx:154
FeldmanCousins.cxx:155
FeldmanCousins.cxx:156
FeldmanCousins.cxx:157
FeldmanCousins.cxx:158
FeldmanCousins.cxx:159
FeldmanCousins.cxx:160
FeldmanCousins.cxx:161
FeldmanCousins.cxx:162
FeldmanCousins.cxx:163
FeldmanCousins.cxx:164
FeldmanCousins.cxx:165
FeldmanCousins.cxx:166
FeldmanCousins.cxx:167
FeldmanCousins.cxx:168
FeldmanCousins.cxx:169
FeldmanCousins.cxx:170
FeldmanCousins.cxx:171
FeldmanCousins.cxx:172
FeldmanCousins.cxx:173
FeldmanCousins.cxx:174
FeldmanCousins.cxx:175
FeldmanCousins.cxx:176
FeldmanCousins.cxx:177
FeldmanCousins.cxx:178
FeldmanCousins.cxx:179
FeldmanCousins.cxx:180
FeldmanCousins.cxx:181
FeldmanCousins.cxx:182
FeldmanCousins.cxx:183
FeldmanCousins.cxx:184
FeldmanCousins.cxx:185
FeldmanCousins.cxx:186
FeldmanCousins.cxx:187
FeldmanCousins.cxx:188
FeldmanCousins.cxx:189
FeldmanCousins.cxx:190
FeldmanCousins.cxx:191
FeldmanCousins.cxx:192
FeldmanCousins.cxx:193
FeldmanCousins.cxx:194
FeldmanCousins.cxx:195
FeldmanCousins.cxx:196
FeldmanCousins.cxx:197
FeldmanCousins.cxx:198
FeldmanCousins.cxx:199
FeldmanCousins.cxx:200
FeldmanCousins.cxx:201
FeldmanCousins.cxx:202
FeldmanCousins.cxx:203
FeldmanCousins.cxx:204
FeldmanCousins.cxx:205
FeldmanCousins.cxx:206
FeldmanCousins.cxx:207
FeldmanCousins.cxx:208
FeldmanCousins.cxx:209
FeldmanCousins.cxx:210
FeldmanCousins.cxx:211
FeldmanCousins.cxx:212
FeldmanCousins.cxx:213
FeldmanCousins.cxx:214
FeldmanCousins.cxx:215
FeldmanCousins.cxx:216
FeldmanCousins.cxx:217
FeldmanCousins.cxx:218
FeldmanCousins.cxx:219
FeldmanCousins.cxx:220
FeldmanCousins.cxx:221
FeldmanCousins.cxx:222
FeldmanCousins.cxx:223
FeldmanCousins.cxx:224
FeldmanCousins.cxx:225
FeldmanCousins.cxx:226
FeldmanCousins.cxx:227
FeldmanCousins.cxx:228
FeldmanCousins.cxx:229
FeldmanCousins.cxx:230
FeldmanCousins.cxx:231
FeldmanCousins.cxx:232
FeldmanCousins.cxx:233
FeldmanCousins.cxx:234
FeldmanCousins.cxx:235
FeldmanCousins.cxx:236
FeldmanCousins.cxx:237
FeldmanCousins.cxx:238
FeldmanCousins.cxx:239
FeldmanCousins.cxx:240
FeldmanCousins.cxx:241
FeldmanCousins.cxx:242
FeldmanCousins.cxx:243
FeldmanCousins.cxx:244
FeldmanCousins.cxx:245
FeldmanCousins.cxx:246
FeldmanCousins.cxx:247
FeldmanCousins.cxx:248
FeldmanCousins.cxx:249
FeldmanCousins.cxx:250
FeldmanCousins.cxx:251
FeldmanCousins.cxx:252
FeldmanCousins.cxx:253
FeldmanCousins.cxx:254
FeldmanCousins.cxx:255
FeldmanCousins.cxx:256
FeldmanCousins.cxx:257
FeldmanCousins.cxx:258
FeldmanCousins.cxx:259
FeldmanCousins.cxx:260
FeldmanCousins.cxx:261
FeldmanCousins.cxx:262
FeldmanCousins.cxx:263
FeldmanCousins.cxx:264
FeldmanCousins.cxx:265
FeldmanCousins.cxx:266
FeldmanCousins.cxx:267
FeldmanCousins.cxx:268
FeldmanCousins.cxx:269
FeldmanCousins.cxx:270
FeldmanCousins.cxx:271
FeldmanCousins.cxx:272
FeldmanCousins.cxx:273
FeldmanCousins.cxx:274
FeldmanCousins.cxx:275
FeldmanCousins.cxx:276
FeldmanCousins.cxx:277
FeldmanCousins.cxx:278
FeldmanCousins.cxx:279
FeldmanCousins.cxx:280
FeldmanCousins.cxx:281
FeldmanCousins.cxx:282
FeldmanCousins.cxx:283