| 53 |
|
|
| 54 |
|
|
| 55 |
|
double simpleProfile2(ConfigFile * config, string Type, |
| 56 |
< |
int dat, double bkg, double bkg_uncertainty, double sig, double |
| 56 |
> |
double dat, double bkg, double bkg_uncertainty, double sig, double |
| 57 |
|
sig_uncertainty, double xsec, double ExpNsigLimit ) //absolute uncertainties! |
| 58 |
|
{ |
| 59 |
|
TStopwatch t; |
| 64 |
|
///////////////////////////////////////// |
| 65 |
|
char * name = new char[1024]; |
| 66 |
|
RooWorkspace* wspace = new RooWorkspace(); |
| 67 |
< |
sprintf(name,"Poisson::countingModel(obs[150,0,300],sum(s[%d,0,%d]*ratioSigEff[1.,0,2.],b[100,0,300]*ratioBkgEff[1.,0.,2.]))", (int)ExpNsigLimit,(int)(ExpNsigLimit*2.5)); |
| 67 |
> |
sprintf(name,"Poisson::countingModel(obs[50,0,200],sum(s[50,0,100]*ratioSigEff[1.,0,2.],b[50,0,100]*ratioBkgEff[1.,0.,2.]))"); |
| 68 |
|
wspace->factory(name); // counting model |
| 69 |
|
//wspace->factory("Gaussian::sigConstraint(1,ratioSigEff,0.257)"); // 5% signal efficiency uncertainty |
| 70 |
|
//wspace->factory("Gaussian::bkgConstraint(1,ratioBkgEff,0.136)"); // 10% background efficiency uncertainty |
| 71 |
|
sprintf(name,"Gaussian::sigConstraint(1,ratioSigEff,%f)",sig_uncertainty/sig); |
| 72 |
|
wspace->factory(name); // 5% signal efficiency uncertainty |
| 73 |
< |
sprintf(name,"Gaussian::bkgConstraint(1,ratioBkgEff,%f)",sig_uncertainty/sig); |
| 73 |
> |
sprintf(name,"Gaussian::bkgConstraint(1,ratioBkgEff,%f)",bkg_uncertainty/bkg); |
| 74 |
|
wspace->factory(name); // 10% background efficiency uncertainty |
| 75 |
|
wspace->factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)"); // product of terms |
| 76 |
|
wspace->Print(); |
| 107 |
|
// First, let's use a Calculator based on the Profile Likelihood Ratio |
| 108 |
|
//ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest); |
| 109 |
|
ProfileLikelihoodCalculator plc(*data, modelConfig); |
| 110 |
< |
plc.SetTestSize(.05); |
| 110 |
> |
plc.SetTestSize(.10); |
| 111 |
|
ConfInterval* lrint = plc.GetInterval(); // that was easy. |
| 112 |
|
|
| 113 |
|
// Let's make a plot |
| 131 |
|
plotInt.Draw(); |
| 132 |
|
|
| 133 |
|
// Second, use a Calculator based on the Feldman Cousins technique |
| 134 |
< |
FeldmanCousins fc(*data, modelConfig); |
| 135 |
< |
fc.UseAdaptiveSampling(true); |
| 136 |
< |
fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed |
| 137 |
< |
fc.SetNBins(100); // number of points to test per parameter |
| 138 |
< |
fc.SetTestSize(.05); |
| 139 |
< |
// fc.SaveBeltToFile(true); // optional |
| 140 |
< |
ConfInterval* fcint = NULL; |
| 141 |
< |
fcint = fc.GetInterval(); // that was easy. |
| 134 |
> |
// FeldmanCousins fc(*data, modelConfig); |
| 135 |
> |
// fc.UseAdaptiveSampling(true); |
| 136 |
> |
// fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed |
| 137 |
> |
// fc.SetNBins(100); // number of points to test per parameter |
| 138 |
> |
// fc.SetTestSize(.10); //95% single sided |
| 139 |
> |
// fc.AdditionalNToysFactor(5); |
| 140 |
> |
// // fc.SaveBeltToFile(true); // optional |
| 141 |
> |
// ConfInterval* fcint = NULL; |
| 142 |
> |
// fcint = fc.GetInterval(); // that was easy. |
| 143 |
|
|
| 144 |
|
RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true)); |
| 145 |
|
|
| 154 |
|
ProposalFunction* pdfProp = ph.GetProposalFunction(); // that was easy |
| 155 |
|
|
| 156 |
|
MCMCCalculator mc(*data, modelConfig); |
| 157 |
< |
mc.SetNumIters(20000); // steps to propose in the chain |
| 158 |
< |
mc.SetTestSize(.05); // 95% CL |
| 159 |
< |
mc.SetNumBurnInSteps(40); // ignore first N steps in chain as "burn in" |
| 157 |
> |
mc.SetNumIters(100000); // steps to propose in the chain |
| 158 |
> |
mc.SetTestSize(.10); // 95% CL single-sided |
| 159 |
> |
mc.SetNumBurnInSteps(100); // ignore first N steps in chain as "burn in" |
| 160 |
|
mc.SetProposalFunction(*pdfProp); |
| 161 |
|
mc.SetLeftSideTailFraction(0.5); // find a "central" interval |
| 162 |
|
MCMCInterval* mcInt = (MCMCInterval*)mc.GetInterval(); // that was easy |
| 170 |
|
config->add("RooSimpleProfile.signal."+Type+"UpperLimit", up); |
| 171 |
|
config->add("RooSimpleProfile.xsec."+Type+"UpperLimit", up/sig * xsec); |
| 172 |
|
|
| 173 |
< |
// Get Lower and Upper limits from FeldmanCousins with profile construction |
| 174 |
< |
if (fcint != NULL) { |
| 175 |
< |
double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s); |
| 176 |
< |
double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s); |
| 177 |
< |
cout << "FC lower limit on s = " << fcll << endl; |
| 178 |
< |
cout << "FC upper limit on s = " << fcul << endl; |
| 179 |
< |
config->add("RooFC.signal."+Type+"UpperLimit", fcul); |
| 180 |
< |
config->add("RooFC.xsec."+Type+"UpperLimit", fcul/sig * xsec); |
| 181 |
< |
//TLine* fcllLine = new TLine(fcll, 0, fcll, 1); |
| 182 |
< |
//TLine* fculLine = new TLine(fcul, 0, fcul, 1); |
| 183 |
< |
//fcllLine->SetLineColor(kRed); |
| 184 |
< |
//fculLine->SetLineColor(kRed); |
| 185 |
< |
//fcllLine->Draw("same"); |
| 186 |
< |
//fculLine->Draw("same"); |
| 187 |
< |
dataCanvas->Update(); |
| 188 |
< |
} |
| 173 |
> |
// // Get Lower and Upper limits from FeldmanCousins with profile construction |
| 174 |
> |
// if (fcint != NULL) { |
| 175 |
> |
// double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s); |
| 176 |
> |
// double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s); |
| 177 |
> |
// cout << "FC lower limit on s = " << fcll << endl; |
| 178 |
> |
// cout << "FC upper limit on s = " << fcul << endl; |
| 179 |
> |
// config->add("RooFC.signal."+Type+"UpperLimit", fcul); |
| 180 |
> |
// config->add("RooFC.xsec."+Type+"UpperLimit", fcul/sig * xsec); |
| 181 |
> |
// //TLine* fcllLine = new TLine(fcll, 0, fcll, 1); |
| 182 |
> |
// //TLine* fculLine = new TLine(fcul, 0, fcul, 1); |
| 183 |
> |
// //fcllLine->SetLineColor(kRed); |
| 184 |
> |
// //fculLine->SetLineColor(kRed); |
| 185 |
> |
// //fcllLine->Draw("same"); |
| 186 |
> |
// //fculLine->Draw("same"); |
| 187 |
> |
// dataCanvas->Update(); |
| 188 |
> |
// } |
| 189 |
|
/* |
| 190 |
|
|
| 191 |
|
// Plot MCMC interval and print some statistics |
| 251 |
|
double ExpNsigLimit = config.read<double>("ExpNsigLimit",100.0); |
| 252 |
|
|
| 253 |
|
simpleProfile2(&config, "Obs", data, bkg, bkgUncert, sig, sigUncert, xsec, ExpNsigLimit); |
| 254 |
< |
simpleProfile2(&config, "Exp", (int)bkg, bkg, bkgUncert, sig, sigUncert, xsec, ExpNsigLimit); |
| 254 |
> |
simpleProfile2(&config, "Exp", bkg, bkg, bkgUncert, sig, sigUncert, xsec, ExpNsigLimit); |
| 255 |
|
|
| 256 |
|
//write stuff: |
| 257 |
|
time_t rawtime; |
