diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc
index a064c4d..2e9f967 100644
--- a/examples/Test_Dpipi.cc
+++ b/examples/Test_Dpipi.cc
@@ -1,324 +1,327 @@
#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
#include "TFile.h"
#include "TH2.h"
#include "TRandom.h"
#include "TString.h"
#include "TSystem.h"
#include "TF1.h"
#include "TCanvas.h"
#include "LauDaughters.hh"
#include "LauDecayTimePdf.hh"
#include "LauEffModel.hh"
#include "LauIsobarDynamics.hh"
#include "LauMagPhaseCoeffSet.hh"
#include "LauRandom.hh"
#include "LauRealImagCoeffSet.hh"
#include "LauTimeDepFitModel.hh"
#include "LauVetoes.hh"
#include "LauFlavTag.hh"
#include "Lau1DHistPdf.hh"
#include "Lau1DCubicSpline.hh"
#include "Test_Dpipi_ProgOpts.hh"
int main(const int argc, const char ** argv)
{
const TestDpipi_ProgramSettings settings{argc,argv};
if ( settings.helpRequested ) {
return EXIT_SUCCESS;
}
if ( ! settings.parsedOK ) {
return EXIT_FAILURE;
}
const Bool_t fixPhiMix{ settings.fixPhiMix || settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS };
const Bool_t useSinCos{kTRUE};
Double_t nSigEvents{0};
switch (settings.dType) {
case LauTimeDepFitModel::CPEigenvalue::CPEven :
nSigEvents = 15000;
break;
case LauTimeDepFitModel::CPEigenvalue::CPOdd :
nSigEvents = 5000;
break;
case LauTimeDepFitModel::CPEigenvalue::QFS :
nSigEvents = 50000;
break;
}
LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0");
LauDaughters* daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar");
// efficiency
LauVetoes* vetoes = new LauVetoes();
//vetoes->addMassVeto( 2, 2.00776, 2.01276 );
LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes);
LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes);
//Args for flavTag: useAveDelta - kFALSE and useEtaPrime - kFALSE
LauFlavTag* flavTag = new LauFlavTag(kFALSE,kFALSE);
flavTag->setTrueTagVarName("trueTag");
+ if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS) {
+ flavTag->setDecayFlvVarName("decayFlv");
+ }
TFile* etaFile = TFile::Open("ft-eta-hist.root");
TH1* etaHist = dynamic_cast<TH1*>(etaFile->Get("ft_eta_hist"));
Lau1DHistPdf* etaHistPDF = new Lau1DHistPdf("eta",etaHist,0.0,0.5,kTRUE,kFALSE);
const Double_t meanEta { etaHistPDF->getMean() };
// if the tagging is perfect then also make it perfectly efficient, otherwise 50% efficient
const Double_t tagEffVal { (meanEta == 0.0) ? 1.0 : 0.5 };
std::pair<Double_t,Double_t> tagEff {tagEffVal, tagEffVal};
// use a null calibration for the time being, so p0 = <eta> and p1 = 1
std::pair<Double_t,Double_t> calib0 {meanEta, meanEta};
std::pair<Double_t,Double_t> calib1 {1.0, 1.0};
flavTag->addTagger("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1);
// signal dynamics
LauIsobarDynamics* sigModelB0bar = new LauIsobarDynamics(daughtersB0bar, effModelB0bar);
sigModelB0bar->setIntFileName("integ_B0bar.dat");
sigModelB0bar->addResonance("D*+_2", 2, LauAbsResonance::RelBW);
sigModelB0bar->addResonance("D*+_0", 2, LauAbsResonance::RelBW);
sigModelB0bar->addResonance("rho0(770)", 3, LauAbsResonance::RelBW);
sigModelB0bar->addResonance("f_0(980)", 3, LauAbsResonance::RelBW);
sigModelB0bar->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW);
LauIsobarDynamics* sigModelB0 = new LauIsobarDynamics(daughtersB0, effModelB0);
sigModelB0->setIntFileName("integ_B0.dat");
sigModelB0->addResonance("D*-_2", 1, LauAbsResonance::RelBW);
sigModelB0->addResonance("D*-_0", 1, LauAbsResonance::RelBW);
sigModelB0->addResonance("rho0(770)", 3, LauAbsResonance::RelBW);
sigModelB0->addResonance("f_0(980)", 3, LauAbsResonance::RelBW);
sigModelB0->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW);
// fit model
LauTimeDepFitModel* fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag);
std::vector<LauAbsCoeffSet*> coeffset;
coeffset.push_back( new LauRealImagCoeffSet("D*+_2", 1.00, 0.00, kTRUE, kTRUE) );
coeffset.push_back( new LauRealImagCoeffSet("D*+_0", 0.53*TMath::Cos( 3.00), 0.53*TMath::Sin( 3.00), kFALSE, kFALSE) );
coeffset.push_back( new LauRealImagCoeffSet("rho0(770)", 1.22*TMath::Cos( 2.25), 1.22*TMath::Sin( 2.25), kFALSE, kFALSE) );
coeffset.push_back( new LauRealImagCoeffSet("f_0(980)", 0.19*TMath::Cos(-2.48), 0.19*TMath::Sin(-2.48), kFALSE, kFALSE) );
coeffset.push_back( new LauRealImagCoeffSet("f_2(1270)", 0.75*TMath::Cos( 2.97), 0.75*TMath::Sin( 2.97), kFALSE, kFALSE) );
for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
fitModel->setAmpCoeffSet(*iter);
}
fitModel->setCPEigenvalue( settings.dType );
fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos );
// production asymmetry
fitModel->setAsymmetries(0.0,kTRUE);
// Delta t PDFs
const Double_t minDt(0.0);
const Double_t maxDt(15.0);
const Double_t minDtErr(0.0);
const Double_t maxDtErr(0.215);
const std::vector<Bool_t> scale {
settings.perEventTimeErr && kTRUE,
settings.perEventTimeErr && kTRUE,
};
const UInt_t nGauss(scale.size());
LauParameter * mean0 = new LauParameter("dt_mean_0", scale[0] ? -1.63e-3 : -1.84e-03, -0.01, 0.01, kTRUE );
LauParameter * mean1 = new LauParameter("dt_mean_1", scale[1] ? -1.63e-3 : -3.62e-03, -0.01, 0.01, kTRUE );
LauParameter * sigma0 = new LauParameter("dt_sigma_0", scale[0] ? 0.991 : 3.05e-02, 0.0, 2.0, kTRUE );
LauParameter * sigma1 = new LauParameter("dt_sigma_1", scale[1] ? 1.80 : 6.22e-02, 0.0, 2.5, kTRUE );
LauParameter * frac1 = new LauParameter("dt_frac_1", scale[0] && scale[1] ? 0.065 : 0.761, 0.0, 1.0, kTRUE);
LauParameter * tau = new LauParameter("dt_tau", 1.520, 0.5, 5.0, settings.fixLifetime);
LauParameter * freq = new LauParameter("dt_deltaM", 0.5064, 0.0, 1.0, settings.fixDeltaM);
std::vector<LauAbsRValue*> dtPars {
mean0,
mean1,
sigma0,
sigma1,
frac1,
tau,
freq
};
// Decay time acceptance histogram
TFile* dtaFile = TFile::Open("dta-hist.root");
TH1* dtaHist = dynamic_cast<TH1*>(dtaFile->Get("dta_hist"));
// Create the spline knot positions and y-values from the histogram contents
std::vector<Double_t> dtvals;
std::vector<Double_t> effvals;
dtvals.push_back(minDt);
effvals.push_back(0.0);
for ( Int_t bin{0}; bin < dtaHist->GetNbinsX(); ++bin ) {
dtvals.push_back( dtaHist->GetBinCenter(bin+1) );
effvals.push_back( dtaHist->GetBinContent(bin+1) );
}
dtvals.push_back(maxDt);
effvals.push_back(effvals.back());
// Decay time error histogram
TFile* dteFile = TFile::Open("dte-hist.root");
TH1* dteHist = dynamic_cast<TH1*>(dteFile->Get("dte_hist"));
LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "decayTime", "decayTimeErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime, settings.timeEffModel );
dtPdf->doSmearing(settings.timeResolution);
if ( settings.perEventTimeErr ) {
dtPdf->setErrorHisto( dteHist );
}
switch(settings.timeEffModel)
{
case LauDecayTimePdf::EfficiencyMethod::Spline:
{
- fitModel->setASqMaxValue(0.34);
+ fitModel->setASqMaxValue(0.06);
Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural);
dtPdf->setEffiSpline(dtEffSpline);
break;
}
case LauDecayTimePdf::EfficiencyMethod::Binned:
{
- fitModel->setASqMaxValue(0.34);
+ fitModel->setASqMaxValue(0.06);
dtPdf->setEffiHist(dtaHist);
break;
}
case LauDecayTimePdf::EfficiencyMethod::Flat:
{
- fitModel->setASqMaxValue(3.4);
+ fitModel->setASqMaxValue(4.1);
break;
}
}
fitModel->setSignalDtPdf( dtPdf );
// set the number of signal events
std::cout<<"nSigEvents = "<<nSigEvents<<std::endl;
LauParameter* nSigPar = new LauParameter("signalEvents", nSigEvents, -2.0*nSigEvents, 2.0*nSigEvents, kTRUE);
fitModel->setNSigEvents(nSigPar);
// set the number of experiments
if (settings.command == Command::Generate) {
fitModel->setNExpts(settings.nExptGen, settings.firstExptGen);
} else {
fitModel->setNExpts(settings.nExptFit, settings.firstExptFit);
}
fitModel->useAsymmFitErrors(kFALSE);
fitModel->useRandomInitFitPars(kFALSE);
fitModel->doPoissonSmearing(kFALSE);
fitModel->doEMLFit(kFALSE);
fitModel->writeLatexTable(kFALSE);
TString dTypeStr;
switch (settings.dType) {
case LauTimeDepFitModel::CPEigenvalue::CPEven :
dTypeStr = "CPEven";
break;
case LauTimeDepFitModel::CPEigenvalue::CPOdd :
dTypeStr = "CPOdd";
break;
case LauTimeDepFitModel::CPEigenvalue::QFS :
dTypeStr = "QFS";
break;
}
TString dataFile("");
TString treeName("fitTree");
TString rootFileName("");
TString tableFileName("");
TString fitToyFileName("");
TString splotFileName("");
dataFile = "TEST-Dpipi";
dataFile += "_"+dTypeStr;
switch(settings.timeEffModel)
{
case LauDecayTimePdf::EfficiencyMethod::Spline:
dataFile += "_Spline";
break;
case LauDecayTimePdf::EfficiencyMethod::Binned:
dataFile += "_Hist";
break;
case LauDecayTimePdf::EfficiencyMethod::Flat:
dataFile += "_Flat";
break;
}
if (settings.timeResolution) {
if (settings.perEventTimeErr) {
dataFile += "_DTRperevt";
} else {
dataFile += "_DTRavg";
}
} else {
dataFile += "_DTRoff";
}
dataFile += "_expts";
dataFile += settings.firstExptGen;
dataFile += "-";
dataFile += settings.firstExptGen+settings.nExptGen-1;
dataFile += ".root";
if (settings.command == Command::Generate) {
rootFileName = "dummy.root";
tableFileName = "genResults";
} else {
rootFileName = "fit"; rootFileName += settings.iFit;
rootFileName += "_Results_"; rootFileName += dTypeStr;
rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1;
rootFileName += ".root";
tableFileName = "fit"; tableFileName += settings.iFit;
tableFileName += "_Results_"; tableFileName += dTypeStr;
tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1;
fitToyFileName = "fit"; fitToyFileName += settings.iFit;
fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr;
fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1;
fitToyFileName += ".root";
splotFileName = "fit"; splotFileName += settings.iFit;
splotFileName += "_sPlot_"; splotFileName += dTypeStr;
splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1;
splotFileName += ".root";
}
// Generate toy from the fitted parameters
//fitModel->compareFitData(1, fitToyFileName);
// Write out per-event likelihoods and sWeights
//fitModel->writeSPlotData(splotFileName, "splot", kFALSE);
// Execute the generation/fit
switch (settings.command) {
case Command::Generate :
fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName );
break;
case Command::Fit :
fitModel->run( "fit", dataFile, treeName, rootFileName, tableFileName );
break;
case Command::SimFit :
fitModel->runTask( dataFile, treeName, rootFileName, tableFileName, "localhost", settings.port );
break;
}
return EXIT_SUCCESS;
}
diff --git a/inc/LauFlavTag.hh b/inc/LauFlavTag.hh
index f9e9191..902f874 100644
--- a/inc/LauFlavTag.hh
+++ b/inc/LauFlavTag.hh
@@ -1,236 +1,236 @@
/*
Copyright 2017 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauFlavTag.hh
\brief File containing declaration of LauFlavTag class.
*/
/*! \class LauFlavTag
\brief Class for defining the flavour tagging approach.
Define the flavour tagging categories and all associated parameters
to be passed to the relevant fit models.
*/
#ifndef LAU_FLAVTAG
#define LAU_FLAVTAG
// TODO - audit these includes, there seem to be a number that are not necessary
#include <vector>
#include <map>
#include <set>
#include <iostream>
#include "TString.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauParameter.hh"
#include "LauFitDataTree.hh"
#include "LauAbsFitModel.hh"
#include "LauDecayTimePdf.hh"
#include "LauAbsPdf.hh"
class LauFlavTag final {
public:
//! Constructor
/*!
\param [in] useAveDelta use average and delta variables for tagging calibration and efficiency
\param [in] useEtaPrime use eta prime rather the eta as the mistag throughout
*/
LauFlavTag(const Bool_t useAveDelta = kFALSE, const Bool_t useEtaPrime = kFALSE);
//! Initialise
// TODO is this needed? Commented for the moment (here and where called in LauTimeDepFitModel)
//void initialise();
// TODO - need to decide which functions need to be public (interface) and which should be private (implementation details)
//! Change the dilutions, delta dilutions and tagCatFrac for signal if needed
/*!
\param [in] name the name of the tagger
\param [in] tagVarName the tagging variable name of the tagger in the ntuple
\param [in] mistagVarName the associated mistag variable name of the same tagger in the ntuple
\param [in] etapdf the mistag distribution for the tagger
\param [in] tagEff tagging efficiency - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag
\param [in] calib_p0 calibration parameter p0 - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag
\param [in] calib_p1 calibration parameter p1 - (particle, antiparticle) or (average, delta) depending on useAveDelta_ flag
*/
// Need to set remember the position in the vector using a map for later reference
//void addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf,
// const Double_t tagEff_b0=1.0, const Double_t calib_p0_b0=1.0, const Double_t calib_p1_b0=1.0,
// const Double_t tagEff_b0bar=-1.0, const Double_t calib_p0_b0bar=-1.0, const Double_t calib_p1_b0bar=-1.0);
void addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf,
const std::pair <Double_t, Double_t> tagEff, const std::pair <Double_t, Double_t> calib_p0,
const std::pair <Double_t, Double_t> calib_p1);
//! Read in the input fit data variables, e.g. m13Sq and m23Sq
void cacheInputFitVars(LauFitDataTree* inputFitData);
void updateEventInfo(const ULong_t iEvt);
const std::vector<TString>& getTagVarNames() const {return tagVarName_;};
const std::vector<TString>& getMistagVarNames() const {return mistagVarName_;};
const TString& getTrueTagVarName() const {return trueTagVarName_;};
- const TString& getDecayVarName() const {return decayVarName_;};
+ const TString& getDecayFlvVarName() const {return decayFlvVarName_;};
Int_t getCurEvtTrueTagFlv() const {return curEvtTrueTagFlv_;};
Int_t getCurEvtDecayFlv() const {return curEvtDecayFlv_;};
const std::vector<Int_t>& getCurEvtTagFlv() const {return curEvtTagFlv_;};
const std::vector<Double_t>& getCurEvtMistag() const {return curEvtMistag_;};
ULong_t getNTaggers() const {return tagVarName_.size();}
//! Get map of calibration parameters for each tagging category
std::vector<LauParameter*> getCalibP0B0(){return calib_p0_B0_;};
std::vector<LauParameter*> getCalibP0B0bar(){return calib_p0_B0bar_;};
std::vector<LauParameter*> getCalibP1B0(){return calib_p1_B0_;};
std::vector<LauParameter*> getCalibP1B0bar(){return calib_p1_B0bar_;};
//! Get map of alternative calibration parameters for each tagging category
std::vector<LauParameter*> getCalibP0Ave(){return calib_p0_ave_;};
std::vector<LauParameter*> getCalibP0Delta(){return calib_p0_delta_;};
std::vector<LauParameter*> getCalibP1Ave(){return calib_p1_ave_;};
std::vector<LauParameter*> getCalibP1Delta(){return calib_p1_delta_;};
//! Get map of alternative calibration parameters for each tagging category
std::vector<LauParameter*> getTagEffAve(){return tagEff_ave_;};
std::vector<LauParameter*> getTagEffDelta(){return tagEff_delta_;};
std::vector<LauParameter*> getTagEffB0(){return tagEff_B0_;};
std::vector<LauParameter*> getTagEffB0bar(){return tagEff_B0bar_;};
const std::vector<Double_t>& getPerEvtAvgMistag() const {return perEvtAvgMistag_;};
Double_t getLittleOmega(const ULong_t position, const Int_t flag) const;
Double_t getCapitalOmega(const ULong_t position, const Int_t flag) const;
Double_t getEtaGen(const ULong_t position);
//! Return the Boolean controlling if we use the alternative tagging calibration parameters
Bool_t getUseAveDelta() const {return useAveDelta_;};
void setTrueTagVarName(TString trueTagVarName);
- void setDecayVarName(TString decayVarName);
+ void setDecayFlvVarName(TString decayFlvVarName);
//! Gaussian constraints for P0 parameters for a given tagger
/*!
param [in] name name of the tagger
param [in] constraint1 the (mean, sigma) for the particle or average parameter
param [in] constraint2 the (mean, sigma) for the antiparticle or delta parameter
*/
void addP0GaussianConstraints(const TString name, const std::pair <Double_t, Double_t> constraint1, const std::pair <Double_t, Double_t> constraint2);
//! Gaussian constraints for P1 parameters for a given tagger
/*!
param [in] name name of the tagger
param [in] constraint1 the (mean, sigma) for the particle or average parameter
param [in] constraint2 the (mean, sigma) for the antiparticle or delta parameter
*/
void addP1GaussianConstraints(const TString name, const std::pair <Double_t, Double_t> constraint1, const std::pair <Double_t, Double_t> constraint2);
//! Gaussian constraints for tagging efficiency parameters for a given tagger
/*!
param [in] name name of the tagger
param [in] constraint1 the (mean, sigma) for the particle or average parameter
param [in] constraint2 the (mean, sigma) for the antiparticle or delta parameter
*/
void addTagEffGaussianConstraints(const TString name, const std::pair <Double_t, Double_t> constraint1, const std::pair <Double_t, Double_t> constraint2);
private:
//! Map to link taggers to their vector position
std::map<TString, Int_t> taggerPosition_;
//! Flavour tagging variable name
std::vector<TString> tagVarName_;
//! Per event mistag variable name
std::vector<TString> mistagVarName_;
//! True tag variable name for normalisation decays
TString trueTagVarName_;
//! Decay flavour tag variable name for normalisation decays
- TString decayVarName_;
+ TString decayFlvVarName_;
//! Vector of flavour tags for each event
std::vector< std::vector<Int_t> > evtTagFlv_;
//! Flavour tag for current event
std::vector<Int_t> curEvtTagFlv_;
//! Vector of mistags for each event
std::vector< std::vector<Double_t> > evtMistag_;
//! Per event mistag for current event
std::vector<Double_t> curEvtMistag_;
//! Vector of true tags for each event
std::vector< Int_t > evtTrueTagFlv_;
//! Vector of decay tags for each event
std::vector< Int_t > evtDecayFlv_;
//! True tag from normalisation mode for current event
Int_t curEvtTrueTagFlv_{0};
//! True tag from normalisation mode for current event
Int_t curEvtDecayFlv_{0};
//! Per-event average mistag value (eta hat)
std::vector<Double_t> perEvtAvgMistag_;
//! Calibration parameters
std::vector<LauParameter*> calib_p0_B0_;
std::vector<LauParameter*> calib_p0_B0bar_;
std::vector<LauParameter*> calib_p1_B0_;
std::vector<LauParameter*> calib_p1_B0bar_;
//! Alternative calibration parameters
std::vector<LauParameter*> calib_p0_ave_;
std::vector<LauParameter*> calib_p0_delta_;
std::vector<LauParameter*> calib_p1_ave_;
std::vector<LauParameter*> calib_p1_delta_;
//! Flag to use alternative calibration parameters
Bool_t useAveDelta_;
//! Flag to use eta prime not eta for the mistag
Bool_t useEtaPrime_;
//! Tagging efficiency parameters
std::vector<LauParameter*> tagEff_B0_;
std::vector<LauParameter*> tagEff_B0bar_;
std::vector<LauParameter*> tagEff_ave_;
std::vector<LauParameter*> tagEff_delta_;
//! Eta PDFs
std::vector<LauAbsPdf*> etaPdfs_;
ClassDef(LauFlavTag,0) // Flavour tagging set up
};
#endif
diff --git a/inc/LauTimeDepFitModel.hh b/inc/LauTimeDepFitModel.hh
index 93f4476..2daf216 100644
--- a/inc/LauTimeDepFitModel.hh
+++ b/inc/LauTimeDepFitModel.hh
@@ -1,712 +1,715 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauTimeDepFitModel.hh
\brief File containing declaration of LauTimeDepFitModel class.
*/
/*! \class LauTimeDepFitModel
\brief Class for defining a time-dependent fit model.
LauTimeDepFitModel is a class that allows the user to define a three-body
Dalitz plot according to the isobar model, i.e. defining a set of
resonances that have complex amplitudes that can interfere with each other.
It extends the LauSimpleFitModel and LauCPFitModel models in that it allows
the fitting of time-dependent particle/antiparticle decays to
flavour-conjugate Dalitz plots, including their interference through mixing.
*/
#ifndef LAU_TIMEDEP_FIT_MODEL
#define LAU_TIMEDEP_FIT_MODEL
#include <vector>
#include <map>
#include <set>
#include <iostream>
#include "TString.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "LauAbsFitModel.hh"
#include "LauConstants.hh"
#include "LauEmbeddedData.hh"
#include "LauParameter.hh"
#include "LauFlavTag.hh"
#include "LauCategoryFlavTag.hh"
class LauAbsBkgndDPModel;
class LauAbsCoeffSet;
class LauAbsPdf;
class LauDecayTimePdf;
class LauIsobarDynamics;
class LauKinematics;
class LauScfMap;
class LauTimeDepFitModel : public LauAbsFitModel {
public:
//! Possible CP eigenvalues (the intrinsic CP of the final state particles)
enum CPEigenvalue {
CPOdd = -1, /*!< CP odd final state */
QFS = 0, /*!< Quasi Flavour Specific final state */
CPEven = 1 /*!< CP even final state */
};
//! Constructor
/*!
\param [in] modelB0bar DP model for the antiparticle
\param [in] modelB0 DP model for the particle
\param [in] flavTag flavour tagging information
*/
LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag);
//! Destructor
virtual ~LauTimeDepFitModel();
//! Set the signal event yield
/*!
\param [in] nSigEvents contains the signal yield and option to fix it
*/
virtual void setNSigEvents(LauParameter* nSigEvents);
//! Set the signal event yield and asymmetry
/*!
\param [in] nSigEvents contains the signal yield and option to fix it
\param [in] sigAsym contains the signal asymmetry and option to fix it
*/
virtual void setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym);
//! Set the number of background events
/*!
The name of the parameter must be that of the corresponding background category (so that it can be correctly assigned)
\param [in] nBkgndEvents contains the name, yield and option to fix the yield of the background
*/
virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents);
//! Set the background event yield and asymmetry
/*!
\param [in] nBkgEvents contains the background yield and option to fix it
\param [in] BkgAsym contains the background asymmetry and option to fix it
*/
virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym);
//! Set the background DP models
/*!
\param [in] bkgndClass the name of the background class
\param [in] model the DP model of the background
*/
void setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* model);
//! Switch on/off storage of amplitude info in generated ntuple
void storeGenAmpInfo(Bool_t storeInfo) { storeGenAmpInfo_ = storeInfo; }
//! Set CP eigenvalue
/*!
The CP eigenvalue can be supplied on an event-by-event
basis, e.g. if the data contains daughters that are D0
mesons that can decay to either K+ K- (CP even) or KS pi0
(CP odd).
This method allows you to set the default value that should
be used if the data does not contain this information as
well as the name of the variable in the data that will
specify this information.
If completely unspecified all events will be assumed to be
CP even.
\param defaultCPEV the default for the eigenvalue
\param evVarName the variable name in the data tree that specifies the CP eigenvalue
*/
inline void setCPEigenvalue( const CPEigenvalue defaultCPEV, const TString& cpevVarName = "" )
{
cpEigenValue_ = defaultCPEV;
cpevVarName_ = cpevVarName;
}
//! Set the DP amplitude coefficients
/*!
\param [in] coeffSet the set of coefficients
*/
void setAmpCoeffSet(LauAbsCoeffSet* coeffSet);
//! Set the decay time PDFs
/*!
\param [in] position the tagger position in the vectors
\param [in] pdf the signal decay time PDF
*/
void setSignalDtPdf(LauDecayTimePdf* pdf);
//! Set the decay time PDFs
/*!
\param [in] tagCat the tagging category for which the PDF should be used
\param [in] pdf the background decay time PDF
*/
void setBackgroundDtPdf(LauDecayTimePdf* pdf);
//! Set the signal PDF for a given variable
/*!
\param [in] tagCat the tagging category for which the PDF should be used
\param [in] pdf the PDF to be added to the signal model
*/
void setSignalPdfs(LauAbsPdf* pdf);
//! Set the background PDF
/*!
\param [in] bkgndClass the name of the background class
\param [in] pdf the PDF to be added to the background model
*/
void setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf);
void setSignalFlavTagPdfs( const Int_t tagCat, LauAbsPdf* pdf);
void setBkgdFlavTagPdfs( const TString name, LauAbsPdf* pdf);
//! Embed full simulation events for the signal, rather than generating toy from the PDFs
/*!
\param [in] tagCat the tagging category for which the file should be used
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble
\param [in] reuseEventsWithinExperiment
\param [in] useReweighting
*/
void embedSignal(const TString& fileName, const TString& treeName,
const Bool_t reuseEventsWithinEnsemble, const Bool_t reuseEventsWithinExperiment = kFALSE);
void embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE);
//! Set the value of the mixing phase
/*!
\param [in] phiMix the value of the mixing phase
\param [in] fixPhiMix whether the value should be fixed or floated
\param [in] useSinCos whether to use the sine and cosine as separate parameters or to just use the mixing phase itself
*/
void setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos = kFALSE);
//! Initialise the fit
virtual void initialise();
//! Initialise the signal DP model
virtual void initialiseDPModels();
//! Recalculate Normalization the signal DP models
virtual void recalculateNormalisation();
//! Update the coefficients
virtual void updateCoeffs();
// Toy MC generation and fitting overloaded functions
virtual Bool_t genExpt();
//! Set the maximum value of A squared to be used in the accept/reject
/*!
\param [in] value the new value
*/
inline void setASqMaxValue(const Double_t value) {aSqMaxSet_ = value;}
//! Weight events based on the DP model
/*!
\param [in] dataFileName the name of the data file
\param [in] dataTreeName the name of the data tree
*/
virtual void weightEvents( const TString& dataFileName, const TString& dataTreeName );
//! Calculate things that depend on the fit parameters after they have been updated by Minuit
virtual void propagateParUpdates();
//! Read in the input fit data variables, e.g. m13Sq and m23Sq
virtual void cacheInputFitVars();
//! Check the initial fit parameters
virtual void checkInitFitParams();
//! Get the fit results and store them
/*!
\param [in] tablePrefixName prefix for the name of the output file
*/
virtual void finaliseFitResults(const TString& tablePrefixName);
//! Save the pdf Plots for all the resonances of experiment number fitExp
/*!
TODO - not working in this model!!
\param [in] label prefix for the file name to be saved
*/
virtual void savePDFPlots(const TString& label);
//! Save the pdf Plots for the sum of ressonances correspondint to "sin" of experiment number fitExp
/*!
TODO - not working in this model!!
\param [in] label prefix for the file name to be saved
\param [in] spin spin of the wave to be saved
*/
virtual void savePDFPlotsWave(const TString& label, const Int_t& spin);
//! Print the fit fractions, total DP rate and mean efficiency
/*!
\param [out] output the stream to which to print
*/
virtual void printFitFractions(std::ostream& output);
//! Print the asymmetries
/*!
\param [out] output the stream to which to print
*/
virtual void printAsymmetries(std::ostream& output);
//! Write the fit results in latex table format
/*!
\param [in] outputFile the name of the output file
*/
virtual void writeOutTable(const TString& outputFile);
//! Store the per event likelihood values
virtual void storePerEvtLlhds();
// Methods to do with calculating the likelihood functions
// and manipulating the fitting parameters.
//! Get the total likelihood for each event
/*!
\param [in] iEvt the event number
*/
virtual Double_t getTotEvtLikelihood(const UInt_t iEvt);
//! Calculate the signal and background likelihoods for the DP for a given event
/*!
\param [in] iEvt the event number
*/
virtual void getEvtDPDtLikelihood(const UInt_t iEvt);
//! Determine the signal and background likelihood for the extra variables for a given event
/*!
\param [in] iEvt the event number
*/
virtual void getEvtExtraLikelihoods(const UInt_t iEvt);
virtual void getEvtFlavTagLikelihood(const UInt_t iEvt);
//! Get the total number of events
/*!
\return the total number of events
*/
virtual Double_t getEventSum() const;
//! Set the fit parameters for the DP model
void setSignalDPParameters();
//! Set the fit parameters for the decay time PDFs
void setDecayTimeParameters();
//! Set the fit parameters for the extra PDFs
void setExtraPdfParameters();
//! Set the initial yields
void setFitNEvents();
//! Set the calibration parameters
void setCalibParams();
//! Set the tagging efficiency parameters
void setTagEffParams();
//! Set the efficiency parameters
void setEffiParams();
//! Set the asymmetry parameters
void setAsymParams();
//! Set the tagging asymmetry parameters
void setFlavTagAsymParams();
//! Set-up other parameters that are derived from the fit results, e.g. fit fractions
void setExtraNtupleVars();
//! Set production and detections asymmetries
void setAsymmetries(const Double_t AProd, const Bool_t AProdFix);
//! Randomise the initial fit parameters
void randomiseInitFitPars();
//! Method to set up the storage for background-related quantities called by setBkgndClassNames
virtual void setupBkgndVectors();
//! Calculate the CP asymmetries
/*!
\param [in] initValues is this before or after the fit
*/
void calcAsymmetries(const Bool_t initValues = kFALSE);
//! Finalise value of mixing phase
void checkMixingPhase();
//! Return the map of signal decay time PDFs
typedef std::map< Int_t, LauDecayTimePdf*> LauTagCatDtPdfMap;
LauDecayTimePdf* getSignalDecayTimePdf(){return signalDecayTimePdf_;}
//! Return the map of background decay time PDFs
std::vector<LauDecayTimePdf*> getBackgroundDecayTimePdfs(){return backgroundDecayTimePdfs_;}
protected:
typedef std::map< Int_t, LauParameter> LauTagCatParamMap;
typedef std::map< Int_t, LauPdfList > LauTagCatPdfListMap;
typedef std::map< Int_t, LauAbsPdf* > LauTagCatPdfMap;
typedef std::map< TString, LauAbsPdf* > LauBkgdPdfMap;
typedef std::map< Int_t, Int_t > LauTagCatEventsMap;
typedef std::map< Int_t, LauEmbeddedData* > LauTagCatEmbDataMap;
//typedef std::map< Int_t, std::pair<Int_t,Double_t> > LauTaggerGenInfo;
//typedef std::map< std::pair<TString,Int_t>, LauTaggerGenInfo > LauGenInfo;
typedef std::map< TString, std::pair<Int_t,Double_t> > LauGenInfo;
typedef std::vector<LauTagCatEmbDataMap> LauTagCatEmbDataMapList;
typedef std::vector<LauAbsBkgndDPModel*> LauBkgndDPModelList;
typedef std::vector<LauPdfList> LauBkgndPdfsList;
typedef std::vector<LauAbsRValue*> LauBkgndYieldList;
typedef std::vector<Bool_t> LauBkgndReuseEventsList;
//! Determine the number of events to generate for each hypothesis
LauGenInfo eventsToGenerate();
//! Generate signal event
Bool_t generateSignalEvent();
//! Generate background event
/*!
\param [in] bgID ID number of the background class
*/
Bool_t generateBkgndEvent(UInt_t bgID);
//! Setup the required ntuple branches
void setupGenNtupleBranches();
//! Store all of the DP and decay time information
void setDPDtBranchValues();
//! Generate from the extra PDFs
/*!
\param [in] extraPdfs the list of extra PDFs
\param [in] embeddedData the embedded data sample
*/
void generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData);
//! Add sPlot branches for the extra PDFs
/*!
\param [in] extraPdfs the list of extra PDFs
\param [in] prefix the list of prefixes for the branch names
*/
void addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix);
//! Set the branches for the sPlot ntuple with extra PDFs
/*!
\param [in] extraPdfs the list of extra PDFs
\param [in] prefix the list of prefixes for the branch names
\param [in] iEvt the event number
*/
Double_t setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt);
//! Update the signal events after Minuit sets background parameters
void updateSigEvents();
//! Add branches to store experiment number and the event number within the experiment
virtual void setupSPlotNtupleBranches();
//! Returns the names of all variables in the fit
virtual LauSPlot::NameSet variableNames() const;
//! Returns the names and yields of species that are free in the fit
virtual LauSPlot::NumbMap freeSpeciesNames() const;
//! Returns the names and yields of species that are fixed in the fit
virtual LauSPlot::NumbMap fixdSpeciesNames() const;
//! Returns the species and variables for all 2D PDFs in the fit
virtual LauSPlot::TwoDMap twodimPDFs() const;
//! Check if the signal is split into well-reconstructed and mis-reconstructed types
virtual Bool_t splitSignal() const { return kFALSE; }
//! Check if the mis-reconstructed signal is to be smeared in the DP
virtual Bool_t scfDPSmear() const { return kFALSE; }
//! Add the parameters from each PDF into the fit
/*!
\param [in] theMap the container of PDFs
*/
UInt_t addParametersToFitList(LauPdfList* theList);
//! Add the parameters from each decay time PDF into the fit
/*!
\param [in] theMap the container of PDFs
*/
UInt_t addParametersToFitList(std::vector<LauDecayTimePdf*> theVector);
//! Calculate the component integrals of the interference term
void calcInterferenceTermIntegrals();
//! Calculate the total integral of the interference term
void calcInterTermNorm();
private:
//! Dalitz plot PDF for the antiparticle
LauIsobarDynamics* sigModelB0bar_;
//! Dalitz plot PDF for the particle
LauIsobarDynamics* sigModelB0_;
//! Kinematics object for antiparticle
LauKinematics* kinematicsB0bar_;
//! Kinematics object for particle
LauKinematics* kinematicsB0_;
//! The background Dalitz plot models
LauBkgndDPModelList BkgndDPModels_;
//! The background PDFs
LauBkgndPdfsList BkgndPdfs_;
//! Background boolean
Bool_t usingBkgnd_;
//! LauFlavTag object for flavour tagging
LauFlavTag* flavTag_;
//! Flavour tag for current event
std::vector<Int_t> curEvtTagFlv_;
//! Per event mistag for current event
std::vector<Double_t> curEvtMistag_;
//! Per event transformed mistag for current event
std::vector<Double_t> curEvtMistagPrime_;
- //! True tag for the current event
+ //! True tag flavour (i.e. flavour at production) for the current event (only relevant for toy generation)
Int_t curEvtTrueTagFlv_;
+ //! Flavour at decay for the current event (only relevant for QFS)
+ Int_t curEvtDecayFlv_;
+
//! Number of signal components
UInt_t nSigComp_;
//! Number of signal DP parameters
UInt_t nSigDPPar_;
//! Number of decay time PDF parameters
UInt_t nDecayTimePar_;
//! Number of extra PDF parameters
UInt_t nExtraPdfPar_;
//! Number of normalisation parameters (yields, asymmetries)
UInt_t nNormPar_;
//! Number of calibration parameters (p0, p1)
UInt_t nCalibPar_;
//! Number of tagging efficneyc parameters
UInt_t nTagEffPar_;
//! Number of efficiency parameters (p0, p1)
UInt_t nEffiPar_;
//! Number of asymmetry parameters
UInt_t nAsymPar_;
//! Number of tagging asymmetry parameters
UInt_t nTagAsymPar_;
//! The complex coefficients for antiparticle
std::vector<LauComplex> coeffsB0bar_;
//! The complex coefficients for particle
std::vector<LauComplex> coeffsB0_;
//! Magnitudes and Phases or Real and Imaginary Parts
std::vector<LauAbsCoeffSet*> coeffPars_;
//! The integrals of the efficiency corrected amplitude cross terms for each pair of amplitude components
/*!
Calculated as the sum of A* x Abar x efficiency
*/
std::vector< std::vector<LauComplex> > fifjEffSum_;
//! The normalisation for the term 2.0*Re(A*Abar*phiMix)
Double_t interTermReNorm_;
//! The normalisation for the term 2.0*Im(A*Abar*phiMix)
Double_t interTermImNorm_;
//! The antiparticle fit fractions
LauParArray fitFracB0bar_;
//! The particle fit fractions
LauParArray fitFracB0_;
//! The fit fraction asymmetries
std::vector<LauParameter> fitFracAsymm_;
//! A_CP parameter
std::vector<LauParameter> acp_;
//! The mean efficiency for the antiparticle
LauParameter meanEffB0bar_;
//! The mean efficiency for the particle
LauParameter meanEffB0_;
//! The average DP rate for the antiparticle
LauParameter DPRateB0bar_;
//! The average DP rate for the particle
LauParameter DPRateB0_;
//! Signal yields
LauParameter* signalEvents_;
//! Signal asymmetry
LauParameter* signalAsym_;
//! Background yield(s)
LauBkgndYieldList bkgndEvents_;
//! Background asymmetries(s)
LauBkgndYieldList bkgndAsym_;
//! CP eigenvalue variable name
TString cpevVarName_;
//! CP eigenvalue for current event
CPEigenvalue cpEigenValue_;
//! Vector to store event CP eigenvalues
std::vector<CPEigenvalue> evtCPEigenVals_;
//! The mass difference between the neutral mass eigenstates
LauParameter deltaM_;
//! The width difference between the neutral mass eigenstates
LauParameter deltaGamma_;
//! The lifetime
LauParameter tau_;
//! The mixing phase
LauParameter phiMix_;
//! The sine of the mixing phase
LauParameter sinPhiMix_;
//! The cosine of the mixing phase
LauParameter cosPhiMix_;
//! Flag whether to use the sine and cosine of the mixing phase or the phase itself as the fit parameters
Bool_t useSinCos_;
//! e^{-i*phiMix}
LauComplex phiMixComplex_;
//! Signal Decay time PDFs (one per tagging category)
LauDecayTimePdf* signalDecayTimePdf_;
//! Background Decay time PDFs (one per tagging category)
std::vector<LauDecayTimePdf*> backgroundDecayTimePdfs_;
//! Decay time for the current event
Double_t curEvtDecayTime_;
//! Decay time error for the current event
Double_t curEvtDecayTimeErr_;
//! PDFs for other variables
LauPdfList* sigExtraPdf_;
//! eta PDFs for each TagCat
LauPdfList sigFlavTagPdf_;
//! eta PDFs for each background
LauBkgdPdfMap bkgdFlavTagPdf_;
//! Production asymmetry between B0 and B0bar
LauParameter AProd_;
// Toy generation stuff
//! The maximum allowed number of attempts when generating an event
Int_t iterationsMax_;
//! The number of unsucessful attempts to generate an event so far
Int_t nGenLoop_;
//! The value of A squared for the current event
Double_t ASq_;
//! The maximum value of A squared that has been seen so far while generating
Double_t aSqMaxVar_;
//! The maximum allowed value of A squared
Double_t aSqMaxSet_;
//! Flag for storage of amplitude info in generated ntuple
Bool_t storeGenAmpInfo_;
//! The signal event tree for embedding fully simulated events
LauEmbeddedData* signalTree_;
//! The background event tree for embedding fully simulated events
std::vector<LauEmbeddedData*> bkgndTree_;
//! Boolean to control reuse of embedded signal events
Bool_t reuseSignal_;
//! Vector of booleans to reuse background events
LauBkgndReuseEventsList reuseBkgnd_;
// Likelihood values
//! Signal DP likelihood value
Double_t sigDPLike_;
//! Signal likelihood from extra PDFs
Double_t sigExtraLike_;
Double_t sigFlavTagLike_;
Double_t bkgdFlavTagLike_;
//! Total signal likelihood
Double_t sigTotalLike_;
//! Background DP likelihood value(s)
std::vector<Double_t> bkgndDPLike_;
//! Background likelihood value(s) from extra PDFs
std::vector<Double_t> bkgndExtraLike_;
//! Total background likelihood(s)
std::vector<Double_t> bkgndTotalLike_;
ClassDef(LauTimeDepFitModel,0) // Time-dependent neutral model
};
#endif
diff --git a/src/LauFlavTag.cc b/src/LauFlavTag.cc
index b475f53..65c1de3 100644
--- a/src/LauFlavTag.cc
+++ b/src/LauFlavTag.cc
@@ -1,467 +1,467 @@
/*
Copyright 2017 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauFlavTag.cc
\brief File containing implementation of LauFlavTag class.
*/
// TODO - audit these includes, there seem to be a number that are not necessary
#include <iostream>
#include <iomanip>
#include <fstream>
#include <map>
#include <vector>
#include "TFile.h"
#include "TMinuit.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TVirtualFitter.h"
#include "TH1D.h"
#include "LauAbsBkgndDPModel.hh"
#include "LauAbsCoeffSet.hh"
#include "LauAbsPdf.hh"
#include "LauAsymmCalc.hh"
#include "LauComplex.hh"
#include "LauConstants.hh"
#include "LauDPPartialIntegralInfo.hh"
#include "LauDaughters.hh"
#include "LauDecayTimePdf.hh"
#include "LauFitNtuple.hh"
#include "LauGenNtuple.hh"
#include "LauIsobarDynamics.hh"
#include "LauKinematics.hh"
#include "LauPrint.hh"
#include "LauRandom.hh"
#include "LauScfMap.hh"
#include "LauFlavTag.hh"
#include "Lau1DHistPdf.hh"
ClassImp(LauFlavTag)
LauFlavTag::LauFlavTag(const Bool_t useAveDelta, const Bool_t useEtaPrime) :
useAveDelta_(useAveDelta),
useEtaPrime_(useEtaPrime)
{
}
void LauFlavTag::addTagger(const TString& name, const TString& tagVarName, const TString& mistagVarName, LauAbsPdf* etapdf, const std::pair <Double_t, Double_t> tagEff, const std::pair <Double_t, Double_t> calib_p0, const std::pair <Double_t, Double_t> calib_p1)
{
// Find how many taggers have already been added
const ULong_t position { tagVarName_.size() };
// Update map to relate tagger name and position in the vectors
taggerPosition_[name]=position;
// Fill vectors
tagVarName_.push_back(tagVarName);
mistagVarName_.push_back(mistagVarName);
if (etapdf){
etaPdfs_.push_back(etapdf);
Lau1DHistPdf* etahistpdf = dynamic_cast<Lau1DHistPdf*>(etapdf);
if (etahistpdf){
perEvtAvgMistag_.push_back(etahistpdf->getMean());
} else {
std::cerr << "WARNING in LauFlavTag::addTagger : Couldn't determine average eta value from PDF. Setting it to 0.4." << std::endl;
perEvtAvgMistag_.push_back(0.4);
}
} else {
std::cerr << "ERROR in LauFlavTag::addTagger : Eta PDF pointer is NULL" << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
//Use particle/antiparticle variables
if (!useAveDelta_){
TString tagEff_b0Name("tagEff_b0_"+name);
TString tagEff_b0barName("tagEff_b0bar_"+name);
TString calib_p0_b0Name("calib_p0_b0_"+name);
TString calib_p0_b0barName("calib_p0_b0bar_"+name);
TString calib_p1_b0Name("calib_p1_b0_"+name);
TString calib_p1_b0barName("calib_p1_b0bar_"+name);
LauParameter* tageffb0 = new LauParameter(tagEff_b0Name,tagEff.first,0.0,1.0,kTRUE);
tagEff_B0_.push_back(tageffb0);
tagEff_B0_[position]->initValue(tagEff.first); tagEff_B0_[position]->genValue(tagEff.first);
tagEff_B0_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp0b0 = new LauParameter(calib_p0_b0Name,calib_p0.first,-10.0,10.0,kTRUE);
calib_p0_B0_.push_back(calibp0b0);
calib_p0_B0_[position]->initValue(calib_p0.first); calib_p0_B0_[position]->genValue(calib_p0.first);
calib_p0_B0_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp1b0 = new LauParameter(calib_p1_b0Name,calib_p1.first,0.0,1.5,kTRUE);
calib_p1_B0_.push_back(calibp1b0);
calib_p1_B0_[position]->initValue(calib_p1.first); calib_p1_B0_[position]->genValue(calib_p1.first);
calib_p1_B0_[position]->fixed(kTRUE); //Update once full code in place
if (tagEff.second==-1.0 && calib_p0.second==-1.0 && calib_p1.second==-1.0){
tagEff_B0bar_.push_back(tagEff_B0_[position]->createClone(tagEff_b0barName));
calib_p0_B0bar_.push_back(calib_p0_B0_[position]->createClone(calib_p0_b0barName));
calib_p1_B0bar_.push_back(calib_p1_B0_[position]->createClone(calib_p1_b0barName));
} else {
LauParameter* tageffb0bar = new LauParameter(tagEff_b0barName,tagEff.second,0.0,1.0,kTRUE);
tagEff_B0bar_.push_back(tageffb0bar);
tagEff_B0bar_[position]->initValue(tagEff.second); tagEff_B0bar_[position]->genValue(tagEff.second);
tagEff_B0bar_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp0b0bar = new LauParameter(calib_p0_b0barName,calib_p0.second,-10.0,10.0,kTRUE);
calib_p0_B0bar_.push_back(calibp0b0bar);
calib_p0_B0bar_[position]->initValue(calib_p0.second); calib_p0_B0bar_[position]->genValue(calib_p0.second);
calib_p0_B0bar_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp1b0bar = new LauParameter(calib_p1_b0barName,calib_p1.second,0.0,1.5,kTRUE);
calib_p1_B0bar_.push_back(calibp1b0bar);
calib_p1_B0bar_[position]->initValue(calib_p1.second); calib_p1_B0bar_[position]->genValue(calib_p1.second);
calib_p1_B0bar_[position]->fixed(kTRUE); //Update once full code in place
}
} else { //Use average and delta variables
TString tagEff_aveName("tagEff_ave_"+name);
TString tagEff_deltaName("tagEff_delta_"+name);
TString calib_p0_aveName("calib_p0_ave_"+name);
TString calib_p0_deltaName("calib_p0_delta_"+name);
TString calib_p1_aveName("calib_p1_ave_"+name);
TString calib_p1_deltaName("calib_p1_delta_"+name);
LauParameter* tageffave = new LauParameter(tagEff_aveName,tagEff.first,0.0,1.0,kTRUE);
tagEff_ave_.push_back(tageffave);
tagEff_ave_[position]->initValue(tagEff.first); tagEff_ave_[position]->genValue(tagEff.first);
tagEff_ave_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp0ave = new LauParameter(calib_p0_aveName,calib_p0.first,-10.0,10.0,kTRUE);
calib_p0_ave_.push_back(calibp0ave);
calib_p0_ave_[position]->initValue(calib_p0.first); calib_p0_ave_[position]->genValue(calib_p0.first);
calib_p0_ave_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp1ave = new LauParameter(calib_p1_aveName,calib_p1.first,0.0,1.5,kTRUE);
calib_p1_ave_.push_back(calibp1ave);
calib_p1_ave_[position]->initValue(calib_p1.first); calib_p1_ave_[position]->genValue(calib_p1.first);
calib_p1_ave_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* tageffdelta = new LauParameter(tagEff_deltaName,tagEff.second,-1.0,1.0,kTRUE);
tagEff_delta_.push_back(tageffdelta);
tagEff_delta_[position]->initValue(tagEff.second); tagEff_delta_[position]->genValue(tagEff.second);
tagEff_delta_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp0delta = new LauParameter(calib_p0_deltaName,calib_p0.second,-10.0,10.0,kTRUE);
calib_p0_delta_.push_back(calibp0delta);
calib_p0_delta_[position]->initValue(calib_p0.second); calib_p0_delta_[position]->genValue(calib_p0.second);
calib_p0_delta_[position]->fixed(kTRUE); //Update once full code in place
LauParameter* calibp1delta = new LauParameter(calib_p1_deltaName,calib_p1.second,-10.0,10.0,kTRUE);
calib_p1_delta_.push_back(calibp1delta);
calib_p1_delta_[position]->initValue(calib_p1.second); calib_p1_delta_[position]->genValue(calib_p1.second);
calib_p1_delta_[position]->fixed(kTRUE); //Update once full code in place
}
std::cout<<"INFO in LauFlavTag::addTagger : Added tagger with name "<< name << std::endl;
}
void LauFlavTag::cacheInputFitVars(LauFitDataTree* inputFitData)
{
evtTagFlv_.clear();
evtMistag_.clear();
evtTrueTagFlv_.clear();
evtDecayFlv_.clear();
// Loop over the taggers to check the branches
for (ULong_t i=0; i < tagVarName_.size(); ++i){
if ( ! inputFitData->haveBranch( tagVarName_[i] ) ) {
std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << tagVarName_[i] << "\"." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( ! inputFitData->haveBranch( mistagVarName_[i] ) ) {
std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << mistagVarName_[i] << "\"." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
if ( ! inputFitData->haveBranch( trueTagVarName_ ) ) {
std::cerr << "ERROR in LauFlavTag::cacheInputFitVars : Input data does not contain branch \"" << trueTagVarName_ << "\"." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
const ULong_t nEvents { inputFitData->nEvents() };
evtTagFlv_.reserve( nEvents );
evtMistag_.reserve( nEvents );
evtTrueTagFlv_.reserve( nEvents );
evtDecayFlv_.reserve( nEvents );
LauFitData::const_iterator fitdata_iter;
for (ULong_t iEvt = 0; iEvt < nEvents; iEvt++) {
const LauFitData& dataValues = inputFitData->getData(iEvt);
//Clear vectors
curEvtTagFlv_.clear();
curEvtMistag_.clear();
// For untagged events see if we have a truth tag for normalisation modes
- curEvtTrueTagFlv_ = static_cast<Int_t>( dataValues.at( trueTagVarName_ ) );
+ curEvtTrueTagFlv_ = ( trueTagVarName_ != "" ) ? static_cast<Int_t>( dataValues.at( trueTagVarName_ ) ) : 0;
if ( curEvtTrueTagFlv_ > 1 ) {
std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid true tagging output " << curEvtTrueTagFlv_ << " for event " << iEvt << ", setting it to +1" << std::endl;
curEvtTrueTagFlv_ = +1;
} else if ( curEvtTrueTagFlv_ < -1 ){
std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid true tagging output " << curEvtTrueTagFlv_ << " for event " << iEvt << ", setting it to -1" << std::endl;
curEvtTrueTagFlv_ = -1;
}
evtTrueTagFlv_.push_back(curEvtTrueTagFlv_);
// Flavour at decay
- curEvtDecayFlv_ = static_cast<Int_t>( dataValues.at( decayVarName_ ) );
+ curEvtDecayFlv_ = ( decayFlvVarName_ != "" ) ? static_cast<Int_t>( dataValues.at( decayFlvVarName_ ) ) : 0;
if ( curEvtDecayFlv_ > 1 ) {
std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid decay tagging output " << curEvtDecayFlv_ << " for event " << iEvt << ", setting it to +1" << std::endl;
curEvtDecayFlv_ = +1;
} else if ( curEvtDecayFlv_ < -1 ){
std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid decay tagging output " << curEvtDecayFlv_ << " for event " << iEvt << ", setting it to -1" << std::endl;
curEvtDecayFlv_ = -1;
}
evtDecayFlv_.push_back(curEvtDecayFlv_);
for (ULong_t i=0; i < tagVarName_.size(); ++i){
curEvtTagFlv_.push_back( static_cast<Int_t>( dataValues.at( tagVarName_[i] ) ) );
if ( curEvtTagFlv_[i] > 1 ) {
std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_[i] << " for event " << iEvt << ", setting it to +1" << std::endl;
curEvtTagFlv_[i] = +1;
} else if ( curEvtTagFlv_[i] < -1 ) {
std::cerr << "WARNING in LauFlavTag::cacheInputFitVars : Invalid tagging output " << curEvtTagFlv_[i] << " for event " << iEvt << ", setting it to -1" << std::endl;
curEvtTagFlv_[i] = -1;
}
curEvtMistag_.push_back( static_cast<Double_t>( dataValues.at( mistagVarName_[i] ) ) );
// Calibrated mistag > 0.5 is just a tag flip - handled automatically in getCapitalOmega function
if (curEvtMistag_[i] > 0.5){
std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<<curEvtMistag_[i]<<" is larger than 0.5, setting to 0.5"<<std::endl;
curEvtMistag_[i] = 0.5;
}
if (curEvtMistag_[i] < 0.0){
std::cerr<<"WARNING in LauFlavTag::cacheInputFitVars : Mistag value "<<curEvtMistag_[i]<<" is less than 0.0, setting to 0.0"<<std::endl;
curEvtMistag_[i] = 0.0;
}
}
evtTagFlv_.push_back(curEvtTagFlv_);
evtMistag_.push_back(curEvtMistag_);
}
}
void LauFlavTag::updateEventInfo(const ULong_t iEvt)
{
//Clear vector
curEvtTagFlv_.clear();
curEvtMistag_.clear();
//Assign current event variables
curEvtTagFlv_ = evtTagFlv_[iEvt];
curEvtMistag_ = evtMistag_[iEvt];
curEvtTrueTagFlv_ = evtTrueTagFlv_[iEvt];
curEvtDecayFlv_ = evtDecayFlv_[iEvt];
}
Double_t LauFlavTag::getLittleOmega(const ULong_t position, const Int_t flag) const
{
if (TMath::Abs(flag) != 1){
std::cerr << "ERROR in LauFlavTag::getLittleOmega : Invalid flag, you must request either omega (+1) or omega bar (-1) to be returned" << std::endl;
return 0.0;
}
Double_t calibp0(0.), calibp1(0.), calibp0bar(0.), calibp1bar(0.);
//If we are floating average omega and delta omega we need to use those parameters instead
if (useAveDelta_){
calibp0 = calib_p0_ave_[position]->unblindValue() + 0.5*calib_p0_delta_[position]->unblindValue();
calibp0bar = calib_p0_ave_[position]->unblindValue() - 0.5*calib_p0_delta_[position]->unblindValue();
calibp1 = calib_p1_ave_[position]->unblindValue() + 0.5*calib_p1_delta_[position]->unblindValue();
calibp1bar = calib_p1_ave_[position]->unblindValue() - 0.5*calib_p1_delta_[position]->unblindValue();
} else {
calibp0 = calib_p0_B0_[position]->unblindValue();
calibp0bar = calib_p0_B0bar_[position]->unblindValue();
calibp1 = calib_p1_B0_[position]->unblindValue();
calibp1bar = calib_p1_B0bar_[position]->unblindValue();
}
if (flag == 1){
return calibp0 + calibp1 * (curEvtMistag_[position] - perEvtAvgMistag_[position]);
}
else{
return calibp0bar + calibp1bar * (curEvtMistag_[position] - perEvtAvgMistag_[position]);
}
std::cerr << "ERROR in LauFlavTag::getLittleOmega : Current event flavour tag not defined" << std::endl;
return 0.0;
}
Double_t LauFlavTag::getCapitalOmega(const ULong_t position, const Int_t flag) const
{
if (TMath::Abs(flag) != 1){
std::cerr << "ERROR in LauFlavTag::getCapitalOmega : Invalid flag, you must request either Omega (+1) or Omega bar (-1) to be returned" << std::endl;
return 0.0;
}
//Delta functions to control which terms contribute
Int_t deltap1(0), deltam1(0), delta0(0);
if (curEvtTagFlv_[position] == -1){
deltam1 = 1;
} else if(curEvtTagFlv_[position] == 1){
deltap1 = 1;
} else{
delta0 = 1;
}
//Efficiency
Double_t eff=0.0;
if (useAveDelta_){
if(flag==1){
eff = tagEff_ave_[position]->unblindValue() + 0.5*tagEff_delta_[position]->unblindValue();
} else {
eff = tagEff_ave_[position]->unblindValue() - 0.5*tagEff_delta_[position]->unblindValue();
}
}else{
if(flag==1){
eff = tagEff_B0_[position]->unblindValue();
}else{
eff = tagEff_B0bar_[position]->unblindValue();
}
}
//Little omega
Double_t omega = this->getLittleOmega(position, flag);
Double_t omegaPrime(0.);
//Transform to omega prime - TODO isn't this the inverse, getLittleOmega is actually giving us omegaPrime and on the next line we convert back to omega?
if (useEtaPrime_){
omegaPrime = (1/(1+TMath::Exp(-1.0*omega)));
}else{
omegaPrime = omega;
}
//little omega must be between 0 and 1. Force this for now, if the fits keep getting stuck can look more closely at it.
if (omegaPrime < 0.0){
std::cerr << "WARNING in LauFlavTag::getCapitalOmega the value of little omega is less than 0, shifting to 0" << std::endl;
omegaPrime = 0.0;
}
if (omegaPrime > 1.0){
std::cerr << "WARNING in LauFlavTag::getCapitalOmega the value of little omega is greater than 1, shifting to 1" << std::endl;
omegaPrime = 1.0;
}
//eta PDF value
std::vector<Double_t> abs;
abs.push_back(curEvtMistag_[position]);
etaPdfs_[position]->calcLikelihoodInfo(abs);
const Double_t h { etaPdfs_[position]->getLikelihood() };
const Double_t u { 2.0 }; // the PDF value for a uniform PDF between 0.0 and 0.5
//Put it together
if (flag == 1){ //Particle
return (deltap1*eff*(1-omegaPrime) + deltam1*eff*omegaPrime)*h + delta0*(1-eff)*u;
} else {
return (deltam1*eff*(1-omegaPrime) + deltap1*eff*omegaPrime)*h + delta0*(1-eff)*u;
}
}
Double_t LauFlavTag::getEtaGen(const ULong_t position)
{
//Clear mistag vector for a new event
if(position==0){
curEvtMistag_.clear();
}
LauFitData data { etaPdfs_[position]->generate(nullptr) }; //TODO Add DP dependence?
Double_t etagen { data.at(etaPdfs_[position]->varName()) };
if (etagen > 0.5){etagen = 0.5;}
if (etagen < 0.0){etagen = 0.0;}
curEvtMistag_.push_back(etagen);
return etagen;
}
void LauFlavTag::setTrueTagVarName(TString trueTagVarName){
trueTagVarName_ = std::move(trueTagVarName);
}
-void LauFlavTag::setDecayVarName(TString decayVarName){
- decayVarName_ = std::move(decayVarName);
+void LauFlavTag::setDecayFlvVarName(TString decayFlvVarName){
+ decayFlvVarName_ = std::move(decayFlvVarName);
}
void LauFlavTag::addP0GaussianConstraints(TString name, std::pair <Double_t, Double_t> constraint1, std::pair <Double_t, Double_t> constraint2){
//Does key exist?
if (taggerPosition_.count(name)==0){
std::cerr << "ERROR in LauFlavTag::addP0GaussianConstraints : Tagger name not recognised please check your options" << std::endl;
std::cerr << "ERROR in LauFlavTag::addP0GaussianConstraints : Constraints have not been applied" << std::endl;
return;
}
//Find position in the vector from the tagger name
Double_t pos = taggerPosition_.at(name);
if (!useAveDelta_){
calib_p0_B0_[pos]->addGaussianConstraint(constraint1.first,constraint1.second);
calib_p0_B0bar_[pos]->addGaussianConstraint(constraint2.first,constraint2.second);
}else{
calib_p0_ave_[pos]->addGaussianConstraint(constraint1.first,constraint1.second);
calib_p0_delta_[pos]->addGaussianConstraint(constraint2.first,constraint2.second);
}
std::cout << "INFO in LauFlavTag::addP0GaussianConstraints : Added Gaussian constraints for the P0 calibration parameters of tagger " << name << std::endl;
}
void LauFlavTag::addP1GaussianConstraints(TString name, std::pair <Double_t, Double_t> constraint1, std::pair <Double_t, Double_t> constraint2){
//Does key exist?
if (taggerPosition_.count(name)==0){
std::cerr << "ERROR in LauFlavTag::addP1GaussianConstraints : Tagger name not recognised please check your options" << std::endl;
std::cerr << "ERROR in LauFlavTag::addP1GaussianConstraints : Constraints have not been applied" << std::endl;
return;
}
//Find position in the vector from the tagger name
Double_t pos = taggerPosition_.at(name);
if (!useAveDelta_){
calib_p1_B0_[pos]->addGaussianConstraint(constraint1.first,constraint1.second);
calib_p1_B0bar_[pos]->addGaussianConstraint(constraint2.first,constraint2.second);
}else{
calib_p1_ave_[pos]->addGaussianConstraint(constraint1.first,constraint1.second);
calib_p1_delta_[pos]->addGaussianConstraint(constraint2.first,constraint2.second);
}
std::cout << "INFO in LauFlavTag::addP1GaussianConstraints : Added Gaussian constraints for the P1 calibration parameters of tagger " << name << std::endl;
}
void LauFlavTag::addTagEffGaussianConstraints(TString name, std::pair <Double_t, Double_t> constraint1, std::pair <Double_t, Double_t> constraint2){
//Does key exist?
if (taggerPosition_.count(name)==0){
std::cerr << "ERROR in LauFlavTag::addTagEffGaussianConstraints : Tagger name not recognised please check your options" << std::endl;
std::cerr << "ERROR in LauFlavTag::addTagEffGaussianConstraints : Constraints have not been applied" << std::endl;
return;
}
//Find position in the vector from the tagger name
Double_t pos = taggerPosition_.at(name);
if (!useAveDelta_){
tagEff_B0_[pos]->addGaussianConstraint(constraint1.first,constraint1.second);
tagEff_B0bar_[pos]->addGaussianConstraint(constraint2.first,constraint2.second);
}else{
tagEff_ave_[pos]->addGaussianConstraint(constraint1.first,constraint1.second);
tagEff_delta_[pos]->addGaussianConstraint(constraint2.first,constraint2.second);
}
std::cout << "INFO in LauFlavTag::addTagEffGaussianConstraints : Added Gaussian constraints for the tagging efficiency parameters of tagger " << name << std::endl;
}
diff --git a/src/LauTimeDepFitModel.cc b/src/LauTimeDepFitModel.cc
index 4bdd10d..deefbea 100644
--- a/src/LauTimeDepFitModel.cc
+++ b/src/LauTimeDepFitModel.cc
@@ -1,2899 +1,2948 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauTimeDepFitModel.cc
\brief File containing implementation of LauTimeDepFitModel class.
*/
#include <iostream>
#include <iomanip>
#include <fstream>
#include <map>
#include <vector>
#include "TFile.h"
#include "TMinuit.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TVirtualFitter.h"
#include "LauAbsBkgndDPModel.hh"
#include "LauAbsCoeffSet.hh"
#include "LauAbsPdf.hh"
#include "LauAsymmCalc.hh"
#include "LauComplex.hh"
#include "LauConstants.hh"
#include "LauDPPartialIntegralInfo.hh"
#include "LauDaughters.hh"
#include "LauDecayTimePdf.hh"
#include "LauFitNtuple.hh"
#include "LauGenNtuple.hh"
#include "LauIsobarDynamics.hh"
#include "LauKinematics.hh"
#include "LauPrint.hh"
#include "LauRandom.hh"
#include "LauScfMap.hh"
#include "LauTimeDepFitModel.hh"
#include "LauFlavTag.hh"
ClassImp(LauTimeDepFitModel)
LauTimeDepFitModel::LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag) : LauAbsFitModel(),
sigModelB0bar_(modelB0bar),
sigModelB0_(modelB0),
kinematicsB0bar_(modelB0bar ? modelB0bar->getKinematics() : 0),
kinematicsB0_(modelB0 ? modelB0->getKinematics() : 0),
usingBkgnd_(kFALSE),
flavTag_(flavTag),
+ curEvtTrueTagFlv_(0),
+ curEvtDecayFlv_(0),
nSigComp_(0),
nSigDPPar_(0),
nDecayTimePar_(0),
nExtraPdfPar_(0),
nNormPar_(0),
nCalibPar_(0),
nTagEffPar_(0),
nEffiPar_(0),
nAsymPar_(0),
coeffsB0bar_(0),
coeffsB0_(0),
coeffPars_(0),
fitFracB0bar_(0),
fitFracB0_(0),
fitFracAsymm_(0),
acp_(0),
meanEffB0bar_("meanEffB0bar",0.0,0.0,1.0),
meanEffB0_("meanEffB0",0.0,0.0,1.0),
DPRateB0bar_("DPRateB0bar",0.0,0.0,100.0),
DPRateB0_("DPRateB0",0.0,0.0,100.0),
signalEvents_(0),
signalAsym_(0),
cpevVarName_(""),
cpEigenValue_(CPEven),
evtCPEigenVals_(0),
deltaM_("deltaM",0.0),
deltaGamma_("deltaGamma",0.0),
tau_("tau",LauConstants::tauB0),
phiMix_("phiMix", 2.0*LauConstants::beta, -LauConstants::threePi, LauConstants::threePi, kFALSE),
sinPhiMix_("sinPhiMix", TMath::Sin(2.0*LauConstants::beta), -1.0, 1.0, kFALSE),
cosPhiMix_("cosPhiMix", TMath::Cos(2.0*LauConstants::beta), -1.0, 1.0, kFALSE),
useSinCos_(kFALSE),
phiMixComplex_(TMath::Cos(-2.0*LauConstants::beta),TMath::Sin(-2.0*LauConstants::beta)),
signalDecayTimePdf_(),
backgroundDecayTimePdfs_(),
curEvtDecayTime_(0.0),
curEvtDecayTimeErr_(0.0),
sigExtraPdf_(),
sigFlavTagPdf_(),
bkgdFlavTagPdf_(),
AProd_("AProd",0.0,-1.0,1.0,kTRUE),
iterationsMax_(100000000),
nGenLoop_(0),
ASq_(0.0),
aSqMaxVar_(0.0),
aSqMaxSet_(1.25),
storeGenAmpInfo_(kFALSE),
signalTree_(),
reuseSignal_(kFALSE),
sigDPLike_(0.0),
sigExtraLike_(0.0),
sigFlavTagLike_(0.0),
bkgdFlavTagLike_(0.0),
sigTotalLike_(0.0)
{
// Set up ftag here?
// Make sure that the integration scheme will be symmetrised
sigModelB0bar_->forceSymmetriseIntegration(kTRUE);
sigModelB0_->forceSymmetriseIntegration(kTRUE);
}
LauTimeDepFitModel::~LauTimeDepFitModel()
{
for (LauPdfList::iterator pdf_iter = sigExtraPdf_->begin(); pdf_iter != sigExtraPdf_->end(); ++pdf_iter) {
delete *(pdf_iter);
}
for (std::vector<LauEmbeddedData*>::iterator iter = bkgndTree_.begin(); iter != bkgndTree_.end(); ++iter){
delete *(iter);
}
}
void LauTimeDepFitModel::setupBkgndVectors()
{
UInt_t nBkgnds = this->nBkgndClasses();
BkgndDPModels_.resize( nBkgnds );
BkgndPdfs_.resize( nBkgnds );
bkgndEvents_.resize( nBkgnds );
bkgndAsym_.resize( nBkgnds );
bkgndTree_.resize( nBkgnds );
reuseBkgnd_.resize( nBkgnds );
bkgndDPLike_.resize( nBkgnds );
bkgndExtraLike_.resize( nBkgnds );
bkgndTotalLike_.resize( nBkgnds );
}
void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents)
{
if ( nSigEvents == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( signalEvents_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
return;
}
if ( signalAsym_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
return;
}
signalEvents_ = nSigEvents;
signalEvents_->name("signalEvents");
Double_t value = nSigEvents->value();
signalEvents_->range(-2.0*(TMath::Abs(value)+1.0),2.0*(TMath::Abs(value)+1.0));
signalAsym_ = new LauParameter("signalAsym",0.0,-1.0,1.0,kTRUE);
}
void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym)
{
if ( nSigEvents == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The event LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( sigAsym == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The asym LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( signalEvents_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
return;
}
if ( signalAsym_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
return;
}
signalEvents_ = nSigEvents;
signalEvents_->name("signalEvents");
Double_t value = nSigEvents->value();
signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
signalAsym_ = sigAsym;
signalAsym_->name("signalAsym");
signalAsym_->range(-1.0,1.0);
}
void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents)
{
if ( nBkgndEvents == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBgkndEvents : The background yield LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
if ( bkgndEvents_[bkgndID] != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
return;
}
if ( bkgndAsym_[bkgndID] != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl;
return;
}
nBkgndEvents->name( nBkgndEvents->name()+"Events" );
if ( nBkgndEvents->isLValue() ) {
Double_t value = nBkgndEvents->value();
LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
}
bkgndEvents_[bkgndID] = nBkgndEvents;
bkgndAsym_[bkgndID] = new LauParameter(nBkgndEvents->name()+"Asym",0.0,-1.0,1.0,kTRUE);
}
void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym)
{
if ( nBkgndEvents == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( bkgndAsym == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background asym LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
if ( bkgndEvents_[bkgndID] != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
return;
}
if ( bkgndAsym_[bkgndID] != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl;
return;
}
bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
if ( nBkgndEvents->isLValue() ) {
Double_t value = nBkgndEvents->value();
LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
}
bkgndEvents_[bkgndID] = nBkgndEvents;
bkgndAsym_[bkgndID]->name( nBkgndEvents->name()+"Asym" );
if ( bkgndAsym->isLValue() ) {
LauParameter* asym = dynamic_cast<LauParameter*>( bkgndAsym );
asym->range(-1.0, 1.0);
}
bkgndAsym_[bkgndID] = bkgndAsym;
}
void LauTimeDepFitModel::setSignalDtPdf(LauDecayTimePdf* pdf)
{
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : The PDF pointer is null, not adding it."<<std::endl;
return;
}
signalDecayTimePdf_ = pdf;
}
void LauTimeDepFitModel::setBackgroundDtPdf(LauDecayTimePdf* pdf)
{
// TODO If these are all histograms shouldn't need to add much more code in other functions
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setBackgroundDtPdf : The PDF pointer is null, not adding it."<<std::endl;
return;
}
backgroundDecayTimePdfs_.push_back(pdf);
usingBkgnd_ = kTRUE;
}
void LauTimeDepFitModel::setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* model)
{
if (model==0) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModels : the model pointer is null." << std::endl;
return;
}
// check that this background name is valid
if ( ! this->validBkgndClass( bkgndClass) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModel : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
return;
}
UInt_t bkgndID = this->bkgndClassID( bkgndClass );
BkgndDPModels_[bkgndID] = model;
usingBkgnd_ = kTRUE;
}
void LauTimeDepFitModel::setSignalPdfs(LauAbsPdf* pdf)
{
// These "extra variables" are assumed to be purely kinematical, like mES and DeltaE
//or making use of Rest of Event information, and therefore independent of whether
//the parent is a B0 or a B0bar. If this assupmtion doesn't hold, do modify this part!
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setSignalPdfs : The PDF pointer is null."<<std::endl;
return;
}
sigExtraPdf_->push_back(pdf);
}
void LauTimeDepFitModel::setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf)
{
if (pdf==0) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndPdf : PDF pointer is null." << std::endl;
return;
}
// check that this background name is valid
if ( ! this->validBkgndClass( bkgndClass ) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
return;
}
UInt_t bkgndID = this->bkgndClassID( bkgndClass );
BkgndPdfs_[bkgndID].push_back(pdf);
usingBkgnd_ = kTRUE;
}
void LauTimeDepFitModel::setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos)
{
phiMix_.value(phiMix); phiMix_.initValue(phiMix); phiMix_.genValue(phiMix); phiMix_.fixed(fixPhiMix);
const Double_t sinPhiMix = TMath::Sin(phiMix);
sinPhiMix_.value(sinPhiMix); sinPhiMix_.initValue(sinPhiMix); sinPhiMix_.genValue(sinPhiMix); sinPhiMix_.fixed(fixPhiMix);
const Double_t cosPhiMix = TMath::Cos(phiMix);
cosPhiMix_.value(cosPhiMix); cosPhiMix_.initValue(cosPhiMix); cosPhiMix_.genValue(cosPhiMix); cosPhiMix_.fixed(fixPhiMix);
useSinCos_ = useSinCos;
phiMixComplex_.setRealPart(cosPhiMix);
phiMixComplex_.setImagPart(-1.0*sinPhiMix);
}
void LauTimeDepFitModel::initialise()
{
// From the initial parameter values calculate the coefficients
// so they can be passed to the signal model
this->updateCoeffs();
// Initialisation
if (this->useDP() == kTRUE) {
this->initialiseDPModels();
}
//Flavour tagging
//flavTag_->initialise();
if (!this->useDP() && sigExtraPdf_->empty()) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Signal model doesn't exist for any variable."<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (this->useDP() == kTRUE) {
// Check that we have all the Dalitz-plot models
if ((sigModelB0bar_ == 0) || (sigModelB0_ == 0)) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialise : the pointer to one (particle or anti-particle) of the signal DP models is null."<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
// Next check that, if a given component is being used we've got the
// right number of PDFs for all the variables involved
// TODO - should probably check variable names and so on as well
//UInt_t nsigpdfvars(0);
//for ( LauPdfList::const_iterator pdf_iter = sigExtraPdf_.begin(); pdf_iter != sigExtraPdf_.end(); ++pdf_iter ) {
// std::vector<TString> varNames = (*pdf_iter)->varNames();
// for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
// if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
// ++nsigpdfvars;
// }
// }
//}
//if (usingBkgnd_) {
// for (LauBkgndPdfsList::const_iterator bgclass_iter = BkgndPdfsB0_.begin(); bgclass_iter != BkgndPdfsB0_.end(); ++bgclass_iter) {
// UInt_t nbkgndpdfvars(0);
// const LauPdfList& pdfList = (*bgclass_iter);
// for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) {
// std::vector<TString> varNames = (*pdf_iter)->varNames();
// for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
// if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
// ++nbkgndpdfvars;
// }
// }
// }
// if (nbkgndpdfvars != nsigpdfvars) {
// std::cerr << "ERROR in LauTimeDepFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl;
// gSystem->Exit(EXIT_FAILURE);
// }
// }
//}
// Clear the vectors of parameter information so we can start from scratch
this->clearFitParVectors();
// Set the fit parameters for signal and background models
this->setSignalDPParameters();
// Set the fit parameters for the decay time models
this->setDecayTimeParameters();
// Set the fit parameters for the extra PDFs
this->setExtraPdfParameters();
// Set the initial bg and signal events
this->setFitNEvents();
// Handle flavour-tagging calibration parameters
this->setCalibParams();
// Add tagging efficiency parameters
this->setTagEffParams();
// Add the efficiency parameters
this->setEffiParams();
//Asymmetry terms AProd and in setAsymmetries()?
//this->setAsymParams();
// Check that we have the expected number of fit variables
const LauParameterPList& fitVars = this->fitPars();
if (fitVars.size() != (nSigDPPar_ + nDecayTimePar_ + nExtraPdfPar_ + nNormPar_ + nCalibPar_ + nTagEffPar_ + nEffiPar_)) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Number of fit parameters not of expected size."<<std::endl;
std::cout<< fitVars.size() << " != " << nSigDPPar_ << nDecayTimePar_ << nExtraPdfPar_ << nNormPar_ << nCalibPar_ << nTagEffPar_ << nEffiPar_ << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
this->setExtraNtupleVars();
}
void LauTimeDepFitModel::recalculateNormalisation()
{
sigModelB0bar_->recalculateNormalisation();
sigModelB0_->recalculateNormalisation();
sigModelB0bar_->modifyDataTree();
sigModelB0_->modifyDataTree();
this->calcInterferenceTermIntegrals();
}
void LauTimeDepFitModel::initialiseDPModels()
{
if (sigModelB0bar_ == 0) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0bar signal DP model doesn't exist"<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (sigModelB0_ == 0) {
std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0 signal DP model doesn't exist"<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// Need to check that the number of components we have and that the dynamics has matches up
const UInt_t nAmpB0bar = sigModelB0bar_->getnTotAmp();
const UInt_t nAmpB0 = sigModelB0_->getnTotAmp();
if ( nAmpB0bar != nAmpB0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Unequal number of signal DP components in the particle and anti-particle models: " << nAmpB0bar << " != " << nAmpB0 << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( nAmpB0bar != nSigComp_ ) {
std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Number of signal DP components in the model (" << nAmpB0bar << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
std::cout<<"INFO in LauTimeDepFitModel::initialiseDPModels : Initialising signal DP model"<<std::endl;
sigModelB0bar_->initialise(coeffsB0bar_);
sigModelB0_->initialise(coeffsB0_);
fifjEffSum_.clear();
fifjEffSum_.resize(nSigComp_);
for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
fifjEffSum_[iAmp].resize(nSigComp_);
}
// calculate the integrals of the A*Abar terms
this->calcInterferenceTermIntegrals();
this->calcInterTermNorm();
// Add backgrounds
if (usingBkgnd_ == kTRUE) {
for (LauBkgndDPModelList::iterator iter = BkgndDPModels_.begin(); iter != BkgndDPModels_.end(); ++iter) {
(*iter)->initialise();
}
}
}
void LauTimeDepFitModel::calcInterferenceTermIntegrals()
{
const std::vector<LauDPPartialIntegralInfo*>& integralInfoListB0bar = sigModelB0bar_->getIntegralInfos();
const std::vector<LauDPPartialIntegralInfo*>& integralInfoListB0 = sigModelB0_->getIntegralInfos();
// TODO should check (first time) that they match in terms of number of entries in the vectors and that each entry has the same number of points, ranges, weights etc.
LauComplex A, Abar, fifjEffSumTerm;
for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) {
fifjEffSum_[iAmp][jAmp].zero();
}
}
const UInt_t nIntegralRegions = integralInfoListB0bar.size();
for ( UInt_t iRegion(0); iRegion < nIntegralRegions; ++iRegion ) {
const LauDPPartialIntegralInfo* integralInfoB0bar = integralInfoListB0bar[iRegion];
const LauDPPartialIntegralInfo* integralInfoB0 = integralInfoListB0[iRegion];
const UInt_t nm13Points = integralInfoB0bar->getnm13Points();
const UInt_t nm23Points = integralInfoB0bar->getnm23Points();
for (UInt_t m13 = 0; m13 < nm13Points; ++m13) {
for (UInt_t m23 = 0; m23 < nm23Points; ++m23) {
const Double_t weight = integralInfoB0bar->getWeight(m13,m23);
const Double_t eff = integralInfoB0bar->getEfficiency(m13,m23);
const Double_t effWeight = eff*weight;
for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
A = integralInfoB0->getAmplitude(m13, m23, iAmp);
for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) {
Abar = integralInfoB0bar->getAmplitude(m13, m23, jAmp);
fifjEffSumTerm = Abar*A.conj();
fifjEffSumTerm.rescale(effWeight);
fifjEffSum_[iAmp][jAmp] += fifjEffSumTerm;
}
}
}
}
}
}
void LauTimeDepFitModel::calcInterTermNorm()
{
const std::vector<Double_t>& fNormB0bar = sigModelB0bar_->getFNorm();
const std::vector<Double_t>& fNormB0 = sigModelB0_->getFNorm();
LauComplex norm;
for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) {
LauComplex coeffTerm = coeffsB0bar_[jAmp]*coeffsB0_[iAmp].conj();
coeffTerm *= fifjEffSum_[iAmp][jAmp];
coeffTerm.rescale(fNormB0bar[jAmp] * fNormB0[iAmp]);
norm += coeffTerm;
}
}
norm *= phiMixComplex_;
interTermReNorm_ = 2.0*norm.re();
interTermImNorm_ = 2.0*norm.im();
}
void LauTimeDepFitModel::setAmpCoeffSet(LauAbsCoeffSet* coeffSet)
{
// Is there a component called compName in the signal models?
TString compName = coeffSet->name();
TString conjName = sigModelB0bar_->getConjResName(compName);
const LauDaughters* daughtersB0bar = sigModelB0bar_->getDaughters();
const LauDaughters* daughtersB0 = sigModelB0_->getDaughters();
const Bool_t conjugate = daughtersB0bar->isConjugate( daughtersB0 );
if ( ! sigModelB0bar_->hasResonance(compName) ) {
if ( ! sigModelB0bar_->hasResonance(conjName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
return;
}
std::cerr<<"WARNING in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\" but does contain the conjugate \""<<conjName<<"\", resetting name to use the conjugate."<<std::endl;
TString tmp = compName;
compName = conjName;
conjName = tmp;
coeffSet->name( compName );
}
if ( conjugate ) {
if ( ! sigModelB0_->hasResonance(conjName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<conjName<<"\"."<<std::endl;
return;
}
} else {
if ( ! sigModelB0_->hasResonance(compName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
return;
}
}
// Do we already have it in our list of names?
for (std::vector<LauAbsCoeffSet*>::const_iterator iter=coeffPars_.begin(); iter!=coeffPars_.end(); ++iter) {
if ((*iter)->name() == compName) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : Have already set coefficients for \""<<compName<<"\"."<<std::endl;
return;
}
}
coeffSet->index(nSigComp_);
coeffPars_.push_back(coeffSet);
TString parName = coeffSet->baseName(); parName += "FitFracAsym";
fitFracAsymm_.push_back(LauParameter(parName, 0.0, -1.0, 1.0));
acp_.push_back(coeffSet->acp());
++nSigComp_;
std::cout<<"INFO in LauTimeDepFitModel::setAmpCoeffSet : Added coefficients for component \""<<compName<<"\" to the fit model."<<std::endl;
}
void LauTimeDepFitModel::calcAsymmetries(Bool_t initValues)
{
// Calculate the CP asymmetries
// Also calculate the fit fraction asymmetries
for (UInt_t i = 0; i < nSigComp_; i++) {
acp_[i] = coeffPars_[i]->acp();
LauAsymmCalc asymmCalc(fitFracB0bar_[i][i].value(), fitFracB0_[i][i].value());
Double_t asym = asymmCalc.getAsymmetry();
fitFracAsymm_[i].value(asym);
if (initValues) {
fitFracAsymm_[i].genValue(asym);
fitFracAsymm_[i].initValue(asym);
}
}
}
void LauTimeDepFitModel::setSignalDPParameters()
{
// Set the fit parameters for the signal model.
nSigDPPar_ = 0;
if ( ! this->useDP() ) {
return;
}
std::cout << "INFO in LauTimeDepFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl;
// Place isobar coefficient parameters in vector of fit variables
LauParameterPList& fitVars = this->fitPars();
for (UInt_t i = 0; i < nSigComp_; ++i) {
LauParameterPList pars = coeffPars_[i]->getParameters();
for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) {
if ( !(*iter)->clone() ) {
fitVars.push_back(*iter);
++nSigDPPar_;
}
}
}
// Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector
// Need to make sure that they are unique because some might appear in both DP models
LauParameterPSet& resVars = this->resPars();
resVars.clear();
LauParameterPList& sigDPParsB0bar = sigModelB0bar_->getFloatingParameters();
LauParameterPList& sigDPParsB0 = sigModelB0_->getFloatingParameters();
for ( LauParameterPList::iterator iter = sigDPParsB0bar.begin(); iter != sigDPParsB0bar.end(); ++iter ) {
if ( resVars.insert(*iter).second ) {
fitVars.push_back(*iter);
++nSigDPPar_;
}
}
for ( LauParameterPList::iterator iter = sigDPParsB0.begin(); iter != sigDPParsB0.end(); ++iter ) {
if ( resVars.insert(*iter).second ) {
fitVars.push_back(*iter);
++nSigDPPar_;
}
}
}
UInt_t LauTimeDepFitModel::addParametersToFitList(std::vector<LauDecayTimePdf*> theVector)
{
UInt_t counter(0);
LauParameterPList& fitVars = this->fitPars();
// loop through the map
for (std::vector<LauDecayTimePdf*>::iterator iter = theVector.begin(); iter != theVector.end(); ++iter) {
// grab the pdf and then its parameters
LauDecayTimePdf* thePdf = *iter; // The first one is the tagging category
LauAbsRValuePList& rvalues = thePdf->getParameters();
// loop through the parameters
for (LauAbsRValuePList::iterator pars_iter = rvalues.begin(); pars_iter != rvalues.end(); ++pars_iter) {
LauParameterPList params = (*pars_iter)->getPars();
for (LauParameterPList::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter) {
// for each "original" parameter add it to the list of fit parameters and increment the counter
if ( !(*params_iter)->clone() && ( !(*params_iter)->fixed() ||
(this->twoStageFit() && (*params_iter)->secondStage()) ) ) {
fitVars.push_back(*params_iter);
++counter;
}
}
}
}
return counter;
}
UInt_t LauTimeDepFitModel::addParametersToFitList(LauPdfList* theList)
{
UInt_t counter(0);
counter += this->addFitParameters(*(theList));
return counter;
}
void LauTimeDepFitModel::setDecayTimeParameters()
{
nDecayTimePar_ = 0;
std::cout << "INFO in LauTimeDepFitModel::setDecayTimeParameters : Setting the initial fit parameters of the DecayTime Pdfs." << std::endl;
LauParameterPList& fitVars = this->fitPars();
// Loop over the Dt PDFs
LauAbsRValuePList& rvalues = signalDecayTimePdf_->getParameters();
// loop through the parameters
for (LauAbsRValuePList::iterator pars_iter = rvalues.begin(); pars_iter != rvalues.end(); ++pars_iter) {
LauParameterPList params = (*pars_iter)->getPars();
for (LauParameterPList::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter) {
// for each "original" parameter add it to the list of fit parameters and increment the counter
if ( !(*params_iter)->clone() && ( !(*params_iter)->fixed() ||
(this->twoStageFit() && (*params_iter)->secondStage()) ) ) {
fitVars.push_back(*params_iter);
++nDecayTimePar_;
}
}
}
if (usingBkgnd_){
nDecayTimePar_ += this->addParametersToFitList(backgroundDecayTimePdfs_);
}
if (useSinCos_) {
fitVars.push_back(&sinPhiMix_);
fitVars.push_back(&cosPhiMix_);
nDecayTimePar_ += 2;
} else {
fitVars.push_back(&phiMix_);
++nDecayTimePar_;
}
}
void LauTimeDepFitModel::setExtraPdfParameters()
{
// Include the parameters of the PDF for each tagging category in the fit
// NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE)
// With the new "cloned parameter" scheme only "original" parameters are passed to the fit.
// Their clones are updated automatically when the originals are updated.
nExtraPdfPar_ = 0;
std::cout << "INFO in LauTimeDepFitModel::setExtraPdfParameters : Setting the initial fit parameters of the extra Pdfs." << std::endl;
if (sigExtraPdf_){
nExtraPdfPar_ += this->addFitParameters((*sigExtraPdf_));
}
if (usingBkgnd_ == kTRUE) {
for (LauBkgndPdfsList::iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) {
nExtraPdfPar_ += this->addFitParameters(*iter);
}
}
}
void LauTimeDepFitModel::setFitNEvents()
{
nNormPar_ = 0;
std::cout << "INFO in LauTimeDepFitModel::setFitNEvents : Setting the initial fit parameters of the signal and ackground yields." << std::endl;
// Initialise the total number of events to be the sum of all the hypotheses
Double_t nTotEvts = signalEvents_->value();
this->eventsPerExpt(TMath::FloorNint(nTotEvts));
LauParameterPList& fitVars = this->fitPars();
// if doing an extended ML fit add the signal fraction into the fit parameters
if (this->doEMLFit()) {
std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for signal and background components..."<<std::endl;
fitVars.push_back(signalEvents_);
++nNormPar_;
} else {
std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..."<<std::endl;
}
// if not using the DP in the model we need an explicit signal asymmetry parameter
if (this->useDP() == kFALSE) {
fitVars.push_back(signalAsym_);
++nNormPar_;
}
// TODO arguably should delegate this
//LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac();
// tagging-category fractions for signal events
//for (LauTagCatParamMap::iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) {
// if (iter == signalTagCatFrac.begin()) {
// continue;
// }
// LauParameter* par = &((*iter).second);
// fitVars.push_back(par);
// ++nNormPar_;
//}
// Backgrounds
if (usingBkgnd_ == kTRUE) {
for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
std::vector<LauParameter*> parameters = (*iter)->getPars();
for ( LauParameter* parameter : parameters ) {
if(!parameter->clone()) {
fitVars.push_back(parameter);
++nNormPar_;
}
}
}
for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) {
std::vector<LauParameter*> parameters = (*iter)->getPars();
for ( LauParameter* parameter : parameters ) {
if(!parameter->clone()) {
fitVars.push_back(parameter);
++nNormPar_;
}
}
}
}
}
void LauTimeDepFitModel::setAsymParams()
{
nAsymPar_ = 0;
LauParameterPList& fitVars = this->fitPars();
if (!AProd_.fixed()){
fitVars.push_back(&AProd_);
nAsymPar_+=1;
}
}
void LauTimeDepFitModel::setTagEffParams()
{
nTagEffPar_ = 0;
Bool_t useAltPars = flavTag_->getUseAveDelta();
std::cout << "INFO in LauTimeDepFitModel::setTagEffParams : Setting the initial fit parameters for flavour tagging efficiencies." << std::endl;
if (useAltPars){
std::vector<LauParameter*> tageff_ave = flavTag_->getTagEffAve();
std::vector<LauParameter*> tageff_delta = flavTag_->getTagEffDelta();
LauParameterPList& fitVars = this->fitPars();
for(std::vector<LauParameter*>::iterator iter = tageff_ave.begin(); iter != tageff_ave.end(); ++iter){
LauParameter* eff = *iter;
if (eff->fixed()){continue;}
fitVars.push_back(eff);
++nTagEffPar_;
}
for(std::vector<LauParameter*>::iterator iter = tageff_delta.begin(); iter != tageff_delta.end(); ++iter){
LauParameter* eff = *iter;
if (eff->fixed()){continue;}
fitVars.push_back(eff);
++nTagEffPar_;
}
} else {
std::vector<LauParameter*> tageff_b0 = flavTag_->getTagEffB0();
std::vector<LauParameter*> tageff_b0bar = flavTag_->getTagEffB0bar();
LauParameterPList& fitVars = this->fitPars();
for(std::vector<LauParameter*>::iterator iter = tageff_b0.begin(); iter != tageff_b0.end(); ++iter){
LauParameter* eff = *iter;
if (eff->fixed()){continue;}
fitVars.push_back(eff);
++nTagEffPar_;
}
for(std::vector<LauParameter*>::iterator iter = tageff_b0bar.begin(); iter != tageff_b0bar.end(); ++iter){
LauParameter* eff = *iter;
if (eff->fixed()){continue;}
fitVars.push_back(eff);
++nTagEffPar_;
}
}
}
void LauTimeDepFitModel::setCalibParams()
{
Bool_t useAltPars = flavTag_->getUseAveDelta();
std::cout << "INFO in LauTimeDepFitModel::setCalibParams : Setting the initial fit parameters of the flavour tagging calibration parameters." << std::endl;
if (useAltPars){
std::vector<LauParameter*> p0pars_ave = flavTag_->getCalibP0Ave();
std::vector<LauParameter*> p0pars_delta = flavTag_->getCalibP0Delta();
std::vector<LauParameter*> p1pars_ave = flavTag_->getCalibP1Ave();
std::vector<LauParameter*> p1pars_delta = flavTag_->getCalibP1Delta();
LauParameterPList& fitVars = this->fitPars();
for(std::vector<LauParameter*>::iterator iter = p0pars_ave.begin(); iter != p0pars_ave.end(); ++iter){
LauParameter* p0 = *iter;
if (p0->fixed()){continue;}
fitVars.push_back(p0);
++nCalibPar_;
}
for(std::vector<LauParameter*>::iterator iter = p0pars_delta.begin(); iter != p0pars_delta.end(); ++iter){
LauParameter* p0 = *iter;
if (p0->fixed()){continue;}
fitVars.push_back(p0);
++nCalibPar_;
}
for(std::vector<LauParameter*>::iterator iter = p1pars_ave.begin(); iter != p1pars_ave.end(); ++iter){
LauParameter* p1 = *iter;
if (p1->fixed()){continue;}
fitVars.push_back(p1);
++nCalibPar_;
}
for(std::vector<LauParameter*>::iterator iter = p1pars_delta.begin(); iter != p1pars_delta.end(); ++iter){
LauParameter* p1 = *iter;
if (p1->fixed()){continue;}
fitVars.push_back(p1);
++nCalibPar_;
}
} else {
std::vector<LauParameter*> p0pars_b0 = flavTag_->getCalibP0B0();
std::vector<LauParameter*> p0pars_b0bar = flavTag_->getCalibP0B0bar();
std::vector<LauParameter*> p1pars_b0 = flavTag_->getCalibP1B0();
std::vector<LauParameter*> p1pars_b0bar = flavTag_->getCalibP1B0bar();
LauParameterPList& fitVars = this->fitPars();
for(std::vector<LauParameter*>::iterator iter = p0pars_b0.begin(); iter != p0pars_b0.end(); ++iter){
LauParameter* p0 = *iter;
if (p0->fixed()){continue;}
fitVars.push_back(p0);
++nCalibPar_;
}
for(std::vector<LauParameter*>::iterator iter = p0pars_b0bar.begin(); iter != p0pars_b0bar.end(); ++iter){
LauParameter* p0 = *iter;
if (p0->fixed()){continue;}
fitVars.push_back(p0);
++nCalibPar_;
}
for(std::vector<LauParameter*>::iterator iter = p1pars_b0.begin(); iter != p1pars_b0.end(); ++iter){
LauParameter* p1 = *iter;
if (p1->fixed()){continue;}
fitVars.push_back(p1);
++nCalibPar_;
}
for(std::vector<LauParameter*>::iterator iter = p1pars_b0bar.begin(); iter != p1pars_b0bar.end(); ++iter){
LauParameter* p1 = *iter;
if (p1->fixed()){continue;}
fitVars.push_back(p1);
++nCalibPar_;
}
}
}
void LauTimeDepFitModel::setEffiParams()
{
nEffiPar_ = 0;
LauParameterPList& fitVars = this->fitPars();
std::vector<LauParameter*>& effiPars = signalDecayTimePdf_->getEffiPars();
for(std::vector<LauParameter*>::iterator iter = effiPars.begin(); iter != effiPars.end(); ++iter){
LauParameter* par = *iter;
if (par->fixed()){continue;}
fitVars.push_back(par);
++nEffiPar_;
}
}
void LauTimeDepFitModel::setExtraNtupleVars()
{
// Set-up other parameters derived from the fit results, e.g. fit fractions.
if (this->useDP() != kTRUE) {
return;
}
// First clear the vectors so we start from scratch
this->clearExtraVarVectors();
LauParameterList& extraVars = this->extraPars();
// Add the B0 and B0bar fit fractions for each signal component
fitFracB0bar_ = sigModelB0bar_->getFitFractions();
if (fitFracB0bar_.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar_[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j = i; j < nSigComp_; j++) {
TString name = fitFracB0bar_[i][j].name();
name.Insert( name.Index("FitFrac"), "B0bar" );
fitFracB0bar_[i][j].name(name);
extraVars.push_back(fitFracB0bar_[i][j]);
}
}
fitFracB0_ = sigModelB0_->getFitFractions();
if (fitFracB0_.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0_[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j = i; j < nSigComp_; j++) {
TString name = fitFracB0_[i][j].name();
name.Insert( name.Index("FitFrac"), "B0" );
fitFracB0_[i][j].name(name);
extraVars.push_back(fitFracB0_[i][j]);
}
}
// Calculate the ACPs and FitFrac asymmetries
this->calcAsymmetries(kTRUE);
// Add the Fit Fraction asymmetry for each signal component
for (UInt_t i = 0; i < nSigComp_; i++) {
extraVars.push_back(fitFracAsymm_[i]);
}
// Add the calculated CP asymmetry for each signal component
for (UInt_t i = 0; i < nSigComp_; i++) {
extraVars.push_back(acp_[i]);
}
// Now add in the DP efficiency values
Double_t initMeanEffB0bar = sigModelB0bar_->getMeanEff().initValue();
meanEffB0bar_.value(initMeanEffB0bar); meanEffB0bar_.initValue(initMeanEffB0bar); meanEffB0bar_.genValue(initMeanEffB0bar);
extraVars.push_back(meanEffB0bar_);
Double_t initMeanEffB0 = sigModelB0_->getMeanEff().initValue();
meanEffB0_.value(initMeanEffB0); meanEffB0_.initValue(initMeanEffB0); meanEffB0_.genValue(initMeanEffB0);
extraVars.push_back(meanEffB0_);
// Also add in the DP rates
Double_t initDPRateB0bar = sigModelB0bar_->getDPRate().initValue();
DPRateB0bar_.value(initDPRateB0bar); DPRateB0bar_.initValue(initDPRateB0bar); DPRateB0bar_.genValue(initDPRateB0bar);
extraVars.push_back(DPRateB0bar_);
Double_t initDPRateB0 = sigModelB0_->getDPRate().initValue();
DPRateB0_.value(initDPRateB0); DPRateB0_.initValue(initDPRateB0); DPRateB0_.genValue(initDPRateB0);
extraVars.push_back(DPRateB0_);
}
void LauTimeDepFitModel::setAsymmetries(const Double_t AProd, const Bool_t AProdFix){
AProd_.value(AProd);
AProd_.fixed(AProdFix);
}
void LauTimeDepFitModel::finaliseFitResults(const TString& tablePrefixName)
{
// Retrieve parameters from the fit results for calculations and toy generation
// and eventually store these in output root ntuples/text files
// Now take the fit parameters and update them as necessary
// i.e. to make mag > 0.0, phase in the right range.
// This function will also calculate any other values, such as the
// fit fractions, using any errors provided by fitParErrors as appropriate.
// Also obtain the pull values: (measured - generated)/(average error)
if (this->useDP() == kTRUE) {
for (UInt_t i = 0; i < nSigComp_; ++i) {
// Check whether we have "a > 0.0", and phases in the right range
coeffPars_[i]->finaliseValues();
}
}
// update the pulls on the event fractions and asymmetries
if (this->doEMLFit()) {
signalEvents_->updatePull();
}
if (this->useDP() == kFALSE) {
signalAsym_->updatePull();
}
// Finalise the pulls on the decay time parameters
signalDecayTimePdf_->updatePulls();
// and for backgrounds if required
if (usingBkgnd_){
for (std::vector<LauDecayTimePdf*>::iterator iter = backgroundDecayTimePdfs_.begin(); iter != backgroundDecayTimePdfs_.end(); ++iter) {
LauDecayTimePdf* pdf = *iter;
pdf->updatePulls();
}
}
if (useSinCos_) {
cosPhiMix_.updatePull();
sinPhiMix_.updatePull();
} else {
this->checkMixingPhase();
}
if (usingBkgnd_ == kTRUE) {
for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
std::vector<LauParameter*> parameters = (*iter)->getPars();
for ( LauParameter* parameter : parameters ) {
parameter->updatePull();
}
}
for (LauBkgndYieldList::iterator iter = bkgndAsym_.begin(); iter != bkgndAsym_.end(); ++iter) {
std::vector<LauParameter*> parameters = (*iter)->getPars();
for ( LauParameter* parameter : parameters ) {
parameter->updatePull();
}
}
}
// Update the pulls on all the extra PDFs' parameters
if (sigExtraPdf_){
this->updateFitParameters(*(sigExtraPdf_));
}
if (usingBkgnd_ == kTRUE) {
for (LauBkgndPdfsList::iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) {
this->updateFitParameters(*iter);
}
}
// Fill the fit results to the ntuple
// update the coefficients and then calculate the fit fractions and ACP's
if (this->useDP() == kTRUE) {
this->updateCoeffs();
sigModelB0bar_->updateCoeffs(coeffsB0bar_); sigModelB0bar_->calcExtraInfo();
sigModelB0_->updateCoeffs(coeffsB0_); sigModelB0_->calcExtraInfo();
LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions();
if (fitFracB0bar.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0bar.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0bar[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
LauParArray fitFracB0 = sigModelB0_->getFitFractions();
if (fitFracB0.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
fitFracB0bar_[i][j].value(fitFracB0bar[i][j].value());
fitFracB0_[i][j].value(fitFracB0[i][j].value());
}
}
meanEffB0bar_.value(sigModelB0bar_->getMeanEff().value());
meanEffB0_.value(sigModelB0_->getMeanEff().value());
DPRateB0bar_.value(sigModelB0bar_->getDPRate().value());
DPRateB0_.value(sigModelB0_->getDPRate().value());
this->calcAsymmetries();
// Then store the final fit parameters, and any extra parameters for
// the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate)
this->clearExtraVarVectors();
LauParameterList& extraVars = this->extraPars();
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
extraVars.push_back(fitFracB0bar_[i][j]);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
extraVars.push_back(fitFracB0_[i][j]);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
extraVars.push_back(fitFracAsymm_[i]);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
extraVars.push_back(acp_[i]);
}
extraVars.push_back(meanEffB0bar_);
extraVars.push_back(meanEffB0_);
extraVars.push_back(DPRateB0bar_);
extraVars.push_back(DPRateB0_);
this->printFitFractions(std::cout);
this->printAsymmetries(std::cout);
}
const LauParameterPList& fitVars = this->fitPars();
const LauParameterList& extraVars = this->extraPars();
LauFitNtuple* ntuple = this->fitNtuple();
ntuple->storeParsAndErrors(fitVars, extraVars);
// find out the correlation matrix for the parameters
ntuple->storeCorrMatrix(this->iExpt(), this->fitStatus(), this->covarianceMatrix());
// Fill the data into ntuple
ntuple->updateFitNtuple();
// Print out the partial fit fractions, phases and the
// averaged efficiency, reweighted by the dynamics (and anything else)
if (this->writeLatexTable()) {
TString sigOutFileName(tablePrefixName);
sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex";
this->writeOutTable(sigOutFileName);
}
}
void LauTimeDepFitModel::printFitFractions(std::ostream& output)
{
// Print out Fit Fractions, total DP rate and mean efficiency
// First for the B0bar events
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
output<<"B0bar FitFraction for component "<<i<<" ("<<compName<<") = "<<fitFracB0bar_[i][i]<<std::endl;
}
output<<"B0bar overall DP rate (integral of matrix element squared) = "<<DPRateB0bar_<<std::endl;
output<<"B0bar average efficiency weighted by whole DP dynamics = "<<meanEffB0bar_<<std::endl;
// Then for the B0 sample
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
const TString conjName(sigModelB0bar_->getConjResName(compName));
output<<"B0 FitFraction for component "<<i<<" ("<<conjName<<") = "<<fitFracB0_[i][i]<<std::endl;
}
output<<"B0 overall DP rate (integral of matrix element squared) = "<<DPRateB0_<<std::endl;
output<<"B0 average efficiency weighted by whole DP dynamics = "<<meanEffB0_<<std::endl;
}
void LauTimeDepFitModel::printAsymmetries(std::ostream& output)
{
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
output<<"Fit Fraction asymmetry for component "<<i<<" ("<<compName<<") = "<<fitFracAsymm_[i]<<std::endl;
}
for (UInt_t i = 0; i < nSigComp_; i++) {
const TString compName(coeffPars_[i]->name());
output<<"ACP for component "<<i<<" ("<<compName<<") = "<<acp_[i].value()<<" +- "<<acp_[i].error()<<std::endl;
}
}
void LauTimeDepFitModel::writeOutTable(const TString& outputFile)
{
// Write out the results of the fit to a tex-readable table
std::ofstream fout(outputFile);
LauPrint print;
std::cout<<"INFO in LauTimeDepFitModel::writeOutTable : Writing out results of the fit to the tex file "<<outputFile<<std::endl;
if (this->useDP() == kTRUE) {
// print the fit coefficients in one table
coeffPars_.front()->printTableHeading(fout);
for (UInt_t i = 0; i < nSigComp_; i++) {
coeffPars_[i]->printTableRow(fout);
}
fout<<"\\hline"<<std::endl;
fout<<"\\end{tabular}"<<std::endl<<std::endl;
// print the fit fractions in another
fout<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
fout<<"\\hline"<<std::endl;
fout<<"Component & \\Bzb\\ Fit Fraction & \\Bz\\ Fit Fraction & Fit Fraction Asymmetry & $A_{\\CP}$ \\\\"<<std::endl;
fout<<"\\hline"<<std::endl;
Double_t fitFracSumB0bar(0.0);
Double_t fitFracSumB0(0.0);
for (UInt_t i = 0; i < nSigComp_; i++) {
TString resName = coeffPars_[i]->name();
resName = resName.ReplaceAll("_", "\\_");
fout<<resName<<" & $";
Double_t fitFracB0bar = fitFracB0bar_[i][i].value();
fitFracSumB0bar += fitFracB0bar;
print.printFormat(fout, fitFracB0bar);
fout << "$ & $" << std::endl;
Double_t fitFracB0 = fitFracB0_[i][i].value();
fitFracSumB0 += fitFracB0;
print.printFormat(fout, fitFracB0);
fout << "$ & $" << std::endl;
Double_t fitFracAsymm = fitFracAsymm_[i].value();
print.printFormat(fout, fitFracAsymm);
fout << "$ & $" << std::endl;
Double_t acp = acp_[i].value();
Double_t acpErr = acp_[i].error();
print.printFormat(fout, acp);
fout<<" \\pm ";
print.printFormat(fout, acpErr);
fout<<"$ \\\\"<<std::endl;
}
fout<<"\\hline"<<std::endl;
// Also print out sum of fit fractions
fout << "Fit Fraction Sum = & $";
print.printFormat(fout, fitFracSumB0bar);
fout << "$ & $";
print.printFormat(fout, fitFracSumB0);
fout << "$ & & \\\\" << std::endl;
fout << "\\hline \n\\hline" << std::endl;
fout << "DP rate = & $";
print.printFormat(fout, DPRateB0bar_.value());
fout << "$ & $";
print.printFormat(fout, DPRateB0_.value());
fout << "$ & & \\\\" << std::endl;
fout << "$< \\varepsilon > =$ & $";
print.printFormat(fout, meanEffB0bar_.value());
fout << "$ & $";
print.printFormat(fout, meanEffB0_.value());
fout << "$ & & \\\\" << std::endl;
if (useSinCos_) {
fout << "$\\sinPhiMix =$ & $";
print.printFormat(fout, sinPhiMix_.value());
fout << " \\pm ";
print.printFormat(fout, sinPhiMix_.error());
fout << "$ & & & & & & & \\\\" << std::endl;
fout << "$\\cosPhiMix =$ & $";
print.printFormat(fout, cosPhiMix_.value());
fout << " \\pm ";
print.printFormat(fout, cosPhiMix_.error());
fout << "$ & & & & & & & \\\\" << std::endl;
} else {
fout << "$\\phiMix =$ & $";
print.printFormat(fout, phiMix_.value());
fout << " \\pm ";
print.printFormat(fout, phiMix_.error());
fout << "$ & & & & & & & \\\\" << std::endl;
}
fout << "\\hline \n\\end{tabular}" << std::endl;
}
if (!sigExtraPdf_->empty()) {
fout<<"\\begin{tabular}{|l|c|}"<<std::endl;
fout<<"\\hline"<<std::endl;
fout<<"\\Extra Signal PDFs' Parameters: & \\\\"<<std::endl;
this->printFitParameters(*(sigExtraPdf_), fout);
if (usingBkgnd_ == kTRUE && !BkgndPdfs_.empty()) {
fout << "\\hline" << std::endl;
fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl;
for (LauBkgndPdfsList::const_iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) {
this->printFitParameters(*iter, fout);
}
}
fout<<"\\hline \n\\end{tabular}"<<std::endl<<std::endl;
}
}
void LauTimeDepFitModel::checkInitFitParams()
{
// Update the number of signal events to be total-sum(background events)
this->updateSigEvents();
// Check whether we want to have randomised initial fit parameters for the signal model
if (this->useRandomInitFitPars() == kTRUE) {
this->randomiseInitFitPars();
}
}
void LauTimeDepFitModel::randomiseInitFitPars()
{
// Only randomise those parameters that are not fixed!
std::cout<<"INFO in LauTimeDepFitModel::randomiseInitFitPars : Randomising the initial values of the coefficients of the DP components (and phiMix)..."<<std::endl;
for (UInt_t i = 0; i < nSigComp_; i++) {
coeffPars_[i]->randomiseInitValues();
}
phiMix_.randomiseValue(-LauConstants::pi, LauConstants::pi);
if (useSinCos_) {
sinPhiMix_.initValue(TMath::Sin(phiMix_.initValue()));
cosPhiMix_.initValue(TMath::Cos(phiMix_.initValue()));
}
}
LauTimeDepFitModel::LauGenInfo LauTimeDepFitModel::eventsToGenerate()
{
// Determine the number of events to generate for each hypothesis
// If we're smearing then smear each one individually
// NB this individual smearing has to be done individually per tagging category as well
LauGenInfo nEvtsGen;
// Signal
// If we're including the DP and decay time we can't decide on the tag
// yet, since it depends on the whole DP+dt PDF, however, if
// we're not then we need to decide.
Double_t evtWeight(1.0);
Double_t nEvts = signalEvents_->genValue();
if ( nEvts < 0.0 ) {
evtWeight = -1.0;
nEvts = TMath::Abs( nEvts );
}
Double_t sigAsym(0.0);
if (this->useDP() == kFALSE) {
sigAsym = signalAsym_->genValue();
//TODO fill in here if we care
} else {
Double_t rateB0bar = sigModelB0bar_->getDPRate().value();
Double_t rateB0 = sigModelB0_->getDPRate().value();
if ( rateB0bar+rateB0 > 1e-30) {
sigAsym = (rateB0bar-rateB0)/(rateB0bar+rateB0);
}
//for (LauTagCatParamMap::const_iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) {
// const LauParameter& par = iter->second;
// Double_t eventsbyTagCat = par.value() * nEvts;
// if (this->doPoissonSmearing()) {
// eventsbyTagCat = LauRandom::randomFun()->Poisson(eventsbyTagCat);
// }
// eventsB0[iter->first] = std::make_pair( TMath::Nint(eventsbyTagCat), evtWeight );
//}
//nEvtsGen[std::make_pair("signal",0)] = eventsB0; // generate signal event, decide tag later.
nEvtsGen["signal"] = std::make_pair( nEvts, evtWeight );
}
std::cout<<"INFO in LauTimeDepFitModel::eventsToGenerate : Generating toy MC with:"<<std::endl;
std::cout<<" : Signal asymmetry = "<<sigAsym<<" and number of signal events = "<<signalEvents_->genValue()<<std::endl;
return nEvtsGen;
}
Bool_t LauTimeDepFitModel::genExpt()
{
// Routine to generate toy Monte Carlo events according to the various models we have defined.
// Determine the number of events to generate for each hypothesis
LauGenInfo nEvts = this->eventsToGenerate();
Bool_t genOK(kTRUE);
Int_t evtNum(0);
const UInt_t nBkgnds = this->nBkgndClasses();
std::vector<TString> bkgndClassNames(nBkgnds);
std::vector<TString> bkgndClassNamesGen(nBkgnds);
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
TString name( this->bkgndClassName(iBkgnd) );
bkgndClassNames[iBkgnd] = name;
bkgndClassNamesGen[iBkgnd] = "gen"+name;
}
// Loop over the hypotheses and generate the appropriate number of
// events for each one
for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {
// find the category of events (e.g. signal)
const TString& evtCategory(iter->first);
// Type
const TString& type(iter->first);
// Number of events
Int_t nEvtsGen( iter->second.first );
// get the event weight for this category
const Double_t evtWeight( iter->second.second );
for (Int_t iEvt(0); iEvt<nEvtsGen; ++iEvt) {
this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
if (evtCategory == "signal") {
this->setGenNtupleIntegerBranchValue("genSig",1);
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 );
}
// All the generate*Event() methods have to fill in curEvtDecayTime_ and curEvtDecayTimeErr_
// In addition, generateSignalEvent has to decide on the tag and fill in curEvtTagFlv_
genOK = this->generateSignalEvent();
} else {
this->setGenNtupleIntegerBranchValue("genSig",0);
UInt_t bkgndID(0);
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
Int_t gen(0);
if ( bkgndClassNames[iBkgnd] == type ) {
gen = 1;
bkgndID = iBkgnd;
}
this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen );
}
genOK = this->generateBkgndEvent(bkgndID);
}
if (!genOK) {
// If there was a problem with the generation then break out and return.
// The problem model will have adjusted itself so that all should be OK next time.
break;
}
if (this->useDP() == kTRUE) {
this->setDPDtBranchValues(); // store DP, decay time and tagging variables in the ntuple
}
// Store the event's tag and tagging category
this->setGenNtupleIntegerBranchValue("cpEigenvalue", cpEigenValue_);
- this->setGenNtupleDoubleBranchValue(flavTag_->getTrueTagVarName(),curEvtTrueTagFlv_);
- std::vector<TString> tagVarName = flavTag_->getTagVarNames();
- std::vector<TString> mistagVarName = flavTag_->getMistagVarNames();
+ const TString& trueTagVarName { flavTag_->getTrueTagVarName() };
+ if ( trueTagVarName != "" ) {
+ this->setGenNtupleIntegerBranchValue(trueTagVarName, curEvtTrueTagFlv_);
+ }
+
+ if ( cpEigenValue_ == QFS ) {
+ const TString& decayFlvVarName { flavTag_->getDecayFlvVarName() };
+ if ( decayFlvVarName != "" ) {
+ this->setGenNtupleIntegerBranchValue(decayFlvVarName, curEvtDecayFlv_);
+ }
+ }
+
+ const std::vector<TString>& tagVarNames { flavTag_->getTagVarNames() };
+ const std::vector<TString>& mistagVarNames { flavTag_->getMistagVarNames() };
// Loop over the taggers - values set via generateSignalEvent
const ULong_t nTaggers {flavTag_->getNTaggers()};
for (ULong_t i=0; i<nTaggers; ++i){
- this->setGenNtupleIntegerBranchValue(tagVarName[i],curEvtTagFlv_[i]);
- this->setGenNtupleDoubleBranchValue(mistagVarName[i],curEvtMistag_[i]);
+ this->setGenNtupleIntegerBranchValue(tagVarNames[i], curEvtTagFlv_[i]);
+ this->setGenNtupleDoubleBranchValue(mistagVarNames[i], curEvtMistag_[i]);
}
// Store the event number (within this experiment)
// and then increment it
this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum);
++evtNum;
// Write the values into the tree
this->fillGenNtupleBranches();
// Print an occasional progress message
if (iEvt%1000 == 0) {std::cout<<"INFO in LauTimeDepFitModel::genExpt : Generated event number "<<iEvt<<" out of "<<nEvtsGen<<" "<<evtCategory<<" events."<<std::endl;}
}
} //end of loop over species and tagFlv.
if (this->useDP() && genOK) {
sigModelB0bar_->checkToyMC(kTRUE);
sigModelB0_->checkToyMC(kTRUE);
std::cout<<"aSqMaxSet = "<<aSqMaxSet_<<" and aSqMaxVar = "<<aSqMaxVar_<<std::endl;
// Get the fit fractions if they're to be written into the latex table
if (this->writeLatexTable()) {
LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions();
if (fitFracB0bar.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0bar.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0bar[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
LauParArray fitFracB0 = sigModelB0_->getFitFractions();
if (fitFracB0.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j(i); j<nSigComp_; ++j) {
fitFracB0bar_[i][j].value(fitFracB0bar[i][j].value());
fitFracB0_[i][j].value(fitFracB0[i][j].value());
}
}
meanEffB0bar_.value(sigModelB0bar_->getMeanEff().value());
meanEffB0_.value(sigModelB0_->getMeanEff().value());
DPRateB0bar_.value(sigModelB0bar_->getDPRate().value());
DPRateB0_.value(sigModelB0_->getDPRate().value());
}
}
// If we're reusing embedded events or if the generation is being
// reset then clear the lists of used events
if (reuseSignal_ || !genOK) {
if (signalTree_) {
signalTree_->clearUsedList();
}
}
for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
LauEmbeddedData* data = bkgndTree_[bkgndID];
if (reuseBkgnd_[bkgndID] || !genOK) {
if (data) {
data->clearUsedList();
}
}
}
return genOK;
}
Bool_t LauTimeDepFitModel::generateSignalEvent()
{
// Generate signal event, including SCF if necessary.
// DP:DecayTime generation follows.
// If it's ok, we then generate mES, DeltaE, Fisher/NN...
Bool_t genOK(kTRUE);
Bool_t generatedEvent(kFALSE);
Bool_t doSquareDP = kinematicsB0bar_->squareDP();
doSquareDP &= kinematicsB0_->squareDP();
LauKinematics* kinematics(kinematicsB0bar_);
if (this->useDP()) {
if (signalTree_) {
signalTree_->getEmbeddedEvent(kinematics);
//curEvtTagFlv_ = TMath::Nint(signalTree_->getValue("tagFlv"));
curEvtDecayTimeErr_ = signalTree_->getValue(signalDecayTimePdf_->varErrName());
curEvtDecayTime_ = signalTree_->getValue(signalDecayTimePdf_->varName());
if (signalTree_->haveBranch("mcMatch")) {
Int_t match = TMath::Nint(signalTree_->getValue("mcMatch"));
if (match) {
this->setGenNtupleIntegerBranchValue("genTMSig",1);
this->setGenNtupleIntegerBranchValue("genSCFSig",0);
} else {
this->setGenNtupleIntegerBranchValue("genTMSig",0);
this->setGenNtupleIntegerBranchValue("genSCFSig",1);
}
}
} else {
nGenLoop_ = 0;
+
// generate the decay time error (NB the kTRUE forces the generation of a new value)
curEvtDecayTimeErr_ = signalDecayTimePdf_->generateError(kTRUE);
// clear vectors
curEvtTagFlv_.clear();
+ curEvtMistag_.clear();
std::vector<LauParameter*> tageffB0 = flavTag_->getTagEffB0();
std::vector<LauParameter*> tageffB0bar = flavTag_->getTagEffB0bar();
std::vector<LauParameter*> tageffave = flavTag_->getTagEffAve();
std::vector<LauParameter*> tageffdelta = flavTag_->getTagEffDelta();
Double_t tagEffB0(0.), tagEffB0bar(0.);
- curEvtMistag_.clear();
curEvtTrueTagFlv_ = 0;
+ curEvtDecayFlv_ = 0;
// First choose the true tag, accounting for the production asymmetry
// CONVENTION WARNING
Double_t random = LauRandom::randomFun()->Rndm();
if (random <= 0.5 * ( 1.0 - AProd_.unblindValue() ) ) {
curEvtTrueTagFlv_ = 1; // B0 tag
} else {
curEvtTrueTagFlv_ = -1; // B0bar tag
}
// Next generate the tag decisions and per-event mistag probabilities
Double_t randNo(0);
const ULong_t nTaggers { flavTag_->getNTaggers() };
for(ULong_t position{0}; position<nTaggers; ++position){
curEvtMistag_.push_back(flavTag_->getEtaGen(position));
if(flavTag_->getUseAveDelta()){
tagEffB0 = tageffave[position]->unblindValue() + 0.5*tageffdelta[position]->unblindValue();
tagEffB0bar = tageffave[position]->unblindValue() - 0.5*tageffdelta[position]->unblindValue();
} else {
tagEffB0 = tageffB0[position]->unblindValue();
tagEffB0bar = tageffB0bar[position]->unblindValue();
}
if (curEvtTrueTagFlv_ == 1){
randNo = LauRandom::randomFun()->Rndm();
// Try to tag in tageff% of cases
if (randNo <= tagEffB0) {
randNo = LauRandom::randomFun()->Rndm();
// Account for (calibrated) mistag
if (randNo > flavTag_->getLittleOmega(position,1)){
curEvtTagFlv_.push_back(1); // B0 tag
} else {
curEvtTagFlv_.push_back(-1); // B0bar tag
}
} else {
curEvtTagFlv_.push_back(0); // Untagged
}
} else {
randNo = LauRandom::randomFun()->Rndm();
// Try to tag in tageff% of cases
if (randNo <= tagEffB0bar) {
randNo = LauRandom::randomFun()->Rndm();
// Account for (calibrated) mistag
if (randNo > flavTag_->getLittleOmega(position,-1)){
curEvtTagFlv_.push_back(-1); // B0bar tag
} else {
curEvtTagFlv_.push_back(1); // B0 tag
}
} else {
curEvtTagFlv_.push_back(0); // Untagged
}
}
}
// Now generate from the combined DP / decay-time PDF
while (generatedEvent == kFALSE && nGenLoop_ < iterationsMax_) {
// Generate the DP position
Double_t m13Sq{0.0}, m23Sq{0.0};
kinematicsB0bar_->genFlatPhaseSpace(m13Sq, m23Sq);
// Next, calculate the total A and Abar for the given DP position
sigModelB0_->calcLikelihoodInfo(m13Sq, m23Sq);
sigModelB0bar_->calcLikelihoodInfo(m13Sq, m23Sq);
- // Retrieve the amplitudes and efficiency from the dynamics
- const LauComplex& Abar { sigModelB0bar_->getEvtDPAmp() };
- const LauComplex& A { sigModelB0_->getEvtDPAmp() };
- const Double_t dpEff { sigModelB0bar_->getEvtEff() };
-
- // Next calculate the DP terms
- const Double_t aSqSum { A.abs2() + Abar.abs2() };
- const Double_t aSqDif { A.abs2() - Abar.abs2() };
-
- Double_t interTermRe { 0.0 };
- Double_t interTermIm { 0.0 };
- if ( cpEigenValue_ != QFS ) {
- const LauComplex inter { Abar * A.conj() * phiMixComplex_ };
- if ( cpEigenValue_ == CPEven ) {
- interTermIm = 2.0 * inter.im();
- interTermRe = 2.0 * inter.re();
- } else {
- interTermIm = -2.0 * inter.im();
- interTermRe = -2.0 * inter.re();
- }
- }
-
// Generate decay time
const Double_t tMin = signalDecayTimePdf_->minAbscissa();
const Double_t tMax = signalDecayTimePdf_->maxAbscissa();
curEvtDecayTime_ = LauRandom::randomFun()->Rndm()*(tMax-tMin) + tMin;
+ // TODO - should the generation of the decay time error be moved in here?
+ // - I think not because it's generated form its own PDF rather than picking a value from a uniform dist like we're doing for the DP and t
+
// Calculate all the decay time info
signalDecayTimePdf_->calcLikelihoodInfo(curEvtDecayTime_,curEvtDecayTimeErr_);
- // Get the decay time acceptance
- const Double_t dtEff { signalDecayTimePdf_->getEffiTerm() };
- // First get all the decay time terms
+ // Retrieve the amplitudes and efficiency from the dynamics
+ const LauComplex& Abar { sigModelB0bar_->getEvtDPAmp() };
+ const LauComplex& A { sigModelB0_->getEvtDPAmp() };
+ const Double_t dpEff { sigModelB0bar_->getEvtEff() };
+
+ // Also retrieve all the decay time terms
const Double_t dtCos { signalDecayTimePdf_->getCosTerm() };
const Double_t dtSin { signalDecayTimePdf_->getSinTerm() };
const Double_t dtCosh { signalDecayTimePdf_->getCoshTerm() };
const Double_t dtSinh { signalDecayTimePdf_->getSinhTerm() };
+ // and the decay time acceptance
+ const Double_t dtEff { signalDecayTimePdf_->getEffiTerm() };
- // Combine DP and decay-time info for all terms
- // Multiplying the cos and sin terms by the true flavour at production
- const Double_t coshTerm { dtCosh * aSqSum };
- const Double_t sinhTerm { dtSinh * interTermRe };
- const Double_t cosTerm { dtCos * aSqDif * curEvtTrueTagFlv_ };
- const Double_t sinTerm { dtSin * interTermIm * curEvtTrueTagFlv_ };
-
- // Sum to obtain the total and multiply by the efficiency
- const Double_t ASq { ( coshTerm + sinhTerm + cosTerm - sinTerm ) * dpEff * dtEff };
- //std::cout << "Total Amplitude Eff: " << ASq << std::endl;
-
- //Finally we throw the dice to see whether this event should be generated
- //We make a distinction between the likelihood of TM and SCF to tag the SCF events as such
- Double_t randNum = LauRandom::randomFun()->Rndm();
- if (randNum <= ASq/aSqMaxSet_ ) {
- generatedEvent = kTRUE;
- nGenLoop_ = 0;
- if (ASq > aSqMaxVar_) {aSqMaxVar_ = ASq;}
+
+ const Double_t ASq { A.abs2() };
+ const Double_t AbarSq { Abar.abs2() };
+
+ if ( cpEigenValue_ == QFS) {
+ // Calculate the total intensities for each flavour-specific final state
+ const Double_t ATotSq { ( ASq * dtCosh + curEvtTrueTagFlv_ * ASq * dtCos ) * dpEff * dtEff };
+ const Double_t AbarTotSq { ( AbarSq * dtCosh - curEvtTrueTagFlv_ * AbarSq * dtCos ) * dpEff * dtEff };
+ const Double_t ASumSq { ATotSq + AbarTotSq };
+
+ //Finally we throw the dice to see whether this event should be generated (and, if so, which final state)
+ const Double_t randNum = LauRandom::randomFun()->Rndm();
+ if (randNum <= ASumSq / aSqMaxSet_ ) {
+ generatedEvent = kTRUE;
+ nGenLoop_ = 0;
+ if (ASumSq > aSqMaxVar_) {aSqMaxVar_ = ASumSq;}
+
+ if ( randNum <= ATotSq / aSqMaxSet_ ) {
+ curEvtDecayFlv_ = +1;
+ } else {
+ curEvtDecayFlv_ = -1;
+ }
+ } else {
+ nGenLoop_++;
+ }
} else {
- nGenLoop_++;
+ // Calculate the DP terms
+ const Double_t aSqSum { ASq + AbarSq };
+ const Double_t aSqDif { ASq - AbarSq };
+ const LauComplex inter { Abar * A.conj() * phiMixComplex_ };
+ const Double_t interTermRe { ( cpEigenValue_ == CPEven ) ? 2.0 * inter.im() : -2.0 * inter.im() };
+ const Double_t interTermIm { ( cpEigenValue_ == CPEven ) ? 2.0 * inter.re() : -2.0 * inter.re() };
+
+ // Combine DP and decay-time info for all terms
+ // Multiplying the cos and sin terms by the true flavour at production
+ const Double_t coshTerm { aSqSum * dtCosh };
+ const Double_t sinhTerm { interTermRe * dtSinh };
+ const Double_t cosTerm { curEvtTrueTagFlv_ * aSqDif * dtCos };
+ const Double_t sinTerm { curEvtTrueTagFlv_ * interTermIm * dtSin };
+
+ // Sum to obtain the total and multiply by the efficiency
+ const Double_t ATotSq { ( coshTerm + sinhTerm + cosTerm - sinTerm ) * dpEff * dtEff };
+
+ //Finally we throw the dice to see whether this event should be generated
+ const Double_t randNum = LauRandom::randomFun()->Rndm();
+ if (randNum <= ATotSq/aSqMaxSet_ ) {
+ generatedEvent = kTRUE;
+ nGenLoop_ = 0;
+ if (ATotSq > aSqMaxVar_) {aSqMaxVar_ = ATotSq;}
+ } else {
+ nGenLoop_++;
+ }
}
} // end of while !generatedEvent loop
} // end of if (signalTree_) else control
} else {
if ( signalTree_ ) {
signalTree_->getEmbeddedEvent(0);
//curEvtTagFlv_ = TMath::Nint(signalTree_->getValue("tagFlv"));
curEvtDecayTimeErr_ = signalTree_->getValue(signalDecayTimePdf_->varErrName());
curEvtDecayTime_ = signalTree_->getValue(signalDecayTimePdf_->varName());
}
}
// Check whether we have generated the toy MC OK.
if (nGenLoop_ >= iterationsMax_) {
aSqMaxSet_ = 1.01 * aSqMaxVar_;
genOK = kFALSE;
std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Hit max iterations: setting aSqMaxSet_ to "<<aSqMaxSet_<<std::endl;
} else if (aSqMaxVar_ > aSqMaxSet_) {
aSqMaxSet_ = 1.01 * aSqMaxVar_;
genOK = kFALSE;
std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Found a larger ASq value: setting aSqMaxSet_ to "<<aSqMaxSet_<<std::endl;
}
if (genOK) {
//Some variables, like Fisher or NN, might use m13Sq and m23Sq from the kinematics
//kinematicsB0bar_ is up to date, update kinematicsB0_
kinematicsB0_->updateKinematics(kinematicsB0bar_->getm13Sq(), kinematicsB0bar_->getm23Sq() );
this->generateExtraPdfValues(sigExtraPdf_, signalTree_);
}
// Check for problems with the embedding
if (signalTree_ && (signalTree_->nEvents() == signalTree_->nUsedEvents())) {
std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events."<<std::endl;
signalTree_->clearUsedList();
}
return genOK;
}
Bool_t LauTimeDepFitModel::generateBkgndEvent([[maybe_unused]] UInt_t bkgndID)
{
// Generate Bkgnd event
Bool_t genOK(kTRUE);
//LauAbsBkgndDPModel* model(0);
//LauEmbeddedData* embeddedData(0);
//LauPdfList* extraPdfs(0);
//LauKinematics* kinematics(0);
//model = BkgndDPModels_[bkgndID];
//if (this->enableEmbedding()) {
// // find the right embedded data for the current tagging category
// LauTagCatEmbDataMap::const_iterator emb_iter = bkgndTree_[bkgndID].find(curEvtTagCat_);
// embeddedData = (emb_iter != bkgndTree_[bkgndID].end()) ? emb_iter->second : 0;
//}
//extraPdfs = &BkgndPdfs_[bkgndID];
//kinematics = kinematicsB0bar_;
//if (this->useDP()) {
// if (embeddedData) {
// embeddedData->getEmbeddedEvent(kinematics);
// } else {
// if (model == 0) {
// const TString& bkgndClass = this->bkgndClassName(bkgndID);
// std::cerr << "ERROR in LauCPFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << bkgndClass << "\"." << std::endl;
// gSystem->Exit(EXIT_FAILURE);
// }
// genOK = model->generate();
// }
//} else {
// if (embeddedData) {
// embeddedData->getEmbeddedEvent(0);
// }
//}
//if (genOK) {
// this->generateExtraPdfValues(extraPdfs, embeddedData);
//}
//// Check for problems with the embedding
//if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
// const TString& bkgndClass = this->bkgndClassName(bkgndID);
// std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded " << bkgndClass << " events used up, clearing the list of used events." << std::endl;
// embeddedData->clearUsedList();
//}
return genOK;
}
void LauTimeDepFitModel::setupGenNtupleBranches()
{
// Setup the required ntuple branches
this->addGenNtupleDoubleBranch("evtWeight");
this->addGenNtupleIntegerBranch("genSig");
this->addGenNtupleIntegerBranch("cpEigenvalue");
- std::vector<TString> tagVarName = flavTag_->getTagVarNames();
+ const TString& trueTagVarName { flavTag_->getTrueTagVarName() };
+ if ( trueTagVarName != "" ) {
+ this->addGenNtupleIntegerBranch(trueTagVarName);
+ }
+
+ if ( cpEigenValue_ == QFS ) {
+ const TString& decayFlvVarName { flavTag_->getDecayFlvVarName() };
+ if ( decayFlvVarName != "" ) {
+ this->addGenNtupleIntegerBranch(decayFlvVarName);
+ }
+ }
+
+ const std::vector<TString>& tagVarNames { flavTag_->getTagVarNames() };
+ const std::vector<TString>& mistagVarNames { flavTag_->getMistagVarNames() };
const ULong_t nTaggers {flavTag_->getNTaggers()};
for (ULong_t position{0}; position<nTaggers; ++position){
- this->addGenNtupleIntegerBranch(tagVarName[position]);
+ this->addGenNtupleIntegerBranch(tagVarNames[position]);
+ this->addGenNtupleDoubleBranch(mistagVarNames[position]);
}
if (this->useDP() == kTRUE) {
// Let's add the decay time variables.
this->addGenNtupleDoubleBranch(signalDecayTimePdf_->varName());
this->addGenNtupleDoubleBranch(signalDecayTimePdf_->varErrName());
this->addGenNtupleDoubleBranch("m12");
this->addGenNtupleDoubleBranch("m23");
this->addGenNtupleDoubleBranch("m13");
this->addGenNtupleDoubleBranch("m12Sq");
this->addGenNtupleDoubleBranch("m23Sq");
this->addGenNtupleDoubleBranch("m13Sq");
this->addGenNtupleDoubleBranch("cosHel12");
this->addGenNtupleDoubleBranch("cosHel23");
this->addGenNtupleDoubleBranch("cosHel13");
if (kinematicsB0bar_->squareDP() && kinematicsB0_->squareDP()) {
this->addGenNtupleDoubleBranch("mPrime");
this->addGenNtupleDoubleBranch("thPrime");
}
// Can add the real and imaginary parts of the B0 and B0bar total
// amplitudes seen in the generation (restrict this with a flag
// that defaults to false)
if ( storeGenAmpInfo_ ) {
this->addGenNtupleDoubleBranch("reB0Amp");
this->addGenNtupleDoubleBranch("imB0Amp");
this->addGenNtupleDoubleBranch("reB0barAmp");
this->addGenNtupleDoubleBranch("imB0barAmp");
}
}
// Let's look at the extra variables for signal in one of the tagging categories
if ( sigExtraPdf_ ) {
for (LauPdfList::const_iterator pdf_iter = sigExtraPdf_->begin(); pdf_iter != sigExtraPdf_->end(); ++pdf_iter) {
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
this->addGenNtupleDoubleBranch( (*var_iter) );
}
}
}
}
}
void LauTimeDepFitModel::setDPDtBranchValues()
{
// Store the decay time variables.
this->setGenNtupleDoubleBranchValue(signalDecayTimePdf_->varName(),curEvtDecayTime_);
this->setGenNtupleDoubleBranchValue(signalDecayTimePdf_->varErrName(),curEvtDecayTimeErr_);
// CONVENTION WARNING
// TODO check - for now use B0 for any tags
//LauKinematics* kinematics(0);
//if (curEvtTagFlv_[position]<0) {
LauKinematics* kinematics = kinematicsB0_;
//} else {
// kinematics = kinematicsB0bar_;
//}
// Store all the DP information
this->setGenNtupleDoubleBranchValue("m12", kinematics->getm12());
this->setGenNtupleDoubleBranchValue("m23", kinematics->getm23());
this->setGenNtupleDoubleBranchValue("m13", kinematics->getm13());
this->setGenNtupleDoubleBranchValue("m12Sq", kinematics->getm12Sq());
this->setGenNtupleDoubleBranchValue("m23Sq", kinematics->getm23Sq());
this->setGenNtupleDoubleBranchValue("m13Sq", kinematics->getm13Sq());
this->setGenNtupleDoubleBranchValue("cosHel12", kinematics->getc12());
this->setGenNtupleDoubleBranchValue("cosHel23", kinematics->getc23());
this->setGenNtupleDoubleBranchValue("cosHel13", kinematics->getc13());
if (kinematics->squareDP()) {
this->setGenNtupleDoubleBranchValue("mPrime", kinematics->getmPrime());
this->setGenNtupleDoubleBranchValue("thPrime", kinematics->getThetaPrime());
}
// Can add the real and imaginary parts of the B0 and B0bar total
// amplitudes seen in the generation (restrict this with a flag
// that defaults to false)
if ( storeGenAmpInfo_ ) {
if ( this->getGenNtupleIntegerBranchValue("genSig")==1 ) {
LauComplex Abar = sigModelB0bar_->getEvtDPAmp();
LauComplex A = sigModelB0_->getEvtDPAmp();
this->setGenNtupleDoubleBranchValue("reB0Amp", A.re());
this->setGenNtupleDoubleBranchValue("imB0Amp", A.im());
this->setGenNtupleDoubleBranchValue("reB0barAmp", Abar.re());
this->setGenNtupleDoubleBranchValue("imB0barAmp", Abar.im());
} else {
this->setGenNtupleDoubleBranchValue("reB0Amp", 0.0);
this->setGenNtupleDoubleBranchValue("imB0Amp", 0.0);
this->setGenNtupleDoubleBranchValue("reB0barAmp", 0.0);
this->setGenNtupleDoubleBranchValue("imB0barAmp", 0.0);
}
}
}
void LauTimeDepFitModel::generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData)
{
// CONVENTION WARNING
LauKinematics* kinematics = kinematicsB0_;
//LauKinematics* kinematics(0);
//if (curEvtTagFlv_<0) {
// kinematics = kinematicsB0_;
//} else {
// kinematics = kinematicsB0bar_;
//}
// Generate from the extra PDFs
if (extraPdfs) {
for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
LauFitData genValues;
if (embeddedData) {
genValues = embeddedData->getValues( (*pdf_iter)->varNames() );
} else {
genValues = (*pdf_iter)->generate(kinematics);
}
for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) {
TString varName = var_iter->first;
if ( varName != "m13Sq" && varName != "m23Sq" ) {
Double_t value = var_iter->second;
this->setGenNtupleDoubleBranchValue(varName,value);
}
}
}
}
}
void LauTimeDepFitModel::propagateParUpdates()
{
// Update the complex mixing phase
if (useSinCos_) {
phiMixComplex_.setRealPart(cosPhiMix_.unblindValue());
phiMixComplex_.setImagPart(-1.0*sinPhiMix_.unblindValue());
} else {
phiMixComplex_.setRealPart(TMath::Cos(-1.0*phiMix_.unblindValue()));
phiMixComplex_.setImagPart(TMath::Sin(-1.0*phiMix_.unblindValue()));
}
// Update the total normalisation for the signal likelihood
if (this->useDP() == kTRUE) {
this->updateCoeffs();
sigModelB0bar_->updateCoeffs(coeffsB0bar_);
sigModelB0_->updateCoeffs(coeffsB0_);
this->calcInterTermNorm();
}
// Update the decay time normalisation
//if ( signalDecayTimePdf_ ) {
// TODO - at present this isn't needed since we're brute force updating everything for every fit iteration anyway
// - should make this intelligent (only update if certain parameters are floating and have changed in the last iteration) - probably put the intelligence inside LauDecayTimePdf
// - probably needs a new function - calcNorm is required to run per-event in some scenarios, so the LauDecayTimePdf object should work this out for itself
// - this function could recalculate and recache everything that depends on the parameters that have changed
// - then calcLikelihoodInfo(iEvt) can revert to just retrieving from the cache
//signalDecayTimePdf_->calcNorm();
//}
// TODO
// - maybe also need to add an update of the background decay time PDFs here
// Update the signal events from the background numbers if not doing an extended fit
// And update the tagging category fractions
this->updateSigEvents();
}
void LauTimeDepFitModel::updateSigEvents()
{
// The background parameters will have been set from Minuit.
// We need to update the signal events using these.
if (!this->doEMLFit()) {
Double_t nTotEvts = this->eventsPerExpt();
Double_t signalEvents = nTotEvts;
signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
LauAbsRValue* nBkgndEvents = (*iter);
if ( nBkgndEvents->isLValue() ) {
LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
yield->range(-2.0*nTotEvts,2.0*nTotEvts);
}
}
// Subtract background events (if any) from signal.
if (usingBkgnd_ == kTRUE) {
for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
signalEvents -= (*iter)->value();
}
}
if ( ! signalEvents_->fixed() ) {
signalEvents_->value(signalEvents);
}
}
}
void LauTimeDepFitModel::cacheInputFitVars()
{
// Fill the internal data trees of the signal and background models.
// Note that we store the events of both charges in both the
// negative and the positive models. It's only later, at the stage
// when the likelihood is being calculated, that we separate them.
LauFitDataTree* inputFitData = this->fitData();
evtCPEigenVals_.clear();
const Bool_t hasCPEV = ( (cpevVarName_ != "") && inputFitData->haveBranch( cpevVarName_ ) );
UInt_t nEvents = inputFitData->nEvents();
evtCPEigenVals_.reserve( nEvents );
LauFitData::const_iterator fitdata_iter;
for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
const LauFitData& dataValues = inputFitData->getData(iEvt);
// if the CP-eigenvalue is in the data use those, otherwise use the default
if ( hasCPEV ) {
fitdata_iter = dataValues.find( cpevVarName_ );
const Int_t cpEV = static_cast<Int_t>( fitdata_iter->second );
if ( cpEV == 1 ) {
cpEigenValue_ = CPEven;
} else if ( cpEV == -1 ) {
cpEigenValue_ = CPOdd;
} else if ( cpEV == 0 ) {
cpEigenValue_ = QFS;
} else {
std::cerr<<"WARNING in LauTimeDepFitModel::cacheInputFitVars : Unknown value: "<<cpEV<<" for CP eigenvalue, setting to CP-even"<<std::endl;
cpEigenValue_ = CPEven;
}
}
evtCPEigenVals_.push_back( cpEigenValue_ );
}
// We'll cache the DP amplitudes at the end because we'll
// append some points that the other PDFs won't deal with.
// Flavour tagging information
flavTag_->cacheInputFitVars(inputFitData);
if (this->useDP() == kTRUE) {
// DecayTime and SigmaDecayTime
signalDecayTimePdf_->cacheInfo(*inputFitData);
}
// ...and then the extra PDFs
if (sigExtraPdf_){
this->cacheInfo((*sigExtraPdf_), *inputFitData);
}
if(usingBkgnd_ == kTRUE){
for (LauBkgndPdfsList::iterator iter = BkgndPdfs_.begin(); iter != BkgndPdfs_.end(); ++iter) {
this->cacheInfo((*iter), *inputFitData);
}
}
if (this->useDP() == kTRUE) {
sigModelB0bar_->fillDataTree(*inputFitData);
sigModelB0_->fillDataTree(*inputFitData);
if (usingBkgnd_ == kTRUE) {
for (LauBkgndDPModelList::iterator iter = BkgndDPModels_.begin(); iter != BkgndDPModels_.end(); ++iter) {
(*iter)->fillDataTree(*inputFitData);
}
}
}
}
Double_t LauTimeDepFitModel::getTotEvtLikelihood(const UInt_t iEvt)
{
-
- // Find out whether the tag-side B was a B0 or a B0bar.
- curEvtTagFlv_ = flavTag_->getCurEvtTagFlv();
-
// Get the CP eigenvalue of the current event
cpEigenValue_ = evtCPEigenVals_[iEvt];
// Get the DP and DecayTime likelihood for signal (TODO and eventually backgrounds)
this->getEvtDPDtLikelihood(iEvt);
// Get the flavour tagging likelihood from eta PDFs (per tagging category - TODO backgrounds to come later)
sigFlavTagLike_ = 1.0;
//this->getEvtFlavTagLikelihood(iEvt);
// Get the combined extra PDFs likelihood for signal (TODO and eventually backgrounds)
this->getEvtExtraLikelihoods(iEvt);
// Construct the total likelihood for signal, qqbar and BBbar backgrounds
Double_t sigLike = sigDPLike_ * sigFlavTagLike_ * sigExtraLike_;
//std::cout << "DP like = " << sigDPLike_ << std::endl;
//std::cout << "flav tag like = " << sigFlavTagLike_ << std::endl;
//std::cout << "extra like = " << sigExtraLike_ << std::endl;
// TODO
Double_t signalEvents = signalEvents_->unblindValue();
if (this->useDP() == kFALSE) {
//signalEvents *= 0.5 * (1.0 + curEvtTagFlv_ * signalAsym_->unblindValue());
}
if ( ! signalEvents_->fixed() ) {
sigLike *= signalEvents;
}
return sigLike;
}
Double_t LauTimeDepFitModel::getEventSum() const
{
Double_t eventSum(0.0);
eventSum += signalEvents_->unblindValue();
return eventSum;
}
void LauTimeDepFitModel::getEvtDPDtLikelihood(const UInt_t iEvt)
{
// Function to return the signal and background likelihoods for the
// Dalitz plot for the given event evtNo.
if ( ! this->useDP() ) {
// There's always going to be a term in the likelihood for the
// signal, so we'd better not zero it.
sigDPLike_ = 1.0;
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
if (usingBkgnd_ == kTRUE) {
bkgndDPLike_[bkgndID] = 1.0;
} else {
bkgndDPLike_[bkgndID] = 0.0;
}
}
return;
}
// Calculate event quantities
+ // Get the DP dynamics, decay time, and flavour tagging to calculate
+ // everything required for the likelihood calculation
- // Get the dynamics to calculate everything required for the likelihood calculation
sigModelB0bar_->calcLikelihoodInfo(iEvt);
sigModelB0_->calcLikelihoodInfo(iEvt);
- // Background part
- // TODO add them into the actual Likelihood calculations
- const UInt_t nBkgnds = this->nBkgndClasses();
- for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
- if (usingBkgnd_ == kTRUE) {
- bkgndDPLike_[bkgndID] = BkgndDPModels_[bkgndID]->getLikelihood(iEvt);
- } else {
- bkgndDPLike_[bkgndID] = 0.0;
- }
- }
+ signalDecayTimePdf_->calcLikelihoodInfo(iEvt);
+
+ flavTag_->updateEventInfo(iEvt);
+
// Retrieve the amplitudes and efficiency from the dynamics
- const LauComplex& Abar { sigModelB0bar_->getEvtDPAmp() };
- const LauComplex& A { sigModelB0_->getEvtDPAmp() };
+ LauComplex Abar { sigModelB0bar_->getEvtDPAmp() };
+ LauComplex A { sigModelB0_->getEvtDPAmp() };
+
const Double_t dpEff { sigModelB0bar_->getEvtEff() };
+
+ // If this is a QFS decay, one of the DP amplitudes needs to be zeroed
+ if (cpEigenValue_ == QFS){
+ curEvtDecayFlv_ = flavTag_->getCurEvtDecayFlv();
+ if ( curEvtDecayFlv_ == +1 ) {
+ Abar.zero();
+ } else if ( curEvtDecayFlv_ == -1 ) {
+ A.zero();
+ }
+ }
+
+
// Next calculate the DP terms
const Double_t aSqSum { A.abs2() + Abar.abs2() };
const Double_t aSqDif { A.abs2() - Abar.abs2() };
Double_t interTermRe { 0.0 };
Double_t interTermIm { 0.0 };
if ( cpEigenValue_ != QFS ) {
const LauComplex inter { Abar * A.conj() * phiMixComplex_ };
if ( cpEigenValue_ == CPEven ) {
interTermIm = 2.0 * inter.im();
interTermRe = 2.0 * inter.re();
} else {
interTermIm = -2.0 * inter.im();
interTermRe = -2.0 * inter.re();
}
}
// First get all the decay time terms
- signalDecayTimePdf_->calcLikelihoodInfo(iEvt);
// TODO Backgrounds
// Get the decay time acceptance
const Double_t dtEff { signalDecayTimePdf_->getEffiTerm() };
// Get all the decay time terms
const Double_t dtCos { signalDecayTimePdf_->getCosTerm() };
const Double_t dtSin { signalDecayTimePdf_->getSinTerm() };
const Double_t dtCosh { signalDecayTimePdf_->getCoshTerm() };
const Double_t dtSinh { signalDecayTimePdf_->getSinhTerm() };
// Get the decay time error term
const Double_t dtErrLike { signalDecayTimePdf_->getErrTerm() };
// Get flavour tagging terms
- flavTag_->updateEventInfo(iEvt);
-
Double_t omega{1.0};
Double_t omegabar{1.0};
const ULong_t nTaggers { flavTag_->getNTaggers() };
for (ULong_t position{0}; position<nTaggers; ++position){
omega *= flavTag_->getCapitalOmega(position,+1);
omegabar *= flavTag_->getCapitalOmega(position,-1);
}
const Double_t prodAsym { AProd_.unblindValue() };
const Double_t ftOmegaHyp { ((1.0 - prodAsym)*omega + (1.0 + prodAsym)*omegabar) };
const Double_t ftOmegaTrig { ((1.0 - prodAsym)*omega - (1.0 + prodAsym)*omegabar) };
- Double_t coshTerm { dtCosh * ftOmegaHyp * aSqSum };
- Double_t sinhTerm { dtSinh * ftOmegaHyp * interTermRe };
- Double_t cosTerm { dtCos * ftOmegaTrig * aSqDif };
- Double_t sinTerm { dtSin * ftOmegaTrig * interTermIm };
-
- curEvtTrueTagFlv_ = flavTag_->getCurEvtTrueTagFlv();
-
- if (curEvtTrueTagFlv_ != 0 && cpEigenValue_ == QFS){
- cosTerm *= curEvtTrueTagFlv_;
- sinTerm *= curEvtTrueTagFlv_;
- }
+ const Double_t coshTerm { ftOmegaHyp * dtCosh * aSqSum };
+ const Double_t sinhTerm { ftOmegaHyp * dtSinh * interTermRe };
+ const Double_t cosTerm { ftOmegaTrig * dtCos * aSqDif };
+ const Double_t sinTerm { ftOmegaTrig * dtSin * interTermIm };
// Combine all terms to get the total amplitude squared
const Double_t ASq { coshTerm + sinhTerm + cosTerm - sinTerm };
// Calculate the DP and time normalisation
const Double_t normASqSum { sigModelB0_->getDPNorm() + sigModelB0bar_->getDPNorm() };
const Double_t normASqDiff { sigModelB0_->getDPNorm() - sigModelB0bar_->getDPNorm() };
Double_t normInterTermRe { 0.0 };
Double_t normInterTermIm { 0.0 };
if ( cpEigenValue_ != QFS ) {
// TODO - double check this sign flipping here (it's presumably right but...)
normInterTermRe = ( cpEigenValue_ == CPOdd ) ? -1.0 * interTermReNorm_ : interTermReNorm_;
normInterTermIm = ( cpEigenValue_ == CPOdd ) ? -1.0 * interTermImNorm_ : interTermImNorm_;
}
- //const Double_t normExpTerm { signalDecayTimePdf_->getNormTermExp() };
const Double_t normCoshTerm { signalDecayTimePdf_->getNormTermCosh() };
const Double_t normSinhTerm { signalDecayTimePdf_->getNormTermSinh() };
const Double_t normCosTerm { signalDecayTimePdf_->getNormTermCos() };
const Double_t normSinTerm { signalDecayTimePdf_->getNormTermSin() };
- Double_t asymPart { - 2.0 * prodAsym * ( normASqDiff * normCosTerm + normInterTermIm * normSinTerm ) };
- // TODO - double check what to do about the true flavour here
- if (curEvtTrueTagFlv_ != 0 && cpEigenValue_ == QFS){
- asymPart *= curEvtTrueTagFlv_;
- }
+ const Double_t normHyp { normASqSum * normCoshTerm + normInterTermRe * normSinhTerm };
+ const Double_t normTrig { - prodAsym * ( normASqDiff * normCosTerm + normInterTermIm * normSinTerm ) };
// Combine all terms to get the total normalisation
- const Double_t norm { normASqSum * normCoshTerm + normInterTermRe * normSinhTerm + asymPart };
+ const Double_t norm { 2.0 * ( normHyp + normTrig ) };
// Multiply the squared-amplitude by the efficiency (DP and decay time) and decay-time error likelihood
// and normalise to obtain the signal likelihood
sigDPLike_ = ( ASq * dpEff * dtEff * dtErrLike ) / norm;
+
+
+ // Background part
+ // TODO add them into the actual Likelihood calculations
+ const UInt_t nBkgnds = this->nBkgndClasses();
+ for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
+ if (usingBkgnd_ == kTRUE) {
+ bkgndDPLike_[bkgndID] = BkgndDPModels_[bkgndID]->getLikelihood(iEvt);
+ } else {
+ bkgndDPLike_[bkgndID] = 0.0;
+ }
+ }
}
void LauTimeDepFitModel::getEvtExtraLikelihoods(const UInt_t iEvt)
{
// Function to return the signal and background likelihoods for the
// extra variables for the given event evtNo.
sigExtraLike_ = 1.0; //There's always a likelihood term for signal, so we better not zero it.
// First, those independent of the tagging of the event:
// signal
if (sigExtraPdf_) {
sigExtraLike_ = this->prodPdfValue( (*sigExtraPdf_), iEvt );
}
// Background
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
if (usingBkgnd_) {
bkgndExtraLike_[bkgndID] = this->prodPdfValue( BkgndPdfs_[bkgndID], iEvt );
} else {
bkgndExtraLike_[bkgndID] = 0.0;
}
}
}
void LauTimeDepFitModel::getEvtFlavTagLikelihood(const UInt_t iEvt)
{
// Function to return the signal and background likelihoods for the
// extra variables for the given event evtNo.
sigFlavTagLike_ = 1.0; //There's always a likelihood term for signal, so we better not zero it.
// Loop over taggers
const ULong_t nTaggers { flavTag_->getNTaggers() };
for (ULong_t position{0}; position<nTaggers; ++position){
if (sigFlavTagPdf_[position]) {
sigFlavTagPdf_[position]->calcLikelihoodInfo(iEvt);
sigFlavTagLike_ = sigFlavTagPdf_[position]->getLikelihood();
}
}
if (sigFlavTagLike_<=0){
std::cout<<"INFO in LauTimeDepFitModel::getEvtFlavTagLikelihood : Event with 0 FlavTag Liklihood"<<std::endl;
}
// TODO Add in the background components too
}
void LauTimeDepFitModel::updateCoeffs()
{
coeffsB0bar_.clear(); coeffsB0_.clear();
coeffsB0bar_.reserve(nSigComp_); coeffsB0_.reserve(nSigComp_);
for (UInt_t i = 0; i < nSigComp_; ++i) {
coeffsB0bar_.push_back(coeffPars_[i]->antiparticleCoeff());
coeffsB0_.push_back(coeffPars_[i]->particleCoeff());
}
}
void LauTimeDepFitModel::checkMixingPhase()
{
Double_t phase = phiMix_.value();
Double_t genPhase = phiMix_.genValue();
// Check now whether the phase lies in the right range (-pi to pi).
Bool_t withinRange(kFALSE);
while (withinRange == kFALSE) {
if (phase > -LauConstants::pi && phase < LauConstants::pi) {
withinRange = kTRUE;
} else {
// Not within the specified range
if (phase > LauConstants::pi) {
phase -= LauConstants::twoPi;
} else if (phase < -LauConstants::pi) {
phase += LauConstants::twoPi;
}
}
}
// A further problem can occur when the generated phase is close to -pi or pi.
// The phase can wrap over to the other end of the scale -
// this leads to artificially large pulls so we wrap it back.
Double_t diff = phase - genPhase;
if (diff > LauConstants::pi) {
phase -= LauConstants::twoPi;
} else if (diff < -LauConstants::pi) {
phase += LauConstants::twoPi;
}
// finally store the new value in the parameter
// and update the pull
phiMix_.value(phase);
phiMix_.updatePull();
}
void LauTimeDepFitModel::embedSignal(const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
{
if (signalTree_) {
std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Already embedding signal from file."<<std::endl;
return;
}
if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
std::cerr<<"WARNING in LauTimeDepFitModel::embedSignal : Conflicting options provided, will not reuse events at all."<<std::endl;
reuseEventsWithinExperiment = kFALSE;
}
signalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
Bool_t dataOK = signalTree_->findBranches();
if (!dataOK) {
delete signalTree_; signalTree_ = 0;
std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Problem creating data tree for embedding."<<std::endl;
return;
}
reuseSignal_ = reuseEventsWithinEnsemble;
}
void LauTimeDepFitModel::embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
{
if ( ! this->validBkgndClass( bkgndClass ) ) {
std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
return;
}
UInt_t bkgndID = this->bkgndClassID( bkgndClass );
LauEmbeddedData* bkgTree = bkgndTree_[bkgndID];
if (bkgTree) {
std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Already embedding background from a file." << std::endl;
return;
}
bkgTree = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
Bool_t dataOK = bkgTree->findBranches();
if (!dataOK) {
delete bkgTree; bkgTree = 0;
std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Problem creating data tree for embedding." << std::endl;
return;
}
reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;
if (this->enableEmbedding() == kFALSE) {
this->enableEmbedding(kTRUE);
}
}
void LauTimeDepFitModel::setupSPlotNtupleBranches()
{
// add branches for storing the experiment number and the number of
// the event within the current experiment
this->addSPlotNtupleIntegerBranch("iExpt");
this->addSPlotNtupleIntegerBranch("iEvtWithinExpt");
// Store the efficiency of the event (for inclusive BF calculations).
if (this->storeDPEff()) {
this->addSPlotNtupleDoubleBranch("efficiency");
}
// Store the total event likelihood for each species.
this->addSPlotNtupleDoubleBranch("sigTotalLike");
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
TString name( this->bkgndClassName(iBkgnd) );
name += "TotalLike";
this->addSPlotNtupleDoubleBranch(name);
}
}
// Store the DP likelihoods
if (this->useDP()) {
this->addSPlotNtupleDoubleBranch("sigDPLike");
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
TString name( this->bkgndClassName(iBkgnd) );
name += "DPLike";
this->addSPlotNtupleDoubleBranch(name);
}
}
}
// Store the likelihoods for each extra PDF
const LauPdfList* pdfList( sigExtraPdf_ );
this->addSPlotNtupleBranches(pdfList, "sig");
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
const TString& bkgndClass = this->bkgndClassName(iBkgnd);
const LauPdfList* pdfList2 = &(BkgndPdfs_[iBkgnd]);
this->addSPlotNtupleBranches(pdfList2, bkgndClass);
}
}
}
void LauTimeDepFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix)
{
if (!extraPdfs) {
return;
}
// Loop through each of the PDFs
for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
// Count the number of input variables that are not
// DP variables (used in the case where there is DP
// dependence for e.g. MVA)
UInt_t nVars(0);
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
++nVars;
}
}
if ( nVars == 1 ) {
// If the PDF only has one variable then
// simply add one branch for that variable
TString varName = (*pdf_iter)->varName();
TString name(prefix);
name += varName;
name += "Like";
this->addSPlotNtupleDoubleBranch(name);
} else if ( nVars == 2 ) {
// If the PDF has two variables then we
// need a branch for them both together and
// branches for each
TString allVars("");
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
allVars += (*var_iter);
TString name(prefix);
name += (*var_iter);
name += "Like";
this->addSPlotNtupleDoubleBranch(name);
}
TString name(prefix);
name += allVars;
name += "Like";
this->addSPlotNtupleDoubleBranch(name);
} else {
std::cerr<<"WARNING in LauTimeDepFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs."<<std::endl;
}
}
}
Double_t LauTimeDepFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt)
{
// Store the PDF value for each variable in the list
Double_t totalLike(1.0);
Double_t extraLike(0.0);
if ( !extraPdfs ) {
return totalLike;
}
for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
// calculate the likelihood for this event
(*pdf_iter)->calcLikelihoodInfo(iEvt);
extraLike = (*pdf_iter)->getLikelihood();
totalLike *= extraLike;
// Count the number of input variables that are not
// DP variables (used in the case where there is DP
// dependence for e.g. MVA)
UInt_t nVars(0);
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
++nVars;
}
}
if ( nVars == 1 ) {
// If the PDF only has one variable then
// simply store the value for that variable
TString varName = (*pdf_iter)->varName();
TString name(prefix);
name += varName;
name += "Like";
this->setSPlotNtupleDoubleBranchValue(name, extraLike);
} else if ( nVars == 2 ) {
// If the PDF has two variables then we
// store the value for them both together
// and for each on their own
TString allVars("");
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
allVars += (*var_iter);
TString name(prefix);
name += (*var_iter);
name += "Like";
Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) );
this->setSPlotNtupleDoubleBranchValue(name, indivLike);
}
TString name(prefix);
name += allVars;
name += "Like";
this->setSPlotNtupleDoubleBranchValue(name, extraLike);
} else {
std::cerr<<"WARNING in LauAllFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs."<<std::endl;
}
}
return totalLike;
}
LauSPlot::NameSet LauTimeDepFitModel::variableNames() const
{
LauSPlot::NameSet nameSet;
if (this->useDP()) {
nameSet.insert("DP");
}
for (LauPdfList::const_iterator pdf_iter = sigExtraPdf_->begin(); pdf_iter != sigExtraPdf_->end(); ++pdf_iter) {
// Loop over the variables involved in each PDF
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
// If they are not DP coordinates then add them
if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
nameSet.insert( (*var_iter) );
}
}
}
return nameSet;
}
LauSPlot::NumbMap LauTimeDepFitModel::freeSpeciesNames() const
{
LauSPlot::NumbMap numbMap;
if (!signalEvents_->fixed() && this->doEMLFit()) {
numbMap["sig"] = signalEvents_->genValue();
}
if ( usingBkgnd_ ) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
const TString& bkgndClass = this->bkgndClassName(iBkgnd);
const LauAbsRValue* par = bkgndEvents_[iBkgnd];
if (!par->fixed()) {
numbMap[bkgndClass] = par->genValue();
if ( ! par->isLValue() ) {
std::cerr << "WARNING in LauTimeDepFitModel::freeSpeciesNames : \"" << par->name() << "\" is a LauFormulaPar, which implies it is perhaps not entirely free to float in the fit, so the sWeight calculation may not be reliable" << std::endl;
}
}
}
}
return numbMap;
}
LauSPlot::NumbMap LauTimeDepFitModel::fixdSpeciesNames() const
{
LauSPlot::NumbMap numbMap;
if (signalEvents_->fixed() && this->doEMLFit()) {
numbMap["sig"] = signalEvents_->genValue();
}
if ( usingBkgnd_ ) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
const TString& bkgndClass = this->bkgndClassName(iBkgnd);
const LauAbsRValue* par = bkgndEvents_[iBkgnd];
if (par->fixed()) {
numbMap[bkgndClass] = par->genValue();
}
}
}
return numbMap;
}
LauSPlot::TwoDMap LauTimeDepFitModel::twodimPDFs() const
{
LauSPlot::TwoDMap twodimMap;
const LauPdfList* pdfList = sigExtraPdf_;
for (LauPdfList::const_iterator pdf_iter = pdfList->begin(); pdf_iter != pdfList->end(); ++pdf_iter) {
// Count the number of input variables that are not DP variables
UInt_t nVars(0);
for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
++nVars;
}
}
if ( nVars == 2 ) {
twodimMap.insert( std::make_pair( "sig", std::make_pair( (*pdf_iter)->varNames()[0], (*pdf_iter)->varNames()[1] ) ) );
}
}
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
const TString& bkgndClass = this->bkgndClassName(iBkgnd);
const LauPdfList& pdfList2 = BkgndPdfs_[iBkgnd];
for (LauPdfList::const_iterator pdf_iter = pdfList2.begin(); pdf_iter != pdfList2.end(); ++pdf_iter) {
// Count the number of input variables that are not DP variables
UInt_t nVars(0);
std::vector<TString> varNames = (*pdf_iter)->varNames();
for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
++nVars;
}
}
if ( nVars == 2 ) {
twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) );
}
}
}
}
return twodimMap;
}
void LauTimeDepFitModel::storePerEvtLlhds()
{
std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Storing per-event likelihood values..."<<std::endl;
LauFitDataTree* inputFitData = this->fitData();
// if we've not been using the DP model then we need to cache all
// the info here so that we can get the efficiency from it
if (!this->useDP() && this->storeDPEff()) {
sigModelB0bar_->initialise(coeffsB0bar_);
sigModelB0_->initialise(coeffsB0_);
sigModelB0bar_->fillDataTree(*inputFitData);
sigModelB0_->fillDataTree(*inputFitData);
}
UInt_t evtsPerExpt(this->eventsPerExpt());
LauIsobarDynamics* sigModel(sigModelB0bar_);
for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) {
// Find out whether we have B0bar or B0
flavTag_->updateEventInfo(iEvt);
curEvtTagFlv_ = flavTag_->getCurEvtTagFlv();
curEvtMistag_ = flavTag_->getCurEvtMistag();
// the DP information
this->getEvtDPDtLikelihood(iEvt);
if (this->storeDPEff()) {
if (!this->useDP()) {
sigModel->calcLikelihoodInfo(iEvt);
}
this->setSPlotNtupleDoubleBranchValue("efficiency",sigModel->getEvtEff());
}
if (this->useDP()) {
sigTotalLike_ = sigDPLike_;
this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_);
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
TString name = this->bkgndClassName(iBkgnd);
name += "DPLike";
this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]);
}
}
} else {
sigTotalLike_ = 1.0;
}
// the signal PDF values
sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigExtraPdf_, "sig", iEvt);
// the background PDF values
LauBkgndPdfsList* bkgndPdfs(0);
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
const TString& bkgndClass = this->bkgndClassName(iBkgnd);
LauPdfList& pdfs = (*bkgndPdfs)[iBkgnd];
bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt);
}
}
// the total likelihoods
this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
if (usingBkgnd_) {
const UInt_t nBkgnds = this->nBkgndClasses();
for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
TString name = this->bkgndClassName(iBkgnd);
name += "TotalLike";
this->setSPlotNtupleDoubleBranchValue(name,bkgndTotalLike_[iBkgnd]);
}
}
// fill the tree
this->fillSPlotNtupleBranches();
}
std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Finished storing per-event likelihood values."<<std::endl;
}
void LauTimeDepFitModel::weightEvents( const TString& /*dataFileName*/, const TString& /*dataTreeName*/ )
{
std::cerr << "ERROR in LauTimeDepFitModel::weightEvents : Method not available for this fit model." << std::endl;
return;
}
void LauTimeDepFitModel::savePDFPlots(const TString& /*label*/)
{
}
void LauTimeDepFitModel::savePDFPlotsWave(const TString& /*label*/, const Int_t& /*spin*/)
{
}
void LauTimeDepFitModel::setBkgdFlavTagPdfs(const TString name, LauAbsPdf* pdf)
{
// TODO - when backgrounds are added - implement a check that "name" really is the name of a background category
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setBkgdFlavTagPdfs : The PDF pointer is null."<<std::endl;
return;
}
bkgdFlavTagPdf_[name] = pdf;
}