D47.1759146693.diff
No OneTemporary
Actions

Size

30 KB

Referenced Files

None

Subscribers

None

D47.1759146693.diff
View Options

	diff --git a/EvtGenBase/EvtKine.hh b/EvtGenBase/EvtKine.hh
	--- a/EvtGenBase/EvtKine.hh
	+++ b/EvtGenBase/EvtKine.hh
	@@ -62,4 +62,9 @@
	EvtComplex wignerD( int j, int m1, int m2, double phi, double theta,
	double gamma );

	+//
	+// Function to calculate momentum of daughters in two body decay in mothers
	+// rest frame.
	+double twoBodyMomentum( const double M, const double m1, const double m2 );
	+
	#endif
	diff --git a/EvtGenModels/EvtFourBodyPhsp.hh b/EvtGenModels/EvtFourBodyPhsp.hh
	new file mode 100644
	--- /dev/null
	+++ b/EvtGenModels/EvtFourBodyPhsp.hh
	@@ -0,0 +1,92 @@
	+
	+/***********************************************************************
	+* 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 EVTFOURBODYPHSP_HH
	+#define EVTFOURBODYPHSP_HH
	+
	+#include <vector>
	+#include <array>
	+#include <utility>
	+
	+#include "EvtGenBase/EvtDecayProb.hh"
	+
	+class EvtParticle;
	+
	+class EvtFourBodyPhsp : public EvtDecayProb {
	+ public:
	+ enum Shape
	+ {
	+ rectangle = 1,
	+ trapezoid = 2,
	+ pentagon = 3,
	+ variable = 4
	+ };
	+
	+ std::string getName() override;
	+ EvtDecayBase* clone() override;
	+
	+ void init() override;
	+ void initProbMax() override;
	+
	+ void decay( EvtParticle* parent ) override;
	+
	+ private:
	+ std::array<double, 4> phspFactor( const double mM, const double m12,
	+ const double m34,
	+ std::array<double, 4>& daughters ) const;
	+ Shape determineBoundaryShape( const double m12Min, const double m12Max,
	+ const double m34Max,
	+ const double mMother ) const;
	+
	+ std::pair<double, double> generatePairMasses(
	+ const double m12Min, const double m12Max, const double m34Min,
	+ const double m34Max, const double mMother,
	+ const EvtFourBodyPhsp::Shape shape ) const;
	+ std::pair<double, double> generateRectangle( const double m12Min,
	+ const double m12Max,
	+ const double m34Min,
	+ const double m34Max ) const;
	+ std::pair<double, double> generateTrapezoid( const double m12Min,
	+ const double m12Max,
	+ const double m34Min,
	+ const double mMother ) const;
	+
	+ std::array<double, 4> m_daughterMasses{ -1, -1, -1, -1 };
	+
	+ double m_m12Min;
	+ double m_m12Max;
	+ double m_m34Min;
	+ double m_m34Max;
	+
	+ double m_trapNorm;
	+ double m_trapCoeff1;
	+ double m_trapCoeff2;
	+
	+ double m_pentagonSplit;
	+ double m_pentagonFraction;
	+
	+ Shape m_boundaryShape;
	+
	+ bool m_stableMother{true};
	+ bool m_stableDaughters{true};
	+ bool m_fixedBoundary{true};
	+};
	+
	+#endif
	diff --git a/src/EvtGenBase/EvtKine.cpp b/src/EvtGenBase/EvtKine.cpp
	--- a/src/EvtGenBase/EvtKine.cpp
	+++ b/src/EvtGenBase/EvtKine.cpp
	@@ -28,7 +28,7 @@
	#include "EvtGenBase/EvtVector4R.hh"
	#include "EvtGenBase/EvtdFunction.hh"

	-#include <math.h>
	+#include <cmath>

	double EvtDecayAngle( const EvtVector4R& p, const EvtVector4R& q,
	const EvtVector4R& d )
	@@ -41,7 +41,7 @@
	double md2 = d.mass2();

	double cost = ( pd * mq2 - pq * qd ) /
	- sqrt( ( pq * pq - mq2 * mp2 ) * ( qd * qd - mq2 * md2 ) );
	+ std::sqrt( ( pq * pq - mq2 * mp2 ) * ( qd * qd - mq2 * md2 ) );

	return cost;
	}
	@@ -80,7 +80,7 @@
	x = d1_perp.dot( h1_perp );
	y = d1_prime.dot( h1_perp );

	- double chi = atan2( y, x );
	+ double chi = std::atan2( y, x );

	if ( chi < 0.0 )
	chi += EvtConst::twoPi;
	@@ -99,7 +99,7 @@
	double pq = p * q;

	return q.mass() * ( p * l ) /
	- sqrt( -( pq * pq - p.mass2() * q.mass2() ) * l.mass2() );
	+ std::sqrt( -( pq * pq - p.mass2() * q.mass2() ) * l.mass2() );
	}

	// Calculate phi using the given 4 vectors (all in the same frame)
	@@ -117,9 +117,9 @@

	double y = p.scalartripler3( z, q, d ) + alpha * p.scalartripler3( z, q, q );
	double x = ( zq * ( qd + alpha * q2 ) - q2 * ( zd + alpha * zq ) ) /
	- sqrt( q2 );
	+ std::sqrt( q2 );

	- double phi = atan2( y, x );
	+ double phi = std::atan2( y, x );

	return phi < 0 ? ( phi + EvtConst::twoPi ) : phi;
	}
	@@ -131,3 +131,15 @@

	return exp( gp ) * EvtdFunction::d( j, m1, m2, theta ) * exp( gm );
	}
	+
	+double twoBodyMomentum( const double M, const double m1, const double m2 )
	+{
	+ const double MSq = M * M;
	+ const double mSum = m1 + m2;
	+ const double mDiff = m1 - m2;
	+ double result = ( MSq - mDiff * mDiff ) * ( MSq - mSum * mSum );
	+ if ( result < 0 ) {
	+ return 0;
	+ }
	+ return std::sqrt( result ) / ( 2 * M );
	+}
	diff --git a/src/EvtGenModels/EvtFourBodyPhsp.cpp b/src/EvtGenModels/EvtFourBodyPhsp.cpp
	new file mode 100644
	--- /dev/null
	+++ b/src/EvtGenModels/EvtFourBodyPhsp.cpp
	@@ -0,0 +1,446 @@
	+
	+/***********************************************************************
	+* 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 <cmath>
	+
	+#include "EvtGenBase/EvtKine.hh"
	+#include "EvtGenBase/EvtPDL.hh"
	+#include "EvtGenBase/EvtParticle.hh"
	+#include "EvtGenBase/EvtRandom.hh"
	+#include "EvtGenBase/EvtReport.hh"
	+#include "EvtGenModels/EvtFourBodyPhsp.hh"
	+
	+std::string EvtFourBodyPhsp::getName()
	+{
	+ return "FOURBODYPHSP";
	+}
	+
	+EvtDecayBase* EvtFourBodyPhsp::clone()
	+{
	+ return new EvtFourBodyPhsp;
	+}
	+
	+std::array<double, 4> EvtFourBodyPhsp::phspFactor(
	+ const double mM, const double m12, const double m34,
	+ std::array<double, 4>& daughters ) const
	+{
	+ std::array<double, 4> result;
	+ result[1] = twoBodyMomentum( mM, m12, m34 );
	+ result[2] = twoBodyMomentum( m12, daughters[0], daughters[1] );
	+ result[3] = twoBodyMomentum( m34, daughters[2], daughters[3] );
	+ result[0] = result[1] * result[2] * result[3];
	+
	+ return result;
	+}
	+
	+void EvtFourBodyPhsp::init()
	+{
	+ // Check that we have right number of daughters
	+ checkNDaug( 4 );
	+
	+ // Check whether mother is quasi-stable
	+ auto parent = getParentId();
	+ double width = EvtPDL::getWidth( parent );
	+ if ( width > 1e-6 ) {
	+ m_stableMother = false;
	+ }
	+
	+ // Check whether all daughters are stable
	+ for ( int i = 0; i < 4; ++i ) {
	+ auto daughter = getDaug( i );
	+ width = EvtPDL::getWidth( daughter );
	+ if ( width > 1e-6 ) {
	+ m_stableDaughters = false;
	+ m_daughterMasses[i] = EvtPDL::getMinMass( daughter );
	+ } else {
	+ m_daughterMasses[i] = EvtPDL::getMass( daughter );
	+ }
	+ }
	+
	+ // check correct number of arguments
	+ checkNArg( 0, 2, 4 );
	+ double mass1 = m_daughterMasses[0];
	+ double mass2 = m_daughterMasses[1];
	+ double mass3 = m_daughterMasses[2];
	+ double mass4 = m_daughterMasses[3];
	+ double mMother = EvtPDL::getMaxMass( parent );
	+ if ( getNArg() > 2 ) {
	+ m_m12Min = getArg( 0 );
	+ m_m12Max = getArg( 1 );
	+ m_m34Min = getArg( 2 );
	+ m_m34Max = getArg( 3 );
	+ } else {
	+ if ( getNArg() > 0 ) {
	+ m_m12Min = getArg( 0 );
	+ m_m12Max = getArg( 1 );
	+ } else {
	+ m_m12Min = mass1 + mass2;
	+ m_m12Max = mMother - mass3 - mass4;
	+ }
	+ m_m34Min = mass3 + mass4;
	+ m_m34Max = mMother - mass1 - mass2;
	+ if ( m_stableDaughters == false \|\| m_stableMother == false ) {
	+ m_fixedBoundary = false;
	+ }
	+ }
	+ // Make sure that we have correct boundaries
	+ if ( m_m12Min < mass1 + mass2 ) {
	+ m_m12Min = mass1 + mass2;
	+ }
	+ if ( m_m12Max > mMother - mass3 - mass4 ) {
	+ m_m12Max = mMother - mass3 - mass4;
	+ }
	+ if ( m_m34Min < mass3 + mass4 ) {
	+ m_m34Min = mass3 + mass4;
	+ }
	+ if ( m_m34Max > mMother - mass1 - mass2 ) {
	+ m_m34Max = mMother - mass1 - mass2;
	+ }
	+
	+ if ( m_stableDaughters && m_stableMother ) {
	+ m_boundaryShape = determineBoundaryShape( m_m12Min, m_m12Max, m_m34Max,
	+ mMother );
	+ } else {
	+ m_boundaryShape = Shape::variable;
	+ }
	+ // If we have fixed boundary, we can precalculate some variables for
	+ // m12 and m34 generation
	+ if ( m_fixedBoundary ) {
	+ if ( m_boundaryShape == Shape::trapezoid ) {
	+ const double m12Diff = m_m12Max - m_m12Min;
	+ const double minSum = m_m12Min + m_m34Min;
	+ m_trapNorm = ( mMother - m_m34Min ) * m12Diff -
	+ 0.5 * ( m12Diff * ( m_m12Max + m_m12Min ) );
	+ m_trapCoeff1 = mMother - m_m34Min;
	+ m_trapCoeff2 = mMother * mMother - 2 * mMother * minSum +
	+ minSum * minSum;
	+ }
	+ if ( m_boundaryShape == Shape::pentagon ) {
	+ m_pentagonSplit = mMother - m_m34Max;
	+ const double area1 = ( m_pentagonSplit - m_m12Min ) *
	+ ( m_m34Max - m_m34Min );
	+ const double pm12Diff = m_m12Max - m_pentagonSplit;
	+ const double area2 = 0.5 * pm12Diff *
	+ ( mMother + m_m34Max - m_m12Max ) -
	+ pm12Diff * m_m34Min;
	+ m_pentagonFraction = area1 / ( area1 + area2 );
	+ const double m12Diff = m_m12Max - m_pentagonSplit;
	+ const double minSum = m_pentagonSplit + m_m34Min;
	+ m_trapNorm = ( mMother - m_m34Min ) * m12Diff -
	+ 0.5 * ( m12Diff * ( m_m12Max + m_pentagonSplit ) );
	+ m_trapCoeff1 = mMother - m_m34Min;
	+ m_trapCoeff2 = mMother * mMother - 2 * mMother * minSum +
	+ minSum * minSum;
	+ }
	+ }
	+}
	+
	+void EvtFourBodyPhsp::initProbMax()
	+{
	+ double startM12 = m_m12Min + ( m_m12Max - m_m12Min ) / 20.;
	+ double startM34 = m_m34Min + ( m_m34Max - m_m34Min ) / 20.;
	+ bool contCond = true;
	+ int iteration = 0;
	+
	+ auto parent = getParentId();
	+ double mMother = EvtPDL::getMaxMass( parent );
	+
	+ double funcValue = 0;
	+ while (contCond){
	+ ++iteration;
	+ double currentM12 = startM12;
	+ double currentM34 = startM34;
	+ funcValue = phspFactor( mMother, currentM12, currentM34,
	+ m_daughterMasses )[0];
	+ // Maximum along m12
	+ double step = ( m_m12Max - m_m12Min ) / 100.;
	+ while ( step > 1e-4 ) {
	+ double point1 = currentM12 + step;
	+ if ( point1 > m_m12Max ) {
	+ point1 = m_m12Max;
	+ }
	+ if ( currentM34 > mMother - point1 ) {
	+ point1 = mMother - currentM34;
	+ }
	+ double point2 = currentM12 - step;
	+ if ( point2 < m_m12Min ) {
	+ point2 = m_m12Min;
	+ }
	+ double value1 = phspFactor( mMother, point1, currentM34,
	+ m_daughterMasses )[0];
	+ double value2 = phspFactor( mMother, point2, currentM34,
	+ m_daughterMasses )[0];
	+ if ( value1 > funcValue && value1 > value2 ) {
	+ currentM12 = point1;
	+ funcValue = value1;
	+ } else if ( value2 > funcValue ) {
	+ currentM12 = point2;
	+ funcValue = value2;
	+ }
	+ step /= 2.;
	+ }
	+ // Maximum along m34
	+ step = ( mMother - currentM12 - m_m34Min ) / 100.;
	+ while ( step > 1e-4 ) {
	+ double point1 = currentM34 + step;
	+ if ( point1 > m_m34Max ) {
	+ point1 = m_m34Max;
	+ }
	+ if ( point1 > mMother - currentM12 ) {
	+ point1 = mMother - currentM12;
	+ }
	+ double point2 = currentM34 - step;
	+ if ( point2 < m_m34Min ) {
	+ point2 = m_m34Min;
	+ }
	+ double value1 = phspFactor( mMother, currentM12, point1,
	+ m_daughterMasses )[0];
	+ double value2 = phspFactor( mMother, currentM12, point2,
	+ m_daughterMasses )[0];
	+ if ( value1 > funcValue && value1 > value2 ) {
	+ currentM34 = point1;
	+ funcValue = value1;
	+ } else if ( value2 > funcValue ) {
	+ currentM34 = point2;
	+ funcValue = value2;
	+ }
	+ step /= 2.;
	+ }
	+
	+ // Check termination condition
	+ double m12Diff = currentM12 - startM12;
	+ double m34Diff = currentM34 - startM34;
	+ double distSq = m12Diff * m12Diff + m34Diff * m34Diff;
	+ if (distSq < 1e-8 \|\| iteration > 50){
	+ contCond = false;
	+ }
	+ startM12 = currentM12;
	+ startM34 = currentM34;
	+ }
	+
	+ setProbMax( funcValue * 1.05 );
	+}
	+
	+void EvtFourBodyPhsp::decay( EvtParticle* parent )
	+{
	+
	+ parent->makeDaughters( getNDaug(), getDaugs() );
	+ bool massTreeStatus = parent->generateMassTree();
	+ if ( !massTreeStatus ) {
	+ EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	+ << "Failed to generate daughters masses in EvtFourBodyPhsp."
	+ << std::endl;
	+ ::abort();
	+ }
	+
	+ double mMother = parent->mass();
	+
	+ // Need to check whether boundaries are OK and whether we need to work
	+ // out boundary shape
	+ double cM12Min, cM12Max;
	+ double cM34Min, cM34Max;
	+ EvtFourBodyPhsp::Shape cShape;
	+ if ( m_fixedBoundary ) {
	+ cM12Min = m_m12Min;
	+ cM12Max = m_m12Max;
	+ cM34Min = m_m34Min;
	+ cM34Max = m_m34Max;
	+ cShape = m_boundaryShape;
	+ } else {
	+ // In this case at least one particle has non-zero width and thus
	+ // boundaries and shape of the region can change
	+ for ( int i = 0; i < 4; ++i ) {
	+ auto daughter = parent->getDaug( i );
	+ m_daughterMasses[i] = daughter->mass();
	+ }
	+ cM12Min = m_m12Min > ( m_daughterMasses[0] + m_daughterMasses[1] )
	+ ? m_m12Min
	+ : m_daughterMasses[0] + m_daughterMasses[1];
	+ cM12Max = m_m12Max <
	+ ( mMother - m_daughterMasses[2] - m_daughterMasses[3] )
	+ ? m_m12Max
	+ : mMother - m_daughterMasses[2] - m_daughterMasses[3];
	+ cM34Min = m_m34Min > ( m_daughterMasses[2] + m_daughterMasses[3] )
	+ ? m_m34Min
	+ : m_daughterMasses[2] + m_daughterMasses[3];
	+ cM34Max = m_m34Max <
	+ ( mMother - m_daughterMasses[0] - m_daughterMasses[1] )
	+ ? m_m34Max
	+ : mMother - m_daughterMasses[0] - m_daughterMasses[1];
	+ cShape = determineBoundaryShape( cM12Min, cM12Max, cM34Max, mMother );
	+ }
	+
	+ // Generate m12 and m34
	+ auto masses = generatePairMasses( cM12Min, cM12Max, cM34Min, cM34Max,
	+ mMother, cShape );
	+ const double m12 = masses.first;
	+ const double m34 = masses.second;
	+
	+ // calculate probability, it will return array with 4 elements with
	+ // probability, q, p1 and p3
	+ auto probEval = phspFactor( mMother, m12, m34, m_daughterMasses );
	+ setProb( probEval[0] );
	+
	+ // initialise kinematics
	+ const double cosTheta1 = EvtRandom::Flat(-1.0, 1.0);
	+ const double sinTheta1 = std::sqrt( 1 - cosTheta1 * cosTheta1 );
	+ const double cosTheta3 = EvtRandom::Flat(-1.0, 1.0);
	+ const double sinTheta3 = std::sqrt( 1 - cosTheta3 * cosTheta3 );
	+ const double phi = EvtRandom::Flat( 0., EvtConst::twoPi );
	+ // m12 and m34 are put along z-axis, 1 and 2 go to x-z plane and 3-4
	+ // plane is rotated by phi compared to 1-2 plane. All momenta are set
	+ // in 12 and 34 rest frames and then boosted to parent rest frame
	+ const double p1x = probEval[2] * sinTheta1;
	+ const double p1z = probEval[2] * cosTheta1;
	+ const double p1Sq = probEval[2] * probEval[2];
	+ const double en1 = std::sqrt( m_daughterMasses[0] * m_daughterMasses[0] +
	+ p1Sq );
	+ const double en2 = std::sqrt( m_daughterMasses[1] * m_daughterMasses[1] +
	+ p1Sq );
	+ const double p3T = probEval[3] * sinTheta3;
	+ const double p3x = p3T * std::cos( phi );
	+ const double p3y = p3T * std::sin( phi );
	+ const double p3z = probEval[3] * cosTheta3;
	+ const double p3Sq = probEval[3] * probEval[3];
	+ const double en3 = std::sqrt( m_daughterMasses[2] * m_daughterMasses[2] +
	+ p3Sq );
	+ const double en4 = std::sqrt( m_daughterMasses[3] * m_daughterMasses[3] +
	+ p3Sq );
	+
	+ EvtVector4R mom1( en1, p1x, 0.0, p1z );
	+ EvtVector4R mom2( en2, -p1x, 0.0, -p1z );
	+ EvtVector4R mom3( en3, p3x, p3y, p3z );
	+ EvtVector4R mom4( en4, -p3x, -p3y, -p3z );
	+
	+ const double qSq = probEval[1] * probEval[1];
	+ const double en12 = std::sqrt( m12 * m12 + qSq );
	+ const double en34 = std::sqrt( m34 * m34 + qSq );
	+ EvtVector4R q12( en12, 0.0, 0.0, probEval[1] );
	+ EvtVector4R q34( en34, 0.0, 0.0, -probEval[1] );
	+ mom1.applyBoostTo( q12 );
	+ mom2.applyBoostTo( q12 );
	+ mom3.applyBoostTo( q34 );
	+ mom4.applyBoostTo( q34 );
	+
	+ // As final step, rotate everything randomly in space
	+ const double euler1 = EvtRandom::Flat( 0., EvtConst::twoPi );
	+ const double euler2 = std::acos( EvtRandom::Flat( -1.0, 1.0 ) );
	+ const double euler3 = EvtRandom::Flat( 0., EvtConst::twoPi );
	+ mom1.applyRotateEuler( euler1, euler2, euler3 );
	+ mom2.applyRotateEuler( euler1, euler2, euler3 );
	+ mom3.applyRotateEuler( euler1, euler2, euler3 );
	+ mom4.applyRotateEuler( euler1, euler2, euler3 );
	+
	+ // Set momenta for daughters
	+ auto daug = parent->getDaug( 0 );
	+ daug->init( daug->getId(), mom1 );
	+ daug = parent->getDaug( 1 );
	+ daug->init( daug->getId(), mom2 );
	+ daug = parent->getDaug( 2 );
	+ daug->init( daug->getId(), mom3 );
	+ daug = parent->getDaug( 3 );
	+ daug->init( daug->getId(), mom4 );
	+}
	+
	+EvtFourBodyPhsp::Shape EvtFourBodyPhsp::determineBoundaryShape(
	+ const double m12Min, const double m12Max, const double m34Max,
	+ const double mMother ) const
	+{
	+ double maxY = mMother - m12Min;
	+ const bool corner1 = m34Max < maxY;
	+ maxY = mMother - m12Max;
	+ const bool corner2 = m34Max < maxY;
	+
	+ if ( corner1 && corner2 ) {
	+ return Shape::rectangle;
	+ } else if ( !corner1 && !corner2 ) {
	+ return Shape::trapezoid;
	+ }
	+ return Shape::pentagon;
	+}
	+
	+std::pair<double, double> EvtFourBodyPhsp::generatePairMasses(
	+ const double m12Min, const double m12Max, const double m34Min,
	+ const double m34Max, const double mMother,
	+ const EvtFourBodyPhsp::Shape shape ) const
	+{
	+ switch ( shape ) {
	+ case EvtFourBodyPhsp::Shape::rectangle:
	+ return generateRectangle( m12Min, m12Max, m34Min, m34Max );
	+ break;
	+ case EvtFourBodyPhsp::Shape::trapezoid:
	+ return generateTrapezoid( m12Min, m12Max, m34Min, mMother );
	+ break;
	+ case EvtFourBodyPhsp::Shape::pentagon:
	+ double split, fraction;
	+ if ( m_fixedBoundary ) {
	+ split = m_pentagonSplit;
	+ fraction = m_pentagonFraction;
	+ } else {
	+ split = mMother - m34Max;
	+ const double area1 = ( split - m12Min ) * ( m34Max - m34Min );
	+ const double pm12Diff = m12Max - split;
	+ const double area2 = 0.5 * pm12Diff *
	+ ( mMother + m34Max - m12Max ) -
	+ pm12Diff * m34Min;
	+ fraction = area1 / ( area1 + area2 );
	+ }
	+ if ( EvtRandom::Flat() < fraction ) {
	+ return generateRectangle( m12Min, split, m34Min, m34Max );
	+ } else {
	+ return generateTrapezoid( split, m12Max, m34Min, mMother );
	+ }
	+ break;
	+ default:
	+ return std::make_pair( m12Min, m34Min );
	+ break;
	+ }
	+}
	+
	+std::pair<double, double> EvtFourBodyPhsp::generateRectangle(
	+ const double m12Min, const double m12Max, const double m34Min,
	+ const double m34Max ) const
	+{
	+ return std::make_pair( EvtRandom::Flat( m12Min, m12Max ),
	+ EvtRandom::Flat( m34Min, m34Max ) );
	+}
	+
	+std::pair<double, double> EvtFourBodyPhsp::generateTrapezoid(
	+ const double m12Min, const double m12Max, const double m34Min,
	+ const double mMother ) const
	+{
	+ double norm, coeff1, coeff2;
	+ if ( m_fixedBoundary ) {
	+ norm = m_trapNorm;
	+ coeff1 = m_trapCoeff1;
	+ coeff2 = m_trapCoeff2;
	+ } else {
	+ const double m12Diff = m12Max - m12Min;
	+ const double minSum = m12Min + m34Min;
	+ norm = ( mMother - m34Min ) * m12Diff -
	+ 0.5 * ( m12Diff * ( m12Max + m12Min ) );
	+ coeff1 = mMother - m34Min;
	+ coeff2 = mMother * mMother - 2 * mMother * minSum + minSum * minSum;
	+ }
	+ const double rnd = EvtRandom::Flat();
	+ const double m12 = coeff1 - std::sqrt( -2.0 * rnd * norm + coeff2 );
	+ const double m34 = EvtRandom::Flat( m34Min, mMother - m12 );
	+ return std::make_pair( m12, m34 );
	+}
	diff --git a/src/EvtGenModels/EvtModelReg.cpp b/src/EvtGenModels/EvtModelReg.cpp
	--- a/src/EvtGenModels/EvtModelReg.cpp
	+++ b/src/EvtGenModels/EvtModelReg.cpp
	@@ -160,6 +160,7 @@
	#include "EvtGenModels/EvtbsToLLLL.hh"
	#include "EvtGenModels/EvtbsToLLLLHyperCP.hh"
	#include "EvtGenModels/EvtThreeBodyPhsp.hh"
	+#include "EvtGenModels/EvtFourBodyPhsp.hh"
	#include "EvtGenModels/EvtEtaLLPiPi.hh"

	#include <assert.h>
	@@ -342,5 +343,6 @@
	modelist.registerModel( new EvtPsi2JpsiPiPi );

	modelist.registerModel( new EvtThreeBodyPhsp );
	+ modelist.registerModel( new EvtFourBodyPhsp );
	modelist.registerModel( new EvtEtaLLPiPi );
	}
	diff --git a/test/evtgenlhc_test1.cc b/test/evtgenlhc_test1.cc
	--- a/test/evtgenlhc_test1.cc
	+++ b/test/evtgenlhc_test1.cc
	@@ -64,6 +64,7 @@
	#include "TH1.h"
	#include "TH2.h"
	#include "TROOT.h"
	+#include "TTree.h"
	#include "TString.h"

	#include <cmath>
	@@ -158,6 +159,7 @@
	void runBaryonic( int nEvent, EvtGen& myGenerator );
	void run3BPhspRegion( int nEvent, EvtGen& myGenerator );
	void runFlatSqDalitz( int nEvent, EvtGen& myGenerator );
	+void runFourBody( int nevent, EvtGen& myGenerator );

	int main( int argc, char* argv[] )
	{
	@@ -541,6 +543,12 @@
	runFlatSqDalitz( nevent, myGenerator );
	}

	+ if ( !strcmp( argv[1], "4bodyPhsp" ) ) {
	+ int nevent = atoi( argv[2] );
	+ EvtRadCorr::setNeverRadCorr();
	+ runFourBody( nevent, myGenerator );
	+ }
	+
	//*******************************************************
	//test of the rotations and boosts performed in EvtGen.
	// Added by Lange and Ryd Jan 5,2000.
	@@ -5819,3 +5827,125 @@
	file->Close();
	EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "SUCCESS\n";
	}
	+
	+void runFourBody( int nevent, EvtGen& myGenerator )
	+{
	+ TFile* file = new TFile( "fourBody.root", "RECREATE" );
	+
	+ double m12, m13, m14, m23, m24, m34, m123, m124, m234;
	+ double mB, m1, m2, m3, m4;
	+ double pBe, pBx, pBy, pBz;
	+ double p1e, p1x, p1y, p1z;
	+ double p2e, p2x, p2y, p2z;
	+ double p3e, p3x, p3y, p3z;
	+ double p4e, p4x, p4y, p4z;
	+ double theta1, theta3, chi;
	+ TTree* tree = new TTree( "tree", "");
	+ tree->Branch("m12", &m12, "m12/D");
	+ tree->Branch("m13", &m13, "m13/D");
	+ tree->Branch("m14", &m14, "m14/D");
	+ tree->Branch("m23", &m23, "m23/D");
	+ tree->Branch("m24", &m24, "m24/D");
	+ tree->Branch("m34", &m34, "m34/D");
	+ tree->Branch("m123", &m123, "m123/D");
	+ tree->Branch("m124", &m124, "m124/D");
	+ tree->Branch("m234", &m234, "m234/D");
	+ tree->Branch("mB", &mB, "mB/D");
	+ tree->Branch("m1", &m1, "m1/D");
	+ tree->Branch("m2", &m2, "m2/D");
	+ tree->Branch("m3", &m3, "m3/D");
	+ tree->Branch("m4", &m4, "m4/D");
	+ tree->Branch("pBe", &pBe, "pBe/D");
	+ tree->Branch("pBx", &pBx, "pBx/D");
	+ tree->Branch("pBy", &pBy, "pBy/D");
	+ tree->Branch("pBz", &pBz, "pBz/D");
	+ tree->Branch("p1e", &p1e, "p1e/D");
	+ tree->Branch("p1x", &p1x, "p1x/D");
	+ tree->Branch("p1y", &p1y, "p1y/D");
	+ tree->Branch("p1z", &p1z, "p1z/D");
	+ tree->Branch("p2e", &p2e, "p2e/D");
	+ tree->Branch("p2x", &p2x, "p2x/D");
	+ tree->Branch("p2y", &p2y, "p2y/D");
	+ tree->Branch("p2z", &p2z, "p2z/D");
	+ tree->Branch("p3e", &p3e, "p3e/D");
	+ tree->Branch("p3x", &p3x, "p3x/D");
	+ tree->Branch("p3y", &p3y, "p3y/D");
	+ tree->Branch("p3z", &p3z, "p3z/D");
	+ tree->Branch("p4e", &p4e, "p4e/D");
	+ tree->Branch("p4x", &p4x, "p4x/D");
	+ tree->Branch("p4y", &p4y, "p4y/D");
	+ tree->Branch("p4z", &p4z, "p4z/D");
	+ tree->Branch("theta1", &theta1, "theta1/D");
	+ tree->Branch("theta3", &theta3, "theta3/D");
	+ tree->Branch("chi", &chi, "chi/D");
	+
	+ int count = 1;
	+
	+ char udecay_name[100];
	+ strcpy( udecay_name, "exampleFiles/4BodyPhsp.DEC" );
	+ myGenerator.readUDecay( udecay_name );
	+
	+ static EvtId B = EvtPDL::getId( std::string( "B+" ) );
	+
	+ do {
	+ EvtVector4R pinit( EvtPDL::getMass( B ), 0.0, 0.0, 0.0 );
	+
	+ EvtParticle* root_part = EvtParticleFactory::particleFactory( B, pinit );
	+
	+ myGenerator.generateDecay( root_part );
	+
	+ mB = root_part->mass();
	+ m1 = root_part->getDaug( 0 )->mass();
	+ m2 = root_part->getDaug( 1 )->mass();
	+ m3 = root_part->getDaug( 2 )->mass();
	+ m4 = root_part->getDaug( 3 )->mass();
	+
	+ EvtParticle* daug1 = root_part->getDaug( 0 );
	+ EvtParticle* daug2 = root_part->getDaug( 1 );
	+ EvtParticle* daug3 = root_part->getDaug( 2 );
	+ EvtParticle* daug4 = root_part->getDaug( 3 );
	+ m12 = ( daug1->getP4() + daug2->getP4() ).mass();
	+ m13 = ( daug1->getP4() + daug3->getP4() ).mass();
	+ m14 = ( daug1->getP4() + daug4->getP4() ).mass();
	+ m23 = ( daug2->getP4() + daug3->getP4() ).mass();
	+ m24 = ( daug2->getP4() + daug4->getP4() ).mass();
	+ m34 = ( daug3->getP4() + daug4->getP4() ).mass();
	+ m123 = ( daug1->getP4() + daug2->getP4() + daug3->getP4() ).mass();
	+ m124 = ( daug1->getP4() + daug2->getP4() + daug4->getP4() ).mass();
	+ m234 = ( daug2->getP4() + daug3->getP4() + daug4->getP4() ).mass();
	+ pBe = root_part->getP4().get( 0 );
	+ pBx = root_part->getP4().get( 1 );
	+ pBy = root_part->getP4().get( 2 );
	+ pBz = root_part->getP4().get( 3 );
	+ p1e = daug1->getP4().get( 0 );
	+ p1x = daug1->getP4().get( 1 );
	+ p1y = daug1->getP4().get( 2 );
	+ p1z = daug1->getP4().get( 3 );
	+ p2e = daug2->getP4().get( 0 );
	+ p2x = daug2->getP4().get( 1 );
	+ p2y = daug2->getP4().get( 2 );
	+ p2z = daug2->getP4().get( 3 );
	+ p3e = daug3->getP4().get( 0 );
	+ p3x = daug3->getP4().get( 1 );
	+ p3y = daug3->getP4().get( 2 );
	+ p3z = daug3->getP4().get( 3 );
	+ p4e = daug4->getP4().get( 0 );
	+ p4x = daug4->getP4().get( 1 );
	+ p4y = daug4->getP4().get( 2 );
	+ p4z = daug4->getP4().get( 3 );
	+
	+ theta1 = EvtDecayAngle( root_part->getP4(),
	+ daug1->getP4() + daug2->getP4(), daug1->getP4() );
	+ theta3 = EvtDecayAngle( root_part->getP4(),
	+ daug3->getP4() + daug4->getP4(), daug3->getP4() );
	+ chi = EvtDecayAngleChi( root_part->getP4(), daug1->getP4(),
	+ daug2->getP4(), daug3->getP4(), daug4->getP4() );
	+ tree->Fill();
	+
	+ root_part->deleteTree();
	+ } while ( count++ < nevent );
	+
	+ file->Write();
	+ file->Close();
	+ EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "SUCCESS\n";
	+}
	diff --git a/test/exampleFiles/4BodyPhsp.DEC b/test/exampleFiles/4BodyPhsp.DEC
	new file mode 100644
	--- /dev/null
	+++ b/test/exampleFiles/4BodyPhsp.DEC
	@@ -0,0 +1,10 @@
	+Decay B+
	+#1.000 mu+ mu- K+ rho0 PHSP;
	+#1.000 mu+ mu- K+ pi0 PHSP;
	+#1.000 mu+ mu- K+ pi0 FOURBODYPHSP;
	+#1.000 mu+ mu- K+ pi0 FOURBODYPHSP 1.0 2.0 1.0 2.0;
	+1.000 mu+ mu- K+ rho0 FOURBODYPHSP;
	+Enddecay
	+CDecay B-
	+
	+End

File Metadata

Mime Type: text/plain
Expires: Mon, Sep 29, 12:51 PM (16 h, 5 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6533766
Default Alt Text: D47.1759146693.diff (30 KB)

D47.1759146693.diffNo OneTemporaryActions

D47.1759146693.diffView Options

File Metadata

Event Timeline

D47.1759146693.diff
No OneTemporary
Actions

D47.1759146693.diff
View Options