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 . *
+***********************************************************************/
+
+#ifndef EVTFOURBODYPHSP_HH
+#define EVTFOURBODYPHSP_HH
+
+#include
+#include
+#include
+
+#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 phspFactor( const double mM, const double m12,
+ const double m34,
+ std::array& daughters ) const;
+ Shape determineBoundaryShape( const double m12Min, const double m12Max,
+ const double m34Max,
+ const double mMother ) const;
+
+ std::pair generatePairMasses(
+ const double m12Min, const double m12Max, const double m34Min,
+ const double m34Max, const double mMother,
+ const EvtFourBodyPhsp::Shape shape ) const;
+ std::pair generateRectangle( const double m12Min,
+ const double m12Max,
+ const double m34Min,
+ const double m34Max ) const;
+ std::pair generateTrapezoid( const double m12Min,
+ const double m12Max,
+ const double m34Min,
+ const double mMother ) const;
+
+ std::array 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
+#include
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 . *
+***********************************************************************/
+
+#include
+
+#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 EvtFourBodyPhsp::phspFactor(
+ const double mM, const double m12, const double m34,
+ std::array& daughters ) const
+{
+ std::array 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 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 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 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
@@ -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
@@ -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