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.
628 lines
28 KiB
628 lines
28 KiB
//Renata Kopecna
|
|
|
|
|
|
#include <pulls.hh>
|
|
|
|
#include <string>
|
|
#include <sstream>
|
|
#include <iostream>
|
|
|
|
#include "folder.hh"
|
|
#include <fitter.hh>
|
|
#include <bu2kstarmumu_plotter.hh>
|
|
#include <bu2kstarmumu_generator.hh>
|
|
#include <bu2kstarmumu_pdf.hh>
|
|
#include <bu2kstarmumu_parameters.hh>
|
|
#include <paths.hh>
|
|
#include <helpers.hh>
|
|
#include <event.hh>
|
|
|
|
#include <TTree.h>
|
|
|
|
#include <spdlog.h>
|
|
|
|
//////////////////////////
|
|
// INIT VALUES FOR FITS //
|
|
//////////////////////////
|
|
|
|
//Do they change anywhere or can I set them in a header?
|
|
//TODO replace by the parameters->init
|
|
double f_m_B[4] = {5282.32, 5282.08, 5278.57, 5279.48};
|
|
double f_lambda1[4] = {-0.00386, -0.00610, -0.00362, -0.00578};
|
|
double f_lambda2[4] = {-0.03320, -0.00120, -0.04600, -0.00100};
|
|
double f_alpha_1[4] = { 1.7910, 1.8000, 1.6510, 1.5400};
|
|
double f_alpha_2[4] = { 2.3200, 2.4300, 2.0100, 2.2500};
|
|
double f_n_1[4] = { 1.4400, 1.4900, 1.6800, 1.6100};
|
|
double f_n_2[4] = { 2.8900, 3.7000, 3.8700, 3.7000};
|
|
double f_sigma_1[4] = {14.8300, 15.4800, 14.8400, 14.4};
|
|
double f_sigma_2[4] = {14.0700, 13.1000, 15.6200, 13.3};
|
|
double f_f_CB[4] = { 0.6360, 0.3600, 0.3040, 0.6200};
|
|
double f_m_scale[4] = { 1.1814, 1.0648, 1.1547, 1.1168};
|
|
|
|
|
|
//TODO posibly to be deleted so ignoring for now!
|
|
void do_pulls(fcnc::options opts, bool doToyPulls, bool doMCPulls, int job_id, double SwaveFraction, UInt_t nPulls, Double_t nToyEvents, UInt_t nMCEvents ,std::string angularsuffix, bool Fit1bin, bool Fit2bins, bool FitAllbins ){ //TODO: change to int and return some errors possibly
|
|
//TODO: figure out jobs_to_do and job_id, but there is no manipulation with it, so so far so good
|
|
//Why is it not opts.job_id?
|
|
if(doToyPulls)std::cout << "[INFO]\t\tGetting pull distributions and test fits from Toy Events!" << std::endl;
|
|
|
|
//disable eps output if job_id is larger than -1. This is needed for HD cluster jobs:
|
|
if(job_id > -1)opts.write_eps = false;
|
|
|
|
//have extra option for the usage of binning:
|
|
const bool UseBinnedFit = false;
|
|
|
|
opts.fit_fl = false;
|
|
opts.fit_afb = false;
|
|
opts.only_angles = false;
|
|
opts.only_Bmass = false;
|
|
if(opts.only_Bmass) opts.only_angles = false;
|
|
|
|
opts.use_angular_acc = true;
|
|
opts.use_weighted_bkg = true;
|
|
opts.fit_full_angular_bkg = true;
|
|
opts.always_generate_full_angular = true;
|
|
opts.weighted_fit = true;
|
|
|
|
opts.swave = true;
|
|
opts.shift_lh = false;
|
|
|
|
opts.squared_hesse = true;
|
|
opts.minos_errors = false;
|
|
|
|
opts.generate_mkpi = true;
|
|
opts.fit_mkpi = false;
|
|
opts.use_mkpi = true;
|
|
opts.flat_bkg = false;
|
|
|
|
//change the label of the produced plots:
|
|
opts.plot_label = "LHCb toys";
|
|
|
|
//SM
|
|
std::vector<double> s1s, s1c, s2s, s2c, s3, s4, s5, s6s, s6c, s7, s8, s9;
|
|
std::vector<double> f_signal, f_bckgnd, f_bckcoeff;
|
|
double f_sig_norm = 0.0;
|
|
std::vector<UInt_t> events_per_bin;
|
|
|
|
//determine number of bins:
|
|
const UInt_t nBins = UseBinnedFit ? opts.TheQ2binsmin.size() : 1;
|
|
if(!UseBinnedFit){
|
|
opts.TheQ2binsmin = {8.68};
|
|
opts.TheQ2binsmax = {10.09};
|
|
}
|
|
|
|
//create TTree to save pull values
|
|
double treeValue, treeStart, treeError, treeErrorUp, treeErrorDown;
|
|
UInt_t treeBin, treeFitNumber, treeFitResult, treeIndex;
|
|
|
|
TTree * T = nullptr;
|
|
if(doToyPulls || doMCPulls){
|
|
T = new TTree("PullTree", "PullTree");
|
|
T->Branch("v", &treeValue, "v/D");
|
|
T->Branch("i", &treeIndex, "i/i");
|
|
T->Branch("s", &treeStart, "s/D");
|
|
T->Branch("e", &treeError, "e/D");
|
|
T->Branch("eu", &treeErrorUp, "eu/D");
|
|
T->Branch("ed", &treeErrorDown, "ed/D");
|
|
T->Branch("b", &treeBin , "b/i");
|
|
T->Branch("f", &treeFitNumber, "f/i");
|
|
T->Branch("r", &treeFitResult, "r/i");
|
|
T->Branch("fs", &SwaveFraction, "fs/D");
|
|
}
|
|
|
|
//TODO: fix the name of the genLvl root file
|
|
basic_params params_genLvl = basic_params();
|
|
params_genLvl.Run = 2;
|
|
params_genLvl.year = 2017;
|
|
params_genLvl.nBins = nBins;
|
|
std::string genPath = final_result_name_genLvlMC(params_genLvl, nBins, false, true, false);
|
|
|
|
s1s = load_param_values_into_vector("S1s",genPath);
|
|
s3 = load_param_values_into_vector("S3", genPath);
|
|
s4 = load_param_values_into_vector("S4", genPath);
|
|
s5 = load_param_values_into_vector("S5", genPath);
|
|
s6s = load_param_values_into_vector("S6s",genPath);
|
|
s7 = load_param_values_into_vector("S7", genPath);
|
|
s8 = load_param_values_into_vector("S8", genPath);
|
|
s9 = load_param_values_into_vector("S9", genPath);
|
|
|
|
if(UseBinnedFit){
|
|
if(nBins == 8){
|
|
f_signal = {101., 61., 62., 96., 125., 124., 129., 69.};
|
|
f_bckgnd = {131., 206., 279., 358., 358., 247., 158., 93.};
|
|
f_bckcoeff = {-0.0052, -0.0084, -0.0088, -0.0031, -0.0028, -0.0158, -0.0079, -0.0058}; //for notation: exp(lambda * x)
|
|
|
|
//results from Signal channel generator level MC fits:
|
|
s1c = { 0.2450, 0.6996, 0.7993, 0.7384, 0.6414, 0.4201, 0.3513, 0.3344};
|
|
s2s = { 0.1448, 0.0749, 0.0514, 0.0662, 0.0900, 0.1449, 0.1620, 0.1662};
|
|
s2c = {-0.1953,-0.6639,-0.7774,-0.7254,-0.6334,-0.4171,-0.3497,-0.3332};
|
|
s3 = { 0.0000, 0.0020,-0.0080,-0.0120,-0.0240,-0.0620,-0.1520,-0.2330};
|
|
s4 = { 0.0870, 0.0040,-0.1080,-0.1950,-0.2510,-0.2770,-0.2900,-0.3040};
|
|
s5 = { 0.2440, 0.0950,-0.1590,-0.3050,-0.4140,-0.4210,-0.3400,-0.2490};
|
|
s6s = {-0.1270,-0.2070,-0.0860, 0.0950, 0.3020, 0.5140, 0.5570, 0.4700};
|
|
s6c = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000};
|
|
s7 = {-0.0060,-0.0100,-0.0030, 0.0020, 0.0030,-0.0020,-0.0020, 0.0050};
|
|
s8 = { 0.0010,-0.0030,-0.0040,-0.0030,-0.0010, 0.0030,-0.0010, 0.0030};
|
|
s9 = {-0.0020, 0.0000, 0.0020, 0.0040,-0.0020,-0.0010,-0.0010, 0.0010};
|
|
|
|
}
|
|
else if(nBins == 4){
|
|
spdlog::error("YOU NEED TO UPDATE THESE FIRTS!");
|
|
assert(0);
|
|
f_signal = {101., 61., 62., 96.};
|
|
f_bckgnd = {131., 206., 279., 358.};
|
|
f_bckcoeff = {-0.0052, -0.0084, -0.0088, -0.0031}; //for notation: exp(lambda * x)
|
|
|
|
//results from Signal channel generator level MC fits:
|
|
s1s = { 0.5120, 0.1800, 0.1370, 0.1930};
|
|
s1c = { 0.2450, 0.6996, 0.7993, 0.7384};
|
|
s2s = { 0.1448, 0.0749, 0.0514, 0.0662};
|
|
s2c = {-0.1953,-0.6639,-0.7774,-0.7254};
|
|
s3 = { 0.0000, 0.0020,-0.0080,-0.0120};
|
|
s4 = { 0.0870, 0.0040,-0.1080,-0.1950};
|
|
s5 = { 0.2440, 0.0950,-0.1590,-0.3050};
|
|
s6s = {-0.1270,-0.2070,-0.0860, 0.0950};
|
|
s6c = { 0.0000, 0.0000, 0.0000, 0.0000};
|
|
s7 = {-0.0060,-0.0100,-0.0030, 0.0020};
|
|
s8 = { 0.0010,-0.0030,-0.0040,-0.0030};
|
|
s9 = {-0.0020, 0.0000, 0.0020, 0.0040};
|
|
}
|
|
else if(nBins == 2){
|
|
f_signal = { 61.+ 62.+ 96.+125., 129.+69.};
|
|
f_bckgnd = {206.+279.+358.+358., 158.+93.};
|
|
f_bckcoeff = {-0.0064, -0.0065};
|
|
|
|
s1s = { 0.1959, 0.4917};
|
|
s1c = { 0.7441, 0.3447};
|
|
s2s = { 0.0646, 0.1636};
|
|
s2c = {-0.7217, -0.3432};
|
|
s3 = {-0.0042, -0.1747};
|
|
s4 = { 0.1033, 0.2929};
|
|
s5 = {-0.1458, -0.3157};
|
|
s6s = { 0.0652, -0.5338};
|
|
s6c = { 0.0000, 0.0000};
|
|
s7 = { 0.0340, -0.0000};
|
|
s8 = {-0.0115, -0.0002};
|
|
s9 = { 0.0003, 0.0003};
|
|
}
|
|
else {
|
|
spdlog::error("No SM values given for binning scheme with {0:d} q2 bins. Exit", nBins);
|
|
assert(0);
|
|
}
|
|
}
|
|
else{
|
|
f_signal = {0.5};
|
|
f_bckgnd = {0.5};
|
|
f_bckcoeff = {-2.835e-3}; //for notation: exp(lambda * x)
|
|
s1s = { 0.2372};
|
|
s1c = { 0.6863};
|
|
s2s = { 0.0};
|
|
s2c = {-0.6760};
|
|
s3 = {-0.0107};
|
|
s4 = { 0.2163};
|
|
s5 = {-0.3806};
|
|
s6s = {-0.1957};
|
|
s6c = { 0.0};
|
|
s7 = { 0.0290};
|
|
s8 = {-0.0113};
|
|
s9 = { 0.0003};
|
|
}
|
|
|
|
assert(f_signal.size() == f_bckgnd.size());
|
|
f_sig_norm = std::accumulate(f_signal.begin(), f_signal.end(), 0.0)
|
|
+ std::accumulate(f_bckgnd.begin(), f_bckgnd.end(), 0.0);
|
|
|
|
for(UInt_t i = 0; i < f_signal.size(); i++){
|
|
f_signal.at(i) /= f_sig_norm;
|
|
f_bckgnd.at(i) /= f_sig_norm;
|
|
}
|
|
|
|
fcnc::folder fldr(&opts);
|
|
fcnc::fitter f(&opts);
|
|
fcnc::bu2kstarmumu_plotter thePlotter(&opts);
|
|
|
|
std::vector< std::vector<fcnc::bu2kstarmumu_parameters*> >theParameters;
|
|
std::vector< std::vector<fcnc::bu2kstarmumu_pdf*> > thePDFs;
|
|
for(UInt_t c = 0; c < nBins; c++){
|
|
//add vector of params to matrix:
|
|
std::vector<fcnc::bu2kstarmumu_parameters*> params_per_bin;
|
|
params_per_bin.clear();
|
|
theParameters.push_back(params_per_bin);
|
|
|
|
//add vector of pdfs to matrix:
|
|
std::vector<fcnc::bu2kstarmumu_pdf*> pdf_per_bin;
|
|
pdf_per_bin.clear();
|
|
thePDFs.push_back(pdf_per_bin);
|
|
}
|
|
assert(theParameters.size() == nBins);
|
|
assert(thePDFs.size() == nBins);
|
|
|
|
//pointers to current parameter, pdf and generator:
|
|
fcnc::bu2kstarmumu_parameters * params[nBins];
|
|
fcnc::bu2kstarmumu_pdf * prob[nBins];
|
|
fcnc::bu2kstarmumu_generator * gen = nullptr;
|
|
|
|
std::vector<int> fitresults;
|
|
std::vector<fcnc::event>selection;
|
|
|
|
//number of fits is either determined by the larger number of pulls or fits requested
|
|
UInt_t nFits = nPulls;
|
|
|
|
//vectors to save values for pull plots:
|
|
std::vector<int>var_indexs;
|
|
std::vector< std::vector< std::vector<double> > >pull_values;
|
|
std::vector< std::vector< std::vector<double> > >pull_errors;
|
|
std::vector< std::vector<double> >pull_starts;
|
|
|
|
|
|
//initialize parameters, events and fitter for every FitTest:
|
|
for(UInt_t n = 0; n < nFits; n++){
|
|
|
|
//if binning is used, loop over the number of q2 bins, else only run loop once!
|
|
for(UInt_t b = 0; b < nBins; b++){
|
|
|
|
std::cout << std::endl;
|
|
std::cout << "****************************************" << std::endl;
|
|
std::cout << "[FIT]\t\tStarting fit >> " << n << " <<" << std::endl;
|
|
std::cout << "[FIT]\t\tIn bin >> " << b << " <<" << std::endl;
|
|
std::cout << "****************************************" << std::endl;
|
|
std::cout << std::endl;
|
|
|
|
opts.swave = true;
|
|
|
|
opts.q2_min = UseBinnedFit ? opts.TheQ2binsmin.at(b) : opts.TheQ2binsmin.front();
|
|
opts.q2_max = UseBinnedFit ? opts.TheQ2binsmax.at(b) : opts.TheQ2binsmax.back();
|
|
|
|
if(n == 0){
|
|
//initiate parameters for mass fit
|
|
params[b] = new fcnc::bu2kstarmumu_parameters(&opts);
|
|
params[b]->f_sig.init(f_signal.at(b)/(f_signal.at(b)+f_bckgnd.at(b)), 0.0, 1.0, 0.01);
|
|
if(params[b]->f_sig.get_value() == 0.0)
|
|
params[b]->m_b.init(PDGMASS_B, B_MASS_LOW, B_MASS_HIGH, 0.0);
|
|
else
|
|
params[b]->m_b.init(PDGMASS_B, B_MASS_LOW, B_MASS_HIGH, 0.01);
|
|
if(opts.twotailedcrystalball)
|
|
params[b]->m_res_1.init(1.0, 0.0, 1.0, 0.0);
|
|
else
|
|
params[b]->m_res_1.init(0.391, 0.0, 1.0, 0.0);// 0.01);
|
|
params[b]->fm_tau.init(1.0, 0.0, 1.0, 0.0);
|
|
if(opts.fit_lambda){
|
|
if(params[b]->f_sig.get_value() == 1.0)
|
|
params[b]->m_lambda.init(f_bckcoeff.at(b), -1.0e-1, -1.0e-6, 0);
|
|
else
|
|
params[b]->m_lambda.init(f_bckcoeff.at(b), -1.0e-1, -1.0e-6, 1.0e-5);
|
|
params[b]->m_lambda_2.init(-1.664e-3, -1.0e-1, -1.0e-6, 0.0);
|
|
params[b]->m_tau.init( 1./2.835e-3, 100.0, 1.0e+4, 0.0);
|
|
params[b]->m_tau_2.init(1./1.664e-3, 0.0, 1.0e+4, 0.0);
|
|
}
|
|
else{
|
|
params[b]->m_lambda.init(-1.0e-3, -1.0e-1, -1.0e-6, 0.0);
|
|
params[b]->m_lambda_2.init(-1.5e-3, -1.0e-1, -1.0e-6, 0.0);
|
|
params[b]->m_tau.init(-1./f_bckcoeff.at(b), 100.0, 1.0e+4, 1.0);
|
|
params[b]->m_tau_2.init(6.01e+2, 0.0, 1.0e+4, 0.0);
|
|
}
|
|
if(opts.twotailedcrystalball)
|
|
params[b]->m_sigma_1.init(f_sigma_1[n], 5.0, 200.0, 0.0);//, 0.1);
|
|
else{
|
|
params[b]->m_sigma_1.init(f_sigma_1[n], 5.0, 200.0, 0.0);//, 0.1);
|
|
params[b]->m_sigma_2.init(f_sigma_2[n], 5.0, 200.0, 0.0);//, 0.1);
|
|
}
|
|
params[b]->alpha_1.init(f_alpha_1[n], 0.1, 10.0, 0.0);//, 0.15);
|
|
params[b]->alpha_2.init(f_alpha_2[n], 0.1, 10.0, 0.0);//, 0.20);
|
|
params[b]->n_1.init(f_n_1[n], 0.1, 15.0, 0.0);//, 2.0);
|
|
params[b]->n_2.init(f_n_2[n], 0.1, 10.0, 0.0);//, 0.16);
|
|
params[b]->m_scale.init(f_m_scale[n], 0.0, 2.0, 0.0);
|
|
|
|
params[b]->load_only_Bmass_param_values("fitresult_SignalFit_MC_SimultaneousFit_2Dfit_bin"+std::to_string(b)+"_pdf"+std::to_string(n)+".txt");
|
|
|
|
//S-wave
|
|
double sWaveStepSize = 0.01;
|
|
params[b]->FS.init(SwaveFraction, 0.0, 1.0, sWaveStepSize);
|
|
|
|
//B0 fit results
|
|
params[b]->SS1.init(-0.231, -1.0, 1.0, (opts.full_angular || opts.folding != 4 ? sWaveStepSize : 0.0));
|
|
params[b]->SS2.init( 0.023, -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? sWaveStepSize : 0.0));
|
|
params[b]->SS3.init( 0.003, -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? sWaveStepSize : 0.0));
|
|
params[b]->SS4.init( 0.001, -1.0, 1.0, (opts.full_angular || opts.folding > 2 ? sWaveStepSize : 0.0));
|
|
params[b]->SS5.init(-0.068, -1.0, 1.0, (opts.full_angular ? sWaveStepSize : 0.0));
|
|
|
|
if(opts.fit_mkpi || opts.use_mkpi){
|
|
params[b]->gammakstar.init(0.0503, 0.01, 0.8, 0.001); //PDG value
|
|
params[b]->mkstar.init(PDGMASS_K_STAR_PLUS / 1000., 0.7, 1.2, 0.001); //PDG value
|
|
params[b]->gammakstarplus.init(0.236, 0.1, 1.0, 0.0); //PDG value
|
|
params[b]->mkstarplus.init(1.41, 1.2, 1.6, 0.0); //PDG value
|
|
params[b]->asphase.init(TMath::Pi(), 0.0, 2.0*TMath::Pi(), 0.0);
|
|
params[b]->a.init(1.95, 0.001, 20.0, 0.0);
|
|
params[b]->r.init(1.78, 0.001, 10.0, 0.0);
|
|
//background parameter
|
|
if(params[b]->f_sig() != 1.0){
|
|
params[b]->cbkgmkpi0.init(1.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgmkpi1.init(0.0, -1.0, 1.0, 0.01);
|
|
params[b]->cbkgmkpi2.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgmkpi3.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgmkpi4.init(0.0, -1.0, 1.0, 0.0);
|
|
opts.mkpi_threshold = false;
|
|
params[b]->nthreshold.init(1.5, 0.0, 15.0, 0.0);
|
|
}
|
|
}
|
|
|
|
//define center of q2bin as effective q2:
|
|
params[b]->eff_q2.init((UseBinnedFit ? 0.5*(opts.TheQ2binsmin.at(b)+opts.TheQ2binsmax.at(b)) : 6.0), opts.TheQ2binsmin.front(), opts.TheQ2binsmax.back(), 0.0);
|
|
|
|
//if only mass fit: do not active angles in fit, i.e. set stepsize to 0.0
|
|
double angleStepSize = opts.only_Bmass ? 0.0 : 0.001;
|
|
if(params[b]->f_sig.get_value() == 0.0)
|
|
angleStepSize = 0.0;
|
|
if(opts.fit_fl)
|
|
//params[b]->Fl.init(3.0/4.0*(s1c[b]-s2c[b]/3.0), 0.0, 1.0, angleStepSize);
|
|
params[b]->Fl.init(1.0-4.0/3.0*s1s[b], 0.0, 1.0, angleStepSize);
|
|
else
|
|
params[b]->S1s.init(s1s[b], -1.0, 1.0, angleStepSize);
|
|
params[b]->S3.init(s3[b], -1.0, 1.0, angleStepSize/10.);
|
|
params[b]->S4.init(s4[b], -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? angleStepSize : 0.0));
|
|
params[b]->S5.init(s5[b], -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? angleStepSize : 0.0));
|
|
if(opts.fit_afb)
|
|
params[b]->Afb.init(3.0/4.0*s6s[b], -0.75, 0.75, (opts.full_angular || opts.folding == 0 ? angleStepSize: 0.0));
|
|
else
|
|
params[b]->S6s.init(s6s[b], -1.0, 1.0, (opts.full_angular || opts.folding == 0 ? angleStepSize : 0.0));
|
|
params[b]->S7.init(s7[b], -1.0, 1.0, (opts.full_angular || opts.folding == 3 ? angleStepSize/10. : 0.0));
|
|
params[b]->S8.init(s8[b], -1.0, 1.0, (opts.full_angular || opts.folding == 4 ? angleStepSize/10. : 0.0));
|
|
params[b]->S9.init(s9[b], -1.0, 1.0, (opts.full_angular || opts.folding == 0 ? angleStepSize/100. : 0.0));
|
|
|
|
if(!opts.flat_bkg){
|
|
//ctl
|
|
params[b]->cbkgctl0.init(1.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgctl2.init(0.0, -1.0, 1.0, 0.01);
|
|
params[b]->cbkgctl4.init(0.0, -1.0, 1.0, 0.0);
|
|
if(opts.full_angular || opts.folding == 0){
|
|
params[b]->cbkgctl1.init(0.0, -1.0, 1.0, 0.01);
|
|
params[b]->cbkgctl3.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
else{
|
|
params[b]->cbkgctl1.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgctl3.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
|
|
//ctk
|
|
params[b]->cbkgctk0.init(1.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgctk2.init(0.05, -1.0, 1.0, 0.01);
|
|
params[b]->cbkgctk4.init(0.0, -1.0, 1.0, 0.0);
|
|
if(opts.full_angular || opts.folding != 4){
|
|
params[b]->cbkgctk1.init(0.1, -1.0, 1.0, 0.01);
|
|
params[b]->cbkgctk3.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
else{
|
|
params[b]->cbkgctk1.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgctk3.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
|
|
//phi
|
|
params[b]->cbkgphi0.init(1.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgphi2.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgphi4.init(0.0, -1.0, 1.0, 0.0);
|
|
if(opts.full_angular || opts.folding == 0 || opts.folding == 3 || opts.folding == 4){
|
|
params[b]->cbkgphi1.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgphi3.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
else{
|
|
params[b]->cbkgphi1.init(0.0, -1.0, 1.0, 0.0);
|
|
params[b]->cbkgphi3.init(0.0, -1.0, 1.0, 0.0);
|
|
}
|
|
|
|
params[b]->load_only_bckgnd_param_values(("fitresult_OnlyBckgnd_bin"+std::to_string(b)+".txt").c_str());
|
|
|
|
}
|
|
}
|
|
else{
|
|
//S-wave
|
|
double sWaveStepSize = 0.01;
|
|
params[b]->FS.init(SwaveFraction, 0.0, 1.0, sWaveStepSize);
|
|
/*
|
|
//JPSI Fit results
|
|
params[b]->SS1.init( 0.549, -1.0, 1.0, (opts.full_angular || opts.folding != 4 ? sWaveStepSize : 0.0));
|
|
params[b]->SS2.init(-0.128, -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? sWaveStepSize : 0.0));
|
|
params[b]->SS3.init( 0.003, -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? sWaveStepSize : 0.0));
|
|
params[b]->SS4.init(-0.006, -1.0, 1.0, (opts.full_angular || opts.folding > 2 ? sWaveStepSize : 0.0));
|
|
params[b]->SS5.init(-0.149, -1.0, 1.0, (opts.full_angular ? sWaveStepSize : 0.0));
|
|
*/
|
|
|
|
//B0 fit results
|
|
params[b]->SS1.init(-0.231, -1.0, 1.0, (opts.full_angular || opts.folding != 4 ? sWaveStepSize : 0.0));
|
|
params[b]->SS2.init( 0.023, -1.0, 1.0, (opts.full_angular || opts.folding == 1 ? sWaveStepSize : 0.0));
|
|
params[b]->SS3.init( 0.003, -1.0, 1.0, (opts.full_angular || opts.folding == 2 ? sWaveStepSize : 0.0));
|
|
params[b]->SS4.init( 0.001, -1.0, 1.0, (opts.full_angular || opts.folding > 2 ? sWaveStepSize : 0.0));
|
|
params[b]->SS5.init(-0.068, -1.0, 1.0, (opts.full_angular ? sWaveStepSize : 0.0));
|
|
params[b]->reset_parameters();
|
|
|
|
}
|
|
|
|
//create vectors to save parameter values for pull plots:
|
|
if(n == 0){
|
|
if(doToyPulls || doMCPulls){
|
|
UInt_t pp = 0;
|
|
for(UInt_t p = 0; p < params[b]->nparameters(); p++){
|
|
if(params[b]->get_parameter(p)->get_step_size() != 0.0){
|
|
if(b == 0){
|
|
std::vector< std::vector<double> >pulls_per_parameter;
|
|
pull_values.push_back(pulls_per_parameter);
|
|
|
|
std::vector< std::vector<double> >errors_per_parameter;
|
|
pull_errors.push_back(errors_per_parameter);
|
|
|
|
std::vector<double>starts_per_parameter;
|
|
pull_starts.push_back(starts_per_parameter);
|
|
|
|
var_indexs.push_back(p);
|
|
}
|
|
//check that vectors have same size and that no alloc error occurs:
|
|
assert(pp < pull_values.size());
|
|
assert(var_indexs.size() == pull_values.size());
|
|
assert(pp < pull_errors.size());
|
|
assert(var_indexs.size() == pull_errors.size());
|
|
assert(pp < pull_starts.size());
|
|
assert(var_indexs.size() == pull_starts.size());
|
|
|
|
//save vectors for values and start values to matrices
|
|
std::vector<double>pulls_per_param_per_bin;
|
|
pull_values.at(pp).push_back(pulls_per_param_per_bin);
|
|
std::vector<double>errors_per_param_per_bin;
|
|
pull_errors.at(pp).push_back(errors_per_param_per_bin);
|
|
pull_starts.at(pp).push_back(params[b]->get_parameter(p)->get_start_value());
|
|
|
|
//increase parameter index:
|
|
pp++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//pdf
|
|
prob[b] = new fcnc::bu2kstarmumu_pdf(&opts, params[b]);
|
|
if(opts.use_angular_acc || opts.weighted_fit){
|
|
prob[b]->load_coeffs_eff_phsp_4d();
|
|
}
|
|
prob[b]->update_cached_normalization(params[b]);
|
|
|
|
//create vector with events according to the requested fits/pulls
|
|
selection.clear();
|
|
|
|
Double_t nEvents = 0.;
|
|
|
|
if(UseBinnedFit){
|
|
nEvents *= f_signal.at(b)+f_bckgnd.at(b);
|
|
}
|
|
else
|
|
nEvents = nToyEvents;
|
|
|
|
std::cout << "Generate " << nEvents << " events!" << std::endl;
|
|
gen = new fcnc::bu2kstarmumu_generator(&opts);
|
|
selection = gen->generate(nEvents, params[b], prob[b]);
|
|
//delete gen;
|
|
if(!opts.full_angular && opts.always_generate_full_angular){
|
|
for(UInt_t ee = 0; ee < selection.size(); ee++){
|
|
fldr.fold(&selection.at(ee));
|
|
}
|
|
}
|
|
|
|
//check if all angles are almost 0.0
|
|
if(false){
|
|
for(UInt_t i = 0; i < selection.size(); i++)
|
|
if((fabs(selection.at(i).costhetak) < 0.001) && (fabs(selection.at(i).costhetal) < 0.001) && (fabs(selection.at(i).phi) < 0.001))
|
|
std::cout << "[WARNING]\tEvent" << i << "\tctk: " << selection.at(i).costhetak << "\tctl: " << selection.at(i).costhetal << "\tphi: " << selection.at(i).phi << std::endl;
|
|
}
|
|
|
|
if(n == 0)events_per_bin.push_back(selection.size());
|
|
|
|
//deactive the S-wave for the fit!
|
|
/*
|
|
opts.swave = false;
|
|
params[b]->FS.init(0.0, 0.0, 1.0, 0.0);
|
|
params[b]->SS1.init( 0.0, -1.0, 1.0, 0.0);
|
|
params[b]->SS2.init( 0.0, -1.0, 1.0, 0.0);
|
|
params[b]->SS3.init( 0.0, -1.0, 1.0, 0.0);
|
|
params[b]->SS4.init( 0.0, -1.0, 1.0, 0.0);
|
|
params[b]->SS5.init( 0.0, -1.0, 1.0, 0.0);
|
|
*/
|
|
//set values to jpsi results for the fit!
|
|
params[b]->FS.init(0.0, 0.0, 1.0, 0.0);
|
|
params[b]->SS1.init( 0.549, -1.0, 1.0, 0.0);
|
|
params[b]->SS2.init(-0.128, -1.0, 1.0, 0.0);
|
|
params[b]->SS3.init( 0.003, -1.0, 1.0, 0.0);
|
|
params[b]->SS4.init(-0.006, -1.0, 1.0, 0.0);
|
|
params[b]->SS5.init(-0.149, -1.0, 1.0, 0.0);
|
|
|
|
//make sure that the correct coefficients are chosen for the fit:
|
|
if(!opts.full_angular && (opts.use_angular_acc || opts.weighted_fit) && opts.always_generate_full_angular){
|
|
prob[b]->load_coeffs_eff_phsp_4d();
|
|
prob[b]->update_cached_normalization(params[b]);
|
|
}
|
|
|
|
//fit the events:
|
|
bool do_fit = true;
|
|
int fitresult = 0;
|
|
if(do_fit){
|
|
fitresult = f.fit(prob[b], params[b], &selection);
|
|
}
|
|
else{ //don't fit, just update the efficiencies
|
|
fitresult = 300;
|
|
if(opts.use_angular_acc || opts.weighted_fit){
|
|
prob[b]->update_cached_efficiencies(params[b], &selection);
|
|
}
|
|
}
|
|
|
|
fitresults.push_back(fitresult);
|
|
//nametag per bin!
|
|
std::string plotname;
|
|
plotname = "ToyGen_"+std::to_string(nToyEvents)+"ToyEvents_FS_"+std::to_string(SwaveFraction)+"_";
|
|
|
|
plotname.append(angularsuffix);
|
|
if(UseBinnedFit){
|
|
plotname.append("_bin"+std::to_string(b));
|
|
if(Fit1bin)
|
|
plotname.append("_1BIN");
|
|
if(Fit2bins)
|
|
plotname.append("_2BINS");
|
|
if(FitAllbins)
|
|
plotname.append("_9BINS");
|
|
}
|
|
|
|
//plot the current pdf with event data:
|
|
if(n == 0 && opts.write_eps){
|
|
std::cout << "PLOT: " << selection.size() << " events" << std::endl;
|
|
thePlotter.plot_data(prob[b], params[b], &selection, get_PullPlot_path(), plotname+"_Fit"+std::to_string(n), true);
|
|
|
|
//add pdfs for all binning:
|
|
thePDFs.at(b).push_back(prob[b]);
|
|
theParameters.at(b).push_back(params[b]);
|
|
|
|
}
|
|
|
|
//save param values to pull_values vector:
|
|
if(doToyPulls){
|
|
for(UInt_t p = 0; p < var_indexs.size(); p++){
|
|
int index = var_indexs.at(p);
|
|
|
|
double val = params[b]->get_parameter(index)->get_value();
|
|
double err = params[b]->get_parameter(index)->get_error();
|
|
|
|
pull_values.at(p).at(b).push_back(val);
|
|
pull_errors.at(p).at(b).push_back(err);
|
|
|
|
//save values to tree:
|
|
treeValue = params[b]->get_parameter(index)->get_value();
|
|
treeStart = params[b]->get_parameter(index)->get_start_value();
|
|
treeError = params[b]->get_parameter(index)->get_error();
|
|
treeErrorUp = params[b]->get_parameter(index)->get_error_up();
|
|
treeErrorDown = params[b]->get_parameter(index)->get_error_down();
|
|
|
|
treeIndex = index;
|
|
treeBin = b;
|
|
treeFitNumber = n;
|
|
treeFitResult = fitresult;
|
|
|
|
//treeVarName = params[b]->get_parameter(index)->get_name();
|
|
//treeVarDesc = params[b]->get_parameter(index)->get_description();
|
|
|
|
T->Fill();
|
|
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//save TTree to .root file:
|
|
if(doToyPulls){
|
|
TFile * F = new TFile(("PullResults_"
|
|
+std::to_string(job_id == -10 ? nMCEvents : nToyEvents)
|
|
+(job_id == -10 ? "_MC_" : "_Toys_")
|
|
+std::to_string(nPulls)+"_Fits_"
|
|
+angularsuffix
|
|
+"_F_S_"+std::to_string(SwaveFraction)
|
|
+(job_id >= 0 ? "_job"+std::to_string(job_id) : "")
|
|
+".root").c_str(), "RECREATE");
|
|
F->cd();
|
|
T->Write();
|
|
F->Close();
|
|
delete T;
|
|
}
|
|
}
|