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.
138 lines
6.1 KiB
138 lines
6.1 KiB
//Renata Kopecna
|
|
|
|
#include <TFile.h>
|
|
#include <TCanvas.h>
|
|
#include <TChain.h>
|
|
#include <TF1.h>
|
|
#include <TH1D.h>
|
|
|
|
#include "GetMeanError.hh"
|
|
#include <design.hh>
|
|
#include <helpers.hh>
|
|
#include "ScriptHelpers.hh"
|
|
|
|
|
|
std::vector<double> getSM(std::string observable){
|
|
|
|
if (observable == "Fl") return {+0.248569008046490070, +0.2581971160208132000, +0.42526014659559220000, +0.49420669353506260000, +0.1850102996703149000}; //s1s
|
|
else if (observable == "S3") return {+0.000732569649322501, -0.0324523859139632900, -0.08563864422457230000, -0.18934424485647136000, -0.0130310503598977030}; //s3
|
|
else if (observable == "S4") return {-0.028317693174889150, -0.2425837027660650600, -0.28116530113868093000, -0.29668531019781386000, -0.1468982883895949300}; //s4
|
|
else if (observable == "S5") return {+0.035479698896675250, -0.3745017565795041000, -0.40726931726860705000, -0.30161750053585110000, -0.1915219985323076700}; //s5
|
|
else if (observable == "Afb")return {-0.122319885918514240, +0.2473816826915200600, +0.52432045165541470000, +0.52106439209368990000, +0.0147578966169414490}; //s6s
|
|
else if (observable == "S7") return {-0.019468600975846337, -0.0105435770924197210, -0.00219285191192886900, -0.00119446595061885200, -0.0160492529928752330}; //s7
|
|
else if (observable == "S8") return {-0.010192773160727949, -0.0043019891737210840, +0.00104987431866559550, +0.00027244645255727460, -0.0070279543261629670}; //s8
|
|
else if (observable == "S9") return {-0.000756965302130655, -0.0006959911482550733, +0.00044741790607562027, +0.00026464193866571387, -0.0007199408885633002}; //s9
|
|
else{
|
|
spdlog::warn("Wrong observable, returning an empty vector");
|
|
return {};
|
|
}
|
|
}
|
|
|
|
std::vector<double> fixSM(std::string observable){
|
|
std::vector<double> tmp = {};
|
|
for (auto val : getSM(observable)){
|
|
if (observable == "Fl") tmp.push_back(round((1-4.0*val/3.0)*100.0)/100.0);
|
|
else if (observable=="Afb") tmp.push_back(round((3.0*val/4.0)*100.0)/100.0);
|
|
else tmp.push_back(round(val*100.0)/100.0);
|
|
}
|
|
return tmp;
|
|
}
|
|
|
|
int loadAllFiles(std::string observable, basic_params params){
|
|
//load chain with all results of bootstrapping
|
|
TChain * ch = new TChain(observable.c_str());
|
|
|
|
std::string files = final_result_name_toys(1111, params.reference, params.nBins, true, params, params.Run, false, false, false, false, false, true);
|
|
replace(files, "1111", "*"); //integer replaced by a string *
|
|
replace(files, "ToysFit/", "ToysFit/"+std::to_string(params.jobID)+"/");
|
|
|
|
spdlog::debug("Loading files: " + files);
|
|
|
|
ch->Add(files.c_str());
|
|
int nEntries = ch->GetEntries();
|
|
spdlog::debug("Loaded {0:d} entries", nEntries);
|
|
if (nEntries == 0){
|
|
spdlog::error("No files found!");
|
|
return 404;
|
|
}
|
|
|
|
//link variables to branches
|
|
|
|
double value = DEFAULT_TREE_VAL;
|
|
double error = DEFAULT_TREE_ERR;
|
|
double errorup = DEFAULT_TREE_ERR;
|
|
double errordown = DEFAULT_TREE_ERR;
|
|
|
|
int migrad = DEFAULT_TREE_INT;
|
|
int cov = DEFAULT_TREE_INT;
|
|
int bin = DEFAULT_TREE_INT;
|
|
int pdf = DEFAULT_TREE_INT;
|
|
|
|
ch->SetBranchStatus("*",1); //TODO optimize
|
|
ch->SetBranchAddress("value", &value);
|
|
ch->SetBranchAddress("error", &error);
|
|
ch->SetBranchAddress("error_up", &errorup);
|
|
ch->SetBranchAddress("error_down", &errordown);
|
|
ch->SetBranchAddress("migrad", &migrad);
|
|
ch->SetBranchAddress("status_cov", &cov);
|
|
ch->SetBranchAddress("bin", &bin);
|
|
ch->SetBranchAddress("pdf", &pdf);
|
|
|
|
const int nBins = params.nBins;
|
|
std::vector<double> errors[nBins]; //Don't fill a histogram just yet, just put the values in a vector first
|
|
|
|
//Loop over the entries and fill the vector
|
|
for (int iter = 0; iter < nEntries; iter++){
|
|
ch->GetEntry(iter);
|
|
//require converged migrad
|
|
if (migrad!=0) continue;
|
|
//require fitresult = 300 or 100 (or 200)
|
|
if(!(cov == 1 || cov == 3 || cov == 2)) continue;
|
|
if (pdf == 1) continue; //Just take one pdf, as the values are shared
|
|
spdlog::trace("Parameter value at entry {0:d}: {1:f}", iter, value);
|
|
spdlog::trace("Parameter error at entry {0:d}: {1:f}", iter, error);
|
|
errors[bin].push_back(error);
|
|
}
|
|
|
|
//Debug print in case
|
|
for (int b = 0; b < nBins; b++){
|
|
spdlog::debug("Got {0:d} entries for bin {1:d}", errors[b].size(),b);
|
|
}
|
|
|
|
//Now init a histogram for each bin and fill it
|
|
std::vector<double> means;
|
|
for (int b = 0; b < nBins; b++){
|
|
//Get the max and the min of the errors
|
|
double max = *max_element(errors[b].begin(), errors[b].end());
|
|
double min = *min_element(errors[b].begin(), errors[b].end());
|
|
if (min == max) continue; //If errors are all the same, it was fixed, so don't plot
|
|
//Easier than checking folding
|
|
if (max>1.0) max = 0.5;
|
|
if (min<0.01) min = 0.01;
|
|
//Create the TH1D for the errors
|
|
TH1D *h_errors = new TH1D(observable.c_str(),
|
|
(observable+" in q^{2} bin #"+std::to_string(b)
|
|
+";"+observable+";Entries").c_str(),
|
|
errors[b].size()/10, min, max);
|
|
//Fill the TH1D from the read vector
|
|
for (auto err: errors[b]) h_errors->Fill(err);
|
|
//Fit the distribution with a gaussia
|
|
TF1 *fit_func = new TF1("fit_func","gaus",min,max);
|
|
fit_func->SetParameters(h_errors->GetMaximum(), h_errors->GetMean(), h_errors->GetRMS() );
|
|
h_errors->Fit("fit_func",spdlog_debug() ? "R" : "RQ");
|
|
//Plot and save it
|
|
plotAndSave(h_errors,observable,"./"+observable+"_"+std::to_string(params.jobID)+"_"+std::to_string(b),"eps");
|
|
//Save the mean to a vector so you can print it later
|
|
means.push_back(fit_func->GetParameter("Mean") );
|
|
h_errors->Clear();
|
|
}
|
|
if (means.size()>0){
|
|
std::vector<double>valSM = fixSM(observable);
|
|
std::cout << "$" << observable << "$\t";
|
|
for_indexed (auto val: means) std::cout << " & $" << valSM[i] << " \\pm " << val << "$";
|
|
std::cout << "\\\\" << std::endl;
|
|
}
|
|
return 0;
|
|
|
|
|
|
}
|