Page MenuHomeHEPForge
Files F19215088

D104.1759092863.diff
No OneTemporary
Actions

Size
23 KB
Referenced Files
None
Subscribers
None

D104.1759092863.diff
View Options

diff --git a/doc/ReleaseNotes.md b/doc/ReleaseNotes.md
--- a/doc/ReleaseNotes.md
+++ b/doc/ReleaseNotes.md
@@ -1,5 +1,13 @@
# Laura++ release notes
+16th Nov 2023 Thomas Latham
+* Add functions to fit models that return:
+ - signal DP amplitudes
+ - signal and background DP likelihood(s)
+* Add example to plot histograms of likelihood values
+* Caveat: SCF signal is not currently handled - an appropriate warning is printed
+* See https://phab.hepforge.org/T227 for more details
+
16th Oct 2023 Thomas Latham
* Remove all uses of deprecated ROOT ClassImp macro
* Generate ROOT dictionary only for classes: LauAbsRValue, LauBlind and LauParameter
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -22,6 +22,7 @@
MergeDataFiles
mixedSampleTest
PlotKMatrixTAmp
+ PlotLikelihood
PlotResults
point2PointTestSample
QuasiFlatSqDalitz
diff --git a/examples/PlotLikelihood.cc b/examples/PlotLikelihood.cc
new file mode 100644
--- /dev/null
+++ b/examples/PlotLikelihood.cc
@@ -0,0 +1,194 @@
+
+/*
+Copyright 2023 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
+*/
+
+#include <cstdlib>
+#include <iostream>
+#include <memory>
+#include <vector>
+
+#include "TFile.h"
+#include "TH2.h"
+#include "TString.h"
+#include "TTree.h"
+
+#include "LauSimpleFitModel.hh"
+#include "LauBkgndDPModel.hh"
+#include "LauDaughters.hh"
+#include "LauEffModel.hh"
+#include "LauIsobarDynamics.hh"
+#include "LauMagPhaseCoeffSet.hh"
+#include "LauResonanceMaker.hh"
+#include "LauVetoes.hh"
+
+int main( /*int argc, char** argv*/ )
+{
+ // If you want to use square DP histograms for efficiency,
+ // backgrounds or you just want the square DP co-ordinates
+ // stored in the toy MC ntuple then set this to kTRUE
+ Bool_t squareDP = kFALSE;
+
+ // This defines the DP => decay is B+ -> pi+ pi+ pi-
+ // Particle 1 = pi+
+ // Particle 2 = pi+
+ // Particle 3 = pi-
+ // The DP is defined in terms of m13Sq and m23Sq
+ LauDaughters* daughters = new LauDaughters("B+", "pi+", "pi+", "pi-", squareDP);
+
+ // Optionally apply some vetoes to the DP
+ // (example syntax given but commented-out)
+ LauVetoes* vetoes = new LauVetoes();
+ //Double_t DMin = 1.70;
+ //Double_t DMax = 1.925;
+ //Double_t JpsiMin = 3.051;
+ //Double_t JpsiMax = 3.222;
+ //Double_t psi2SMin = 3.676;
+ //Double_t psi2SMax = 3.866;
+ //vetoes->addMassVeto(1, DMin, DMax); // D0 veto, m23 (and automatically m13)
+ //vetoes->addMassVeto(1, JpsiMin, JpsiMax); // J/psi veto, m23 (and automatically m13)
+ //vetoes->addMassVeto(1, psi2SMin, psi2SMax); // psi(2S) veto, m23 (and automatically m13)
+
+ // Define the efficiency model (defaults to unity everywhere)
+ LauEffModel* effModel = new LauEffModel(daughters, vetoes);
+
+ // Can optionally provide a histogram to model variation over DP
+ // (example syntax given but commented-out)
+ //TFile *effHistFile = TFile::Open("histoFiles/B3piNRDPEff.root", "read");
+ //TH2* effHist = dynamic_cast<TH2*>(effHistFile->Get("effHist"));
+ //Bool_t useInterpolation = kTRUE;
+ //Bool_t fluctuateBins = kFALSE;
+ //Bool_t useUpperHalf = kTRUE;
+ //effModel->setEffHisto(effHist, useInterpolation, fluctuateBins, 0.0, 0.0, useUpperHalf, squareDP);
+
+ // Create the isobar model
+
+ // Set the values of the Blatt-Weisskopf barrier radii and whether they are fixed or floating
+ LauResonanceMaker& resMaker = LauResonanceMaker::get();
+ resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 5.0 );
+ resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 );
+ resMaker.fixBWRadius( LauBlattWeisskopfFactor::Parent, kTRUE );
+ resMaker.fixBWRadius( LauBlattWeisskopfFactor::Light, kTRUE );
+
+ LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel);
+ // Add various components to the isobar model,
+ // modifying some resonance parameters
+ LauAbsResonance* reson(0);
+ // addResonance arguments: resName, resPairAmpInt, resType
+ reson = sigModel->addResonance("rho0(770)", 1, LauAbsResonance::GS); // resPairAmpInt = 1 => resonance mass is m23.
+ reson = sigModel->addResonance("rho0(1450)", 1, LauAbsResonance::RelBW);
+ reson = sigModel->addResonance("f_0(980)", 1, LauAbsResonance::Flatte);
+ reson->setResonanceParameter("g1",0.2);
+ reson->setResonanceParameter("g2",1.0);
+ reson = sigModel->addResonance("f_2(1270)", 1, LauAbsResonance::RelBW);
+ const TString nrName = "BelleNR_Swave";
+ reson = sigModel->addResonance(nrName, 0, LauAbsResonance::BelleSymNRNoInter);
+ reson->setResonanceParameter("alpha", 0.2);
+
+ // Set the maximum signal DP ASq value
+ // If you do not provide a value, one will be determined automatically,
+ // which should be close to the true maximum but is not guaranteed to
+ // be optimal.
+ // Any value, whether manually provided or automatically determined,
+ // will be automatically adjusted to avoid bias or extreme inefficiency
+ // but it is best to set this by hand once you've found the right value
+ // through some trial and error.
+ sigModel->setASqMaxValue(0.35);
+
+ // Create the fit model
+ LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel);
+
+ // Create the complex coefficients for the isobar model
+ // Here we're using the magnitude and phase form:
+ // c_j = a_j exp(i*delta_j)
+ std::vector<LauAbsCoeffSet*> coeffset;
+ coeffset.push_back( new LauMagPhaseCoeffSet("rho0(770)", 1.00, 0.00, kTRUE, kTRUE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("rho0(1450)", 0.37, 1.99, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("f_0(980)", 0.27, -1.59, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet("f_2(1270)", 0.53, 1.39, kFALSE, kFALSE) );
+ coeffset.push_back( new LauMagPhaseCoeffSet(nrName, 0.54, -0.84, kFALSE, kFALSE) );
+ for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
+ fitModel->setAmpCoeffSet(*iter);
+ }
+
+ // Set the signal yield and define whether it is fixed or floated
+ const Double_t nSigEvents = 1500.0;
+ LauParameter * signalEvents = new LauParameter("signalEvents",nSigEvents,-1.0*nSigEvents,2.0*nSigEvents,kFALSE);
+ fitModel->setNSigEvents(signalEvents);
+
+ // Set up a background model
+
+ // First declare the names of the background class(es)
+ std::vector<TString> bkgndNames(1);
+ bkgndNames[0] = "qqbar";
+ fitModel->setBkgndClassNames( bkgndNames );
+
+ // Define and set the yield parameter for the background
+ const Double_t nBkg = 1250.0;
+ LauParameter* nBkgndEvents = new LauParameter("qqbar",nBkg,-2.0*nBkg,2.0*nBkg,kFALSE);
+ fitModel->setNBkgndEvents( nBkgndEvents );
+
+ // Create the background DP model
+ LauBkgndDPModel* qqbarModel = new LauBkgndDPModel(daughters, vetoes);
+
+ // Load in background DP model histogram
+ // (example syntax given but commented-out - the background will be treated as being uniform in the DP in the absence of a histogram)
+ //TString qqFileName("histoFiles/offResDP.root");
+ //TFile* qqFile = TFile::Open(qqFileName.Data(), "read");
+ //TH2* qqDP = dynamic_cast<TH2*>(qqFile->Get("AllmTheta")); // m', theta'
+ //qqbarModel->setBkgndHisto(qqDP, useInterpolation, fluctuateBins, useUpperHalf, squareDP);
+
+ // Add the background DP model into the fit model
+ fitModel->setBkgndDPModel( "qqbar", qqbarModel );
+
+
+ // Set the names of the files to read/write
+ const TString histFileName("plot-3pi.root");
+ std::unique_ptr<TFile> histFile { TFile::Open( histFileName, "recreate" ) };
+
+ const Int_t nBins{100};
+ const Double_t xMin{0.0};
+ const Double_t xMax{28.0};
+ const Double_t yMin{0.0};
+ const Double_t yMax{28.0};
+ TH2* signalLike = new TH2D( "signalLike", "Signal Likelihood", nBins, xMin, xMax, nBins, yMin, yMax );
+ TH2* qqbarLike = new TH2D( "qqbarLike", "qqbar Likelihood", nBins, xMin, xMax, nBins, yMin, yMax );
+
+ for ( Int_t iX{1}; iX < nBins; ++iX ) {
+
+ const Double_t m13Sq { signalLike->GetXaxis()->GetBinCenter( iX ) };
+
+ for ( Int_t iY{1}; iY < nBins; ++iY ) {
+ const Double_t m23Sq { signalLike->GetYaxis()->GetBinCenter( iY ) };
+
+ const auto likelihoods { fitModel->getDPLikelihoods( m13Sq, m23Sq ) };
+
+ signalLike->SetBinContent( iX, iY, likelihoods.at("signal") );
+ qqbarLike ->SetBinContent( iX, iY, likelihoods.at("qqbar") );
+ }
+ }
+
+ histFile->Write();
+ histFile->Close();
+
+ return EXIT_SUCCESS;
+}
diff --git a/inc/LauAbsBkgndDPModel.hh b/inc/LauAbsBkgndDPModel.hh
--- a/inc/LauAbsBkgndDPModel.hh
+++ b/inc/LauAbsBkgndDPModel.hh
@@ -57,6 +57,24 @@
*/
virtual Bool_t generate() = 0;
+ //! Get likelihood for a given DP position
+ /*!
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return the likelihood value
+ */
+ virtual Double_t getLikelihood( const Double_t m13Sq,
+ const Double_t m23Sq ) = 0;
+
+ //! Get unnormalised likelihood for a given DP position
+ /*!
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return the unnormalised likelihood value
+ */
+ virtual Double_t getUnNormValue( const Double_t m13Sq,
+ const Double_t m23Sq ) = 0;
+
//! Get likelihood for a given event
/*!
\param [in] iEvt the event number
diff --git a/inc/LauAbsFitModel.hh b/inc/LauAbsFitModel.hh
--- a/inc/LauAbsFitModel.hh
+++ b/inc/LauAbsFitModel.hh
@@ -71,6 +71,7 @@
#include <set>
#include <vector>
+#include "LauComplex.hh"
#include "LauFitObject.hh"
#include "LauFormulaPar.hh"
#include "LauSimFitTask.hh"
@@ -294,6 +295,24 @@
*/
void setParametersFileFallback(const TString& fileName, const TString& treeName, const std::map<TString, Double_t>& parameters, const Bool_t fix);
+ //! Calculate the DP amplitude(s) for a given DP position
+ /*!
+ If not already done, this function will initialise the DP models
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return a container of complex amplitudes, labelled to indicate to which component they belong
+ */
+ virtual std::map<TString, LauComplex> getDPAmps( const Double_t m13Sq, const Double_t m23Sq ) = 0;
+
+ //! Calculate the DP likelihood(s) for a given DP position
+ /*!
+ If not already done, this function will initialise the DP models
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return a container of likelihood values, labelled to indicate to which component they belong
+ */
+ virtual std::map<TString, Double_t> getDPLikelihoods( const Double_t m13Sq, const Double_t m23Sq ) = 0;
+
protected:
// Some typedefs
diff --git a/inc/LauBkgndDPModel.hh b/inc/LauBkgndDPModel.hh
--- a/inc/LauBkgndDPModel.hh
+++ b/inc/LauBkgndDPModel.hh
@@ -91,6 +91,31 @@
*/
virtual Bool_t generate();
+ //! Get likelihood for a given DP position
+ /*!
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return the likelihood value
+ */
+ virtual Double_t getLikelihood( const Double_t m13Sq,
+ const Double_t m23Sq );
+
+ //! Get unnormalised likelihood for a given DP position
+ /*!
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return the unnormalised likelihood value
+ */
+ virtual Double_t getUnNormValue( const Double_t m13Sq,
+ const Double_t m23Sq );
+
+ //! Get likelihood for a given event
+ /*!
+ \param [in] iEvt the event number
+ \return the likelihood value
+ */
+ virtual Double_t getLikelihood(UInt_t iEvt);
+
//! Get unnormalised likelihood for a given event
/*!
\param [in] iEvt the event number
@@ -104,13 +129,6 @@
*/
virtual Double_t getPdfNorm() const {return pdfNorm_;}
- //! Get likelihood for a given event
- /*!
- \param [in] iEvt the event number
- \return the likelihood value
- */
- virtual Double_t getLikelihood(UInt_t iEvt);
-
//! Cache the input data and (if appropriate) the per-event likelihood values
/*!
\param [in] fitDataTree the input data
diff --git a/inc/LauCPFitModel.hh b/inc/LauCPFitModel.hh
--- a/inc/LauCPFitModel.hh
+++ b/inc/LauCPFitModel.hh
@@ -224,6 +224,24 @@
*/
virtual void setAmpCoeffSet(LauAbsCoeffSet* coeffSet);
+ //! Calculate the DP amplitude(s) for a given DP position
+ /*!
+ If not already done, this function will initialise the DP models
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return a container of complex amplitudes, labelled to indicate to which component they belong
+ */
+ virtual std::map<TString, LauComplex> getDPAmps( const Double_t m13Sq, const Double_t m23Sq );
+
+ //! Calculate the DP likelihood(s) for a given DP position
+ /*!
+ If not already done, this function will initialise the DP models
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return a container of likelihood values, labelled to indicate to which component they belong
+ */
+ virtual std::map<TString, Double_t> getDPLikelihoods( const Double_t m13Sq, const Double_t m23Sq );
+
protected:
//! Define a map to be used to store a category name and numbers
typedef std::map< std::pair<TString,Int_t>, std::pair<Int_t,Double_t> > LauGenInfo;
diff --git a/inc/LauSimpleFitModel.hh b/inc/LauSimpleFitModel.hh
--- a/inc/LauSimpleFitModel.hh
+++ b/inc/LauSimpleFitModel.hh
@@ -173,6 +173,24 @@
*/
virtual void setAmpCoeffSet(LauAbsCoeffSet* coeffSet);
+ //! Calculate the DP amplitude(s) for a given DP position
+ /*!
+ If not already done, this function will initialise the DP models
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return a container of complex amplitudes, labelled to indicate to which component they belong
+ */
+ virtual std::map<TString, LauComplex> getDPAmps( const Double_t m13Sq, const Double_t m23Sq );
+
+ //! Calculate the DP likelihood(s) for a given DP position
+ /*!
+ If not already done, this function will initialise the DP models
+ \param [in] m13Sq the invariant mass squared of children 1 and 3
+ \param [in] m23Sq the invariant mass squared of children 2 and 3
+ \return a container of likelihood values, labelled to indicate to which component they belong
+ */
+ virtual std::map<TString, Double_t> getDPLikelihoods( const Double_t m13Sq, const Double_t m23Sq );
+
protected:
//! Define a map to be used to store a category name and numbers
typedef std::map< TString, std::pair<Int_t,Double_t> > LauGenInfo;
diff --git a/src/LauBkgndDPModel.cc b/src/LauBkgndDPModel.cc
--- a/src/LauBkgndDPModel.cc
+++ b/src/LauBkgndDPModel.cc
@@ -246,36 +246,45 @@
return;
}
- Double_t m13Sq(0.0), m23Sq(0.0), mPrime(0.0), thetaPrime(0.0);
-
- UInt_t nEvents = inputFitTree.nEvents();
- LauKinematics* kinematics = this->getKinematics();
+ const UInt_t nEvents { inputFitTree.nEvents() };
// clear and resize the data vector
bgData_.clear();
bgData_.resize(nEvents);
+ Double_t m13Sq { 0.0 }, m23Sq { 0.0 };
+
for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
const LauFitData& dataValues = inputFitTree.getData(iEvt);
- LauFitData::const_iterator iter = dataValues.find("m13Sq");
- m13Sq = iter->second;
- iter = dataValues.find("m23Sq");
- m23Sq = iter->second;
+ m13Sq = dataValues.at("m13Sq");
+ m23Sq = dataValues.at("m23Sq");
- // Update the kinematics. This will also update m' and theta' if squareDP = true
- kinematics->updateKinematics(m13Sq, m23Sq);
+ bgData_[iEvt] = this->getUnNormValue( m13Sq, m23Sq );
+ }
+}
- if (squareDP_ == kTRUE) {
- mPrime = kinematics->getmPrime();
- thetaPrime = kinematics->getThetaPrime();
- curEvtHistVal_ = this->calcHistValue(mPrime, thetaPrime);
- } else {
- curEvtHistVal_ = this->calcHistValue(m13Sq, m23Sq);
- }
+Double_t LauBkgndDPModel::getUnNormValue(const Double_t m13Sq, const Double_t m23Sq)
+{
+ // Update the kinematics. This will also update m' and theta' if squareDP = true
+ LauKinematics* kinematics = this->getKinematics();
+ kinematics->updateKinematics(m13Sq, m23Sq);
+
+ if (squareDP_) {
+ const Double_t mPrime { kinematics->getmPrime() };
+ const Double_t thetaPrime { kinematics->getThetaPrime() };
+ curEvtHistVal_ = this->calcHistValue(mPrime, thetaPrime);
+ } else {
+ curEvtHistVal_ = this->calcHistValue(m13Sq, m23Sq);
+ }
+
+ return curEvtHistVal_;
+}
- bgData_[iEvt] = curEvtHistVal_;
- }
+Double_t LauBkgndDPModel::getLikelihood(const Double_t m13Sq, const Double_t m23Sq)
+{
+ this->getUnNormValue( m13Sq, m23Sq );
+ return curEvtHistVal_ / this->getPdfNorm();
}
Double_t LauBkgndDPModel::getUnNormValue(UInt_t iEvt)
diff --git a/src/LauCPFitModel.cc b/src/LauCPFitModel.cc
--- a/src/LauCPFitModel.cc
+++ b/src/LauCPFitModel.cc
@@ -2915,6 +2915,80 @@
std::cout << "INFO in LauCPFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl;
}
+std::map<TString, LauComplex> LauCPFitModel::getDPAmps(const Double_t m13Sq, const Double_t m23Sq)
+{
+ // Initialise the DP model, if not already done
+ if ( negCoeffs_.empty() || posCoeffs_.empty() ) {
+ this->updateCoeffs();
+ this->initialiseDPModels();
+ }
+
+ if ( ! posKinematics_->withinDPLimits( m13Sq, m23Sq ) ) {
+ return { { "signal_" + negParent_, {} },
+ { "signal_" + posParent_, {} } };
+ }
+
+ negSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
+ posSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
+
+ return { { "signal_" + negParent_, negSigModel_->getEvtDPAmp() },
+ { "signal_" + posParent_, posSigModel_->getEvtDPAmp() } };
+}
+
+std::map<TString, Double_t> LauCPFitModel::getDPLikelihoods(const Double_t m13Sq, const Double_t m23Sq)
+{
+ // Initialise the DP model, if not already done
+ if ( negCoeffs_.empty() || posCoeffs_.empty() ) {
+ this->updateCoeffs();
+ this->initialiseDPModels();
+ }
+
+ std::map<TString, Double_t> likelihoods;
+
+ if ( ! posKinematics_->withinDPLimits( m13Sq, m23Sq ) ) {
+ likelihoods.emplace( "signal_" + negParent_, 0.0 );
+ likelihoods.emplace( "signal_" + posParent_, 0.0 );
+ if ( usingBkgnd_ ) {
+ const UInt_t nBkgnds { this->nBkgndClasses() };
+ for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
+ likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + negParent_, 0.0 );
+ likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + posParent_, 0.0 );
+ }
+ }
+ return likelihoods;
+ }
+
+ negSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
+ posSigModel_->calcLikelihoodInfo( m13Sq, m23Sq );
+
+ // See comment in getEvtDPLikelihood for explanation of the 2.0 factor
+ const Double_t norm { 2.0 / ( negSigModel_->getDPNorm() + posSigModel_->getDPNorm() ) };
+
+ likelihoods.emplace( "signal_" + negParent_, negSigModel_->getEvtIntensity() * norm );
+ likelihoods.emplace( "signal_" + posParent_, posSigModel_->getEvtIntensity() * norm );
+
+ // TODO - SCF signal
+ static bool warningIssued { false };
+ if ( useSCF_ && ! warningIssued ) {
+ warningIssued = true;
+ std::cerr << "WARNING in LauCPFitModel::getDPLikelihoods : calculation of SCF likelihoods not currently implemented in this function\n";
+ std::cerr << " : signal likelihood will just be the truth-matched value";
+ std::cerr << std::endl;
+ }
+
+ if ( usingBkgnd_ ) {
+ const UInt_t nBkgnds { this->nBkgndClasses() };
+ for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
+ likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + negParent_,
+ negBkgndDPModels_[bkgndID]->getLikelihood( m13Sq, m23Sq ) );
+ likelihoods.emplace( this->bkgndClassName( bkgndID ) + "_" + posParent_,
+ posBkgndDPModels_[bkgndID]->getLikelihood( m13Sq, m23Sq ) );
+ }
+ }
+
+ return likelihoods;
+}
+
void LauCPFitModel::embedNegSignal(const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
Bool_t useReweighting)
diff --git a/src/LauSimpleFitModel.cc b/src/LauSimpleFitModel.cc
--- a/src/LauSimpleFitModel.cc
+++ b/src/LauSimpleFitModel.cc
@@ -2015,6 +2015,67 @@
std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl;
}
+std::map<TString, LauComplex> LauSimpleFitModel::getDPAmps(const Double_t m13Sq, const Double_t m23Sq)
+{
+ // Initialise the DP model, if not already done
+ if ( coeffs_.empty() ) {
+ this->updateCoeffs();
+ this->initialiseDPModels();
+ }
+
+ if ( ! kinematics_->withinDPLimits( m13Sq, m23Sq ) ) {
+ return { { "signal", {} } };
+ }
+
+ sigDPModel_->calcLikelihoodInfo( m13Sq, m23Sq );
+
+ return { { "signal", sigDPModel_->getEvtDPAmp() } };
+}
+
+std::map<TString, Double_t> LauSimpleFitModel::getDPLikelihoods(const Double_t m13Sq, const Double_t m23Sq)
+{
+ // Initialise the DP model, if not already done
+ if ( coeffs_.empty() ) {
+ this->updateCoeffs();
+ this->initialiseDPModels();
+ }
+
+ std::map<TString, Double_t> likelihoods;
+
+ if ( ! kinematics_->withinDPLimits( m13Sq, m23Sq ) ) {
+ likelihoods.emplace( "signal", 0.0 );
+ if ( usingBkgnd_ ) {
+ const UInt_t nBkgnds { this->nBkgndClasses() };
+ for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
+ likelihoods.emplace( this->bkgndClassName( bkgndID ), 0.0 );
+ }
+ }
+ return likelihoods;
+ }
+
+ sigDPModel_->calcLikelihoodInfo( m13Sq, m23Sq );
+ likelihoods.emplace( "signal", sigDPModel_->getEvtLikelihood() );
+
+ // TODO - SCF signal
+ static bool warningIssued { false };
+ if ( useSCF_ && ! warningIssued ) {
+ warningIssued = true;
+ std::cerr << "WARNING in LauSimpleFitModel::getDPLikelihoods : calculation of SCF likelihoods not currently implemented in this function\n";
+ std::cerr << " : signal likelihood will just be the truth-matched value";
+ std::cerr << std::endl;
+ }
+
+ if ( usingBkgnd_ ) {
+ const UInt_t nBkgnds { this->nBkgndClasses() };
+ for ( UInt_t bkgndID( 0 ); bkgndID < nBkgnds; ++bkgndID ) {
+ likelihoods.emplace( this->bkgndClassName( bkgndID ),
+ bkgndDPModels_[bkgndID]->getLikelihood( m13Sq, m23Sq ) );
+ }
+ }
+
+ return likelihoods;
+}
+
void LauSimpleFitModel::embedSignal(const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
Bool_t useReweighting)

File Metadata

Mime Type
text/plain
Expires
Sun, Sep 28, 9:54 PM (10 h, 57 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6547348
Default Alt Text
D104.1759092863.diff (23 KB)

Event Timeline

HEPForge.org · .