| 1 |
#include "RooPlot.h"
|
| 2 |
#include "RooRandom.h"
|
| 3 |
#include "RooRealVar.h"
|
| 4 |
#include "RooGaussian.h"
|
| 5 |
#include "RooPolynomial.h"
|
| 6 |
#include "RooArgSet.h"
|
| 7 |
#include "RooAddPdf.h"
|
| 8 |
#include "RooDataSet.h"
|
| 9 |
#include "RooExtendPdf.h"
|
| 10 |
#include "RooConstVar.h"
|
| 11 |
#include "RooMsgService.h"
|
| 12 |
|
| 13 |
#ifndef __CINT__ // problem including this file with CINT
|
| 14 |
#include "RooGlobalFunc.h"
|
| 15 |
#endif
|
| 16 |
|
| 17 |
#include "RooStats/FeldmanCousins.h"
|
| 18 |
#include "RooStats/PointSetInterval.h"
|
| 19 |
|
| 20 |
#include "TCanvas.h"
|
| 21 |
#include "TH1F.h"
|
| 22 |
#include "TFile.h"
|
| 23 |
#include "TChain.h"
|
| 24 |
#include "TGraphErrors.h"
|
| 25 |
#include "TF1.h"
|
| 26 |
|
| 27 |
#include "RooWorkspace.h"
|
| 28 |
|
| 29 |
#include "./_tdrstyle.C"
|
| 30 |
|
| 31 |
#include "RooStats/RooStatsUtils.h"
|
| 32 |
#include "RooStats/ProfileLikelihoodTestStat.h"
|
| 33 |
|
| 34 |
void median(int ntoys = 10000){
|
| 35 |
using namespace RooFit;
|
| 36 |
using namespace RooStats;
|
| 37 |
setTDRStyle();
|
| 38 |
RooRandom::randomGenerator()->SetSeed(10);
|
| 39 |
|
| 40 |
TFile file("wsGI05_20.root");
|
| 41 |
RooWorkspace *wspace = (RooWorkspace*)file.Get("myWS");
|
| 42 |
// ModelConfig *modelConfig = (ModelConfig*)wspace->genobj("sigBg");
|
| 43 |
|
| 44 |
RooAbsPdf* model = wspace->pdf("model");
|
| 45 |
RooAbsPdf* null = wspace->pdf("null");
|
| 46 |
RooRealVar* qT = wspace->var("qT");
|
| 47 |
RooRealVar* dy_slope = wspace->var("dy_slope");
|
| 48 |
RooRealVar* luminosity= wspace->var("luminosity");
|
| 49 |
RooRealVar* sig_yield = wspace->var("sig_yield");
|
| 50 |
|
| 51 |
*dy_slope = 0.6555;
|
| 52 |
*luminosity = 35.53;
|
| 53 |
RooArgSet POI(*sig_yield);
|
| 54 |
|
| 55 |
TCanvas *c1 = new TCanvas();
|
| 56 |
ProfileLikelihoodTestStat* testStatistic = new ProfileLikelihoodTestStat(*model);
|
| 57 |
TFile *f = new TFile("SamplingDistributions_gi05_35pb_NuisLumiSlope_20.root");
|
| 58 |
const int nBins=100;
|
| 59 |
TH1F *step[nBins];
|
| 60 |
TH1F *cl90 = new TH1F("cl90","",100,0,10), *cl95 = new TH1F("cl95","",100,0,10);
|
| 61 |
char buff[64];
|
| 62 |
for(int bin=0; bin<nBins; bin++){
|
| 63 |
sprintf(buff,"step_%d",bin+1);
|
| 64 |
step[bin] = new TH1F(buff,buff,500,0,5);
|
| 65 |
sprintf(buff,"hist_temp;%d",bin+1);
|
| 66 |
TH1 *tmp = (TH1*)f->Get(buff);
|
| 67 |
if( tmp==0 ) cout<<"No toys for "<<bin<<endl;
|
| 68 |
else {
|
| 69 |
sprintf(buff,"hist_temp_%d",bin+1);
|
| 70 |
tmp->SetName(buff);
|
| 71 |
step[bin]->Add(tmp);
|
| 72 |
//cout<<step[bin]->GetName()<<endl;
|
| 73 |
}
|
| 74 |
}
|
| 75 |
for(int iter=0; iter<100; iter++){
|
| 76 |
cout<<"Iteration: "<<iter<<endl;
|
| 77 |
RooRandom::randomGenerator()->SetSeed(iter+100);
|
| 78 |
*sig_yield = 0;
|
| 79 |
RooDataSet* data = null->generate(*qT,RooFit::Extended());
|
| 80 |
TGraphErrors *ge = new TGraphErrors(nBins);
|
| 81 |
// RooAbsReal* nll = model->createNLL(*data, RooFit::CloneData(kFALSE));
|
| 82 |
RooNLLVar *nll = new RooNLLVar("nll","",*model,*data, RooFit::Extended());
|
| 83 |
for(int bin=0; bin<nBins; bin++){
|
| 84 |
*sig_yield = (bin+0.5) * sig_yield->getMax()/float(nBins);
|
| 85 |
//cout<<"NLL("<<norm->getVal()<<")="<<nll->evaluate()<<" ";//endl;//((RooProfileLL*)testStatistic->GetTestStatistic())->bestFitParams().getRealValue("norm")<<endl;
|
| 86 |
cout<<"NLL("<<sig_yield->getVal()<<")="<<nll->getVal()<<" ";//endl;//((RooProfileLL*)testStatistic->GetTestStatistic())->bestFitParams().getRealValue("norm")<<endl;
|
| 87 |
Double_t thisTestStatistic = testStatistic->Evaluate(*data,POI);
|
| 88 |
int llr_bin = step[bin]->GetXaxis()->FindBin(thisTestStatistic);
|
| 89 |
double pValue = step[bin]->Integral(llr_bin,step[bin]->GetNbinsX())/step[bin]->Integral();
|
| 90 |
double pUncert = sqrt(step[bin]->Integral(llr_bin,step[bin]->GetNbinsX()))/step[bin]->Integral();
|
| 91 |
cout<<"TestStatistic at "<<(bin+1)<<"/"<<100<<"("<<sig_yield->getVal()<<") = "<<thisTestStatistic<<" pValue="<<pValue<<" nDataEvents="<<data->numEntries()<<endl;
|
| 92 |
//ge->SetPoint(bin,(bin+0.5)*10./nBins,pValue);
|
| 93 |
ge->SetPoint(bin,(bin+0.5) * sig_yield->getMax()/float(nBins) ,pValue);
|
| 94 |
ge->SetPointError(bin,0,pUncert);
|
| 95 |
}
|
| 96 |
//sprintf(buff,"m=500 GeV/c^{2} profile (iteration=%d)",iter);
|
| 97 |
sprintf(buff,"gi05_20_iteration_%d",iter);
|
| 98 |
ge->SetName(buff);
|
| 99 |
ge->Draw("AP");
|
| 100 |
sprintf(buff,"gi05_20_iteration_%d.eps",iter);
|
| 101 |
c1->Print(buff);
|
| 102 |
model->fitTo(*data);
|
| 103 |
cout<<"Best fit's yield="<<sig_yield->getVal()<<" LLR="<<testStatistic->Evaluate(*data,POI)<<" NLL="<<nll->evaluate()<<endl;
|
| 104 |
|
| 105 |
//*norm = 3.1;
|
| 106 |
//cout<<"For norm="<<norm->getVal()<<" LLR="<<testStatistic->Evaluate(*data,POI)<<" NLL="<<nll->evaluate()<<endl;
|
| 107 |
//RooPlot* frame = qT->frame();
|
| 108 |
//data->plotOn(frame);
|
| 109 |
//model->plotOn(frame);
|
| 110 |
//frame->Draw();
|
| 111 |
//sprintf(buff,"./tmp/data_%d.png",iter);
|
| 112 |
//c1->Print(buff);
|
| 113 |
}
|
| 114 |
}
|
| 115 |
|
