You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
629 lines
30 KiB
629 lines
30 KiB
/**
|
|
* @file feldman_cousins.hh
|
|
* @author Christoph Langenbruch, David Gerick, Renata Kopecna
|
|
* @date 2021-01-12
|
|
*
|
|
*/
|
|
|
|
#include <feldman_cousins.hh>
|
|
|
|
//TODO: removed what all is unused
|
|
#include <iostream>
|
|
#include <complex>
|
|
#include <assert.h>
|
|
|
|
#include <event.hh>
|
|
#include <parameters.hh>
|
|
#include <funcs.hh>
|
|
#include <pdf.hh>
|
|
#include <plotter.hh>
|
|
#include <options.hh>
|
|
#include <bu2kstarmumu_generator.hh>
|
|
#include <fitter.hh>
|
|
#include "folder.hh"
|
|
|
|
#include <thread>
|
|
#include <mutex>
|
|
|
|
#include <TMinuit.h>
|
|
#include <TFile.h>
|
|
#include <TTree.h>
|
|
#include <TRandom3.h>
|
|
#include <TMatrix.h>
|
|
#include <TH1D.h>
|
|
#include <TCanvas.h>
|
|
#include <TROOT.h>
|
|
#include <TStyle.h>
|
|
#include <TVector.h>
|
|
#include <TMatrixD.h>
|
|
#include <TMatrixDSym.h>
|
|
#include <TMatrixTSym.h>
|
|
#include <TMatrixT.h>
|
|
#include <TDecompChol.h>
|
|
|
|
|
|
#include <spdlog.h>
|
|
|
|
bool fcnc::feldman_cousins::fc_1d(std::string parname, unsigned int nsteps, double parmin, double parmax, unsigned int ntoys, unsigned int pdfidx,
|
|
std::vector<pdf*> probs, std::vector<parameters*> params, generator* gen, fitter* f,
|
|
std::vector<std::vector<event> *> events, unsigned int bin, int from, int to){
|
|
std::ostringstream sout;
|
|
bool TwoBins = false;
|
|
if(opts->TheQ2binsmin.size() == 8) TwoBins = false;
|
|
else if(opts->TheQ2binsmin.size() == 2) TwoBins = true;
|
|
else assert(opts->TheQ2binsmin.size() == 2 || opts->TheQ2binsmin.size() == 8);
|
|
|
|
sout << "FCresult" << (TwoBins ? "_2BINS_" : "") << bin << "_" << parname << "_" << from << "_" << to << ".root";
|
|
struct stat buffer;
|
|
if(stat (sout.str().c_str(), &buffer) == 0){
|
|
spdlog::warn("[SKIP]\t\tFC job already run as root-file '" + sout.str() + "' exists, continue to next job!");
|
|
return false;
|
|
}
|
|
else{
|
|
spdlog::info("Feldman-Cousins study for parameter="+ parname + " in q2bin={0:d}", bin);
|
|
}
|
|
// if(!opts->full_angular)
|
|
// opts->always_generate_full_angular = false;
|
|
bool cache_minos = opts->minos_errors;
|
|
bool cache_hesse = opts->hesse_postrun;
|
|
bool cache_sq_hesse = opts->squared_hesse;
|
|
bool cache_shift = opts->shift_lh;
|
|
//int cache_print = spdlog::default_logger_raw()->level();
|
|
|
|
const unsigned int nPDF = probs.size();
|
|
assert(params.size() == nPDF);
|
|
assert(events.size() == nPDF);
|
|
assert(pdfidx < nPDF);
|
|
|
|
//temp save parameter set
|
|
unsigned int nparameters = params.at(pdfidx)->nparameters();
|
|
std::vector<std::vector<double> > smvalues (nPDF, std::vector<double> (nparameters, 0.0));
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for (unsigned int k=0; k<nparameters; k++){
|
|
smvalues.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
|
|
}
|
|
}
|
|
|
|
parameter* par = params.at(pdfidx)->get_parameter(parname);
|
|
if (par != 0){
|
|
opts->minos_errors = false;
|
|
opts->squared_hesse = false;
|
|
opts->hesse_postrun = false;
|
|
opts->simplex_prerun = false;
|
|
opts->shift_lh = false;
|
|
|
|
opts->minuit_strategy = 2.0;
|
|
//fitter f(opts);
|
|
int statusdatafloated = f->fit(probs, params, events);
|
|
double lhdatafloated = f->likelihood();
|
|
double valuedatafloated = par->get_value();
|
|
spdlog::info("Data fit finished with -2logLh = {0:f}",lhdatafloated );
|
|
std::vector<std::vector<double> > floatingvalues(nPDF, std::vector<double> (nparameters, 0.0));
|
|
std::vector<std::vector<double> > prev_error_low(nPDF, std::vector<double> (nparameters, 0.0));
|
|
std::vector<std::vector<double> > prev_error_high(nPDF, std::vector<double> (nparameters, 0.0));
|
|
std::vector<std::vector<double> > prev_meas(nPDF, std::vector<double> (nparameters, 0.0));
|
|
std::vector<std::vector<bool> > isConstraint(nPDF, std::vector<bool> (nparameters, false));
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for (unsigned int k=0; k<nparameters; k++){
|
|
floatingvalues.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
|
|
prev_error_low.at(n).at(k) = params.at(n)->get_parameter(k)->get_previous_error_low();
|
|
prev_error_high.at(n).at(k) = params.at(n)->get_parameter(k)->get_previous_error_high();
|
|
prev_meas.at(n).at(k) = params.at(n)->get_parameter(k)->get_previous_measurement();
|
|
isConstraint.at(n).at(k) = params.at(n)->get_parameter(k)->get_gaussian_constraint();
|
|
}
|
|
}
|
|
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for (unsigned int k=0; k<nparameters; k++){
|
|
if(isConstraint.at(n).at(k)){
|
|
spdlog::info("[PDF{0:d}][PAR{1:d}]\tpar={2:s}: prev_measurement={3:f} prev_low={4:f}prev_high={5:f}",
|
|
n, k, params.at(n)->get_parameter(k)->get_name(),
|
|
prev_meas.at(n).at(k),prev_error_low.at(n).at(k) ,prev_error_high.at(n).at(k));
|
|
}
|
|
}
|
|
}
|
|
|
|
//opts->minuit_strategy = 1.0;//this should be 2 for the fixed data fit I think
|
|
if (from < 0) from = 0;
|
|
if (to > int(nsteps) || to < 0) to = nsteps;
|
|
std::vector<double> lhdatafixed(to-from, 0.0);
|
|
std::vector<int> statusdatafixed(to-from, 0);
|
|
std::vector<std::vector<double> > lhtoysfixed;
|
|
std::vector<std::vector<double> > lhtoysfloated;
|
|
|
|
std::vector<std::vector<int> > statustoysfixed;
|
|
std::vector<std::vector<int> > statustoysfloated;
|
|
|
|
std::vector<std::vector< double> > params_dfloated (nPDF, std::vector< double> (nparameters, 0.0));
|
|
std::vector<std::vector<std::vector<double>>>params_dfixed(nPDF,std::vector<std::vector<double>>
|
|
(to-from, std::vector<double>(nparameters, 0.0)));
|
|
std::vector<std::vector<std::vector< std::vector< double> > > > params_tfixed (nPDF,
|
|
std::vector<std::vector< std::vector< double> > > (to - from,
|
|
std::vector< std::vector< double> > (ntoys,
|
|
std::vector< double> (nparameters, 0.0))));
|
|
std::vector<std::vector<std::vector< std::vector< double> > > > params_tfloated (nPDF,
|
|
std::vector<std::vector< std::vector< double> > > (to - from,
|
|
std::vector< std::vector< double> > (ntoys,
|
|
std::vector< double> (nparameters, 0.0))));
|
|
|
|
//buffer floated data values
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for (unsigned int k=0; k<nparameters; k++){
|
|
params_dfloated.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
|
|
}
|
|
}
|
|
|
|
std::string pardesc = par->get_description();
|
|
for (int j=from; j<to; j++){//loop over steps
|
|
double parcur = parmin+(j+0.5)*(parmax-parmin)/nsteps;
|
|
spdlog::info("[STEP]\tEvaluating step {0:d}out of {1:d} steps!",j-from+1, to-from);
|
|
spdlog::info("[INFO]\tStudying point "+parname+"={0:f}",parcur);
|
|
|
|
//init parameters from floating fit values
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for(unsigned int k=0; k<nparameters; k++){
|
|
if(isConstraint.at(n).at(k)){
|
|
params.at(n)->get_parameter(k)->init(floatingvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
|
|
params.at(n)->get_parameter(k)->get_step_size(), prev_error_low.at(n).at(k), prev_error_high.at(n).at(k));
|
|
params.at(n)->get_parameter(k)->set_previous_measurement(prev_meas.at(n).at(k));
|
|
}
|
|
else{
|
|
params.at(n)->get_parameter(k)->init(floatingvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
|
|
params.at(n)->get_parameter(k)->get_step_size());
|
|
}
|
|
}
|
|
params.at(n)->reset_parameters();
|
|
}
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
params.at(n)->get_parameter(parname)->init(parcur, parmin, parmax, 0.0);
|
|
}
|
|
opts->minuit_strategy = 2.0;//more precise for fixed data fit (used for lh)
|
|
opts->hesse_postrun = false;
|
|
opts->simplex_prerun = false;
|
|
spdlog::info("Start fit with fixed parameter");
|
|
statusdatafixed.at(j-from) = f->fit(probs, params, events);
|
|
lhdatafixed.at(j-from) = f->likelihood();
|
|
spdlog::info("[DONE]\tFinished fit with fixed parameter");
|
|
opts->minuit_strategy = 1.0;//less precise for toys -> speed more important
|
|
|
|
opts->hesse_postrun = false;
|
|
opts->simplex_prerun = false;
|
|
|
|
std::vector<std::vector<double> > fitvalues (nPDF, std::vector<double> (nparameters, 0.0));
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
//parameters from fixed point
|
|
for (unsigned int k=0; k<nparameters; k++){
|
|
fitvalues.at(n).at(k) = params.at(n)->get_parameter(k)->get_value();
|
|
}
|
|
//buffer fixed data values
|
|
for (unsigned int k=0; k<nparameters; k++){
|
|
params_dfixed.at(n).at(j-from).at(k) = params.at(n)->get_parameter(k)->get_value();
|
|
}
|
|
}
|
|
|
|
//make all toys with these params!
|
|
std::vector<std::vector<std::vector<event> > >toys(ntoys, std::vector<std::vector<event> > ());
|
|
spdlog::info("Start generating toys");
|
|
for (unsigned int k = 0; k < ntoys; k++){
|
|
spdlog::debug("[TOYS]\tGenerating toy sample {0:d}/{1:d}",k+1,ntoys );
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
std::vector<event> toy = gen->generate(events.at(n)->size(), params.at(n), probs.at(n));
|
|
if(!opts->full_angular){
|
|
fcnc::folder fldr(opts);
|
|
for(UInt_t e = 0; e < toy.size(); e++){
|
|
fldr.fold(&toy.at(e));
|
|
if(toy.at(e).m < 1000.){
|
|
spdlog::warn("[CORRUPT]\tEvent={0:d} pdf={0:f} toy={0:f}",e, n, k );
|
|
spdlog::warn("[CORRUPT]\tctl={0:f} ctk={0:f} phi={0:f} q2={0:f}",
|
|
toy.at(e).costhetal, toy.at(e).costhetak, toy.at(e).phi, toy.at(e).q2 );
|
|
}
|
|
}
|
|
}
|
|
toys.at(k).push_back(toy);
|
|
}
|
|
}
|
|
spdlog::info("[DONE]\tGenerating toys");
|
|
|
|
//init parameters from fixed point
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for(unsigned int k=0; k<nparameters; k++){
|
|
if(isConstraint.at(n).at(k)){
|
|
params.at(n)->get_parameter(k)->init(fitvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
|
|
params.at(n)->get_parameter(k)->get_step_size(), prev_error_low.at(n).at(k), prev_error_high.at(n).at(k));
|
|
params.at(n)->get_parameter(k)->set_previous_measurement(prev_meas.at(n).at(k));
|
|
}
|
|
else{
|
|
params.at(n)->get_parameter(k)->init(fitvalues.at(n).at(k), params.at(n)->get_parameter(k)->get_min(), params.at(n)->get_parameter(k)->get_max(),
|
|
params.at(n)->get_parameter(k)->get_step_size());
|
|
}
|
|
}
|
|
}
|
|
//fit toys repeatedly
|
|
std::vector<double> toysfixed;
|
|
std::vector<double> toysfloated;
|
|
std::vector<int> statusfixed;
|
|
std::vector<int> statusfloated;
|
|
spdlog::info("Start fitting all toys");
|
|
for (unsigned int k=0; k<ntoys; k++){
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
params.at(n)->reset_parameters();//these are not the fitted/generated parameters->fix, they should float however, so should be ok?
|
|
}
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
params.at(n)->get_parameter(parname)->init(parcur, parmin, parmax, 0.0);
|
|
}
|
|
|
|
//save pointers to toy events into vector:
|
|
std::vector<std::vector<event> * > leToys;
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
leToys.push_back(&toys.at(k).at(n));
|
|
}
|
|
|
|
//fit with fixed parameter
|
|
spdlog::info("[TOYS]\tFitting fixed toy sample {0:d}/{1:d} at step {2:d}/{3:d}",
|
|
k+1, ntoys, j-from+1, to-from );
|
|
int resultfixed = f->fit(probs, params, leToys);
|
|
toysfixed.push_back(f->likelihood());
|
|
statusfixed.push_back(resultfixed);
|
|
|
|
//buffer fixed toy values
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for (unsigned int l=0; l<nparameters; l++){
|
|
params_tfixed.at(n).at(j-from).at(k).at(l) = params.at(n)->get_parameter(l)->get_value();
|
|
}
|
|
params.at(n)->reset_parameters();
|
|
|
|
}
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
params.at(n)->get_parameter(parname)->init(parcur, parmin, parmax, 0.01);
|
|
}
|
|
|
|
//fit with floated parameter
|
|
spdlog::info("[TOYS]\tFitting freed toy sample0:d}/{1:d} at step {2:d}/{3:d}",
|
|
k+1, ntoys, j-from+1, to-from );
|
|
int resultfloated = f->fit(probs, params, leToys);
|
|
toysfloated.push_back(f->likelihood());
|
|
statusfloated.push_back(resultfloated);
|
|
|
|
//buffer fixed toy values
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
for (unsigned int l=0; l<nparameters; l++){
|
|
params_tfloated.at(n).at(j-from).at(k).at(l) = params.at(n)->get_parameter(l)->get_value();
|
|
}
|
|
}
|
|
}
|
|
lhtoysfixed.push_back(toysfixed);
|
|
lhtoysfloated.push_back(toysfloated);
|
|
statustoysfixed.push_back(statusfixed);
|
|
statustoysfloated.push_back(statusfloated);
|
|
}// end loop over steps
|
|
spdlog::info("[DONE]\tFitting all toys");
|
|
spdlog::debug("tLoop over all {0:d} steps completed. Writing results to file!",nsteps );
|
|
|
|
//save to root file
|
|
TFile* fout = new TFile(sout.str().c_str(), "RECREATE");
|
|
fout->cd();
|
|
TTree* t = new TTree(parname.c_str(), pardesc.c_str());
|
|
int pdf=0, step=0, statustfixed=0, statustfloated=0, toy=0, statusdfixed=0, statusdfloated=0;
|
|
double value=0.0, lhdfloated=lhdatafloated, lhdfixed=0.0, lhtfixed=0.0, lhtfloated=0.0;
|
|
t->Branch("bin",&bin,"bin/I");
|
|
t->Branch("pdf",&pdf,"pdf/I");
|
|
t->Branch("step",&step,"step/I");
|
|
t->Branch("value",&value,"value/D");
|
|
t->Branch("toy",&toy,"toy/I");
|
|
t->Branch("statusdatafixed",&statusdfixed,"statusdatafixed/I");
|
|
t->Branch("statusdatafloated",&statusdfloated,"statusdatafloated/I");
|
|
t->Branch("statustoyfixed",&statustfixed,"statustoyfixed/I");
|
|
t->Branch("statustoyfloated",&statustfloated,"statustoyfloated/I");
|
|
t->Branch("lhdatafloated",&lhdfloated,"lhdatafloated/D");
|
|
t->Branch("valuedatafloated",&valuedatafloated,"valuedatafloated/D");
|
|
t->Branch("lhdatafixed",&lhdfixed,"lhdatafixed/D");
|
|
t->Branch("lhtoyfixed",&lhtfixed,"lhtoyfixed/D");
|
|
t->Branch("lhtoyfloated",&lhtfloated,"lhtoyfloated/D");
|
|
//also save the fitted values
|
|
std::vector<double> tfloated(nparameters, 0.0), tfixed(nparameters, 0.0), dfloated(nparameters, 0.0), dfixed(nparameters, 0.0);
|
|
for (unsigned int k = 0; k < nparameters; k++){
|
|
fcnc::parameter* param = params.at(pdfidx)->get_parameter(k);
|
|
if (param->get_step_size() != 0.0){
|
|
std::string parname(param->get_name());
|
|
std::string pardesc(param->get_description());
|
|
//TTree* t = new TTree(parname.c_str(), pardesc.c_str());
|
|
t->Branch((std::string("datafixed_")+parname).c_str(),&dfixed.at(k),(std::string("datafixed_")+parname+std::string("/D")).c_str());
|
|
t->Branch((std::string("datafloated_")+parname).c_str(),&dfloated.at(k),(std::string("datafloated_")+parname+std::string("/D")).c_str());
|
|
t->Branch((std::string("toyfixed_")+parname).c_str(),&tfixed.at(k),(std::string("toyfixed_")+parname+std::string("/D")).c_str());
|
|
t->Branch((std::string("toyfloated_")+parname).c_str(),&tfloated.at(k),(std::string("toyfloated_")+parname+std::string("/D")).c_str());
|
|
}
|
|
}
|
|
for(unsigned int n = 0; n < nPDF; n++){
|
|
pdf = n; //TODO: what is this supposed to do?
|
|
for (int j=from; j<to; j++){
|
|
step = j;
|
|
for (unsigned int k=0; k<ntoys; k++){
|
|
toy = k;
|
|
value = parmin+(j+0.5)*(parmax-parmin)/nsteps;
|
|
statustfixed = statustoysfixed.at(j-from).at(k);
|
|
statustfloated = statustoysfloated.at(j-from).at(k);
|
|
statusdfloated = statusdatafloated;
|
|
statusdfixed = statusdatafixed.at(j-from);
|
|
lhtfixed = lhtoysfixed.at(j-from).at(k);
|
|
lhtfloated = lhtoysfloated.at(j-from).at(k);
|
|
lhdfixed = lhdatafixed.at(j-from);
|
|
//lhdfloated = lhdatafloated; //always the same
|
|
for (unsigned int l = 0; l < nparameters; l++){
|
|
tfloated.at(l) = params_tfloated.at(n).at(j-from).at(k).at(l);
|
|
tfixed.at(l) = params_tfixed.at(n).at(j-from).at(k).at(l);
|
|
dfloated.at(l) = params_dfloated.at(n).at(l);
|
|
dfixed.at(l) = params_dfixed.at(n).at(j-from).at(l);
|
|
}
|
|
t->Fill();
|
|
}
|
|
}
|
|
}
|
|
t->Write();
|
|
fout->Write();
|
|
fout->Close();
|
|
delete fout;
|
|
/*
|
|
//analyse results
|
|
std::string suffix = "";
|
|
TH1D* cl = new TH1D((std::string("cl")+parname+suffix).c_str(), (std::string(";")+pardesc+std::string(";Confidence level")).c_str(), nsteps, parmin, parmax);
|
|
TH1D* lh = new TH1D((std::string("lh")+parname+suffix).c_str(), (std::string(";")+pardesc+std::string(";Confidence level")).c_str(), nsteps, parmin, parmax);
|
|
for (unsigned int j=0; j<nsteps; j++)
|
|
{
|
|
unsigned int nlarger = 0;
|
|
for (unsigned int k=0; k<ntoys; k++)
|
|
{
|
|
double dtoy = lhtoysfixed.at(j).at(k)-lhtoysfloated.at(j).at(k);
|
|
double ddata = lhdatafixed.at(j)-lhdatafloated;
|
|
if (dtoy > ddata)
|
|
nlarger++;
|
|
}
|
|
lh->SetBinContent(j+1, TMath::Erf(sqrt(lhdatafixed.at(j)-lhdatafloated)/sqrt(2.0)));
|
|
cl->SetBinContent(j+1, 1.0-double(nlarger)/ntoys);
|
|
}
|
|
double low = parmin;
|
|
const double thecl = 0.683;
|
|
if (cl->GetBinContent(1) > thecl)//???
|
|
{
|
|
for (int i=0; i<nsteps; i++)
|
|
if (cl->GetBinContent(i+1) > thecl && cl->GetBinContent(i+2) < thecl)
|
|
{
|
|
low = cl->GetBinCenter(i+1) + (thecl-cl->GetBinContent(i+1))*(cl->GetBinCenter(i+2)-cl->GetBinCenter(i+1))/(cl->GetBinContent(i+2)-cl->GetBinContent(i+1));
|
|
break;
|
|
}
|
|
}
|
|
double high = parmax;
|
|
if (cl->GetBinContent(nsteps) > thecl)
|
|
{
|
|
for (int i=nsteps-1; i>=0; i--)
|
|
if (cl->GetBinContent(i+1) < thecl && cl->GetBinContent(i+2) > thecl)
|
|
{
|
|
high = cl->GetBinCenter(i+1) + (thecl-cl->GetBinContent(i+1))*(cl->GetBinCenter(i+2)-cl->GetBinCenter(i+1))/(cl->GetBinContent(i+2)-cl->GetBinContent(i+1));
|
|
break;
|
|
}
|
|
}
|
|
TCanvas* c1 = new TCanvas("c1", "c1", 1600, 1200);
|
|
TLine* l = new TLine();
|
|
c1->cd();
|
|
cl->SetMinimum(0.0);
|
|
cl->SetMaximum(1.1);
|
|
cl->Draw();
|
|
lh->SetLineColor(4);
|
|
lh->Draw("lsame");
|
|
l->SetLineColor(2);
|
|
l->DrawLine(low, 0.0, low, 1.1);
|
|
l->DrawLine(high, 0.0, high, 1.1);
|
|
c1->Print((std::string("plots/cl")+parname+suffix+std::string(".eps")).c_str(),"eps");
|
|
//delete cl;
|
|
delete l;
|
|
delete c1;
|
|
delete cl;
|
|
*/
|
|
}
|
|
else{
|
|
spdlog::critical("Could not find parameter "+parname);
|
|
assert(0);
|
|
}
|
|
spdlog::info("Finished Feldman-Cousins study for parameter "+ parname);
|
|
opts->hesse_postrun = cache_hesse;
|
|
opts->squared_hesse = cache_sq_hesse;
|
|
opts->minos_errors = cache_minos;
|
|
opts->shift_lh = cache_shift;
|
|
return true;
|
|
}
|
|
|
|
bool fcnc::feldman_cousins::fc_1d(std::string parname, unsigned int nsteps, double parmin, double parmax, unsigned int ntoys,
|
|
pdf* prob, parameters* params, generator* gen, fitter* f,
|
|
std::vector<event> * events, unsigned int bin, int from, int to){
|
|
std::vector<pdf *> the_probs;
|
|
the_probs.push_back(prob);
|
|
|
|
std::vector<parameters *> the_params;
|
|
the_params.push_back(params);
|
|
|
|
std::vector<std::vector<event>* > the_events;
|
|
the_events.push_back(events);
|
|
|
|
return fc_1d(parname, nsteps, parmin, parmax, ntoys, 0,
|
|
the_probs, the_params, gen, f,
|
|
the_events, bin, from, to);
|
|
}
|
|
|
|
bool fcnc::feldman_cousins::fc_2d(std::string parxname, unsigned int nstepsx, double parxmin, double parxmax,
|
|
std::string paryname, unsigned int nstepsy, double parymin, double parymax,
|
|
unsigned int ntoys, unsigned int pdfidx, //this parameter just changes output file
|
|
std::vector<pdf*> probs, std::vector<parameters*> params, generator* gen, fitter* f,
|
|
std::vector<std::vector<event> *> events, unsigned int bin, int from, int to){
|
|
spdlog::info("Feldman-Cousins study for parameters" + parxname +", " +paryname );
|
|
/*
|
|
bool cache_minos = opts->minos_errors;
|
|
bool cache_hesse = opts->hesse_postrun;
|
|
bool cache_shift = opts->shift_lh;
|
|
int cache_print = opts->print_level;
|
|
|
|
const unsigned int nPDF = probs.size();
|
|
assert(params.size() == nPDF);
|
|
assert(events.size() == nPDF);
|
|
assert(pdfidx < nPDF);
|
|
|
|
parameter* parx = params.at(pdfidx)->get_parameter(parxname);
|
|
parameter* pary = params.at(pdfidx)->get_parameter(paryname);
|
|
if (parx != 0 && pary != 0)
|
|
{
|
|
opts->minos_errors = false;
|
|
opts->hesse_postrun = false;
|
|
opts->shift_lh = false;
|
|
|
|
//fitter f(opts);
|
|
f->fit(probs, params, events);
|
|
double lhdatafloated = f.likelihood();
|
|
double valuexdatafloated = parx->get_value();
|
|
double valueydatafloated = pary->get_value();
|
|
spdlog::info("Data fit finished with -2logLh = {0:f}",lhdatafloated );
|
|
|
|
opts->print_level = -1;
|
|
//std::cout << from << " {0:f}",to << " {0:f}",from - to << std::endl;
|
|
if (from < 0)
|
|
from = 0;
|
|
if (to > int(nstepsx*nstepsy) || to < 0)
|
|
to = nstepsx*nstepsy;
|
|
std::vector<double> lhdatafixed(to-from, 0.0);
|
|
std::vector<std::vector<double> > lhtoysfixed;
|
|
std::vector<std::vector<double> > lhtoysfloated;
|
|
|
|
std::vector<std::vector<int> > statustoysfixed;
|
|
std::vector<std::vector<int> > statustoysfloated;
|
|
|
|
std::string parxdesc = parx->get_description();
|
|
std::string parydesc = pary->get_description();
|
|
for (int j=from; j<to; j++)
|
|
{
|
|
//xidx = j%nstepsx
|
|
//yidx = j/nstepsx
|
|
double parxcur = parxmin+((j%nstepsx)+0.5)*(parxmax-parxmin)/nstepsx;
|
|
double parycur = parymin+((j/nstepsx)+0.5)*(parymax-parymin)/nstepsy;
|
|
|
|
spdlog::info("Studying point" +parxname +"="+ parxcur + paryname + "=" +parycur );
|
|
params.reset_parameters();
|
|
|
|
parx->init(parxname, parxdesc, parxcur, parxmin, parxmax, 0.0);
|
|
pary->init(paryname, parydesc, parycur, parymin, parymax, 0.0);
|
|
f->fit(probs, params, events);
|
|
lhdatafixed.at(j-from) = f->likelihood();
|
|
|
|
//make all toys with these params!
|
|
std::vector<std::vector<event> > toys;
|
|
for (unsigned int k=0; k<ntoys; k++)
|
|
{
|
|
std::vector<event> toy = gen->generate(events.size(), params, prob);
|
|
toys.push_back(toy);
|
|
}
|
|
//fit toys repeatedly
|
|
std::vector<double> toysfixed;
|
|
std::vector<double> toysfloated;
|
|
std::vector<int> statusfixed;
|
|
std::vector<int> statusfloated;
|
|
for (unsigned int k=0; k<ntoys; k++)
|
|
{
|
|
params.reset_parameters();//these are not the fitted/generated parameters->fix, they should float however, so should be ok?
|
|
parx->init(parxname, parxdesc, parxcur, parxmin, parxmax, 0.0);
|
|
pary->init(paryname, parydesc, parycur, parymin, parymax, 0.0);
|
|
int resultfixed = f->fit(probs, params, &toys.at(k));
|
|
toysfixed.push_back(f->likelihood());
|
|
statusfixed.push_back(resultfixed);
|
|
|
|
params.reset_parameters();
|
|
parx->init(parxname, parxdesc, parxcur, parxmin, parxmax, 0.01);
|
|
pary->init(paryname, parydesc, parycur, parymin, parymax, 0.01);
|
|
int resultfloated = f->fit(probs, params, &toys.at(k));
|
|
toysfloated.push_back(f->likelihood());
|
|
statusfloated.push_back(resultfloated);
|
|
}
|
|
lhtoysfixed.push_back(toysfixed);
|
|
lhtoysfloated.push_back(toysfloated);
|
|
statustoysfixed.push_back(statusfixed);
|
|
statustoysfloated.push_back(statusfloated);
|
|
}
|
|
//std::cout <<"after loop" );
|
|
|
|
//save to root file
|
|
std::ostringstream sout;
|
|
sout << "FCresult_" << bin << "_" << parxname << "_" << paryname << "_" << from << "_" << to << ".root";
|
|
TFile* fout = new TFile(sout.str().c_str(), "RECREATE");
|
|
fout->cd();
|
|
TTree* t = new TTree((parxname+std::string("_")+paryname).c_str(), (std::string(";")+parxdesc+std::string(";")+parydesc).c_str());
|
|
int step=0, stepx=0, stepy=0, statustfixed=0, statustfloated=0, toy=0;
|
|
double valuex=0.0, valuey, lhdfloated=lhdatafloated, lhdfixed=0.0, lhtfixed=0.0, lhtfloated=0.0;
|
|
t->Branch("bin",&bin,"bin/I");
|
|
t->Branch("step",&step,"step/I");
|
|
t->Branch("stepx",&stepx,"stepx/I");
|
|
t->Branch("stepy",&stepy,"stepy/I");
|
|
t->Branch("valuex",&valuex,"valuex/D");
|
|
t->Branch("valuey",&valuey,"valuey/D");
|
|
t->Branch("valuexdatafloated",&valuexdatafloated,"valuexdatafloated/D");
|
|
t->Branch("valueydatafloated",&valueydatafloated,"valueydatafloated/D");
|
|
t->Branch("toy",&toy,"toy/I");
|
|
t->Branch("statustoyfixed",&statustfixed,"statustoyfixed/I");
|
|
t->Branch("statustoyfloated",&statustfloated,"statustoyfloated/I");
|
|
t->Branch("lhdatafloated",&lhdfloated,"lhdatafloated/D");
|
|
t->Branch("lhdatafixed",&lhdfixed,"lhdatafixed/D");
|
|
t->Branch("lhtoyfixed",&lhtfixed,"lhtoyfixed/D");
|
|
t->Branch("lhtoyfloated",&lhtfloated,"lhtoyfloated/D");
|
|
for (int j=from; j<to; j++)
|
|
{
|
|
for (unsigned int k=0; k<ntoys; k++)
|
|
{
|
|
step = j;
|
|
stepx = j%nstepsx;
|
|
stepy = j/nstepsx;
|
|
toy = k;
|
|
valuex = parxmin+((j%nstepsx)+0.5)*(parxmax-parxmin)/nstepsx;
|
|
valuey = parymin+((j/nstepsx)+0.5)*(parymax-parymin)/nstepsy;
|
|
statustfixed = statustoysfixed.at(j-from).at(k);
|
|
statustfloated = statustoysfloated.at(j-from).at(k);
|
|
lhtfixed = lhtoysfixed.at(j-from).at(k);
|
|
lhtfloated = lhtoysfloated.at(j-from).at(k);
|
|
lhdfixed = lhdatafixed.at(j-from);
|
|
lhdfloated = lhdatafloated;
|
|
t->Fill();
|
|
}
|
|
}
|
|
t->Write();
|
|
fout->Write();
|
|
fout->Close();
|
|
delete fout;
|
|
}
|
|
else{
|
|
if (parx == 0) spdlog::error("Could not find parameter" +parxname );
|
|
if (pary == 0) spdlog::error("Could not find parameter" +paryname );
|
|
assert(0);
|
|
}
|
|
spdlog::info("Finished Feldman-Cousins study for parameters "+parxname + "," +paryname );
|
|
opts->hesse_postrun = cache_hesse;
|
|
opts->minos_errors = cache_minos;
|
|
opts->shift_lh = cache_shift;
|
|
opts->print_level = cache_print;
|
|
*/
|
|
return true;
|
|
}
|
|
|
|
bool fcnc::feldman_cousins::fc_2d(std::string parxname, unsigned int nstepsx, double parxmin, double parxmax,
|
|
std::string paryname, unsigned int nstepsy, double parymin, double parymax,
|
|
unsigned int ntoys, //this parameter just changes output file
|
|
pdf* prob, parameters* params, generator* gen, fitter* f,
|
|
std::vector<event> * events, unsigned int bin, int from, int to){
|
|
std::vector<pdf *> the_probs;
|
|
the_probs.push_back(prob);
|
|
|
|
std::vector<parameters *> the_params;
|
|
the_params.push_back(params);
|
|
|
|
std::vector<std::vector<event>* > the_events;
|
|
the_events.push_back(events);
|
|
|
|
return fc_2d(parxname, nstepsx, parxmin, parxmax,
|
|
paryname, nstepsy, parymin, parymax,
|
|
ntoys, 0,
|
|
the_probs, the_params, gen, f,
|
|
the_events, bin, from, to);
|
|
}
|
|
|