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