Differential D36 Diff 134 src/EvtGenModels/EvtFlatSqDalitz.cpp

Changeset View

Standalone View

src/EvtGenModels/EvtFlatSqDalitz.cpp

	Show All 14 Lines
	* GNU General Public License for more details. *			* GNU General Public License for more details. *
	* *			* *
	* You should have received a copy of the GNU General Public License *			* You should have received a copy of the GNU General Public License *
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *			* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/			***********************************************************************/

	#include "EvtGenModels/EvtFlatSqDalitz.hh"			#include "EvtGenModels/EvtFlatSqDalitz.hh"

	#include "EvtGenBase/EvtDiracSpinor.hh"
	#include "EvtGenBase/EvtGenKine.hh"			#include "EvtGenBase/EvtGenKine.hh"
	#include "EvtGenBase/EvtPDL.hh"			#include "EvtGenBase/EvtPDL.hh"
	#include "EvtGenBase/EvtParticle.hh"			#include "EvtGenBase/EvtParticle.hh"
	#include "EvtGenBase/EvtPatches.hh"			#include "EvtGenBase/EvtPatches.hh"
	#include "EvtGenBase/EvtReport.hh"			#include "EvtGenBase/EvtReport.hh"
	#include "EvtGenBase/EvtTensor4C.hh"			#include "EvtGenBase/EvtConst.hh"
	#include "EvtGenBase/EvtVector4C.hh"			#include "EvtGenBase/EvtRandom.hh"

	#include <fstream>			#include <fstream>
	#include <stdio.h>			#include <stdio.h>
	#include <stdlib.h>			#include <stdlib.h>
	#include <string>			#include <string>
	using std::fstream;			using std::fstream;

	EvtFlatSqDalitz::~EvtFlatSqDalitz()			EvtFlatSqDalitz::~EvtFlatSqDalitz()
	{			{
	}			}

	std::string EvtFlatSqDalitz::getName()			std::string EvtFlatSqDalitz::getName()
	{			{
	return "FLATSQDALITZ";			return "FLATSQDALITZ";
	}			}

	EvtDecayBase* EvtFlatSqDalitz::clone()			EvtDecayBase* EvtFlatSqDalitz::clone()
	{			{
	return new EvtFlatSqDalitz;			return new EvtFlatSqDalitz;
	}			}

	void EvtFlatSqDalitz::initProbMax()			void EvtFlatSqDalitz::initProbMax()
	{			{
	setProbMax( 1. );			noProbMax();
	}			}

	void EvtFlatSqDalitz::init()			void EvtFlatSqDalitz::init()
	{			{
	// check that there are 0 arguments
	checkNArg( 0 );

	//check there are 3 daughters			//check there are 3 daughters
	checkNDaug( 3 );			checkNDaug( 3 );

				// check that there are 0 arguments
				checkNArg( 0, 2, 4 );

				if ( getNArg() > 0 ) {
				m_mPrimeMin = getArg( 0 );
				m_mPrimeMax = getArg( 1 );
				}
				if ( getNArg() > 2 ) {
				m_thetaPrimeMin = getArg( 2 );
				m_thetaPrimeMax = getArg( 3 );
				}
	}			}

	void EvtFlatSqDalitz::decay( EvtParticle* p )			void EvtFlatSqDalitz::decay( EvtParticle* p )
	{			{
	p->initializePhaseSpace( getNDaug(), getDaugs() );			p->makeDaughters( getNDaug(), getDaugs() );
				p->generateMassTree();
				const double mParent = p->mass();
				EvtParticle* daug1 = p->getDaug( 0 );
				EvtParticle* daug2 = p->getDaug( 1 );
				EvtParticle* daug3 = p->getDaug( 2 );
				const double mDaug1 = daug1->mass();
				const double mDaug2 = daug2->mass();
				const double mDaug3 = daug3->mass();

				// Generate m' and theta'
				const double mPrime = EvtRandom::Flat( m_mPrimeMin, m_mPrimeMax );
				const double thetaPrime = EvtRandom::Flat( m_thetaPrimeMin, m_thetaPrimeMax );

				// calculate m12 and m23
				const double m12 = 0.5 * ( std::cos( mPrime * EvtConst::pi ) + 1 ) *
				( mParent - ( mDaug1 + mDaug2 + mDaug3 ) ) +
				mDaug1 + mDaug2;
				const double m12Sq = m12 * m12;

				const double part1 = std::sqrt(
				-4 * mDaug1 * mDaug1 * m12Sq +
				jbackUnsubmitted Not Done Inline Actions Perhaps reduce some of the code here by introducing variables we can reuse and add functions that calculate equivalent expressions: `const double mDaug1Sq = mDaug1 * mDaug1;` `const double part1 = getMtm( m12Sq, mDaug1Sq, mDaug2Sq );` `const double part2 = getMtm( m12Sq, mDaug3Sq, mParentSq );` jback: Perhaps reduce some of the code here by introducing variables we can reuse and add functions…
				std::pow( mDaug1 * mDaug1 + m12Sq - mDaug2 * mDaug2, 2 ) );
				const double part2 = std::sqrt(
				-4 * mDaug3 * mDaug3 * m12Sq +
				std::pow( -mDaug3 * mDaug3 - m12Sq + mParent * mParent, 2 ) );
				const double m23Sq =
				0.5 * ( ( ( mDaug2 * mDaug2 - mDaug1 * mDaug1 ) *
				( mParent * mParent - mDaug3 * mDaug3 ) +
				part1 * part2 * std::cos( EvtConst::pi * thetaPrime ) ) /
				m12Sq +
				mParent * mParent + mDaug1 * mDaug1 + mDaug2 * mDaug2 +
				mDaug3 * mDaug3 - m12Sq );
				tlathamUnsubmitted Not Done Inline Actions Looking at LauKinematics:: updateSqDPMassSquares, this part looks maybe a bit more complicated than it needs to be. tlatham: Looking at [[ https://laura.hepforge.org/doc/doxygen/latest/LauKinematics_8cc_source.
				krepsAuthorUnsubmitted Done Inline Actions I strongly suspect that what I did was basically the same just dumped to a single expression, but I appreciate that it was not very transparent and easy to understand. The new version should fix that. kreps: I strongly suspect that what I did was basically the same just dumped to a single expression…

	double mB = p->mass();			// Turn m12 and m23 into momenta
	double m1 = p->getDaug( 0 )->mass();			EvtGenKine::ThreeBodyKine(m12Sq, m23Sq, p);
	double m2 = p->getDaug( 1 )->mass();
	double m3 = p->getDaug( 2 )->mass();

	EvtVector4R p4_1 = p->getDaug( 0 )->getP4();
	EvtVector4R p4_2 = p->getDaug( 1 )->getP4();
	EvtVector4R p4_3 = p->getDaug( 2 )->getP4();

	EvtVector4R p4_12 = p4_1 + p4_2;
	EvtVector4R p4_13 = p4_1 + p4_3;
	// do not compute p4_23 to avoid breaking p4 conservation ???
	EvtVector4R p4_23 = p4_2 + p4_3;

	double m12 = p4_12.mass();
	double m13 = p4_13.mass();
	double m23 = p4_23.mass();

	double m12norm = 2 * ( ( m12 - ( m1 + m2 ) ) / ( mB - ( m1 + m2 + m3 ) ) ) -
	1;
	double mPrime = acos( m12norm ) / EvtConst::pi;
	double thPrime = acos( ( m12 * m12 * ( m23 * m23 - m13 * m13 ) -
	( m2 * m2 - m1 * m1 ) * ( mB * mB - m3 * m3 ) ) /
	( sqrt( pow( m12 * m12 + m1 * m1 - m2 * m2, 2 ) -
	4 * m12 * m12 * m1 * m1 ) *
	sqrt( pow( -m12 * m12 + mB * mB - m3 * m3, 2 ) -
	4 * m12 * m12 * m3 * m3 ) ) ) /
	EvtConst::pi;

	double p3st = sqrt( pow( mB * mB - m3 * m3 - m12 * m12, 2 ) -
	pow( 2 * m12 * m3, 2 ) ) /
	( 2 * m12 );
	double p1st = sqrt( pow( m2 * m2 - m1 * m1 - m12 * m12, 2 ) -
	pow( 2 * m12 * m1, 2 ) ) /
	( 2 * m12 );

	double jacobian = 2 * pow( EvtConst::pi, 2 ) * sin( EvtConst::pi * mPrime ) *
	sin( EvtConst::pi * thPrime ) * p1st * p3st * m12 *
	( mB - ( m1 + m2 + m3 ) );

	double prob = 1. / jacobian; //pow(1./(jacobian),2);

	// std::cout << mB << " " << mPrime << " " << thPrime << " " << prob << std::endl;

	setProb( prob );

	if ( prob < 1 )
	setProb( prob );
	else
	setProb( 1. );

	return;			return;
	}			}