|
|
//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;
}
|
