No OneTemporary
Actions

Size

38 KB

Subscribers

None

View Options

	diff --git a/EvtGenModels/EvtDalitzTable.hh b/EvtGenModels/EvtDalitzTable.hh
	index acb95a7..8dc3a4d 100644
	--- a/EvtGenModels/EvtDalitzTable.hh
	+++ b/EvtGenModels/EvtDalitzTable.hh
	@@ -1,80 +1,80 @@

	/***********************************************************************
	* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
	* *
	* This file is part of EvtGen. *
	* *
	* EvtGen is free software: you can redistribute it and/or modify *
	* it under the terms of the GNU General Public License as published by *
	* the Free Software Foundation, either version 3 of the License, or *
	* (at your option) any later version. *
	* *
	* EvtGen is distributed in the hope that it will be useful, *
	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	* GNU General Public License for more details. *
	* *
	* You should have received a copy of the GNU General Public License *
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/

	#ifndef EVTDALITZTABLE_HPP
	#define EVTDALITZTABLE_HPP

	#include "EvtGenBase/EvtCyclic3.hh"
	#include "EvtGenBase/EvtDalitzPlot.hh"
	#include "EvtGenBase/EvtDalitzReso.hh"
	#include "EvtGenBase/EvtId.hh"
	#include "EvtGenBase/EvtSpinType.hh"

	#include "EvtGenModels/EvtDalitzDecayInfo.hh"

	#include <map>
	#include <string>
	#include <vector>

	// Description: Model to describe a generic dalitz decay

	class EvtDalitzTable {
	public:
	- static EvtDalitzTable* getInstance( const std::string dec_name = "",
	- bool verbose = true );
	+ static const EvtDalitzTable& getInstance( const std::string dec_name = "",
	+ bool verbose = true );

	+ std::vector<EvtDalitzDecayInfo> getDalitzTable( const EvtId& parent ) const;
	+
	+ protected:
	+ EvtDalitzTable();
	+ ~EvtDalitzTable();
	+
	+ private:
	bool fileHasBeenRead( const std::string dec_name ) const;
	void readXMLDecayFile( const std::string dec_name, bool verbose = true );
	void checkParticle( std::string particle ) const;

	void addDecay( EvtId parent, const EvtDalitzDecayInfo& dec );
	void copyDecay( EvtId parent, EvtId* daughters, EvtId copy, EvtId* copyd );

	- std::vector<EvtDalitzDecayInfo> getDalitzTable( const EvtId& parent );
	-
	- protected:
	- EvtDalitzTable();
	- ~EvtDalitzTable();
	-
	- private:
	EvtDalitzReso getResonance( std::string shape, EvtDalitzPlot dp,
	EvtCyclic3::Pair angPair,
	EvtCyclic3::Pair resPair,
	EvtSpinType::spintype spinType, double mass,
	double width, double FFp, double FFr,
	double alpha, double aLass, double rLass,
	double BLass, double phiBLass, double RLass,
	double phiRLass, double cutoffLass );
	int getDaughterPairs(
	EvtId* resDaughter, EvtId* daughter,
	std::vector<std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair>>& angAndResPairs );

	std::map<EvtId, std::vector<EvtDalitzDecayInfo>> m_dalitztable;
	std::vector<std::string> m_readFiles;

	EvtDalitzTable( const EvtDalitzTable& );
	EvtDalitzTable& operator=( const EvtDalitzTable& );

	//to calculate probMax
	double calcProbMax( EvtDalitzPlot dp, EvtDalitzDecayInfo* model );
	double calcProb( EvtDalitzPoint point, EvtDalitzDecayInfo* model );
	};

	#endif
	diff --git a/src/EvtGenModels/EvtDalitzTable.cpp b/src/EvtGenModels/EvtDalitzTable.cpp
	index bf447d0..99c5943 100644
	--- a/src/EvtGenModels/EvtDalitzTable.cpp
	+++ b/src/EvtGenModels/EvtDalitzTable.cpp
	@@ -1,673 +1,672 @@

	/***********************************************************************
	* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
	* *
	* This file is part of EvtGen. *
	* *
	* EvtGen is free software: you can redistribute it and/or modify *
	* it under the terms of the GNU General Public License as published by *
	* the Free Software Foundation, either version 3 of the License, or *
	* (at your option) any later version. *
	* *
	* EvtGen is distributed in the hope that it will be useful, *
	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	* GNU General Public License for more details. *
	* *
	* You should have received a copy of the GNU General Public License *
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/

	#include "EvtGenModels/EvtDalitzTable.hh"

	#include "EvtGenBase/EvtCyclic3.hh"
	#include "EvtGenBase/EvtDalitzPlot.hh"
	#include "EvtGenBase/EvtPDL.hh"
	#include "EvtGenBase/EvtParserXml.hh"
	#include "EvtGenBase/EvtReport.hh"
	#include "EvtGenBase/EvtSpinType.hh"

	#include <sstream>
	#include <stdlib.h>

	using std::endl;
	using std::fstream;
	using std::ifstream;

	EvtDalitzTable::EvtDalitzTable()
	{
	m_dalitztable.clear();
	m_readFiles.clear();
	}

	EvtDalitzTable::~EvtDalitzTable()
	{
	m_dalitztable.clear();
	m_readFiles.clear();
	}

	-EvtDalitzTable* EvtDalitzTable::getInstance( const std::string dec_name,
	- bool verbose )
	+const EvtDalitzTable& EvtDalitzTable::getInstance( const std::string dec_name,
	+ bool verbose )
	{
	- static thread_local EvtDalitzTable* theDalitzTable = nullptr;
	+ static thread_local EvtDalitzTable theDalitzTable;

	- if ( theDalitzTable == nullptr ) {
	- theDalitzTable = new EvtDalitzTable();
	- }
	-
	- if ( !theDalitzTable->fileHasBeenRead( dec_name ) ) {
	- theDalitzTable->readXMLDecayFile( dec_name, verbose );
	+ if ( !theDalitzTable.fileHasBeenRead( dec_name ) ) {
	+ theDalitzTable.readXMLDecayFile( dec_name, verbose );
	}

	return theDalitzTable;
	}

	bool EvtDalitzTable::fileHasBeenRead( const std::string dec_name ) const
	{
	std::vector<std::string>::const_iterator i = m_readFiles.begin();
	for ( ; i != m_readFiles.end(); i++ ) {
	if ( ( *i ).compare( dec_name ) == 0 ) {
	return true;
	}
	}
	return false;
	}

	void EvtDalitzTable::readXMLDecayFile( const std::string dec_name, bool verbose )
	{
	if ( verbose ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "EvtDalitzTable: Reading in xml parameter file " << dec_name
	<< endl;
	}

	m_readFiles.push_back( dec_name );

	EvtDalitzDecayInfo* dalitzDecay = nullptr;
	double probMax = 0;
	EvtId ipar;
	std::string decayParent = "";
	std::string daugStr = "";
	EvtId daughter[3];

	EvtDalitzPlot dp;
	EvtComplex cAmp;
	std::vector<std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair>> angAndResPairs;
	std::string shape( "" );
	EvtSpinType::spintype spinType( EvtSpinType::SCALAR );
	double mass( 0. ), width( 0. ), FFp( 0. ), FFr( 0. );
	std::vector<EvtFlatteParam> flatteParams;
	//Nonres parameters
	double alpha( 0. );
	//LASS parameters
	double aLass( 0. ), rLass( 0. ), BLass( 0. ), phiBLass( 0. ), RLass( 0. ),
	phiRLass( 0. ), cutoffLass( -1. );

	EvtParserXml parser;
	parser.open( dec_name );

	bool endReached = false;

	while ( parser.readNextTag() ) {
	//TAGS FOUND UNDER DATA
	if ( parser.getParentTagTitle() == "data" ) {
	if ( parser.getTagTitle() == "dalitzDecay" ) {
	int nDaughters = 0;

	decayParent = parser.readAttribute( "particle" );
	daugStr = parser.readAttribute( "daughters" );
	probMax = parser.readAttributeDouble( "probMax", -1 );

	checkParticle( decayParent );
	ipar = EvtPDL::getId( decayParent );

	std::istringstream daugStream( daugStr );

	std::string daugh;
	while ( std::getline( daugStream, daugh, ' ' ) ) {
	checkParticle( daugh );
	daughter[nDaughters++] = EvtPDL::getId( daugh );
	}

	if ( nDaughters != 3 ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Expected to find three daughters for dalitzDecay of "
	<< decayParent << " near line "
	<< parser.getLineNumber() << ", "
	<< "found " << nDaughters << endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Will terminate execution!" << endl;
	::abort();
	}

	double m_d1 = EvtPDL::getMass( daughter[0] ),
	m_d2 = EvtPDL::getMass( daughter[1] ),
	m_d3 = EvtPDL::getMass( daughter[2] ),
	M = EvtPDL::getMass( ipar );

	dp = EvtDalitzPlot( m_d1, m_d2, m_d3, M );

	dalitzDecay = new EvtDalitzDecayInfo( daughter[0], daughter[1],
	daughter[2] );

	} else if ( parser.getTagTitle() == "copyDalitz" ) {
	int nDaughters = 0;
	int nCopyDaughters = 0;
	EvtId copyDaughter[3];

	decayParent = parser.readAttribute( "particle" );
	daugStr = parser.readAttribute( "daughters" );

	std::string copyParent = parser.readAttribute( "copy" );
	std::string copyDaugStr = parser.readAttribute( "copyDaughters" );

	checkParticle( decayParent );
	ipar = EvtPDL::getId( decayParent );

	checkParticle( copyParent );
	EvtId copypar = EvtPDL::getId( copyParent );

	std::istringstream daugStream( daugStr );
	std::istringstream copyDaugStream( copyDaugStr );

	std::string daugh;
	while ( std::getline( daugStream, daugh, ' ' ) ) {
	checkParticle( daugh );
	daughter[nDaughters++] = EvtPDL::getId( daugh );
	}
	while ( std::getline( copyDaugStream, daugh, ' ' ) ) {
	checkParticle( daugh );
	copyDaughter[nCopyDaughters++] = EvtPDL::getId( daugh );
	}

	if ( nDaughters != 3 \|\| nCopyDaughters != 3 ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Expected to find three daughters for copyDecay of "
	<< decayParent << " from " << copyParent
	<< " near line " << parser.getLineNumber() << ", "
	<< "found " << nDaughters << " and " << nCopyDaughters
	<< endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Will terminate execution!" << endl;
	::abort();
	}

	copyDecay( ipar, daughter, copypar, copyDaughter );

	} else if ( parser.getTagTitle() == "/data" ) { //end of data
	endReached = true;
	parser.close();
	break;
	}
	//TAGS FOUND UNDER DALITZDECAY
	} else if ( parser.getParentTagTitle() == "dalitzDecay" ) {
	if ( parser.getTagTitle() == "resonance" ) {
	flatteParams.clear();

	//Amplitude
	EvtComplex ampFactor(
	parser.readAttributeDouble( "ampFactorReal", 1. ),
	parser.readAttributeDouble( "ampFactorImag", 0. ) );
	double mag = parser.readAttributeDouble( "mag", -999. );
	double phase = parser.readAttributeDouble( "phase", -999. );
	double real = parser.readAttributeDouble( "real", -999. );
	double imag = parser.readAttributeDouble( "imag", -999. );

	if ( ( real != -999. \|\| imag != -999. ) && mag == -999. &&
	phase == -999. ) {
	if ( real == -999. ) {
	real = 0;
	}
	if ( imag == -999. ) {
	imag = 0;
	}
	mag = sqrt( real * real + imag * imag );
	phase = atan2( imag, real ) * EvtConst::radToDegrees;
	}
	if ( mag == -999. ) {
	mag = 1.;
	}
	if ( phase == -999. ) {
	phase = 0.;
	}

	cAmp = ampFactor *
	EvtComplex( cos( phase * 1.0 / EvtConst::radToDegrees ),
	sin( phase * 1.0 / EvtConst::radToDegrees ) ) *
	mag;

	//Resonance particle properties
	mass = 0.;
	width = 0.;
	spinType = EvtSpinType::SCALAR;

	std::string particle = parser.readAttribute( "particle" );
	if ( particle != "" ) {
	EvtId resId = EvtPDL::getId( particle );
	if ( resId == EvtId( -1, -1 ) ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Unknown particle name:" << particle.c_str()
	<< endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Will terminate execution!" << endl;
	::abort();
	} else {
	mass = EvtPDL::getMeanMass( resId );
	width = EvtPDL::getWidth( resId );
	spinType = EvtPDL::getSpinType( resId );
	}
	}

	width = parser.readAttributeDouble( "width", width );
	mass = parser.readAttributeDouble( "mass", mass );
	switch ( parser.readAttributeInt( "spin", -1 ) ) {
	case -1: //not set here
	break;
	case 0:
	spinType = EvtSpinType::SCALAR;
	break;
	case 1:
	spinType = EvtSpinType::VECTOR;
	break;
	case 2:
	spinType = EvtSpinType::TENSOR;
	break;
	default:
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Unsupported spin near line "
	<< parser.getLineNumber() << " of XML file." << endl;
	::abort();
	}

	//Shape and form factors
	shape = parser.readAttribute( "shape" );
	FFp = parser.readAttributeDouble( "BlattWeisskopfFactorParent",
	0.0 );
	FFr = parser.readAttributeDouble( "BlattWeisskopfFactorResonance",
	1.5 );

	//Shape specific attributes
	if ( shape == "NonRes_Exp" ) {
	alpha = parser.readAttributeDouble( "alpha", 0.0 );
	}
	if ( shape == "LASS" ) {
	aLass = parser.readAttributeDouble( "a", 2.07 );
	rLass = parser.readAttributeDouble( "r", 3.32 );
	BLass = parser.readAttributeDouble( "B", 1.0 );
	phiBLass = parser.readAttributeDouble( "phiB", 0.0 );
	RLass = parser.readAttributeDouble( "R", 1.0 );
	phiRLass = parser.readAttributeDouble( "phiR", 0.0 );
	cutoffLass = parser.readAttributeDouble( "cutoff", -1.0 );
	}

	//Daughter pairs for resonance
	angAndResPairs.clear();

	std::string resDaugStr = parser.readAttribute( "resDaughters" );
	if ( resDaugStr != "" ) {
	std::istringstream resDaugStream( resDaugStr );
	std::string resDaug;
	int nResDaug( 0 );
	EvtId resDaughter[2];
	while ( std::getline( resDaugStream, resDaug, ' ' ) ) {
	checkParticle( resDaug );
	resDaughter[nResDaug++] = EvtPDL::getId( resDaug );
	}
	if ( nResDaug != 2 ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Resonance must have exactly 2 daughters near line "
	<< parser.getLineNumber() << " of XML file." << endl;
	::abort();
	}
	int nRes = getDaughterPairs( resDaughter, daughter,
	angAndResPairs );
	if ( nRes == 0 ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Resonance daughters must match decay daughters near line "
	<< parser.getLineNumber() << " of XML file." << endl;
	::abort();
	}
	if ( parser.readAttributeBool( "normalise", true ) )
	cAmp /= sqrt( nRes );
	}
	if ( angAndResPairs.empty() ) {
	switch ( parser.readAttributeInt( "daughterPair" ) ) {
	case 1:
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::BC, EvtCyclic3::AB ) );
	break;
	case 2:
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::CA, EvtCyclic3::BC ) );
	break;
	case 3:
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::AB, EvtCyclic3::CA ) );
	break;
	default:
	if ( shape ==
	"NonRes" ) { //We don't expect a pair for non-resonant terms but add dummy values for convenience
	angAndResPairs.push_back( std::make_pair(
	EvtCyclic3::AB, EvtCyclic3::AB ) );
	break;
	}
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Daughter pair must be 1, 2 or 3 near line "
	<< parser.getLineNumber() << " of XML file."
	<< endl;
	::abort();
	}
	}

	if ( parser.isTagInline() ) {
	std::vector<std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair>>::iterator
	it = angAndResPairs.begin();
	for ( ; it != angAndResPairs.end(); it++ ) {
	std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair> pairs = *it;
	EvtDalitzReso resonance = getResonance(
	shape, dp, pairs.first, pairs.second, spinType,
	mass, width, FFp, FFr, alpha, aLass, rLass, BLass,
	phiBLass, RLass, phiRLass, cutoffLass );
	dalitzDecay->addResonance( cAmp, resonance );
	}
	}
	} else if ( parser.getTagTitle() == "/dalitzDecay" ) {
	if ( probMax < 0 ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "probMax is not defined for " << decayParent
	<< " -> " << daugStr << endl;
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Will now estimate probMax. This may take a while. Once probMax is calculated, update the XML file to skip this step in future."
	<< endl;
	probMax = calcProbMax( dp, dalitzDecay );
	}
	dalitzDecay->setProbMax( probMax );
	addDecay( ipar, *dalitzDecay );
	delete dalitzDecay;
	dalitzDecay = nullptr;
	} else if ( verbose ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Unexpected tag " << parser.getTagTitle()
	<< " found in XML decay file near line "
	<< parser.getLineNumber() << ". Tag will be ignored."
	<< endl;
	}
	//TAGS FOUND UNDER RESONANCE
	} else if ( parser.getParentTagTitle() == "resonance" ) {
	if ( parser.getTagTitle() == "flatteParam" ) {
	EvtFlatteParam param( parser.readAttributeDouble( "mass1" ),
	parser.readAttributeDouble( "mass2" ),
	parser.readAttributeDouble( "g" ) );
	flatteParams.push_back( param );
	} else if ( parser.getTagTitle() == "/resonance" ) {
	std::vector<std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair>>::iterator
	it = angAndResPairs.begin();
	for ( ; it != angAndResPairs.end(); it++ ) {
	std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair> pairs = *it;
	EvtDalitzReso resonance = getResonance(
	shape, dp, pairs.first, pairs.second, spinType, mass,
	width, FFp, FFr, alpha, aLass, rLass, BLass, phiBLass,
	RLass, phiRLass, cutoffLass );

	std::vector<EvtFlatteParam>::iterator flatteIt =
	flatteParams.begin();
	for ( ; flatteIt != flatteParams.end(); flatteIt++ ) {
	resonance.addFlatteParam( ( *flatteIt ) );
	}

	dalitzDecay->addResonance( cAmp, resonance );
	}
	}
	}
	}

	if ( !endReached ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Either the decay file ended prematurely or the file is badly formed.\n"
	<< "Error occured near line" << parser.getLineNumber() << endl;
	::abort();
	}
	}

	void EvtDalitzTable::checkParticle( std::string particle ) const
	{
	if ( EvtPDL::getId( particle ) == EvtId( -1, -1 ) ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Unknown particle name:" << particle.c_str() << endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Will terminate execution!" << endl;
	::abort();
	}
	}

	void EvtDalitzTable::addDecay( EvtId parent, const EvtDalitzDecayInfo& dec )
	{
	if ( m_dalitztable.find( parent ) != m_dalitztable.end() ) {
	m_dalitztable[parent].push_back( dec );
	} else {
	m_dalitztable[parent].push_back( dec );
	}
	}

	void EvtDalitzTable::copyDecay( EvtId parent, EvtId* daughters, EvtId copy,
	EvtId* copyd )
	{
	EvtDalitzDecayInfo decay( daughters[0], daughters[1], daughters[2] );
	std::vector<EvtDalitzDecayInfo> copyTable = getDalitzTable( copy );
	std::vector<EvtDalitzDecayInfo>::iterator i = copyTable.begin();
	for ( ; i != copyTable.end(); i++ ) {
	EvtId daughter1 = ( *i ).daughter1();
	EvtId daughter2 = ( *i ).daughter2();
	EvtId daughter3 = ( *i ).daughter3();

	if ( ( copyd[0] == daughter1 && copyd[1] == daughter2 &&
	copyd[2] == daughter3 ) \|\|
	( copyd[0] == daughter1 && copyd[1] == daughter3 &&
	copyd[2] == daughter2 ) \|\|
	( copyd[0] == daughter2 && copyd[1] == daughter1 &&
	copyd[2] == daughter3 ) \|\|
	( copyd[0] == daughter2 && copyd[1] == daughter3 &&
	copyd[2] == daughter1 ) \|\|
	( copyd[0] == daughter3 && copyd[1] == daughter1 &&
	copyd[2] == daughter2 ) \|\|
	( copyd[0] == daughter3 && copyd[1] == daughter2 &&
	copyd[2] == daughter1 ) ) {
	decay.setProbMax( ( *i ).getProbMax() );
	std::vector<std::pair<EvtComplex, EvtDalitzReso>>::const_iterator j =
	( *i ).getResonances().begin();
	for ( ; j != ( *i ).getResonances().end(); j++ ) {
	decay.addResonance( ( *j ) );
	}
	addDecay( parent, decay );
	return;
	}
	}
	//if we get here then there was no match
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Did not find dalitz decays for particle:" << copy << "\n";
	}

	-std::vector<EvtDalitzDecayInfo> EvtDalitzTable::getDalitzTable( const EvtId& parent )
	+std::vector<EvtDalitzDecayInfo> EvtDalitzTable::getDalitzTable(
	+ const EvtId& parent ) const
	{
	std::vector<EvtDalitzDecayInfo> table;
	- if ( m_dalitztable.find( parent ) != m_dalitztable.end() ) {
	- table = m_dalitztable[parent];
	+
	+ auto iter = m_dalitztable.find( parent );
	+ if ( iter != m_dalitztable.end() ) {
	+ table = iter->second;
	}

	if ( table.empty() ) {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Did not find dalitz decays for particle:" << parent << "\n";
	}

	return table;
	}

	EvtDalitzReso EvtDalitzTable::getResonance(
	std::string shape, EvtDalitzPlot dp, EvtCyclic3::Pair angPair,
	EvtCyclic3::Pair resPair, EvtSpinType::spintype spinType, double mass,
	double width, double FFp, double FFr, double alpha, double aLass,
	double rLass, double BLass, double phiBLass, double RLass, double phiRLass,
	double cutoffLass )
	{
	if ( shape == "RBW" \|\| shape == "RBW_CLEO" ) {
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::RBW_CLEO, FFp, FFr );
	} else if ( shape == "RBW_CLEO_ZEMACH" ) {
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::RBW_CLEO_ZEMACH, FFp, FFr );
	} else if ( shape == "GS" \|\| shape == "GS_CLEO" ) {
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::GS_CLEO, FFp, FFr );
	} else if ( shape == "GS_CLEO_ZEMACH" ) {
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::GS_CLEO_ZEMACH, FFp, FFr );
	} else if ( shape == "GAUSS" \|\| shape == "GAUSS_CLEO" ) {
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::GAUSS_CLEO, FFp, FFr );
	} else if ( shape == "GAUSS_CLEO_ZEMACH" ) {
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::GAUSS_CLEO_ZEMACH, FFp, FFr );
	} else if ( shape == "Flatte" ) {
	return EvtDalitzReso( dp, resPair, mass );
	} else if ( shape == "LASS" ) {
	return EvtDalitzReso( dp, resPair, mass, width, aLass, rLass, BLass,
	phiBLass, RLass, phiRLass, cutoffLass, true );
	} else if ( shape == "NonRes" ) {
	return EvtDalitzReso();
	} else if ( shape == "NonRes_Linear" ) {
	return EvtDalitzReso( dp, resPair, EvtDalitzReso::NON_RES_LIN );
	} else if ( shape == "NonRes_Exp" ) {
	return EvtDalitzReso( dp, resPair, EvtDalitzReso::NON_RES_EXP, alpha );
	} else { //NBW
	if ( shape != "NBW" )
	EvtGenReport( EVTGEN_WARNING, "EvtGen" )
	<< "EvtDalitzTable: shape " << shape
	<< " is unknown. Defaulting to NBW." << endl;
	return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width,
	EvtDalitzReso::NBW, FFp, FFr );
	}
	}

	int EvtDalitzTable::getDaughterPairs(
	EvtId* resDaughter, EvtId* daughter,
	std::vector<std::pair<EvtCyclic3::Pair, EvtCyclic3::Pair>>& angAndResPairs )
	{
	int n( 0 );
	if ( resDaughter[0] == daughter[0] && resDaughter[1] == daughter[1] ) {
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::BC, EvtCyclic3::AB ) );
	n++;
	} else if ( resDaughter[0] == daughter[1] && resDaughter[1] == daughter[0] ) {
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::CA, EvtCyclic3::AB ) );
	n++;
	}

	if ( resDaughter[0] == daughter[1] && resDaughter[1] == daughter[2] ) {
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::CA, EvtCyclic3::BC ) );
	n++;
	} else if ( resDaughter[0] == daughter[2] && resDaughter[1] == daughter[1] ) {
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::AB, EvtCyclic3::BC ) );
	n++;
	}

	if ( resDaughter[0] == daughter[2] && resDaughter[1] == daughter[0] ) {
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::AB, EvtCyclic3::CA ) );
	n++;
	} else if ( resDaughter[0] == daughter[0] && resDaughter[1] == daughter[2] ) {
	angAndResPairs.push_back(
	std::make_pair( EvtCyclic3::BC, EvtCyclic3::CA ) );
	n++;
	}

	return n;
	}

	double EvtDalitzTable::calcProbMax( EvtDalitzPlot dp, EvtDalitzDecayInfo* model )
	{
	double factor = 1.2; //factor to increase our final answer by
	int nStep( 1000 ); //number of steps - total points will be 3nStepnStep

	double maxProb( 0 );
	double min( 0 ), max( 0 ), step( 0 ), min2( 0 ), max2( 0 ), step2( 0 );

	//first do AB, BC
	min = dp.qAbsMin( EvtCyclic3::AB );
	max = dp.qAbsMax( EvtCyclic3::AB );
	step = ( max - min ) / nStep;
	for ( int i = 0; i < nStep; ++i ) {
	double qAB = min + i * step;
	min2 = dp.qMin( EvtCyclic3::BC, EvtCyclic3::AB, qAB );
	max2 = dp.qMax( EvtCyclic3::BC, EvtCyclic3::AB, qAB );
	step2 = ( max2 - min2 ) / nStep;
	for ( int j = 0; j < nStep; ++j ) {
	double qBC = min2 + j * step2;
	EvtDalitzCoord coord( EvtCyclic3::AB, qAB, EvtCyclic3::BC, qBC );
	EvtDalitzPoint point( dp, coord );
	double prob = calcProb( point, model );
	if ( prob > maxProb )
	maxProb = prob;
	}
	}

	//next do BC, CA
	min = dp.qAbsMin( EvtCyclic3::BC );
	max = dp.qAbsMax( EvtCyclic3::BC );
	step = ( max - min ) / nStep;
	for ( int i = 0; i < nStep; ++i ) {
	double qBC = min + i * step;
	min2 = dp.qMin( EvtCyclic3::CA, EvtCyclic3::BC, qBC );
	max2 = dp.qMax( EvtCyclic3::CA, EvtCyclic3::BC, qBC );
	step2 = ( max2 - min2 ) / nStep;
	for ( int j = 0; j < nStep; ++j ) {
	double qCA = min2 + j * step2;
	EvtDalitzCoord coord( EvtCyclic3::BC, qBC, EvtCyclic3::CA, qCA );
	EvtDalitzPoint point( dp, coord );
	double prob = calcProb( point, model );
	if ( prob > maxProb )
	maxProb = prob;
	}
	}

	//finally do CA, AB
	min = dp.qAbsMin( EvtCyclic3::CA );
	max = dp.qAbsMax( EvtCyclic3::CA );
	step = ( max - min ) / nStep;
	for ( int i = 0; i < nStep; ++i ) {
	double qCA = min + i * step;
	min2 = dp.qMin( EvtCyclic3::AB, EvtCyclic3::CA, qCA );
	max2 = dp.qMax( EvtCyclic3::AB, EvtCyclic3::CA, qCA );
	step2 = ( max2 - min2 ) / nStep;
	for ( int j = 0; j < nStep; ++j ) {
	double qAB = min2 + j * step2;
	EvtDalitzCoord coord( EvtCyclic3::CA, qCA, EvtCyclic3::AB, qAB );
	EvtDalitzPoint point( dp, coord );
	double prob = calcProb( point, model );
	if ( prob > maxProb )
	maxProb = prob;
	}
	}
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Largest probability found was " << maxProb << endl;
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Setting probMax to " << factor * maxProb << endl;
	return factor * maxProb;
	}

	double EvtDalitzTable::calcProb( EvtDalitzPoint point, EvtDalitzDecayInfo* model )
	{
	std::vector<std::pair<EvtComplex, EvtDalitzReso>> resonances =
	model->getResonances();

	EvtComplex amp( 0, 0 );
	std::vector<std::pair<EvtComplex, EvtDalitzReso>>::iterator i =
	resonances.begin();
	for ( ; i != resonances.end(); i++ ) {
	std::pair<EvtComplex, EvtDalitzReso> res = ( *i );
	amp += res.first * res.second.evaluate( point );
	}
	return abs2( amp );
	}
	diff --git a/src/EvtGenModels/EvtGenericDalitz.cpp b/src/EvtGenModels/EvtGenericDalitz.cpp
	index 8ea73b9..919a7e0 100644
	--- a/src/EvtGenModels/EvtGenericDalitz.cpp
	+++ b/src/EvtGenModels/EvtGenericDalitz.cpp
	@@ -1,129 +1,129 @@

	/***********************************************************************
	* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
	* *
	* This file is part of EvtGen. *
	* *
	* EvtGen is free software: you can redistribute it and/or modify *
	* it under the terms of the GNU General Public License as published by *
	* the Free Software Foundation, either version 3 of the License, or *
	* (at your option) any later version. *
	* *
	* EvtGen is distributed in the hope that it will be useful, *
	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	* GNU General Public License for more details. *
	* *
	* You should have received a copy of the GNU General Public License *
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/

	#include "EvtGenModels/EvtGenericDalitz.hh"

	#include "EvtGenBase/EvtDalitzPoint.hh"
	#include "EvtGenBase/EvtPDL.hh"
	#include "EvtGenBase/EvtParticle.hh"

	#include "EvtGenModels/EvtDalitzTable.hh"

	std::string EvtGenericDalitz::getName() const
	{
	return "GENERIC_DALITZ";
	}

	EvtDecayBase* EvtGenericDalitz::clone() const
	{
	return new EvtGenericDalitz();
	}

	void EvtGenericDalitz::init()
	{
	checkNArg( 1 );

	EvtId parnum = getParentId();
	EvtId d1 = getDaug( 0 );
	EvtId d2 = getDaug( 1 );
	EvtId d3 = getDaug( 2 );

	std::vector<EvtDalitzDecayInfo> decays =
	- EvtDalitzTable::getInstance( getArgStr( 0 ) )->getDalitzTable( parnum );
	+ EvtDalitzTable::getInstance( getArgStr( 0 ) ).getDalitzTable( parnum );

	std::vector<EvtDalitzDecayInfo>::iterator i = decays.begin();
	for ( ; i != decays.end(); i++ ) {
	EvtId daughter1 = ( *i ).daughter1();
	EvtId daughter2 = ( *i ).daughter2();
	EvtId daughter3 = ( *i ).daughter3();

	if ( d1 == daughter1 && d2 == daughter2 && d3 == daughter3 ) {
	m_d1 = 0;
	m_d2 = 1;
	m_d3 = 2;
	} else if ( d1 == daughter1 && d2 == daughter3 && d3 == daughter2 ) {
	m_d1 = 0;
	m_d2 = 2;
	m_d3 = 1;
	} else if ( d1 == daughter2 && d2 == daughter1 && d3 == daughter3 ) {
	m_d1 = 1;
	m_d2 = 0;
	m_d3 = 2;
	} else if ( d1 == daughter2 && d2 == daughter3 && d3 == daughter1 ) {
	m_d1 = 1;
	m_d2 = 2;
	m_d3 = 0;
	} else if ( d1 == daughter3 && d2 == daughter1 && d3 == daughter2 ) {
	m_d1 = 2;
	m_d2 = 0;
	m_d3 = 1;
	} else if ( d1 == daughter3 && d2 == daughter2 && d3 == daughter1 ) {
	m_d1 = 2;
	m_d2 = 1;
	m_d3 = 0;
	} else {
	continue;
	}

	m_resonances = ( *i ).getResonances();
	setProbMax( ( *i ).getProbMax() );
	return;
	}
	}

	void EvtGenericDalitz::decay( EvtParticle* p )
	{
	p->initializePhaseSpace( getNDaug(), getDaugs() );

	EvtVector4R p4_d1 = p->getDaug( m_d1 )->getP4();
	EvtVector4R p4_d2 = p->getDaug( m_d2 )->getP4();
	EvtVector4R p4_d3 = p->getDaug( m_d3 )->getP4();

	double mA = p->getDaug( m_d1 )->mass();
	double mB = p->getDaug( m_d2 )->mass();
	double mC = p->getDaug( m_d3 )->mass();

	double m2AB = ( p4_d1 + p4_d2 ).mass2();
	double m2CA = ( p4_d1 + p4_d3 ).mass2();
	double m2BC = ( p4_d2 + p4_d3 ).mass2();

	EvtDalitzPoint point( mA, mB, mC, m2AB, m2BC, m2CA );

	EvtComplex amp( 0, 0 );
	std::vector<std::pair<EvtComplex, EvtDalitzReso>>::iterator i =
	m_resonances.begin();
	for ( ; i != m_resonances.end(); i++ ) {
	std::pair<EvtComplex, EvtDalitzReso> res = ( *i );
	amp += res.first * res.second.evaluate( point );
	}

	vertex( amp );
	return;
	}

	std::string EvtGenericDalitz::getParamName( int i )
	{
	switch ( i ) {
	case 0:
	return "xmlFile";
	default:
	return "";
	}
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:01 PM (1 h, 44 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4461189
Default Alt Text: (38 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions