kopecna/EWP-BplusToKstMuMu-AngAna
kopecna/EWP-BplusToKstMuMu-AngAna
2
1
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
8d454c69bd
EWP-BplusToKstMuMu-AngAna/Code/FCNCFitter/sources/Core/bu2kstarmumu_pdf.cc

5381 lines
251 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

//Renata Kopecna
#include <fstream>
#include "folder.hh"
#include <funcs.hh>
#include <constants.hh>
#include <bu2kstarmumu_pdf.hh>
#include <pdf.hh>
#include <options.hh>
#include <values.hh>
#include <design.hh>
#include <integrals.hh>
#include <event.hh>
#include <helpers.hh>
#include <paths.hh>
#include <TMinuit.h>
#include <spdlog.h>
#include <TROOT.h>
#include <TFile.h>
#include <TStyle.h>
//Helper function to create a polynomial for the background
std::vector<double> init_ch_ctl(const fcnc::bu2kstarmumu_parameters* params){
std::vector<double> bkg_ctl;
bkg_ctl.push_back(params->cbkgctl0());
bkg_ctl.push_back(params->cbkgctl1());
bkg_ctl.push_back(params->cbkgctl2());
bkg_ctl.push_back(params->cbkgctl3());
bkg_ctl.push_back(params->cbkgctl4());
return bkg_ctl;
}
std::vector<double> init_ch_ctk(const fcnc::bu2kstarmumu_parameters* params){
std::vector<double> bkg_ctk;
bkg_ctk.push_back(params->cbkgctk0());
bkg_ctk.push_back(params->cbkgctk1());
bkg_ctk.push_back(params->cbkgctk2());
bkg_ctk.push_back(params->cbkgctk3());
bkg_ctk.push_back(params->cbkgctk4());
bkg_ctk.push_back(params->cbkgctk5());
bkg_ctk.push_back(params->cbkgctk6());
return bkg_ctk;
}
std::vector<double> init_ch_phi(const fcnc::bu2kstarmumu_parameters* params){
std::vector<double> bkg_phi;
bkg_phi.push_back(params->cbkgphi0());
bkg_phi.push_back(params->cbkgphi1());
bkg_phi.push_back(params->cbkgphi2());
bkg_phi.push_back(params->cbkgphi3());
bkg_phi.push_back(params->cbkgphi4());
return bkg_phi;
}
//using namespace std;
using namespace fcnc;
//this is a wrapper for the complex error function
//which is needed for a gaussian convolution of exp*cos/sin
//this fortran implementation however is not thread safe!
extern "C" void wwerfwrap_(double *zr, double *zi, double *wr, double* wi);
//boost::mutex global_mutex_2;
std::mutex global_mutex_2;
bu2kstarmumu_pdf* bu2kstarmumu_pdf::current_pdf;
bu2kstarmumu_pdf::bu2kstarmumu_pdf(options* o, const bu2kstarmumu_parameters* params):
opts(o)
{
//vector of efficiency coefficients containing 1260 elements, coefficients stored in bu2kstarmumu_coefficients.cc
//coeffs_eff_4d = fcnc::coefficients_nominal;
//load_coeffs_eff_phsp_4d();
fldr = new fcnc::folder(opts);
};
bu2kstarmumu_pdf::~bu2kstarmumu_pdf()
{
};
void bu2kstarmumu_pdf::get_moments(std::vector<fcnc::event> events, std::vector<double>& moments, std::vector<double>& covariance) const
{
if(!opts->full_angular){
spdlog::error("No folding implemented for moments. Debug: 1");
assert(0);
}
std::vector<double> cov(NTERMS*NTERMS, 0.0);//the actual covariance matrix for the moments
std::vector<double> m(NTERMS, 0.0);//need to first calculate the mean
double nweight = 0.0; //sum of weights
double nweightsquared = 0.0; //sum of weights squared, not used anywhere
for (auto meas: events){
double weight = meas.weight*meas.sweight;
nweight += weight;
nweightsquared += weight*weight;
}
for (auto meas: events){
double weight = meas.weight*meas.sweight;
weight /= nweight;//nominal
//Calculate the angular terms for each parameter
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
fj(meas.costhetal, meas.costhetak, meas.phi,
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
double ffs=0.0, fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0;
if(opts->swave){
swave_fj(meas.costhetal, meas.costhetak, meas.phi,
ffs, fs1, fs2, fs3, fs4, fs5);
}
//TODO when bored: this is very annoying and f_whatever should be a struct or so and one should loop
m.at(0) += f1*weight;
m.at(1) += f2*weight;
m.at(2) += f3*weight;
m.at(3) += f4*weight;
m.at(4) += f5*weight;
m.at(5) += f6*weight;
m.at(6) += f7*weight;
m.at(7) += f8*weight;
m.at(8) += f9*weight;
m.at(9) += f10*weight;
m.at(10) += f11*weight;
m.at(11) += f12*weight;
if(opts->swave){
m.at(12) += ffs*weight;
m.at(13) += fs1*weight;
m.at(14) += fs2*weight;
m.at(15) += fs3*weight;
m.at(16) += fs4*weight;
m.at(17) += fs5*weight;
}
}
//calculate rms
for (auto meas: events){
double weight = meas.weight*meas.sweight;//dangerous?
weight /= nweight;//nominal
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
fj(meas.costhetal, meas.costhetak, meas.phi,
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
double ffs=0.0, fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0;
if(opts->swave){
swave_fj(meas.costhetal, meas.costhetak, meas.phi,
ffs, fs1, fs2, fs3, fs4, fs5);
}
std::vector<double> diff(NTERMS*NTERMS, 0.0);
diff.at(0) = f1 - m.at(0); //aaaaaaaarrrrggghhh
diff.at(1) = f2 - m.at(1);
diff.at(2) = f3 - m.at(2);
diff.at(3) = f4 - m.at(3);
diff.at(4) = f5 - m.at(4);
diff.at(5) = f6 - m.at(5);
diff.at(6) = f7 - m.at(6);
diff.at(7) = f8 - m.at(7);
diff.at(8) = f9 - m.at(8);
diff.at(9) = f10 - m.at(9);
diff.at(10) = f11 - m.at(10);
diff.at(11) = f12 - m.at(11);
if(opts->swave){
diff.at(12) = ffs - m.at(12);
diff.at(13) = fs1 - m.at(13);
diff.at(14) = fs2 - m.at(14);
diff.at(15) = fs3 - m.at(15);
diff.at(16) = fs4 - m.at(16);
diff.at(17) = fs5 - m.at(17);
}
//for swave==false, only loop over the first 12x12 elements
//the matrix however stays in the 18x18 format, with 0.0 for all j > 11 || k > 11
UInt_t mm = opts->swave ? NTERMS : NTERMS_P;
for (unsigned int j = 0; j<mm; j++){
for (unsigned int k = 0; k<mm; k++){
cov.at(j*NTERMS+k) += (diff.at(j))*(diff.at(k))*sqr(weight);
}
}
}
moments = m;
covariance = cov;
return;
};
double corr_from_obs(std::vector<double> obs, std::vector<double> obscov, int i, int k){
return obscov.at(i+NTERMS*k)/sqrt(obscov.at(k+NTERMS*k))/sqrt(obscov.at(i+NTERMS*i));
}
void bu2kstarmumu_pdf::save_moments_to_obs(bu2kstarmumu_parameters* params,
std::vector<double> obs,
std::vector<double> obscov) const
{
assert(obs.size() == NTERMS);
assert(obscov.size() == NTERMS*NTERMS);
//**** Fl ****//
params->Fl.set_value(1.0-4.0/3.0*obs.at(0));
params->Fl.set_error(4.0/3.0*sqrt(obscov.at(0+NTERMS*0)));
if(opts->fit_fl) params->Fl.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++){
if(i == 0) params->Fl.set_correlation(1.0, 0); //self correlation is one by definition
else if(i == 7 && opts->fit_afb){
//correlation to Afb is equal to correlation of S1s and S6s
params->Fl.set_correlation(corr_from_obs(obs,obscov,7,0), 7);
}
else{
params->Fl.set_correlation(corr_from_obs(obs,obscov,i,0)/sqrt(3.0/4.0), i);
}
}
//**** S1s ****//
params->S1s.set_value(obs.at(0));
params->S1s.set_error(sqrt(obscov.at(0+NTERMS*0)));
if(!opts->fit_fl) params->S1s.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++){
params->S1s.set_correlation(corr_from_obs(obs,obscov,i,0), i);
}
//**** S3 ****//
params->S3.set_value(obs.at(4));
params->S3.set_error(sqrt(obscov.at(4+NTERMS*4)));
params->S3.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++){
params->S3.set_correlation(corr_from_obs(obs,obscov,i,4), i);
}
//**** S4 ****//
params->S4.set_value(obs.at(5));
params->S4.set_error(sqrt(obscov.at(5+NTERMS*5)));
params->S4.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++)
params->S4.set_correlation(corr_from_obs(obs,obscov,i,5), i);
//**** S5 ****//
params->S5.set_value(obs.at(6));
params->S5.set_error(sqrt(obscov.at(6+NTERMS*6)));
params->S5.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++)
params->S5.set_correlation(corr_from_obs(obs,obscov,i,6), i);
//**** Afb ****//
params->Afb.set_value(3.0/4.0*obs.at(7));
params->Afb.set_error(3.0/4.0*sqrt(obscov.at(7+NTERMS*7)));
if(opts->fit_afb)
params->Afb.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++){
if(i == 7) params->Afb.set_correlation(1.0, 7);//self correlation is not scaled
else if(i == 0 && opts->fit_afb){
//correlation to Fl is equal to correlation of S1s and S6s
params->Afb.set_correlation(corr_from_obs(obs,obscov,0,7), 0);
}
else params->Afb.set_correlation(corr_from_obs(obs,obscov,i,7)/sqrt(4.0/3.0), i);
}
//**** S6s ****//
params->S6s.set_value(obs.at(7));
params->S6s.set_error(sqrt(obscov.at(7+NTERMS*7)));
if(!opts->fit_afb) params->S6s.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->S6s.set_correlation(corr_from_obs(obs,obscov,i,7), i);
//**** S7 ****//
params->S7.set_value(obs.at(9));
params->S7.set_error(sqrt(obscov.at(9+NTERMS*9)));
params->S7.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->S7.set_correlation(corr_from_obs(obs,obscov,i,9), i);
//**** S8 ****//
params->S8.set_value(obs.at(10));
params->S8.set_error(sqrt(obscov.at(10+NTERMS*10)));
params->S8.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->S8.set_correlation(corr_from_obs(obs,obscov,i,10), i);
//**** S9 ****//
params->S9.set_value(obs.at(11));
params->S9.set_error(sqrt(obscov.at(11+NTERMS*11)));
params->S9.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->S9.set_correlation(corr_from_obs(obs,obscov,i,11), i);
//**** FS ****//
params->FS.set_value(obs.at(12));
params->FS.set_error(sqrt(obscov.at(12+NTERMS*12)));
if(opts->swave) params->FS.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->FS.set_correlation(corr_from_obs(obs,obscov,i,12), i);
//**** SS1 ****//
params->SS1.set_value(obs.at(13));
params->SS1.set_error(sqrt(obscov.at(13+NTERMS*13)));
if(opts->swave) params->SS1.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->SS1.set_correlation(corr_from_obs(obs,obscov,i,13), i);
//**** SS2 ****//
params->SS2.set_value(obs.at(14));
params->SS2.set_error(sqrt(obscov.at(14+NTERMS*14)));
if(opts->swave) params->SS2.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->SS2.set_correlation(corr_from_obs(obs,obscov,i,14), i);
//**** SS3 ****//
params->SS3.set_value(obs.at(15));
params->SS3.set_error(sqrt(obscov.at(15+NTERMS*15)));
if(opts->swave) params->SS3.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->SS3.set_correlation(corr_from_obs(obs,obscov,i,15), i);
//**** SS4 ****//
params->SS4.set_value(obs.at(16));
params->SS4.set_error(sqrt(obscov.at(16+NTERMS*16)));
if(opts->swave)params->SS4.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->SS4.set_correlation(corr_from_obs(obs,obscov,i,16), i);
//**** SS5 ****//
params->SS5.set_value(obs.at(17));
params->SS5.set_error(sqrt(obscov.at(17+NTERMS*17)));
if(opts->swave) params->SS5.set_step_size(0.0001);
for(unsigned int i = 0; i < NTERMS; i++) params->SS5.set_correlation(corr_from_obs(obs,obscov,i,17), i);
}
double bu2kstarmumu_pdf::prob(const bu2kstarmumu_parameters* params, const event& meas) const
{
//safety checks:
if(opts->fit_mkpi && opts->use_mkpi){
spdlog::critical("Do not use the direct fit of m(Kpi) and the included fit of m(Kpi) in the angular s-Wave at the same time!");
spdlog::critical("5D fit:\t'use_mkpi'");
spdlog::critical("4D+1D fit:\t'fit_mkpi'");
assert(0);
}
if(opts->only_mass2DFit && opts->only_Bmass){
spdlog::critical("Cannot fit m(Kpi) and m(B) with option 'only_Bmass'!");
assert(0);
}
if(opts->only_mass2DFit && !opts->fit_mkpi){
spdlog::critical("When using the 2D mass fit, enable 'fit_mkpi' as well!");
assert(0);
}
if(opts->only_mkpi && !opts->fit_mkpi){
spdlog::critical("When using the m(Kpi) mass fit, enable 'fit_mkpi' as well!");
assert(0);
}
if((opts->only_angles && opts->only_Bmass) || (opts->only_mkpi && opts->only_Bmass) || (opts->only_angles && opts->only_mkpi)){
spdlog::critical("Do not use two or more of 'only_Bmass', 'only_mkpi' and 'only_angles' at the same time!");
spdlog::critical("only_Bmass = " + boolToString(opts->only_Bmass));
spdlog::critical("only_mkpi = " + boolToString(opts->only_mkpi));
spdlog::critical("only_angles = " + boolToString(opts->only_angles));
assert(0);
}
if(opts->use_angular_acc && opts->weighted_fit){
spdlog::critical("Do not use 'use_angular_acc' and 'weighted_fit' at the same time!");
spdlog::critical("'use_angular_acc'\tConvolute acceptance function into fit pdf");
spdlog::critical("'weighted_fit'\tAssign per event weight according to acceptance");
assert(0);
}
if(opts->twotailedcrystalball && opts->crystalball){
spdlog::critical("Cannot use CB and two-tailed CB option. Choose only one and chose wisely:");
spdlog::critical("Either 'twotailedcrystalball' or 'crystalball'");
assert(0);
}
if(opts->full_angular && opts->folding > -1){
spdlog::critical("Cannot run full angular fit and folding at the same time");
assert(0);
}
//determine signal and background probability
double fsig = (opts->extended_ml ? params->n_sig()/(params->n_sig()+params->n_bkg()) : params->f_sig());
double sig_part = 0.0;
double bkg_part = 0.0;
if(fsig != 0.0)//signal part:
sig_part += fsig
* angular_sig_prob(params, meas)
* m_sig_prob(params, meas)
* mkpi_sig_prob(params, meas);
if(fsig != 1.0)//background part:
bkg_part += (1.0-fsig)
* angular_bkg_prob(params, meas)
* m_bkg_prob(params, meas)
* mkpi_bkg_prob(params, meas);
//force individual probabilities for signal and background to be positive
if(opts->individual_penalties){
if(sig_part < 0.0){
sig_part = 1.0e-12;
bkg_part = 0.0;
}
if(bkg_part < 0.0){
bkg_part = 1.0e-12;
sig_part = 0.0;
}
}
//add background and signal parts:
double result = sig_part + bkg_part;
//check final probability
if(result < 0.0){
std::lock_guard<std::mutex> lock(global_mutex_2);
spdlog::debug("Event has negative probability {0:e}", result);
spdlog::debug("Signal_part:\t{0:e}", sig_part);
spdlog::debug("m(B):\t\t{0:e}",m_sig_prob(params, meas)); //THIS one crashes and idk why
spdlog::debug("m(Kpi):\t{0:e}",mkpi_sig_prob(params, meas));
spdlog::debug("angles:\t{0:e}",angular_sig_prob(params, meas));
spdlog::debug("FS:\t\t{0:e}", params->FS.get_value());
spdlog::debug("Backgrnd_part:\t{0:e}", bkg_part);
if (fsig != 1.0){
spdlog::debug("m(B):\t{0:e}",m_bkg_prob(params, meas));
spdlog::debug("m(Kpi):\t{0:e}",mkpi_bkg_prob(params, meas));
spdlog::debug("angles:\t{0:e}",angular_bkg_prob(params, meas));
}
const bool variablepushback = false; //OMFG
if(variablepushback){
double baseline = 1.0e-12;
double absres = fabs(result) > 1.0e-6 ? fabs(result) : 1.0e-6;
double logres = log10(absres)+6.0;
double factor = pow(10.0, -logres);
result = baseline*factor;
}
else{
result = 1.0e-12;
}
}
else if(std::isinf(result) || std::isnan(result)){
std::lock_guard<std::mutex> lock(global_mutex_2);
params->print_parameters(false);
spdlog::error("Event has probability {0:f}", result);
spdlog::error("Signal_part:\t{0:f}", sig_part);
spdlog::error("fsig:\t\t{0:f}", fsig);
spdlog::error("m(B):\t\t{0:f}", m_sig_prob(params, meas));
spdlog::error("m(Kpi):\t{0:e}", mkpi_sig_prob(params, meas));
spdlog::error("angles:\t{0:e}",angular_sig_prob(params, meas));
spdlog::error("FS:\t\t{0:f}", params->FS.get_value());
spdlog::error("m(K*+) s-Wave NORM:\t{0:f}", mkpi_swave_2_norm);
spdlog::error("m(K*+) p-Wave NORM:\t{0:f}", mkpi_pwave_2_norm);
spdlog::error("Backgrnd_part:\t{0:f}", bkg_part);
assert(0);
}
return result;
};
void bu2kstarmumu_pdf::init(parameters* params) {
init(dynamic_cast<bu2kstarmumu_parameters*>(params));
//update_cached_normalization(params);
};
struct bin_edges{ //This could've been a class in case someone is bored in the future
//limits for all variable (range) (4D):
double q2min, q2max;
double ctkmin, ctkmax;
double ctlmin, ctlmax;
double phimin, phimax;
unsigned int nq2bins, nctlbins, nctkbins, nphibins;
//keep this as functions in case of modyfing the ranges during running the code!
//range of each dimension (4D):
double range_q2(){
return q2max - q2min;
}
double range_ctk(){
return ctkmax - ctkmin;
}
double range_ctl(){
return ctlmax - ctlmin;
}
double range_phi(){
return phimax - phimin;
}
//bin width for each variable
double binWidth_q2(){
return range_q2()/nq2bins;
}
double binWidth_ctl(){
return range_ctl()/nctlbins;
}
double binWidth_ctk(){
return range_ctk()/nctkbins;
}
double binWidth_phi(){
return range_phi()/nphibins;
}
//bin edges for each variable
std::vector<double> edges_q2(){
std::vector<double> tmp;
for (unsigned int h = 0; h<=nq2bins; h++){
tmp.push_back(q2min + h * binWidth_q2());
}
return tmp;
}
std::vector<double> edges_ctl(){
std::vector<double> tmp;
for (unsigned int h = 0; h<=nctlbins; h++){
tmp.push_back(ctlmin + h * binWidth_ctl());
}
return tmp;
}
std::vector<double> edges_ctk(){
std::vector<double> tmp;
for (unsigned int h = 0; h<=nctkbins; h++){
tmp.push_back(ctkmin + h * binWidth_ctk());
}
return tmp;
}
std::vector<double> edges_phi(){
std::vector<double> tmp;
for (unsigned int h = 0; h<=nphibins; h++){
tmp.push_back(phimin + h * binWidth_phi());
}
return tmp;
}
//Total number of bins in 4D
unsigned int getNbins(){
return nq2bins*nctlbins*nctkbins*nphibins;
}
//constructor
bin_edges(fcnc::options *opts, unsigned int nq2,unsigned int nctl,unsigned int nctk,unsigned int nphi){
//get the limits for all variables (4D):
q2min = opts->TheQ2binsmin.front(); q2max = opts->TheQ2binsmax.back();
ctkmin = opts->ctk_min; ctkmax = opts->ctk_max;
ctlmin = opts->ctl_min; ctlmax = opts->ctl_max;
phimin = opts->phi_min; phimax = opts->phi_max;
nq2bins = nq2;
nctlbins = nctl;
nctkbins = nctk;
nphibins = nphi;
}
};
unsigned int get_bin_in4D(int h, int i, int j, int k, bin_edges edges){
/// h = q2, i = ctl, j = ctk, k = phi
//Return the bin number for given event evt in 4D (q2, thetal, thetak, phi)
//Starting from 0!
unsigned int bin = edges.nctlbins*edges.nctkbins*edges.nphibins*h //q2
+edges.nctkbins*edges.nphibins*i //thetaL
+edges.nphibins*j //thetaK
+k; //phi
if (bin >= edges.getNbins()){
spdlog::critical("Bin {0:d} is larger than the dimension of {1:d}!",bin,edges.getNbins());
spdlog::error("h {0:d},\ti {1:d},\tj {2:d},\tk {3:d}", h, i, j, k);
spdlog::error("q2 {0:d},\tctl {1:d},\tctk {2:d},\tphi {3:d}",
edges.nq2bins, edges.nctlbins, edges.nctkbins, edges.nphibins);
assert(0);
}
else return bin;
}
unsigned int get_bin_in4D(event evt, bin_edges edges ){
//Return the bin number for given event evt in 4D (q2, thetal, thetak, phi)
//Starting from 0!
unsigned int h = (evt.q2 - edges.q2min ) / edges.binWidth_q2(); //division here truncates the result"chi
unsigned int i = (evt.costhetal - edges.ctlmin) / edges.binWidth_ctl();//eg 18/10 = 1
unsigned int j = (evt.costhetak - edges.ctkmin) / edges.binWidth_ctk();
unsigned int k = (evt.phi - edges.phimin) / edges.binWidth_phi();
return get_bin_in4D(h,i,j,k, edges);
}
unsigned int npolynom::get_bin_in4D(int h, int i, int j, int k){
/// h = q2, i = ctl, j = ctk, k = phi
//Return the bin number for given event evt in 4D (q2, thetal, thetak, phi)
//Starting from 0!
return ctl*ctk*phi*h //q2
+ctk*phi*i //thetaL
+phi*j //thetaK
+k; //phi
}
unsigned int get_bin_in3D(int i, int j, int k, npolynom npol){
/// h = q2, i = ctl, j = ctk, k = phi
//Return the bin number for given event evt in 4D (q2, thetal, thetak, phi)
//Starting from 0!
return npol.ctk*npol.phi*i //thetaL
+npol.phi*j //thetaK
+k; //phi
}
unsigned int get_bin_in3D(int i, int j, int k, bin_edges edges){
/// h = q2, i = ctl, j = ctk, k = phi
//Return the bin number for given event evt in 4D (q2, thetal, thetak, phi)
//Starting from 0!
return edges.nctkbins*edges.nphibins*i //thetaL
+edges.nphibins*j //thetaK
+k; //phi
}
void get_chi2(bin_edges edges, npolynom npol,
std::vector<event> eventsInQ2bin, std::vector<double> coeffs_moments,
values* globalvalues,
bool assumectleven, bool assumePhiEven, bool write_eps,
std::string subfolder, std::string appendix){
spdlog::info( "Starting chi2 calculation" );
std::vector<double> nsel_hist(edges.getNbins(), 0.0);
std::vector<double> nsel_hist_errsq(edges.getNbins(), 0.0);
//this is actually not the correct error here, due to the weights!
//should be sqrt(sum_i weights_i^2)
double nsel = 0.0;
int nEvents = eventsInQ2bin.size();
for (auto &meas: eventsInQ2bin) {
unsigned int bin = get_bin_in4D(meas,edges);
nsel_hist.at(bin) += meas.weight;
nsel += meas.weight;
nsel_hist_errsq.at(bin) += sqr(meas.weight);
}
//Set errors
for (auto err: nsel_hist_errsq){
if (err == 0.0) err = 1.0;
}
//compare with prediction
double neps = 0.0;
std::vector<double> neps_hist(edges.getNbins(), 0.0);
std::vector<double> q2_edges = edges.edges_q2();
std::vector<double> ctl_edges = edges.edges_ctl();
std::vector<double> ctk_edges = edges.edges_ctk();
std::vector<double> phi_edges = edges.edges_phi();
for (unsigned int h = 0; h<edges.nq2bins; h++){
for (unsigned int i = 0; i<edges.nctlbins; i++){
for (unsigned int j = 0; j<edges.nctkbins; j++){
for (unsigned int k = 0; k<edges.nphibins; k++){
double integral = 0.0;
for (unsigned int l = 0; l<npol.q2; l++){
for (unsigned int m = 0; m<npol.ctl; m++){
for (unsigned int n = 0; n<npol.ctk; n++){
for (unsigned int o = 0; o<npol.phi; o++){
double result = coeffs_moments.at(npol.get_bin_in4D(l,m,n,o))
*(pow(q2_edges[h+1],l+1) - pow(q2_edges[h],l+1))/(l+1.0)
*(pow(ctl_edges[i+1],m+1) - pow(ctl_edges[i],m+1))/(m+1.0)
*(pow(ctk_edges[j+1],n+1) - pow(ctk_edges[j],n+1))/(n+1.0)
*(pow(phi_edges[k+1],o+1) - pow(phi_edges[k],o+1))/(o+1.0);
integral += result;
}
}
}
}
neps_hist.at(get_bin_in4D(h,i,j,k,edges)) = integral;
neps += integral;
}
}
}
}
spdlog::debug( "nsel {0:f}\n\t nEvents {1:d}\n\t neps {2:f}\n\t nsel/neps {3:f} ", nsel, nEvents, neps, nsel/neps );
TH1D* hpull = new TH1D("hpull", ";(N_{sel}-N_{fit})/#sqrt{N_{sel}};",100, -5.0, +5.0);
TH1D* hnsel = new TH1D("hnsel", ";N_{sel};",100, 0.0, +25.0);
TH1D* hnfit = new TH1D("hnfit", ";N_{fit};",100, 0.0, +25.0);
//rescale fitted efficiency to MC statistics
//This would be nicer in a vector of vectors or so, but oh well
for (unsigned int h = 0; h<edges.nq2bins; h++)
for (unsigned int i = 0; i<edges.nctlbins; i++)
for (unsigned int j = 0; j<edges.nctkbins; j++)
for (unsigned int k = 0; k<edges.nphibins; k++)
neps_hist.at(get_bin_in4D(h,i,j,k, edges)) *= nsel/neps;
//calculate chi2 and perform "pull" plot
double chi2 = 0.0;
//double errorscale = 1.01;
for (unsigned int h = 0; h<edges.nq2bins; h++){
for (unsigned int i = 0; i<edges.nctlbins; i++){
for (unsigned int j = 0; j<edges.nctkbins; j++){
for (unsigned int k = 0; k<edges.nphibins; k++){
unsigned int bin = get_bin_in4D(h,i,j,k, edges);
//nominal
chi2 += sqr(nsel_hist.at(bin)-neps_hist.at(bin))/fabs(neps_hist.at(bin));
hpull->Fill((nsel_hist.at(bin)-neps_hist.at(bin))/sqrt(fabs(neps_hist.at(bin))));
spdlog::trace("nsel_hist.at({0:d})\t{1:f}",bin,nsel_hist.at(bin));
spdlog::trace("neps_hist.at({0:d})\t{1:f}",bin,neps_hist.at(bin));
hnsel->Fill(nsel_hist.at(bin));
hnfit->Fill(neps_hist.at(bin));
}
}
}
}
spdlog::debug( "chi2\t {0:f}", chi2 );
//Calculate the number of parameters
int nparams = npol.q2*npol.ctl*npol.ctk*npol.phi-1;
spdlog::debug( "Nparams\t {0:d}", nparams );
if (assumectleven && assumePhiEven){
nparams = npol.q2*((npol.ctl+1)/2)*npol.ctk*((npol.phi+1)/2)-1;
}
else{
if (assumectleven) nparams = npol.q2*((npol.ctl+1)/2)*npol.ctk*npol.phi-1;
if (assumePhiEven) nparams = npol.q2*npol.ctl*npol.ctk*((npol.phi+1)/2)-1;
}
spdlog::debug( "Nparams (accounting for even ctl and phi)\t {0:d}", nparams );
//Get degrees of freedom
int nbins = edges.getNbins();
spdlog::debug( "Nbins\t {0:d}", nbins);
int ndof = nbins - nparams - 1;
spdlog::debug( "ndof\t {0:d}", ndof );
spdlog::debug( "chi2/ndof (accounting for even ctl and phi)\t {0:f}",chi2/ndof );
spdlog::debug( "chi2/Nbins\t\t {0:f}",chi2/nbins );
spdlog::debug( "pvalue using ndof\t {0:f}",TMath::Prob(chi2, ndof) );
spdlog::debug( "pvalue using nbins\t {0:f}", TMath::Prob(chi2, nbins) );
//=========================================
//
// save chi2 to global values:
//
//=========================================
globalvalues->angacccorr_chi2 = chi2;
globalvalues->angacccorr_ndof = ndof;
//recount of nparams
int ncount = 0;
for (unsigned int l = 0; l<npol.q2; l++)
for (unsigned int m = 0; m<npol.ctl; m++)
for (unsigned int n = 0; n<npol.ctk; n++)
for (unsigned int o = 0; o<npol.phi; o++)
if (coeffs_moments.at(npol.get_bin_in4D(l,m,n,o)) != 0.0) ncount++;
spdlog::debug( "Recount of Nparams gives:\t {0:d}",ncount );
if(write_eps){
plotAndSave(hpull,"cpull",subfolder + "pull"+appendix, "eps");
plotAndSave(hnsel,"csel",subfolder + "nsel"+appendix, "eps");
plotAndSave(hnfit,"cfit",subfolder + "nfit"+appendix, "eps");
}
delete hpull;
delete hnsel;
delete hnfit;
return;
}
void get_legendre(bin_edges edges, npolynom npol,
std::vector<event> eventsInQ2bin, double nWeightedEvents,
std::vector<double> &poly_moments,
bool fejerkernel, bool assumectleven, bool assumePhiEven,
bool checkSignificance, bool calculateCovariance){
//use Legendre polynomials
spdlog::info("DOLEGENDRE start");
//should def use this:
//calculates all legendre polynomials up to including degree n
//void fcnc::legendre(double x, int n, std::vector<double>& results)
//this is not the best way to do things
//should rather
//loop over events
//determine, for every event the legendre polynomials order 0 to nq2/nctl/nctk/nphi
//loop over all orders
//add weight*leg_q2.at(h)*leg_ctl.at(i)*leg_ctk.at(j)*leg_phi.at(k) in correct bin
//optimized version
std::vector<double> chijk_moments = std::vector<double>(npol.getSize(), 0.0);
for (auto & meas: eventsInQ2bin){ //This is all computed then in the range of -1 to +1
double q2 = 2.0 * (meas.q2 - edges.q2min ) / edges.range_q2() - 1.0;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin) / edges.range_ctl()- 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin) / edges.range_ctk()- 1.0;
double phi = 2.0 * (meas.phi - edges.phimin) / edges.range_phi()- 1.0;
double weight = meas.weight;
//calculates all legendre polynomials up to including degree n
std::vector<double> leg_q2(npol.q2, 0.0);
fcnc::legendre(q2, npol.q2, leg_q2);
std::vector<double> leg_ctl(npol.ctl, 0.0);
fcnc::legendre(ctl, npol.ctl, leg_ctl);
std::vector<double> leg_ctk(npol.ctk, 0.0);
fcnc::legendre(ctk, npol.ctk, leg_ctk);
std::vector<double> leg_phi(npol.phi, 0.0);
fcnc::legendre(phi, npol.phi, leg_phi);
double result = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
unsigned int bin = npol.get_bin_in4D(h,i,j,k);
result = 0.0;
//result
if (assumectleven && i%2 != 0 ) result = 0.0;
else if (assumePhiEven && k%2 != 0) result = 0.0;
else{
result =
weight
*leg_q2.at(h)*(2.0*h+1.0)/2.0
*leg_ctl.at(i)*(2.0*i+1.0)/2.0
*leg_ctk.at(j)*(2.0*j+1.0)/2.0
*leg_phi.at(k)*(2.0*k+1.0)/2.0
/nWeightedEvents;
if (fejerkernel)
result *= (1.0 - h/double(npol.q2-1))
*(1.0 - i/double(npol.ctl-1))
*(1.0 - j/double(npol.ctk-1))
*(1.0 - k/double(npol.phi-1));
}
chijk_moments.at(bin) += result;
}
}
}
}
}
//check significance
//can be optimized as above
if (checkSignificance){
spdlog::info("starting CHECK SIGNIFICANT for LEGENDRE" );
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
unsigned int bin = npol.get_bin_in4D(h,i,j,k);
double result = chijk_moments.at(bin);
double sigma = 0.0;
for (unsigned int e=0; e<eventsInQ2bin.size(); e++){
const event & meas = eventsInQ2bin.at(e);
double q2 = 2.0 * (meas.q2 - edges.q2min ) / edges.range_q2() - 1.0;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin ) / edges.range_ctl() - 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin ) / edges.range_ctk() - 1.0;
double phi = 2.0 * (meas.phi - edges.phimin ) / edges.range_phi() - 1.0;
sigma += sqr(result - meas.weight // /nWeightedEvents
*fcnc::legendre(q2, h)*(2.0*h+1.0)/2.0
*fcnc::legendre(ctl, i)*(2.0*i+1.0)/2.0
*fcnc::legendre(ctk, j)*(2.0*j+1.0)/2.0
*fcnc::legendre(phi, k)*(2.0*k+1.0)/2.0
);// /eventsInQ2bin.size();
}
sigma = sqrt(sigma);
sigma /= nWeightedEvents;
if (fabs(result/sigma) < 1.0) chijk_moments.at(bin) = 0.0;
}
}
}
}
spdlog::info( "finished CHECK SIGNIFICANT for LEGENDRE" );
}
//Log
spdlog::debug( "//vector of legendre efficiency coefficients containing {0:d} elements", chijk_moments.size());
spdlog::trace( "legendre_coeffs_eff_4d = { " );
if (spdlog_trace()){
for (unsigned int i=0; i<chijk_moments.size(); i++){
std::cout<< std::setprecision(15) << std::scientific << chijk_moments.at(i) << (i<chijk_moments.size()-1 ? ", " : " ");
}
}
//calculate covariance
if (calculateCovariance){
const unsigned int nsize = npol.getSize();
spdlog::info( "Start COVARIANCE matrix for LEGENDRE polynoms" );
spdlog::debug( "calculating {0:d} by {1:d} convariance matrix ({2:d} doubles)" , nsize, nsize, nsize*nsize );
//optimized
std::vector<double> covariance(nsize*nsize, 0.0);
//startB = clock();
for (unsigned int e=0; e<eventsInQ2bin.size(); e++){
if (e % (eventsInQ2bin.size()/100) == 0)spdlog::debug( "{0:f}%", round(double(e*100)/eventsInQ2bin.size()));
const event & meas = eventsInQ2bin.at(e);
double q2 = 2.0 * (meas.q2 - edges.q2min ) / edges.range_q2() - 1.0;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin) / edges.range_ctl()- 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin) / edges.range_ctk()- 1.0;
double phi = 2.0 * (meas.phi - edges.phimin) / edges.range_phi()- 1.0;
double weight = meas.weight;
//calculates all legendre polynomials up to including degree n
std::vector<double> leg_q2(npol.q2, 0.0);
legendre(q2, npol.q2, leg_q2);
std::vector<double> leg_ctl(npol.ctl, 0.0);
legendre(ctl, npol.ctl, leg_ctl);
std::vector<double> leg_ctk(npol.ctk, 0.0);
legendre(ctk, npol.ctk, leg_ctk);
std::vector<double> leg_phi(npol.phi, 0.0);
legendre(phi, npol.phi, leg_phi);
double wfactor = weight*weight/(nWeightedEvents*nWeightedEvents);
double resulti=0.0, resultj=0.0, cov=0.0;
unsigned int bini=0, binj=0;
int hfactor=0, ifactor=0, jfactor=0, kfactor=0,
lfactor=0, mfactor=0, nfactor=0, ofactor=0,
hifactor=0, hijfactor=0, hijkfactor=0,
lmfactor=0, lmnfactor=0, lmnofactor=0;
double aq2=0.0,aq2ctl=0.0,aq2ctlctk=0.0,aq2ctlctkphi=0.0;
double bq2=0.0,bq2ctl=0.0,bq2ctlctk=0.0,bq2ctlctkphi=0.0;
for (unsigned int h = 0; h<npol.q2; h++){
hfactor = (2*h+1);
aq2 = 0.0625*leg_q2[h];
for (unsigned int i = 0; i<npol.ctl; i++){
ifactor = (2*i+1);
hifactor = hfactor*ifactor;
aq2ctl = aq2*leg_ctl[i];
for (unsigned int j = 0; j<npol.ctk; j++){
jfactor = (2*j+1);
hijfactor = hifactor*jfactor;
aq2ctlctk = aq2ctl*leg_ctk[j];
for (unsigned int k = 0; k<npol.phi; k++){
kfactor = (2*k+1);
hijkfactor = hijfactor*kfactor;
aq2ctlctkphi = aq2ctlctk*leg_phi[k];
bini = npol.get_bin_in4D(h,i,j,k);
resulti = chijk_moments[bini];
if (resulti != 0.0){
for (unsigned int l = 0; l<npol.q2; l++){
lfactor = (2*l+1);
bq2 = 0.0625*leg_q2[l];
for (unsigned int m = 0; m<npol.ctl; m++)
{
mfactor = (2*m+1);
lmfactor = lfactor*mfactor;
bq2ctl = bq2*leg_ctl[m];
for (unsigned int n = 0; n<npol.ctk; n++){
nfactor = (2*n+1);
lmnfactor = lmfactor*nfactor;
bq2ctlctk = bq2ctl*leg_ctk[n];
for (unsigned int o = 0; o<npol.phi; o++){
ofactor = (2*o+1);
lmnofactor = lmnfactor*ofactor;
bq2ctlctkphi = bq2ctlctk*leg_phi[o];
binj = npol.get_bin_in4D(l,m,n,o);
resultj = chijk_moments[binj];
if (binj <= bini && resultj != 0.0){
cov = (wfactor
*(resulti-hijkfactor*aq2ctlctkphi)
*(resultj-lmnofactor*bq2ctlctkphi)
);
covariance[nsize*bini+binj] += cov;
}
}
}
}
}
}
}
}
}
}
}
for (unsigned int bini = 0; bini<nsize; bini++){
for (unsigned int binj = 0; binj<nsize; binj++){
//covariance matrix is symmetric
if (binj < bini) covariance.at(nsize*binj+bini) = covariance.at(nsize*bini+binj);
}
}
//Print uncertanities in a matrix
spdlog::info("uncertainties[{0:0.8e}]", nsize);
for (unsigned int i =0; i<nsize; i++) std::cout << covariance.at(nsize*i+i) << (i < nsize-1 ? ", " : " };");
//Print covariance matrix
spdlog::info("Covariance matrix cov[{0:0.8e}]", nsize*nsize);
for (unsigned int i =0; i<nsize; i++){
for (unsigned int j =0; j<nsize; j++){
std::cout <<covariance.at(nsize*i+j) << ", ";
}
std::cout << std::endl;
}
std::cout << "};" << std::endl;//needs to be converted into correct range and translated to monomial coefficients? no, should do these transformations after every acceptance randomization
double correlation = 0.0;
for (unsigned int i =0; i<nsize; i++){
for (unsigned int j =0; j<nsize; j++){
double corr = covariance.at(nsize*i+j)/sqrt(covariance.at(nsize*i+i))/sqrt(covariance.at(nsize*j+j));
if (fabs(corr) > fabs(correlation) && i != j) correlation = corr;
}
}
spdlog::info("Largest correlation found is {0:f}", correlation );
}
spdlog::debug("start CONVERT legendre to monomial coefficients" );
//transform legendre coefficients to monomial coefficients
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double tmp = chijk_moments.at(npol.get_bin_in4D(h,i,j,k));
std::vector<double> leg_coeffs_q2(npol.q2, 0.0);
std::vector<double> leg_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> leg_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> leg_coeffs_phi(npol.phi, 0.0);
std::vector<double> poly_coeffs_q2(npol.q2, 0.0);
std::vector<double> poly_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> poly_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> poly_coeffs_phi(npol.phi, 0.0);
leg_coeffs_q2.at(h) = 1.0;
leg_coeffs_ctl.at(i) = 1.0;
leg_coeffs_ctk.at(j) = 1.0;
leg_coeffs_phi.at(k) = 1.0;
legendre_to_poly(leg_coeffs_q2, poly_coeffs_q2);
legendre_to_poly(leg_coeffs_ctl, poly_coeffs_ctl);
legendre_to_poly(leg_coeffs_ctk, poly_coeffs_ctk);
legendre_to_poly(leg_coeffs_phi, poly_coeffs_phi);
for (unsigned int l = 0; l<=h; l++)
for (unsigned int m = 0; m<=i; m++)
for (unsigned int n = 0; n<=j; n++)
for (unsigned int o = 0; o<=k; o++){
unsigned int bin = npol.get_bin_in4D(l,m,n,o);
poly_moments.at(bin) += tmp*poly_coeffs_q2.at(l)*poly_coeffs_ctl.at(m)*poly_coeffs_ctk.at(n)*poly_coeffs_phi.at(o);
}
}//tested, works
spdlog::debug( "end CONVERT legendre to monomial coefficients" );
spdlog::debug( "poly_moments size: {0:d}", poly_moments.size() );
spdlog::info( "DOLEGENDRE end" );
}
void get_chebyschev(bin_edges edges, npolynom npol,
std::vector<event> eventsInQ2bin, double nWeightedEvents,
std::vector<double> &poly_moments,
bool fejerkernel, bool assumectleven, bool assumePhiEven,
bool checkSignificance, bool calculateCovariance
){
//use Legendre polynomials
spdlog::info("CHEBYSCHEV start");
std::vector<double> chijk_moments = std::vector<double>(npol.getSize(), 0.0);
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double result = 0.0;
if (assumectleven && i%2 != 0)result = 0.0;
else if (assumePhiEven && k%2 != 0)result = 0.0;
else{
for (auto &meas: eventsInQ2bin){
double q2 = 2.0 * (meas.q2 - edges.q2min ) / edges.range_q2() - 1.0;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin) / edges.range_ctl()- 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin) / edges.range_ctk()- 1.0;
double phi = 2.0 * (meas.phi - edges.phimin) / edges.range_phi()- 1.0;
result += meas.weight
*fcnc::chebyshev(q2, h)/sqrt(1.0-q2*q2)*(h == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctl, i)/sqrt(1.0-ctl*ctl)*(i == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctk, j)/sqrt(1.0-ctk*ctk)*(j == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(phi, k)/sqrt(1.0-phi*phi)*(k == 0 ? 1.0/MY_PI : 2.0/MY_PI)
/nWeightedEvents;// /eventsInQ2bin.size();
}
if (fejerkernel){
result *= (1.0 - h/double(npol.q2-1))
*(1.0 - i/double(npol.ctl-1))
*(1.0 - j/double(npol.ctk-1))
*(1.0 - k/double(npol.phi-1));
}
}
//now have the "mean" in result
//mu = sum_i w_i/N
//sigma(mu) = sqrt(sum_i (w_i-mu)^2)/N
//else
chijk_moments.push_back(result);
}
//check significance
if (checkSignificance){
spdlog::info( "start CHECK SIGNIFICANCE for CHEBYSHEV polynoms" );
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
unsigned int bin = npol.get_bin_in4D(h,i,j,k);
double result = chijk_moments.at(bin);
double sigma = 0.0;
for (auto &meas : eventsInQ2bin){
double q2 = 2.0 * (meas.q2 - edges.q2min ) / edges.range_q2() - 1.0;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin) / edges.range_ctl()- 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin) / edges.range_ctk()- 1.0;
double phi = 2.0 * (meas.phi - edges.phimin) / edges.range_phi()- 1.0;
sigma += sqr(result - meas.weight // /nWeightedEvents
*fcnc::chebyshev(q2, h)/sqrt(1-q2*q2)*(h == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctl, i)/sqrt(1-ctl*ctl)*(i == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctk, j)/sqrt(1-ctk*ctk)*(j == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(phi, k)/sqrt(1-phi*phi)*(k == 0 ? 1.0/MY_PI : 2.0/MY_PI));// /eventsInQ2bin.size();
}
sigma = sqrt(sigma);
sigma /= nWeightedEvents;
if (fabs(result/sigma) < 1.0) chijk_moments.at(bin) = 0.0;
}
spdlog::info( "end CHECK SIGNIFICANCE for CHEBYSHEV polynoms" );
}
//printchebyshev
spdlog::debug( "//vector of chebyshev efficiency coefficients containing {0:d} elements",chijk_moments.size());
spdlog::trace( "chebyshev_coeffs_eff_4d = { " );
if (spdlog_trace()){
for (unsigned int i=0; i<chijk_moments.size(); i++){
std::cout<< std::setprecision(15) << std::scientific << chijk_moments.at(i) << (i<chijk_moments.size()-1 ? ", " : " ");
}
}
//calculate covariance
if (calculateCovariance){
spdlog::debug( "calculating {0:d} by {1:d} convariance matrix ({2:d} doubles)" , npol.getSize(), npol.getSize(), npol.getSize2() );
std::vector<double> covariance(npol.getSize2(), 0.0);
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
unsigned int bini = npol.get_bin_in4D(h,i,j,k);
double resulti = chijk_moments.at(bini);
if (resulti != 0.0){
for (unsigned int l = 0; l<npol.q2; l++){
for (unsigned int m = 0; m<npol.ctl; m++){
for (unsigned int n = 0; n<npol.ctk; n++){
for (unsigned int o = 0; o<npol.phi; o++){
unsigned int binj = npol.get_bin_in4D(l,m,n,o);
if (binj > bini) continue;//covariance matrix is symmetric
double resultj = chijk_moments.at(binj);
if (((bini*npol.getSize()+binj) * 100) % (npol.getSize2()) == 0)
spdlog::debug( "covariance matrix calculation {0:f}%", ((bini*npol.getSize()+binj) * 100) / (npol.getSize2()));
if (resultj != 0.0){
double cov = 0.0;
for (auto &meas : eventsInQ2bin){
double q2 = 2.0 * (meas.q2 - edges.q2min ) / edges.range_q2() - 1.0;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin) / edges.range_ctl()- 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin) / edges.range_ctk()- 1.0;
double phi = 2.0 * (meas.phi - edges.phimin) / edges.range_phi()- 1.0;
cov += ((resulti - meas.weight
*fcnc::chebyshev(q2, h)/sqrt(1-q2*q2)*(h == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctl, i)/sqrt(1-ctl*ctl)*(i == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctk, j)/sqrt(1-ctk*ctk)*(j == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(phi, k)/sqrt(1-phi*phi)*(k == 0 ? 1.0/MY_PI : 2.0/MY_PI))
*(resultj - meas.weight
*fcnc::chebyshev(q2, l)/sqrt(1-q2*q2)*(l == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctl, m)/sqrt(1-ctl*ctl)*(m == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctk, n)/sqrt(1-ctk*ctk)*(n == 0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(phi, o)/sqrt(1-phi*phi)*(o == 0 ? 1.0/MY_PI : 2.0/MY_PI))
);
}
//cov = sqrt(cov); //to get sigma on diagonal
cov /= sqr(nWeightedEvents);
covariance.at(npol.getSize()*bini+binj) = cov;
covariance.at(npol.getSize()*binj+bini) = cov;
}
}
}
}
}
}
}
}
}
}
spdlog::info("uncertainties[{0:0.8e}]", npol.getSize());
if (spdlog_info()){
for (unsigned int i =0; i<npol.getSize(); i++){
std::cout <<covariance.at(npol.getSize()*i+i) << (i < npol.getSize()-1 ? ", " : " };");
}
std::cout << std::endl;//needs to be converted into correct range and translated to monomial coefficients?
}
spdlog::info("Covariance matrix cov[{0:0.8e}]", npol.getSize2());
for (unsigned int i =0; i<npol.getSize(); i++){
for (unsigned int j =0; j<npol.getSize(); j++){
std::cout << covariance.at(npol.getSize()*i+j) << ", ";
}
std::cout << std::endl;
}
std::cout <<"};" << std::endl;//needs to be converted into correct range and translated to monomial coefficients? no, should do these transformations after every acceptance randomization
double correlation = 0.0;
for (unsigned int i =0; i<npol.getSize(); i++)
for (unsigned int j =0; j<npol.getSize(); j++)
{
double corr = covariance.at(npol.getSize()*i+j)/sqrt(covariance.at(npol.getSize()*i+i))/sqrt(covariance.at(npol.getSize()*j+j));
if (fabs(corr) > fabs(correlation) && i != j)
correlation = corr;
}
spdlog::info( "Largest correlation found is {0:f}", correlation );
}
//transform chebyshev coefficients to monomial coefficients
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
double c = chijk_moments.at(npol.get_bin_in4D(h,i,j,k));
std::vector<double> ch_coeffs_q2(npol.q2, 0.0);
std::vector<double> ch_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> ch_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> ch_coeffs_phi(npol.phi, 0.0);
std::vector<double> poly_coeffs_q2(npol.q2, 0.0);
std::vector<double> poly_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> poly_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> poly_coeffs_phi(npol.phi, 0.0);
ch_coeffs_q2.at(h) = 1.0;
ch_coeffs_ctl.at(i) = 1.0;
ch_coeffs_ctk.at(j) = 1.0;
ch_coeffs_phi.at(k) = 1.0;
chebyshev_to_poly(ch_coeffs_q2, poly_coeffs_q2);
chebyshev_to_poly(ch_coeffs_ctl, poly_coeffs_ctl);
chebyshev_to_poly(ch_coeffs_ctk, poly_coeffs_ctk);
chebyshev_to_poly(ch_coeffs_phi, poly_coeffs_phi);
for (unsigned int l = 0; l<=h; l++){
for (unsigned int m = 0; m<=i; m++){
for (unsigned int n = 0; n<=j; n++){
for (unsigned int o = 0; o<=k; o++){
unsigned int bin = npol.get_bin_in4D(l,m,n,o);
poly_moments.at(bin) += c*poly_coeffs_q2.at(l)*poly_coeffs_ctl.at(m)*poly_coeffs_ctk.at(n)*poly_coeffs_phi.at(o);
}
}
}
}
}
}
}
}
}
double integrate_4D(bin_edges edges, int h, int i, int j, int k){
return integrate_x_to_n(edges.q2min, edges.q2max,h)//careful with q2max
*integrate_x_to_n(edges.ctlmin, edges.ctlmax,i)
*integrate_x_to_n(edges.ctkmin, edges.ctkmax,j)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
}
double integrate_3D(bin_edges edges, int i, int j, int k){
return integrate_x_to_n(edges.ctlmin, edges.ctlmax,i)
*integrate_x_to_n(edges.ctkmin, edges.ctkmax,j)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
}
void fill_moments_hist_1D(TH1D *hist_moments, std::vector<double> coeffs_moments,
bin_edges edges, npolynom npol, double q2low, double q2high){
bool do_ctl = (std::string(hist_moments->GetName()).find("ctl") != std::string::npos);
bool do_ctk = (std::string(hist_moments->GetName()).find("ctk") != std::string::npos);
bool do_phi = (std::string(hist_moments->GetName()).find("phi") != std::string::npos);
bool do_q2 = (std::string(hist_moments->GetName()).find("q2") != std::string::npos);
for (int l=1; l<=hist_moments->GetNbinsX(); l++){
double a = hist_moments->GetBinLowEdge(l);
double b = hist_moments->GetBinLowEdge(l+1);
double integral = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
double result = coeffs_moments.at(npol.get_bin_in4D(h,i,j,k))
*(pow(do_q2 ? b : q2high ,h+1)/(h+1.0) - pow(do_q2 ? a : q2low ,h+1)/(h+1.0)) //Possibly replace by get_result;
*(pow(do_ctl ? b : edges.ctlmax,i+1)/(i+1.0) - pow(do_ctl ? a : edges.ctlmin,i+1)/(i+1.0))
*(pow(do_ctk ? b : edges.ctkmax,j+1)/(j+1.0) - pow(do_ctk ? a : edges.ctkmin,j+1)/(j+1.0))
*(pow(do_phi ? b : edges.phimax,k+1)/(k+1.0) - pow(do_phi ? a : edges.phimin,k+1)/(k+1.0));
integral += result;
}
}
}
}
hist_moments->SetBinContent(l, integral);
}
return;
}
void fill_moments_hist_2D(TH2D *hist_moments, std::vector<double> coeffs_moments,
bin_edges edges, npolynom npol, double q2low, double q2high){
bool do_ctl = (std::string(hist_moments->GetName()).find("ctl") != std::string::npos);
bool do_ctk = (std::string(hist_moments->GetName()).find("ctk") != std::string::npos);
bool do_phi = (std::string(hist_moments->GetName()).find("phi") != std::string::npos);
bool do_q2 = (std::string(hist_moments->GetName()).find("q2") != std::string::npos);
spdlog::trace("do_ctl\t"+boolToString(do_ctl));
spdlog::trace("do_ctk\t"+boolToString(do_ctk));
spdlog::trace("do_phi\t"+boolToString(do_phi));
spdlog::trace("do_q2\t"+boolToString(do_q2));
bool is_ctl_First = !do_q2; //ctl is second if q2
bool is_ctk_First = !(do_q2 || do_ctl); //ctk is second if ctl or q2
//phi is always second and q2 always first
spdlog::trace("is_ctl_First\t"+boolToString(is_ctl_First));
spdlog::trace("is_ctk_First\t"+boolToString(is_ctk_First));
for (int l=1; l<=hist_moments->GetNbinsX(); l++){
double a = hist_moments->GetXaxis()->GetBinLowEdge(l);
double b = hist_moments->GetXaxis()->GetBinLowEdge(l+1);
for (int m=1; m<=hist_moments->GetNbinsY(); m++){
double c = hist_moments->GetYaxis()->GetBinLowEdge(m);
double d = hist_moments->GetYaxis()->GetBinLowEdge(m+1);
double integral = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
double result = coeffs_moments.at(npol.get_bin_in4D(h,i,j,k))
*integrate_x_to_n(do_q2 ? a : q2low, do_q2 ? b : q2high ,h)
*integrate_x_to_n(do_ctl ? (is_ctl_First ? a : c) : edges.ctlmin,
do_ctl ? (is_ctl_First ? b : d) : edges.ctlmax,i)
*integrate_x_to_n(do_ctk ? (is_ctk_First ? a : c) : edges.ctkmin,
do_ctk ? (is_ctk_First ? b : d) : edges.ctkmax, j)
*integrate_x_to_n(do_phi ? c : edges.phimin,do_phi ? d : edges.phimax,k);
integral += result;
}
}
}
}
spdlog::trace("\t*** fill_moments_hist_2D: ***\t´");
spdlog::trace("l: {0:d}\tm: {1:d}\tint: {2:f}", l, m, integral);
hist_moments->SetBinContent(l, m, integral);
}
}
return;
}
int get_correct_index(std::string whichVar, int l, int m, int n, int o){
if (whichVar == "q2") return l;
if (whichVar == "ctl") return m;
if (whichVar == "ctk") return n;
if (whichVar == "phi") return o;
return -1;
}
int change_index(std::string whichVar, npolynom npol, int p, int l, int m, int n, int o){
if (whichVar == "q2") return npol.get_bin_in4D(p,m,n,o);
if (whichVar == "ctl") return npol.get_bin_in4D(l,p,n,o);
if (whichVar == "ctk") return npol.get_bin_in4D(l,m,p,o);
if (whichVar == "phi") return npol.get_bin_in4D(l,m,n,p);
return -1;
}
void correct_coefficients(npolynom npol, bin_edges edges, std::vector<double> &out_poly_moments, std::vector<double> in_poly_moments, std::string whichVar){
//Just rebin
unsigned int npol_var;
double intvar_min;
double intvar_max;
if (whichVar == "q2"){ //If this would be na efin array, my life would be orders of magnitute easier
npol_var = npol.q2;
intvar_min = edges.q2min;
intvar_max = edges.q2max;
}
else if (whichVar == "ctl") {
npol_var = npol.ctl;
intvar_min = edges.ctlmin;
intvar_max = edges.ctlmax;
}
else if (whichVar == "ctk") {
npol_var = npol.ctk;
intvar_min = edges.ctkmin;
intvar_max = edges.ctkmax;
}
else if (whichVar == "phi") {
npol_var = npol.phi;
intvar_min = edges.phimin;
intvar_max = edges.phimax;
}
else{
spdlog::error("Wrong index name! Either q2, ctl, ctk or phi!");
assert (0);
}
spdlog::debug("npol_var = {0:d}\tintvar_min = {1:f}\tintvar_max = {2:f}", npol_var, intvar_min, intvar_max);
for (unsigned int l = 0; l<npol.q2; l++){
for (unsigned int m = 0; m<npol.ctl; m++){
for (unsigned int n = 0; n<npol.ctk; n++){
for (unsigned int o = 0; o<npol.phi; o++){
unsigned int bin = npol.get_bin_in4D(l,m,n,o);
spdlog::trace("l={0:d}\tm={1:d}\tn={2:d}\to={3:d}\t bin ={4:d}",l,m,n,o,bin);
std::vector<double> poly(npol_var, 0.0);
poly.at(get_correct_index(whichVar,l,m,n,o)) = in_poly_moments.at(bin);
std::vector<double> correct(npol_var, 0.0);
correct_poly(poly, correct, intvar_min, intvar_max);
for (unsigned int p = 0; p<npol_var; p++){
unsigned int tobin = change_index(whichVar,npol,p,l,m,n,o);
spdlog::trace("p={0:d}\t tobin ={1:d}",p,tobin);
out_poly_moments.at(tobin) += correct.at(p);
}
}
}
}
}
return;
}
void bu2kstarmumu_pdf::parametrize_eff_phsp_4d(const std::vector<event>& events, values* globalvalues, int tagNumber, bool assumePhiEven, bool checkSignificance, bool runMinuit, bool checkFactorization, bool do3Dmoments){
spdlog::debug("----> start parametrization of angular acceptance correction coefficients with {0:d} events.", events.size());
const UInt_t nBinsQ2 = 40;
const UInt_t nBinsAngles = 40;
bool ignore_q2weights = false;
if (ignore_q2weights) spdlog::warn("Not using weights from genLvl PHSP to flatten the q2 distribution!!!");
//Get the pointer to this class of bu2kstarmumu_pdf
current_pdf = this;
std::vector<double> q2binsmin = opts->TheQ2binsmin; //Check TheQ2binsmax and q2binsmin
std::vector<double> q2binsmax = opts->TheQ2binsmax;
spdlog::debug("Vector q2binsmin:\t"+convert_vector_to_string(q2binsmin));
spdlog::debug("Vector q2binsmax:\t"+convert_vector_to_string(q2binsmax));
//name tag for all plots:
std::string appendix = "_"+DECAY_NAME
+ (opts->KS ? opts->DDLL+"_" : "_")
+ "Run" + std::to_string(opts->run)
+ ((opts->angacccorrperyear && opts->year != -1) ? "_"+std::to_string(opts->year) : "")
+ ((opts->folding != -1) ? "_folding"+std::to_string(opts->folding) : "")
+ ((opts->systematic == 3) ? "_"+std::to_string(opts->orderincrease)+"moreOrders" : "")
+ std::string(tagNumber != 0 ? "_"+std::to_string(tagNumber) : "");
std::string subfolder = "plots/angular/" + std::string(tagNumber != 0 ? "test/" : "");
spdlog::debug("Saving all plots into folder " + subfolder);
//Set orders of chebyshev/legendre polynomials
npolynom npol(opts); //Should be a const, but then all the future shuffling gets difficult
//fold the events if needed:
spdlog::debug("Init folder....");
//TODO remove fcnc::folder fldr(opts);
//only accept events inside q2min and q2max
std::vector<event> filtered_events;
//std::vector<event> fullrange_events;
unsigned int ee;
unsigned int N = events.size();
TRandom3 * rnd = new TRandom3(15*opts->job_id);
spdlog::debug("Filtering events...");
for (unsigned int e=0; e<N; e++){
//choose random events for bootstrapping
if (opts->systematic == 1) ee = int(rnd->Rndm() * N);
else ee = e; //or the correct, nominal order
fcnc::event evt = events.at(ee);
if (!isEvtInAngleRange(&(evt), opts)) continue;
filtered_events.push_back(evt);
if(!opts->full_angular) fldr->fold(&filtered_events.back());
}
if (!ignore_q2weights){
spdlog::debug("Getting q2 weights...");
//weight for all events, reweight from PHSP q2 model to flat q2 distribution using histogram
TFile* fphspweightq2 = new TFile(CONFIG_PHSP_WEIGHT.c_str());
TH1D* phspweightq2 = dynamic_cast<TH1D*>(fphspweightq2->Get(opts->IsFlatQ2 ? "flatq2weights" : "phspweightq2"));
bool keep_weights = opts->systematic != 2;
bool vary_weights = opts->systematic == 4;
//reweight filtered events
for (auto evt: filtered_events){
double q2 = evt.q2;
if(vary_weights){
evt.weight += evt.delta_weight * rnd->Gaus(0., 1.);
}
if (q2 < 0.2){
if(keep_weights) evt.weight *= phspweightq2->GetBinContent(2);//keep loaded totalweight!
else evt.weight = phspweightq2->GetBinContent(2);
}
else{
int bin = phspweightq2->FindFixBin(q2);
if (bin >=1 && bin <= phspweightq2->GetNbinsX()){
if(keep_weights) evt.weight *= phspweightq2->GetBinContent(bin);//keep loaded totalweight!
else evt.weight = phspweightq2->GetBinContent(bin);
}
}
}
delete phspweightq2;
fphspweightq2->Close();
delete fphspweightq2;
}
//for systematic 9: trigger weights, multiply existing weights by 1/L0Muon_eff
if(opts->systematic == 9){
spdlog::debug("Getting systematic 9");
//get years from current data-set:
std::vector<int> years; //empty vector
for (auto event: filtered_events){
bool year_found = false;
for(auto yr: years){
if(yr == event.year){
year_found = true;
break;
}
}
if(!year_found) years.push_back(event.year); //First event always gets filled
//Technically one should make an if to check it is all years, but given the size of the event sample this is okayish
}
assert(years.size());
std::sort(years.begin(), years.end());
//get L0Muon efficiency histo(s)
TFile * f_L0_w[years.size()];
TH1D * h_L0_w[years.size()];
for(unsigned int y = 0; y < years.size(); y++){
f_L0_w[y] = new TFile((std::string("config/")+(opts->KS?"KshortPiplus":"KplusPi0Resolved")+"_L0MuonEfficiency_"+std::to_string(years.at(y))+".root").c_str(), "READ");
if(!f_L0_w[y]->IsOpen()){
spdlog::critical( "Could not find file with L0Muon efficiencies: 'config/" + std::string(opts->KS ? "KshortPiplus" : "KplusPi0Resolved") + "_L0MuonEfficiency_{0:d}.root'", years.at(y) );
assert(0);
}
h_L0_w[y] = dynamic_cast<TH1D*>(f_L0_w[y]->Get("L0MuonEffRatioPHSP"));
spdlog::info("Succesfully loaded L0Muon efficiency reweighting histogram for year={0:d}", years.at(y));
}
for (unsigned int i=0; i<filtered_events.size(); i++){
double max_mu_pt = filtered_events.at(i).max_mu_pt;
unsigned int year_idx = 0;
for(unsigned int y = 0; y < years.size(); y++){
if(years.at(y) == filtered_events.at(i).year){
year_idx = y;
break;
}
}
if(max_mu_pt > 20000.){
filtered_events.at(i).weight *= h_L0_w[year_idx]->GetBinContent(h_L0_w[year_idx]->GetNbinsX());
}
else{
int bin = h_L0_w[year_idx]->FindFixBin(max_mu_pt);
filtered_events.at(i).weight *= h_L0_w[year_idx]->GetBinContent(bin);
}
}
for(unsigned int y = 0; y < years.size(); y++){
delete h_L0_w[y];
f_L0_w[y]->Close();
delete f_L0_w[y];
}
}
//set the pointer that is accessible from the minuit fcn
eff_events = filtered_events;
//make plot to check weight distribution if running a test
if(tagNumber >0){
spdlog::debug("Plotting weights");
TH1D* hweight = new TH1D("hweight", ";w;", 100, 0.0, 1.0/EFF_CUTOFF);
for (unsigned int i=0; i<filtered_events.size(); i++){
spdlog::trace("Weight of event {0:d}\t\t{1:f}", i , filtered_events.at(i).weight);
hweight->Fill(filtered_events.at(i).weight);
}
TCanvas * cweights = new TCanvas();
cweights->cd();
hweight->Draw("");
if(opts->write_eps)cweights->Print((subfolder + "weights"+appendix+".eps").c_str(), "eps");
if(opts->write_eps)cweights->Print((subfolder + "weights"+appendix+".root").c_str(), "root");
delete hweight;
delete cweights;
}
//---------------------------------------
//DETERMINATION FROM ACCEPTANCE HISTOGRAM
//---------------------------------------
unsigned int nq2bins = 18;
unsigned int nctlbins = 10;
unsigned int nctkbins = 10;
unsigned int nphibins = 5;
bool adapt_bins_to_folding = false;
if(adapt_bins_to_folding){
nphibins /= 2;
if(opts->folding > 0)nctlbins /= 2;
}
spdlog::debug("Options: ");
spdlog::debug("q2min={0:f}\tq2max={1:f}", opts->TheQ2binsmin.front(),opts->TheQ2binsmax.back());
spdlog::debug("ctkmin={0:f}\tctkmax={1:f}", opts->ctk_min, opts->ctk_max);
spdlog::debug("ctlmin={0:f}\tctlmax={1:f}", opts->ctl_min, opts->ctl_max);
spdlog::debug("phimin={0:f}\tphimax={1:f}", opts->phi_min, opts->phi_max);
bin_edges edges (opts, nq2bins, nctlbins, nctkbins, nphibins);
const std::vector<double> q2_edges = edges.edges_q2();
const std::vector<double> ctl_edges = edges.edges_ctl();
const std::vector<double> ctk_edges = edges.edges_ctk();
const std::vector<double> phi_edges = edges.edges_phi(); //This could've been a struct/class
spdlog::trace("q2 edges: " + convert_vector_to_string(q2_edges));
spdlog::trace("ctl edges: " + convert_vector_to_string(ctl_edges));
spdlog::trace("ctk edges: " + convert_vector_to_string(ctk_edges));
spdlog::trace("phi edges: " + convert_vector_to_string(phi_edges));
//should compare full 4D with
//3D in q2 bins
//3D factorizing in bins of q2
spdlog::debug("Getting number of weighted events...");
std::vector<double> eps_hist(edges.getNbins(), 0.0);
double nWeightedEvents = 0.0;
for (auto &meas: filtered_events){
unsigned int bin = get_bin_in4D(meas,edges);
eps_hist.at(bin) += meas.weight;
nWeightedEvents += meas.weight;
}
//rescale efficiency
spdlog::debug("Rescaling efficiency...");
for (unsigned int h = 0; h<edges.nq2bins; h++){
for (unsigned int i = 0; i<edges.nctlbins; i++){
for (unsigned int j = 0; j<edges.nctkbins; j++){
for (unsigned int k = 0; k<edges.nphibins; k++){
eps_hist.at(get_bin_in4D(h,i,j,k,edges)) *= edges.getNbins()/nWeightedEvents;
}
}
}
}
//determine xis in bins of q2
spdlog::debug("Getting xis in bins of q2...");
std::vector<TH1D*> hxi_hist;
for (auto angObs: get_angObser_withTeX_vec()){
hxi_hist.push_back(new TH1D(("hxi"+angObs[0]+"hist").c_str(), (";q^{2};"+angObs[1]).c_str(),
edges.nq2bins, edges.q2min, edges.q2max));
}
std::vector<int>int_xi(hxi_hist.size(),0);
for (unsigned int h = 0; h<edges.nq2bins; h++){
std::fill(int_xi.begin(), int_xi.end(), 0); //fill int_xi with zeroes
for (unsigned int i = 0; i<edges.nctlbins; i++)
for (unsigned int j = 0; j<edges.nctkbins; j++)
for (unsigned int k = 0; k<edges.nphibins; k++){
unsigned int bin = get_bin_in4D(h,i,j,k,edges);
double eps = eps_hist.at(bin);
double f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12;
if(opts->full_angular){
integrated_fj_noacc(ctl_edges[i], ctl_edges[i+1],
ctk_edges[j], ctk_edges[j+1],
phi_edges[k], phi_edges[k+1],
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
}
else{
folded_integrated_fj_noacc(ctl_edges[i], ctl_edges[i+1],
ctk_edges[j], ctk_edges[j+1],
phi_edges[k], phi_edges[k+1],
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
}
int_xi[0] += eps*f1; //KILL ME
int_xi[1] += eps*f2;
int_xi[2] += eps*f3;
int_xi[3] += eps*f4;
int_xi[4] += eps*f5;
int_xi[5] += eps*f6;
int_xi[6] += eps*f7;
int_xi[7] += eps*f8;
int_xi[8] += eps*f9;
int_xi[9] += eps*f10;
int_xi[10] += eps*f11;
int_xi[11] += eps*f12;
}
for_indexed(auto hist: hxi_hist) hist->SetBinContent(h+1,int_xi[i]);
}
//-----------------------------------
//DETERMINE XIS USING UNBINNED METHOD
//-----------------------------------
spdlog::debug("Getting xis using unbinned method...");
std::vector<TH1D*> hxi_ub;
for (auto angObs: get_angObser_withTeX_vec()){
hxi_ub.push_back(new TH1D(("hxi"+angObs[0]+"ub").c_str(), (";q^{2};"+angObs[1]).c_str(),
edges.nq2bins, edges.q2min, edges.q2max));
}
//make it pretty
for (auto hist: hxi_ub) hist->SetLineColor(4);
for (unsigned int h = 0; h<edges.nq2bins; h++){
std::fill(int_xi.begin(), int_xi.end(), 0); //fill int_xi with zeroes, no need to initialize (see line 961)
//double q2a = hxi_ub[0]->GetBinLowEdge(h); //check this is the same as q2_edges!!!!! should be, but to be sure
//double q2b = hxi_ub[0]->GetBinLowEdge(h+1);
for (auto &&meas: filtered_events){
if (meas.q2 < q2_edges[h] || meas.q2 > q2_edges[h]) continue;
double f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12;
if(opts->full_angular){
fj(meas.costhetal, meas.costhetak, meas.phi,
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
}
else{
folded_fj(meas.costhetal, meas.costhetak, meas.phi,
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
}
double weight = meas.weight/nWeightedEvents*edges.nq2bins*(edges.range_ctl()*edges.range_ctk()*edges.range_phi());
int_xi[0] += f1*weight; //KILL ME
int_xi[1] += f2*weight;
int_xi[2] += f3*weight;
int_xi[3] += f4*weight;
int_xi[4] += f5*weight;
int_xi[5] += f6*weight;
int_xi[6] += f7*weight;
int_xi[7] += f8*weight;
int_xi[8] += f9*weight;
int_xi[9] += f10*weight;
int_xi[10] += f11*weight;
int_xi[11] += f12*weight;
}
for_indexed(auto hist: hxi_ub) hist->SetBinContent(h,int_xi[i]);
}
//------------------------------------------------------
//DETERMINE POLYNOMIAL ACCEPTANCE WITH METHOD OF MOMENTS
//------------------------------------------------------
const bool calculateCovariance = false;
const bool assumectleven = opts->systematic == 2;
const bool fejerkernel = false; //see Alexander Weisse, Holger Feshke, Chebyshev expansion techniques
const bool doLegendre = true; //somehow, chebyshev is broken now, great job me
std::vector<double> poly_moments(npol.getSize(), 0.0);
if (doLegendre){
get_legendre(edges, npol,
filtered_events, nWeightedEvents,
poly_moments,
fejerkernel, assumectleven, assumePhiEven,
checkSignificance, calculateCovariance);
}
else{
get_chebyschev(edges, npol,
filtered_events, nWeightedEvents,
poly_moments,
fejerkernel, assumectleven, assumePhiEven,
checkSignificance, calculateCovariance);
}
spdlog::debug("Getting correct polynomial coefficients in cos(thetak) and cos(thetal)");
spdlog::debug("poly_moments size:{0:d}", poly_moments.size());
//correct coefficients from -1 +1 to ctlmin .. ctlmax for ctl, bit difficult, requires binomial coefficients
std::vector<double> poly_moments_ctl(npol.getSize(), 0.0);
std::vector<double> poly_moments_ctk(npol.getSize(), 0.0);
std::vector<double> poly_moments_phi(npol.getSize(), 0.0);
std::vector<double> poly_moments_q2(npol.getSize(), 0.0);
correct_coefficients(npol,edges,poly_moments_ctl,poly_moments, "ctl");
correct_coefficients(npol,edges,poly_moments_ctk,poly_moments_ctl, "ctk");
correct_coefficients(npol,edges,poly_moments_phi,poly_moments_ctk, "phi");
correct_coefficients(npol,edges,poly_moments_q2, poly_moments_phi, "q2");
spdlog::debug("Getting moments....");
std::vector<double> coeffs_moments(poly_moments_q2);
//calculate norm after full integration
double full_norm = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
full_norm += coeffs_moments.at(npol.get_bin_in4D(h,i,j,k))*integrate_4D(edges, h, i, j, k);
}
}
}
}
//assume a flat efficiency of 1.0 -> scale by edges.range_q2()*RANGE_3D()/full_norm
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
coeffs_moments.at(npol.get_bin_in4D(h,i,j,k)) *= edges.range_q2()*RANGE_3D()/full_norm;
}
}
}
}
//print to std::out
if (spdlog_trace()){
spdlog::trace( "//vector of efficiency coefficients containing {0:d} elements", coeffs_moments.size());
spdlog::trace( "coeffs_eff_4d = { " );
for (unsigned int i=0; i<coeffs_moments.size(); i++){
std::cout << std::setprecision(5) << std::scientific << coeffs_moments.at(i)<< (i<coeffs_moments.size()-1 ? ", " : " ");
}
}
//set this as default coefficients now
coeffs_eff_4d = coeffs_moments;
//calculate the integrated angular terms int f_i dOmega (q2) = xi(q2)
std::vector<TH1D*> hxi_moments;
for (auto &angObs: get_angObser_withTeX_vec()){
hxi_moments.push_back(new TH1D(("hxi"+angObs[0]+"_moments").c_str(), (";q^{2};"+angObs[1]).c_str(),
edges.nq2bins, edges.q2min, edges.q2max));
}
for (int l=1; l<=hxi_moments[0]->GetNbinsX(); l++) {
double a = hxi_moments[0]->GetXaxis()->GetBinLowEdge(l);
double b = hxi_moments[0]->GetXaxis()->GetBinLowEdge(l+1);
std::vector<double> integral_xi(hxi_moments.size(),0.0);
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
unsigned int bin =npol.get_bin_in4D(h, i, j, k);
double coeff = C_PWAVE*coeffs_moments.at(bin)*pow(0.5*(a+b), h);
//f1 = c * sinthetak2; //j1s
integral_xi[0] += coeff
* integrate_x_to_n(edges.ctlmin, edges.ctlmax, i)
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j) - integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
* integrate_x_to_n(edges.phimin, edges.phimax,k);
//f2 = c * costhetak2; //j1c
integral_xi[1] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i)
*integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f3 = c * sinthetak2*cos2thetal; //j2s
integral_xi[2] += coeff
*(2.0*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2)-integrate_x_to_n(edges.ctlmin, edges.ctlmax, i))
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f4 = c * costhetak2*cos2thetal; //j2c
integral_xi[3] += coeff
*(2.0*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2)-integrate_x_to_n(edges.ctlmin, edges.ctlmax, i))
*integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f5 = c * sinthetak2*sinthetal2*cos(2.0*phi); //j3
integral_xi[4] += coeff
*(integrate_x_to_n(edges.ctlmin, edges.ctlmax, i)-integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2))
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n_times_cos_2x(edges.phimin, edges.phimax,k);
//f6 = c * sin2thetak*sin2thetal*cos(phi); //j4
integral_xi[5] += coeff
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i+1)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_cos_x(edges.phimin, edges.phimax,k);
//f7 = c * sin2thetak*sinthetal*cos(phi); //j5
integral_xi[6] += coeff
*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_cos_x(edges.phimin, edges.phimax,k);
//f8 = c * sinthetak2*costhetal; //j6s
integral_xi[7] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+1)
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f9 = c * costhetak2*costhetal; //j6c
integral_xi[8] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+1)
*integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f10 = c * sin2thetak*sinthetal*sin(phi); //j7
integral_xi[9] += coeff
*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_sin_x(edges.phimin, edges.phimax,k);
//f11 = c * sin2thetak*sin2thetal*sin(phi); //j8
integral_xi[10] += coeff
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i+1)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_sin_x(edges.phimin, edges.phimax,k);
//f12 = c * sinthetak2*sinthetal2*sin(2.0*phi); //j9
integral_xi[11] += coeff
*(integrate_x_to_n(edges.ctlmin, edges.ctlmax, i) - integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2))
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j) - integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n_times_sin_2x(edges.phimin, edges.phimax,k);
}
}
}
}
for_indexed(auto hist: hxi_moments) hist->SetBinContent(l, is_param_folded(i,opts) ? 0.0 : integral_xi[i]);
}
//-----------------------------------------------------
//DETERMINE POLYNOMIAL ACCEPTANCE USING UNBINNED ML FIT
//-----------------------------------------------------
std::vector<TH1D*> hxi_ml_hist;
for (auto &angObs: get_angObser_withTeX_vec()){
hxi_ml_hist.push_back(new TH1D(("hxi"+angObs[0]+"_ml").c_str(), (";q^{2};"+angObs[1]).c_str(),
nBinsQ2, edges.q2min, edges.q2max));
}
for (auto &hist: hxi_ml_hist){
hist->SetLineColor(3);
}
//possible TODO for when you're bored: move into it's own function
//Running minuit is expensive, use only if necesary!
if (runMinuit){
spdlog::info("Running Minuit....");
//setup minuit
TMinuit* minuit_eff = new TMinuit(10000);
int errorcode;
minuit_eff->SetFCN(&(eff_fcn_phsp_4d));
double strategy_list[1] = {2.0};
minuit_eff->mnexcm("SET STR", strategy_list, 1, errorcode);
minuit_eff->SetMaxIterations(1000000);
minuit_eff->SetErrorDef(1.0);
const bool second_run = false;
//define params (could also feed with starting values from method of moments)
unsigned int q2_limit = 6;
unsigned int ctl_limit = 3;
unsigned int ctk_limit = 4;
unsigned int phi_limit = 1;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
std::ostringstream sout;
sout << "c" << h << i << j << k;
if (h==0 && i==0 && j==0 && k==0){
minuit_eff->DefineParameter(npol.get_bin_in4D(h,i,j,k), sout.str().c_str(), 1.0, 0.0, 0.0, 0.0);
}
else if (h >= q2_limit || i >= ctl_limit || j >= ctk_limit || k >= phi_limit){
//this is to be able to do higher orders with the moments method and lower ones with minuit
minuit_eff->DefineParameter(npol.get_bin_in4D(h,i,j,k), sout.str().c_str(), 0.0, 0.0, 0.0, 0.0);
}
else if ((second_run && i == 0 && k == 0)|| (!second_run && i % 2 == 0 && k % 2 == 0)){
minuit_eff->DefineParameter(npol.get_bin_in4D(h,i,j,k), sout.str().c_str(), 0.0, 0.001, 0.0, 0.0);
}
else{
minuit_eff->DefineParameter(npol.get_bin_in4D(h,i,j,k), sout.str().c_str(), 0.0, 0.0, 0.0, 0.0);
}
}
}
}
}
double migrad_options[2] = {100000000.0, 0.1}; //{max calls, tolerance}
spdlog::info( "Running Migrad to determine eff shape" );
minuit_eff->mnset();
minuit_eff->mnexcm("MIG", migrad_options, 2, errorcode);
if (second_run){
spdlog::info( "Running Migrad a second time releasing parameters" );
//second run releasing some parameters
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
std::ostringstream sout;
sout << "c" << h << i << j << k;
if ((i > 0 || k > 0) && i % 2 == 0 && k % 2 == 0
&& h < q2_limit && i < ctl_limit && j < ctk_limit && k < phi_limit
)// && (i+j+k)<=4)
minuit_eff->DefineParameter(npol.get_bin_in4D(h,i,j,k), sout.str().c_str(), 0.0, 0.001, 0.0, 0.0);
}
spdlog::info( "Running Migrad to determine eff shape" );
minuit_eff->mnexcm("MIG", migrad_options, 2, errorcode);
}
//extract fitted parameters
spdlog::debug("Extracting fitted parameters.");
std::vector<double> chijk_ml;
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double v,e;
minuit_eff->GetParameter(npol.get_bin_in4D(h,i,j,k), v, e);
chijk_ml.push_back(v);
}
delete minuit_eff;
spdlog::info( "Minuit finished" );
//transform chebyshev coefficients to monomial coefficients
std::vector<double> poly_ml(npol.getSize(), 0.0);
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double c = chijk_ml.at(npol.get_bin_in4D(h,i,j,k));
std::vector<double> ch_coeffs_q2(npol.q2, 0.0);
std::vector<double> ch_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> ch_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> ch_coeffs_phi(npol.phi,0.0);
ch_coeffs_q2.at(h) = 1.0;
ch_coeffs_ctl.at(i) = 1.0;
ch_coeffs_ctk.at(j) = 1.0;
ch_coeffs_phi.at(k) = 1.0;
std::vector<double> poly_coeffs_q2(npol.q2, 0.0);
std::vector<double> poly_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> poly_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> poly_coeffs_phi(npol.phi,0.0);
chebyshev_to_poly(ch_coeffs_q2, poly_coeffs_q2);
chebyshev_to_poly(ch_coeffs_ctl, poly_coeffs_ctl);
chebyshev_to_poly(ch_coeffs_ctk, poly_coeffs_ctk);
chebyshev_to_poly(ch_coeffs_phi, poly_coeffs_phi);
for (unsigned int l = 0; l<=h; l++)
for (unsigned int m = 0; m<=i; m++)
for (unsigned int n = 0; n<=j; n++)
for (unsigned int o = 0; o<=k; o++){
unsigned int bin = npol.get_bin_in4D(l,m,n,o);
poly_ml.at(bin) += c*poly_coeffs_q2.at(l)*poly_coeffs_ctl.at(m)*poly_coeffs_ctk.at(n)*poly_coeffs_phi.at(o);
}
}//tested, works
//correct coefficients from -1 +1 to ctlmin .. ctlmax for ctl, bit difficult, requires binomial coefficients
std::vector<double> poly_ml_ctl(npol.getSize(), 0.0);
std::vector<double> poly_ml_ctk(npol.getSize(), 0.0);
std::vector<double> poly_ml_phi(npol.getSize(), 0.0);
std::vector<double> poly_ml_q2(npol.getSize(), 0.0);
correct_coefficients(npol, edges, poly_ml_ctl, poly_ml, "ctl");
correct_coefficients(npol, edges, poly_ml_ctk, poly_ml_ctl, "ctk");
correct_coefficients(npol, edges, poly_ml_phi, poly_ml_ctk, "phi");
correct_coefficients(npol, edges, poly_ml_q2, poly_ml_phi, "q2");
std::vector<double> coeffs_ml(poly_ml_q2);
//calculate norm after full integration
double full_norm = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
full_norm += coeffs_ml.at(npol.get_bin_in4D(h,i,j,k))*integrate_4D(edges,h,i,j,k);
}
}
}
}
//assume a flat efficiency of 1.0 -> scale by edges.range_q2()*RANGE_3D()/full_norm
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
coeffs_ml.at(npol.get_bin_in4D(h,i,j,k)) *= edges.range_q2()*RANGE_3D()/full_norm;
// delta_q2*delta_ctl*delta_ctk*delta_phi/full_norm;
}
}
}
}
//calculate the integrated angular terms int f_i dOmega (q2) = xi(q2)
for (int l=1; l<=hxi_ml_hist[0]->GetNbinsX(); l++){
double a = hxi_ml_hist[0]->GetXaxis()->GetBinLowEdge(l);
double b = hxi_ml_hist[0]->GetXaxis()->GetBinLowEdge(l+1);
std::vector<double> integral_xi(hxi_ml_hist.size(), 0.0);
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
unsigned int bin = npol.get_bin_in4D(h,i,j,k);
double coeff = C_PWAVE*coeffs_ml.at(bin)*pow(0.5*(a+b), h);
//f1 = c * sinthetak2; //j1s
integral_xi[0] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i)
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f2 = c * costhetak2; //j1c
integral_xi[1] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i)
*integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f3 = c * sinthetak2*cos2thetal; //j2s
integral_xi[2] += coeff
*(2.0*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2)-integrate_x_to_n(edges.ctlmin, edges.ctlmax, i))
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f4 = c * costhetak2*cos2thetal; //j2c
integral_xi[3] += coeff
*(2.0*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2)-integrate_x_to_n(edges.ctlmin, edges.ctlmax, i))
*integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f5 = c * sinthetak2*sinthetal2*cos(2.0*phi); //j3
integral_xi[4] += coeff
*(integrate_x_to_n(edges.ctlmin, edges.ctlmax, i)-integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2))
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n_times_cos_2x(edges.phimin, edges.phimax,k);
//f6 = c * sin2thetak*sin2thetal*cos(phi); //j4
integral_xi[5] += coeff
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i+1)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_cos_x(edges.phimin, edges.phimax,k);
//f7 = c * sin2thetak*sinthetal*cos(phi); //j5
integral_xi[6] += coeff
*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_cos_x(edges.phimin, edges.phimax,k);
//f8 = c * sinthetak2*costhetal; //j6s
integral_xi[7] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+1)
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j)-integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f9 = c * costhetak2*costhetal; //j6c
integral_xi[8] += coeff
*integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+1)
*integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2)
*integrate_x_to_n(edges.phimin, edges.phimax,k);
//f10 = c * sin2thetak*sinthetal*sin(phi); //j7
integral_xi[9] += coeff
*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_sin_x(edges.phimin, edges.phimax,k);
//f11 = c * sin2thetak*sin2thetal*sin(phi); //j8
integral_xi[10] += coeff
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctlmin, edges.ctlmax, i+1)
*2.0*integrate_x_to_n_times_sqrt_1_minus_x2(edges.ctkmin, edges.ctkmax, j+1)
*integrate_x_to_n_times_sin_x(edges.phimin, edges.phimax,k);
//f12 = c * sinthetak2*sinthetal2*sin(2.0*phi); //j9
integral_xi[11] += coeff
*(integrate_x_to_n(edges.ctlmin, edges.ctlmax, i) - integrate_x_to_n(edges.ctlmin, edges.ctlmax, i+2))
*(integrate_x_to_n(edges.ctkmin, edges.ctkmax, j) - integrate_x_to_n(edges.ctkmin, edges.ctkmax, j+2))
*integrate_x_to_n_times_sin_2x(edges.phimin, edges.phimax,k);
}
}
}
}
for_indexed(auto hist: hxi_ml_hist) hist->SetBinContent(l, is_param_folded(i,opts) ? 0.0 : integral_xi[i]);
}
}
//-----------------------
// Factorization in q2 bins
//-----------------------
//std::vector<double> nfact();
if (checkFactorization){ //possible TODO for when you're bored: move into it's own function
spdlog::info( "CHECKFACTORIZATION start " );
//fill histo for comparison
std::vector<double> nsel_hist(edges.getNbins(), 0.0);
std::vector<double> nsel_hist_errsq(edges.getNbins(), 0.0);
for (auto &meas: filtered_events){
unsigned int h = (meas.q2 - edges.q2min ) / edges.binWidth_q2();
unsigned int i = (meas.costhetal - edges.ctlmin ) / edges.binWidth_ctl();
unsigned int j = (meas.costhetak - edges.ctkmin ) / edges.binWidth_ctk();
unsigned int k = (meas.phi - edges.phimin ) / edges.binWidth_phi();
unsigned int bin = get_bin_in4D(h,i,j,k,edges);
nsel_hist.at(bin) += meas.weight;
nsel_hist_errsq.at(bin) += sqr(meas.weight);
}
for (unsigned int i=0; i<nsel_hist_errsq.size(); i++){
if (nsel_hist_errsq.at(i) == 0.0) nsel_hist_errsq.at(i) = 1.0;
}
//loop over all q2 bins
double chi2 = 0.0;
for (unsigned int h=0; h<q2_edges.size()-1; h++){
double q2a = q2_edges.at(h);
double q2b = q2_edges.at(h+1);
std::vector<double> ctl_moments(npol.ctl, 0.0);
std::vector<double> ctk_moments(npol.ctk, 0.0);
std::vector<double> phi_moments(npol.phi,0.0);
double nfactweighted = 0.0;
//determine weighted number of events in this bin and fill control histos
std::ostringstream sout;
sout << h;
sout << "_" << DECAY_NAME << (opts->KS ? opts->DDLL+"_" : "_")
<< "Run" << opts->run
<< ((opts->angacccorrperyear && opts->year != -1) ? "_"+std::to_string(opts->year) : "")
<< "_" << tagNumber;
TH1D* hctl = new TH1D("hctl", ";cos#Theta_{L};", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* hctk = new TH1D("hctk", ";cos#Theta_{K};", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* hphi = new TH1D("hphi", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
for (auto &meas : filtered_events){
if (meas.q2 > q2a && meas.q2 < q2b){
nfactweighted += meas.weight;
hctl->Fill(meas.costhetal, meas.weight);
hctk->Fill(meas.costhetak, meas.weight);
hphi->Fill(meas.phi, meas.weight);
}
}
//loop over events in this bin, determine factorizing moments
for (auto & meas: filtered_events){
if (meas.q2 > q2a && meas.q2 < q2b){
double weight = meas.weight;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin ) / edges.range_ctl() - 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin ) / edges.range_ctk() - 1.0;
double phi = 2.0 * (meas.phi - edges.phimin ) / edges.range_phi() - 1.0;
for (unsigned int i = 0; i<npol.ctl; i++){
if (assumectleven && i%2 != 0){
ctl_moments.at(i) = 0.0;
}
else{
//ctl_moments.at(i)+=weight*fcnc::chebyshev(ctl,i)/sqrt(1-ctl*ctl)*(i==0?1.0/pi:2.0/pi)/nfactweighted;
ctl_moments.at(i) += weight*fcnc::legendre(ctl, i)*(2.0*i+1.0)/2.0/nfactweighted;
}
}
for (unsigned int j = 0; j<npol.ctk; j++){
//ctk_moments.at(j) += weight*fcnc::chebyshev(ctk, j)/sqrt(1-ctk*ctk)*(j == 0 ? 1.0/pi : 2.0/pi)/nfactweighted;
ctk_moments.at(j) += weight*fcnc::legendre(ctk, j)*(2.0*j+1.0)/2.0/nfactweighted;
}
for (unsigned int k = 0; k<npol.phi; k++){
if (assumePhiEven && k%2 != 0){
phi_moments.at(k) = 0.0;
}
else{
//phi_moments.at(k) += weight*fcnc::chebyshev(phi, k)/sqrt(1-phi*phi)*(k == 0 ? 1.0/pi : 2.0/pi)/nfactweighted;
phi_moments.at(k) += weight*fcnc::legendre(phi, k)*(2.0*k+1.0)/2.0/nfactweighted;
}
}
}
}
//make polynomial
std::vector<double> ctl_polys(npol.ctl, 0.0);
std::vector<double> ctk_polys(npol.ctk, 0.0);
std::vector<double> phi_polys(npol.phi,0.0);
if (doLegendre){
legendre_to_poly(ctl_moments, ctl_polys);
legendre_to_poly(ctk_moments, ctk_polys);
legendre_to_poly(phi_moments, phi_polys);
}
else{
chebyshev_to_poly(ctl_moments, ctl_polys);
chebyshev_to_poly(ctk_moments, ctk_polys);
chebyshev_to_poly(phi_moments, phi_polys);
}
//transform angles from -1 +1 to angle_min angle_max
std::vector<double> ctl_correct(npol.ctl, 0.0);
std::vector<double> ctk_correct(npol.ctk, 0.0);
std::vector<double> phi_correct(npol.phi,0.0);
for (unsigned int i=0; i<ctk_polys.size(); i++) ctk_correct.at(i) = ctk_polys.at(i);
for (unsigned int i=0; i<ctl_polys.size(); i++) ctl_correct.at(i) = ctl_polys.at(i);
correct_poly(phi_polys, phi_correct, -TMath::Pi(), +TMath::Pi());
//determine chi2 in this bin
std::vector<double> neps_hist(nctlbins*nctkbins*nphibins, 0.0);
for (unsigned int i = 0; i<nctlbins; i++){
double ctla = edges.ctlmin + edges.range_ctl() * i /double(edges.nctlbins);
double ctlb = edges.ctlmin + edges.range_ctl() * (i+1.0) /double(edges.nctlbins);
for (unsigned int j = 0; j<nctkbins; j++){
double ctka = edges.ctkmin + edges.range_ctk() * j /double(edges.nctkbins);
double ctkb = edges.ctkmin + edges.range_ctk() * (j+1.0) /double(edges.nctkbins);
for (unsigned int k = 0; k<nphibins; k++){
double phi_a = edges.phimin + edges.range_phi() * k / double(edges.nphibins);
double phi_b = edges.phimin + edges.range_phi() * (k+1.0) / double(edges.nphibins);
double result_ctl = 0.0, result_ctk = 0.0, result_phi = 0.0;
for (unsigned int m = 0; m<npol.ctl; m++) result_ctl += ctl_correct.at(m)*integrate_x_to_n(ctla,ctlb,m);
for (unsigned int n = 0; n<npol.ctk; n++) result_ctk += ctk_correct.at(n)*integrate_x_to_n(ctka,ctkb,n);
for (unsigned int o = 0; o<npol.phi; o++) result_phi += phi_correct.at(o)*integrate_x_to_n(phi_a,phi_b,o);
neps_hist.at(get_bin_in3D(i,j,k,edges)) = result_ctl*result_ctk*result_phi;
}
}
}
//rescale to be sure
double nnorm_ctl = 0.0, nnorm_ctk = 0.0, nnorm_phi = 0.0;
for (unsigned int m = 0; m<npol.ctl; m++) nnorm_ctl += ctl_correct.at(m)*integrate_x_to_n(edges.ctlmin, edges.ctlmax, m);
for (unsigned int n = 0; n<npol.ctk; n++) nnorm_ctk += ctk_correct.at(n)*integrate_x_to_n(edges.ctkmin, edges.ctkmax, n);
for (unsigned int o = 0; o<npol.phi; o++) nnorm_phi += phi_correct.at(o)*integrate_x_to_n(edges.phimin, edges.phimax, o);
for (unsigned int i = 0; i<nctlbins; i++){
for (unsigned int j = 0; j<nctkbins; j++){
for (unsigned int k = 0; k<nphibins; k++){
neps_hist.at(get_bin_in3D(i,j,k,edges)) *= (nfactweighted/nnorm_ctl/nnorm_ctk/nnorm_phi);
}
}
}
//should scale nfactweighted/nfactfit
//determine chi2 //TODO: not working?
double curr_chi2 = 0.0;
TH1D* hpull = new TH1D("hpull", ";(N_{sel}-N_{fit})/#sqrt{N_{sel}};", 100, -5.0, 5.0);
for (unsigned int i = 0; i<nctlbins; i++){
for (unsigned int j = 0; j<nctkbins; j++){
for (unsigned int k = 0; k<nphibins; k++){
unsigned int bin = get_bin_in4D(h,i,j,k,edges);
curr_chi2 += sqr(nsel_hist.at(bin)-neps_hist.at(get_bin_in3D(i,j,k,edges)))/fabs(neps_hist.at(get_bin_in3D(i,j,k,edges))); //TODO: h is a wrong bin!
hpull->Fill((nsel_hist.at(bin)-neps_hist.at(get_bin_in3D(i,j,k,edges)))/sqrt(fabs(neps_hist.at(get_bin_in3D(i,j,k,edges)))));//TODO: h is a wrong bin!
}
}
}
//do plots
TH1D* ectl = new TH1D("ectl", ";cos#Theta_{L};", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* ectk = new TH1D("ectk", ";cos#Theta_{K};", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* ephi = new TH1D("ephi", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
for (int l=1; l<=ectl->GetNbinsX(); l++){
double integral = 0.0;
for (unsigned int i = 0; i<npol.ctl; i++){
integral += ctl_correct.at(i)*integrate_x_to_n(ectl->GetBinLowEdge(l),ectl->GetBinLowEdge(l+1),i);
}
ectl->SetBinContent(l, integral);
}
for (int l=1; l<=ectk->GetNbinsX(); l++){
double integral = 0.0;
for (unsigned int j = 0; j<npol.ctk; j++){
integral += ctk_correct.at(j)*integrate_x_to_n(ectl->GetBinLowEdge(l),ectl->GetBinLowEdge(l+1),j);
}
ectk->SetBinContent(l, integral);
}
for (int l=1; l<=ephi->GetNbinsX(); l++){
double integral = 0.0;
for (unsigned int k = 0; k<npol.phi; k++){
integral += phi_correct.at(k)*integrate_x_to_n(ectl->GetBinLowEdge(l),ectl->GetBinLowEdge(l+1),k);
}
ephi->SetBinContent(l, integral);
}
TLatex *tex = getPrettyTex(0.06,31);
std::ostringstream sregion;
sregion << std::fixed << std::setprecision(2) << q2a << "<q^{2}<" << q2b;
TCanvas* c_pull = new TCanvas("c_pull", "c_pull", 1600, 1200);
c_pull->cd();
hpull->Draw();
tex->DrawLatex(0.85, 0.15, sregion.str().c_str());
if(opts->write_eps)c_pull->Print((std::string(subfolder + "fact_pull")+sout.str()+std::string(".eps")).c_str(),"eps");
plotAngular(hctl,ectl,opts->write_eps,"fact_ctleff",sout.str(),sregion.str(),tex,subfolder);
plotAngular(hctk,ectk,opts->write_eps,"fact_ctkeff",sout.str(),sregion.str(),tex,subfolder);
plotAngular(hphi,ephi,opts->write_eps,"fact_phieff",sout.str(),sregion.str(),tex,subfolder);
delete hpull;
delete ectl;
delete ectk;
delete ephi;
delete hctl;
delete hctk;
delete hphi;
delete c_pull;
spdlog::debug( "fact chi2 in bin {0:d}: {1:f}", h, curr_chi2 );
chi2 += curr_chi2;
}
spdlog::info( "fact total chi2: {0:f}",chi2 );
spdlog::info( "uses {0:d}*{1:d}*{2:d}*{3:d}={4:d} bins",edges.nq2bins, edges.nctlbins, edges.nctkbins, edges.nphibins, edges.getNbins());
spdlog::info( "Chi2/Nbins = {0:f}", chi2/(edges.getNbins()) );
spdlog::info( "CHECKFACTORIZATION end " );
}
//-----------------------
// 3D MOMENTS in q2 bins
//-----------------------
if (do3Dmoments){
spdlog::info("DOTHREEDMOMENTS start ");
//fill histo for comparison
std::vector<double> nsel_hist(edges.getNbins(), 0.0);
std::vector<double> nsel_hist_errsq(edges.getNbins(), 0.0);
for (unsigned int e=0; e<filtered_events.size(); e++){
const event & meas = filtered_events.at(e);
unsigned int bin = get_bin_in4D(meas,edges);
nsel_hist.at(bin) += meas.weight;
nsel_hist_errsq.at(bin) += sqr(meas.weight);
}
for (unsigned int i=0; i<nsel_hist_errsq.size(); i++){
if (nsel_hist_errsq.at(i) == 0.0) nsel_hist_errsq.at(i) = 1.0;
}
//loop over all q2 bins
double chi2 = 0.0;
for (unsigned int h=0; h<q2_edges.size()-1; h++){
double q2a = q2_edges.at(h);
double q2b = q2_edges.at(h+1);
double nthreedweighted = 0.0;
//determine weighted number of events in this bin and fill control histos
std::ostringstream sout;
sout << h;
sout << "_" << DECAY_NAME << (opts->KS ? opts->DDLL+"_" : "_")
<< "Run" << opts->run
<< ((opts->angacccorrperyear && opts->year != -1) ? "_"+std::to_string(opts->year) : "");
TH1D* hctl = new TH1D("hctl", ";cos#Theta_{L};", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* hctk = new TH1D("hctk", ";cos#Theta_{K};", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* hphi = new TH1D("hphi", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
for (auto &meas: filtered_events){
if (meas.q2 > q2a && meas.q2 < q2b){
nthreedweighted += meas.weight;
hctl->Fill(meas.costhetal, meas.weight);
hctk->Fill(meas.costhetak, meas.weight);
hphi->Fill(meas.phi, meas.weight);
}
}
//loop over events in this bin, determine factorizing moments
spdlog::debug("Getting factorizing moments");
std::vector<double> threed_moments(npol.ctl*npol.ctk*npol.phi,0.0);
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
for (auto &meas: filtered_events){
if (meas.q2 > q2a && meas.q2 < q2b){
double weight = meas.weight;
double ctl = 2.0 * (meas.costhetal - edges.ctlmin ) / edges.range_ctl() - 1.0;
double ctk = 2.0 * (meas.costhetak - edges.ctkmin ) / edges.range_ctk() - 1.0;
double phi = 2.0 * (meas.phi - edges.phimin ) / edges.range_phi() - 1.0;
double result = 0.0;
if (assumectleven && i%2 != 0) result = 0.0;
else if (assumePhiEven && k%2 != 0) result = 0.0;
else {
if (doLegendre){
result = weight
*fcnc::legendre(ctl, i)*(2.0*i+1.0)/2.0
*fcnc::legendre(ctk, j)*(2.0*j+1.0)/2.0
*fcnc::legendre(phi, k)*(2.0*k+1.0)/2.0
/nthreedweighted;// /filtered_events.size();
}
else{
result = weight
*fcnc::chebyshev(ctl, i)/sqrt(1-ctl*ctl)*(i==0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(ctk, j)/sqrt(1-ctk*ctk)*(j==0 ? 1.0/MY_PI : 2.0/MY_PI)
*fcnc::chebyshev(phi, k)/sqrt(1-phi*phi)*(k==0 ? 1.0/MY_PI : 2.0/MY_PI)
/nthreedweighted;// /filtered_events.size();
}
}
threed_moments.at( get_bin_in3D(i,j,k,npol)) += result;
}
}
}
}
}
//make polynomial
//transform legendre coefficients to monomial coefficients
spdlog::debug("Transforming legendre coefficients to monomial coefficients....");
std::vector<double> threedpoly_moments(npol.ctl*npol.ctk*npol.phi,0.0);
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
double tmp = threed_moments.at(get_bin_in3D(i,j,k,npol));
std::vector<double> leg_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> leg_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> leg_coeffs_phi(npol.phi,0.0);
leg_coeffs_ctl.at(i) = 1.0;
leg_coeffs_ctk.at(j) = 1.0;
leg_coeffs_phi.at(k) = 1.0;
std::vector<double> poly_coeffs_ctl(npol.ctl, 0.0);
std::vector<double> poly_coeffs_ctk(npol.ctk, 0.0);
std::vector<double> poly_coeffs_phi(npol.phi,0.0);
legendre_to_poly(leg_coeffs_ctl, poly_coeffs_ctl);
legendre_to_poly(leg_coeffs_ctk, poly_coeffs_ctk);
legendre_to_poly(leg_coeffs_phi, poly_coeffs_phi);
for (unsigned int m = 0; m<=i; m++){
for (unsigned int n = 0; n<=j; n++){
for (unsigned int o = 0; o<=k; o++){
unsigned int bin = get_bin_in3D(m,n,o,npol);
threedpoly_moments.at(bin) += tmp*poly_coeffs_ctl.at(m)*poly_coeffs_ctk.at(n)*poly_coeffs_phi.at(o);
}
}
}
}
}
}
//correct coefficients from -1 +1 to -pi .. +pi for phi, this is easy
std::vector<double> threedpoly_moments_phi(threedpoly_moments);
for (unsigned int m = 0; m<npol.ctl; m++){
for (unsigned int n = 0; n<npol.ctk; n++){
for (unsigned int o = 0; o<npol.phi; o++){
threedpoly_moments_phi.at(get_bin_in3D(m,n,o,npol)) *= pow(1.0/MY_PI, o);
}
}
}
//calculate norm after full integration
double full_norm = 0.0;
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
double result = threedpoly_moments_phi.at(get_bin_in3D(i,j,k,npol))*integrate_3D(edges,i,j,k);
full_norm += result;
}
}
}
spdlog::debug("Rescaling...");
//rescale to be sure
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
threedpoly_moments_phi.at(get_bin_in3D(i,j,k,npol)) *= nthreedweighted/full_norm;
}
}
}
//determine nfit in this bin
std::vector<double> neps_hist(nctlbins*nctkbins*nphibins, 0.0);
for (unsigned int i = 0; i<nctlbins; i++){
double ctla = edges.ctlmin + i * edges.range_ctl();
double ctlb = edges.ctlmin + (i+1) * edges.range_ctl();
for (unsigned int j = 0; j<nctkbins; j++){
double ctka = edges.ctkmin + j * edges.range_ctk();
double ctkb = edges.ctkmin + (j+1) * edges.range_ctk();
for (unsigned int k = 0; k<nphibins; k++){
double phi_a = edges.phimin + k * edges.range_phi();
double phi_b = edges.phimin + (k+1) * edges.range_phi();
double result = 0.0;
for (unsigned int m = 0; m<npol.ctl; m++)
for (unsigned int n = 0; n<npol.ctk; n++)
for (unsigned int o = 0; o<npol.phi; o++)
result += threedpoly_moments_phi.at(get_bin_in3D(m,n,o,npol))
*(pow(ctlb,m+1)/(m+1.0) - pow(ctla, m+1)/(m+1.0))
*(pow(ctkb,n+1)/(n+1.0) - pow(ctka, n+1)/(n+1.0))
*(pow(phi_b,o+1)/(o+1.0) - pow(phi_a, o+1)/(o+1.0));
unsigned int bin = get_bin_in3D(i,j,k,edges);
neps_hist.at(bin) = result;
}
}
}
spdlog::debug("Getting chi2...");
//determine chi2
double curr_chi2 = 0.0;
TH1D* hpull = new TH1D("hpull", ";(N_{sel}-N_{fit})/#sqrt{N_{sel}};", 100, -5.0, 5.0);
for (unsigned int i = 0; i<nctlbins; i++)
for (unsigned int j = 0; j<nctkbins; j++)
for (unsigned int k = 0; k<nphibins; k++){
unsigned int bin =get_bin_in4D(h,i,j,k,edges);
curr_chi2 += sqr(nsel_hist.at(bin)-neps_hist.at(get_bin_in3D(i,j,k,edges)))/fabs(neps_hist.at(get_bin_in3D(i,j,k,edges)));
hpull->Fill((nsel_hist.at(bin)-neps_hist.at(get_bin_in3D(i,j,k,edges)))/sqrt(fabs(neps_hist.at(get_bin_in3D(i,j,k,edges)))));
}
//spdlog::debug( "ok" );
//do plots
TH1D* ectl = new TH1D("ectl", ";cos#Theta_{L};", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* ectk = new TH1D("ectk", ";cos#Theta_{K};", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* ephi = new TH1D("ephi", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
for (int l=1; l<=ectl->GetNbinsX(); l++){
double a = ectl->GetBinLowEdge(l);
double b = ectl->GetBinLowEdge(l+1);
double integral = 0.0;
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double result = threedpoly_moments_phi.at(get_bin_in3D(i,j,k,npol))
*(pow(b ,i+1)/(i+1.0) - pow(a ,i+1)/(i+1.0))
*(pow(edges.ctkmax,j+1)/(j+1.0) - pow(edges.ctkmin,j+1)/(j+1.0))
*(pow(edges.phimax,k+1)/(k+1.0) - pow(edges.phimin,k+1)/(k+1.0));
integral += result;
}
ectl->SetBinContent(l, integral);
}
for (int l=1; l<=ectk->GetNbinsX(); l++){
double a = ectk->GetBinLowEdge(l);
double b = ectk->GetBinLowEdge(l+1);
double integral = 0.0;
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double result = threedpoly_moments_phi.at(get_bin_in3D(i,j,k,npol))
*(pow(edges.ctlmax,i+1)/(i+1.0) - pow(edges.ctlmin,i+1)/(i+1.0))
*(pow(b ,j+1)/(j+1.0) - pow(a ,j+1)/(j+1.0))
*(pow(edges.phimax,k+1)/(k+1.0) - pow(edges.phimin,k+1)/(k+1.0));
integral += result;
}
ectk->SetBinContent(l, integral);
}
for (int l=1; l<=ephi->GetNbinsX(); l++){
double a = ephi->GetBinLowEdge(l);
double b = ephi->GetBinLowEdge(l+1);
double integral = 0.0;
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++){
double result = threedpoly_moments_phi.at(get_bin_in3D(i,j,k,npol))
*(pow(edges.ctlmax,i+1)/(i+1.0) - pow(edges.ctlmin,i+1)/(i+1.0))
*(pow(edges.ctkmax,j+1)/(j+1.0) - pow(edges.ctkmin,j+1)/(j+1.0))
*(pow(b ,k+1)/(k+1.0) - pow(a ,k+1)/(k+1.0));
integral += result;
}
ephi->SetBinContent(l, integral);
}
spdlog::debug("Plotting...");
TLatex *tex = getPrettyTex(0.06,31);
std::ostringstream sregion;
sregion << std::fixed << std::setprecision(2) << q2a << "<q^{2}<" << q2b;
TCanvas* c_zero = new TCanvas("c_zero", "c_zero", 1600, 1200);
c_zero->cd();
hpull->Draw();
tex->DrawLatex(0.85, 0.15, sregion.str().c_str());
if(opts->write_eps)c_zero->Print((std::string(subfolder + "3d_pull")+sout.str()+std::string(".eps")).c_str(),"eps");
plotAngular(hctl,ectl,opts->write_eps,"3d_ctleff", sout.str(), sregion.str(),tex,subfolder);
plotAngular(hctk,ectk,opts->write_eps,"3d_ctkeff", sout.str(), sregion.str(),tex,subfolder);
plotAngular(hphi,ephi,opts->write_eps,"3d_phieff", sout.str(), sregion.str(),tex,subfolder);
delete hpull;
delete ectl;
delete ectk;
delete ephi;
delete hctl;
delete hctk;
delete hphi;
delete c_zero;
spdlog::debug( "3d chi2 in bin {0:d}: {1:f}", h, curr_chi2 );
chi2 += curr_chi2;
}
spdlog::debug( "3d total chi2: {0:f}",chi2 );
spdlog::debug( "3d uses {0:d}*{1:d}*{2:d}*{3:d}={4_d} bins",edges.nq2bins, edges.nctlbins, edges.nctkbins, edges.nphibins, edges.getNbins());
spdlog::debug( "3d chi2/Nbins = {0:f}", chi2/(edges.getNbins()) );
spdlog::info( "DOTHREEDMOMENTS end " );
}
//-----------------------
// CHI2 TEST FOR MOMENTS
//-----------------------
bool dochi2test = true;
if(opts->only_4_1D_chi2)dochi2test = false;
spdlog::debug("dochi2: "+boolToString(dochi2test));
if (dochi2test) {
get_chi2(edges, npol, filtered_events, coeffs_moments, globalvalues,
assumectleven, assumePhiEven, opts->write_eps, subfolder, appendix);
}
//-----------------------
//MAKE NICE CONTROL PLOTS
//-----------------------
const bool do_plots = opts->write_eps || opts->only_4_1D_chi2;
if (do_plots){
spdlog::info( "Start PLOTS.." );
gROOT->SetStyle("Plain");
TPad foo;
const int NRGBs = 5;
const int NCont = 250;
double stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
double red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
double green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
double blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
gStyle->SetNumberContours(NCont);
gStyle->SetLabelFont(132, "xyz");
gStyle->SetTitleFont(132, "xyz");
gStyle->SetLegendFont(132);
gStyle->SetStatFont(132);
//new histograms binned in q2
//compare the three angular projections for data and moments
for (unsigned int u=0; u<q2_edges.size()-1; u++){
std::ostringstream sout;
sout << u;
std::string binno = sout.str();
TH1D* hctl = new TH1D("hctl", ";cos#Theta_{L};", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* hctk = new TH1D("hctk", ";cos#Theta_{K};", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* hphi = new TH1D("hphi", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctlctk = new TH2D("hctlctk", ";cos#Theta_{L};cos#Theta_{K}", nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.ctkmin, edges.ctkmax);
TH2D* hctlphi = new TH2D("hctlphi", ";cos#Theta_{L};#phi [rad]" , nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctkphi = new TH2D("hctkphi", ";cos#Theta_{K};#phi [rad]" , nBinsAngles, edges.ctkmin, edges.ctkmax, nBinsAngles, edges.phimin, edges.phimax);
double q2a = q2_edges.at(u);
double q2b = q2_edges.at(u+1);
for (auto & meas: filtered_events){
if (meas.q2 > q2a && meas.q2 < q2b){
hctl->Fill(meas.costhetal, meas.weight);
hctk->Fill(meas.costhetak, meas.weight);
hphi->Fill(meas.phi, meas.weight);
hctlctk->Fill(meas.costhetal, meas.costhetak, meas.weight);
hctlphi->Fill(meas.costhetal, meas.phi, meas.weight);
hctkphi->Fill(meas.costhetak, meas.phi, meas.weight);
}
}
TH1D* hctl_moments = new TH1D("hctl_moments", ";cos#Theta_{L};", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* hctk_moments = new TH1D("hctk_moments", ";cos#Theta_{K};", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* hphi_moments = new TH1D("hphi_moments", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctlctk_moments = new TH2D("hctlctk_moments", ";cos#Theta_{L};cos#Theta_{K}", nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.ctkmin, edges.ctkmax);
TH2D* hctlphi_moments = new TH2D("hctlphi_moments", ";cos#Theta_{L};#phi [rad]" , nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctkphi_moments = new TH2D("hctkphi_moments", ";cos#Theta_{K};#phi [rad]" , nBinsAngles, edges.ctkmin, edges.ctkmax, nBinsAngles, edges.phimin, edges.phimax);
//1d plots
fill_moments_hist_1D(hctl_moments, coeffs_moments, edges, npol, q2a, q2b);
fill_moments_hist_1D(hctk_moments, coeffs_moments, edges, npol, q2a, q2b);
fill_moments_hist_1D(hphi_moments, coeffs_moments, edges, npol, q2a, q2b);
//2d plots
fill_moments_hist_2D(hctlctk_moments, coeffs_moments, edges, npol, q2a, q2b);
fill_moments_hist_2D(hctlphi_moments, coeffs_moments, edges, npol, q2a, q2b);
fill_moments_hist_2D(hctkphi_moments, coeffs_moments, edges, npol, q2a, q2b);
TLatex *tex = getPrettyTex(0.06,31);
std::ostringstream sregion;
sregion << std::fixed << std::setprecision(2) << q2a << "<q^{2}<" << q2b;
plotAngular(hctl,hctl_moments,opts->write_eps,"ctleff",binno+appendix,sregion.str(), tex, subfolder);
plotAngular(hctk,hctk_moments,opts->write_eps,"ctkeff",binno+appendix,sregion.str(), tex, subfolder);
plotAngular(hphi,hphi_moments,opts->write_eps,"phieff",binno+appendix,sregion.str(), tex, subfolder);
if (tagNumber < 0){ //Do not plot when testing quickly
plotAngular(hctlctk,hctlctk_moments,opts->write_eps,"ctlctkeff",binno+appendix,
sregion.str(), tex, subfolder);
plotAngular(hctlphi,hctlphi_moments,opts->write_eps,"ctlphieff",binno+appendix,
sregion.str(), tex, subfolder);
plotAngular(hctkphi,hctkphi_moments,opts->write_eps,"ctkphieff",binno+appendix,
sregion.str(), tex, subfolder);
}
delete hctl; delete hctl_moments;
delete hctk; delete hctk_moments;
delete hphi; delete hphi_moments;
delete hctlctk; delete hctlctk_moments;
delete hctkphi; delete hctkphi_moments;
delete hctlphi; delete hctlphi_moments;
}//end extra plots in bins of q2
//data histograms
TH1D* hq2 = new TH1D("hq2" , ";q^{2} [GeV^{2}];" , nBinsQ2 ,edges.q2min, edges.q2max);
TH1D* hctl = new TH1D("hctl", ";cos(#Theta_{L});", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* hctk = new TH1D("hctk", ";cos(#Theta_{K});", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* hphi = new TH1D("hphi", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctlctk = new TH2D("hctlctk", ";cos#Theta_{L};cos#Theta_{K}", nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.ctkmin, edges.ctkmax);
TH2D* hctlphi = new TH2D("hctlphi", ";cos#Theta_{L};#phi [rad]" , nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctkphi = new TH2D("hctkphi", ";cos#Theta_{K};#phi [rad]" , nBinsAngles, edges.ctkmin, edges.ctkmax, nBinsAngles, edges.phimin, edges.phimax);
TH2D* hq2ctl = new TH2D("hq2ctl", ";q^{2};cos#Theta_{L}", nBinsQ2, edges.q2min, edges.q2max, nBinsAngles, edges.ctlmin, edges.ctlmax);
TH2D* hq2ctk = new TH2D("hq2ctk", ";q^{2};cos#Theta_{K}", nBinsQ2, edges.q2min, edges.q2max, nBinsAngles, edges.ctkmin, edges.ctkmax);
TH2D* hq2phi = new TH2D("hq2phi", ";q^{2};#phi [rad]" , nBinsQ2, edges.q2min, edges.q2max, nBinsAngles, edges.phimin, edges.phimax);
hq2->Sumw2();
hctl->Sumw2();
hctk->Sumw2();
hphi->Sumw2();
for (unsigned int i=0; i<filtered_events.size(); i++){
const event& meas = filtered_events.at(i);
hq2->Fill(meas.q2, meas.weight);
hctl->Fill(meas.costhetal, meas.weight);
hctk->Fill(meas.costhetak, meas.weight);
hphi->Fill(meas.phi, meas.weight);
hctlctk->Fill(meas.costhetal, meas.costhetak, meas.weight);
hctlphi->Fill(meas.costhetal, meas.phi, meas.weight);
hctkphi->Fill(meas.costhetak, meas.phi, meas.weight);
hq2ctl->Fill(meas.q2, meas.costhetal, meas.weight);
hq2ctk->Fill(meas.q2, meas.costhetak, meas.weight);
hq2phi->Fill(meas.q2, meas.phi, meas.weight);
}
//fitted histograms, moments method
TH1D* hq2_moments = new TH1D("hq2_moments" , ";q^{2} [GeV^{2}];" , nBinsQ2, edges.q2min, edges.q2max);
TH1D* hctl_moments = new TH1D("hctl_moments", ";cos(#Theta_{L});", nBinsAngles, edges.ctlmin, edges.ctlmax);
TH1D* hctk_moments = new TH1D("hctk_moments", ";cos(#Theta_{K});", nBinsAngles, edges.ctkmin, edges.ctkmax);
TH1D* hphi_moments = new TH1D("hphi_moments", ";#phi [rad];" , nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctlctk_moments = new TH2D("hctlctk_moments", ";cos#Theta_{L};cos#Theta_{K}", nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.ctkmin, edges.ctkmax);
TH2D* hctlphi_moments = new TH2D("hctlphi_moments", ";cos#Theta_{L};#phi [rad]" , nBinsAngles, edges.ctlmin, edges.ctlmax, nBinsAngles, edges.phimin, edges.phimax);
TH2D* hctkphi_moments = new TH2D("hctkphi_moments", ";cos#Theta_{K};#phi [rad]" , nBinsAngles, edges.ctkmin, edges.ctkmax, nBinsAngles, edges.phimin, edges.phimax);
TH2D* hq2ctl_moments = new TH2D("hq2ctl_moments", ";q^{2};cos#Theta_{L}", nBinsQ2, edges.q2min, edges.q2max, nBinsAngles, edges.ctlmin, edges.ctlmax);
TH2D* hq2ctk_moments = new TH2D("hq2ctk_moments", ";q^{2};cos#Theta_{K}", nBinsQ2, edges.q2min, edges.q2max, nBinsAngles, edges.ctkmin, edges.ctkmax);
TH2D* hq2phi_moments = new TH2D("hq2phi_moments", ";q^{2};#phi [rad]" , nBinsQ2, edges.q2min, edges.q2max, nBinsAngles, edges.phimin, edges.phimax);
//1d plots
fill_moments_hist_1D(hq2_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_1D(hctl_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_1D(hctk_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_1D(hphi_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
//2d plots
fill_moments_hist_2D(hctlctk_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_2D(hctlphi_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_2D(hctkphi_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_2D(hq2ctl_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_2D(hq2ctk_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
fill_moments_hist_2D(hq2phi_moments,coeffs_moments,edges,npol,edges.q2min,edges.q2max);
if(opts->only_4_1D_chi2){
//reset globalvalues
globalvalues->chi2_4_1D = 0.0;
globalvalues->chi2_4_1D_ctk = 0.0;
globalvalues->chi2_4_1D_ctl = 0.0;
globalvalues->chi2_4_1D_phi = 0.0;
globalvalues->chi2_4_1D_q2 = 0.0;
globalvalues->ndof_4_1D = 0;
globalvalues->ndof_4_1D_ctk = 0;
globalvalues->ndof_4_1D_ctl = 0;
globalvalues->ndof_4_1D_phi = 0;
globalvalues->ndof_4_1D_q2 = 0;
//first scale moments to match integral of event histograms
hctk_moments->Scale(hctk->Integral()/hctk_moments->Integral());
hctl_moments->Scale(hctl->Integral()/hctl_moments->Integral());
hphi_moments->Scale(hphi->Integral()/hphi_moments->Integral());
hq2_moments ->Scale(hq2->Integral()/hq2_moments->Integral());
//get chi2 from all 4 variables
for(UInt_t b = 1; b <= nBinsAngles; b++){
globalvalues->chi2_4_1D_ctk += sqr((hctk->GetBinContent(b) - hctk_moments->GetBinContent(b))) / fabs(hctk_moments->GetBinContent(b));
globalvalues->chi2_4_1D_ctl += sqr((hctl->GetBinContent(b) - hctl_moments->GetBinContent(b))) / fabs(hctl_moments->GetBinContent(b));
globalvalues->chi2_4_1D_phi += sqr((hphi->GetBinContent(b) - hphi_moments->GetBinContent(b))) / fabs(hphi_moments->GetBinContent(b));
}
for(UInt_t b = 1; b <= nBinsQ2; b++){
globalvalues->chi2_4_1D_q2 += sqr((hq2->GetBinContent(b) - hq2_moments->GetBinContent(b))) / fabs(hq2_moments->GetBinContent(b));
}
//sum of all four 1D histograms
globalvalues->chi2_4_1D = globalvalues->chi2_4_1D_ctk + globalvalues->chi2_4_1D_ctl + globalvalues->chi2_4_1D_phi + globalvalues->chi2_4_1D_q2;
//determine ndof for all angles and q2
//NOT sure if these ndofs are 100% correct, but good enough for a relative measure
globalvalues->ndof_4_1D_ctk = nBinsAngles - (npol.ctk-1) - 1;
globalvalues->ndof_4_1D_ctl = nBinsAngles - (npol.ctl-1) - 1;
globalvalues->ndof_4_1D_phi = nBinsAngles - (npol.phi-1) - 1;
globalvalues->ndof_4_1D_q2 = nBinsQ2 - (npol.q2-1) - 1;
globalvalues->ndof_4_1D = 3*nBinsAngles + nBinsQ2 - 1 - (npol.ctk + npol.ctl + npol.phi + npol.q2 - 1);
}
//plot
plotAngular(hq2,hq2_moments,opts->write_eps,"q2eff",appendix, subfolder);
plotAngular(hctl,hctl_moments,opts->write_eps,"ctleff",appendix, subfolder);
plotAngular(hctk,hctk_moments,opts->write_eps,"ctkeff",appendix, subfolder);
plotAngular(hphi,hphi_moments,opts->write_eps,"phieff",appendix, subfolder);
if (tagNumber < 0){ //do not plot when testing
plotAngular(hctlctk,hctlctk_moments,opts->write_eps,"ctlctkeff",appendix, subfolder);
plotAngular(hctlphi,hctlphi_moments,opts->write_eps,"ctlphieff",appendix, subfolder);
plotAngular(hctkphi,hctkphi_moments,opts->write_eps,"ctkphieff",appendix, subfolder);
plotAngular(hq2ctl,hq2ctl_moments,opts->write_eps,"q2ctleff",appendix, subfolder);
plotAngular(hq2ctk,hq2ctk_moments,opts->write_eps,"q2ctkeff",appendix, subfolder);
plotAngular(hq2phi,hq2phi_moments,opts->write_eps,"q2phieff",appendix, subfolder);
}
delete hq2; delete hq2_moments;
delete hctl; delete hctl_moments;
delete hctk; delete hctk_moments;
delete hphi; delete hphi_moments;
delete hctlctk; delete hctlctk_moments;
delete hctkphi; delete hctkphi_moments;
delete hctlphi; delete hctlphi_moments;
delete hq2ctl; delete hq2ctl_moments;
delete hq2ctk; delete hq2ctk_moments;
delete hq2phi; delete hq2phi_moments;
hxi_moments[0]->SetMinimum(1.0); hxi_moments[0]->SetMaximum(+2.0);
hxi_moments[1]->SetMinimum(0.25); hxi_moments[1]->SetMaximum(+1.25);
hxi_moments[2]->SetMinimum(-1.0); hxi_moments[2]->SetMaximum(+0.0);
hxi_moments[3]->SetMinimum(-0.75); hxi_moments[3]->SetMaximum(+0.25);
for (int i = 4; i<12; i++){
hxi_moments[i]->SetMinimum(-0.5);
hxi_moments[i]->SetMaximum(+0.5);
}
for (auto hist: hxi_moments) hist->GetYaxis()->SetTitleOffset(1.4);
TCanvas* c10 = new TCanvas("c10", "c10", 1600, 1200);
c10->Divide(4,3);
plotAngularInOne(c10,1,hxi_moments[0],hxi_ub[0],hxi_hist[0],hxi_ml_hist[0],runMinuit,edges.q2min,edges.q2max,1.5,false);
plotAngularInOne(c10,2,hxi_moments[1],hxi_ub[1],hxi_hist[1],hxi_ml_hist[1],runMinuit,edges.q2min,edges.q2max,0.75,false);
plotAngularInOne(c10,3,hxi_moments[2],hxi_ub[2],hxi_hist[2],hxi_ml_hist[2],runMinuit,edges.q2min,edges.q2max,-0.5,false);
plotAngularInOne(c10,4,hxi_moments[3],hxi_ub[3],hxi_hist[3],hxi_ml_hist[3],runMinuit,edges.q2min,edges.q2max,-0.25,false);
for (int i = 4; i<12; i++){
plotAngularInOne(c10,i+1,hxi_moments[i],hxi_ub[i],hxi_hist[i],hxi_ml_hist[i],runMinuit,
edges.q2min,edges.q2max,0.0,is_param_folded(i,opts));
}
if(opts->write_eps)c10->Print((subfolder + "/xis"+appendix+".eps").c_str(),"eps");
delete c10;
spdlog::info( "end PLOTS" );
}//end make plots
//delete histograms
for (auto hist: hxi_moments) delete hist;
for (auto hist: hxi_ml_hist) delete hist;
for (auto hist: hxi_hist) delete hist;
for (auto hist: hxi_ub) delete hist;
}; //end parametrize_eff_phsp_4d()
double bu2kstarmumu_pdf::get_ang_eff(fcnc::options *opts,
double q2, double ctl, double ctk, double phi) const{
npolynom npol(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
//Posibly can be faster, but with more memory consumption and unless I calculate this for like a million events it doesn't make a difference. This way the code is readable and usable at many places
double eff = 0.0;
for(unsigned int h = 0; h<npol.q2; h++){
for(unsigned int i = 0; i<npol.ctl; i++){
for(unsigned int j = 0; j<npol.ctk; j++){
for(unsigned int k = 0; k<npol.phi; k++){
eff += coeffs_eff_4d[npol.get_bin_in4D(h,i,j,k)]
*pow(q2,h)*pow(ctl,i)*pow(ctk,j)*pow(phi,k);
}
}
}
}
//spdlog::trace("ctl={0:f}\t ctk={1:f}\t phi={2:f}\t eff={3:f}", ctl, ctk, phi, eff);
return eff;
}
double bu2kstarmumu_pdf::get_ang_eff(fcnc::options *opts, fcnc::event evt, bool fa) const{
//fa=true meansusing the evt. _fa values
if (fa) return get_ang_eff(opts, evt.q2, evt.costhetal_fa, evt.costhetak_fa, evt.phi_fa);
else return get_ang_eff(opts, evt.q2, evt.costhetal, evt.costhetak, evt.phi);
}
void bu2kstarmumu_pdf::eff_fcn(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag)
{
//TODO full angular fit
npolynom npol(current_pdf->opts);
std::vector<double> chebyshev_costhetal(npol.ctl, 0.0);
std::vector<double> chebyshev_costhetak(npol.ctk, 0.0);
std::vector<double> chebyshev_phi(npol.phi, 0.0);
for (unsigned int j = 0; j<npol.ctl; j++) chebyshev_costhetal.at(j) = params[j];
for (unsigned int j = 0; j<npol.ctk; j++) chebyshev_costhetak.at(j) = params[npol.ctl+j];
for (unsigned int j = 0; j<npol.phi; j++) chebyshev_phi.at(j) = params[npol.ctl+npol.ctk+j];
std::vector<double> coeffs_costhetal(npol.ctl, 0.0);
std::vector<double> coeffs_costhetak(npol.ctk, 0.0);
std::vector<double> coeffs_phi(npol.phi, 0.0);
std::vector<double> temp_ctl(npol.ctl, 0.0);
std::vector<double> temp_ctk(npol.ctk, 0.0);
std::vector<double> temp_phi(npol.phi, 0.0);
chebyshev_to_poly(chebyshev_costhetal, temp_ctl);
chebyshev_to_poly(chebyshev_costhetak, temp_ctk);
chebyshev_to_poly(chebyshev_phi, temp_phi);
for (unsigned int i=0; i<coeffs_costhetak.size(); i++){
coeffs_costhetak.at(i) = temp_ctk.at(i);
}
for (unsigned int i=0; i<coeffs_costhetal.size(); i++){
coeffs_costhetal.at(i) = temp_ctl.at(i);
}
correct_poly(temp_phi, coeffs_phi, -TMath::Pi(), +TMath::Pi());
double cijk = 0.0;
double f1 = 0.0, f2 = 0.0, f3 = 0.0, f4 = 0.0, f5 = 0.0,
f6 = 0.0, f7 = 0.0, f8 = 0.0, f9 = 0.0, f10 = 0.0, f11 = 0.0, f12 = 0.0;
const double pi = MY_PI;
double c = C_PWAVE;
if(current_pdf->opts->full_angular)c *= 2.0; //folding in phi
if (current_pdf->opts->full_angular)
{
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++)
{
cijk = c*coeffs_costhetal.at(i)*coeffs_costhetak.at(j)*coeffs_phi.at(k);
f1 += cijk * (pow(-1,i)/(i+1.0)+1.0/(i+1.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))*(pow(-1,k)*pow(pi,k+1)/(k+1.0)+pow(pi,k+1)/(k+1.0));
f2 += cijk * (pow(-1,i)/(i+1.0)+1.0/(i+1.0))*(pow(-1,j)/(j+3.0)+1.0/(j+3.0))*(pow(-1,k)*pow(pi,k+1)/(k+1.0)+pow(pi,k+1)/(k+1.0));
f3 += cijk * (((i-1.0)*pow(-1,i)/(i*i+4.0*i+3.0)+(i-1.0)/(i*i+4.0*i+3.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))*(pow(-1,k)*pow(pi,k+1)/(k+1.0)+pow(pi,k+1)/(k+1.0)));
f4 += cijk * (((i-1.0)*pow(-1,i)/(i*i+4.0*i+3.0)+(i-1.0)/(i*i+4.0*i+3.0))*(pow(-1,j)/(j+3.0)+1/(j+3.0))*(pow(-1,k)*pow(pi,k+1.0)/(k+1.0)+pow(pi,k+1)/(k+1.0)));
f5 += cijk * (2.0*pow(-1,i)/(i*i+4.0*i+3.0)+2.0/(i*i+4.0*i+3.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))
*(integral_x_to_n_times_cos_2x(pi,k)-integral_x_to_n_times_cos_2x(-pi,k));
f6 += cijk * 4.0
*(integral_x_to_n_times_sqrt_1_minus_x2(1,j+1)-integral_x_to_n_times_sqrt_1_minus_x2(-1,j+1))
*(integral_x_to_n_times_sqrt_1_minus_x2(1,i+1)-integral_x_to_n_times_sqrt_1_minus_x2(-1,i+1))
*(integral_x_to_n_times_cos_x(pi,k)-integral_x_to_n_times_cos_x(-pi,k));
f7 += cijk * 2.0
*(integral_x_to_n_times_sqrt_1_minus_x2(1,j+1)-integral_x_to_n_times_sqrt_1_minus_x2(-1,j+1))
*(integral_x_to_n_times_sqrt_1_minus_x2(1,i)-integral_x_to_n_times_sqrt_1_minus_x2(-1,i))
*(integral_x_to_n_times_cos_x(pi,k)-integral_x_to_n_times_cos_x(-pi,k));
f8 += cijk * (1.0/(i+2.0)-pow(-1,i)/(i+2.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))*(pow(-1,k)*pow(pi,k+1)/(k+1.0)+pow(pi,k+1)/(k+1.0));
f9 += cijk * (1.0/(i+2.0)-pow(-1,i)/(i+2.0))*(pow(-1,j)/(j+3.0)+1.0/(j+3.0))*(pow(-1,k)*pow(pi,k+1)/(k+1.0)+pow(pi,k+1)/(k+1.0));
f10 += cijk * 2.0
*(integral_x_to_n_times_sqrt_1_minus_x2(1,j+1)-integral_x_to_n_times_sqrt_1_minus_x2(-1,j+1))
*(integral_x_to_n_times_sqrt_1_minus_x2(1,i)-integral_x_to_n_times_sqrt_1_minus_x2(-1,i))
*(integral_x_to_n_times_sin_x(pi,k)-integral_x_to_n_times_sin_x(-pi,k));
f11 += cijk * 4.0
*(integral_x_to_n_times_sqrt_1_minus_x2(1,j+1)-integral_x_to_n_times_sqrt_1_minus_x2(-1,j+1))
*(integral_x_to_n_times_sqrt_1_minus_x2(1,i+1)-integral_x_to_n_times_sqrt_1_minus_x2(-1,i+1))
*(integral_x_to_n_times_sin_x(pi,k)-integral_x_to_n_times_sin_x(-pi,k));
f12 += cijk * (2.0*pow(-1,i)/(i*i+4.0*i+3.0)+2.0/(i*i+4.0*i+3.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))
*(integral_x_to_n_times_sin_2x(pi,k)-integral_x_to_n_times_sin_2x(-pi,k));
}
}
else
{
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++)
{
cijk = c * coeffs_costhetal.at(i)*coeffs_costhetak.at(j)*coeffs_phi.at(k);
/*
f1 += cijk * 2.0*(2.0*pow(-1,i)/(i*i+4.0*i+3.0)+2.0/(i*i+4.0*i+3.0))*(pow(-1,j)/(j+3.0)+1.0/(j+3.0))*pow(MY_PI,k+1)/(k+1.0);
f2 += cijk * ((i+2.0)*pow(-1,i)/(i*i+4.0*i+3.0)+(i+2.0)/(i*i+4.0*i+3.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))*pow(MY_PI,k+1)/(k+1.0);
f3 += cijk * (2.0*pow(-1,i)/(i*i+4.0*i+3.0)+2.0/(i*i+4.0*i+3.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))
*(integral_x_to_n_times_cos_2x(MY_PI, k) - integral_x_to_n_times_cos_2x(0.0, k));
f4 += cijk * 4.0*(1.0/(i+2.0)-pow(-1,i)/(i+2.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))*pow(MY_PI,k+1)/(3.0*(k+1.0));
f5 += cijk * (2.0*pow(-1,i)/(i*i+4.0*i+3.0)+2.0/(i*i+4.0*i+3.0))*(2.0*pow(-1,j)/(j*j+4.0*j+3.0)+2.0/(j*j+4.0*j+3.0))
*(integral_x_to_n_times_sin_2x(MY_PI, k) - integral_x_to_n_times_sin_2x(0.0, k));
*/
//TODO add f1 to f12 for new folding!
spdlog::error("Function 'eff_fcn' not yet implemented for folding");
assert(0);
}
}
//calculate eps()
std::vector<double> ch_costhetal;
std::vector<double> ch_costhetak;
std::vector<double> ch_phi;
double lh_costhetal = 0.0;
double lh_costhetak = 0.0;
double lh_phi = 0.0;
unsigned int nsize = current_pdf->eff_events.size();
for (unsigned int i=0; i<nsize; i++)
{
event meas = current_pdf->eff_events.at(i);
chebyshev(meas.costhetal, npol.ctl, ch_costhetal);
chebyshev(meas.costhetak, npol.ctk, ch_costhetak);
if (current_pdf->opts->full_angular)
chebyshev(meas.phi/MY_PI, npol.phi, ch_phi);
else
chebyshev(2.0*meas.phi/MY_PI-1.0, npol.phi, ch_phi);
//double result = 0.0;
double result_costhetal = 0.0;
for (unsigned int j = 0; j<npol.ctl; j++){
result_costhetal += params[j]*ch_costhetal.at(j);// /norm_costhetal;
}
lh_costhetal += -2.0*TMath::Log(result_costhetal);
double result_costhetak = 0.0;
for (unsigned int j = 0; j<npol.ctk; j++){
result_costhetak += params[npol.ctl+j]*ch_costhetak.at(j);// /norm_costhetak;
}
lh_costhetak += -2.0*TMath::Log(result_costhetak);
double result_phi = 0.0;
for (unsigned int j = 0; j<npol.phi; j++){
result_phi += params[npol.ctl+npol.ctk+j]*ch_phi.at(j);// /norm_phi;
}
lh_phi += -2.0*TMath::Log(result_phi);
}
//multiply with the pdf
const double j1s = current_pdf->eff_params->J1s();
const double j1c = current_pdf->eff_params->J1c();
const double j2s = current_pdf->eff_params->J2s();
const double j2c = current_pdf->eff_params->J2c();
const double j3 = current_pdf->eff_params->J3();
const double j4 = current_pdf->eff_params->J4();
const double j5 = current_pdf->eff_params->J5();
const double j6s = current_pdf->eff_params->J6s();
const double j6c = current_pdf->eff_params->J6c();
const double j7 = current_pdf->eff_params->J7();
const double j8 = current_pdf->eff_params->J8();
const double j9 = current_pdf->eff_params->J9();
double sig_norm = 0.0;
sig_norm = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9);
double lh_sig = 0.0, prob = 0.0;
for (unsigned int i=0; i<nsize; i++)
{
event meas = current_pdf->eff_events.at(i);
if (current_pdf->opts->full_angular)
{
current_pdf->fj(meas.costhetal, meas.costhetak, meas.phi,
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
prob = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9)/sig_norm;
}
else
{
current_pdf->folded_fj(meas.costhetal, meas.costhetak, meas.phi, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
/*
prob = 9.0/16.0/MY_PI * (fl * f1
+ft * f2
+s3 * f3
+afb * f4
+aim * f5
)/sig_norm;
*/
prob = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9)/sig_norm;
}
lh_sig += -2.0*TMath::Log(prob);
}
//is the integral correct here?
//need to do a nice integration over pdf_sig(costhetal,costhetak,phi)*eps(costhetal,costhetak,phi)
lh = lh_sig + lh_costhetal + lh_costhetak + lh_phi;
};
//this does not assume factorisation of the three dependencies, but assumes phsp MC in 3D
void bu2kstarmumu_pdf::eff_fcn_phsp(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag)
{
//TODO
unsigned int ncosthetal = 5;
unsigned int ncosthetak = 5;
unsigned int nphi = 5;
//make sure we are doing full angular acceptance, no phi folding
assert(current_pdf->opts->full_angular);
//need to account for the norm given the current parameters
double norm = 0.0;
for (unsigned int i = 0; i<ncosthetal; i++)
for (unsigned int j = 0; j<ncosthetak; j++)
for (unsigned int k = 0; k<nphi; k++)
{
double cijk = params[ncosthetak*nphi*i+nphi*j+k];
norm += cijk
*(integral_chebyshev(+1.0, i)-integral_chebyshev(-1.0, i)) //TODO: check the -1,+1
*(integral_chebyshev(+1.0, j)-integral_chebyshev(-1.0, j))
*(integral_chebyshev(+1.0, k)-integral_chebyshev(-1.0, k));//there surely is a more efficient way
}
//this gives Sum_events Log(eps(event))
lh = 0.0;
unsigned int nsize = current_pdf->eff_events.size();
//loop over all events
for (unsigned int l=0; l<nsize; l++) {
const event& meas = current_pdf->eff_events.at(l);
std::vector<double> ch_costhetal;
std::vector<double> ch_costhetak;
std::vector<double> ch_phi;
chebyshev(meas.costhetal, ncosthetal, ch_costhetal);
chebyshev(meas.costhetak, ncosthetak, ch_costhetak);
chebyshev(meas.phi/MY_PI, nphi, ch_phi);
double result = 0.0;
for (unsigned int i = 0; i<ncosthetal; i++)
for (unsigned int j = 0; j<ncosthetak; j++)
for (unsigned int k = 0; k<nphi; k++)
{
double cijk = params[ncosthetak*nphi*i+nphi*j+k];
result += cijk*ch_costhetal.at(i)*ch_costhetak.at(j)*ch_phi.at(k);
}
if (result > 0.0 && !std::isnan(result) && !std::isinf(result))
lh += -2.0*TMath::Log(result/norm) - (-2.0*TMath::Log(1.0/8.0));//can subtract
else
lh += (-100.0*2.0*TMath::Log(1.0/8.0));
}
spdlog::debug( "eff_fcn_phsp lh {0:f}\t npar {1:d}\tiflag {2:d}\tnorm {3:f}", lh, npar, iflag, norm);
};
void bu2kstarmumu_pdf::eff_fcn_phsp_4d(Int_t &npar, Double_t *grad, Double_t &lh, Double_t *params, Int_t iflag)
{
npolynom npol(current_pdf->opts);
double q2min = current_pdf->opts->TheQ2binsmin.front(), q2max = current_pdf->opts->TheQ2binsmax.back();
double ctkmin = current_pdf->opts->ctk_min, ctkmax = current_pdf->opts->ctk_max;
double ctlmin = current_pdf->opts->ctl_min, ctlmax = current_pdf->opts->ctl_max;
double phimin = current_pdf->opts->phi_min, phimax = current_pdf->opts->phi_max;
spdlog::warn("The TMinuit version was not yet tested for folded analysis!");
assert(current_pdf->opts->full_angular);
//need to account for the norm given the current parameters
double norm = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
double chijk = params[npol.get_bin_in4D(h,i,j,k)];
if (chijk == 0.0) continue;
norm += chijk
*(integral_chebyshev(+1.0, h)-integral_chebyshev(-1.0, h))
*(integral_chebyshev(+1.0, i)-integral_chebyshev(-1.0, i))
*(integral_chebyshev(+1.0, j)-integral_chebyshev(-1.0, j))
*(integral_chebyshev(+1.0, k)-integral_chebyshev(-1.0, k));//there surely is a more efficient way
}
}
}
}
//this gives Sum_events Log(eps(event))
lh = 0.0;
unsigned int nsize = current_pdf->eff_events.size();
//loop over all events
for (unsigned int l=0; l<nsize; l++){
const event& meas = current_pdf->eff_events.at(l);
if (meas.q2 < q2min || meas.q2 > q2max) continue;
std::vector<double> ch_q2;
std::vector<double> ch_costhetal;
std::vector<double> ch_costhetak;
std::vector<double> ch_phi;
chebyshev(2.0*(meas.q2 - q2min ) / (q2max - q2min ) - 1.0, npol.q2, ch_q2);
chebyshev(2.0*(meas.costhetal - ctlmin) / (ctlmax - ctlmin) - 1.0, npol.ctl, ch_costhetal);
chebyshev(2.0*(meas.costhetak - ctkmin) / (ctkmax - ctkmin) - 1.0, npol.ctk, ch_costhetak);
chebyshev(2.0*(meas.phi - phimin) / (phimax - phimin) - 1.0, npol.phi, ch_phi);
/*
chebyshev(meas.costhetal, ncosthetal, ch_costhetal);
chebyshev(meas.costhetak, ncosthetak, ch_costhetak);
chebyshev(meas.phi/pi, nphi, ch_phi);
*/
double result = 0.0;
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.q2; k++)
{
double chijk = params[npol.get_bin_in4D(h,i,j,k)];
result += chijk*ch_q2.at(h)*ch_costhetal.at(i)*ch_costhetak.at(j)*ch_phi.at(k);
}
if (result > 0.0 && !std::isnan(result) && !std::isinf(result)){
lh += -2.0*meas.weight*TMath::Log(result/norm) - (-2.0*TMath::Log(1.0/16.0));
}
else{
lh += (-100.0*2.0*TMath::Log(1.0/8.0));
}
}
spdlog::debug( "eff_fcn_phsp_4d lh {0:f}\t npar {1:d}\tiflag {2:d}\tnorm {3:f}", lh, npar, iflag, norm);
};
void bu2kstarmumu_pdf::update_cached_normalization(const bu2kstarmumu_parameters* params)
{
spdlog::trace("updating cached normalisation");
//this caches the overall normalisation for both signal and background
double q2 = params->eff_q2();
//swave norm
if (opts->use_mkpi || opts->fit_mkpi){
double mkpi_min = opts->mkpi_min/1.0e+3;
double mkpi_max = opts->mkpi_max/1.0e+3;
//double q2 = params->eff_q2();
double asphase = params->asphase();
double a = params->a();
double r = params->r();
double gammakstar = params->gammakstar();
double mkstar = params->mkstar();
double gammakstarplus = params->gammakstarplus();
double mkstarplus = params->mkstarplus();
double R = params->R();
double gammaf800 = params->gammaf800();
double mf800 = params->mf800();
double f800mag = params->f800mag();
double f800phase = params->f800phase();
//integrate over mkpi2
const double mk = PDGMASS_KST_KAON/1.0e+3;
const double mpi = PDGMASS_KST_PION/1.0e+3;
const double mb = PDGMASS_B/1.0e+3;
const double from = opts->mkpi_full_range_norm ? (mk+mpi) : mkpi_min;
const double to = opts->mkpi_full_range_norm ? mb : mkpi_max;
std::complex<double> mkpi_swavepwave_norm;
if (opts->simple_mkpi){
mkpi_swave_2_norm = int_mkpi_simple_kstarzero_amp_squared(from, to, q2,
f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus);
mkpi_pwave_2_norm = int_mkpi_simple_bw_kstar_amp_squared(from, to, q2, gammakstar, mkstar);
mkpi_swavepwave_norm = int_mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(from, to, q2,
gammakstar, mkstar, asphase,
f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus);
mkpi_re_swavepwave_norm = mkpi_swavepwave_norm.real();
mkpi_im_swavepwave_norm = mkpi_swavepwave_norm.imag();
}
else if (opts->isobar){
mkpi_swave_2_norm = int_mkpi_kstarzero_isobar_amp_squared(from, to, q2,
f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R);
mkpi_pwave_2_norm = int_mkpi_bw_kstar_amp_squared(from, to, q2, gammakstar, mkstar, R);
mkpi_swavepwave_norm = int_mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(from, to, q2,
gammakstar, mkstar, asphase,
f800mag, f800phase, gammaf800, mf800, gammakstarplus, mkstarplus, R);
mkpi_re_swavepwave_norm = mkpi_swavepwave_norm.real();
mkpi_im_swavepwave_norm = mkpi_swavepwave_norm.imag();
}
else{//nominal lass
mkpi_swave_2_norm = int_mkpi_kstarzero_lass_amp_squared(from, to, q2, asphase, a, r, gammakstarplus, mkstarplus, R);
mkpi_pwave_2_norm = int_mkpi_bw_kstar_amp_squared(from, to, q2, gammakstar, mkstar, R);
mkpi_swavepwave_norm = int_mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(from, to, q2,
gammakstar, mkstar,
asphase, a, r, gammakstarplus, mkstarplus, R);
mkpi_re_swavepwave_norm = mkpi_swavepwave_norm.real();
mkpi_im_swavepwave_norm = mkpi_swavepwave_norm.imag();
}
if (opts->use_p2){
//mkpi_bkg(
std::vector<double> ch_p2(5, 0.0);
ch_p2.at(0) = params->cbkgp20();
ch_p2.at(1) = params->cbkgp21();
ch_p2.at(2) = params->cbkgp22();
ch_p2.at(3) = params->cbkgp23();
ch_p2.at(4) = params->cbkgp24();
std::vector<double> poly_p2(ch_p2.size(), 0.0);
chebyshev_to_poly(ch_p2, poly_p2);
std::vector<double> poly_corr_p2(ch_p2.size(), 0.0);
correct_poly(poly_p2, poly_corr_p2, mkpi_min*mkpi_min, mkpi_max*mkpi_max);
mkpi_bkg_norm = 0.0;
for (unsigned int i=0; i<ch_p2.size(); i++)
mkpi_bkg_norm += poly_corr_p2.at(i)*(pow(mkpi_max*mkpi_max,i+1)-pow(mkpi_min*mkpi_min,i+1))/(i+1.0);
}
else{
std::vector<double> ch_mkpi(5, 0.0);
ch_mkpi.at(0) = params->cbkgmkpi0();
ch_mkpi.at(1) = params->cbkgmkpi1();
ch_mkpi.at(2) = params->cbkgmkpi2();
ch_mkpi.at(3) = params->cbkgmkpi3();
ch_mkpi.at(4) = params->cbkgmkpi4();
std::vector<double> poly_mkpi(ch_mkpi.size(), 0.0);
chebyshev_to_poly(ch_mkpi, poly_mkpi);
std::vector<double> poly_corr_mkpi(ch_mkpi.size(), 0.0);
correct_poly(poly_mkpi, poly_corr_mkpi, mkpi_min, mkpi_max);
mkpi_bkg_norm = 0.0;
if (opts->mkpi_threshold){
mkpi_bkg_norm = fcnc::int_threshold(mkpi_min, mkpi_max, 0.633, params->nthreshold());
}
else{
for(unsigned int i=0; i<poly_mkpi.size(); i++)
mkpi_bkg_norm += poly_corr_mkpi.at(i)*(pow(mkpi_max,i+1)-pow(mkpi_min,i+1))/(i+1.0);
}
}
}
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
double fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0, fs6=0.0;
const double j1s = params->J1s();
const double j1c = params->J1c();
const double j2s = params->J2s();
const double j2c = params->J2c();
const double j3 = params->J3();
const double j4 = params->J4();
const double j5 = params->J5();
const double j6s = params->J6s();
const double j6c = params->J6c();
const double j7 = params->J7();
const double j8 = params->J8();
const double j9 = params->J9();
const double fs = params->FS();
const double js1 = params->SS1();
const double js2 = params->SS2();
const double js3 = params->SS3();
const double js4 = params->SS4();
const double js5 = params->SS5();
if(spdlog_trace()){ //print out the integrated angular moments. with and without weights, as well as their ratios
load_coeffs_eff_phsp_4d();
integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12, q2);
double pwave_weighted = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9);
integrated_fj_noacc(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
double pwave_noweights = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9);
swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
spdlog::trace("Integrated P-wave:");
spdlog::trace("P-wave no weights: {0:0.5f}",pwave_noweights);
spdlog::trace("P-wave weighted: {0:0.5f}",pwave_weighted);
spdlog::trace("P-wave ratio: {0:0.5f}\n",pwave_weighted / pwave_noweights);
if (opts->swave){
double swave_noweights = (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
double swave_weighted = (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
spdlog::trace("Integrated S-wave:");
spdlog::trace("S-wave weighted: {0:0.5f}{0:0.5f}",swave_weighted);
spdlog::trace("S-wave no weights: {0:0.5f}",swave_noweights);
spdlog::trace("S-wave ratio: {0:0.5f}\n",swave_weighted / swave_noweights);
spdlog::trace("Double ratio (s/p):{0:0.5f}",(swave_weighted / swave_noweights) / (pwave_weighted / pwave_noweights));
}
//assert(0);
}
if(opts->full_angular){
if(!opts->cache_fis){
if(opts->use_angular_acc){
integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12, q2);
}
else{//this is also what the weighted fit uses!
integrated_fj_noacc(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
}
fbuffer_f1 = f1;
fbuffer_f2 = f2;
fbuffer_f3 = f3;
fbuffer_f4 = f4;
fbuffer_f5 = f5;
fbuffer_f6 = f6;
fbuffer_f7 = f7;
fbuffer_f8 = f8;
fbuffer_f9 = f9;
fbuffer_f10 = f10;
fbuffer_f11 = f11;
fbuffer_f12 = f12;
}
else{
f1 = fbuffer_f1;
f2 = fbuffer_f2;
f3 = fbuffer_f3;
f4 = fbuffer_f4;
f5 = fbuffer_f5;
f6 = fbuffer_f6;
f7 = fbuffer_f7;
f8 = fbuffer_f8;
f9 = fbuffer_f9;
f10 = fbuffer_f10;
f11 = fbuffer_f11;
f12 = fbuffer_f12;
}
fnorm_sig = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9);
if (opts->fit_asymmetries){
fnorm_sig = (f1*j1s + f2*j1c + f3*j2s + f4*j2c);
}
bool swave = opts->swave;
if (swave){
if (!opts->cache_fis){
if (opts->use_angular_acc){
swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
}
else{//weighted fit and fit without acceptance
swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
fbuffer_fs1 = fs1;
fbuffer_fs2 = fs2;
fbuffer_fs3 = fs3;
fbuffer_fs4 = fs4;
fbuffer_fs5 = fs5;
fbuffer_fs6 = fs6;
}
else{
fs1 = fbuffer_fs1;
fs2 = fbuffer_fs2;
fs3 = fbuffer_fs3;
fs4 = fbuffer_fs4;
fs5 = fbuffer_fs5;
fs6 = fbuffer_fs6;
}
fnorm_sig = (1.0-fs) * (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9)
+ (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
if (opts->fit_asymmetries){
fnorm_sig = (1.0-fs)*(f1*j1s + f2*j1c + f3*j2s + f4*j2c)
+ (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
}
}
spdlog::trace("Fnormsig (full_angular) {0:f}",fnorm_sig );
spdlog::trace("Normalization integrals P-wave");
spdlog::trace("\tj1s {0:0.4f}",f1/fnorm_sig);
spdlog::trace("\tj1c {0:0.4f}",f2/fnorm_sig);
spdlog::trace("\tj2s {0:0.4f}",f3/fnorm_sig);
spdlog::trace("\tj2c {0:0.4f}",f4/fnorm_sig);
spdlog::trace("\tj3 {0:0.4f}",f5/fnorm_sig);
spdlog::trace("\tj4 {0:0.4f}",f6/fnorm_sig);
spdlog::trace("\tj5 {0:0.4f}",f7/fnorm_sig);
spdlog::trace("\tj6s {0:0.4f}",f8/fnorm_sig);
spdlog::trace("\tj6c {0:0.4f}",f9/fnorm_sig);
spdlog::trace("\tj7 {0:0.4f}",f10/fnorm_sig);
spdlog::trace("\tj8 {0:0.4f}",f11/fnorm_sig);
spdlog::trace("\tj9 {0:0.4f}",f12/fnorm_sig);
if(opts->swave){
spdlog::trace("Normalization integrals S-wave");
spdlog::trace("fs {0:0.4f}",fs1/fnorm_sig);
spdlog::trace("\tas1 {0:0.4f}",fs2/fnorm_sig);
spdlog::trace("\tas2 {0:0.4f}",fs3/fnorm_sig);
spdlog::trace("\tas3 {0:0.4f}",fs4/fnorm_sig);
spdlog::trace("\tas4 {0:0.4f}",fs5/fnorm_sig);
spdlog::trace("\tas5 {0:0.4f}",fs6/fnorm_sig);
}
spdlog::trace("Angular amplitudes: ");
spdlog::trace("\tj1s {0:0.4f}",j1s);
spdlog::trace("\tj1c {0:0.4f}",j1c);
spdlog::trace("\tj2s {0:0.4f}",j2s);
spdlog::trace("\tj2c {0:0.4f}",j2c);
spdlog::trace("\tj3 {0:0.4f}",j3);
spdlog::trace("\tj4 {0:0.4f}",j4);
spdlog::trace("\tj5 {0:0.4f}",j5);
spdlog::trace("\tj6s {0:0.4f}",j6s);
spdlog::trace("\tj6c {0:0.4f}",j6c);
spdlog::trace("\tj7 {0:0.4f}",j7);
spdlog::trace("\tj8 {0:0.4f}",j8);
spdlog::trace("\tj9 {0:0.4f}",j9);
}
else //folded normalization!
{
if (!opts->cache_fis) {
if (opts->use_angular_acc){
folded_integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12, q2);
}
else{//this is also what the weighted fit uses!
folded_integrated_fj_noacc(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12);
}
fbuffer_f1 = f1;
fbuffer_f2 = f2;
fbuffer_f3 = f3;
fbuffer_f4 = f4;
fbuffer_f5 = f5;
fbuffer_f6 = f6;
fbuffer_f7 = f7;
fbuffer_f8 = f8;
fbuffer_f9 = f9;
fbuffer_f10 = f10;
fbuffer_f11 = f11;
fbuffer_f12 = f12;
}
else{
f1 = fbuffer_f1;
f2 = fbuffer_f2;
f3 = fbuffer_f3;
f4 = fbuffer_f4;
f5 = fbuffer_f5;
f6 = fbuffer_f6;
f7 = fbuffer_f7;
f8 = fbuffer_f8;
f9 = fbuffer_f9;
f10 = fbuffer_f10;
f11 = fbuffer_f11;
f12 = fbuffer_f12;
}
fnorm_sig = (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9);
if (opts->fit_asymmetries){
fnorm_sig = (f1*j1s + f2*j1c + f3*j2s + f4*j2c);
}
bool swave = opts->swave;
if (swave){
if (!opts->cache_fis) {
if (opts->use_angular_acc ){
folded_swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
}
else{//weighted fit and fit without acceptance
folded_swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
fbuffer_fs1 = fs1;
fbuffer_fs2 = fs2;
fbuffer_fs3 = fs3;
fbuffer_fs4 = fs4;
fbuffer_fs5 = fs5;
fbuffer_fs6 = fs6;
}
else {
fs1 = fbuffer_fs1;
fs2 = fbuffer_fs2;
fs3 = fbuffer_fs3;
fs4 = fbuffer_fs4;
fs5 = fbuffer_fs5;
fs6 = fbuffer_fs6;
}
fnorm_sig = (1-fs) * (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9)
+ (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
if (opts->fit_asymmetries)
fnorm_sig = (1-fs)*(f1*j1s + f2*j1c + f3*j2s + f4*j2c)
+(fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
}
spdlog::trace("Fnormsig (folding {0:d}) {1:f}", opts->folding, fnorm_sig );
spdlog::trace("Normalization integrals P-wave");
spdlog::trace("\tj1s {0:0.4f}",f1/fnorm_sig);
spdlog::trace("\tj1c {0:0.4f}",f2/fnorm_sig);
spdlog::trace("\tj2s {0:0.4f}",f3/fnorm_sig);
spdlog::trace("\tj2c {0:0.4f}",f4/fnorm_sig);
spdlog::trace("\tj3 {0:0.4f}",f5/fnorm_sig);
spdlog::trace("\tj4 {0:0.4f}",f6/fnorm_sig);
spdlog::trace("\tj5 {0:0.4f}",f7/fnorm_sig);
spdlog::trace("\tj6s {0:0.4f}",f8/fnorm_sig);
spdlog::trace("\tj6c {0:0.4f}",f9/fnorm_sig);
spdlog::trace("\tj7 {0:0.4f}",f10/fnorm_sig);
spdlog::trace("\tj8 {0:0.4f}",f11/fnorm_sig);
spdlog::trace("\tj9 {0:0.4f}",f12/fnorm_sig);
if(opts->swave){
spdlog::trace("Normalization integrals S-wave");
spdlog::trace("fs {0:0.4f}",fs1/fnorm_sig);
spdlog::trace("\tas1 {0:0.4f}",fs2/fnorm_sig);
spdlog::trace("\tas2 {0:0.4f}",fs3/fnorm_sig);
spdlog::trace("\tas3 {0:0.4f}",fs4/fnorm_sig);
spdlog::trace("\tas4 {0:0.4f}",fs5/fnorm_sig);
spdlog::trace("\tas5 {0:0.4f}",fs6/fnorm_sig);
}
spdlog::trace("Angular amplitudes: ");
spdlog::trace("\tj1s {0:0.4f}",j1s);
spdlog::trace("\tj1c {0:0.4f}",j1c);
spdlog::trace("\tj2s {0:0.4f}",j2s);
spdlog::trace("\tj2c {0:0.4f}",j2c);
spdlog::trace("\tj3 {0:0.4f}",j3);
spdlog::trace("\tj4 {0:0.4f}",j4);
spdlog::trace("\tj5 {0:0.4f}",j5);
spdlog::trace("\tj6s {0:0.4f}",j6s);
spdlog::trace("\tj6c {0:0.4f}",j6c);
spdlog::trace("\tj7 {0:0.4f}",j7);
spdlog::trace("\tj8 {0:0.4f}",j8);
spdlog::trace("\tj9 {0:0.4f}",j9);
}
//bkg norm
//this is used to convolute the acceptance into the pdf
if(opts->use_angular_acc && opts->use_weighted_bkg){ //weighted background pdf
if (opts->full_angular){
fnorm_bkg = integral_bkg_chebyshev(params, CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX, q2);
}
else if(opts->fit_full_angular_bkg){
fnorm_bkg = integral_bkg_chebyshev(params, CTL_MIN, CTL_MAX, opts->ctk_min, opts->ctk_max, PHI_MIN, PHI_MAX, q2); //TODO: check here whether this makes any sense (https://gitlab.cern.ch/LHCb-RD/ewp-Bplus2Kstmumu-AngAna/-/merge_requests/21)
}
else{
fnorm_bkg = integral_bkg_chebyshev(params, opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max, q2);
}
}
else{ //non weighted background pdf
if (opts->full_angular){
if (opts->flat_bkg) fnorm_bkg = RANGE_3D();
else fnorm_bkg = integral_bkg(params, opts->ctl_min, opts->ctl_max,
opts->ctk_min, opts->ctk_max,
opts->phi_min, opts->phi_max);
}
else{ //folded angles
if (opts->flat_bkg){
fnorm_bkg = (opts->ctk_max - opts->ctk_min) * (opts->ctl_max - opts->ctl_min) * (opts->phi_max - opts->phi_min);
}
else{
if(opts->fit_full_angular_bkg){
fnorm_bkg = integral_bkg(params, CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX);
}
else{
fnorm_bkg = integral_bkg(params, opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max);
}
}
}
}
}
void bu2kstarmumu_pdf::update_cached_integrated_fis(const bu2kstarmumu_parameters* params)
{
if(opts->cache_fis){
double q2 = params->eff_q2();
if(opts->full_angular){
if(opts->use_angular_acc){
integrated_fj_chebyshev(fbuffer_f1, fbuffer_f2, fbuffer_f3, fbuffer_f4, fbuffer_f5, fbuffer_f6,
fbuffer_f7, fbuffer_f8, fbuffer_f9, fbuffer_f10, fbuffer_f11, fbuffer_f12, q2);
swave_integrated_fj_chebyshev(fbuffer_fs1, fbuffer_fs2, fbuffer_fs3, fbuffer_fs4, fbuffer_fs5, fbuffer_fs6, q2);
}
else{//this is also what the weighted fit uses!
integrated_fj_noacc(fbuffer_f1, fbuffer_f2, fbuffer_f3, fbuffer_f4, fbuffer_f5, fbuffer_f6,
fbuffer_f7, fbuffer_f8, fbuffer_f9, fbuffer_f10, fbuffer_f11, fbuffer_f12);
swave_integrated_fj_noacc(fbuffer_fs1, fbuffer_fs2, fbuffer_fs3, fbuffer_fs4, fbuffer_fs5, fbuffer_fs6);
}
}
else{ //folded
if(opts->use_angular_acc){
folded_integrated_fj_chebyshev(fbuffer_f1, fbuffer_f2, fbuffer_f3, fbuffer_f4, fbuffer_f5, fbuffer_f6,
fbuffer_f7, fbuffer_f8, fbuffer_f9, fbuffer_f10, fbuffer_f11, fbuffer_f12, q2);
folded_swave_integrated_fj_chebyshev(fbuffer_fs1, fbuffer_fs2, fbuffer_fs3, fbuffer_fs4, fbuffer_fs5, fbuffer_fs6, q2);
}
else{//this is also what the weighted fit uses!
folded_integrated_fj_noacc(fbuffer_f1, fbuffer_f2, fbuffer_f3, fbuffer_f4, fbuffer_f5, fbuffer_f6,
fbuffer_f7, fbuffer_f8, fbuffer_f9, fbuffer_f10, fbuffer_f11, fbuffer_f12);
folded_swave_integrated_fj_noacc(fbuffer_fs1, fbuffer_fs2, fbuffer_fs3, fbuffer_fs4, fbuffer_fs5, fbuffer_fs6);
}
}
}
}
double bu2kstarmumu_pdf::mkpi_pwave_2(const bu2kstarmumu_parameters* params, const event& meas) const
{
double mkpi = meas.mkpi/1.0e+3;
double q2 = params->eff_q2();
//if (opts->weighted_fit) TODO
// q2 = meas.q2;
double gammakstar = params->gammakstar();
double mkstar = params->mkstar();
double R = params->R();
if (opts->simple_mkpi) return mkpi_simple_bw_kstar_amp_squared(mkpi, q2, gammakstar, mkstar);
else return mkpi_bw_kstar_amp_squared(mkpi, q2, gammakstar, mkstar, R);
}
double bu2kstarmumu_pdf::mkpi_swave_2(const bu2kstarmumu_parameters* params, const event& meas) const
{
double mkpi = meas.mkpi/1.0e+3;
double q2 = params->eff_q2();
//if (opts->weighted_fit) TODO
// q2 = meas.q2;
double asphase = params->asphase();
double a = params->a();
double r = params->r();
double R = params->R();
double gammakstarplus = params->gammakstarplus();
double mkstarplus = params->mkstarplus();
double gammaf800 = params->gammaf800();
double mf800 = params->mf800();
double f800mag = params->f800mag();
double f800phase = params->f800phase();
if (opts->simple_mkpi){
return mkpi_simple_kstarzero_amp_squared(mkpi, q2,
f800mag, f800phase, gammaf800, mf800,
gammakstarplus, mkstarplus);
}
else if (opts->isobar){
return mkpi_kstarzero_isobar_amp_squared(mkpi, q2,
f800mag, f800phase, gammaf800, mf800,
gammakstarplus, mkstarplus, R);
}
else{//LASS
return mkpi_kstarzero_lass_amp_squared(mkpi, q2, asphase, a, r, gammakstarplus, mkstarplus, R);
}
}
double bu2kstarmumu_pdf::mkpi_re_swavepwave(const bu2kstarmumu_parameters* params, const event& meas) const
{
double mkpi = meas.mkpi/1.0e+3;
double q2 = params->eff_q2();
// if (opts->weighted_fit) //TODO
// q2 = meas.q2;
double gammakstar = params->gammakstar();
double mkstar = params->mkstar();
double asphase = params->asphase();
double a = params->a();
double r = params->r();
double gammakstarplus = params->gammakstarplus();
double mkstarplus = params->mkstarplus();
double R = params->R();
double gammaf800 = params->gammaf800();
double mf800 = params->mf800();
double f800mag = params->f800mag();
double f800phase = params->f800phase();
std::complex<double> swavepwave;
if (opts->simple_mkpi){
swavepwave = mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(mkpi, q2,
gammakstar, mkstar, asphase,
f800mag, f800phase, gammaf800, mf800,
gammakstarplus, mkstarplus);
}
else if (opts->isobar){
swavepwave = mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(mkpi, q2,
gammakstar, mkstar, asphase,
f800mag, f800phase, gammaf800, mf800,
gammakstarplus, mkstarplus, R);
}
else{//LASS
swavepwave = mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(mkpi, q2, gammakstar, mkstar, asphase,
a, r, gammakstarplus, mkstarplus, R);
}
return swavepwave.real();
}
double bu2kstarmumu_pdf::mkpi_im_swavepwave(const bu2kstarmumu_parameters* params, const event& meas) const
{
double mkpi = meas.mkpi/1.0e+3;
double q2 = params->eff_q2();
//if (opts->weighted_fit) //TODO
// q2 = meas.q2;
double gammakstar = params->gammakstar();
double mkstar = params->mkstar();
double asphase = params->asphase();
double a = params->a();
double r = params->r();
double gammakstarplus = params->gammakstarplus();
double mkstarplus = params->mkstarplus();
double R = params->R();
double gammaf800 = params->gammaf800();
double mf800 = params->mf800();
double f800mag = params->f800mag();
double f800phase = params->f800phase();
std::complex<double> swavepwave;
if (opts->simple_mkpi){
swavepwave = mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(mkpi, q2,
gammakstar, mkstar, asphase,
f800mag, f800phase, gammaf800, mf800,
gammakstarplus, mkstarplus);
}
else if (opts->isobar){
swavepwave = mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(mkpi, q2,
gammakstar, mkstar, asphase,
f800mag, f800phase, gammaf800, mf800,
gammakstarplus, mkstarplus, R);
}
else{//LASS
swavepwave = mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(mkpi, q2, gammakstar, mkstar, asphase,
a, r, gammakstarplus, mkstarplus, R);
}
return swavepwave.imag();
}
double bu2kstarmumu_pdf::mkpi_bkg(const bu2kstarmumu_parameters* params, const event& meas) const
{
double mkpi = meas.mkpi/1.0e+3;
double mkpi_min = opts->mkpi_min/1.0e+3;
double mkpi_max = opts->mkpi_max/1.0e+3;
double probbkg = 0.0;
if (opts->use_p2) {
std::vector<double> ch_p2(5, 0.0);
ch_p2.at(0) = params->cbkgp20();
ch_p2.at(1) = params->cbkgp21();
ch_p2.at(2) = params->cbkgp22();
ch_p2.at(3) = params->cbkgp23();
ch_p2.at(4) = params->cbkgp24();
std::vector<double> poly_p2(ch_p2.size(), 0.0);
chebyshev_to_poly(ch_p2, poly_p2);
std::vector<double> poly_corr_p2(ch_p2.size(), 0.0);
correct_poly(poly_p2, poly_corr_p2, mkpi_min*mkpi_min, mkpi_max*mkpi_max);
for (unsigned int i=0; i<ch_p2.size(); i++)
probbkg += poly_corr_p2.at(i)*pow(mkpi*mkpi,i);
}
else
{
std::vector<double> ch_mkpi(5, 0.0);
ch_mkpi.at(0) = params->cbkgmkpi0();
ch_mkpi.at(1) = params->cbkgmkpi1();
ch_mkpi.at(2) = params->cbkgmkpi2();
ch_mkpi.at(3) = params->cbkgmkpi3();
ch_mkpi.at(4) = params->cbkgmkpi4();
std::vector<double> poly_mkpi(ch_mkpi.size(), 0.0);
chebyshev_to_poly(ch_mkpi, poly_mkpi);
std::vector<double> poly_corr_mkpi(ch_mkpi.size(), 0.0);
correct_poly(poly_mkpi, poly_corr_mkpi, mkpi_min, mkpi_max);
if (opts->mkpi_threshold){
//threshold(double x, double thr, double n);
probbkg = fcnc::threshold(mkpi, 0.633, params->nthreshold());
}
else{
for (unsigned int i=0; i<ch_mkpi.size(); i++) probbkg += poly_corr_mkpi.at(i)*pow(mkpi,i);
}
}
return probbkg;
}
void bu2kstarmumu_pdf::update_cached_xis(const bu2kstarmumu_parameters* params, std::vector<event>* events)
{
spdlog::info("Updating cached xis");
//this is used if per_event_norm is on
//it stores the 18 xis for every event
//otherwise using a different norm per event is computationally not possible
//double q2 = params->eff_q2();
assert(opts->use_event_norm);
/*
for (unsigned int e=0; e<events->size(); e++)
{
if ((e*100) % events->size() == 0)
std::cout << (e*100)/events->size() << "% " << std::endl;
event& meas = events->at(e);D
double q2 = meas.q2;
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12, q2);
meas.xis[0] = f1;
meas.xis[1] = f2;
meas.xis[2] = f3;
meas.xis[3] = f4;
meas.xis[4] = f5;
meas.xis[5] = f6;
meas.xis[6] = f7;
meas.xis[7] = f8;
meas.xis[8] = f9;
meas.xis[9] = f10;
meas.xis[10] = f11;
meas.xis[11] = f12;
bool swave = opts->swave;
if (swave)
{
double fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0, fs6=0.0;
swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
meas.xis[12] = fs1;
meas.xis[13] = fs2;
meas.xis[14] = fs3;
meas.xis[15] = fs4;
meas.xis[16] = fs5;
meas.xis[17] = fs6;
}
}
*/
/*
for (unsigned int e=0; e<events->size(); e++)
{
event& meas = events->at(e);
std::cout << meas.xis[0] << " "
<< meas.xis[1] << " "
<< meas.xis[2] << " "
<< meas.xis[3] << " "
<< meas.xis[4] << " "
<< meas.xis[5] << " "
<< meas.xis[6] << " "
<< meas.xis[7] << " "
<< meas.xis[8] << " "
<< meas.xis[9] << " "
<< meas.xis[10] << " "
<< meas.xis[11] << " "
<< meas.xis[12] << " "
<< meas.xis[13] << " "
<< meas.xis[14] << " "
<< meas.xis[15] << " "
<< meas.xis[16] << " "
<< meas.xis[17] << " q2 " << meas.q2 << std::endl;
}
*/
//std::cout << std::endl;
}
void bu2kstarmumu_pdf::update_cached_normalization(const parameters* params){
update_cached_normalization(static_cast<const bu2kstarmumu_parameters*>(params));
}
void bu2kstarmumu_pdf::update_cached_integrated_fis(const parameters* params){
update_cached_integrated_fis(static_cast<const bu2kstarmumu_parameters*>(params));
}
void bu2kstarmumu_pdf::update_cached_efficiencies(const parameters* params, std::vector<event>* events)
{
update_cached_efficiencies(static_cast<const bu2kstarmumu_parameters*>(params), events);
}
void bu2kstarmumu_pdf::update_cached_xis(const parameters* params, std::vector<event>* events)
{
update_cached_xis(static_cast<const bu2kstarmumu_parameters*>(params), events);
}
double bu2kstarmumu_pdf::m_sig_prob(const bu2kstarmumu_parameters* params, const event& meas) const{
if (opts->only_angles || opts->only_mkpi) return 1.0;
double mmin=params->m_b.get_min();
double mmax=params->m_b.get_max();
if(meas.m > mmax || meas.m < mmin){
spdlog::error("m(B) outside of allowed range: mmin={0:f}, mmax={1:f} m={2:f}", mmin, mmax, meas.m);
assert(0);
}
double result = 0.0;
if(opts->crystalball){
result = params->m_res_1()
* crystalball(meas.m, params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->n_1())
/ integrate_crystalball(params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->n_1(), mmin, mmax);
if(params->m_res_1() != 1.0)
result += (1.0-params->m_res_1())
* crystalball(meas.m, params->m_b(), params->m_scale()*params->m_sigma_2(), params->alpha_2(), params->n_2())
/ integrate_crystalball(params->m_b(), params->m_scale()*params->m_sigma_2(), params->alpha_2(), params->n_2(), mmin, mmax);
}
else if(opts->twotailedcrystalball){
result = twotailedcrystalball(meas.m, params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->alpha_2(), params->n_1(), params->n_2())
/ integrate_twotailedcrystalball(params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->alpha_2(), params->n_1(), params->n_2(), mmin, mmax);
}
else{
result = params->m_res_1()
* normed_gauss(meas.m, params->m_scale()*params->m_sigma_1(), params->m_b(), mmin, mmax);
if(params->m_res_1() != 1.0)
result += (1.0-params->m_res_1())
* normed_gauss(meas.m, params->m_scale()*params->m_sigma_2(), params->m_b(), mmin, mmax);
}
return result;
};
double bu2kstarmumu_pdf::mkpi_sig_prob(const bu2kstarmumu_parameters* params, const event& meas) const
{
if(!opts->fit_mkpi || opts->only_angles || opts->only_Bmass) return 1.0;
assert(meas.mkpi <= opts->mkpi_max);
assert(meas.mkpi >= opts->mkpi_min);
double result = 0.0;
double prob_mkpi_pwave = mkpi_pwave_2(params, meas);
double prob_mkpi_swave = mkpi_swave_2(params, meas);
double prob_mkpi_pwaveswave_re = mkpi_re_swavepwave(params, meas);
double prob_mkpi_pwaveswave_im = mkpi_im_swavepwave(params, meas);
double fs = params->FS();
//only m(B) and m(K*+) fitted
if(opts->only_mass2DFit){
if(opts->swave){
result = (1.0-fs)
* prob_mkpi_pwave/mkpi_pwave_2_norm
+ fs
* prob_mkpi_swave/mkpi_swave_2_norm;
}
else{
result = prob_mkpi_pwave/mkpi_pwave_2_norm;
}
}//end only m(B) and m(K*+)
else{//angular dimensions included in the fit
double q2 = meas.q2;
//integrated S_i
double xi1, xi2, xi3, xi4, xi5, xi6, xi7, xi8, xi9, xi10, xi11, xi12;
if(opts->use_angular_acc) integrated_fj_chebyshev(xi1, xi2, xi3, xi4, xi5, xi6, xi7, xi8, xi9, xi10, xi11, xi12, q2);
else integrated_fj_noacc(xi1, xi2, xi3, xi4, xi5, xi6, xi7, xi8, xi9, xi10, xi11, xi12);
//integrated SS_i
double sxi1, sxi2, sxi3, sxi4, sxi5, sxi6;
if(opts->use_angular_acc) swave_integrated_fj_chebyshev(sxi1, sxi2, sxi3, sxi4, sxi5, sxi6, q2);
else swave_integrated_fj_noacc(sxi1, sxi2, sxi3, sxi4, sxi5, sxi6);
//load values for required parameters:
const double s1s = params->S1s();
const double s2s = s1s/3.0;
const double s1c = 1.0 - 4.0/3.0 * s1s;
const double s2c = -s1c;
const double s3 = params->S3();
const double s4 = params->S4();
const double s5 = params->S5();
const double s6s = params->S6s();
const double s6c = 0.;
const double s7 = params->S7();
const double s8 = params->S8();
const double s9 = params->S9();
const double ss1 = params->SS1();
const double ss2 = params->SS2();
const double ss3 = params->SS3();
const double ss4 = params->SS4();
const double ss5 = params->SS5();
double norm = (1 - fs) * (xi1*s1s + xi2*s1c + xi3*s2s + xi4*s2c + xi5*s3 + xi6*s4
+ xi7*s5 + xi8*s6s + xi9*s6c + xi10*s7 + xi11*s8 + xi12*s9)
+ (sxi1*fs)
+ (sxi2*ss1 + sxi3*ss2 + sxi4*ss3)
+ (sxi5*ss4 + sxi6*ss5);
double A = (xi1*s1s + xi2*s1c + xi3*s2s + xi4*s2c + xi5*s3 + xi6*s4
+ xi7*s5 + xi8*s6s + xi9*s6c + xi10*s7 + xi11*s8 + xi12*s9)/norm;
double B = sxi1/norm;
double C = (sxi2*ss1 + sxi3*ss2 + sxi4*ss3)/norm;
double D = (sxi5*ss4 + sxi6*ss5)/norm;
result = (1.0-fs) * A * prob_mkpi_pwave/mkpi_pwave_2_norm
+ fs * B * prob_mkpi_swave/mkpi_swave_2_norm
+ C * prob_mkpi_pwaveswave_re/mkpi_re_swavepwave_norm
+ D * prob_mkpi_pwaveswave_im/mkpi_im_swavepwave_norm;
}//end of prob with angular parts included
return result;
};
double bu2kstarmumu_pdf::integral_m_sig_prob(const bu2kstarmumu_parameters* params, double ma, double mb) const
{
double result = 0.0;
if (opts->crystalball){
double mmin=params->m_b.get_min(), mmax=params->m_b.get_max();
result = params->m_res_1()*integrate_crystalball(params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->n_1(), ma, mb)
/ integrate_crystalball(params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->n_1(), mmin, mmax);
if (params->m_res_1() != 1.0)
result += (1.0-params->m_res_1()) * integrate_crystalball(params->m_b(), params->m_scale()*params->m_sigma_2(), params->alpha_2(), params->n_2(), ma, mb)
/ integrate_crystalball(params->m_b(), params->m_scale()*params->m_sigma_2(), params->alpha_2(), params->n_2(), mmin, mmax);
}
else if(opts->twotailedcrystalball){
double mmin=params->m_b.get_min(), mmax=params->m_b.get_max();
result = integrate_twotailedcrystalball(params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->alpha_2(), params->n_1(), params->n_2(), ma, mb)
/ integrate_twotailedcrystalball(params->m_b(), params->m_scale()*params->m_sigma_1(), params->alpha_1(), params->alpha_2(), params->n_1(), params->n_2(), mmin, mmax);
}
else{
double mean = params->m_b();
double min = params->m_b.get_min();
double max = params->m_b.get_max();
double sigma1 = params->m_scale()*params->m_sigma_1();
double sigma2 = params->m_scale()*params->m_sigma_2();
double fres = params->m_res_1();
double fullrange = (fres)*1.0/2.0*(TMath::Erf((max-mean)/(sqrt(2.0)*sigma1))-TMath::Erf((min-mean)/(sqrt(2.0)*sigma1)));
if (params->m_res_1() != 1.0)
fullrange += (1.0-fres)*1.0/2.0*(TMath::Erf((max-mean)/(sqrt(2.0)*sigma2))-TMath::Erf((min-mean)/(sqrt(2.0)*sigma2)));
double integral = (fres)*1.0/2.0*(TMath::Erf((mb-mean)/(sqrt(2.0)*sigma1))-TMath::Erf((ma-mean)/(sqrt(2.0)*sigma1)));
if (params->m_res_1() != 1.0)
integral += (1.0-fres)*1.0/2.0*(TMath::Erf((mb-mean)/(sqrt(2.0)*sigma2))-TMath::Erf((ma-mean)/(sqrt(2.0)*sigma2)));
result = integral/fullrange;
}
return result;
};
double bu2kstarmumu_pdf::m_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const
{
if (opts->only_angles || opts->only_mkpi) return 1.0;
double mmin=params->m_b.get_min(), mmax=params->m_b.get_max();
assert(meas.m <= mmax);
assert(meas.m >= mmin);
if(opts->cutsignalwindow){
assert(meas.m <= opts->m_min || meas.m >= opts->m_max);
}
double fm1 = params->fm_tau();
double tau1 = params->m_tau();
double tau2 = params->m_tau_2();
double lambda1 = params->m_lambda();
double lambda2 = params->m_lambda_2();
double result = 0.;
if(opts->fit_lambda){
if(lambda1 == 0.0) return 1.0;
double norm1 = 1.0/lambda1*(exp(mmax*lambda1) - exp(mmin*lambda1));
if(opts->cutsignalwindow){
norm1 -= 1.0/lambda1*(exp(opts->m_max*lambda1) - exp(opts->m_min*lambda1));
}
result = exp(meas.m*lambda1)/norm1;
if (fm1 != 1.0){
if(lambda2 == 0.0) return result;
result *= fm1;
double norm2 = 1.0/lambda2*(exp(mmax*lambda2) - exp(mmin*lambda2));
if(opts->cutsignalwindow)
norm2 -= 1.0/lambda2*(exp(opts->m_max*lambda2) - exp(opts->m_min*lambda2));
result += (1.0-fm1)*exp(meas.m*lambda2)/norm2;
}
}
else{
double norm1 = tau1*(exp(-mmin/tau1) - exp(-mmax/tau1));
if(opts->cutsignalwindow)
norm1 -= tau1*(exp(-opts->m_min/tau1) - exp(-opts->m_max/tau1));
result = exp(-meas.m/tau1)/norm1;
if (fm1 != 1.0){
result *= fm1;
double norm2 = tau2*(exp(-mmin/tau2) - exp(-mmax/tau2));
if(opts->cutsignalwindow)
norm2 -= tau2*(exp(-opts->m_min/tau2) - exp(-opts->m_max/tau2));
result += (1.0-fm1)*exp(-meas.m/tau2)/norm2;
}
}
return result;
};
double bu2kstarmumu_pdf::mkpi_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const
{
if(!(opts->fit_mkpi || opts->use_mkpi) || opts->only_angles || opts->only_Bmass)
return 1.0;
assert(meas.mkpi <= opts->mkpi_max);
assert(meas.mkpi >= opts->mkpi_min);
double result = 0.0;
double mkpi_min = opts->mkpi_min/1.0e+3;
double mkpi_max = opts->mkpi_max/1.0e+3;
double mkpi = meas.mkpi/1.0e+3;
double normbkg=0.0;
std::vector<double> ch_mkpi(5, 0.0);
ch_mkpi.at(0) = params->cbkgmkpi0();
ch_mkpi.at(1) = params->cbkgmkpi1();
ch_mkpi.at(2) = params->cbkgmkpi2();
ch_mkpi.at(3) = params->cbkgmkpi3();
ch_mkpi.at(4) = params->cbkgmkpi4();
std::vector<double> poly_mkpi(ch_mkpi.size(), 0.0);
chebyshev_to_poly(ch_mkpi, poly_mkpi);
std::vector<double> poly_corr_mkpi(ch_mkpi.size(), 0.0);
correct_poly(poly_mkpi, poly_corr_mkpi, mkpi_min, mkpi_max);
for(unsigned int i=0; i<poly_mkpi.size(); i++)
normbkg += poly_corr_mkpi.at(i)*(pow(mkpi_max,i+1)-pow(mkpi_min,i+1))/(i+1.0);
double probbkg = 0.0;
for (unsigned int i=0; i<ch_mkpi.size(); i++)
probbkg += poly_corr_mkpi.at(i)*pow(mkpi,i);
result = probbkg/normbkg;
return result;
};
double bu2kstarmumu_pdf::integral_m_bkg_prob(const bu2kstarmumu_parameters* params, double ma, double mb) const
{
double mmin=params->m_b.get_min(), mmax=params->m_b.get_max();
double bkg_frac = params->fm_tau();
double tau1 = params->m_tau();
double tau2 = params->m_tau_2();
double lambda1 = params->m_lambda();
double lambda2 = params->m_lambda_2();
double result = 0.0;
if(opts->fit_lambda){
assert(lambda1 <= 0.);
double fullrange1 = 1.0/lambda1*(exp(mmax*lambda1) - exp(mmin*lambda1));
if(opts->cutsignalwindow){
fullrange1 -= 1.0/lambda1*(exp(opts->m_max*lambda1) - exp(opts->m_min*lambda1));
}
double integral1 = 1.0/lambda1*(exp(mb*lambda1) - exp(ma*lambda1));
result = integral1/fullrange1;
if (bkg_frac != 1.0){
assert(lambda2 <= 0.);
result *= bkg_frac;
double fullrange2 = 1.0/lambda2*(exp(mmax*lambda2) - exp(mmin*lambda2));
if(opts->cutsignalwindow){
fullrange2 -= 1.0/lambda2*(exp(opts->m_max*lambda2) - exp(opts->m_min*lambda2));
}
double integral2 = 1.0/lambda2*(exp(mb*lambda2) - exp(ma*lambda2));
result += (1.0-bkg_frac)*integral2/fullrange2;
}
}
else{
double fullrange1 = tau1*(exp(-mmin/tau1) - exp(-mmax/tau1));
if(opts->cutsignalwindow){
fullrange1 -= tau1*(exp(opts->m_max/tau1) - exp(opts->m_min/tau1));
}
double integral1 = tau1*(exp(-ma/tau1) - exp(-mb/tau1));
result = integral1/fullrange1;
if (bkg_frac != 1.0){
result *= bkg_frac;
double fullrange2 = tau2*(exp(-mmin/tau2) - exp(-mmax/tau2));
if(opts->cutsignalwindow){
fullrange2 -= tau2*(exp(opts->m_max/tau2) - exp(opts->m_min/tau2));
}
double integral2 = tau2*(exp(-ma/tau2) - exp(-mb/tau2));
result += (1.0-bkg_frac)*integral2/fullrange2;
}
}
return result;
};
double bu2kstarmumu_pdf::angular_sig_prob(const bu2kstarmumu_parameters* params, const event& meas) const
{
//Don't when only fitting the mass
if (opts->only_Bmass || opts->only_mass2DFit || opts->only_mkpi) return 1.0;
double ctl = meas.costhetal;
double ctk = meas.costhetak;
double phi = meas.phi;
double result_swave = 0.0, result_pwave = 0.0;
double result = 0.0;
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
const double j1s = params->J1s();
const double j1c = params->J1c();
const double j2s = params->J2s();
const double j2c = params->J2c();
const double j3 = params->J3();
const double j4 = params->J4();
const double j5 = params->J5();
const double j6s = params->J6s();
const double j6c = params->J6c();
const double j7 = params->J7();
const double j8 = params->J8();
const double j9 = params->J9();
if (opts->full_angular) fj(ctl, ctk, phi, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
else folded_fj(ctl, ctk, phi, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
if (opts->fit_asymmetries && meas.cp_conjugate){
result_pwave += (f1*j1s + f2*j1c + f3*j2s + f4*j2c - f5*j3 - f6*j4
- f7*j5 - f8*j6s - f9*j6c - f10*j7 - f11*j8 - f12*j9);
}
else{
result_pwave += (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4
+ f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9);
}
if (opts->swave){
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0;
if(opts->full_angular) swave_fj(ctl, ctk, phi, f1, f2, f3, f4, f5, f6);
else folded_swave_fj(ctl, ctk, phi, f1, f2, f3, f4, f5, f6);
const double fs = params->FS();
const double js1 = params->SS1();
const double js2 = params->SS2();
const double js3 = params->SS3();
const double js4 = params->SS4();
const double js5 = params->SS5();
result_pwave *= (1.0 - fs);
result_swave += (fs*f1 + js1*f2 + js2*f3 + js3*f4 + js4*f5 + js5*f6);
if (opts->use_mkpi){
double prob_mkpi_pwave = mkpi_pwave_2(params, meas);
double prob_mkpi_swave = mkpi_swave_2(params, meas);
double prob_mkpi_re_swavepwave = mkpi_re_swavepwave(params, meas);
double prob_mkpi_im_swavepwave = mkpi_im_swavepwave(params, meas);
result_pwave *= prob_mkpi_pwave/mkpi_pwave_2_norm;
if (std::isnan(result_pwave)) spdlog::warn("prob_mkpi_pwave {0:f}",prob_mkpi_pwave );
result_swave = (fs*f1*prob_mkpi_swave/mkpi_swave_2_norm
+ (js1*f2 + js2*f3 + js3*f4)*prob_mkpi_re_swavepwave/mkpi_re_swavepwave_norm
+ (js4*f5 + js5*f6)*prob_mkpi_im_swavepwave/mkpi_im_swavepwave_norm);
}
if (std::isnan(result_swave) || std::isinf(result_swave)){
spdlog::warn("\tresult_swave {0:0.6f}",result_swave);
spdlog::warn("\t mkpi_swave_2_norm {0:0.6f}",mkpi_swave_2_norm);
spdlog::warn("\t mkpi_pwave_2_norm {0:0.6f}",mkpi_pwave_2_norm);
spdlog::warn("\t prob_mkpi_swave {0:0.6f}",mkpi_swave_2(params, meas));
spdlog::warn("\t prob_mkpi_re_swavepwave {0:0.6f}",mkpi_re_swavepwave(params, meas));
spdlog::warn("\t prob_mkpi_im_swavepwave {0:0.6f}",mkpi_im_swavepwave(params, meas));
result_swave = 1.0e-12;
}
}
result = result_pwave + result_swave;
if(std::isnan(result) || std::isinf(result)){
spdlog::error("result={0:0f}",result );
spdlog::error("pwave={0:0f}",result_pwave );
spdlog::error("swave={0:0f}",result_swave );
assert(0);
}
//spdlog::trace("result={0:0f}",result );
//spdlog::trace("pwave={0:0f}",result_pwave );
//spdlog::trace("swave={0:0f}",result_swave );
double eff = 1.0;
if (opts->use_angular_acc){
//for fixed effq2 the efficiency should have been calculated beforehand
//for floating q2 should recalculate as done below
//for event_norm, the efficiency should also have been calculated before
bool q2fixed = (params->eff_q2.get_step_size() == 0.0);
if ((q2fixed || opts->use_event_norm) && (!opts->mcweight)) { //use event norm is only used when convoluting the acceptance into the pdf
if (meas.weight > 0.0) eff = 1.0/meas.weight;
spdlog::trace("eff={0:f}",eff);
}
else{
double q2 = params->eff_q2();
eff = this->get_ang_eff(opts,q2,meas.costhetal_fa, meas.costhetal_fa,meas.phi_fa);
spdlog::trace("eff={0:f}",eff);
spdlog::trace("q2={0:f}",q2);
}
}
if(std::isnan(eff) || std::isinf(eff)){
spdlog::error("eff={0:f}",eff );
assert(0);
}
if(eff/fnorm_sig < 0.0){
spdlog::info("negative acc corr and/or normalizing");
spdlog::info("eff={0:f}",eff );
spdlog::info("norm={0:f}",fnorm_sig );
assert(eff > 0.0);
assert(fnorm_sig > 0.0);
}
//spdlog::trace("eff={0:f},\t fnorm_sig={1:f}",eff, fnorm_sig);
result *= eff/fnorm_sig;
return result;
};
double bu2kstarmumu_pdf::angular_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const{
//Don't when only fitting the mass
if (opts->only_Bmass || opts->only_mass2DFit || opts->only_mkpi) return 1.0;
//calculate bkg prob for all 2(4) possible inversed folding angles:
if(opts->fit_full_angular_bkg && !opts->flat_bkg && !opts->full_angular){
fcnc::event e[3];
fldr->invers_fold(&meas, &e[0], &e[1], &e[2]);
double p1, p2, p3, p4;
p1 = simple_angular_bkg_prob(params, meas);
p4 = simple_angular_bkg_prob(params, e[2]);
if(opts->folding > 0){
p2 = simple_angular_bkg_prob(params, e[0]);
p3 = simple_angular_bkg_prob(params, e[1]);
return p1+p2+p3+p4;
}
else{
return p1+p4;
}
}
else{
return simple_angular_bkg_prob(params, meas);
}
}
double bu2kstarmumu_pdf::simple_angular_bkg_prob(const bu2kstarmumu_parameters* params, const event& meas) const{
if (opts->only_Bmass) return 1.0;
if (opts->only_mass2DFit) return 1.0; //Don't when fitting only mass
double result = 1.0;
double costhetal = meas.costhetal;
double costhetak = meas.costhetak;
double phi = meas.phi;
npolynom npol(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
//determine weight/efficiency for this event
double eff = 1.0;
if(opts->use_angular_acc && opts->use_weighted_bkg){
//for fixed effq2 the efficiency should have been calculated beforehand
//for floating q2 should recalculate as done below
//for event_norm, the efficiency should also have been calculated before
bool q2fixed = (params->eff_q2.get_step_size() == 0.0);
if ((q2fixed || opts->use_event_norm) && (!opts->mcweight)){//use event norm is only used when convoluting the acceptance into the pdf
if (meas.weight > 0.0) eff = 1.0/meas.weight;
}
else eff = this->get_ang_eff(opts,meas,true);
}
if (opts->flat_bkg){ //for a flat background, only - if at all - the angular acceptance correction is needed
result *= eff/fnorm_bkg;
}
else{ //Chebyshev background
//read angular bkg coefficients from parameter object
std::vector<double> ch_ctl = init_ch_ctl(params);
std::vector<double> ch_ctk = init_ch_ctk(params);
std::vector<double> ch_phi = init_ch_phi(params);
//reserve vectors for the (non-)Chebyshev polynomial coefficients
std::vector<double> poly_ctl(ch_ctl.size(), 0.0);
std::vector<double> poly_ctk(ch_ctk.size(), 0.0);
std::vector<double> poly_phi(ch_phi.size(), 0.0);
//these temporary coefficients are defined in the range from -1 to +1
std::vector<double> temp_ctl(ch_ctl.size(), 0.0);
std::vector<double> temp_ctk(ch_ctk.size(), 0.0);
std::vector<double> temp_phi(ch_phi.size(), 0.0);
//transform chebyshev polynomial coefficients to Standard Form polynomials
chebyshev_to_poly(ch_ctl, temp_ctl);
chebyshev_to_poly(ch_ctk, temp_ctk);
chebyshev_to_poly(ch_phi, temp_phi);
//stretch or squeeze the definition of the polynomials in phi
//for (unsigned int i=0; i<poly_ctk.size(); i++) poly_ctk.at(i) = temp_ctk.at(i);
for (unsigned int i=0; i<poly_ctk.size(); i++) poly_ctk.at(i) = temp_ctk.at(i);
for (unsigned int i=0; i<poly_ctl.size(); i++) poly_ctl.at(i) = temp_ctl.at(i);
correct_poly(temp_phi, poly_phi, -TMath::Pi(), +TMath::Pi());
//add up all power terms of the polynomials
double p_costhetal = 0.0;
for (unsigned int i=0; i<ch_ctl.size(); i++) p_costhetal += poly_ctl.at(i)*pow(costhetal,i);
double p_costhetak = 0.0;
for (unsigned int i=0; i<ch_ctk.size(); i++) p_costhetak += poly_ctk.at(i)*pow(costhetak,i);
double p_phi = 0.0;
for (unsigned int i=0; i<ch_phi.size(); i++) p_phi += poly_phi.at(i)*pow(phi,i);
//multiply the indidivual angular dimensions. divide by the normalisation to obtain the angular background probability
result = p_costhetal*p_costhetak*p_phi*eff/fnorm_bkg;
}
return result;
};
double bu2kstarmumu_pdf::integral_bkg(const bu2kstarmumu_parameters* params, double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b) const
{
if (opts->flat_bkg){
return (ctl_b-ctl_a)*(ctk_b-ctk_a)*(phi_b-phi_a); // RANGE_3D() for non-folded
}
else{
std::vector<double> ch_ctl = init_ch_ctl(params);
std::vector<double> ch_ctk = init_ch_ctk(params);
std::vector<double> ch_phi = init_ch_phi(params);
std::vector<double> poly_ctl(ch_ctl.size(), 0.0);
std::vector<double> poly_ctk(ch_ctk.size(), 0.0);
std::vector<double> poly_phi(ch_phi.size(), 0.0);
std::vector<double> temp_ctl(ch_ctl.size(), 0.0);
std::vector<double> temp_ctk(ch_ctk.size(), 0.0);
std::vector<double> temp_phi(ch_phi.size(), 0.0);
chebyshev_to_poly(ch_ctl, temp_ctl);
chebyshev_to_poly(ch_ctk, temp_ctk);
chebyshev_to_poly(ch_phi, temp_phi);
for (unsigned int i=0; i<poly_ctk.size(); i++) poly_ctk.at(i) = temp_ctk.at(i);
for (unsigned int i=0; i<poly_ctl.size(); i++)poly_ctl.at(i) = temp_ctl.at(i);
correct_poly(temp_phi, poly_phi, -TMath::Pi(), +TMath::Pi());
double norm_costhetal = 0.0;
for (unsigned int i=0; i<poly_ctl.size(); i++)norm_costhetal += poly_ctl.at(i)*(pow(ctl_b,i+1)-pow(ctl_a,i+1))/(i+1.0);
double norm_costhetak = 0.0;
for (unsigned int i=0; i<poly_ctk.size(); i++)norm_costhetak += poly_ctk.at(i)*(pow(ctk_b,i+1)-pow(ctk_a,i+1))/(i+1.0);
double norm_phi = 0.0;
for (unsigned int i=0; i<poly_phi.size(); i++)norm_phi += poly_phi.at(i)*(pow(phi_b,i+1)-pow(phi_a,i+1))/(i+1.0);
return norm_costhetal*norm_costhetak*norm_phi;
}
}
double bu2kstarmumu_pdf::integral_bkg_chebyshev(const bu2kstarmumu_parameters* params, double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b, double q2) const{
double result = 0.0;
if (opts->flat_bkg){
double c = 1.0;
//new improved non-factorizing acceptance
npolynom npol = npolynom(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
double chijk = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
const unsigned int bin = npol.get_bin_in4D(h,i,j,k);
chijk = c*pow(q2, h)*coeffs_eff_4d.at(bin);
result += chijk
* (pow(ctl_b, i+1)/(i+1.0) - pow(ctl_a, i+1)/(i+1.0))
* (pow(ctk_b, j+1)/(j+1.0) - pow(ctk_a, j+1)/(j+1.0))
* (pow(phi_b, k+1)/(k+1.0) - pow(phi_a, k+1)/(k+1.0));
}
}
}
}
}
else{
//calculate integral of background with efficiency correction
std::vector<double> ch_ctl = init_ch_ctl(params);
std::vector<double> ch_ctk = init_ch_ctk(params);
std::vector<double> ch_phi = init_ch_phi(params);
std::vector<double> poly_ctl(ch_ctl.size(), 0.0);
std::vector<double> poly_ctk(ch_ctk.size(), 0.0);
std::vector<double> poly_phi(ch_phi.size(), 0.0);
std::vector<double> temp_ctl(ch_ctl.size(), 0.0);
std::vector<double> temp_ctk(ch_ctk.size(), 0.0);
std::vector<double> temp_phi(ch_phi.size(), 0.0);
chebyshev_to_poly(ch_ctl, temp_ctl);
chebyshev_to_poly(ch_ctk, temp_ctk);
chebyshev_to_poly(ch_phi, temp_phi);
for (unsigned int i=0; i<poly_ctk.size(); i++) poly_ctk.at(i) = temp_ctk.at(i);
for (unsigned int i=0; i<poly_ctl.size(); i++) poly_ctl.at(i) = temp_ctl.at(i);
correct_poly(temp_phi, poly_phi, -TMath::Pi(), +TMath::Pi());
double c = 1.0;
//new improved non-factorizing acceptance
//double q2 = params->eff_q2();
npolynom npol(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
double chijk = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
const unsigned int bin = npol.get_bin_in4D(h,i,j,k);
chijk = c*pow(q2, h)*coeffs_eff_4d.at(bin);
for (unsigned int l=0; l<poly_ctl.size(); l++){
for (unsigned int m=0; m<poly_ctk.size(); m++){
for (unsigned int n=0; n<poly_phi.size(); n++){
result += chijk * poly_ctl.at(l) * poly_ctk.at(m) * poly_phi.at(n)
* (pow(ctl_b, i+l+1) - pow(ctl_a, i+l+1)) / (i+l+1.0)
* (pow(ctk_b, j+m+1) - pow(ctk_a, j+m+1)) / (j+m+1.0)
* (pow(phi_b, k+n+1) - pow(phi_a, k+n+1)) / (k+n+1.0);
}
}
}
}
}
}
}
}
return result;
}
double bu2kstarmumu_pdf::angular_bkg_norm(const bu2kstarmumu_parameters* params, bool weighted){
npolynom npol(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
//bkg norm
if(weighted){ //weighted background pdf
double q2 = params->eff_q2();
if (opts->full_angular){
if (opts->flat_bkg){
double integral = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
integral += pow(q2, h)*coeffs_eff_4d.at(npol.get_bin_in4D(h,i,j,k))
* integrate_x_to_n(CTL_MIN,CTL_MAX,i)
* integrate_x_to_n(CTK_MIN,CTK_MAX,j)
* integrate_x_to_n(PHI_MIN,PHI_MAX,k);
}
}
}
}
return integral;
}
else{
return integral_bkg_chebyshev(params, CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX, q2);
}
}
else{ //folded angles
if (opts->flat_bkg){
double integral = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
integral += pow(q2, h)*coeffs_eff_4d.at(npol.get_bin_in4D(h,i,j,k))
* integrate_x_to_n(opts->ctl_min,opts->ctl_max,i)
* integrate_x_to_n(opts->ctk_min,opts->ctk_max,j)
* integrate_x_to_n(opts->phi_min,opts->phi_max,k);
}
}
}
}
return integral;
}
else{
if(opts->fit_full_angular_bkg){
return integral_bkg_chebyshev(params, CTL_MIN, CTL_MAX, opts->ctk_min, opts->ctk_max, PHI_MIN, PHI_MAX, q2); //TODO: check if that makes sense (https://gitlab.cern.ch/LHCb-RD/ewp-Bplus2Kstmumu-AngAna/-/merge_requests/21)
}
else{
return integral_bkg_chebyshev(params, opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max, q2);
}
}
}
}
else{ //non weighted background pdf
if (opts->full_angular){
if (opts->flat_bkg) return RANGE_3D();
else return integral_bkg(params, CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX);
}
else{ //folded angles
if (opts->flat_bkg){
return (opts->ctk_max - opts->ctk_min) * (opts->ctl_max - opts->ctl_min) * (opts->phi_max - opts->phi_min);
}
else{
if(opts->fit_full_angular_bkg){
return integral_bkg(params, CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX);
}
else{
return integral_bkg(params, opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max);
}
}
}
}
}
double bu2kstarmumu_pdf::angular_sig_norm(const bu2kstarmumu_parameters* params, bool weighted){
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
double fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0, fs6=0.0;
const double j1s = params->J1s();
const double j1c = params->J1c();
const double j2s = params->J2s();
const double j2c = params->J2c();
const double j3 = params->J3();
const double j4 = params->J4();
const double j5 = params->J5();
const double j6s = params->J6s();
const double j6c = params->J6c();
const double j7 = params->J7();
const double j8 = params->J8();
const double j9 = params->J9();
const double fs = params->FS();
const double js1 = params->SS1();
const double js2 = params->SS2();
const double js3 = params->SS3();
const double js4 = params->SS4();
const double js5 = params->SS5();
//sig norm
double q2 = params->eff_q2();
if (opts->full_angular){
//get signal pdf
if (weighted) integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, q2);
else integrated_fj_noacc(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
//swave
if(opts->swave){
if (weighted) swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
else swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
//if asymmetries are fitted:
if(opts->fit_asymmetries){
if(opts->swave) return (1.0-fs) * (f1*j1s + f2*j1c + f3*j2s + f4*j2c) + (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
else return f1*j1s + f2*j1c + f3*j2s + f4*j2c;
}
if (opts->swave){
return (1.0-fs) * (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9)
+ (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
}
//signal pdf without swave, without asymmetries and without m(Kpi)
return f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9;
}
else{ //folded angles
if (weighted) folded_integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, q2);
else folded_integrated_fj_noacc(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
if(opts->swave){
if (weighted) folded_swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
else folded_swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
if(opts->fit_asymmetries){
if(opts->swave) return (1-fs) * (f1*j1s + f2*j1c + f3*j2s + f4*j2c) + (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
else return (f1*j1s + f2*j1c + f3*j2s + f4*j2c);
}
if(opts->swave){
return (1-fs) * (f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9)
+ (fs1*fs + fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5);
}
//signal pdf without swave, without asymmetries and without m(Kpi)
return f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9;
}
}
double bu2kstarmumu_pdf::angular_swave_norm(const bu2kstarmumu_parameters* params, bool weighted){
double fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0, fs6=0.0;
const double fs = params->FS();
//sig norm
double q2 = params->eff_q2();
if (opts->full_angular){
if (weighted) swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
else swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
else{ //folded angles
if (weighted) folded_swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
else folded_swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
return fs1*fs;
}
double bu2kstarmumu_pdf::angular_pwave_norm(const bu2kstarmumu_parameters* params, bool weighted){
double f1=0.0, f2=0.0, f3=0.0, f4=0.0, f5=0.0, f6=0.0,
f7=0.0, f8=0.0, f9=0.0, f10=0.0, f11=0.0, f12=0.0;
const double j1s = params->J1s();
const double j1c = params->J1c();
const double j2s = params->J2s();
const double j2c = params->J2c();
const double j3 = params->J3();
const double j4 = params->J4();
const double j5 = params->J5();
const double j6s = params->J6s();
const double j6c = params->J6c();
const double j7 = params->J7();
const double j8 = params->J8();
const double j9 = params->J9();
//sig norm
double q2 = params->eff_q2();
if (opts->full_angular){
if (weighted) integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, q2);
else integrated_fj_noacc(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
}
else{ //folded angles
if (weighted) folded_integrated_fj_chebyshev(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, q2);
else folded_integrated_fj_noacc(f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
}
//if asymmetries are fitted:
if(opts->fit_asymmetries)return f1*j1s + f2*j1c + f3*j2s + f4*j2c;
else return f1*j1s + f2*j1c + f3*j2s + f4*j2c + f5*j3 + f6*j4 + f7*j5 + f8*j6s + f9*j6c + f10*j7 + f11*j8 + f12*j9;
}
double bu2kstarmumu_pdf::angular_pswave_norm(const bu2kstarmumu_parameters* params, bool weighted){
double fs1=0.0, fs2=0.0, fs3=0.0, fs4=0.0, fs5=0.0, fs6=0.0;
const double js1 = params->SS1();
const double js2 = params->SS2();
const double js3 = params->SS3();
const double js4 = params->SS4();
const double js5 = params->SS5();
//sig norm
double q2 = params->eff_q2();
if (opts->full_angular){
if (weighted) swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
else swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
else{ //folded angles
if (weighted) folded_swave_integrated_fj_chebyshev(fs1, fs2, fs3, fs4, fs5, fs6, q2);
else folded_swave_integrated_fj_noacc(fs1, fs2, fs3, fs4, fs5, fs6);
}
return fs2*js1 + fs3*js2 + fs4*js3 + fs5*js4 + fs6*js5;
}
void bu2kstarmumu_pdf::update_cached_efficiencies(const bu2kstarmumu_parameters* params, std::vector<event>* events)
{
if (opts->update_efficiencies){
spdlog::info("Updating cached efficiencies");
//assert(0);
//three ways this is usefull
//-q2_eff is fixed and you want to avoid recalculating this every time
//-weighted fit to determine the event weight that needs to be used
//-when using per event norm this can also be cached
npolynom npol(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
double q2 = params->eff_q2();
unsigned int N = events->size();
for (unsigned int e=0; e<N; e++){
double eff = 0.0;
event& meas = events->at(e);
if (opts->weighted_fit || (opts->use_angular_acc && opts->use_event_norm)){
q2 = meas.q2;
}
eff = this->get_ang_eff(opts,q2, meas.costhetal_fa, meas.costhetak_fa, meas.phi_fa);
if (eff > EFF_CUTOFF){
if (opts->multiply_eff) meas.weight *= 1.0/eff;
else meas.weight = 1.0/eff;
}
else{
spdlog::warn("Event with too low efficiency: {0:f}", eff);
spdlog::warn("\t q2: {0:f} \tq2_eff {1:f} \tctl: {2:f} \tctk: {3:f} \tphi: {4:f} ",
meas.q2, q2, meas.costhetal_fa, meas.costhetak_fa , meas.phi_fa);
if (opts->multiply_eff) meas.weight *= 1.0;
else meas.weight = 1.0;
}
}
return;
}
else return;
}
////////////////////////
/// Full angular ///
////////////////////////
//--------------------//
// S-wave
//--------------------//
void bu2kstarmumu_pdf::swave_fj(double ctl, double ctk, double phi,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const{
f1 = C_SWAVE * sintheta2(ctl);
f2 = C_SWAVE * sintheta2(ctl)*ctk;
f3 = C_SWAVE * sin2theta(ctl)*sintheta_abs(ctk)*cos(phi);
f4 = C_SWAVE * sintheta_abs(ctl)*sintheta_abs(ctk)*cos(phi);
f5 = C_SWAVE * sintheta_abs(ctl)*sintheta_abs(ctk)*sin(phi);
f6 = C_SWAVE * sin2theta(ctl)*sintheta_abs(ctk)*sin(phi);
}
void bu2kstarmumu_pdf::swave_integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const{
f1 = 1.0;
f2 = 0.0;
f3 = 0.0;
f4 = 0.0;
f5 = 0.0;
f6 = 0.0;
swave_integrated_fj_noacc(opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max,
f1, f2, f3, f4, f5, f6);
}
void bu2kstarmumu_pdf::swave_integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const{
f1 = C_SWAVE * integrate_s_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f2 = C_SWAVE * integrate_s_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f3 = C_SWAVE * integrate_s_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f4 = C_SWAVE * integrate_s_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f5 = C_SWAVE * integrate_s_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f6 = C_SWAVE * integrate_s_f6(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
return;
}
void bu2kstarmumu_pdf::swave_integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const //TODO: modify this according to the normDebug branch
{
f1 = 0.0; f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0; f6 = 0.0;
npolynom npol(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
double chijk = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
const unsigned int bin = npol.get_bin_in4D(h,i,j,k);
chijk = C_SWAVE*pow(q2, h)*coeffs_eff_4d.at(bin);
f1 += chijk * integrate_s_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f2 += chijk * integrate_s_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f3 += chijk * integrate_s_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f4 += chijk * integrate_s_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f5 += chijk * integrate_s_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f6 += chijk * integrate_s_f6(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
}
}
}
}
return;
};
void bu2kstarmumu_pdf::swave_integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const{
f1 = 0.0; f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0; f6 = 0.0;
swave_integrated_fj_chebyshev(CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX,
f1, f2, f3, f4, f5, f6, q2);
return;
};
//--------------------//
// P-wave
//--------------------//
void bu2kstarmumu_pdf::fj(double ctl, double ctk, double phi,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const{
f1 = C_PWAVE * sintheta2(ctk); //j1s = 3/4*(1-FL)
f2 = C_PWAVE * costheta2(ctk); //j1c = FL
f3 = C_PWAVE * sintheta2(ctk)*cos2theta(ctl); //j2s = 1/4*(1-FL)
f4 = C_PWAVE * costheta2(ctk)*cos2theta(ctl); //j2c = -FL
f5 = C_PWAVE * sintheta2(ctk)*sintheta2(ctl)*cos(2.0*phi); //j3
f6 = C_PWAVE * sin2theta(ctk)*sin2theta(ctl)*cos(phi); //j4
f7 = C_PWAVE * sin2theta(ctk)*sintheta_abs(ctl)*cos(phi); //j5
f8 = C_PWAVE * sintheta2(ctk)*ctl; //j6s
f9 = C_PWAVE * costheta2(ctk)*ctl; //j6c
f10 = C_PWAVE * sin2theta(ctk)*sintheta_abs(ctl)*sin(phi); //j7
f11 = C_PWAVE * sin2theta(ctk)*sin2theta(ctl)*sin(phi); //j8
f12 = C_PWAVE * sintheta2(ctk)*sintheta2(ctl)*sin(2.0*phi); //j9
};
void bu2kstarmumu_pdf::integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const{
//Normalization of the p-wave in principle
f1 = 0.0; f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0; f6 = 0.0;
f7 = 0.0; f8 = 0.0; f9 = 0.0; f10 = 0.0; f11 = 0.0; f12 = 0.0;
integrated_fj_noacc(CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX,
f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
return;
};
void bu2kstarmumu_pdf::integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const{
f1 = C_PWAVE * integrate_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f2 = C_PWAVE * integrate_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f3 = C_PWAVE * integrate_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f4 = C_PWAVE * integrate_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f5 = C_PWAVE * integrate_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f6 = C_PWAVE * integrate_f6(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f7 = C_PWAVE * integrate_f7(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f8 = C_PWAVE * integrate_f8(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f9 = C_PWAVE * integrate_f9(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f10 = C_PWAVE * integrate_f10(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f11 = C_PWAVE * integrate_f11(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f12 = C_PWAVE * integrate_f12(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
return;
};
void bu2kstarmumu_pdf::integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const{
f1 = 0.0;
f2 = 0.0;
f3 = 0.0;
f4 = 0.0;
f5 = 0.0;
f6 = 0.0;
f7 = 0.0;
f8 = 0.0;
f9 = 0.0;
f10 = 0.0;
f11 = 0.0;
f12 = 0.0;
//new improved non-factorizing acceptance
npolynom npol = npolynom(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
const unsigned int bin = npol.get_bin_in4D(h,i,j,k);
double chijk = C_PWAVE*coeffs_eff_4d.at(bin)*pow(q2, h); //9/32pi * angCorr(bin) * q2^h (not integrating over q2)
f1 += chijk * integrate_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f2 += chijk * integrate_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f3 += chijk * integrate_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f4 += chijk * integrate_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f5 += chijk * integrate_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f6 += chijk * integrate_f6(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f7 += chijk * integrate_f7(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f8 += chijk * integrate_f8(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f9 += chijk * integrate_f9(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f10 += chijk * integrate_f10(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f11 += chijk * integrate_f11(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f12 += chijk * integrate_f12(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
}
}
}
}
}
void bu2kstarmumu_pdf::integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6,
double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const{
f1 = 0.0; f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0; f6 = 0.0;
f7 = 0.0; f8 = 0.0; f9 = 0.0; f10 = 0.0; f11 = 0.0; f12 = 0.0;
integrated_fj_chebyshev(CTL_MIN, CTL_MAX, CTK_MIN, CTK_MAX, PHI_MIN, PHI_MAX,
f1, f2, f3, f4, f5, f6,
f7, f8, f9, f10, f11, f12, q2);
return;
};
/////////////////////
///Angular folding///
/////////////////////
//S-wave
void bu2kstarmumu_pdf::folded_swave_fj(double ctl, double ctk, double phi,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const
{
double c = C_SWAVE*2.0; // folding in phi
if(opts->folding > 0) c *= 2.0; //folding in cos(thetaL)
f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0; f6 = 0.0;
f1 = c * sintheta2(ctl);
if(opts->folding != 4) f2 = c * sintheta2(ctl)*(ctk);
if(opts->folding == 1) f3 = c * sintheta_abs(ctk)*sin2theta(ctl)*cos(phi);
if(opts->folding == 2) f4 = c * sintheta_abs(ctk)*sintheta_abs(ctl)*cos(phi);
if(opts->folding > 2) f5 = c * sintheta_abs(ctk)*sintheta_abs(ctl)*sin(phi);
//f6 always zero no matter what folding
}
void bu2kstarmumu_pdf::folded_swave_integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const
{
f1 = 1.0;
f2 = 0.0;
f3 = 0.0;
f4 = 0.0;
f5 = 0.0;
f6 = 0.0;
folded_swave_integrated_fj_noacc(opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max,
f1, f2, f3, f4, f5, f6);
}
void bu2kstarmumu_pdf::folded_swave_integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6) const{
double c = C_SWAVE*2.0; // folding in phi
if(opts->folding > 0) c *= 2.0; //folding in cos(thetaL)
f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0;
f1 = c * integrate_s_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
if(opts->folding != 4) f2 = c * integrate_s_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
if(opts->folding == 1) f3 = c * integrate_s_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
if(opts->folding == 2) f4 = c * integrate_s_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
if(opts->folding > 2) f5 = c * integrate_s_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f6 = 0.0;
return;
}
void bu2kstarmumu_pdf::folded_swave_integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const{
double c = C_SWAVE*2.0; // folding in phi
if(opts->folding > 0) c *= 2.0; //folding in cos(thetaL)
f1 = 0.0; f2 = 0.0; f3 = 0.0; f4 = 0.0; f5 = 0.0; f6 = 0.0;
npolynom npol = npolynom(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
double chijk = 0.0;
for (unsigned int h = 0; h<npol.q2; h++){
for (unsigned int i = 0; i<npol.ctl; i++){
for (unsigned int j = 0; j<npol.ctk; j++){
for (unsigned int k = 0; k<npol.phi; k++){
const unsigned int bin = npol.get_bin_in4D(h,i,j,k);
chijk = c*pow(q2, h)*coeffs_eff_4d.at(bin);
f1 += chijk * integrate_s_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
if(opts->folding != 4) f2 += chijk * integrate_s_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
if(opts->folding == 1) f3 += chijk * integrate_s_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
if(opts->folding == 2) f4 += chijk * integrate_s_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
if(opts->folding > 2) f5 += chijk * integrate_s_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f6 = 0.0;
}
}
}
}
}
void bu2kstarmumu_pdf::folded_swave_integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double q2) const{
folded_swave_integrated_fj_chebyshev(opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max,
f1, f2, f3, f4, f5, f6, q2);
return;
}
//P-wave
void bu2kstarmumu_pdf::folded_fj(double ctl, double ctk, double phi, double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const
{
//reset all coefficients that might not be calculated due to folding!
f6 = 0.;
f7 = 0.;
f8 = 0.;
f9 = 0.;
f10 = 0.;
f11 = 0.;
f12 = 0.;
//constant coefficient
double c = 2.0*C_PWAVE; //multiply by 2 due to folding in phi
if(opts->folding > 0) c *= 2.0; //folding in cos(thetaL)
f1 = c * sintheta2(ctk); //j1s = 3/4*(1-FL)
f2 = c * costheta2(ctk); //j1c = FL
f3 = c * sintheta2(ctk)*cos2theta(ctl); //j2s = 1/4*(1-FL)
f4 = c * costheta2(ctk)*cos2theta(ctl); //j2c = -FL
f5 = c * sintheta2(ctk)*sintheta2(ctl)*cos(2.0*phi); //j3
if(opts->folding == 1) f6 = c * sin2theta(ctk)*sin2theta(ctl)*cos(phi); //j4
if(opts->folding == 2) f7 = c * sin2theta(ctk)*sintheta_abs(ctl)*cos(phi); //j5
if(opts->folding == 0){
f8 = c * sintheta2(ctk)*ctl; //j6s
f9 = c * costheta2(ctk)*ctl; //j6c
}
if(opts->folding == 3) f10 = c * sin2theta(ctk)*sintheta_abs(ctl)*sin(phi); //j7
if(opts->folding == 4) f11 = c * sin2theta(ctk)*sin2theta(ctl)*sin(phi); //j8
if(opts->folding == 0) f12 = c * sintheta2(ctk)*sintheta2(ctl)*sin(2.0*phi); //j9
};
void bu2kstarmumu_pdf::folded_integrated_fj_noacc(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const
{
folded_integrated_fj_noacc(opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max,
f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12);
};
void bu2kstarmumu_pdf::folded_integrated_fj_noacc(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12) const
{
if(opts->full_angular){
spdlog::error("Cannot use folded integrals for full angular fit");
assert(0);
}
f1 = 0.0;
f2 = 0.0;
f3 = 0.0;
f4 = 0.0;
f5 = 0.0;
f6 = 0.0;
f7 = 0.0;
f8 = 0.0;
f9 = 0.0;
f10 = 0.0;
f11 = 0.0;
f12 = 0.0;
//copied from full angular integral:
double c = C_PWAVE*2.; //folding in phi
if(opts->folding > 0)c *= 2.; //folding in cos(theta_L)
f1 = c * integrate_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f2 = c * integrate_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f3 = c * integrate_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f4 = c * integrate_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f5 = c * integrate_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
if(opts->folding == 0){
f8 = c * integrate_f8(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f9 = c * integrate_f9(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
f12 = c * integrate_f12(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
}
else if(opts->folding == 1) f6 = c * integrate_f6(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
else if(opts->folding == 2) f7 = c * integrate_f7(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
else if(opts->folding == 3) f10 = c * integrate_f10(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
else if(opts->folding == 4) f11 = c * integrate_f11(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, 0, 0, 0);
return;
};
void bu2kstarmumu_pdf::folded_integrated_fj_chebyshev(double ctl_a, double ctl_b, double ctk_a, double ctk_b, double phi_a, double phi_b,
double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const{
if(opts->full_angular){
spdlog::error("Cannot use folded integrals for full angular fit");
assert(0);
}
double c = C_PWAVE* 2.; //folding in phi
if(opts->folding > 0)c *= 2.; //folding in cos(theta_L)
f1 = 0.0;
f2 = 0.0;
f3 = 0.0;
f4 = 0.0;
f5 = 0.0;
f6 = 0.0;
f7 = 0.0;
f8 = 0.0;
f9 = 0.0;
f10 = 0.0;
f11 = 0.0;
f12 = 0.0;
//new improved non-factorizing acceptance
npolynom npol = npolynom(opts);
assert(coeffs_eff_4d.size() == npol.getSize());
double chijk = 0.0;
for (unsigned int h = 0; h<npol.q2; h++)
for (unsigned int i = 0; i<npol.ctl; i++)
for (unsigned int j = 0; j<npol.ctk; j++)
for (unsigned int k = 0; k<npol.phi; k++)
{
const unsigned int bin = npol.get_bin_in4D(h,i,j,k);
chijk = c*pow(q2, h)*coeffs_eff_4d.at(bin);
f1 += chijk*integrate_f1(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f2 += chijk*integrate_f2(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f3 += chijk*integrate_f3(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f4 += chijk*integrate_f4(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f5 += chijk*integrate_f5(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
if(opts->folding == 0){
f8 += chijk*integrate_f8(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f9 += chijk*integrate_f9(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
f12 += chijk*integrate_f12(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
}
else if(opts->folding == 1) f6 += chijk*integrate_f6(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
else if(opts->folding == 2) f7 += chijk*integrate_f7(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
else if(opts->folding == 3) f10 += chijk*integrate_f10(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
else if(opts->folding == 4) f11 += chijk*integrate_f11(ctl_a, ctl_b, ctk_a, ctk_b, phi_a, phi_b, i, j, k);
}
};
void bu2kstarmumu_pdf::folded_integrated_fj_chebyshev(double& f1, double& f2, double& f3, double& f4, double& f5, double& f6, double& f7, double& f8, double& f9, double& f10, double& f11, double& f12, double q2) const{
folded_integrated_fj_chebyshev(opts->ctl_min, opts->ctl_max, opts->ctk_min, opts->ctk_max, opts->phi_min, opts->phi_max,
f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, q2);
return;
};
void bu2kstarmumu_pdf::calculate_sweights(const bu2kstarmumu_parameters* parms, std::vector<event>* ev) const
{
assert(opts->extended_ml);
parameter* psig = parms->get_parameter("n_sig");
parameter* pbkg = parms->get_parameter("n_bkg");
if (psig != 0 && pbkg != 0){
double nsig = psig->get_value();
double nbkg = pbkg->get_value();
double sweight = 0.0, prob_sig_m = 0.0, prob_bkg_m = 0.0;
//obtain inverse variant matrix (V_{nj}^{-1})
double inv_cov_matrix[2*2] = {0.0, 0.0, 0.0, 0.0};
for (unsigned int i = 0; i<ev->size(); i++){
const event& meas = ev->at(i);
prob_sig_m = m_sig_prob(parms, meas);
prob_bkg_m = m_bkg_prob(parms, meas);
inv_cov_matrix[0] += prob_sig_m*prob_sig_m/sqr(nsig*prob_sig_m + nbkg*prob_bkg_m);
inv_cov_matrix[1] += prob_sig_m*prob_bkg_m/sqr(nsig*prob_sig_m + nbkg*prob_bkg_m);
inv_cov_matrix[2] += prob_bkg_m*prob_sig_m/sqr(nsig*prob_sig_m + nbkg*prob_bkg_m);
inv_cov_matrix[3] += prob_bkg_m*prob_bkg_m/sqr(nsig*prob_sig_m + nbkg*prob_bkg_m);
}
spdlog::info("Inverse covariance matrix");
spdlog::info("{0:0.8f} {1:.8f}", inv_cov_matrix[0],inv_cov_matrix[1] );
spdlog::info("{0:0.8f} {1:.8f}", inv_cov_matrix[2],inv_cov_matrix[3] );
//invert inserve matrix to get V_{nj}
double cov_matrix[2*2] = {0.0, 0.0, 0.0, 0.0};
cov_matrix[0] = +inv_cov_matrix[3]/(inv_cov_matrix[0]*inv_cov_matrix[3]-inv_cov_matrix[1]*inv_cov_matrix[2]);
cov_matrix[1] = -inv_cov_matrix[1]/(inv_cov_matrix[0]*inv_cov_matrix[3]-inv_cov_matrix[1]*inv_cov_matrix[2]);
cov_matrix[2] = -inv_cov_matrix[2]/(inv_cov_matrix[0]*inv_cov_matrix[3]-inv_cov_matrix[1]*inv_cov_matrix[2]);
cov_matrix[3] = +inv_cov_matrix[0]/(inv_cov_matrix[0]*inv_cov_matrix[3]-inv_cov_matrix[1]*inv_cov_matrix[2]);
spdlog::info("Covariance matrix");
spdlog::info("{0:0.8f} {1:.8f}", cov_matrix[0], cov_matrix[1] );
spdlog::info("{0:0.8f} {1:.8f}", cov_matrix[2], cov_matrix[3] );
//get sWeights for all events
for (unsigned int i = 0; i<ev->size(); i++){
const event& meas = ev->at(i);
prob_sig_m = m_sig_prob(parms, meas);
prob_bkg_m = m_bkg_prob(parms, meas);
sweight = (prob_sig_m * cov_matrix[0] + prob_bkg_m * cov_matrix[1]) / ( prob_sig_m * nsig + prob_bkg_m * nbkg );
ev->at(i).sweight = sweight;
}
}
};
void bu2kstarmumu_pdf::save_coeffs_eff_phsp_4d(){
std::ofstream ostr(coeffs_eff_phsp_4D(opts));
for(unsigned int c = 0; c < this->coeffs_eff_4d.size(); c++){
ostr << std::setprecision(15) << std::scientific << this->coeffs_eff_4d.at(c) << "\n";
}
spdlog::info("[SAVE]\t\tSaved a total of {0:d} angular acceptance correction coefficients",this->coeffs_eff_4d.size());
}
std::vector<double> bu2kstarmumu_pdf::read_coeffs_eff_phsp_4d(){
unsigned int nCoeffs = opts->eff_order_q2 * opts->eff_order_costhetal * opts->eff_order_costhetak * opts->eff_order_phi;
std::string filename = coeffs_eff_phsp_4D(opts);
std::ifstream file(filename);
if(!file.is_open())spdlog::error("File " + filename + " not found to load coefficients!");
else spdlog::info("Loading a total of {0:d} coefficients from file " + filename + ".",nCoeffs);
std::string line;
double coeff;
std::vector<double>coeffs;
assert(file.is_open());
while(1){
getline(file, line);
if(file.eof())break;
std::istringstream istr(line);
istr >> coeff;
coeffs.push_back(coeff);
}
if(coeffs.size() != nCoeffs){
spdlog::critical("Number of coefficients loaded from file ({0:d}) does not match the number of configurated coefficients ({1:d}). Abort", coeffs.size(), nCoeffs);
assert(0);
}
return coeffs;
}
void bu2kstarmumu_pdf::load_coeffs_eff_phsp_4d(){
this->coeffs_eff_4d = read_coeffs_eff_phsp_4d();
}
Reference in New Issue View Git Blame Copy Permalink