1218 lines
44 KiB
C++
1218 lines
44 KiB
C++
//Renata Kopecna
|
|
|
|
#include <helpers.hh>
|
|
#include <fstream> //Needed to remove/copy files
|
|
#include <bu2kstarmumu_parameters.hh>
|
|
|
|
#include <sstream>
|
|
#include <iostream>
|
|
|
|
#include <options.hh>
|
|
#include <parameters.hh>
|
|
#include <constants.hh>
|
|
#include <spdlog.h>
|
|
#include <fmt/ostr.h>
|
|
#include <TFile.h>
|
|
#include <TTree.h>
|
|
|
|
#include <paths.hh>
|
|
#include <numeric> //includes accumulate
|
|
#include <parse.hh>
|
|
|
|
#include <event.hh>
|
|
#include <spdlog.h>
|
|
|
|
//----------------------//
|
|
// Print utils //
|
|
//----------------------//
|
|
std::string boolToString(bool isTrue){
|
|
if (isTrue) return "TRUE";
|
|
else return "FALSE";
|
|
}
|
|
|
|
void reset_spdlog(){
|
|
spdlog::set_pattern("[%^%l%$]\t %v");
|
|
}
|
|
|
|
void set_spdlog_level(int verboseLvl){
|
|
switch(verboseLvl) {
|
|
case 0:
|
|
spdlog::set_level(spdlog::level::trace);
|
|
break;
|
|
case 1:
|
|
spdlog::set_level(spdlog::level::debug);
|
|
break;
|
|
case 2:
|
|
spdlog::set_level(spdlog::level::info);
|
|
break;
|
|
case 3:
|
|
spdlog::set_level(spdlog::level::warn);
|
|
break;
|
|
case 4:
|
|
spdlog::set_level(spdlog::level::err);
|
|
break;
|
|
case 5:
|
|
spdlog::set_level(spdlog::level::critical);
|
|
break;
|
|
}
|
|
return;
|
|
|
|
}
|
|
|
|
bool spdlog_trace(){
|
|
return (spdlog::default_logger_raw()->level()<=spdlog::level::trace);
|
|
}
|
|
bool spdlog_debug(){
|
|
return (spdlog::default_logger_raw()->level()<=spdlog::level::debug);
|
|
}
|
|
bool spdlog_info(){
|
|
return (spdlog::default_logger_raw()->level()<=spdlog::level::info);
|
|
}
|
|
bool spdlog_warn(){
|
|
return (spdlog::default_logger_raw()->level()<=spdlog::level::warn);
|
|
}
|
|
bool spdlog_error(){
|
|
return (spdlog::default_logger_raw()->level()<=spdlog::level::err);
|
|
}
|
|
|
|
std::string get_sample_from_jobID(int job_id){
|
|
if(job_id == 0) return "DATA";
|
|
else if(job_id == 1) return "SIGNAL MC";
|
|
else if(job_id == 2) return "REFERENCE MC ";
|
|
else if(job_id == 3) return "PHSP MC ";
|
|
else if(job_id == 4) return "GenLevel PHSP MC";
|
|
else if(job_id == 5) return "GenLevel MC";
|
|
else {
|
|
spdlog::error("Invalid index for conversion function given: {0:d}", job_id);
|
|
assert(0);
|
|
}
|
|
}
|
|
|
|
std::string format_double(double value, unsigned int digits){
|
|
std::stringstream out;
|
|
out << std::fixed << std::setprecision(digits) << value;
|
|
return out.str();
|
|
};
|
|
|
|
std::string format_value(double value, double error, unsigned int digits){
|
|
//format value such that the relevant number of digits is given back
|
|
//determine first "digits" significant places of error
|
|
if (error == 0.0){
|
|
std::stringstream out; //to_string(value) returns something else than this.... I hate c++
|
|
out << value;
|
|
return out.str();
|
|
}
|
|
else{
|
|
int first_significant = floor(log10(error));
|
|
if (first_significant > 0) first_significant = 0;
|
|
//ln(e^x)=x log10(c*10^x)=x+log10(c) with log10(c) < 1
|
|
return format_double(value,fabs(first_significant)+digits-1);
|
|
}
|
|
};
|
|
|
|
std::string format_error(double error, unsigned int digits){
|
|
return (format_value(error,error,digits));
|
|
};
|
|
|
|
|
|
//----------------------//
|
|
// Year utils //
|
|
//----------------------//
|
|
//generate vector with years for every run option
|
|
//notation for usage of runs:
|
|
//1 : Run 1
|
|
//2 : Run 2
|
|
//12: Run 1+2
|
|
//21: 'Run 2.1' = 2015+2016
|
|
//22: 'Run 2.2' = 2017+2018
|
|
|
|
std::vector<int> get_years(int Run, bool MCsig, bool MCref){
|
|
std::vector<int>years;
|
|
|
|
if(Run == 1 || Run == 12){
|
|
years.push_back(2011);
|
|
years.push_back(2012);
|
|
}
|
|
|
|
if(Run== 2 || Run == 12 || Run== 21){
|
|
if (!MCsig) years.push_back(2015); //No MC for 2015
|
|
years.push_back(2016);
|
|
}
|
|
if (!MCref){ //No RefMC for 2017+2018
|
|
if(Run == 2 || Run == 12 || Run == 22){
|
|
years.push_back(2017);
|
|
years.push_back(2018);
|
|
}
|
|
}
|
|
return years;
|
|
}
|
|
|
|
//Check if year is a run year
|
|
void check_year(int year){
|
|
std::set<int> yearSet = {2011,2012,2015,2016,2017,2018};
|
|
//Check if year is in a set, if not, throw exception
|
|
assert(yearSet.find(year)!= yearSet.end());
|
|
return;
|
|
}
|
|
|
|
int MC_dataset(bool Reference, bool PHSP){
|
|
if (PHSP) return 3; //PHSP
|
|
else if (Reference) return 2; //Ref MC
|
|
else return 1; //Sig MC
|
|
}
|
|
|
|
int get_yearFromRun(int year){
|
|
if (year == 2011 || year == 2012) return 1;
|
|
if (year == 2015 || year == 2016 || year == 2017 || year == 2018) return 2;
|
|
spdlog::error("Year not recognized!");
|
|
assert(0);
|
|
}
|
|
|
|
//--------------------------------------//
|
|
// General helpers //
|
|
//--------------------------------------//
|
|
|
|
//Copy file
|
|
int copyFile(std::string from, std::string to){
|
|
std::ifstream src(from);
|
|
std::ofstream dst(to);
|
|
dst << src.rdbuf();
|
|
return 0;
|
|
}
|
|
|
|
//Check if file exists
|
|
//https://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c
|
|
bool existsTest (const std::string& name) {
|
|
struct stat buffer;
|
|
return (stat (name.c_str(), &buffer) == 0);
|
|
}
|
|
|
|
//https://stackoverflow.com/questions/50960492/creating-folders-in-c
|
|
//Create parent folders if needed
|
|
void makeFolder (const std::string& name) {
|
|
std::string parentfolder = name.substr(0, name.find_last_of("/\\"));
|
|
if(!existsTest(parentfolder))makeFolder(parentfolder);
|
|
spdlog::info("Creating folder '"+name+"'");
|
|
mkdir(name.c_str(), 0755);
|
|
return;
|
|
}
|
|
|
|
//Convert vector into number-space-number... string
|
|
std::string convert_vector_to_string(std::vector<int> myVector){
|
|
std::stringstream ss;
|
|
copy(myVector.begin(), myVector.end(), std::ostream_iterator<int>(ss, " "));
|
|
std::string s = ss.str();
|
|
return s.substr(0, s.length()-1); // get rid of the trailing space
|
|
}
|
|
|
|
std::string convert_vector_to_string(std::vector<double> myVector){
|
|
std::stringstream ss;
|
|
copy(myVector.begin(), myVector.end(), std::ostream_iterator<double>(ss, " "));
|
|
std::string s = ss.str();
|
|
return s.substr(0, s.length()-1); // get rid of the trailing space
|
|
}
|
|
std::string convert_vector_to_string(std::vector<std::string> myVector){
|
|
std::stringstream ss;
|
|
copy(myVector.begin(), myVector.end(), std::ostream_iterator<std::string>(ss, " "));
|
|
std::string s = ss.str();
|
|
return s.substr(0, s.length()-1); // get rid of the trailing space
|
|
}
|
|
|
|
std::vector<double> merge_two_vecs(std::vector<double> one, std::vector<double> two){
|
|
std::vector<double> merged_vec;
|
|
merged_vec.insert(merged_vec.end(), one.begin(), one.end());
|
|
merged_vec.insert(merged_vec.end(), two.begin(), two.end());
|
|
return merged_vec;
|
|
}
|
|
std::vector<int> merge_two_vecs(std::vector<int> one, std::vector<int> two){
|
|
std::vector<int> merged_vec;
|
|
merged_vec.insert(merged_vec.end(), one.begin(), one.end());
|
|
merged_vec.insert(merged_vec.end(), two.begin(), two.end());
|
|
return merged_vec;
|
|
}
|
|
std::vector<std::string> merge_two_vecs(std::vector<std::string> one, std::vector<std::string> two){ //TODO: I bet there is a way to avoid overloading
|
|
std::vector<std::string> merged_vec;
|
|
merged_vec.insert(merged_vec.end(), one.begin(), one.end());
|
|
merged_vec.insert(merged_vec.end(), two.begin(), two.end());
|
|
return merged_vec;
|
|
}
|
|
|
|
bool replace(std::string& str, const std::string& from, const std::string& to) {
|
|
size_t start_pos = str.find(from);
|
|
if(start_pos == std::string::npos)
|
|
return false;
|
|
while ( (start_pos = str.find(from,start_pos)) != std::string::npos){
|
|
str.replace(start_pos, from.length(), to);
|
|
start_pos += to.length(); // Handles case where 'to' is a substring of 'from'
|
|
}
|
|
return true;
|
|
}
|
|
|
|
//Sum a vector of ints
|
|
int sum_vector(std::vector<int> vec){
|
|
return std::accumulate(vec.begin(), vec.end(), 0);
|
|
}
|
|
//Sum a vector of doubles
|
|
double sum_vector(std::vector<double> vec){
|
|
return std::accumulate(vec.begin(), vec.end(), 0.0);
|
|
}
|
|
|
|
//Return the center of a bin
|
|
double bin_center(double low, double high){
|
|
return (high+low)/2.0;
|
|
}
|
|
|
|
//Return the center of a q2 bin
|
|
double bin_center_q2(fcnc::options opts, int b){
|
|
return bin_center(opts.TheQ2binsmin.at(b),opts.TheQ2binsmax.at(b));
|
|
}
|
|
double bin_center_q2(fcnc::options *opts, int b){
|
|
return bin_center(opts->TheQ2binsmin.at(b),opts->TheQ2binsmax.at(b));
|
|
}
|
|
double bin_center_q2(std::vector<double> q2min, std::vector<double> q2max , int b){
|
|
return bin_center(q2min.at(b),q2max.at(b));
|
|
}
|
|
|
|
//Return half of the width (useful for errors)
|
|
double bin_halfWidth(double low, double high){
|
|
return (high-low)/2.0;
|
|
}
|
|
double bin_halfWidth_q2(std::vector<double> q2min, std::vector<double> q2max, int b){
|
|
return (q2max.at(b)-q2min.at(b))/2.0;
|
|
}
|
|
double bin_halfWidth_q2(fcnc::options opts, int b){
|
|
return (opts.TheQ2binsmax.at(b)-opts.TheQ2binsmin.at(b))/2.0;
|
|
}
|
|
double bin_halfWidth_q2(fcnc::options *opts, int b){
|
|
return (opts->TheQ2binsmax.at(b)-opts->TheQ2binsmin.at(b))/2.0;
|
|
}
|
|
|
|
|
|
//Check if the angles are in the correct range given by options
|
|
bool isEvtInAngleRange(fcnc::event *evt, fcnc::options *opts){
|
|
if (evt->costhetal > opts->ctl_max) return false;
|
|
if (evt->costhetal < opts->ctl_min) return false;
|
|
if (evt->costhetak > opts->ctk_max) return false;
|
|
if (evt->costhetak < opts->ctk_min) return false;
|
|
if (evt->phi > opts->phi_max) return false;
|
|
if (evt->phi < opts->phi_min) return false;
|
|
if (evt->q2 > opts->TheQ2binsmax.back()) return false;
|
|
if (evt->q2 < opts->TheQ2binsmin.front()) return false;
|
|
return true;
|
|
|
|
}
|
|
//Check if the angles are in the correct range given by default
|
|
bool isEvtInAngleRange(fcnc::event *evt){
|
|
if (evt->costhetal > CTL_MAX){
|
|
spdlog::trace("ctl of event:\t{0:f}",evt->costhetal);
|
|
return false;
|
|
}
|
|
if (evt->costhetal < CTL_MIN){
|
|
spdlog::trace("ctl of event:\t{0:f}",evt->costhetal);
|
|
return false;
|
|
}
|
|
if (evt->costhetak > CTK_MAX) {
|
|
spdlog::trace("ctk of event:\t{0:f}",evt->costhetak);
|
|
return false;
|
|
}
|
|
if (evt->costhetak < CTK_MIN) {
|
|
spdlog::trace("ctk of event:\t{0:f}",evt->costhetak);
|
|
return false;
|
|
}
|
|
if (evt->phi > PHI_MAX){
|
|
spdlog::trace("phi of event:\t{0:f}",evt->phi);
|
|
return false;
|
|
}
|
|
if (evt->phi < PHI_MIN) {
|
|
spdlog::trace("phi of event:\t{0:f}",evt->phi);
|
|
return false;
|
|
}
|
|
if (evt->q2 < Q2_MIN_RANGE) return false;
|
|
if (evt->q2 > Q2_MAX_RANGE) return false;
|
|
return true;
|
|
|
|
}
|
|
|
|
//For now filter out events that are dumb dumb for folding four
|
|
bool filterFldFour(fcnc::event *evt, fcnc::options *opts){
|
|
if (opts->folding==4){
|
|
//TODO: this is idiotic and will cause trouble if one cuts on ctl. Unfortunatelly, otherwise it cuts the folded parameters too, which is not very beneficial
|
|
// SO this if (!isEvtInAngleRange(&evt,&opts)) continue; is not possible
|
|
if (evt->costhetak > opts->ctk_max) return false;
|
|
if (evt->costhetak < opts->ctk_min) return false;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
//Check if value is in a vector
|
|
bool isInVec(int key, std::vector<int>vec){
|
|
return (std::find(vec.begin(), vec.end(), key)!= vec.end());
|
|
}
|
|
bool isInVec(double key, std::vector<double>vec){
|
|
return (std::find(vec.begin(), vec.end(), key)!= vec.end());
|
|
}
|
|
bool isInVec(std::string key, std::vector<std::string>vec){
|
|
return (std::find(vec.begin(), vec.end(), key)!= vec.end());
|
|
}
|
|
|
|
|
|
//--------------------------------------//
|
|
// Helpers for converting the tuples //
|
|
//--------------------------------------//
|
|
|
|
//Returns the correct tree name as Data/MC and genLvl have different tree names
|
|
std::string get_inputTree_name(int job_id){
|
|
if(job_id == 0) return "DecayTree";
|
|
else if(job_id == 1) return "DecayTreeTruthMatched";
|
|
else if(job_id == 2) return "DecayTreeTruthMatched";
|
|
else if(job_id == 3) return "DecayTreeTruthMatched";
|
|
else if(job_id == 4) return "DecayTree";
|
|
else if(job_id == 5) return "DecayTree";
|
|
else {
|
|
spdlog::error("Invalid index for conversion function given: {0:d}", job_id);
|
|
assert(0);
|
|
}
|
|
}
|
|
|
|
//--------------------------------------//
|
|
// Helpers for angular corrections //
|
|
//--------------------------------------//
|
|
|
|
TF1 *fitStraightLine(TH1D* hist, std::string lineName, double lowEdge, double highEdge){
|
|
//Define the fit line
|
|
TF1 *fitLine= new TF1(lineName.c_str(), "[0]*x+[1]", lowEdge, highEdge);
|
|
|
|
fitLine->SetParameters(0.0, 1.0);
|
|
fitLine->SetParLimits(0,-0.5,0.5);
|
|
fitLine->SetParLimits(1,0.5,1.5);
|
|
fitLine->SetParNames("a","b");
|
|
|
|
std::string fitOption = spdlog_debug() ? "R" : "RQ"; //C turns off chi2 calculation, faster, Q=quiet
|
|
|
|
hist->Fit(lineName.c_str(), fitOption.c_str());
|
|
|
|
spdlog::debug("Fitted value is:\t{0:f}+-{1:f}",fitLine->GetParameter(0), fitLine->GetParError(0));
|
|
|
|
return fitLine;
|
|
}
|
|
|
|
//--------------------------------------//
|
|
// Helpers for bu2kstarmumu_pdf //
|
|
//--------------------------------------//
|
|
|
|
const std::vector<std::vector<std::string>> get_angObser_withTeX_vec(){
|
|
return {
|
|
{"1s","#xi_{1}^{s}"},
|
|
{"1c","#xi_{1}^{c}"},
|
|
{"2s","#xi_{2}^{s}"},
|
|
{"2c","#xi_{2}^{c}"},
|
|
{"3","#xi_{3}"},
|
|
{"4", "#xi_{4}"},
|
|
{"5","#xi_{5}"},
|
|
{"6s","#xi_{6}^{s}"},
|
|
{"6c","#xi_{6}^{c}"},
|
|
{"7","#xi_{7}"},
|
|
{"8","#xi_{8}"},
|
|
{"9","#xi_{9}"}
|
|
};
|
|
}
|
|
|
|
std::vector<std::string> get_angObser_vec(){
|
|
std::vector<std::string> tmp;
|
|
for (auto obs: get_angObser_withTeX_vec()){
|
|
tmp.push_back(obs[0]);
|
|
}
|
|
return tmp;
|
|
}
|
|
|
|
std::vector<std::string> get_angObserTeX_vec(){
|
|
std::vector<std::string> tmp;
|
|
for (auto obs: get_angObser_withTeX_vec()){
|
|
tmp.push_back(obs[1]);
|
|
}
|
|
return tmp;
|
|
}
|
|
|
|
//--------------------------------------//
|
|
// Helpers for time measuring //
|
|
//--------------------------------------//
|
|
|
|
void runTime::start(){
|
|
sw->Reset();
|
|
sw->Start();
|
|
return;
|
|
}
|
|
|
|
void runTime::stop(time_t startTime){
|
|
sw->Stop();
|
|
real_times.push_back(sw->RealTime());
|
|
cpu_times.push_back(sw->CpuTime());
|
|
cpp_times.push_back(Double_t(time(0) - startTime));
|
|
return;
|
|
}
|
|
|
|
void runTime::print(unsigned int nBins){
|
|
spdlog::info("");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| TIME OF FIT ||");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| BIN\t\tCPU time[s]\treal time[s]\tC++ time[s]\t||");
|
|
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
Int_t CpuHours = Int_t(cpu_times.at(b) / 3600);
|
|
cpu_times.at(b) -= CpuHours * 3600;
|
|
Int_t CpuMinutes = Int_t(cpu_times.at(b) / 60);
|
|
cpu_times.at(b) -= CpuMinutes * 60;
|
|
Int_t CpuSeconds = Int_t(cpu_times.at(b));
|
|
|
|
Int_t RealHours = Int_t(real_times.at(b) / 3600);
|
|
real_times.at(b) -= RealHours * 3600;
|
|
Int_t RealMinutes = Int_t(real_times.at(b) / 60);
|
|
real_times.at(b) -= RealMinutes * 60;
|
|
Int_t RealSeconds = Int_t(real_times.at(b));
|
|
|
|
Int_t CppHours = Int_t(cpp_times.at(b) / 3600);
|
|
cpp_times.at(b) -= CppHours * 3600;
|
|
Int_t CppMinutes = Int_t(cpp_times.at(b) / 60);
|
|
cpp_times.at(b) -= CppMinutes * 60;
|
|
Int_t CppSeconds = Int_t(cpp_times.at(b));
|
|
|
|
spdlog::info("|| {0:d}\t\t{1:d}:{2:02d}:{3:02d}\t\t{4:d}:{5:02d}:{6:02d}\t\t{7:d}:{8:02d}:{9:02d} \t||", b, CpuHours, CpuMinutes, CpuSeconds, RealHours, RealMinutes, RealSeconds, CppHours, CppMinutes, CppSeconds);
|
|
}
|
|
|
|
}
|
|
|
|
|
|
//--------------------------------------//
|
|
// Helpers for fit scripts //
|
|
//--------------------------------------//
|
|
std::vector<std::string> param_string_pPrimes(bool fitFL){
|
|
std::vector<std::string> params;
|
|
if(fitFL) params.push_back("Fl");
|
|
else params.push_back("S1s");
|
|
params.push_back("P1");
|
|
params.push_back("P4");
|
|
params.push_back("P5");
|
|
params.push_back("P2");
|
|
params.push_back("P6");
|
|
params.push_back("P8");
|
|
params.push_back("P3");
|
|
return params;
|
|
}
|
|
std::vector<std::string> param_string_pPrimes(bool fitFL, int folding){
|
|
std::vector<std::string> params;
|
|
if(fitFL) params.push_back("Fl");
|
|
else params.push_back("S1s");
|
|
if(folding == 0){
|
|
params.push_back("P2");
|
|
params.push_back("P3");
|
|
return params;
|
|
}
|
|
else if(folding == 1)params.push_back("P4"); //could be a switch, but meh
|
|
else if(folding == 2)params.push_back("P5");
|
|
else if(folding == 3)params.push_back("P6");
|
|
else if(folding == 4)params.push_back("P8");
|
|
return params;
|
|
}
|
|
|
|
std::vector<std::string> param_string_s(bool fitFL, bool fitAFB){
|
|
std::vector<std::string> params;
|
|
if(fitFL) params.push_back("Fl");
|
|
else params.push_back("S1s");
|
|
params.push_back("S3");
|
|
params.push_back("S4");
|
|
params.push_back("S5");
|
|
if (fitAFB) params.push_back("Afb");
|
|
else params.push_back("S6s");
|
|
params.push_back("S7");
|
|
params.push_back("S8");
|
|
params.push_back("S9");
|
|
return params;
|
|
}
|
|
std::vector<std::string> param_string_s(bool fitFL, bool fitAFB, int folding){
|
|
std::vector<std::string> params;
|
|
if(fitFL) params.push_back("Fl");
|
|
else params.push_back("S1s");
|
|
params.push_back("S3");
|
|
if(folding == 0){
|
|
if(fitAFB) params.push_back("Afb");
|
|
else params.push_back("S6s");
|
|
|
|
params.push_back("S9");
|
|
return params;
|
|
}
|
|
else if(folding == 1)params.push_back("S4");
|
|
else if(folding == 2)params.push_back("S5");
|
|
else if(folding == 3)params.push_back("S7");
|
|
else if(folding == 4)params.push_back("S8");
|
|
return params;
|
|
}
|
|
|
|
std::vector<std::string> param_string_sWave(){
|
|
std::vector<std::string> params;
|
|
params.push_back("FS");
|
|
params.push_back("SS1");
|
|
params.push_back("SS2");
|
|
params.push_back("SS3");
|
|
params.push_back("SS4");
|
|
params.push_back("SS5");
|
|
return params;
|
|
}
|
|
std::vector<std::string> param_string_sWave(int folding){
|
|
std::vector<std::string> params;
|
|
params.push_back("FS");
|
|
if(folding == 1) params.push_back("SS2");
|
|
if(folding == 2) params.push_back("SS3");
|
|
if(folding > 2) params.push_back("SS4");
|
|
if(folding < 4) params.push_back("SS1");
|
|
return params;
|
|
}
|
|
|
|
//TODO: merge with parNamesWithTex
|
|
std::vector<std::string> param_string_bkg(){
|
|
std::vector<std::string> str;
|
|
for (auto ang: ANGLES){
|
|
for (int i = 1; i < 5; i++){
|
|
str.push_back("cbkg"+ang+std::to_string(i));
|
|
}
|
|
}
|
|
return str;
|
|
}
|
|
|
|
std::vector<std::string> param_string_bkg_mkpi(){
|
|
std::vector<std::string> str;
|
|
for (int i = 1; i < 2; i++){
|
|
str.push_back("cbkgmkpi"+std::to_string(i));
|
|
}
|
|
return str;
|
|
}
|
|
|
|
std::vector<std::string> param_string(fcnc::options opts, bool MC){
|
|
std::vector<std::string> params;
|
|
//This is suboptimal but I feel it is easier for the reader
|
|
if (opts.full_angular){
|
|
if(opts.fit_pprimes) params = param_string_pPrimes(opts.fit_fl);
|
|
else params = param_string_s(opts.fit_fl, opts.fit_afb);
|
|
if(opts.swave){
|
|
std::vector<std::string> sWave = param_string_sWave();
|
|
params.insert(params.end(), sWave.begin(), sWave.end());
|
|
}
|
|
if(!MC){
|
|
std::vector<std::string> bkg = PAR_BKG_STRING(-1,opts.bkg_order_costhetal,opts.bkg_order_costhetak); //Returns only varied parameters
|
|
params.insert(params.end(), bkg.begin(), bkg.end());
|
|
}
|
|
return params;
|
|
}
|
|
else{
|
|
if(opts.fit_pprimes) params = param_string_pPrimes(opts.fit_fl,opts.folding);
|
|
else params = param_string_s(opts.fit_fl, opts.fit_afb,opts.folding);
|
|
if(opts.swave){
|
|
std::vector<std::string> sWave = param_string_sWave(opts.folding);
|
|
params.insert(params.end(), sWave.begin(), sWave.end());
|
|
}
|
|
if(!MC){
|
|
std::vector<std::string> bkg = PAR_BKG_STRING(opts.folding,opts.bkg_order_costhetal,opts.bkg_order_costhetak);
|
|
params.insert(params.end(), bkg.begin(), bkg.end());
|
|
}
|
|
return params;
|
|
}
|
|
}
|
|
|
|
std::vector<std::string> params_string_mass(fcnc::options opts){
|
|
std::vector<std::string> params;
|
|
params.push_back("m_b");
|
|
if (!opts.twotailedcrystalball) params.push_back("m_res_1");
|
|
params.push_back("m_sigma_1");
|
|
if (!opts.twotailedcrystalball) params.push_back("m_sigma_2");
|
|
params.push_back("alpha_1");
|
|
params.push_back("alpha_2");
|
|
params.push_back("n_1");
|
|
params.push_back("n_2");
|
|
return params;
|
|
}
|
|
|
|
fcnc::parameter* par_with_correct_name(int pos, std::string parname, fcnc::parameters* theParams){
|
|
//pos is the position of the parameter I want to save
|
|
|
|
fcnc::parameter* par = theParams->get_parameter(pos);
|
|
|
|
//I am not sure, but from the code it seems the parameters should be at the same spot for all pdfs
|
|
//Therefore first check if the name on the spot is correct
|
|
if (par->get_name() == parname.c_str()) return par;
|
|
|
|
unsigned int paridx = 0;
|
|
while(paridx != theParams->nparameters()){
|
|
par = theParams->get_parameter(paridx);
|
|
if (par->get_name() == parname.c_str()) return par;
|
|
paridx++;
|
|
}
|
|
spdlog::warn("Parameter "+parname+" not found in PDF {0:d}. Continue with next PDF.");
|
|
return NULL;
|
|
}
|
|
|
|
|
|
double get_sigmaRatio_fromMC(basic_params params, int nBins, int bin, int pdf){
|
|
|
|
std::string sigma_name = "m_sigma_1"; //In case one needs m_sigma or m_sigma_2 or so
|
|
std::string fileName_sig = final_result_name_MC(params, nBins, false, false, true, false, false);
|
|
std::string fileName_ref = final_result_name_MC(params, 1, true, false, true, false, true);
|
|
|
|
return double(get_param_value_from_rootfile(fileName_sig,sigma_name,pdf,bin))/ double(get_param_value_from_rootfile(fileName_ref,sigma_name,pdf,0));
|
|
}
|
|
|
|
//--------------------------------------//
|
|
// Helpers for printing fit results //
|
|
//--------------------------------------//
|
|
|
|
void print_all_parameters(unsigned int nBins, fcnc::bu2kstarmumu_parameters *theParams[],
|
|
int spdlog_level, std::string savePath){
|
|
//Check whether the desired level >= current level
|
|
if (spdlog_level >= spdlog::default_logger_raw()->level()){
|
|
//If yes, print the parameters at info level
|
|
|
|
//Save current level of spdlog
|
|
spdlog::level::level_enum tmp_log_level = spdlog::default_logger_raw()->level();
|
|
|
|
//Set the level to info so it is printed for sure
|
|
spdlog::set_level(spdlog::level::info);
|
|
spdlog::info("");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| PARAMETERS ||");
|
|
spdlog::info("==========================================");
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
spdlog::info("|| BIN {0:d}", b);
|
|
if(spdlog::default_logger_raw()->level()<=spdlog_level){
|
|
theParams[b]->print_parameters(false);
|
|
}
|
|
if (savePath != ""){
|
|
std::string saveBinPath = savePath + "_bin" + std::to_string(b) + ".txt";
|
|
theParams[b]->save_param_values(saveBinPath + "");
|
|
}
|
|
theParams[b]->print_parameters(false);
|
|
}
|
|
//Reset the level back where it was
|
|
spdlog::set_level(tmp_log_level);
|
|
}
|
|
else return;
|
|
|
|
|
|
}
|
|
|
|
|
|
void print_all_parameters(unsigned int nBins, std::vector<UInt_t> pdf_idx,
|
|
std::vector<fcnc::parameters*> theParams[],
|
|
int spdlog_level){
|
|
//Check whether the desired level >= current level
|
|
if (spdlog_level >= spdlog::default_logger_raw()->level()){
|
|
//If yes, print the parameters at info level
|
|
|
|
//Save current level of spdlog
|
|
spdlog::level::level_enum tmp_log_level = spdlog::default_logger_raw()->level();
|
|
|
|
//Set the level to trace so it is printed for sure
|
|
//Pritn_parameters needs trace
|
|
spdlog::set_level(spdlog::level::trace);
|
|
|
|
spdlog::info("");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| PARAMETERS ||");
|
|
spdlog::info("==========================================");
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
spdlog::info("|| BIN {0:d}\tPDF {1:d}", b, pdf_idx.at(i));
|
|
theParams[b].at(i)->print_parameters(false);
|
|
}
|
|
}
|
|
//Reset the level back where it was
|
|
spdlog::set_level(tmp_log_level);
|
|
}
|
|
else return;
|
|
}
|
|
|
|
void tex_all_parameters(unsigned int nBins, std::vector<UInt_t> pdf_idx,
|
|
std::vector<fcnc::parameters*> theParams[]){
|
|
std::ofstream myFile;
|
|
open_Latex_noteFile(latex_params(), myFile); //open file with parameters
|
|
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
myFile << "\\captionof{table}{PDF " << i << " bin " << b<< "}" << std::endl;
|
|
theParams[b].at(i)->print_parameters(true);
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
void tex_all_parameters(unsigned int nBins, int n_pdf, fcnc::bu2kstarmumu_parameters **theParams){
|
|
std::ofstream myFile;
|
|
open_Latex_noteFile(latex_params(), myFile); //open file with parameters
|
|
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
myFile << "\\captionof{table}{PDF " << n_pdf<< " bin " << b<< "}" << std::endl;
|
|
theParams[b]->print_parameters(true);
|
|
}
|
|
return;
|
|
}
|
|
|
|
|
|
void print_fit_results(unsigned int nBins,std::vector<Int_t> fitresults){
|
|
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| FIT RESULTS ||");
|
|
spdlog::info("==========================================");
|
|
if(spdlog_info()){
|
|
spdlog::info("|| BIN\t\tfit result ||");
|
|
for(UInt_t b = 0; b < nBins; b++){
|
|
spdlog::info("|| {0:d}\t\t{1:d}\t||", b, fitresults.at(b));
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
void print_fit_results(unsigned int nBins, std::vector<UInt_t> pdf_idx, std::vector<Int_t> fitresults[], bool simFit){
|
|
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| FIT RESULTS ||");
|
|
spdlog::info("==========================================");
|
|
if(spdlog_info()){
|
|
spdlog::info("|| BIN\t\tPDF\t\tfit result ||");
|
|
for(UInt_t b = 0; b < nBins; b++){
|
|
if(!simFit){
|
|
for_indexed(auto idx: pdf_idx){
|
|
spdlog::info("|| {0:d}\t\t{1:d}\t\t{2:d} \t||\n", b, idx, fitresults[b].at(i));
|
|
}
|
|
}
|
|
else{
|
|
spdlog::info("|| {0:d}\t\t All PDFs\t\t{1:d} \t||\n", b, fitresults[b].at(0)); //TODO
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void print_sig_yields(unsigned int nBins, std::vector<UInt_t> pdf_idx,
|
|
std::vector<UInt_t> *evts_cntr, std::vector<double> *f_sigs, std::vector<double> *f_sigserr){
|
|
spdlog::info("");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| SIGNAL YIELDS ||");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| BIN\t\tPDF\tsignal yield ||");
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
spdlog::info("|| {0:d}\t\t{1:d}\t{2:f}+/-{3:f}\t||", b, pdf_idx.at(i), evts_cntr[b].at(i)*f_sigs[b].at(i), evts_cntr[b].at(i)*f_sigserr[b].at(i));
|
|
}
|
|
}
|
|
spdlog::info("");
|
|
}
|
|
|
|
void print_bkg_yields(unsigned int nBins, std::vector<UInt_t> pdf_idx,
|
|
std::vector<UInt_t> *evts_cntr, std::vector<double> *f_sigs, std::vector<double> *f_sigserr){
|
|
spdlog::info("");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| BACKGROUND YIELDS ||");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| BIN\t\tPDF\tbkg yield ||");
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
spdlog::info("|| {0:d}\t\t{1:d}\t{2:.1f}+/-{3:.1f}\t ||", b, pdf_idx.at(i), evts_cntr[b].at(i)*(1 - f_sigs[b].at(i)), evts_cntr[b].at(i)*f_sigserr[b].at(i));
|
|
}
|
|
}
|
|
spdlog::info("");
|
|
|
|
}
|
|
|
|
|
|
void print_bkgOnly_yields(unsigned int nBins, std::vector<UInt_t> pdf_idx,
|
|
std::vector<fcnc::bu2kstarmumu_parameters*> theParams[],
|
|
std::vector<double> bkg_int_full_range[], std::vector<UInt_t> *evts_cntr){
|
|
spdlog::info("");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| BACKGROUND YIELDS ||");
|
|
spdlog::info("==========================================");
|
|
spdlog::info("|| BIN\t\tPDF\tbkg yield ||");
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
spdlog::info("[BIN{0:d}]\tBckgnd est.: {1:.4f}\tBkg prob: {2:.4f}\tBkg iac: {3:.4f}\tExpCoeff: {4:.4f}",
|
|
i,
|
|
evts_cntr[b].at(i) * bkg_int_full_range[i].at(b) * (1 - theParams[b].at(i)->f_sig.get_value()),
|
|
bkg_int_full_range[b].at(i),
|
|
1 - theParams[b].at(i)->f_sig.get_value(),
|
|
theParams[b].at(i)->m_lambda.get_value());
|
|
}
|
|
}
|
|
spdlog::info("");
|
|
|
|
}
|
|
|
|
|
|
void print_sig_yields_tex(std::string texFiletag, unsigned int nBins, std::vector<UInt_t> pdf_idx,
|
|
fcnc::options *theOptions,
|
|
std::vector<UInt_t> *evts_cntr, std::vector<double> *f_sigs, std::vector<double> *f_sigserr){
|
|
|
|
//Open texFile
|
|
std::ofstream myFile;
|
|
open_Latex_noteFile(texFiletag, myFile);
|
|
|
|
double yield[nBins], yielderr[nBins];
|
|
|
|
myFile << "\\begin{tabular}{|l|"; //TODO: have a function for this
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++) myFile << "r";
|
|
myFile << "|r|}\\hline" << std::endl;
|
|
myFile << "\\qsq bin";
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
myFile << "\t& Run" << theOptions[i].run;
|
|
}
|
|
myFile << "\\\\" << std::endl;
|
|
myFile << "\\hline\\hline" << std::endl;
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
myFile << b << "\t&";
|
|
yield[b] = 0.0; yielderr[b] = 0.0;
|
|
for(UInt_t i = 0; i < pdf_idx.size(); i++){
|
|
myFile << "$" << std::setprecision(1) << std::fixed << evts_cntr[b].at(i)*f_sigs[b].at(i) << "\\pm" << evts_cntr[b].at(i)*f_sigserr[b].at(i) << "$ \t&";
|
|
yield[b] += evts_cntr[b].at(i)*f_sigs[b].at(i);
|
|
yielderr[b] += evts_cntr[b].at(i)*f_sigserr[b].at(i)*evts_cntr[b].at(i)*f_sigserr[b].at(i);
|
|
}
|
|
myFile << "$" << yield[b] << "\\pm" << TMath::Sqrt(yielderr[b]) << "$\\\\" << std::endl;
|
|
}
|
|
myFile << "\\hline\\end{tabular}" << std::endl << std::endl << std::endl << std::endl;
|
|
|
|
myFile << "\\begin{tabular}{|l|r|}\\hline" << std::endl;
|
|
myFile << "\\qsq bin \t& signal yield\\\\" << std::endl;
|
|
myFile << "\\hline\\hline" << std::endl;
|
|
for(unsigned int b = 0; b < nBins; b++) myFile << b << "\t& $" << yield[b] << "\\pm" << TMath::Sqrt(yielderr[b]) << "$\\\\" << std::endl;
|
|
myFile << "\\hline\\end{tabular}" << std::endl << std::endl;
|
|
myFile.close();
|
|
}
|
|
|
|
|
|
int save_results(std::string results_file, unsigned int nBins, std::vector<UInt_t> pdf_idx, std::vector<int>*fit_results, std::vector<fcnc::parameters*> *theParams, bool simFit, fcnc::options *opts){
|
|
|
|
//Open the root file
|
|
spdlog::info("[SAVE]\t\tSaving result values to file " + results_file);
|
|
TFile* fout = new TFile(results_file.c_str(), "RECREATE");
|
|
fout->cd();
|
|
|
|
int param_index = 0;
|
|
|
|
//loop over parameters:
|
|
for(UInt_t p = 0; p < theParams[0].at(0)->nparameters(); p++){
|
|
fcnc::parameter* param = theParams[0].at(0)->get_parameter(p);
|
|
//if(param->get_step_size() == 0.0) continue; //skip fixed parameters
|
|
std::string parname(param->get_name());
|
|
std::string pardesc(param->get_description());
|
|
spdlog::trace("Creating new TTree for parameter:" + parname);
|
|
TTree* t = new TTree(parname.c_str(), pardesc.c_str());
|
|
int pdf = DEFAULT_TREE_INT;
|
|
int bin = DEFAULT_TREE_INT;
|
|
int migrad = DEFAULT_TREE_INT;
|
|
int status_cov = DEFAULT_TREE_INT;
|
|
int totBins = DEFAULT_TREE_INT;
|
|
|
|
double q2a = DEFAULT_TREE_VAL;
|
|
double q2b = DEFAULT_TREE_VAL;
|
|
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;
|
|
|
|
const unsigned int corr_max = param->correlations.size();
|
|
double tmp_corr[corr_max];
|
|
|
|
t->Branch("pdf",&pdf,"pdf/I");
|
|
t->Branch("bin",&bin,"bin/I");
|
|
t->Branch("totBins",&totBins,"totBins/I");
|
|
t->Branch("migrad",&migrad,"migrad/I");
|
|
t->Branch("status_cov",&status_cov,"status_cov/I");
|
|
|
|
t->Branch("value",&value,"value/D");
|
|
t->Branch("q2min",&q2a,"q2min/D");
|
|
t->Branch("q2max",&q2b,"q2max/D");
|
|
t->Branch("error",&error,"error/D");
|
|
t->Branch("error_up",&error_up,"error_up/D");
|
|
t->Branch("error_down",&error_down,"error_down/D");
|
|
t->Branch("start_value",&start_value,"start_value/D");
|
|
t->Branch("prev_value",&prev_value,"prev_value/D");
|
|
t->Branch("prev_error",&prev_error,"prev_error/D");
|
|
|
|
t->Branch("index",¶m_index, "index/I");
|
|
t->Branch("correlations",&tmp_corr,("correlations["+std::to_string(corr_max)+"]/D").c_str());
|
|
|
|
//loop over PDFs
|
|
for_indexed(auto n_pdf: pdf_idx){
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
fcnc::parameter* par = par_with_correct_name(p, parname, theParams[b].at(i));
|
|
if (par == NULL) break;
|
|
totBins = nBins;
|
|
bin = b;
|
|
pdf = n_pdf;
|
|
migrad = fit_results[b].at(simFit ? 0 : i) % 100;
|
|
status_cov = fit_results[b].at(simFit ? 0 : i) / 100;
|
|
|
|
q2a = opts->TheQ2binsmin[b];
|
|
q2b = opts->TheQ2binsmax[b];
|
|
value = par->get_value();
|
|
error = par->get_error();
|
|
error_up = par->get_error_up();
|
|
error_down = par->get_error_down();
|
|
if (opts->minos_errors){
|
|
if (error_up==0.0 && value < par->get_max()) error_up = par->get_max() - value;
|
|
if (error_down==0.0 && value > par->get_min()) error_down = par->get_min() - value;//needs to be negative
|
|
}
|
|
start_value = par->get_start_value();
|
|
prev_value = par->get_previous_measurement();
|
|
prev_error = par->get_previous_error();
|
|
|
|
for (unsigned int l = 0; l < corr_max; l++) tmp_corr[l] = 0.0; //Fill with zeroes
|
|
for (unsigned int l = 0; l < par->correlations.size(); l++){
|
|
tmp_corr[l] = par->correlations.at(l);
|
|
}
|
|
t->Fill();
|
|
}
|
|
} //end loop over pdfs
|
|
|
|
//Write the parameter tree
|
|
t->Write();
|
|
delete t;
|
|
param_index++;
|
|
}
|
|
|
|
//Close the root file
|
|
spdlog::info("Finished saving the parameters into the file {0:s}.",fout->GetName());
|
|
fout->Close();
|
|
delete fout;
|
|
|
|
//Delete the tex file with parameters
|
|
clear_Latex_noteFile(latex_params());
|
|
//Copy the parameter in TeX format into the same place as the final results file
|
|
tex_all_parameters(nBins,pdf_idx,theParams);
|
|
replace(results_file,".root",".tex");
|
|
copyFile(get_Latex_noteFile(latex_params()),results_file);
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
int save_results(std::string results_file, unsigned int nBins, int n_pdf, std::vector<int> fit_results, fcnc::bu2kstarmumu_parameters **theParams, fcnc::options *opts){
|
|
|
|
//Open the root file
|
|
spdlog::info("[SAVE]\t\tSaving result values to file " + results_file);
|
|
TFile* fout = new TFile(results_file.c_str(), "RECREATE");
|
|
fout->cd();
|
|
|
|
int param_index = 0;
|
|
|
|
//loop over parameters:
|
|
for(UInt_t p = 0; p < theParams[0]->nparameters(); p++){
|
|
fcnc::parameter* param = theParams[0]->get_parameter(p);
|
|
//if(param->get_step_size() == 0.0) continue; //skip fixed parameters
|
|
std::string parname(param->get_name());
|
|
std::string pardesc(param->get_description());
|
|
spdlog::trace("Creating new TTree for parameter:" + parname);
|
|
TTree* t = new TTree(parname.c_str(), pardesc.c_str());
|
|
int pdf = DEFAULT_TREE_INT;
|
|
int bin = DEFAULT_TREE_INT;
|
|
int migrad = DEFAULT_TREE_INT;
|
|
int status_cov = DEFAULT_TREE_INT;
|
|
int totBins = DEFAULT_TREE_INT;
|
|
|
|
double q2a = DEFAULT_TREE_VAL;
|
|
double q2b = DEFAULT_TREE_VAL;
|
|
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;
|
|
|
|
const unsigned int corr_max = param->correlations.size();
|
|
double tmp_corr[corr_max];
|
|
|
|
t->Branch("pdf",&pdf,"pdf/I");
|
|
t->Branch("bin",&bin,"bin/I");
|
|
t->Branch("totBins",&totBins,"totBins/I");
|
|
t->Branch("migrad",&migrad,"migrad/I");
|
|
t->Branch("status_cov",&status_cov,"status_cov/I");
|
|
|
|
t->Branch("value",&value,"value/D");
|
|
t->Branch("q2min",&q2a,"q2min/D");
|
|
t->Branch("q2max",&q2b,"q2max/D");
|
|
t->Branch("error",&error,"error/D");
|
|
t->Branch("error_up",&error_up,"error_up/D");
|
|
t->Branch("error_down",&error_down,"error_down/D");
|
|
t->Branch("start_value",&start_value,"start_value/D");
|
|
t->Branch("prev_value",&prev_value,"prev_value/D");
|
|
t->Branch("prev_error",&prev_error,"prev_error/D");
|
|
|
|
t->Branch("index",¶m_index, "index/I");
|
|
t->Branch("correlations",&tmp_corr,("correlations["+std::to_string(corr_max)+"]/D").c_str());
|
|
|
|
//loop over nBins
|
|
for(unsigned int b = 0; b < nBins; b++){
|
|
fcnc::parameter* par = par_with_correct_name(p, parname, theParams[b]);
|
|
if (par == NULL) break;
|
|
totBins = nBins;
|
|
bin = b;
|
|
pdf = n_pdf;
|
|
migrad = fit_results[b] % 100;
|
|
status_cov = fit_results[b] / 100;
|
|
|
|
q2a = opts->TheQ2binsmin[b];
|
|
q2b = opts->TheQ2binsmax[b];
|
|
value = par->get_value();
|
|
error = par->get_error();
|
|
error_up = par->get_error_up();
|
|
error_down = par->get_error_down();
|
|
if (opts->minos_errors){
|
|
if (error_up==0.0 && value < par->get_max()) error_up = par->get_max() - value;
|
|
if (error_down==0.0 && value > par->get_min()) error_down = par->get_min() - value;//needs to be negative
|
|
}
|
|
start_value = par->get_start_value();
|
|
prev_value = par->get_previous_measurement();
|
|
prev_error = par->get_previous_error();
|
|
|
|
for (unsigned int l = 0; l < corr_max; l++) tmp_corr[l] = 0.0; //Fill with zeroes
|
|
for (unsigned int l = 0; l < par->correlations.size(); l++){
|
|
tmp_corr[l] = par->correlations.at(l);
|
|
}
|
|
t->Fill();
|
|
|
|
}
|
|
//Write the parameter tree
|
|
t->Write();
|
|
delete t;
|
|
param_index++;
|
|
} //End loop over names
|
|
|
|
//Close the root file
|
|
spdlog::info("Finished saving the parameters into the file {0:s}.",fout->GetName());
|
|
fout->Close();
|
|
|
|
//Delete the tex file with parameters
|
|
clear_Latex_noteFile(latex_params());
|
|
//Copy the parameter in TeX format into the same place as the final results file
|
|
tex_all_parameters(nBins,n_pdf,theParams);
|
|
replace(results_file,".root",".tex");
|
|
copyFile(get_Latex_noteFile(latex_params()),results_file);
|
|
|
|
delete fout;
|
|
return 0;
|
|
}
|
|
|
|
|
|
//--------------------------------------//
|
|
// Helpers for writting latex stuff //
|
|
//--------------------------------------//
|
|
int clear_Latex_noteFile(std::string tag){
|
|
if( remove(get_Latex_noteFile(tag).c_str()) != 0 ) spdlog::warn( "Unable to delete file " + get_Latex_noteFile(tag) ); //Do not return an error status cause the file might just not exist, which is ok
|
|
else spdlog::trace( "Latex noteFile successfully deleted" );
|
|
return 0;
|
|
}
|
|
|
|
int open_Latex_noteFile(std::string tag, std::ofstream &myfile){
|
|
//This is not working cause this thing is for whatever efin reason stuck with gcc 4.
|
|
myfile.open (get_Latex_noteFile(tag), std::ios::out | std::ios::app);
|
|
if (!myfile.is_open()){
|
|
spdlog::error("Latex noteFile" + get_Latex_noteFile(tag) + " is not opened.");
|
|
return 404;
|
|
}
|
|
else{
|
|
spdlog::trace("Latex noteFile successfully deleted");
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
//--------------------------------------//
|
|
// Helpers for reading the parameters //
|
|
//--------------------------------------//
|
|
|
|
//Checks if the observable is in the TFile
|
|
//If required Fl and there is S1s, it returns 11, ...
|
|
//Modulo 10 then returns whether the observable actually exists
|
|
//Example usage: if (try_getObservable(paramName, file)%10==0) assert(0)
|
|
int try_getObservable(std::string observable, TFile* file){
|
|
if (!file->GetListOfKeys()->Contains(observable.c_str())){
|
|
//Also check if Fl<->S1s and A_FB<->S6s is available
|
|
if (observable == "Fl") return 10+try_getObservable("S1s",file);
|
|
if (observable == "S1s") return 20+try_getObservable("Fl",file);
|
|
if (observable == "Afb") return 30+try_getObservable("S6s",file);
|
|
if (observable == "S6s") return 40+try_getObservable("Afb",file);
|
|
return 0;
|
|
}
|
|
else return 1;
|
|
}
|
|
|
|
std::vector<double> load_param_values_into_vector(std::string paramName, std::string fileName){
|
|
//There is a function for this in bu2kstarmumu_parameters.cc, but that is inherited from parameters, therefore a parameter class needs to be created for this
|
|
//That would be tedious, so there is a new class that opens the root file with results and reads it from a tree directly into an array
|
|
|
|
//Create the vector that will be returned
|
|
std::vector<double> output;
|
|
|
|
//Open file
|
|
spdlog::info("Opening " + fileName);
|
|
TFile* file = new TFile(fileName.c_str(), "READ");
|
|
|
|
int isFoundFlag = try_getObservable(paramName, file);
|
|
if (!(isFoundFlag%10)){
|
|
spdlog::critical("Wrong tree name " + paramName + "! Returning an empty vector!");
|
|
//It should not crash here as it might be desired to init some observables but not use them
|
|
return {};
|
|
}
|
|
else{
|
|
switch(isFoundFlag){
|
|
case 11: //Need Fl, have S1s
|
|
paramName = "S1s";
|
|
break;
|
|
case 21: //Need S1s, have FL
|
|
paramName = "Fl";
|
|
break;
|
|
case 31: //Need Afb, have S6s
|
|
paramName = "S6s";
|
|
break;
|
|
case 41: //Need S6s, have Afb
|
|
paramName = "Afb";
|
|
break;
|
|
}
|
|
//case 1: all good, move one
|
|
}
|
|
|
|
//Get the tree
|
|
TTree* tree = (TTree*)file->Get(paramName.c_str());
|
|
//Activate only needed branches
|
|
tree->SetBranchStatus("*",0);
|
|
tree->SetBranchStatus("value",1);
|
|
double value = DEFAULT_TREE_VAL;
|
|
tree->SetBranchAddress("value", &value);
|
|
|
|
//Loop over nBins that are saved in the tree
|
|
for (int b = 0; b < tree->GetEntries(); b++){
|
|
tree->GetEntry(b);
|
|
//There are not many bins so just do a switch and fill accordingly
|
|
switch(isFoundFlag){
|
|
case 1:
|
|
output.push_back(value);
|
|
case 11: //Need Fl, have S1s
|
|
output.push_back(1-4.0/3.0*value);
|
|
break;
|
|
case 21: //Need S1s, have FL
|
|
output.push_back(3.0/4.0*(1-value));
|
|
break;
|
|
case 31: //Need Afb, have S6s
|
|
output.push_back(3.0/4.0*value);
|
|
break;
|
|
case 41: //Need S6s, have Afb
|
|
output.push_back(4.0/3.0*value);
|
|
break;
|
|
}
|
|
}
|
|
|
|
spdlog::debug(convert_vector_to_string(output));
|
|
return output;
|
|
}
|