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.
953 lines
38 KiB
953 lines
38 KiB
//Renata Kopecna
|
|
|
|
#include <fstream>
|
|
#include <iostream>
|
|
#include <sstream> // std::istringstream
|
|
|
|
#include <bu2kstarmumu_parameters.hh>
|
|
#include <fitter.hh>
|
|
#include <helpers.hh>
|
|
#include <paths.hh>
|
|
|
|
#include <TFile.h>
|
|
#include <TMath.h>
|
|
#include <TChain.h>
|
|
|
|
#include <spdlog.h>
|
|
|
|
//fixConstr constructor
|
|
fixConstr::fixConstr(bool b_fix,bool b_constr){
|
|
if (b_fix && b_constr){ //Check if both constrain and fix are set to true
|
|
//If yes, set constrain to false without modifying the input bools
|
|
spdlog::warn("Cannot both fix and constrain a parameter.\n The parameter is fixed, setting constrain to false.");
|
|
spdlog::warn("Check your options.");
|
|
fix = true;
|
|
constrain = false;
|
|
return;
|
|
}
|
|
fix = b_fix;
|
|
constrain = b_constr;
|
|
return;
|
|
}
|
|
|
|
using namespace fcnc;
|
|
|
|
void bu2kstarmumu_parameters::use_default_bkg(){
|
|
//Set all backgrounds to be flat
|
|
cbkgctl0.init(1.0, -1.0, 1.0, 0.0);
|
|
cbkgctl1.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctl2.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctl3.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctl4.init(0.0, -1.0, 1.0, 0.0);
|
|
|
|
cbkgctk0.init(1.0, -1.0, 1.0, 0.0);
|
|
cbkgctk1.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctk2.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctk3.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctk4.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctk5.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgctk6.init(0.0, -1.0, 1.0, 0.0);
|
|
|
|
cbkgphi0.init(1.0, -1.0, 1.0, 0.0);
|
|
cbkgphi1.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgphi2.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgphi3.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgphi4.init(0.0, -1.0, 1.0, 0.0);
|
|
|
|
cbkgmkpi0.init(1.0, -1.0, 1.0, 0.0);
|
|
cbkgmkpi1.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgmkpi2.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgmkpi3.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgmkpi4.init(0.0, -1.0, 1.0, 0.0);
|
|
|
|
cswavemkpi0.init(1.0, -1.0, 1.0, 0.0);
|
|
cswavemkpi1.init(0.0, -1.0, 1.0, 0.0);
|
|
cswavemkpi2.init(0.0, -1.0, 1.0, 0.0);
|
|
cswavemkpi3.init(0.0, -1.0, 1.0, 0.0);
|
|
cswavemkpi4.init(0.0, -1.0, 1.0, 0.0);
|
|
|
|
cbkgp20.init(1.0, -1.0, 1.0, 0.0);
|
|
cbkgp21.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgp22.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgp23.init(0.0, -1.0, 1.0, 0.0);
|
|
cbkgp24.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::use_default_observables(){
|
|
if (opts->fit_pprimes){
|
|
Fl.init (0.7, 0.0, 2.0*PAR_ANG_RANGE, 0.0);
|
|
P1.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P2.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P3.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P4.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P5.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P6.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P8.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
}
|
|
else{
|
|
Fl.init (0.7, 0.0, 2.0*PAR_ANG_RANGE, 0.0);
|
|
S1s.init(0.7, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S3.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S4.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S5.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
Afb.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S6s.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S7.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S8.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S9.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
}
|
|
return;
|
|
}
|
|
//This function is used in the parameters constructor!
|
|
void bu2kstarmumu_parameters::use_default(){
|
|
Fl.init (0.7, 0.0, 2.0*PAR_ANG_RANGE, 0.0);
|
|
S1s.init(0.7, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S3.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S4.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S5.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
Afb.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S6s.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S7.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S8.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
S9.init (0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
|
|
P1.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P2.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P3.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P4.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P5.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P6.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
P8.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
|
|
//S-wave parameters
|
|
FS.init (0.0, 0.0, 2*PAR_ANG_RANGE, 0.0);
|
|
SS1.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
SS2.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
SS3.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
SS4.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
SS5.init(0.0, -PAR_ANG_RANGE, PAR_ANG_RANGE, 0.0);
|
|
|
|
//signal fraction
|
|
f_sig.init(PAR_N_SIG/(PAR_N_SIG+PAR_N_BKG), 0.0, 1.0, 0.0);
|
|
n_sig.init(PAR_N_SIG, 0.0, 1.0e+6, 0.0);
|
|
n_bkg.init(PAR_N_BKG, 0.0, 1.0e+6, 0.0);
|
|
//mass parameters
|
|
m_b.init (PDGMASS_B, B_MASS_LOW, B_MASS_HIGH, 0.0);
|
|
m_res_1.init (DOUBLE_CB ? 1.0 : 0.5, 0.0, 1.0, 0.0);
|
|
m_sigma_1.init(PAR_SIGMA, PAR_SIGMA_LOW, PAR_SIGMA_HIGH, 0.0);
|
|
m_sigma_2.init(DOUBLE_CB ? 0.0 : PAR_SIGMA, PAR_SIGMA_LOW, PAR_SIGMA_HIGH, 0.0);
|
|
//crystal ball
|
|
m_scale.init(1.0, 0.0, 2.0, 0.0);
|
|
alpha_1.init(1.0, 0.1, 10.0, 0.0);
|
|
n_1.init (1.0, 0.1, 10.0, 0.0);
|
|
alpha_2.init(1.0, 0.1, 10.0, 0.0);
|
|
n_2.init (1.0, 0.1, 10.0, 0.0);
|
|
|
|
fm_tau.init(FIX_FM ? 1.0 : 0.5, 0.0, 1.0, 0.0);
|
|
m_tau.init(PAR_TAU, PAR_TAU/PAR_TAU_SCALE, PAR_TAU*PAR_TAU_SCALE, 0.0);
|
|
m_tau_2.init(PAR_TAU, PAR_TAU/PAR_TAU_SCALE, PAR_TAU*PAR_TAU_SCALE, 0.0);
|
|
m_lambda.init(PAR_LAMBDA, PAR_LAMBDA/PAR_LAMBDA_SCALE, PAR_LAMBDA*PAR_LAMBDA_SCALE, 0.0);
|
|
m_lambda_2.init(PAR_LAMBDA, PAR_LAMBDA/PAR_LAMBDA_SCALE, PAR_LAMBDA*PAR_LAMBDA_SCALE, 0.0);
|
|
|
|
//Welp, this will end up in Jpsi, but whatever :)
|
|
eff_q2.init(bin_center(Q2_MIN_RANGE,Q2_MAX_RANGE), Q2_MIN_RANGE, Q2_MAX_RANGE, 0.0);
|
|
|
|
//TODO: figure out what these are
|
|
asphase.init(MY_PI, 0.0, 2.0*MY_PI, 0.0);
|
|
a.init(1.0, 0.0, 10.0, 0.0);
|
|
r.init(1.0, 0.0, 10.0, 0.0);
|
|
gammakstar.init(PAR_KSTAR_WIDTH/1.0e3, 0.01, PAR_KSTAR_WIDTH/0.5e3, 0.0); //With of K*
|
|
mkstar.init(PDGMASS_K_STAR_PLUS/1.0e3, (K_ONE_PLUS-PAR_KSTAR_WIDTH)/1.0e3, (K_ONE_PLUS+PAR_KSTAR_WIDTH)/1.0e3, 0.0);
|
|
//This is a K1 or K something
|
|
mkstarplus.init(K_ONE_PLUS/1.0e3, (K_ONE_PLUS-PAR_K1_WIDTH)/1.0e3, (K_ONE_PLUS+PAR_K1_WIDTH)/1.0e3, 0.0);
|
|
gammakstarplus.init(PAR_K1_WIDTH/1.0e3, 0.01, PAR_K1_WIDTH/0.5e3, 0.0); //Width of K1
|
|
|
|
|
|
nthreshold.init(PAR_NTRESHOLD, 1.0, 10.0, 0.0);
|
|
R.init(1.6, 0.0, 10.0, 0.0);
|
|
|
|
//Some other fancy resonances //TODO ask Eluned
|
|
mf800.init(0.682, 0.1, 1.0, 0.0);
|
|
gammaf800.init(0.547, 0.1, 1.0, 0.0);
|
|
f800mag.init(0.1, 0.0, 1.0, 0.0);
|
|
f800phase.init(0.5*MY_PI, -2.0*MY_PI, +2.0*MY_PI, 0.0);
|
|
|
|
//Set all backgrounds to be flat
|
|
use_default_bkg();
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::load_param_values(std::string filename){
|
|
|
|
std::ifstream file(filename.c_str());;
|
|
std::string line = "";
|
|
if (file.is_open()){
|
|
spdlog::info("Parameter values are being read from file "+filename+".");
|
|
}
|
|
else{
|
|
spdlog::error("Could not load values from file " + filename);;
|
|
assert(0);
|
|
}
|
|
|
|
UInt_t p_index;
|
|
std::string p_name;
|
|
Float_t p_value;
|
|
Float_t p_error;
|
|
//loop over lines in file and read values and errors of all parameters
|
|
while(true){ //Jesus christ, another while(true). Shoot me, please.
|
|
getline(file, line);
|
|
spdlog::debug("Line: "+line);;
|
|
if(file.eof()) break;
|
|
std::istringstream istr(line);
|
|
istr >> p_index;
|
|
istr >> p_name;
|
|
istr >> p_value;
|
|
istr >> p_error;
|
|
|
|
//sanity check to make sure the parameter index fits to the parameter name:
|
|
if(p_name != this->get_parameter(p_index)->get_name()){
|
|
spdlog::critical("Trying to assign values of var=" + p_name + " to parameter=" + this->get_parameter(p_index)->get_name());;
|
|
assert(p_name == this->get_parameter(p_index)->get_name());
|
|
}
|
|
|
|
//assign value and error from file to parameters
|
|
this->get_parameter(p_index)->set_values_n_errors(p_value, p_error);
|
|
}
|
|
file.close();
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::load_param(std::string filename, std::string parname){
|
|
|
|
std::ifstream file(filename.c_str());;
|
|
std::string line = "";
|
|
if (!file.is_open()){
|
|
spdlog::error("\tCould not open file " + filename);
|
|
assert(0);
|
|
}
|
|
|
|
UInt_t p_index;
|
|
std::string p_name;
|
|
Float_t p_value;
|
|
Float_t p_error;
|
|
//loop over lines in file and read values and errors of all parameters
|
|
while(true){ //Dear kids,
|
|
getline(file, line);
|
|
spdlog::debug("Line: " + line);
|
|
if(file.eof()) break; //never ever write code like this. Even at gunpoint. Don't.
|
|
std::istringstream istr(line);
|
|
istr >> p_index;
|
|
istr >> p_name;
|
|
istr >> p_value;
|
|
istr >> p_error;
|
|
|
|
//continue if param name is identical
|
|
if(p_name != parname) continue;
|
|
|
|
spdlog::debug("Load parameter: " + p_name + "={0:f}", p_value);
|
|
|
|
//sanity check to make sure the parameter index fits to the parameter name:
|
|
assert(p_name == this->get_parameter(p_index)->get_name());
|
|
|
|
//assign value and error from file to parameters
|
|
this->get_parameter(p_index)->set_values_n_errors(p_value, p_error);
|
|
}
|
|
file.close();
|
|
return;
|
|
}
|
|
|
|
double get_param_valueError_from_rootfile(std::string fileName, std::string name, int PDF, int bin, bool getError){
|
|
//Retuns a value for given parametr from a given rootfile
|
|
|
|
//Open file
|
|
spdlog::info("Opening " + fileName);
|
|
TFile* file = new TFile(fileName.c_str(), "READ");
|
|
|
|
if (!file->GetListOfKeys()->Contains(name.c_str())){
|
|
spdlog::critical("Wrong tree name " + name + "! Abort.");
|
|
assert(0);
|
|
}
|
|
|
|
//Get the tree and set the branches to active
|
|
TTree* tree = (TTree*)file->Get(name.c_str());
|
|
tree->SetBranchStatus("*",1); //Activate all branches
|
|
|
|
//The root file is saved in a way that first bins in first pdf is saved
|
|
//then the next pdf is saved, bin by bin
|
|
//Meaning if I want to load event from bin 2 and pdf 1
|
|
//tree->GetEntry(2)
|
|
//If I want to load event from bin 2 and pdf 2, one needs to do
|
|
//tree->GetEntry(totBins+2)
|
|
//PDFs are named according to the Run! So if only Run 2 is fitted, the returned pdf is 2
|
|
//This implementation is actually useful for checks
|
|
//Therefore, starting entry is equal to int bin
|
|
int entry_position = bin;
|
|
|
|
//Read the q2 bins
|
|
int totBins = DEFAULT_TREE_INT;
|
|
int pdf_idx = DEFAULT_TREE_INT;
|
|
int bin_idx = DEFAULT_TREE_INT;
|
|
double value = DEFAULT_TREE_VAL;
|
|
double error = DEFAULT_TREE_ERR;
|
|
int migrad = DEFAULT_TREE_INT;
|
|
|
|
tree->SetBranchAddress("totBins", &totBins);
|
|
tree->SetBranchAddress("pdf", &pdf_idx);
|
|
tree->SetBranchAddress("bin", &bin_idx);
|
|
tree->SetBranchAddress("value",&value);
|
|
tree->SetBranchAddress("migrad",&migrad);
|
|
tree->SetBranchAddress("error",&error);
|
|
|
|
tree->GetEntry(entry_position);
|
|
|
|
//Check whether we got the right bin position
|
|
if (bin_idx != bin){
|
|
spdlog::critical("Something went very wrong for "+name);
|
|
spdlog::critical("The read pdf {0:d} from position {1:d} doesn't agree with the wanted bin {2:d}", pdf_idx, entry_position, bin);
|
|
assert(0);
|
|
}
|
|
while (pdf_idx != PDF){
|
|
entry_position += totBins;
|
|
if (entry_position >= tree->GetEntries()){
|
|
spdlog::critical("Required PDF {0:d} for " + name + " not found!",PDF);
|
|
assert(0);
|
|
}
|
|
tree->GetEntry(entry_position);
|
|
}
|
|
|
|
|
|
//Check if migrad makes sense
|
|
if (migrad != 0){
|
|
spdlog::warn("Careful! You are reading a fit that had failed migrad for " + name + "!");
|
|
printMigradStatus(migrad);
|
|
}
|
|
|
|
|
|
|
|
file->Close();
|
|
spdlog::debug("Loaded " + name + " value is {0:f}", value);
|
|
|
|
return getError ? error : value;
|
|
}
|
|
|
|
double get_param_value_from_rootfile(std::string fileName, std::string name, int PDF, int bin){
|
|
return get_param_valueError_from_rootfile(fileName, name, PDF, bin, false);
|
|
}
|
|
double get_param_error_from_rootfile(std::string fileName, std::string name, int PDF, int bin){
|
|
return get_param_valueError_from_rootfile(fileName, name, PDF, bin, true);
|
|
}
|
|
|
|
int bu2kstarmumu_parameters::get_param_from_rootfile(std::string fileName, std::vector<std::string> names, int PDF, int bin, fixConstr FC){
|
|
|
|
//Open file
|
|
spdlog::info("Opening " + fileName);
|
|
TFile* file = new TFile(fileName.c_str(), "READ");
|
|
|
|
for (auto name: names){
|
|
if (!file->GetListOfKeys()->Contains(name.c_str())){
|
|
spdlog::critical("Wrong tree name " + name + "! Abort.");
|
|
assert(0);
|
|
}
|
|
|
|
//Get the tree and set the branches to active
|
|
TTree* tree = (TTree*)file->Get(name.c_str());
|
|
tree->SetBranchStatus("*",1); //Activate all branches
|
|
|
|
//The root file is saved in a way that first bins in first pdf is saved
|
|
//then the next pdf is saved, bin by bin
|
|
//Meaning if I want to load event from bin 2 and pdf 1
|
|
//tree->GetEntry(2)
|
|
//If I want to load event from bin 2 and pdf 2, one needs to do
|
|
//tree->GetEntry(totBins+2)
|
|
//PDFs are named according to the Run! So if only Run 2 is fitted, the returned pdf is 2
|
|
//This implementation is actually useful for checks
|
|
//Therefore, starting entry is equal to int bin
|
|
int entry_position = bin;
|
|
|
|
//Read the q2 bins
|
|
int totBins = DEFAULT_TREE_INT;
|
|
int pdf_idx = DEFAULT_TREE_INT;
|
|
int bin_idx = DEFAULT_TREE_INT;
|
|
tree->SetBranchAddress("totBins", &totBins);
|
|
tree->SetBranchAddress("pdf", &pdf_idx);
|
|
tree->SetBranchAddress("bin", &bin_idx);
|
|
|
|
//Now get all the parameters
|
|
int migrad = DEFAULT_TREE_INT;
|
|
int param_index = DEFAULT_TREE_INT;
|
|
|
|
double value = DEFAULT_TREE_VAL;
|
|
double error = DEFAULT_TREE_ERR;
|
|
double error_up = DEFAULT_TREE_ERR;
|
|
double error_down = DEFAULT_TREE_ERR;
|
|
double start_value = DEFAULT_TREE_VAL;
|
|
double prev_value = DEFAULT_TREE_VAL;
|
|
double prev_error = DEFAULT_TREE_ERR;
|
|
|
|
tree->SetBranchAddress("value",&value);
|
|
tree->SetBranchAddress("error",&error);
|
|
tree->SetBranchAddress("error_up",&error_up);
|
|
tree->SetBranchAddress("error_down",&error_down);
|
|
|
|
tree->SetBranchAddress("start_value",&start_value); //Not really used: TODO
|
|
tree->SetBranchAddress("prev_value",&prev_value); //Not really used: TODO
|
|
tree->SetBranchAddress("prev_error",&prev_error); //Not really used: TODO
|
|
|
|
tree->SetBranchAddress("index",¶m_index);
|
|
tree->SetBranchAddress("migrad",&migrad);
|
|
|
|
tree->GetEntry(entry_position);
|
|
|
|
//Check whether the file contains correct numbers of bins
|
|
//If totBins are 1, let's assume this is fine and we are reading from ReferenceChannel
|
|
if (unsigned (totBins) != opts->TheQ2binsmin.size() && totBins != 1){
|
|
spdlog::warn("Wrong number of bins for " + name + "!");
|
|
spdlog::warn("totBins = {0:d}\tnQ2bins= {1:d}",totBins,opts->TheQ2binsmin.size());
|
|
}
|
|
if (bin_idx != bin){
|
|
spdlog::critical("Something went very wrong for "+name);
|
|
spdlog::critical("The read pdf {0:d} from position {1:d} doesn't agree with the wanted bin {2:d}", pdf_idx, entry_position, bin);
|
|
assert(0);
|
|
}
|
|
while (pdf_idx != PDF){
|
|
entry_position += totBins;
|
|
if (entry_position >= tree->GetEntries()){
|
|
spdlog::critical("Required PDF {0:d} for " + name + " not found!",PDF);
|
|
assert(0);
|
|
}
|
|
tree->GetEntry(entry_position);
|
|
}
|
|
|
|
//Check if migrad makes sense
|
|
if (migrad != 0){
|
|
spdlog::warn("Careful! You are reading a fit that had failed migrad for " + name + "!");
|
|
printMigradStatus(migrad);
|
|
}
|
|
|
|
//if opts->hesse is non-active, use error for both error_up and error_down
|
|
if(error_down == DEFAULT_TREE_ERR) error_down = error;
|
|
if(error_up == DEFAULT_TREE_ERR) error_up = error;
|
|
|
|
//Init parameter
|
|
double stepsize = GetParamStepsize(name); //TODO check
|
|
if (FC.fix) stepsize = 0.0;
|
|
else if (stepsize < 0.001) stepsize = 0.01; //Protect zeroes //TODO: possibly add some protection against very very small numbers, but if it works as it is, no need to make it fancier
|
|
this->get_parameter(param_index)->set_previous_measurement(prev_value);
|
|
if (FC.constrain){
|
|
if (error_down==0) error_down = 0.01;
|
|
if (error_up==0) error_up = 0.01;
|
|
}
|
|
std::vector<double> param_range = GetParamRange(name);
|
|
this->get_parameter(param_index)->init(value, param_range.front(), param_range.back(), stepsize,error_down,error_up,FC.constrain);
|
|
}
|
|
|
|
file->Close();
|
|
|
|
return 0;
|
|
}
|
|
|
|
int bu2kstarmumu_parameters::fix_param_from_rootfile(std::string fileName,std::vector<std::string> names, int PDF, int bin){
|
|
return get_param_from_rootfile(fileName, names, PDF, bin, fixConstr(true,false));
|
|
}
|
|
|
|
int bu2kstarmumu_parameters::constrain_param_from_rootfile(std::string fileName,std::vector<std::string> names, int PDF, int bin){
|
|
return get_param_from_rootfile(fileName, names, PDF, bin, fixConstr(false,true));
|
|
}
|
|
|
|
int bu2kstarmumu_parameters::load_param_from_rootfile(std::string fileName,std::vector<std::string> names, int PDF, int bin){
|
|
return get_param_from_rootfile(fileName, names, PDF, bin, fixConstr(false, false));
|
|
}
|
|
|
|
//TODO: the following ones are not needed, as the string vector can be used instead....
|
|
void bu2kstarmumu_parameters::load_only_bckgnd_param_values(std::string filename){
|
|
this->load_param(filename, "cbkgmkpi0");
|
|
this->load_param(filename, "cbkgmkpi1");
|
|
this->load_param(filename, "cbkgmkpi2");
|
|
this->load_param(filename, "cbkgmkpi3");
|
|
this->load_param(filename, "cbkgmkpi4");
|
|
this->load_param(filename, "cbkgctl0");
|
|
this->load_param(filename, "cbkgctl1");
|
|
this->load_param(filename, "cbkgctl2");
|
|
this->load_param(filename, "cbkgctl3");
|
|
this->load_param(filename, "cbkgctl4");
|
|
this->load_param(filename, "cbkgctk0");
|
|
this->load_param(filename, "cbkgctk1");
|
|
this->load_param(filename, "cbkgctk2");
|
|
this->load_param(filename, "cbkgctk3");
|
|
this->load_param(filename, "cbkgctk4");
|
|
this->load_param(filename, "cbkgctk5");
|
|
this->load_param(filename, "cbkgctk6");
|
|
this->load_param(filename, "cbkgphi0");
|
|
this->load_param(filename, "cbkgphi1");
|
|
this->load_param(filename, "cbkgphi2");
|
|
this->load_param(filename, "cbkgphi3");
|
|
this->load_param(filename, "cbkgphi4");
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::load_only_Bmass_param_values(std::string filename){
|
|
this->load_param(filename, "sigma_1");
|
|
this->load_param(filename, "alpha_1");
|
|
this->load_param(filename, "n_1");
|
|
this->load_param(filename, "sigma_2");
|
|
this->load_param(filename, "alpha_2");
|
|
this->load_param(filename, "n_2");
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::save_param_values(std::string filename){
|
|
|
|
std::ostringstream sout;
|
|
spdlog::info("Parameter values are being saved to file "+ filename);;
|
|
|
|
//write all parameters to oss
|
|
for(UInt_t p = 0; p < this->nparameters(); p++){
|
|
sout << std::endl;
|
|
sout << p;
|
|
sout << "\t" << this->get_parameter(p)->get_name();
|
|
if (this->get_parameter(p)->get_error()==0) continue; //Print only not-fixed parameters
|
|
sout << "\t" << this->get_parameter(p)->get_value();
|
|
sout << "\t" << this->get_parameter(p)->get_error();
|
|
if (TMath::Abs(this->get_parameter(p)->get_error()/this->get_parameter(p)->get_value()) > 0.5){
|
|
sout << "\tLARGE ERROR"; //c++s abs returned an int for whatever reason
|
|
spdlog::warn("Par " + this->get_parameter(p)->get_name() + " has a large error:\t {0:f} +- {1:f}", this->get_parameter(p)->get_value() ,this->get_parameter(p)->get_error() );
|
|
}
|
|
else{
|
|
spdlog::info("\tPar " + this->get_parameter(p)->get_name() + ":\t {0:f} +- {1:f}", this->get_parameter(p)->get_value() ,this->get_parameter(p)->get_error() );
|
|
}
|
|
}
|
|
//save oss to file
|
|
std::ofstream file(filename.c_str());
|
|
file << sout.str();
|
|
file.close();
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::add_parameters(){
|
|
//physics parameters
|
|
add_parameter(&Fl);
|
|
add_parameter(&S1s);
|
|
add_parameter(&S3);
|
|
add_parameter(&S4);
|
|
add_parameter(&S5);
|
|
add_parameter(&Afb);
|
|
add_parameter(&S6s);
|
|
add_parameter(&S7);
|
|
add_parameter(&S8);
|
|
add_parameter(&S9);
|
|
add_parameter(&P1);
|
|
add_parameter(&P2);
|
|
add_parameter(&P3);
|
|
add_parameter(&P4);
|
|
add_parameter(&P5);
|
|
add_parameter(&P6);
|
|
add_parameter(&P8);
|
|
add_parameter(&FS);
|
|
add_parameter(&SS1);
|
|
add_parameter(&SS2);
|
|
add_parameter(&SS3);
|
|
add_parameter(&SS4);
|
|
add_parameter(&SS5);
|
|
|
|
//had this at the start
|
|
//signal fraction
|
|
add_parameter(&f_sig);
|
|
add_parameter(&n_sig);
|
|
add_parameter(&n_bkg);
|
|
//mass parameters
|
|
add_parameter(&m_b);
|
|
add_parameter(&m_res_1);
|
|
add_parameter(&m_sigma_1);
|
|
add_parameter(&m_sigma_2);
|
|
add_parameter(&m_scale);
|
|
add_parameter(&alpha_1);
|
|
add_parameter(&alpha_2);
|
|
add_parameter(&n_1);
|
|
add_parameter(&n_2);
|
|
|
|
add_parameter(&fm_tau);
|
|
add_parameter(&m_tau);
|
|
add_parameter(&m_tau_2);
|
|
add_parameter(&m_lambda);
|
|
add_parameter(&m_lambda_2);
|
|
add_parameter(&eff_q2);
|
|
|
|
add_parameter(&asphase);
|
|
add_parameter(&a);
|
|
add_parameter(&r);
|
|
add_parameter(&gammakstar);
|
|
add_parameter(&mkstar);
|
|
add_parameter(&gammakstarplus);
|
|
add_parameter(&mkstarplus);
|
|
|
|
add_parameter(&mf800);
|
|
add_parameter(&gammaf800);
|
|
add_parameter(&f800mag);
|
|
add_parameter(&f800phase);
|
|
|
|
add_parameter(&cbkgctl0);
|
|
add_parameter(&cbkgctl1);
|
|
add_parameter(&cbkgctl2);
|
|
add_parameter(&cbkgctl3);
|
|
add_parameter(&cbkgctl4);
|
|
|
|
add_parameter(&cbkgctk0);
|
|
add_parameter(&cbkgctk1);
|
|
add_parameter(&cbkgctk2);
|
|
add_parameter(&cbkgctk3);
|
|
add_parameter(&cbkgctk4);
|
|
add_parameter(&cbkgctk5);
|
|
add_parameter(&cbkgctk6);
|
|
|
|
add_parameter(&cbkgphi0);
|
|
add_parameter(&cbkgphi1);
|
|
add_parameter(&cbkgphi2);
|
|
add_parameter(&cbkgphi3);
|
|
add_parameter(&cbkgphi4);
|
|
|
|
add_parameter(&cbkgmkpi0);
|
|
add_parameter(&cbkgmkpi1);
|
|
add_parameter(&cbkgmkpi2);
|
|
add_parameter(&cbkgmkpi3);
|
|
add_parameter(&cbkgmkpi4);
|
|
|
|
add_parameter(&cswavemkpi0);
|
|
add_parameter(&cswavemkpi1);
|
|
add_parameter(&cswavemkpi2);
|
|
add_parameter(&cswavemkpi3);
|
|
add_parameter(&cswavemkpi4);
|
|
|
|
add_parameter(&cbkgp20);
|
|
add_parameter(&cbkgp21);
|
|
add_parameter(&cbkgp22);
|
|
add_parameter(&cbkgp23);
|
|
add_parameter(&cbkgp24);
|
|
|
|
add_parameter(&nthreshold);
|
|
add_parameter(&R);
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_mass_parameters(int PDF, int nBins, int bin, double defaultStepSize){
|
|
|
|
//Initialize the mass fit parameters based on the bin, as pi0 screws everyhting up
|
|
std::vector<std::vector<double>> sMassPar = init_mass_params_MC(nBins,PDF);
|
|
//Not the most efficient to always get the full vector with all bins, but oh well
|
|
if(opts->twotailedcrystalball){
|
|
m_res_1.init(1.0, 0.0, 1.0, 0.0);
|
|
m_sigma_1.init(sMassPar[0][bin],15.0, 50.0, defaultStepSize);
|
|
}
|
|
else{
|
|
m_res_1.init(0.5, 0.0, 1.0, defaultStepSize);
|
|
m_sigma_1.init(sMassPar[0][bin],15.0, 40.0, defaultStepSize);
|
|
m_sigma_2.init(sMassPar[1][bin],15.0, 40.0, defaultStepSize);
|
|
}
|
|
|
|
alpha_1.init(sMassPar[2][bin], 0.5, 3.0, defaultStepSize);
|
|
alpha_2.init(sMassPar[3][bin], 0.5, 3.0, defaultStepSize);
|
|
n_1.init(sMassPar[4][bin], 2.0, 10.0, defaultStepSize);
|
|
n_2.init(sMassPar[5][bin], 2.0, 10.0, defaultStepSize);
|
|
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_Bmass(std::string fileName, int PDF, double stepSize, fixConstr FC){
|
|
double massInit = get_param_value_from_rootfile(fileName,"m_b", PDF, 0);
|
|
//fixing the Bmass outweights constraining it:
|
|
if (FC.fix) m_b.init(massInit,B_MASS_LOW,B_MASS_HIGH, 0.0);
|
|
else if (FC.constrain) m_b.init(massInit,B_MASS_LOW,B_MASS_HIGH, stepSize, 1.0, true);
|
|
else m_b.init(massInit,B_MASS_LOW,B_MASS_HIGH, stepSize);
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_angular_background_parameters(bool fitReference, double stepsize){
|
|
//Order of the polynomial is 0-6, each parameter represents x^n
|
|
//TODO: CHECK THE FOLDINGS
|
|
|
|
//fitReference is not used, but in case this needs to be different for signal/reference!
|
|
//First set everything to be flat so one can only change the parameters that actually vary
|
|
use_default_bkg();
|
|
|
|
std::vector<double> init_values = init_bkg(opts->folding, opts->bkg_order_costhetal,opts->bkg_order_costhetak);
|
|
for_indexed(auto name: PAR_BKG_STRING(opts->folding, opts->bkg_order_costhetal,opts->bkg_order_costhetak)){
|
|
std::vector<double> range = GetParamRange(name);
|
|
this->get_parameter(name)->init(init_values[i],range.front(), range.back(),stepsize);
|
|
spdlog::trace("Initialized " + this->get_parameter(name)->get_name());
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_mass_background_parameters(int nBins, int bin, bool useLambda){
|
|
std::vector<double> f_lambda = init_f_lambda(nBins);
|
|
if(useLambda){
|
|
m_lambda.init(f_lambda.at(bin), PAR_LAMBDA*PAR_LAMBDA_SCALE, PAR_LAMBDA/PAR_LAMBDA_SCALE, 1.0e-5);
|
|
//Turn on/off second exponential possibly
|
|
m_lambda_2.init(f_lambda.at(bin), PAR_LAMBDA*PAR_LAMBDA_SCALE, PAR_LAMBDA/PAR_LAMBDA_SCALE, FIX_FM ? 0.0 : 1.0e-5);
|
|
m_tau.init_fixed(0.0);
|
|
m_tau_2.init_fixed(0.0);
|
|
}
|
|
else{
|
|
//Not really used, so if someone ever needs this, modify yourselves!!!
|
|
m_tau.init(-1.0/f_lambda.at(bin), 100, 100000, FIX_FM ? 0.0 : 10.0);
|
|
//Turn on/off second exponential possibly
|
|
m_tau_2.init(-1.0/f_lambda.at(bin), 100, 100000, FIX_FM ? 0.0 : 10.0);
|
|
m_lambda.init_fixed(0.0);
|
|
m_lambda_2.init_fixed(0.0);
|
|
}
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_mkpi_pWave_parameters(bool fitReference, double stepsize){
|
|
//fitReference is not used, but in case this needs to be different for signal/reference!
|
|
gammakstar.init(0.065, 0.055, 0.07, stepsize);
|
|
mkstar.init(PDGMASS_K_STAR_PLUS / 1000., 0.882, 0.902, 0.0); //PDG value
|
|
R.init(1.6, 0.0, 10.0, 0.00);
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_mkpi_sWave_parameters(bool fitReference, double stepsize){
|
|
//fitReference is not used, but in case this needs to be different for signal/reference!
|
|
gammakstarplus.init(0.236, 0.1, 0.5, stepsize); //PDG value is 0.236
|
|
mkstarplus.init(1.41,1.2,1.6, 0.00); //PDG value is 1.41
|
|
asphase.init(TMath::Pi(), 0.0, 2.0*TMath::Pi(), 0.0);
|
|
a.init(1.95, 0.0, 20.0, 0.0);
|
|
r.init(1.78, 0.0, 10.0, 0.0);
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_kpi_background_parameters(bool fitReference, double stepsize){
|
|
//fitReference is not used, but in case this needs to be different for signal/reference!
|
|
cbkgmkpi1.init(-0.99, -1.25, 0.25, stepsize);
|
|
//cbkgmkpi2.init(0.3, -1.0, 1.0, stepsize);
|
|
return;
|
|
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_angular_parameters(int nBins, int bin, double stepsize, double scale, bool blind){
|
|
//Scale sets the range
|
|
|
|
//First, make sure everything is at the starting point
|
|
use_default_observables();
|
|
|
|
//Initialize angular parameters from constants
|
|
//s1s //s3 //s4 //s5 //s6s //s7 //s8 //s9
|
|
std::vector<std::vector<double>> sParameters;
|
|
|
|
if (opts->initSM){
|
|
spdlog::debug("Loading angular parameters from flavio SM preditions");
|
|
sParameters = init_angular_params(nBins,false,false);
|
|
}
|
|
else{
|
|
spdlog::debug("Loading angular parameters from MC result file");
|
|
sParameters = init_angular_params_MC(nBins);
|
|
}
|
|
|
|
double blind_scale = 1.0; //Don't know what it does, but it is set to 1.0 everywhere
|
|
|
|
if (!opts->fit_pprimes){
|
|
|
|
if(opts->fit_fl) Fl.init(1.0-4.0/3.0*sParameters.at(0).at(bin), 0.0, 2.0*scale, stepsize);
|
|
else S1s.init(sParameters.at(0).at(bin), -scale, scale, stepsize);
|
|
|
|
S3.init(sParameters.at(1).at(bin), -scale, scale, stepsize);
|
|
S4.init(sParameters.at(2).at(bin), -scale, scale, (opts->full_angular || opts->folding == 1 ? stepsize : 0.0));
|
|
S5.init(sParameters.at(3).at(bin), -scale, scale, (opts->full_angular || opts->folding == 2 ? stepsize : 0.0));
|
|
|
|
if(opts->fit_afb) Afb.init(3.0/4.0*sParameters.at(4).at(bin), -0.75*scale, 0.75*scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
else S6s.init(sParameters.at(4).at(bin), -scale, scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
|
|
S7.init(sParameters.at(5).at(bin), -scale, scale, (opts->full_angular || opts->folding == 3 ? stepsize : 0.0));
|
|
S8.init(sParameters.at(6).at(bin), -scale, scale, (opts->full_angular || opts->folding == 4 ? stepsize : 0.0));
|
|
S9.init(sParameters.at(7).at(bin), -scale, scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
if(blind){
|
|
Fl.set_blinding(true, blind_scale, true);
|
|
P1.set_blinding(true, blind_scale, true);
|
|
P2.set_blinding(true, blind_scale, true);
|
|
P3.set_blinding(true, blind_scale, true);
|
|
P4.set_blinding(true, blind_scale, true);
|
|
P5.set_blinding(true, blind_scale, true);
|
|
P6.set_blinding(true, blind_scale, true);
|
|
P8.set_blinding(true, blind_scale, true);
|
|
}
|
|
}
|
|
else{
|
|
double fl = 1.0-4.0/3.0*sParameters.at(0).at(bin);
|
|
Fl.init(fl, 0.0, 2.0*scale, stepsize);
|
|
P1.init(2*sParameters.at(0).at(bin)/(1.-fl) , -scale, scale, stepsize);
|
|
P2.init(sParameters.at(4).at(bin)/2/(1.-fl) , -scale, scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
P3.init(-sParameters.at(1).at(bin)/(1.-fl) , -scale, scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
P4.init(sParameters.at(2).at(bin)/sqrt(fl*(1.-fl)), -scale, scale, (opts->full_angular || opts->folding == 1 ? stepsize : 0.0));
|
|
P5.init(sParameters.at(3).at(bin)/sqrt(fl*(1.-fl)), -scale, scale, (opts->full_angular || opts->folding == 2 ? stepsize : 0.0));
|
|
P6.init(sParameters.at(5).at(bin)/sqrt(fl*(1.-fl)), -scale, scale, (opts->full_angular || opts->folding == 3 ? stepsize : 0.0));
|
|
P8.init(sParameters.at(6).at(bin)/sqrt(fl*(1.-fl)), -scale, scale, (opts->full_angular || opts->folding == 4 ? stepsize : 0.0));
|
|
|
|
if(blind){
|
|
if(opts->fit_fl) Fl.set_blinding(true, blind_scale, true);
|
|
else S1s.set_blinding(true, blind_scale, true);
|
|
S3.set_blinding(true, blind_scale, true);
|
|
S4.set_blinding(true, blind_scale, true);
|
|
S5.set_blinding(true, blind_scale, true);
|
|
if(opts->fit_afb) Afb.set_blinding(true, blind_scale, true);
|
|
else S6s.set_blinding(true, blind_scale, true);
|
|
S7.set_blinding(true, blind_scale, true);
|
|
S8.set_blinding(true, blind_scale, true);
|
|
S9.set_blinding(true, blind_scale, true);
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_ang_parameters_fromRefDavid(double stepsize, double scale, bool blind){
|
|
|
|
std::vector<double> angPars = init_angular_params_RefFromDavid();
|
|
//Not the most efficient to always get the full vector, but oh well
|
|
|
|
if (opts->fit_pprimes){
|
|
spdlog::critical("init_ang_parameters_fromRefDavid not implemented yet for Pprimes!");
|
|
assert(0);
|
|
}
|
|
|
|
if(opts->fit_fl) Fl.init(1.0-4.0/3.0*angPars.at(0), 0.0, 2.0*scale, stepsize);
|
|
else S1s.init(angPars.at(0), -scale, scale, stepsize);
|
|
|
|
S3.init(angPars.at(1), -scale, scale, stepsize);
|
|
S4.init(angPars.at(2), -scale, scale, (opts->full_angular || opts->folding == 1 ? stepsize : 0.0));
|
|
S5.init(angPars.at(3), -scale, scale, (opts->full_angular || opts->folding == 2 ? stepsize : 0.0));
|
|
|
|
if(opts->fit_afb) Afb.init(3.0/4.0*angPars.at(4), -0.75*scale, 0.75*scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
else S6s.init(angPars.at(4), -scale, scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
|
|
S7.init(angPars.at(5), -scale, scale, (opts->full_angular || opts->folding == 3 ? stepsize : 0.0));
|
|
S8.init(angPars.at(6), -scale, scale, (opts->full_angular || opts->folding == 4 ? stepsize : 0.0));
|
|
S9.init(angPars.at(7), -scale, scale, (opts->full_angular || opts->folding == 0 ? stepsize : 0.0));
|
|
|
|
int blind_scale = 1.0; //No clue what this is good for, but it is 1.0 everywhere, so
|
|
if(blind){
|
|
if(opts->fit_fl) Fl.set_blinding(true, blind_scale, true);
|
|
else S1s.set_blinding(true, blind_scale, true);
|
|
S3.set_blinding(true, blind_scale, true);
|
|
S4.set_blinding(true, blind_scale, true);
|
|
S5.set_blinding(true, blind_scale, true);
|
|
if(opts->fit_afb) Afb.set_blinding(true, blind_scale, true);
|
|
else S6s.set_blinding(true, blind_scale, true);
|
|
S7.set_blinding(true, blind_scale, true);
|
|
S8.set_blinding(true, blind_scale, true);
|
|
S9.set_blinding(true, blind_scale, true);
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
void bu2kstarmumu_parameters::init_sWave_parameters(double stepsize){
|
|
//TODO: very basic for now
|
|
double scale = 1.0;
|
|
SS1.init( -0.153, -scale, scale, opts->full_angular || opts->folding != 4 ? stepsize : 0.0);
|
|
SS2.init( 0.032, -scale, scale, opts->full_angular || opts->folding == 1 ? stepsize : 0.0);
|
|
SS3.init( -0.001, -scale, scale, opts->full_angular || opts->folding == 2 ? stepsize : 0.0);
|
|
SS4.init( 0.001, -scale, scale, opts->full_angular || opts->folding > 2 ? stepsize : 0.0);
|
|
SS5.init( -0.070, -scale, scale, opts->full_angular ? stepsize : 0.0);
|
|
return;
|
|
}
|
|
|
|
bu2kstarmumu_parameters::bu2kstarmumu_parameters(options* o):
|
|
opts(o)
|
|
{
|
|
use_default();
|
|
add_parameters();
|
|
}
|
|
|
|
|
|
double eventsInBin_fraction(int bin, int run, int nBins, bool fromRef){
|
|
//From the requested total number of events, get the number of events in the given bin for given pdf
|
|
//Get the number of events from the dataFile and from the fraction of nSig/(nSig+nBkg)
|
|
|
|
assert(bin >= 0 && bin < nBins);
|
|
assert(run == 1 || run == 2 || run == 12);
|
|
|
|
TChain * chain = new TChain("Events");
|
|
spdlog::trace("Run 1 data file ='" + get_theFCNCpath(0,1) + "'");
|
|
chain->Add(get_theFCNCpath(0,1).c_str());
|
|
spdlog::trace("Run 2 data file ='" + get_theFCNCpath(0,2) + "'");
|
|
chain->Add(get_theFCNCpath(0,2).c_str());
|
|
spdlog::debug("Chain entries: {0:d}", chain->GetEntries());
|
|
|
|
//cut on the q2
|
|
std::vector<double> q2BinsMin = get_TheQ2binsmin(nBins, fromRef); //No need to have specifically 1 for Ref
|
|
std::vector<double> q2BinsMax = get_TheQ2binsmax(nBins, fromRef); //as if reference, always returns one bin
|
|
std::string whichYear = "";
|
|
switch(run){
|
|
case 1: whichYear = "(year<2013) &&"; break;
|
|
case 2: whichYear = "(year>2013) &&"; break;
|
|
case 12: whichYear = "(year<2020) &&"; break;
|
|
}
|
|
std::string cutOnQ2 = " ";
|
|
for (unsigned int b = 0; b < q2BinsMin.size(); b++){
|
|
cutOnQ2.append( "(" + std::to_string(q2BinsMin[b]) + " < q2 && q2 < "+ std::to_string(q2BinsMax[b]) +") || ");
|
|
}
|
|
cutOnQ2.erase(cutOnQ2.length()-3);
|
|
spdlog::debug("Using cut: '" + cutOnQ2 + "'");
|
|
int n_all = chain->GetEntries(cutOnQ2.c_str()); //number of events in the rare q2 region
|
|
if(n_all == -1) spdlog::debug("Getting number of all entries from TChain failed. Using hardcoded numbers.");
|
|
spdlog::debug("All entries: {0:d}", n_all);
|
|
|
|
cutOnQ2 = whichYear + std::to_string(q2BinsMin[bin]) + " < q2 && q2 < "+ std::to_string(q2BinsMax[bin]);
|
|
spdlog::debug("Using cut: '" + cutOnQ2 + "'");
|
|
int n_bin = chain->GetEntries(cutOnQ2.c_str()); //number of events in the desired q2 bin
|
|
if(n_bin == -1)spdlog::debug("Getting number of entries for bin {0:d} from TChain failed. Using hardcoded numbers.", bin);
|
|
spdlog::debug("Entries in bin: {0:d}", n_bin);
|
|
|
|
//Running on Condor does not allow usage of TTreePlayer, which is needed for the GetEntries with selection cut
|
|
//In that case, -1 is returned. If so, use hardcoded fractions for the 4 or 5 q2 binning scheme.
|
|
if(n_all == -1 || n_bin == -1){
|
|
spdlog::warn("Could not determine event fraction in bin={0:d} for run={1:d} from TChain. Use hardcoded numbers.", bin, run);
|
|
if(nBins == 4 || nBins == 5){
|
|
std::vector<std::vector<double> > frac = fracs(nBins);
|
|
if(run == 12)return frac.at(0).at(bin)+frac.at(1).at(bin);
|
|
else return frac.at(run-1).at(bin);
|
|
}
|
|
spdlog::error("No hardcoded numbers available for {0:d} q2 bins.", nBins);
|
|
return 0.0;
|
|
}
|
|
double fraction = double(n_bin)/double(n_all);
|
|
spdlog::debug("Returning fraction of {0:f}", fraction);
|
|
chain->Clear();
|
|
return fraction;
|
|
}
|
|
|
|
void EventNumbers(unsigned int bin, unsigned int pdf, double & fraction, unsigned int & eventnumbers,
|
|
unsigned int TotalEvents, unsigned int bins, unsigned int pdfs){
|
|
|
|
bool use_unblinded_values = false;//TODO
|
|
|
|
//return the fraction of nSig/(nSig+nBkg) *** in the requested q2bin ***
|
|
//return the number of events (nSig+nBkg) *** in this q2bin *** and in this sub-set given a total number of events (all sub-sets and all q2bins)
|
|
fraction = 0;
|
|
eventnumbers = 0;
|
|
|
|
std::vector<double> f_signal = init_n_signal(bins);
|
|
std::vector<double> f_bckgnd = init_n_bckgnd(bins);
|
|
std::vector<double> f_subset = get_f_subset(pdfs);
|
|
|
|
assert(f_signal.size() == bins);
|
|
assert(f_bckgnd.size() == bins);
|
|
assert(f_subset.size() == pdfs);
|
|
|
|
//assign signal fraction
|
|
fraction = f_signal.at(bin) / (f_signal.at(bin) + f_bckgnd.at(bin));
|
|
|
|
//get total event number in this bin and for this sub-set:
|
|
double f_sig_norm = sum_vector(f_signal) + sum_vector(f_bckgnd);
|
|
double f_event_norm = sum_vector(f_subset);
|
|
|
|
for(auto sig: f_signal) sig /= f_sig_norm;
|
|
for(auto bkg: f_bckgnd) bkg /= f_sig_norm;
|
|
for(auto f_sub: f_subset) f_sub /= f_event_norm;
|
|
|
|
eventnumbers = (unsigned int) (0.5 + TotalEvents * (f_signal.at(bin) + f_bckgnd.at(bin)) * f_subset.at(pdf));
|
|
|
|
if(use_unblinded_values){ //Number of signal and background events for each PDF //TODO: de-hardcode this one
|
|
fraction = sig_unblinded.at(bin).at(pdf) / (sig_unblinded.at(bin).at(pdf) + bkg_unblinded.at(bin).at(pdf));
|
|
eventnumbers = (unsigned int) (0.5 + sig_unblinded.at(bin).at(pdf) + bkg_unblinded.at(bin).at(pdf));
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
|