diff --git a/EvtGenBase/EvtGenKine.hh b/EvtGenBase/EvtGenKine.hh
--- a/EvtGenBase/EvtGenKine.hh
+++ b/EvtGenBase/EvtGenKine.hh
@@ -22,6 +22,7 @@
#define EVTGENKINE_HH
class EvtVector4R;
+class EvtParticle;
class EvtGenKine {
public:
@@ -30,6 +31,14 @@
static double PhaseSpacePole( double M, double m1, double m2, double m3,
double a, EvtVector4R p4[10] );
+
+ /*
+ * Function which takes two invariant masses squared in 3-body decay and
+ * parent after makeDaughters() and generateMassTree() and
+ * calculates/generates momenta of daughters and sets those.
+ */
+ static void ThreeBodyKine( const double m12Sq, const double m23Sq,
+ EvtParticle* p );
};
#endif
diff --git a/EvtGenModels/EvtFlatSqDalitz.hh b/EvtGenModels/EvtFlatSqDalitz.hh
--- a/EvtGenModels/EvtFlatSqDalitz.hh
+++ b/EvtGenModels/EvtFlatSqDalitz.hh
@@ -40,6 +40,12 @@
void initProbMax() override;
void decay( EvtParticle* p ) override;
+
+ private:
+ double m_mPrimeMin{ 0. };
+ double m_mPrimeMax{ 1. };
+ double m_thetaPrimeMin{ 0. };
+ double m_thetaPrimeMax{ 1. };
};
#endif
diff --git a/src/EvtGenBase/EvtGenKine.cpp b/src/EvtGenBase/EvtGenKine.cpp
--- a/src/EvtGenBase/EvtGenKine.cpp
+++ b/src/EvtGenBase/EvtGenKine.cpp
@@ -25,6 +25,7 @@
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4R.hh"
+#include "EvtGenBase/EvtParticle.hh"
#include <iostream>
#include <math.h>
@@ -340,3 +341,60 @@
return 1.0 + a / ( m12sq * m12sq );
}
+
+/*
+ * Function which takes two invariant masses squared in 3-body decay and
+ * parent after makeDaughters() and generateMassTree() and
+ * calculates/generates momenta of daughters and sets those.
+ */
+void EvtGenKine::ThreeBodyKine( const double m12Sq, const double m23Sq,
+ EvtParticle* p )
+{
+ 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();
+ const double mParentSq{ mParent * mParent };
+ const double mDaug1Sq{ mDaug1 * mDaug1 };
+ const double mDaug2Sq{ mDaug2 * mDaug2 };
+ const double mDaug3Sq{ mDaug3 * mDaug3 };
+ const double invMParent{ 1. / mParent };
+
+ const double En1 = 0.5 * ( mParentSq + mDaug1Sq - m23Sq ) * invMParent;
+ const double En3 = 0.5 * ( mParentSq + mDaug3Sq - m12Sq ) * invMParent;
+ const double En2 = mParent - En1 - En3;
+ const double p1mag = std::sqrt( En1 * En1 - mDaug1Sq );
+ const double p2mag = std::sqrt( En2 * En2 - mDaug2Sq );
+ double cosPhi = 0.5 * ( mDaug1Sq + mDaug2Sq + 2 * En1 * En2 - m12Sq ) /
+ ( p1mag * p2mag );
+
+ double sinPhi = std::sqrt( 1 - cosPhi * cosPhi );
+ if ( EvtRandom::Flat( 0., 1. ) > 0.5 ) {
+ sinPhi *= -1;
+ }
+ const double p2x = p2mag * cosPhi;
+ const double p2y = p2mag * sinPhi;
+ const double p3x = -p1mag - p2x;
+ const double p3y = -p2y;
+
+ // Construct 4-momenta and rotate them randomly in space
+ EvtVector4R p1( En1, p1mag, 0., 0. );
+ EvtVector4R p2( En2, p2x, p2y, 0. );
+ EvtVector4R p3( En3, p3x, p3y, 0. );
+ 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 );
+ p1.applyRotateEuler( euler1, euler2, euler3 );
+ p2.applyRotateEuler( euler1, euler2, euler3 );
+ p3.applyRotateEuler( euler1, euler2, euler3 );
+
+ // set momenta for daughters
+ daug1->init( daug1->getId(), p1 );
+ daug2->init( daug2->getId(), p2 );
+ daug3->init( daug3->getId(), p3 );
+
+ return;
+}
diff --git a/src/EvtGenModels/EvtFlatSqDalitz.cpp b/src/EvtGenModels/EvtFlatSqDalitz.cpp
--- a/src/EvtGenModels/EvtFlatSqDalitz.cpp
+++ b/src/EvtGenModels/EvtFlatSqDalitz.cpp
@@ -20,14 +20,13 @@
#include "EvtGenModels/EvtFlatSqDalitz.hh"
-#include "EvtGenBase/EvtDiracSpinor.hh"
+#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh"
+#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
-#include "EvtGenBase/EvtTensor4C.hh"
-#include "EvtGenBase/EvtVector4C.hh"
#include <fstream>
#include <stdio.h>
@@ -51,72 +50,66 @@
void EvtFlatSqDalitz::initProbMax()
{
- setProbMax( 1. );
+ noProbMax();
}
void EvtFlatSqDalitz::init()
{
- // check that there are 0 arguments
- checkNArg( 0 );
-
//check there are 3 daughters
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 )
{
- p->initializePhaseSpace( getNDaug(), getDaugs() );
-
- double mB = p->mass();
- double m1 = p->getDaug( 0 )->mass();
- 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. );
+ 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();
+ const double mParentSq = mParent * mParent;
+ const double mDaug1Sq = mDaug1 * mDaug1;
+ const double mDaug2Sq = mDaug2 * mDaug2;
+ const double mDaug3Sq = mDaug3 * mDaug3;
+
+ // 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 en1 = ( m12Sq - mDaug2Sq + mDaug1Sq ) / ( 2. * m12 );
+ const double en3 = ( mParentSq - m12Sq - mDaug3Sq ) / ( 2. * m12 );
+
+ const double p1 = std::sqrt( en1 * en1 - mDaug1Sq );
+ const double p3 = std::sqrt( en3 * en3 - mDaug3Sq );
+ const double m13Sq =
+ mDaug1Sq + mDaug3Sq +
+ 2.0 * ( en1 * en3 - p1 * p3 * std::cos( EvtConst::pi * thetaPrime ) );
+ const double m23Sq = mParentSq - m12Sq - m13Sq + mDaug1Sq + mDaug2Sq +
+ mDaug3Sq;
+
+ // Turn m12 and m23 into momenta
+ EvtGenKine::ThreeBodyKine( m12Sq, m23Sq, p );
return;
}
diff --git a/src/EvtGenModels/EvtThreeBodyPhsp.cpp b/src/EvtGenModels/EvtThreeBodyPhsp.cpp
--- a/src/EvtGenModels/EvtThreeBodyPhsp.cpp
+++ b/src/EvtGenModels/EvtThreeBodyPhsp.cpp
@@ -20,14 +20,15 @@
#include "EvtGenModels/EvtThreeBodyPhsp.hh"
-#include "EvtGenBase/EvtParticle.hh"
-#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtConst.hh"
+#include "EvtGenBase/EvtGenKine.hh"
+#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
+#include "EvtGenBase/EvtReport.hh"
-#include <iostream>
#include <algorithm>
#include <cmath>
+#include <iostream>
std::string EvtThreeBodyPhsp::getName()
{
@@ -67,11 +68,11 @@
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();
+ 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();
const double m12borderMin = mDaug1 + mDaug2;
@@ -115,42 +116,7 @@
}
// At this moment we have valid invariant masses squared
- const double En1 = 0.5 * ( mParent * mParent + mDaug1 * mDaug1 - m23Sq ) / mParent;
- const double En3 = 0.5 * ( mParent * mParent + mDaug3 * mDaug3 - m12Sq ) / mParent;
- const double En2 = mParent - En1 - En3;
- const double p1mag = std::sqrt( En1 * En1 - mDaug1 * mDaug1 );
- const double p2mag = std::sqrt( En2 * En2 - mDaug2 * mDaug2 );
- double cosPhi = 0.5 * ( mDaug1 * mDaug1 + mDaug2 * mDaug2 +
- 2 * En1 * En2 - m12Sq )/p1mag/p2mag;
- if ( cosPhi > 1.0 ) {
- EvtGenReport( EVTGEN_WARNING, "EvtThreeBodyPhsp" )
- << " Model: cosine larger than one: " << cosPhi << std::endl;
- cosPhi = 1.0;
- }
- double sinPhi = std::sqrt( 1 - cosPhi * cosPhi );
- if ( EvtRandom::Flat( 0., 1. ) > 0.5 ) {
- sinPhi *= -1;
- }
- const double p2x = p2mag * cosPhi;
- const double p2y = p2mag * sinPhi;
- const double p3x = -p1mag - p2x;
- const double p3y = -p2y;
-
- // Construct 4-momenta and rotate them randomly in space
- EvtVector4R p1( En1, p1mag, 0., 0. );
- EvtVector4R p2( En2, p2x, p2y, 0. );
- EvtVector4R p3( En3, p3x, p3y, 0. );
- 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 );
- p1.applyRotateEuler(euler1, euler2, euler3);
- p2.applyRotateEuler(euler1, euler2, euler3);
- p3.applyRotateEuler(euler1, euler2, euler3);
-
- // set momenta for daughters
- daug1->init( getDaug( 0 ), p1 );
- daug2->init( getDaug( 1 ), p2 );
- daug3->init( getDaug( 2 ), p3 );
+ EvtGenKine::ThreeBodyKine( m12Sq, m23Sq, p );
return;
}
diff --git a/test/do_tests b/test/do_tests
--- a/test/do_tests
+++ b/test/do_tests
@@ -98,3 +98,4 @@
time ./evtgenlhc_test1 baryonic 1000
time ./evtgenlhc_test1 phspdecaytimecut 10000
time ./evtgenlhc_test1 3bodyPhsp 1000000
+time ./evtgenlhc_test1 flatSqDalitz 1000000
diff --git a/test/evtgenlhc_test1.cc b/test/evtgenlhc_test1.cc
--- a/test/evtgenlhc_test1.cc
+++ b/test/evtgenlhc_test1.cc
@@ -157,6 +157,7 @@
TH1F* mom = 0 );
void runBaryonic( int nEvent, EvtGen& myGenerator );
void run3BPhspRegion( int nEvent, EvtGen& myGenerator );
+void runFlatSqDalitz( int nEvent, EvtGen& myGenerator );
int main( int argc, char* argv[] )
{
@@ -531,7 +532,13 @@
if ( !strcmp( argv[1], "3bodyPhsp" ) ) {
int nevent = atoi( argv[2] );
EvtRadCorr::setNeverRadCorr();
- run3BPhspRegion( nevent, myGenerator);
+ run3BPhspRegion( nevent, myGenerator );
+ }
+
+ if ( !strcmp( argv[1], "flatSqDalitz" ) ) {
+ int nevent = atoi( argv[2] );
+ EvtRadCorr::setNeverRadCorr();
+ runFlatSqDalitz( nevent, myGenerator );
}
//*******************************************************
@@ -5753,3 +5760,62 @@
file->Close();
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "SUCCESS\n";
}
+
+void runFlatSqDalitz( int nevent, EvtGen& myGenerator )
+{
+ TFile* file = new TFile( "flatSqDalitz.root", "RECREATE" );
+
+ TH2F* dalitz = new TH2F( "h4", "Dalitz", 50, 0.0, 1.0, 50, 0.0, 1.0 );
+
+ int count = 1;
+
+ char udecay_name[100];
+ strcpy( udecay_name, "exampleFiles/flatSqDalitz.dec" );
+ myGenerator.readUDecay( udecay_name );
+
+ static EvtId B = EvtPDL::getId( std::string( "Lambda_b0" ) );
+
+ do {
+ EvtVector4R pinit( EvtPDL::getMass( B ), 0.0, 0.0, 0.0 );
+
+ EvtParticle* root_part = EvtParticleFactory::particleFactory( B, pinit );
+
+ myGenerator.generateDecay( root_part );
+
+ double mB = root_part->mass();
+ double m1 = root_part->getDaug( 0 )->mass();
+ double m2 = root_part->getDaug( 1 )->mass();
+ double m3 = root_part->getDaug( 2 )->mass();
+ double mBSq{ mB * mB };
+ double m1Sq{ m1 * m1 };
+ double m2Sq{ m2 * m2 };
+ double m3Sq{ m3 * m3 };
+
+ EvtParticle* daug1 = root_part->getDaug( 0 );
+ EvtParticle* daug2 = root_part->getDaug( 1 );
+ EvtParticle* daug3 = root_part->getDaug( 2 );
+ double m12 = ( daug1->getP4() + daug2->getP4() ).mass();
+ double m13 = ( daug1->getP4() + daug3->getP4() ).mass();
+ double m12Sq{ m12 * m12 };
+ double m13Sq{ m13 * m13 };
+
+ double m12norm =
+ 2 * ( ( m12 - ( m1 + m2 ) ) / ( mB - ( m1 + m2 + m3 ) ) ) - 1;
+ double mPrime = acos( m12norm ) / EvtConst::pi;
+ double en1 = ( m12Sq - m2Sq + m1Sq ) / ( 2.0 * m12 );
+ double en3 = ( mBSq - m12Sq - m3Sq ) / ( 2.0 * m12 );
+ double p1 = std::sqrt( en1 * en1 - m1Sq );
+ double p3 = std::sqrt( en3 * en3 - m3Sq );
+ double cosTheta = ( -m13Sq + m1Sq + m3Sq + 2. * en1 * en3 ) /
+ ( 2. * p1 * p3 );
+ double thPrime = acos( cosTheta ) / EvtConst::pi;
+
+ dalitz->Fill( mPrime, thPrime );
+
+ root_part->deleteTree();
+ } while ( count++ < nevent );
+
+ file->Write();
+ file->Close();
+ EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "SUCCESS\n";
+}
diff --git a/test/exampleFiles/flatSqDalitz.dec b/test/exampleFiles/flatSqDalitz.dec
new file mode 100644
--- /dev/null
+++ b/test/exampleFiles/flatSqDalitz.dec
@@ -0,0 +1,17 @@
+Alias MyLambda Lambda0
+Alias MyantiLambda anti-Lambda0
+ChargeConj MyLambda MyantiLambda
+#
+Decay Lambda_b0
+ 1.000 K+ K- MyLambda FLATSQDALITZ;
+Enddecay
+CDecay anti-Lambda_b0
+#
+Decay MyLambda
+ 1.000 p+ pi- PHSP;
+Enddecay
+CDecay MyantiLambda
+#
+End
+
+