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.
1260 lines
55 KiB
1260 lines
55 KiB
//Renata Kopecna
|
|
|
|
#include <fitter.hh>
|
|
#include <paths.hh>
|
|
|
|
#include <event.hh>
|
|
#include <parameters.hh>
|
|
#include <pdf.hh>
|
|
#include <plotter.hh>
|
|
#include <options.hh>
|
|
#include <helpers.hh>
|
|
#include <funcs.hh>
|
|
#include <design.hh>
|
|
#include <constants.hh>
|
|
|
|
#include <TH1D.h>
|
|
#include <TMatrixD.h>
|
|
#include <TMatrixDSym.h>
|
|
#include <TRandom3.h>
|
|
#include <TVector.h>
|
|
#include <TDecompChol.h>
|
|
|
|
#include <spdlog.h>
|
|
|
|
#include <fstream>
|
|
#include <iostream>
|
|
#include <complex>
|
|
#include <assert.h>
|
|
#include <thread>
|
|
#include <mutex>
|
|
|
|
//Log of the covariance matrix status
|
|
void printCovMatrixStatus(int status, bool squared){
|
|
std::string hesse = "Hesse";
|
|
if (squared) hesse = "Squared Hesse";
|
|
if (status==0) spdlog::warn(hesse + " returns 0: Hesse not calculated.");
|
|
else if (status==1) spdlog::warn(hesse + " returns 1: Approximation only.");
|
|
else if (status==2) spdlog::warn(hesse + " returns 2: Full, forced pos def");
|
|
else if (status==3) spdlog::info(hesse + " returns 3: ALL GOOD");
|
|
else spdlog::warn(hesse + " returns {0:d} -> SOMETHING VERY WRONG", status);
|
|
return;
|
|
}
|
|
|
|
//Migrad status log so one doesn't have to look it up all the time
|
|
void printMigradStatus(int status){
|
|
if (status==0) spdlog::info("MIGRAD returns 0: ALL GOOD");
|
|
else if (status==1) spdlog::warn("MIGRAD returns 1: Blank command, ignored.");
|
|
else if (status==2) spdlog::warn("MIGRAD returns 2: Command line unreadable, ignored.");
|
|
else if (status==3) spdlog::warn("MIGRAD returns 3: Unknown command, ignored");
|
|
else if (status==4) spdlog::warn("MIGRAD returns 4: Abnormal termination, MIGRAD not converged or something.");
|
|
else if (status==9) spdlog::warn("MIGRAD returns 9: Reserved.");
|
|
else if (status==10) spdlog::warn("MIGRAD returns 10: End command.");
|
|
else if (status==11) spdlog::warn("MIGRAD returns 11: Exit/Stop command.");
|
|
else if (status==12) spdlog::warn("MIGRAD returns 12: Return command.");
|
|
else spdlog::warn("MIGRAD returns {0:d} -> SOMETHING VERY WRONG", status);
|
|
return;
|
|
}
|
|
|
|
|
|
//set plotter
|
|
void fcnc::fitter::set_plotters(plotter* plot) {
|
|
std::vector<plotter* > p;
|
|
p.push_back(plot);
|
|
plots = p;
|
|
return;
|
|
}
|
|
|
|
///set plotters
|
|
void fcnc::fitter::set_plotters(std::vector<plotter*> plotters) {
|
|
plots = plotters;
|
|
return;
|
|
}
|
|
///define common parameters for all pdfs
|
|
void fcnc::fitter::set_common_parameters(std::vector<std::string> common_pars) {
|
|
common_params = common_pars;
|
|
return;
|
|
}
|
|
///is the parameter with name one of the common parameters
|
|
bool fcnc::fitter::is_common(std::string name) {
|
|
for (unsigned int i=0; i<common_params.size(); i++){
|
|
if (common_params.at(i) == name) return true;
|
|
}
|
|
return false;
|
|
}
|
|
///is this the first occurence of the parameter?
|
|
bool fcnc::fitter::is_first(std::string name, unsigned int idx_i, unsigned int idx_j) {
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (name == params.at(i)->get_parameter(j)->get_name()){
|
|
if (i==idx_i && j==idx_j) return true;
|
|
else return false;
|
|
}
|
|
}
|
|
}
|
|
spdlog::critical("[FATAL]\t\tThe parameter '"+name+"' was not found!");
|
|
assert(0);
|
|
return true;
|
|
}
|
|
|
|
///returns the nth minuit index for parameter with name
|
|
int fcnc::fitter::nth_minuit_index(std::string name, unsigned int n) {
|
|
assert(n);
|
|
assert(n <= params.size());
|
|
unsigned int start_param = 0;
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (name == params.at(i)->get_parameter(j)->get_name()){
|
|
if (n == 1) return start_param+j;
|
|
else n--;
|
|
}
|
|
}
|
|
start_param += minuit_one_fitter->params.at(i)->nparameters();
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
///returns the first minuit index for parameter with name
|
|
int fcnc::fitter::first_minuit_index(std::string name) {
|
|
return nth_minuit_index(name, 1);
|
|
}
|
|
|
|
///returns wether this is an analysis with blinded parameters
|
|
bool fcnc::fitter::is_blind() {
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->is_blind()){
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
//WHAT THE FUCK IS THIS DOING HERE?
|
|
int fcnc::fitter::fit(pdf* probs, parameters* parms, std::vector<event>* ev, std::string index){
|
|
std::vector<pdf *> the_probs;
|
|
the_probs.push_back(probs);
|
|
|
|
std::vector<parameters *> the_params;
|
|
the_params.push_back(parms);
|
|
|
|
std::vector<std::vector<event>* > the_events;
|
|
the_events.push_back(ev);
|
|
|
|
return fit(the_probs, the_params, the_events, index);
|
|
}
|
|
|
|
void fcnc::fitter::init(pdf* probs, parameters* parms, std::vector<event>* ev, std::string index){
|
|
std::vector<pdf *> the_probs;
|
|
the_probs.push_back(probs);
|
|
|
|
std::vector<parameters *> the_params;
|
|
the_params.push_back(parms);
|
|
|
|
std::vector<std::vector<event>* > the_events;
|
|
the_events.push_back(ev);
|
|
|
|
return init(the_probs, the_params, the_events, index);
|
|
}
|
|
|
|
void fcnc::fitter::minuit_one_fcn(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *parameters, Int_t iflag){
|
|
//set the parameter sets according to minuits current parameter values stored in the params pointer
|
|
//do not forget to set common parameters correctly. these have no directly corresponding minuit index, use first_index instead
|
|
unsigned int start_param = 0;
|
|
for (unsigned int i = 0; i < minuit_one_fitter->params.size(); i++){
|
|
for (unsigned int j = 0; j < minuit_one_fitter->params.at(i)->nparameters(); j++){
|
|
unsigned int index = start_param + j;
|
|
std::string par_name = minuit_one_fitter->params.at(i)->get_parameter(j)->get_name();
|
|
if (minuit_one_fitter->is_common(par_name))
|
|
minuit_one_fitter->params.at(i)->get_parameter(j)->set_value(parameters[minuit_one_fitter->first_minuit_index(par_name)]);
|
|
else
|
|
minuit_one_fitter->params.at(i)->get_parameter(j)->set_value(parameters[index]);
|
|
}
|
|
start_param += minuit_one_fitter->params.at(i)->nparameters();
|
|
}
|
|
lh = minuit_one_fitter->likelihood();
|
|
return;
|
|
}
|
|
|
|
fcnc::fitter::fitter(options* o):
|
|
opts(o)
|
|
{
|
|
minuit_one = new TMinuit(1000);//params->nparameters());
|
|
minuit_one->SetFCN(&(minuit_one_fcn));
|
|
//this makes minuit use highest precision
|
|
int errorcode;
|
|
double strategy_list[1];// = {2.0}
|
|
strategy_list[0] = opts->minuit_strategy;
|
|
minuit_one->mnexcm("SET STR", strategy_list, 1, errorcode);
|
|
minuit_one->SetMaxIterations(40000);
|
|
//need to set this value to 0.5 if you do not use -2*ln(LH)
|
|
minuit_one->SetErrorDef(1.0);//default
|
|
spdlog::info("Machine precision: {0:f}", minuit_one->fEpsmac); //8.88178e-16
|
|
//verbose minuit
|
|
//minuit_one->SetPrintLevel(o->print_level);
|
|
//minuit_one->SetPrintLevel(-1);
|
|
empirical_constant = 0.0;
|
|
}
|
|
|
|
fcnc::fitter::~fitter(){
|
|
delete minuit_one;
|
|
}
|
|
|
|
void fcnc::fitter::init(std::vector<pdf*> probs, std::vector<parameters*> parms, std::vector< std::vector<event> * > ev, std::string index){
|
|
square_weights = false;
|
|
//set all important pointers and initialize
|
|
minuit_one_fitter = this;
|
|
//this could be different for every fit (feldman-cousins)
|
|
int errorcode;
|
|
double strategy_list[1];// = {2.0}
|
|
strategy_list[0] = opts->minuit_strategy;
|
|
minuit_one->mnexcm("SET STR", strategy_list, 1, errorcode);
|
|
|
|
//prob = probability;
|
|
pdfs = probs;
|
|
params = parms;
|
|
events = ev;
|
|
|
|
reset_param_values();
|
|
//this was originally before the reset... put here to catch reset of swave
|
|
for (unsigned int i=0; i<pdfs.size(); i++)
|
|
pdfs.at(i)->init(params.at(i));
|
|
|
|
//update cached efficiencies in three situations
|
|
//1. q2 fixed and regular acceptance
|
|
//2. unfolded fit to determine the weight=1/eff
|
|
//3. per event norm
|
|
//cache fis
|
|
if (opts->cache_fis){
|
|
spdlog::info("Caching fis");
|
|
for (unsigned int i = 0; i < pdfs.size(); i++)
|
|
pdfs.at(i)->update_cached_integrated_fis(params.at(i));
|
|
}
|
|
for (unsigned int i = 0; i < pdfs.size(); i++){
|
|
parameter* effq2 = params.at(i)->get_parameter("eff_q2");
|
|
bool q2fixed = false;
|
|
if (effq2) q2fixed = (effq2->get_step_size() == 0.0);
|
|
if ((q2fixed && opts->use_angular_acc) || opts->weighted_fit || opts->use_event_norm){
|
|
pdfs.at(i)->update_cached_efficiencies(params.at(i), events.at(i));
|
|
}
|
|
}
|
|
//per event xis
|
|
if (opts->use_event_norm){
|
|
for (unsigned int i = 0; i < pdfs.size(); i++){
|
|
pdfs.at(i)->update_cached_xis(params.at(i), events.at(i));
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
int fcnc::fitter::fit(std::vector<pdf*> probs, std::vector<parameters*> parms, std::vector< std::vector<event> * > ev, std::string index){
|
|
//can only use either squared hesse or minos
|
|
if(opts->minos_errors && opts->squared_hesse){
|
|
spdlog::critical("Cannot use both MINOS and squared HESSE for uncertainty deterimination. Choose one!");
|
|
assert(0);
|
|
}
|
|
square_weights = false;
|
|
//set all important pointers and initialize
|
|
minuit_one_fitter = this;
|
|
//this could be different for every fit (feldman-cousins)
|
|
int errorcode;
|
|
double strategy_list[1];// = {2.0}
|
|
strategy_list[0] = opts->minuit_strategy;
|
|
minuit_one->mnexcm("SET STR", strategy_list, 1, errorcode);
|
|
|
|
//prob = probability;
|
|
pdfs = probs;
|
|
params = parms;
|
|
events = ev;
|
|
|
|
//make sure to reduce output for blind analysis
|
|
if (is_blind()) minuit_one->SetPrintLevel(-1);
|
|
else if (spdlog_debug()) minuit_one->SetPrintLevel(0);
|
|
else minuit_one->SetPrintLevel(-1);
|
|
|
|
reset_param_values();
|
|
//this was originally before the reset... put here to catch reset of swave
|
|
for (unsigned int i=0; i<pdfs.size(); i++){
|
|
pdfs.at(i)->init(params.at(i));
|
|
}
|
|
|
|
//update cached efficiencies in three situations
|
|
//1. q2 fixed and regular acceptance
|
|
//2. unfolded fit to determine the weight=1/eff
|
|
//3. per event norm
|
|
//cache fis
|
|
if (opts->cache_fis){
|
|
spdlog::info("Caching fis");
|
|
for (unsigned int i = 0; i < pdfs.size(); i++){
|
|
pdfs.at(i)->update_cached_integrated_fis(params.at(i));
|
|
}
|
|
}
|
|
for (unsigned int i = 0; i < pdfs.size(); i++){
|
|
parameter* effq2 = params.at(i)->get_parameter("eff_q2");
|
|
bool q2fixed = false;
|
|
if (effq2) q2fixed = (effq2->get_step_size() == 0.0);
|
|
if ((q2fixed && opts->use_angular_acc) || opts->weighted_fit || opts->use_event_norm){
|
|
pdfs.at(i)->update_cached_efficiencies(params.at(i), events.at(i));
|
|
}
|
|
}
|
|
//per event xis
|
|
if (opts->use_event_norm){
|
|
for (unsigned int i = 0; i < pdfs.size(); i++){
|
|
pdfs.at(i)->update_cached_xis(params.at(i), events.at(i));
|
|
}
|
|
}
|
|
|
|
//get empirical factor
|
|
empirical_constant = 0.0;
|
|
if (opts->shift_lh){
|
|
unsigned int the_size=0;
|
|
for (unsigned int i=0; i<ev.size(); i++) the_size += ev.at(i)->size();
|
|
empirical_constant = likelihood()/the_size;
|
|
}
|
|
|
|
spdlog::debug("Determined empirical constant of: {0:f}", empirical_constant);
|
|
|
|
int result = 0;
|
|
unsigned int varied_params = 0;
|
|
//determine the number of varied parameters
|
|
//count only parameters that are not common among the PDFs!
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0 && !is_common(params.at(i)->get_parameter(j)->get_name())){
|
|
varied_params++;
|
|
spdlog::debug("[PDF{0:d}]\tPAR='"+params.at(i)->get_parameter(j)->get_name()+"'\trange=[{1:f} - {2:f}]", i+1, params.at(i)->get_parameter(j)->get_min(), params.at(i)->get_parameter(j)->get_max());
|
|
}
|
|
}
|
|
}
|
|
|
|
spdlog::debug("Count number of varied parameters: {0:d} PDF(s)", params.size());
|
|
spdlog::debug("Number of varied, non-common parameters: {0:d}", varied_params);
|
|
//assert(!(varied_params % params.size()));
|
|
|
|
//add the parameters to the counter, that are common AND varied
|
|
for (unsigned int i = 0; i < common_params.size(); i++){
|
|
bool is_varied = false;
|
|
for (unsigned int j = 0; j < params.size(); j++){
|
|
for (unsigned int k = 0; k < params.at(j)->nparameters(); k++){
|
|
if (common_params.at(i) == params.at(j)->get_parameter(k)->get_name() && params.at(j)->get_parameter(k)->get_step_size() != 0.0){
|
|
is_varied = true;
|
|
spdlog::debug("[COMMON]\tPAR='"+params.at(j)->get_parameter(k)->get_name()+"'\trange=[{0:f} - {1:f}]", params.at(j)->get_parameter(k)->get_min(), params.at(j)->get_parameter(k)->get_max());
|
|
break;
|
|
}
|
|
}
|
|
if(is_varied)break;
|
|
}
|
|
if (is_varied){
|
|
varied_params++;
|
|
}
|
|
}
|
|
|
|
spdlog::info("Including common parameters, number of varied parameters is: {0:d}", varied_params);
|
|
|
|
//open texFile
|
|
std::ofstream myFile;
|
|
if (index != ""){
|
|
open_Latex_noteFile(latex_fitterFile(index), myFile);
|
|
}
|
|
|
|
|
|
//Use Markov Chain Monte Carlo //TODO when bored, move this to a separate function
|
|
if (opts->use_mcmc){
|
|
unsigned int nlength = 10000;
|
|
//just run a single chain
|
|
std::vector<TVectorD> chain;//(nlength);
|
|
for (unsigned int i=0; i<nlength; i++){
|
|
chain.push_back(TVectorD(varied_params));
|
|
}
|
|
TVectorD current(varied_params);
|
|
//fill in starting values
|
|
//TMatrixDSym S1(nparams);
|
|
//S1.IdentityMatrix();
|
|
//S1 *= 0.1;
|
|
TMatrixDSym SNminusone(varied_params);
|
|
SNminusone.UnitMatrix();
|
|
//SNminusone *= 0.1; //might be more appropriate for the problem
|
|
|
|
unsigned int idx = 0;
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0 && !is_common(params.at(i)->get_parameter(j)->get_name())){
|
|
current[idx] = params.at(i)->get_parameter(j)->get_value();
|
|
SNminusone(idx, idx) = params.at(i)->get_parameter(j)->get_step_size();//this is actually quite clever
|
|
idx++;
|
|
}
|
|
}
|
|
}
|
|
spdlog::debug("nrows {0:d}\tvaried_params {1:d}", current.GetNrows(), varied_params);
|
|
for(int i=0; i<current.GetNrows(); i++) spdlog::debug("{0:d}", current[i]);
|
|
|
|
TVectorD last(current);
|
|
|
|
double alphastar = 0.234;
|
|
double last_llh = likelihood();
|
|
double minimum_llh = last_llh;
|
|
TVectorD minimum(current);
|
|
|
|
TMatrixDSym SN(SNminusone);
|
|
TRandom3 * rnd = new TRandom3();
|
|
spdlog::info("Starting Markov chain");
|
|
//produce markov chain of the desired length
|
|
for (unsigned int i=0; i<nlength; i++){
|
|
if ((i*100)%nlength == 0) std::cout << i*100/double(nlength) << "%" << std::endl;
|
|
//should optionally repeat this, until values are in parameter range
|
|
TVectorD WN(varied_params);
|
|
bool finished = false;
|
|
while (!finished){//repeat until we are inside the boundaries
|
|
//get random vector
|
|
for (int j=0; j<WN.GetNrows(); j++){
|
|
WN[j] = rnd->Gaus(0.0, 1.0);//could also use students t distribution
|
|
}
|
|
//shift depends on SNminusone
|
|
TVectorD SW(SNminusone*WN);
|
|
current = last + SW;
|
|
finished = true;
|
|
//this will disallow parameters outside the allowed range
|
|
const bool constrain_parameters = true;
|
|
if (constrain_parameters){
|
|
idx = 0;
|
|
for (unsigned int j = 0; j < params.size(); j++){
|
|
for (unsigned int k = 0; k < params.at(j)->nparameters(); k++){
|
|
if (params.at(j)->get_parameter(k)->get_step_size() != 0.0 && !is_common(params.at(j)->get_parameter(k)->get_name())){
|
|
if (current[idx] > params.at(j)->get_parameter(k)->get_max() || current[idx] < params.at(j)->get_parameter(k)->get_min()){
|
|
finished = false;
|
|
}
|
|
idx++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
spdlog::debug("event {0:d}", i);
|
|
//convert index back to parameter* and set to current point
|
|
idx = 0;
|
|
for (unsigned int j = 0; j < params.size(); j++){
|
|
for (unsigned int k = 0; k < params.at(j)->nparameters(); k++){
|
|
if (params.at(j)->get_parameter(k)->get_step_size() != 0.0 && !is_common(params.at(j)->get_parameter(k)->get_name())){
|
|
params.at(j)->get_parameter(k)->set_value(current[idx]);
|
|
idx++;
|
|
}
|
|
}
|
|
}
|
|
for (unsigned int j=0; j<pdfs.size(); j++){//check if this is needed!
|
|
pdfs.at(j)->init(params.at(j));
|
|
}
|
|
//compute likelihood of current point
|
|
double current_llh = likelihood();
|
|
double r = rnd->Rndm();
|
|
double alpha = std::min(1.0, exp(0.5*(last_llh-current_llh)));//llhs are in fact Sum_i -2*log(P_i)
|
|
//debug output
|
|
if (spdlog_debug()){
|
|
for (int i=0; i<current.GetNrows(); i++) std::cout << current[i] << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
//check if current point is accepted
|
|
spdlog::debug("accept event?");
|
|
if (r < alpha){
|
|
chain.at(i) = current;
|
|
last = current;
|
|
last_llh = current_llh;
|
|
}
|
|
else{
|
|
chain.at(i) = last;
|
|
//last stays last
|
|
}
|
|
//debugging output
|
|
spdlog::debug("current_llh {0:f}\t last_llh {1:f}" + std::string((r<alpha) ? " accepted " : " rejected "), current_llh, last_llh);
|
|
//check wether a new likelihood minimum is found and set if this is the case
|
|
if (current_llh < minimum_llh){
|
|
minimum = current;
|
|
minimum_llh = current_llh;
|
|
}
|
|
//update SN
|
|
TMatrixDSym SNminusoneT(SNminusone);
|
|
SNminusoneT.T();
|
|
double etan = std::min(1.0, varied_params*pow(double(i), -2.0/3.0));
|
|
TMatrixDSym WNWNT(varied_params);
|
|
WNWNT.Zero();
|
|
for (int row = 0; row < WNWNT.GetNrows(); row++){
|
|
for (int col = 0; col < WNWNT.GetNcols(); col++){
|
|
WNWNT[row][col] = WN[row]*WN[col]/WN.Norm2Sqr();
|
|
}
|
|
}
|
|
//TMatrixDSym SNSNT(identity + WNWNT*etan*(alpha-alphastar));
|
|
TMatrixDSym SNSNT(varied_params);
|
|
SNSNT.UnitMatrix();
|
|
SNSNT += WNWNT*etan*(alpha-alphastar);
|
|
SNSNT = SNSNT.Similarity(SNminusone);
|
|
TDecompChol chol(SNSNT);
|
|
bool success = chol.Decompose();
|
|
assert(success);
|
|
TMatrixD SNT = chol.GetU();
|
|
TMatrixD SN(SNT);
|
|
SN.T();
|
|
for (int row = 0; row < SN.GetNrows(); row++){
|
|
for (int col = 0; col < SN.GetNcols(); col++){
|
|
SNminusone(row,col) = SN(row,col);
|
|
}
|
|
}
|
|
}
|
|
//best spot found
|
|
spdlog::info("lowest point found at llh of {0:f}\n", minimum_llh);
|
|
if (spdlog_info()){
|
|
for (int i=0; i<minimum.GetNrows(); i++) std::cout << minimum[i] << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
//mean
|
|
spdlog::info("mean point at");
|
|
TVectorD mean(varied_params);
|
|
mean.Zero();
|
|
for (unsigned int i=0; i<nlength; i++){
|
|
TVectorD v(chain.at(i));
|
|
v *= 1.0/double(nlength);
|
|
mean += v;
|
|
}
|
|
if (spdlog_info()){
|
|
for (int i=0; i<mean.GetNrows(); i++) std::cout << mean[i] << std::endl;
|
|
std::cout << std::endl;
|
|
}
|
|
//rms
|
|
spdlog::info("rms is ");
|
|
TVectorD rms(varied_params);
|
|
rms.Zero();
|
|
for (unsigned int i=0; i<nlength; i++){
|
|
TVectorD v(chain.at(i)-mean);
|
|
v.Sqr();
|
|
v *= 1.0/double(nlength);
|
|
rms += v;
|
|
}
|
|
rms.Sqrt();
|
|
if (spdlog_info()){
|
|
for (int i=0; i<rms.GetNrows(); i++) std::cout << rms[i] << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
//get covariance matrix estimate from matrix SN
|
|
spdlog::debug("matrix SNminousone");
|
|
if (spdlog_debug()){
|
|
for (int row = 0; row < SNminusone.GetNrows(); row++){
|
|
for (int col = 0; col < SNminusone.GetNcols(); col++) std::cout << std::setw(8) << SNminusone(row,col) << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
}
|
|
TMatrixD SNminusoneT(SNminusone);
|
|
SNminusoneT.T();
|
|
TMatrixD cov = SNminusone*SNminusoneT;
|
|
spdlog::debug("matrix cov");
|
|
if (spdlog_debug()){
|
|
for (int row = 0; row < cov.GetNrows(); row++){
|
|
for (int col = 0; col < cov.GetNcols(); col++) std::cout << std::setw(8) << cov(row,col) << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
}
|
|
spdlog::debug("errors from matrix cov");
|
|
|
|
if (spdlog_debug()){
|
|
for (int row = 0; row < cov.GetNrows(); row++)std::cout << sqrt(cov(row,row)) << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
//proper estimate of best point and uncertainties from projections
|
|
unsigned int nbins = 1000;
|
|
TH1D* histos[varied_params];
|
|
idx = 0;
|
|
for (unsigned int j = 0; j < params.size(); j++){
|
|
for (unsigned int k = 0; k < params.at(j)->nparameters(); k++){
|
|
if (params.at(j)->get_parameter(k)->get_step_size() != 0.0 && !is_common(params.at(j)->get_parameter(k)->get_name())){
|
|
std::ostringstream hname;
|
|
hname << "hist" << idx;
|
|
histos[idx] = new TH1D(hname.str().c_str(), ";;", nbins,
|
|
params.at(j)->get_parameter(k)->get_min(),
|
|
params.at(j)->get_parameter(k)->get_max());
|
|
idx++;
|
|
}
|
|
}
|
|
}
|
|
//loop over chain and fill histos
|
|
for (unsigned int i=0; i<nlength; i++){
|
|
TVectorD v(chain.at(i));
|
|
for (int row = 0; row < v.GetNrows(); row++) histos[row]->Fill(v[row]);
|
|
}
|
|
if (true){//smooth all histos
|
|
const unsigned int ntimes = 1;
|
|
for (unsigned int i=0; i<varied_params; i++) histos[i]->Smooth(ntimes);
|
|
}
|
|
if (spdlog_debug()){ //If in debug mode, just save the histos for now
|
|
for (unsigned int i=0; i<varied_params; i++){
|
|
plotAndSave(histos[i],"c"+std::to_string(i),std::string(histos[i]->GetName()),"eps");
|
|
}
|
|
}
|
|
//get maximum from histos (for now no smoothing) and add up highes bins until 1 sigma
|
|
//probably want to add from highest bin, no?
|
|
TVectorD max(varied_params);
|
|
TVectorD eup(varied_params);
|
|
TVectorD edown(varied_params);
|
|
for (unsigned int i=0; i<varied_params; i++){
|
|
int maxbin = histos[i]->GetMaximumBin();
|
|
double x = histos[i]->GetBinCenter(maxbin);
|
|
double left = histos[i]->GetBinLowEdge(maxbin);
|
|
double right = histos[i]->GetBinLowEdge(maxbin+1);
|
|
double sum = histos[i]->GetBinContent(maxbin);
|
|
histos[i]->SetBinContent(maxbin, 0.0);
|
|
bool finished = sum/nlength > ONE_SIGMA;
|
|
while (!finished){
|
|
maxbin = histos[i]->GetMaximumBin();
|
|
double nx = histos[i]->GetBinCenter(maxbin);
|
|
sum += histos[i]->GetBinContent(maxbin);
|
|
histos[i]->SetBinContent(maxbin, 0.0);
|
|
if (nx > right) right = histos[i]->GetBinLowEdge(maxbin+1);
|
|
if (nx < left) left = histos[i]->GetBinLowEdge(maxbin);
|
|
finished = sum/nlength > ONE_SIGMA;
|
|
}
|
|
max[i] = x;
|
|
eup[i] = right-x;
|
|
edown[i] = left-x;
|
|
delete histos[i];
|
|
}
|
|
//output of the newly determined mean and errors
|
|
|
|
if (spdlog_info()){
|
|
spdlog::info("max");
|
|
for (int i=0; i<max.GetNrows(); i++) std::cout << max[i] << " ";
|
|
std::cout << std::endl;
|
|
spdlog::info("eup");
|
|
for (int i=0; i<eup.GetNrows(); i++) std::cout << eup[i] << " ";
|
|
std::cout << std::endl;
|
|
spdlog::info("edown");
|
|
for (int i=0; i<edown.GetNrows(); i++) std::cout << edown[i] << " ";
|
|
std::cout << std::endl;
|
|
}
|
|
idx = 0;
|
|
for (unsigned int j = 0; j < params.size(); j++){
|
|
for (unsigned int k = 0; k < params.at(j)->nparameters(); k++){
|
|
if (params.at(j)->get_parameter(k)->get_step_size() != 0.0 && !is_common(params.at(j)->get_parameter(k)->get_name())){
|
|
params.at(j)->get_parameter(k)->set_value(max[idx]);
|
|
params.at(j)->get_parameter(k)->set_error(0.5*(fabs(eup[idx])+fabs(edown[idx])));
|
|
params.at(j)->get_parameter(k)->set_error_up(eup[idx]);
|
|
params.at(j)->get_parameter(k)->set_error_down(edown[idx]);
|
|
//params.at(j)->get_parameter(k)->set_value(mean[idx]);
|
|
//params.at(j)->get_parameter(k)->set_error(sqrt(cov(idx,idx)));
|
|
idx++;
|
|
}
|
|
}
|
|
}
|
|
//MCMC always returns status OK
|
|
return 300;//this means all ok
|
|
}
|
|
//When using no Markov Chain MC: Nominal fit procedure
|
|
else{ //TODO when bored, move this to a separate function
|
|
spdlog::debug("Starting Fit Procedure");
|
|
int errorcode;
|
|
//double eps[1] = {1.0e-10};
|
|
//minuit_one->mnexcm("SET EPS", eps, 1, errorcode);
|
|
//double strategy_list[2] = {0,0.1};
|
|
|
|
double migrad_options[2] = {MIG_MAX_CALLS, MIG_TOLERANCE}; //Maxcalls and tolerance
|
|
// The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP(see SET ERR)
|
|
if (opts->simplex_prerun){
|
|
spdlog::info("Running Simplex");
|
|
minuit_one->mnexcm("SIM", migrad_options, 2, errorcode);
|
|
}
|
|
minuit_one->mnexcm("MIG", migrad_options, 2, errorcode);
|
|
result = errorcode;
|
|
printMigradStatus(result);
|
|
spdlog::warn("EDM: {0:f}", minuit_one->fEDM);
|
|
if (opts->hesse_postrun){
|
|
spdlog::info("Starting hesse");
|
|
minuit_one->mnhess();
|
|
if (spdlog_debug()) minuit_one->mnmatu(1); //Print the correlation matrix, they are annoying so turned off for now
|
|
}
|
|
int nfree, ntot, status_cov;
|
|
double m_fcn, m_edm, up;
|
|
minuit_one->mnstat(m_fcn, m_edm, up, nfree, ntot, status_cov);
|
|
//status_cov 0=not calculated, 1=approximation, 2=full but forced pos. def., 3=full accurate
|
|
result += 100*status_cov;
|
|
if (opts->hesse_postrun){ //Hesse calculates the 2nd derivative of the likelihood profile to assign the correct systematics
|
|
printCovMatrixStatus(status_cov,false);
|
|
}
|
|
double tmp_cov[varied_params * varied_params];
|
|
//spdlog::info("nvaried params = " << varied_params);;
|
|
minuit_one->mnemat(tmp_cov, varied_params);
|
|
if (opts->weighted_fit && opts->asymptotic){//perform correction according to 1911.01303 Eq. 18
|
|
spdlog::debug("Starting an assymptotic treatment.");
|
|
assert(!(opts->squared_hesse));//TODOCL i think hesse_postrun is fine (and will lead to better estimate for the weighted Hessian)
|
|
|
|
//weighted hesse matrix is available as tmp_cov
|
|
//need to iterate over all events
|
|
//need to evaluate first derivatives at central values
|
|
//events in different pdfs are independently distributed
|
|
//just should count every event separately
|
|
TMatrixDSym V(varied_params, tmp_cov);//weighted covariance matrix
|
|
TMatrixDSym C(varied_params);//, 0.0);//tmp_cov_sq);
|
|
for (unsigned int k=0; k<varied_params; k++){
|
|
for (unsigned int l=0; l<varied_params; l++){
|
|
C(k,l) = 0.0;
|
|
}
|
|
}
|
|
//save best fit values
|
|
spdlog::debug("Saving best fit values.");
|
|
std::vector<std::string> param_names;
|
|
std::vector<double> param_values;
|
|
std::vector<double> eps(varied_params, 1.0e-5);
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0){
|
|
//if the parameter is common, only add name once
|
|
if (!is_common(params.at(i)->get_parameter(j)->get_name()) || is_first(params.at(i)->get_parameter(j)->get_name(), i, j)){
|
|
double v,e;
|
|
std::string par_name = params.at(i)->get_parameter(j)->get_name();
|
|
unsigned int idx = first_minuit_index(par_name);
|
|
minuit_one->GetParameter(idx, v, e);
|
|
param_names.push_back(par_name);
|
|
param_values.push_back(v);
|
|
spdlog::debug(par_name+": {0:f}", v);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//the optimized version varies the parameter in the outer loop and performs the update_cached_normalisation only once for + and -
|
|
//loop over data samples
|
|
|
|
spdlog::debug("Looping over data samples.");
|
|
|
|
std::vector<double> lh_extended_nsig(events.size(),0.0);
|
|
std::vector<double> lh_extended_nbkg(events.size(),0.0);
|
|
std::vector<double> lh_extended_nsigbar(events.size(),0.0);
|
|
std::vector<double> lh_extended_nbkgbar(events.size(),0.0);
|
|
for (unsigned int i = 0; i<events.size(); i++){
|
|
//save derivative for every event, every param
|
|
if (i%10==0) spdlog::trace("Event {0:d}",i);
|
|
std::vector<std::vector<double> > derivs; //first idx varied param, second index event
|
|
for (unsigned int k=0; k<varied_params; k++){
|
|
//only want to do this for params.at(i), other data samples use the other param set
|
|
fcnc::parameter* varied = params.at(i)->get_parameter(param_names.at(k));
|
|
if (!varied){
|
|
std::vector<double> derivs_k(events.at(i)->size(),0.0);
|
|
derivs.push_back(derivs_k);
|
|
continue;
|
|
}
|
|
//set varied parameter + everywhere
|
|
varied->set_value(param_values.at(k)+eps.at(k));
|
|
//cache integrals
|
|
pdfs.at(i)->update_cached_normalization(params.at(i));
|
|
if (opts->use_event_norm) pdfs.at(i)->update_cached_normalization(params.at(i), events.at(i));
|
|
|
|
//TODOCL
|
|
|
|
spdlog::debug("Getting likelihood.");
|
|
likelihood();//this is supposed to only get the extended term, for the + varied parameters, should be fast enough //This doesn't do anything!
|
|
//loop over lh_plus events
|
|
std::vector<double> loglh_plus_k(events.at(i)->size(),0.0);
|
|
for (unsigned int j = 0; j < events.at(i)->size(); j++){
|
|
//for every event do finite differences
|
|
const fcnc::event meas = events.at(i)->at(j);
|
|
double probability = 0.0;
|
|
probability = pdfs.at(i)->prob(params.at(i), meas);
|
|
if (probability > 0.0) loglh_plus_k.at(j) = -TMath::Log(probability);
|
|
//TODOCL
|
|
if (meas.cp_conjugate){ //TODOCL CHECK!! //LEONS COMMENT meas.cp_conjugate corresponds to n_sig !!!!!!!
|
|
loglh_plus_k.at(j) += -TMath::Log(lh_extended_nsig.at(i)+lh_extended_nbkg.at(i));
|
|
}
|
|
else{
|
|
loglh_plus_k.at(j) += -TMath::Log(lh_extended_nsigbar.at(i)+lh_extended_nbkgbar.at(i));
|
|
}
|
|
}
|
|
|
|
//set varied parameter - everywhere
|
|
varied->set_value(param_values.at(k)-eps.at(k));
|
|
//cache integrals
|
|
pdfs.at(i)->update_cached_normalization(params.at(i));
|
|
if (opts->use_event_norm) pdfs.at(i)->update_cached_normalization(params.at(i), events.at(i));
|
|
|
|
//TODOCL
|
|
likelihood();//this is supposed to only get the extended term, for the - varied parameters, should be fast enough
|
|
|
|
//loop over lh_minus events
|
|
std::vector<double> loglh_minus_k(events.at(i)->size(),0.0);
|
|
// for (unsigned int j = 0; j < events.at(i)->size(); j++){
|
|
// //for every event do finite differences
|
|
// const fcnc::event meas = events.at(i)->at(j);
|
|
// double probability = 0.0;
|
|
// probability = pdfs.at(i)->prob(params.at(i), meas);
|
|
// if (probability > 0.0) loglh_minus_k.at(j) = -TMath::Log(probability);
|
|
|
|
// //TODOCL
|
|
// if (meas.cp_conjugate){ //TODOCL CHECK!! //LEONS COMMENT meas.cp_conjugate corresponds to n_sig !!!!!!!
|
|
// loglh_minus_k.at(j) += -TMath::Log(lh_extended_nsig.at(i)+lh_extended_nbkg.at(i));
|
|
// }
|
|
// else{
|
|
// loglh_minus_k.at(j) += -TMath::Log(lh_extended_nsigbar.at(i)+lh_extended_nbkgbar.at(i));
|
|
// }
|
|
// }
|
|
|
|
//calculate derivatives for all events
|
|
spdlog::debug("Calculating derivatives for all events");
|
|
std::vector<double> derivs_k(events.at(i)->size(),0.0);
|
|
for (unsigned int j = 0; j < events.at(i)->size(); j++){
|
|
derivs_k.at(j) = (loglh_plus_k.at(j)-loglh_minus_k.at(j))/(2.0*eps.at(k));
|
|
}
|
|
spdlog::trace("Derivatives:+ "+convert_vector_to_string(derivs_k));
|
|
derivs.push_back(derivs_k);
|
|
|
|
//reset varied parameter everywhere
|
|
varied->set_value(param_values.at(k));
|
|
//cache integrals
|
|
pdfs.at(i)->update_cached_normalization(params.at(i));
|
|
if (opts->use_event_norm){
|
|
pdfs.at(i)->update_cached_normalization(params.at(i), events.at(i));
|
|
}
|
|
|
|
}
|
|
|
|
//afterwards determine derivative via finite differences
|
|
|
|
spdlog::debug("Calculating derivatives for finite differences");
|
|
//and save matrix of squared first derivatives
|
|
for (unsigned int j = 0; j < events.at(i)->size(); j++){
|
|
//for every event do finite differences
|
|
const fcnc::event meas = events.at(i)->at(j);
|
|
for (unsigned int k=0; k<varied_params; k++){
|
|
for (unsigned int l=0; l<varied_params; l++){
|
|
C(k,l) += meas.weight*meas.weight*derivs.at(k).at(j)*derivs.at(l).at(j);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
TMatrixD VCV(V,TMatrixD::kMult,TMatrixD(C,TMatrixD::kMult,V));
|
|
TMatrixDSym VCVsym(varied_params);
|
|
spdlog::debug("printing corrected covariance matrix");
|
|
for (unsigned int i=0; i<varied_params; i++){
|
|
for (unsigned int j=0; j<varied_params; j++) {
|
|
if (spdlog_debug()) std::cout << std::scientific << std::setw(9) << std::setprecision(2) << VCV(i,j) << " ";
|
|
VCVsym(i,j) = 0.5*(VCV(i,j)+VCV(j,i));//no check for now
|
|
tmp_cov[varied_params*j + i] = VCVsym(i,j);//save corrected covariance matrix
|
|
}
|
|
if (spdlog_debug()) std::cout << std::endl;
|
|
}
|
|
}
|
|
|
|
if (opts->weighted_fit && opts->squared_hesse){
|
|
square_weights = true;
|
|
spdlog::info("Starting hesse with squared weights");
|
|
minuit_one->mnhess();
|
|
//check that this went fine
|
|
double m_sq_fcn, m_sq_edm, sq_up;
|
|
int status_sq_cov;
|
|
minuit_one->mnstat(m_sq_fcn, m_sq_edm, sq_up, nfree, ntot, status_sq_cov);
|
|
printCovMatrixStatus(status_sq_cov,true);
|
|
//change status to invalied if squared weighted fit failed
|
|
if (status_cov != status_sq_cov){
|
|
result -= 100*status_cov;
|
|
result += 100*((status_cov < status_sq_cov) ? status_cov : status_sq_cov);
|
|
}
|
|
// minuit_one->mnmatu(1);//this also sets something internal in minuit! getparameter seems to get its error from here!
|
|
double tmp_cov_sq[varied_params * varied_params];
|
|
minuit_one->mnemat(tmp_cov_sq, varied_params);
|
|
//lets use root classes for the matrix inversion
|
|
TMatrixDSym V(varied_params, tmp_cov);//tmp_cov
|
|
TMatrixDSym C(varied_params, tmp_cov_sq);
|
|
double det=0.0;
|
|
C.Invert(&det);
|
|
TMatrixD VCV(V,TMatrixD::kMult,TMatrixD(C,TMatrixD::kMult,V));
|
|
TMatrixDSym VCVsym(varied_params);
|
|
spdlog::info("printing corrected covariance matrix");
|
|
for (unsigned int i=0; i<varied_params; i++){
|
|
for (unsigned int j=0; j<varied_params; j++){
|
|
if (spdlog_debug()){
|
|
std::cout << std::scientific << std::setw(9) << std::setprecision(2) << VCV(i,j) << " ";
|
|
}
|
|
VCVsym(i,j) = 0.5*(VCV(i,j)+VCV(j,i));//no check for now
|
|
tmp_cov[varied_params*j + i] = VCVsym(i,j);//save corrected covariance matrix
|
|
}
|
|
if (spdlog_debug()) std::cout << std::endl;
|
|
}
|
|
square_weights = false;
|
|
}
|
|
|
|
//to get the correct latex names
|
|
std::vector<std::string> varied_names;
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0){//if the parameter is common, only add name once
|
|
if (!is_common(params.at(i)->get_parameter(j)->get_name()) || is_first(params.at(i)->get_parameter(j)->get_name(), i, j)){
|
|
varied_names.push_back(params.at(i)->get_parameter(j)->get_description());
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//printing error matrix
|
|
if (spdlog_debug()){
|
|
spdlog::debug("printing error matrix");
|
|
for (unsigned int x = 0; x < varied_params; x++) std::cout << "\tvar" << x << " ";
|
|
std::cout << std::endl;
|
|
for (unsigned int y = 0; y < varied_params; y++){
|
|
std::cout << "var" << y;
|
|
for (unsigned int x = 0; x < varied_params; x++){
|
|
std::cout << " " << std::scientific << std::setw(9) << std::setprecision(2) << tmp_cov[varied_params*y + x];
|
|
}
|
|
std::cout << std::endl;
|
|
}
|
|
}
|
|
//calculating correlation matrix from error matrix
|
|
double cor[varied_params*varied_params];
|
|
for (unsigned int x = 0; x < varied_params; x++){
|
|
for (unsigned int y = 0; y < varied_params; y++){
|
|
cor[varied_params*y + x] = tmp_cov[varied_params*y + x]
|
|
/(sqrt(fabs(tmp_cov[varied_params*x + x]))*sqrt(fabs(tmp_cov[varied_params*y + y])));
|
|
}
|
|
}
|
|
//printing correlation matrix
|
|
|
|
if (index != ""){
|
|
spdlog::info("printing correlation matrix");
|
|
myFile << "Run: " << opts->run << std::endl;
|
|
myFile << "\\begin{tabular}{c";
|
|
for (unsigned int x = 0; x < varied_params; x++){
|
|
myFile << "c";
|
|
}
|
|
myFile << "} \\hline" << std::endl;
|
|
myFile << " ";
|
|
for (unsigned int x = 0; x < varied_params; x++){
|
|
myFile << " & $" << varied_names.at(x) << "$";
|
|
}
|
|
myFile << "\\\\ \\hline" << std::endl;
|
|
|
|
bool upper_half = true;
|
|
for (unsigned int y = 0; y < varied_params; y++){
|
|
myFile << std::setw(20) << "$" << varied_names.at(y) << "$";
|
|
for (unsigned int x = 0; x < varied_params; x++){
|
|
if (x >= y || !upper_half){
|
|
if (fabs(cor[varied_params*y + x]) >= 0.005){
|
|
myFile << " & " << std::fixed << std::setw(5) << std::setprecision(2) << cor[varied_params*y + x];
|
|
}
|
|
else myFile << " & " << std::setw(5) << "-";
|
|
}
|
|
else myFile << " & " << std::setw(5) << " ";
|
|
}
|
|
myFile << " \\\\ " << std::endl;
|
|
}
|
|
myFile <<"\\hline \\end{tabular}" << std::endl;
|
|
}
|
|
|
|
|
|
//saves correlations for all parameters
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j=0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0){
|
|
params.at(i)->get_parameter(j)->correlations.clear();
|
|
unsigned int apary = 0;
|
|
for (unsigned int y=0; y < varied_params; y++){
|
|
if (varied_names.at(y) == params.at(i)->get_parameter(j)->get_description()){
|
|
apary = y;
|
|
break;
|
|
}
|
|
}
|
|
for (unsigned int x=0; x < varied_params; x++)
|
|
params.at(i)->get_parameter(j)->correlations.push_back(cor[varied_params*apary + x]);
|
|
}
|
|
}
|
|
}
|
|
|
|
//minos error analysis
|
|
if (opts->minos_errors && (result%100) == 0){//migrad is successfull
|
|
spdlog::info("Starting Minos error analysis");
|
|
std::vector<int> indices;
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0 && params.at(i)->get_parameter(j)->get_minos()){
|
|
std::string par_name = params.at(i)->get_parameter(j)->get_name();
|
|
unsigned int idx = 0;
|
|
if(is_common(par_name) || opts->use_individual_param_names) idx = first_minuit_index(par_name);
|
|
else idx = nth_minuit_index(par_name, i+1);
|
|
indices.push_back(idx);
|
|
}
|
|
}
|
|
}
|
|
if (indices.size() > 0){//run minos just for subset of parameters
|
|
spdlog::info("MINOS SUBSET");
|
|
int errorcode;
|
|
/*
|
|
double migrad_options[indices.size()+1];
|
|
migrad_options[0] = 100000000.0;
|
|
for (unsigned int i = 0; i < indices.size(); i++)
|
|
//migrad_options[i+1] = indices.at(i);
|
|
migrad_options[i+1] = indices.at(i)+1;
|
|
minuit_one->mnexcm("MINO", migrad_options, indices.size()+1, errorcode);
|
|
*/
|
|
for (unsigned int i = 0; i < indices.size(); i++){
|
|
double migrad_options[2];
|
|
migrad_options[0] = 100000000.0;
|
|
migrad_options[1] = indices.at(i)+1;
|
|
minuit_one->mnexcm("MINO", migrad_options, 2, errorcode);
|
|
}
|
|
}
|
|
else minuit_one->mnmnos();//do minos error analysis on all parameters
|
|
|
|
//write MINOS errors to parameters
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
double e_up, e_down, parab, gcc;
|
|
if (params.at(i)->get_parameter(j)->get_step_size() != 0.0){
|
|
std::string par_name = params.at(i)->get_parameter(j)->get_name();
|
|
unsigned int idx = 0;
|
|
if(is_common(par_name) || opts->use_individual_param_names) idx = first_minuit_index(par_name);
|
|
else idx = nth_minuit_index(par_name, i+1);
|
|
minuit_one->mnerrs(idx, e_up, e_down, parab, gcc);
|
|
params.at(i)->get_parameter(j)->set_error_up(e_up);
|
|
params.at(i)->get_parameter(j)->set_error_down(e_down);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//finally extract fitted parameters from minuit
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
double v,e;
|
|
std::string par_name = params.at(i)->get_parameter(j)->get_name();
|
|
unsigned int idx;
|
|
if(is_common(par_name) || opts->use_individual_param_names) idx = first_minuit_index(par_name);
|
|
else idx = nth_minuit_index(par_name, i + 1);
|
|
|
|
spdlog::trace("[PDF{0:d}]\tPAR{1:d}:\t{2:s}\tTMinuit_idx: {3:d}", i, j, par_name, idx);
|
|
|
|
minuit_one->GetParameter(idx, v, e);
|
|
params.at(i)->get_parameter(j)->set_value(v);
|
|
params.at(i)->get_parameter(j)->set_error(e);
|
|
|
|
|
|
if (opts->weighted_fit){
|
|
//this method is approximate for parameters with limited range!
|
|
int iidx = minuit_one->fNiofex[idx] - 1;
|
|
if (iidx>=0){
|
|
e = sqrt(fabs(tmp_cov[iidx*varied_params+iidx]));
|
|
params.at(i)->get_parameter(j)->set_error(e);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (spdlog_trace()){
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
params.at(i)->print_parameters();
|
|
}
|
|
}
|
|
|
|
if (opts->project){
|
|
for (unsigned int i = 0; i<plots.size(); i++)
|
|
plots.at(i)->plot_data(pdfs.at(i), params.at(i), events.at(i), index);
|
|
}
|
|
spdlog::info("Fit procedure finished");
|
|
}
|
|
|
|
//Close Latex file
|
|
if (index != "") myFile.close();
|
|
|
|
return result;
|
|
}
|
|
|
|
void fcnc::fitter::define_parameter(int i, const parameter & p){
|
|
if (p.get_unlimited()){
|
|
minuit_one->DefineParameter(i, p.get_mn_name(), p(), p.get_step_size(), 0.0, 0.0);
|
|
}
|
|
else{
|
|
minuit_one->DefineParameter(i, p.get_mn_name(), p(), p.get_step_size(), p.get_min(), p.get_max());
|
|
}
|
|
}
|
|
|
|
void fcnc::fitter::reset_param_values(){
|
|
spdlog::debug("Resetting parameter values");
|
|
minuit_one->mncler();
|
|
|
|
unsigned int start_param = 0;
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if (opts->reset_start){
|
|
params.at(i)->get_parameter(j)->reset_start();
|
|
}
|
|
std::string par_name = params.at(i)->get_parameter(j)->get_name();
|
|
if (!is_common(par_name) || is_first(par_name, i, j)){
|
|
define_parameter(start_param+j, *(params.at(i)->get_parameter(j)));
|
|
}
|
|
}
|
|
start_param += params.at(i)->nparameters();
|
|
}
|
|
}
|
|
|
|
double fcnc::fitter::likelihood(){
|
|
lh_tot = 0.0;
|
|
//this is called on every parameter change
|
|
//this can never hurt, also needed in case of weighted/unfolded fit
|
|
for (unsigned int i = 0; i < pdfs.size(); i++){
|
|
pdfs.at(i)->update_cached_normalization(params.at(i));//should do this, but should buffer f1, update f1 on fit command
|
|
if(opts->use_event_norm){
|
|
pdfs.at(i)->update_cached_normalization(params.at(i), events.at(i));
|
|
}
|
|
}
|
|
|
|
//threading the likelihood determination on multiple CPU cores
|
|
if(opts->ncores > 1){
|
|
std::vector<std::thread*> threads;
|
|
for(unsigned int i = 0; i < opts->ncores; i++){
|
|
std::thread* t = new std::thread(std::bind( &fitter::likelihood_thread, this, i));
|
|
threads.push_back(t);
|
|
}
|
|
for(unsigned int i = 0; i < opts->ncores; i++) threads[i]->join();
|
|
for(unsigned int i = 0; i < opts->ncores; i++) delete threads[i];
|
|
}
|
|
else{
|
|
likelihood_thread(0);
|
|
}
|
|
|
|
if (spdlog_trace()){
|
|
unsigned int start_param = 0;
|
|
for(unsigned int i = 0; i < params.size(); i++){
|
|
for(unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
if(params.at(i)->get_parameter(j)->get_step_size() != 0.0){
|
|
std::string par_name = params.at(i)->get_parameter(j)->get_name();
|
|
if(!is_common(par_name) || is_first(par_name, i, j)){
|
|
if(!params.at(i)->get_parameter(j)->is_blind()) std::cout << par_name << ":" << params.at(i)->get_parameter(j)->get_value() << " ";
|
|
else std::cout << par_name << ":" << "XXX" << " ";
|
|
}
|
|
}
|
|
}
|
|
start_param += params.at(i)->nparameters();
|
|
}
|
|
std::cout << std::endl;
|
|
}
|
|
|
|
//there is an additional term for parameters which have been measured previously
|
|
for (unsigned int i = 0; i < params.size(); i++){
|
|
for (unsigned int j = 0; j < params.at(i)->nparameters(); j++){
|
|
const parameter* param = params.at(i)->get_parameter(j);
|
|
if (param->get_gaussian_constraint()){
|
|
if (!is_common(param->get_name()) || is_first(param->get_name(), i, j)){
|
|
double mean = param->get_previous_measurement();
|
|
double x = param->get_value();
|
|
double dup = fabs(param->get_previous_error_high());
|
|
double ddown = fabs(param->get_previous_error_low());
|
|
double sigma = (x>mean ? dup : ddown);
|
|
//spdlog::trace("mean {0:f}\tx {1:f}\tdup {2:f}\tddown {3:f}\tsigma {4:f}",mean,x,dup,ddown,sigma);
|
|
|
|
if (x>mean+dup) sigma = dup;
|
|
else if (x<mean-ddown) sigma = ddown;
|
|
else sigma = ddown + (dup-ddown)*(x-(mean-ddown))/(dup+ddown);
|
|
lh_tot += (-2.0*(-(x-mean)*(x-mean)/(2.0*sigma*sigma)));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (opts->chisquaredstudy){
|
|
double chi2 = 0.0;
|
|
for (unsigned int i = 0; i < pdfs.size(); i++) chi2 += pdfs.at(i)->chi2(params.at(i));
|
|
lh_tot += chi2;
|
|
}
|
|
if (opts->extended_ml){
|
|
double lh_ext = 0.0;
|
|
for(unsigned int i = 0; i < params.size(); i++){
|
|
double nobs = 0.0;
|
|
parameter* pnsig = params.at(i)->get_parameter("n_sig");
|
|
parameter* pnbkg = params.at(i)->get_parameter("n_bkg");
|
|
assert(pnsig && pnbkg);
|
|
double nsig = pnsig->get_value();
|
|
double nbkg = pnbkg->get_value();
|
|
if(opts->weighted_fit){
|
|
for(unsigned int j = 0; j<events.at(i)->size(); j++){
|
|
nobs += events.at(i)->at(j).weight;
|
|
}
|
|
}
|
|
else{
|
|
nobs = events.at(i)->size();
|
|
}
|
|
double prob = TMath::Poisson(nsig+nbkg, nobs);
|
|
if(prob == 0.0) prob = 1e-100;
|
|
lh_ext += -2.0*TMath::Log( prob );
|
|
if (std::isinf(lh_ext) || std::isnan(lh_ext)){
|
|
std::cout << "[PDF" << i << "]\t\tlh_ext = " << lh_ext << "\t prob = " << prob << "\t -2log(prob) = " << -2.0*TMath::Log( prob ) << std::endl;
|
|
std::cout << "[PDF" << i << "]\t\tnsig = " << nsig << "\t n_bkg = " << nbkg << "\t nobs = " << nobs << std::endl;
|
|
assert(0);
|
|
}
|
|
}
|
|
if (std::isinf(lh_ext)){
|
|
spdlog::error("{0:0.16f}", lh_ext);
|
|
spdlog::error("lh (w/o lh_ext) = {0:0.16f}", lh_tot);
|
|
assert(0);
|
|
}
|
|
lh_tot += lh_ext;
|
|
}
|
|
|
|
if (std::isinf(lh_tot || std::isnan(lh_tot))){
|
|
spdlog::error("lh = {0:0.16f}", lh_tot);
|
|
assert(0);
|
|
}
|
|
if (spdlog_trace()) std::cout << "[debug]\t" << std::setprecision(16) << "lh = " << lh_tot << "\r" << std::flush;
|
|
return lh_tot;
|
|
}
|
|
|
|
void fcnc::fitter::likelihood_thread(int thread_id){
|
|
|
|
double probability = 0.0;
|
|
event meas;
|
|
long double result = 0.0;
|
|
|
|
for (unsigned int i = 0; i<events.size(); i++){
|
|
unsigned int min, max;
|
|
unsigned int evts_per_thread = events.at(i)->size()/opts->ncores;
|
|
//is too small for last vector!
|
|
min = evts_per_thread * thread_id;
|
|
max = evts_per_thread * (thread_id+1);
|
|
int last_thread = opts->ncores - 1;
|
|
if (thread_id == last_thread){//last thread needs to take all posibly remaining events
|
|
max = events.at(i)->size();
|
|
}
|
|
std::vector<long double> vec;
|
|
vec.reserve(max - min);
|
|
for (unsigned int j = min; j < max; j++){
|
|
meas = events.at(i)->at(j);
|
|
probability = pdfs.at(i)->prob(params.at(i), meas);
|
|
if (probability < 0.0 || std::isnan(probability) || std::isinf(probability)){
|
|
std::lock_guard<std::mutex> lock(fitter_mutex);
|
|
spdlog::info("Event probability is {0:d}", probability);
|
|
print_event(meas);
|
|
params.at(i)->print_parameters();
|
|
assert(0);
|
|
}
|
|
if (probability == 0.0){
|
|
probability = 1.0;//ignore event
|
|
}
|
|
if (probability != 0.0){
|
|
if (opts->weighted_fit){
|
|
if (square_weights) vec.push_back(-2.0 * meas.weight * meas.weight * TMath::Log(probability) - empirical_constant);//todo: does the empirical constant pose a problem here?
|
|
else vec.push_back(-2.0 * meas.weight * TMath::Log(probability) - empirical_constant);
|
|
}
|
|
else vec.push_back(-2.0 * TMath::Log(probability) - empirical_constant);
|
|
}
|
|
else if (!std::isnan(probability)){
|
|
spdlog::error("Event probability is {0:d}", probability);
|
|
print_event(meas);
|
|
assert(0);
|
|
}
|
|
else{//is inf?
|
|
vec.push_back(8.0);//is inf?//use something like prob max? log(p_min)
|
|
}
|
|
}
|
|
result += add_results<std::vector<long double>::iterator>(vec.begin(), vec.end());
|
|
}
|
|
if (opts->ncores > 1){
|
|
std::lock_guard<std::mutex> lock(fitter_mutex);
|
|
lh_tot += result;
|
|
}
|
|
else{
|
|
lh_tot += result;
|
|
}
|
|
}
|
|
|
|
//boost::mutex fitter::fitter_mutex;
|
|
std::mutex fcnc::fitter::fitter_mutex;
|
|
|
|
fcnc::fitter* fcnc::fitter::minuit_one_fitter;
|