Page MenuHomeHEPForge
Files F7879026

No OneTemporary
Actions

View Options

diff --git a/src/LauBkgndDPModel.cc b/src/LauBkgndDPModel.cc
index 9a89a77..fba3b84 100644
--- a/src/LauBkgndDPModel.cc
+++ b/src/LauBkgndDPModel.cc
@@ -1,288 +1,288 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauBkgndDPModel.cc
\brief File containing implementation of LauBkgndDPModel class.
*/
#include <cstdlib>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DHistDPPdf.hh"
#include "Lau2DSplineDPPdf.hh"
#include "LauBkgndDPModel.hh"
#include "LauDaughters.hh"
#include "LauFitDataTree.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
#include "LauVetoes.hh"
ClassImp(LauBkgndDPModel)
LauBkgndDPModel::LauBkgndDPModel(LauDaughters* daughters, LauVetoes* vetoes) :
LauAbsBkgndDPModel(daughters, vetoes),
symmetricalDP_(kFALSE),
squareDP_(kFALSE),
bgHistDPPdf_(0),
curEvtHistVal_(0.0),
maxPdfHeight_(1.0),
pdfNorm_(1.0),
doneGenWarning_(kFALSE),
lowBinWarningIssued_(kFALSE)
{
if (daughters != 0) {
symmetricalDP_ = daughters->gotSymmetricalDP();
}
}
LauBkgndDPModel::~LauBkgndDPModel()
{
if (bgHistDPPdf_ != 0) {
delete bgHistDPPdf_; bgHistDPPdf_ = 0;
}
}
void LauBkgndDPModel::setBkgndHisto(const TH2* histo, Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
{
Bool_t upperHalf = kFALSE;
if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
LauKinematics* kinematics = this->getKinematics();
const LauVetoes* vetoes = this->getVetoes();
bgHistDPPdf_ = new Lau2DHistDPPdf(histo, kinematics, vetoes, useInterpolation, fluctuateBins, upperHalf, squareDP);
maxPdfHeight_ = bgHistDPPdf_->getMaxHeight();
pdfNorm_ = bgHistDPPdf_->getHistNorm();
}
void LauBkgndDPModel::setBkgndSpline(const TH2* histo, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
{
Bool_t upperHalf = kFALSE;
if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
LauKinematics* kinematics = this->getKinematics();
const LauVetoes* vetoes = this->getVetoes();
bgHistDPPdf_ = new Lau2DSplineDPPdf(histo, kinematics, vetoes, fluctuateBins, upperHalf, squareDP);
maxPdfHeight_ = bgHistDPPdf_->getMaxHeight();
pdfNorm_ = bgHistDPPdf_->getHistNorm();
}
Double_t LauBkgndDPModel::calcHistValue(Double_t xVal, Double_t yVal)
{
// Get the likelihood value of the background in the Dalitz plot.
// Check that we have a valid histogram PDF
if (bgHistDPPdf_ == 0) {
cerr << "WARNING in LauBkgndDPModel::calcHistValue : We don't have a histogram so assuming the likelihood is flat in the Dalitz plot." << endl;
this->setBkgndHisto( 0, kFALSE, kFALSE, kFALSE, kFALSE );
}
// Find out the un-normalised PDF value
Double_t value = bgHistDPPdf_->interpolateXY(xVal, yVal);
// Check that the value is greater than zero
// If we're using a spline then negative values can be caused by adjacent bins that all contain a value of zero.
// The spline requires the value, its first derivatives and the mixed second derivative to be continuous and to match the input histogram
// at the bin centres. Derivatives are calculated using a finite difference approximation taking the difference between the neighbouring bins.
// If two bins are zero but the third is not then the second bin will have a positive first derivative causing the spline to dip below zero
// between the two zero bins to remain smooth. Such dips are unavoidable but are correctly removed here.
if ( value < 0.0 ) {
if(!lowBinWarningIssued_) {
std::cerr << "WARNING in LauBkgndDPModel::calcHistValue : Value " << value << " is less than 0 - setting to 0. You may want to check your histogram!" << std::endl
- << "If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
+ << " : If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
lowBinWarningIssued_=kTRUE;
}
return 0.0;
}
LauKinematics* kinematics = this->getKinematics();
// For square DP co-ordinates, need to divide by Jacobian
if (squareDP_ == kTRUE) {
// Make sure that square DP kinematics are correct, then get Jacobian
kinematics->updateSqDPKinematics(xVal, yVal);
Double_t jacobian = kinematics->calcSqDPJacobian();
value /= jacobian;
}
return value;
}
void LauBkgndDPModel::initialise()
{
}
Bool_t LauBkgndDPModel::generate()
{
// Routine to generate the background, using data provided by an
// already defined histogram.
LauKinematics* kinematics = this->getKinematics();
Bool_t gotBG(kFALSE);
while (gotBG == kFALSE) {
if (squareDP_ == kTRUE) {
// Generate a point in m', theta' space. By construction, this point
// is already within the DP region.
Double_t mPrime(0.0), thetaPrime(0.0);
kinematics->genFlatSqDP(mPrime, thetaPrime);
// If we're in a symmetrical DP then we should only generate events in one half
if ( symmetricalDP_ && thetaPrime > 0.5 ) {
thetaPrime = 1.0 - thetaPrime;
}
// Calculate histogram height for DP point and
// compare with the maximum height
if ( bgHistDPPdf_ != 0 ) {
Double_t bgContDP = bgHistDPPdf_->interpolateXY(mPrime, thetaPrime)/maxPdfHeight_;
if (LauRandom::randomFun()->Rndm() < bgContDP) {
kinematics->updateSqDPKinematics(mPrime, thetaPrime);
gotBG = kTRUE;
}
} else {
if ( !doneGenWarning_ ) {
cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the square DP, which won't be flat in the conventional DP!" << endl;
cerr << "WARNING in LauBkgndDPModel::generate : This should never happen!! What have you done?!" << endl;
doneGenWarning_ = kTRUE;
}
kinematics->updateSqDPKinematics(mPrime, thetaPrime);
gotBG = kTRUE;
}
} else {
// Generate a point in the Dalitz plot (phase-space).
Double_t m13Sq(0.0), m23Sq(0.0);
kinematics->genFlatPhaseSpace(m13Sq, m23Sq);
// If we're in a symmetrical DP then we should only generate events in one half
if ( symmetricalDP_ && m13Sq > m23Sq ) {
Double_t tmpSq = m13Sq;
m13Sq = m23Sq;
m23Sq = tmpSq;
}
// Calculate histogram height for DP point and
// compare with the maximum height
if ( bgHistDPPdf_ != 0 ) {
Double_t bgContDP = bgHistDPPdf_->interpolateXY(m13Sq, m23Sq)/maxPdfHeight_;
if (LauRandom::randomFun()->Rndm() < bgContDP) {
kinematics->updateKinematics(m13Sq, m23Sq);
gotBG = kTRUE;
}
} else {
if ( !doneGenWarning_ ) {
cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the DP." << endl;
doneGenWarning_ = kTRUE;
}
kinematics->updateKinematics(m13Sq, m23Sq);
gotBG = kTRUE;
}
}
}
// Implement veto
Bool_t vetoOK(kTRUE);
const LauVetoes* vetoes = this->getVetoes();
if (vetoes) {
vetoOK = vetoes->passVeto(kinematics);
}
// Call this function recusively until we pass the veto.
if (vetoOK == kFALSE) {
this->generate();
}
return kTRUE;
}
void LauBkgndDPModel::fillDataTree(const LauFitDataTree& inputFitTree)
{
// In LauFitDataTree, the first two variables should always be
// m13^2 and m23^2. Other variables follow thus: charge.
Int_t nBranches = inputFitTree.nBranches();
if (nBranches < 2) {
cout<<"Error in LauBkgndDPModel::fillDataTree. Expecting at least 2 variables "
<<"in input data tree, but there are "<<nBranches<<"! Make sure you have "
<<"the right number of variables in your input data file!"<<endl;
return;
}
Double_t m13Sq(0.0), m23Sq(0.0), mPrime(0.0), thetaPrime(0.0);
UInt_t nEvents = inputFitTree.nEvents();
LauKinematics* kinematics = this->getKinematics();
// clear and resize the data vector
bgData_.clear();
bgData_.resize(nEvents);
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;
// Update the kinematics. This will also update m' and theta' if squareDP = true
kinematics->updateKinematics(m13Sq, m23Sq);
if (squareDP_ == kTRUE) {
mPrime = kinematics->getmPrime();
thetaPrime = kinematics->getThetaPrime();
curEvtHistVal_ = this->calcHistValue(mPrime, thetaPrime);
} else {
curEvtHistVal_ = this->calcHistValue(m13Sq, m23Sq);
}
bgData_[iEvt] = curEvtHistVal_;
}
}
Double_t LauBkgndDPModel::getUnNormValue(UInt_t iEvt)
{
// Retrieve the likelihood for the given event
this->setDataEventNo(iEvt);
return curEvtHistVal_;
}
Double_t LauBkgndDPModel::getLikelihood(UInt_t iEvt)
{
// Retrieve the likelihood for the given event
this->setDataEventNo(iEvt);
Double_t llhd = curEvtHistVal_ / this->getPdfNorm();
return llhd;
}
void LauBkgndDPModel::setDataEventNo(UInt_t iEvt)
{
// Retrieve the data for event iEvt
if (bgData_.size() > iEvt) {
curEvtHistVal_ = bgData_[iEvt];
} else {
cerr<<"ERROR in LauBkgndDPModel::setDataEventNo : Event index too large: "<<iEvt<<" >= "<<bgData_.size()<<"."<<endl;
}
}
diff --git a/src/LauEffModel.cc b/src/LauEffModel.cc
index cff25dc..464d101 100644
--- a/src/LauEffModel.cc
+++ b/src/LauEffModel.cc
@@ -1,191 +1,191 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauEffModel.cc
\brief File containing implementation of LauEffModel class.
*/
#include <cstdlib>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include "TSystem.h"
#include "Lau2DHistDP.hh"
#include "Lau2DSplineDP.hh"
#include "LauDaughters.hh"
#include "LauEffModel.hh"
#include "LauKinematics.hh"
#include "LauVetoes.hh"
ClassImp(LauEffModel)
LauEffModel::LauEffModel(const LauDaughters* daughters, const LauVetoes* vetoes) :
daughters_( daughters ),
vetoes_( vetoes ),
effHisto_( 0 ),
squareDP_( kFALSE ),
fluctuateEffHisto_( kFALSE ),
lowBinWarningIssued_( kFALSE ),
highBinWarningIssued_( kFALSE )
{
if ( daughters_ == 0 ) {
cerr << "ERROR in LauEffModel Constructor : invalid pointer to daughters object supplied." << endl;
gSystem->Exit(EXIT_FAILURE);
}
}
/*LauEffModel::LauEffModel( const LauEffModel& rhs ) :
daughters_( rhs.daughters_ ),
vetoes_( rhs.vetoes_ ),
effHisto_( rhs.effHisto_ ? new Lau2DHistDP( *rhs.effHisto_ ) : 0 ),
squareDP_( rhs.squareDP_ ),
fluctuateEffHisto_( rhs.fluctuateEffHisto_ )
{
}*/
LauEffModel::~LauEffModel()
{
if (effHisto_ != 0) {
delete effHisto_; effHisto_ = 0;
}
}
void LauEffModel::setEffHisto(const TH2* effHisto, Bool_t useInterpolation,
Bool_t fluctuateBins, Double_t avEff, Double_t absError,
Bool_t useUpperHalfOnly, Bool_t squareDP)
{
// Set the efficiency across the Dalitz plot using a predetermined 2D histogram
// with x = m_13^2, y = m_23^2.
Bool_t upperHalf( kFALSE );
if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
if(effHisto_) {
delete effHisto_;
effHisto_=0;
}
// Copy the histogram.
effHisto_ = new Lau2DHistDP(effHisto, daughters_,
useInterpolation, fluctuateBins,
avEff, absError, upperHalf, squareDP);
fluctuateEffHisto_ = fluctuateBins;
if (avEff > 0.0 && absError > 0.0) {
fluctuateEffHisto_ = kTRUE;
}
}
void LauEffModel::setEffSpline(const TH2* effHisto,
Bool_t fluctuateBins, Double_t avEff, Double_t absError,
Bool_t useUpperHalfOnly, Bool_t squareDP)
{
// Set the efficiency across the Dalitz plot using a predetermined 2D histogram
// with x = m_13^2, y = m_23^2.
Bool_t upperHalf( kFALSE );
if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
if(effHisto_) {
delete effHisto_;
effHisto_=0;
}
// Copy the histogram.
effHisto_ = new Lau2DSplineDP(effHisto, daughters_,
fluctuateBins, avEff, absError, upperHalf, squareDP);
fluctuateEffHisto_ = fluctuateBins;
if (avEff > 0.0 && absError > 0.0) {
fluctuateEffHisto_ = kTRUE;
}
}
Double_t LauEffModel::getEffHistValue(Double_t xVal, Double_t yVal) const
{
// Get the efficiency from the 2D histo.
Double_t eff(1.0);
if (effHisto_ != 0) {
eff = effHisto_->interpolateXY(xVal, yVal);
}
return eff;
}
Double_t LauEffModel::calcEfficiency( const LauKinematics* kinematics ) const
{
// Routine to calculate the efficiency for the given event/point in
// the Dalitz plot. This routine uses the 2D histogram set by the
// setEffHisto() funciton.
Double_t eff(1.0);
// Check for vetoes
Bool_t vetoOK(kTRUE);
if ( vetoes_ != 0 ) {
vetoOK = vetoes_->passVeto(kinematics);
}
if (vetoOK == kFALSE) {
// We failed the veto.
eff = 0.0;
} else {
// We are using a 2D histogram representation of the efficiency across the Dalitz plot.
// First find out which bin we are in given the x and y Dalitz plot mass values
// Here, x = m13^2, y = m23^2 or if using square DP x = m', y = theta'.
if (squareDP_ == kTRUE) {
eff = this->getEffHistValue(kinematics->getmPrime(), kinematics->getThetaPrime());
} else {
eff = this->getEffHistValue(kinematics->getm13Sq(), kinematics->getm23Sq());
}
}
// Check that the efficiency is in the allowed range (0-1)
// If we're using a spline then out-of-range efficiencies can be caused by adjacent bins that all contain a value of either zero or one.
// The spline requires the efficiency, its first derivatives and the mixed second derivative to be continuous and to match the input histogram
// at the bin centres. Derivatives are calculated using a finite difference approximation taking the difference between the neighbouring bins.
// If two bins are zero but the third is not then the second bin will have a positive first derivative causing the spline to dip below zero
// between the two zero bins to remain smooth. The analogous case with adjacent maximised bins will cause peaks above one. Such dips are
// unavoidable but are correctly removed here.
if ( eff < 0.0 ) {
if(!lowBinWarningIssued_) {
std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is less than 0 - setting to 0. You may want to check your histogram!" << std::endl
- << "If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
+ << " : If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
lowBinWarningIssued_=kTRUE;
}
eff = 0.0;
} else if ( eff > 1.0 ) {
if(!highBinWarningIssued_) {
std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is greater than 1 - setting to 1. You may want to check your histogram!" << std::endl
- << "If you are using a spline then this could be caused by adjacent full bins. Further warnings will be suppressed." << std::endl;
+ << " : If you are using a spline then this could be caused by adjacent full bins. Further warnings will be suppressed." << std::endl;
highBinWarningIssued_=kTRUE;
}
eff = 1.0;
}
return eff;
}
Bool_t LauEffModel::passVeto( const LauKinematics* kinematics ) const
{
Bool_t pass = kTRUE;
if ( vetoes_ != 0 ) {
pass = vetoes_->passVeto( kinematics );
}
return pass;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:18 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805799
Default Alt Text
(16 KB)

Event Timeline