diff --git a/EvtGenBase/EvtAmp.hh b/EvtGenBase/EvtAmp.hh --- a/EvtGenBase/EvtAmp.hh +++ b/EvtGenBase/EvtAmp.hh @@ -52,29 +52,35 @@ EvtAmp& operator=( const EvtAmp& amp ); /** - * sets the amplitudes calculated in the decay objects - */ - void vertex( const EvtComplex& amp ); + * sets the amplitudes calculated in the decay objects + */ + void vertex( const EvtComplex& amp ) { _amp[0] = amp; } /** - * sets the amplitudes calculated in the decay objects - */ - void vertex( int i1, const EvtComplex& amp ); + * sets the amplitudes calculated in the decay objects + */ + void vertex( int i1, const EvtComplex& amp ) { _amp[i1] = amp; } /** - * sets the amplitudes calculated in the decay objects - */ - void vertex( int i1, int i2, const EvtComplex& amp ); + * sets the amplitudes calculated in the decay objects + */ + void vertex( int i1, int i2, const EvtComplex& amp ) + { + _amp[i1 + i2 * _nstate[0]] = amp; + } /** - * sets the amplitudes calculated in the decay objects - */ - void vertex( int i1, int i2, int i3, const EvtComplex& amp ); + * sets the amplitudes calculated in the decay objects + */ + void vertex( int i1, int i2, int i3, const EvtComplex& amp ) + { + _amp[i1 + ( i2 + i3 * _nstate[1] ) * _nstate[0]] = amp; + } /** - * sets the amplitudes calculated in the decay objects - */ - void vertex( int* i1, const EvtComplex& amp ); + * sets the amplitudes calculated in the decay objects + */ + void vertex( int* i1, const EvtComplex& amp ) { setAmp( i1, amp ); } void dump(); @@ -86,8 +92,10 @@ void setNState( int parent_states, int* daug_states ); // the amplitudes - EvtComplex _amp[125]; - + union { + double _dummybuf[125 * 2]; + EvtComplex _amp[125]; + }; // the number of daughters int _ndaug; diff --git a/EvtGenBase/EvtComplex.hh b/EvtGenBase/EvtComplex.hh --- a/EvtGenBase/EvtComplex.hh +++ b/EvtGenBase/EvtComplex.hh @@ -48,6 +48,8 @@ inline friend double imag( const EvtComplex& c ); inline friend EvtComplex exp( const EvtComplex& c ); friend std::ostream& operator<<( std::ostream& s, const EvtComplex& c ); + inline friend EvtComplex operator&( const EvtComplex&, const EvtComplex& ); + inline friend EvtComplex operator|( const EvtComplex&, const EvtComplex& ); public: EvtComplex() : _rpart( 0.0 ), _ipart( 0.0 ) {} @@ -60,8 +62,8 @@ } inline EvtComplex& operator*=( double d ); inline EvtComplex& operator/=( double d ); - EvtComplex& operator*=( EvtComplex c ); - EvtComplex& operator/=( EvtComplex c ); + inline EvtComplex& operator*=( EvtComplex c ); + inline EvtComplex& operator/=( EvtComplex c ); inline EvtComplex& operator=( const EvtComplex& c ); inline EvtComplex& operator+=( const EvtComplex& c ); inline EvtComplex& operator-=( const EvtComplex& c ); @@ -168,6 +170,22 @@ c1._rpart * c2._ipart + c1._ipart * c2._rpart ); } +// conj(c1)*c2 -- don't forget that the operators "&" and "|" has +// priority less than the arithmetic operations -- explicit brackets +// around are necessary! +EvtComplex operator&( const EvtComplex& c1, const EvtComplex& c2 ) +{ + return EvtComplex( c1._rpart * c2._rpart + c1._ipart * c2._ipart, + c1._rpart * c2._ipart - c1._ipart * c2._rpart ); +} + +// c1*conj(c2) +EvtComplex operator|( const EvtComplex& c1, const EvtComplex& c2 ) +{ + return EvtComplex( c1._rpart * c2._rpart + c1._ipart * c2._ipart, + c1._ipart * c2._rpart - c1._rpart * c2._ipart ); +} + EvtComplex operator-( const EvtComplex& c1, const EvtComplex& c2 ) { return EvtComplex( c1._rpart - c2._rpart, c1._ipart - c2._ipart ); @@ -242,4 +260,27 @@ return exp( c._rpart ) * EvtComplex( cos( c._ipart ), sin( c._ipart ) ); } +inline EvtComplex& EvtComplex::operator*=( EvtComplex c ) +{ + double r = _rpart * c._rpart - _ipart * c._ipart; + double i = _rpart * c._ipart + _ipart * c._rpart; + + _rpart = r; + _ipart = i; + + return *this; +} + +inline EvtComplex& EvtComplex::operator/=( EvtComplex c ) +{ + double inv = 1.0 / ( c._rpart * c._rpart + c._ipart * c._ipart ); + + double r = inv * ( _rpart * c._rpart + _ipart * c._ipart ); + double i = inv * ( _ipart * c._rpart - _rpart * c._ipart ); + + _rpart = r; + _ipart = i; + + return *this; +} #endif diff --git a/EvtGenBase/EvtDiracSpinor.hh b/EvtGenBase/EvtDiracSpinor.hh --- a/EvtGenBase/EvtDiracSpinor.hh +++ b/EvtGenBase/EvtDiracSpinor.hh @@ -62,28 +62,34 @@ public: inline EvtDiracSpinor(); - EvtDiracSpinor( const EvtComplex& sp0, const EvtComplex& sp1, - const EvtComplex& sp2, const EvtComplex& sp3 ); + inline EvtDiracSpinor( const EvtComplex& sp0, const EvtComplex& sp1, + const EvtComplex& sp2, const EvtComplex& sp3 ); inline EvtDiracSpinor( const EvtDiracSpinor& dspinor ); inline EvtDiracSpinor& operator=( const EvtDiracSpinor& dspinor ); inline EvtDiracSpinor& operator+=( const EvtDiracSpinor& u2 ); inline EvtDiracSpinor& operator-=( const EvtDiracSpinor& u2 ); - void set( const EvtComplex& sp0, const EvtComplex& sp1, - const EvtComplex& sp2, const EvtComplex& sp3 ); - void set_spinor( int i, const EvtComplex& sp ); - const EvtComplex& get_spinor( int i ) const; + inline void set( const EvtComplex& sp0, const EvtComplex& sp1, + const EvtComplex& sp2, const EvtComplex& sp3 ); + inline void set_spinor( int i, const EvtComplex& sp ); + inline EvtComplex get_spinor( int i ) const; EvtDiracSpinor conj() const; void applyRotateEuler( double alpha, double beta, double gamma ); void applyBoostTo( const EvtVector4R& p4 ); void applyBoostTo( const EvtVector3R& boost ); EvtDiracSpinor adjoint() const; + inline EvtComplex operator[]( int i ) const { return spinor[i]; } + inline EvtComplex& operator[]( int i ) { return spinor[i]; } + private: EvtComplex spinor[4]; }; +void EvtLeptonVandACurrents( EvtVector4C&, EvtVector4C&, const EvtDiracSpinor&, + const EvtDiracSpinor& ); + EvtDiracSpinor::EvtDiracSpinor() { spinor[0] = EvtComplex(); @@ -100,6 +106,31 @@ spinor[3] = dspinor.spinor[3]; } +EvtDiracSpinor::EvtDiracSpinor( const EvtComplex& sp0, const EvtComplex& sp1, + const EvtComplex& sp2, const EvtComplex& sp3 ) +{ + set( sp0, sp1, sp2, sp3 ); +} + +void EvtDiracSpinor::set( const EvtComplex& sp0, const EvtComplex& sp1, + const EvtComplex& sp2, const EvtComplex& sp3 ) +{ + spinor[0] = sp0; + spinor[1] = sp1; + spinor[2] = sp2; + spinor[3] = sp3; +} + +void EvtDiracSpinor::set_spinor( int i, const EvtComplex& sp ) +{ + spinor[i] = sp; +} + +EvtComplex EvtDiracSpinor::get_spinor( int i ) const +{ + return spinor[i]; +} + EvtDiracSpinor& EvtDiracSpinor::operator=( const EvtDiracSpinor& dspinor ) { spinor[0] = dspinor.spinor[0]; diff --git a/EvtGenBase/EvtGenKine.hh b/EvtGenBase/EvtGenKine.hh --- a/EvtGenBase/EvtGenKine.hh +++ b/EvtGenBase/EvtGenKine.hh @@ -23,15 +23,19 @@ class EvtVector4R; class EvtParticle; +class EvtLinSample; class EvtGenKine { public: - static double PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30], - double mp ); + static double PhaseSpace( int ndaug, const double mass[30], + EvtVector4R p4[30], double mp ); static double PhaseSpacePole( double M, double m1, double m2, double m3, double a, EvtVector4R p4[10] ); + static double PhaseSpacePole2( double M, double m1, double m2, double m3, + EvtVector4R p4[10], const EvtLinSample& ); + /* * Function which takes two invariant masses squared in 3-body decay and * parent after makeDaughters() and generateMassTree() and diff --git a/EvtGenBase/EvtLinSample.hh b/EvtGenBase/EvtLinSample.hh new file mode 100644 --- /dev/null +++ b/EvtGenBase/EvtLinSample.hh @@ -0,0 +1,20 @@ +#ifndef EVTLINSAMPLE_HH +#define EVTLINSAMPLE_HH + +#include +struct EvtPointf { + float x, y; +}; + +struct EvtLinSample { + EvtLinSample() {} + EvtLinSample( const char* ); + EvtLinSample( const std::vector& _v ) { init( _v ); } + void init( const std::vector& _v ); + std::pair operator()( double ) const; + + std::vector v; + std::vector I; +}; + +#endif diff --git a/EvtGenBase/EvtMTRandomEngine.hh b/EvtGenBase/EvtMTRandomEngine.hh --- a/EvtGenBase/EvtMTRandomEngine.hh +++ b/EvtGenBase/EvtMTRandomEngine.hh @@ -31,11 +31,10 @@ virtual double random(); - private: - std::mt19937 engine_; + virtual unsigned long urandom(); - typedef std::uniform_real_distribution URDist; - URDist distribution_; + private: + std::mt19937_64 engine_; }; #endif diff --git a/EvtGenBase/EvtParticle.hh b/EvtGenBase/EvtParticle.hh --- a/EvtGenBase/EvtParticle.hh +++ b/EvtGenBase/EvtParticle.hh @@ -365,7 +365,7 @@ /** * Get forward spin density matrix. */ - EvtSpinDensity getSpinDensityForward() { return _rhoForward; } + const EvtSpinDensity& getSpinDensityForward() const { return _rhoForward; } /** * Set backward spin density matrix. diff --git a/EvtGenBase/EvtPdf.hh b/EvtGenBase/EvtPdf.hh --- a/EvtGenBase/EvtPdf.hh +++ b/EvtGenBase/EvtPdf.hh @@ -144,7 +144,7 @@ public: typedef T result_type; - EvtPdfGen() : _pdf( 0 ) {} + EvtPdfGen() : _pdf( nullptr ) {} EvtPdfGen( const EvtPdfGen& other ) : _pdf( other._pdf ? other._pdf->clone() : nullptr ) { diff --git a/EvtGenBase/EvtRandom.hh b/EvtGenBase/EvtRandom.hh --- a/EvtGenBase/EvtRandom.hh +++ b/EvtGenBase/EvtRandom.hh @@ -34,6 +34,8 @@ static double random(); + static unsigned long urandom(); + //This class does not take ownership of the random engine; //the caller needs to make sure that the engine is not //destroyed. diff --git a/EvtGenBase/EvtRandomEngine.hh b/EvtGenBase/EvtRandomEngine.hh --- a/EvtGenBase/EvtRandomEngine.hh +++ b/EvtGenBase/EvtRandomEngine.hh @@ -31,6 +31,8 @@ virtual double random() = 0; + virtual unsigned long urandom() = 0; + private: }; diff --git a/EvtGenBase/EvtSimpleRandomEngine.hh b/EvtGenBase/EvtSimpleRandomEngine.hh --- a/EvtGenBase/EvtSimpleRandomEngine.hh +++ b/EvtGenBase/EvtSimpleRandomEngine.hh @@ -31,6 +31,8 @@ double random() override; + unsigned long urandom() override; + private: unsigned long int _next; }; diff --git a/EvtGenBase/EvtSpinDensity.hh b/EvtGenBase/EvtSpinDensity.hh --- a/EvtGenBase/EvtSpinDensity.hh +++ b/EvtGenBase/EvtSpinDensity.hh @@ -22,29 +22,77 @@ #define EVTSPINDENSITY_HH #include "EvtGenBase/EvtComplex.hh" +#include + // Description: This class holds a spin density matrix, it is // a complex nxn matrix. class EvtSpinDensity { public: + EvtSpinDensity() : dim( 0 ), rho( nullptr ) {} + EvtSpinDensity( int n ) : dim( 0 ), rho( nullptr ) { setDim( n ); } EvtSpinDensity( const EvtSpinDensity& density ); EvtSpinDensity& operator=( const EvtSpinDensity& density ); virtual ~EvtSpinDensity(); - EvtSpinDensity(); void setDim( int n ); - int getDim() const; - void set( int i, int j, const EvtComplex& rhoij ); - const EvtComplex& get( int i, int j ) const; - double normalizedProb( const EvtSpinDensity& d ); + int getDim() const { return dim; } + void set( int i, int j, const EvtComplex& rhoij ) + { + assert( i < dim && j < dim ); + rho[i * dim + j] = rhoij; + } + EvtComplex get( int i, int j ) const + { + assert( i < dim && j < dim ); + return rho[i * dim + j]; + } + EvtComplex operator()( int i, int j ) const + { + assert( i < dim && j < dim ); + return rho[i * dim + j]; + } + EvtComplex& operator()( int i, int j ) + { + assert( i < dim && j < dim ); + return rho[i * dim + j]; + } + EvtSpinDensity& operator-=( const EvtSpinDensity& d ) + { + assert( dim == d.dim ); + for ( int i = 0; i < dim * dim; i++ ) + rho[i] -= d.rho[i]; + return *this; + } + EvtSpinDensity& operator+=( const EvtSpinDensity& d ) + { + assert( dim == d.dim ); + for ( int i = 0; i < dim * dim; i++ ) + rho[i] += d.rho[i]; + return *this; + } + double normalizedProb( const EvtSpinDensity& d ) const; friend std::ostream& operator<<( std::ostream& s, const EvtSpinDensity& d ); void setDiag( int n ); - int check(); private: - EvtComplexPtrPtr rho; int dim; + EvtComplex* rho; }; +inline EvtSpinDensity operator-( const EvtSpinDensity& a, const EvtSpinDensity& b ) +{ + EvtSpinDensity t( a ); + t -= b; + return t; +} + +inline EvtSpinDensity operator+( const EvtSpinDensity& a, const EvtSpinDensity& b ) +{ + EvtSpinDensity t( a ); + t -= b; + return t; +} + #endif diff --git a/EvtGenBase/EvtStreamInputIterator.hh b/EvtGenBase/EvtStreamInputIterator.hh --- a/EvtGenBase/EvtStreamInputIterator.hh +++ b/EvtGenBase/EvtStreamInputIterator.hh @@ -40,7 +40,7 @@ typedef const Point* pointer; typedef const Point& reference; - EvtStreamInputIterator() : _counter( 0 ) {} + EvtStreamInputIterator() : _counter( nullptr ) {} EvtStreamInputIterator( const EvtStreamInputIterator& other ) : _counter( other._counter ? other._counter->clone() : nullptr ), diff --git a/EvtGenBase/EvtTensor4C.hh b/EvtGenBase/EvtTensor4C.hh --- a/EvtGenBase/EvtTensor4C.hh +++ b/EvtGenBase/EvtTensor4C.hh @@ -33,6 +33,11 @@ EvtTensor4C directProd( const EvtVector4R& c1, const EvtVector4R& c2 ); EvtTensor4C directProd( const EvtVector4C& c1, const EvtVector4C& c2 ); EvtTensor4C directProd( const EvtVector4C& c1, const EvtVector4R& c2 ); + EvtTensor4C asymProd( const EvtVector4R&, const EvtVector4R& ); + EvtVector4R asymProd( const EvtVector4R& a, const EvtVector4R& b, + const EvtVector4R& c ); + EvtVector4C asymProd( const EvtVector4C& a, const EvtVector4R& b, + const EvtVector4R& c ); } // namespace EvtGenFunctions class EvtTensor4C final { @@ -42,6 +47,8 @@ const EvtVector4C& c2 ); friend EvtTensor4C EvtGenFunctions::directProd( const EvtVector4C& c1, const EvtVector4R& c2 ); + friend EvtTensor4C EvtGenFunctions::asymProd( const EvtVector4R&, + const EvtVector4R& ); friend EvtTensor4C rotateEuler( const EvtTensor4C& e, double alpha, double beta, double gamma ); @@ -89,6 +96,7 @@ EvtVector4C cont2( const EvtVector4C& v4 ) const; EvtVector4C cont1( const EvtVector4R& v4 ) const; EvtVector4C cont2( const EvtVector4R& v4 ) const; + EvtTensor4C& addScaled( EvtComplex, const EvtTensor4C& ); private: EvtComplex t[4][4]; diff --git a/EvtGenBase/EvtVector4C.hh b/EvtGenBase/EvtVector4C.hh --- a/EvtGenBase/EvtVector4C.hh +++ b/EvtGenBase/EvtVector4C.hh @@ -45,9 +45,9 @@ friend EvtVector4C operator-( const EvtVector4C& v1, const EvtVector4C& v2 ); public: - EvtVector4C(); - EvtVector4C( const EvtComplex&, const EvtComplex&, const EvtComplex&, - const EvtComplex& ); + inline EvtVector4C(){}; + inline EvtVector4C( const EvtComplex&, const EvtComplex&, const EvtComplex&, + const EvtComplex& ); inline void set( int, const EvtComplex& ); inline void set( const EvtComplex&, const EvtComplex&, const EvtComplex&, const EvtComplex& ); @@ -70,6 +70,15 @@ EvtComplex v[4]; }; +inline EvtVector4C::EvtVector4C( const EvtComplex& e0, const EvtComplex& e1, + const EvtComplex& e2, const EvtComplex& e3 ) +{ + v[0] = e0; + v[1] = e1; + v[2] = e2; + v[3] = e3; +} + EvtVector4C rotateEuler( const EvtVector4C& e, double alpha, double beta, double gamma ); EvtVector4C boostTo( const EvtVector4C& e, const EvtVector4R p4 ); diff --git a/EvtGenBase/EvtVectorParticle.hh b/EvtGenBase/EvtVectorParticle.hh --- a/EvtGenBase/EvtVectorParticle.hh +++ b/EvtGenBase/EvtVectorParticle.hh @@ -33,12 +33,12 @@ void init( EvtId part_n, double e, double px, double py, double pz ); void init( EvtId part_n, const EvtVector4R& p ) override; - void init( EvtId part_n, const EvtVector4R& p, const EvtVector4C&, - const EvtVector4C&, const EvtVector4C& ); - EvtVector4C epsParent( int i ) const override - { - return boostTo( _eps[i], this->getP4() ); - } + // void init( EvtId part_n, const EvtVector4R& p, const EvtVector4C&, + // const EvtVector4C&, const EvtVector4C& ); + EvtVector4C epsParent( int i ) const override; + // { + // return boostTo( _eps[i], this->getP4() ); + // } EvtVector4C eps( int i ) const override { return _eps[i]; } EvtSpinDensity rotateToHelicityBasis() const override; EvtSpinDensity rotateToHelicityBasis( double alpha, double beta, diff --git a/EvtGenModels/EvtPi0Dalitz.hh b/EvtGenModels/EvtPi0Dalitz.hh --- a/EvtGenModels/EvtPi0Dalitz.hh +++ b/EvtGenModels/EvtPi0Dalitz.hh @@ -36,7 +36,7 @@ void decay( EvtParticle* p ) override; private: - double m_poleSize{ 0.00000002 }; + double m_poleSize{ 1 }; // Following are rho mass and width, but in order to keep consistency // with what was done before do not use data from particle table. diff --git a/EvtGenModels/EvtPi0Dalitz.hh b/EvtGenModels/EvtbTosllBSZFF.hh copy from EvtGenModels/EvtPi0Dalitz.hh copy to EvtGenModels/EvtbTosllBSZFF.hh --- a/EvtGenModels/EvtPi0Dalitz.hh +++ b/EvtGenModels/EvtbTosllBSZFF.hh @@ -18,30 +18,25 @@ * along with EvtGen. If not, see . * ***********************************************************************/ -#ifndef EVTPI0DALITZ_HH -#define EVTPI0DALITZ_HH +#ifndef EVTBTOSLLBSZFF_HH +#define EVTBTOSLLBSZFF_HH -#include "EvtGenBase/EvtDecayProb.hh" +#include "EvtGenModels/EvtbTosllFF.hh" -class EvtParticle; +class EvtId; -class EvtPi0Dalitz : public EvtDecayProb { - public: - std::string getName() override; - EvtDecayBase* clone() override; - - void init() override; - void initProbMax() override; +// Description: Form factors for b->sll according to Aoife Bharucha, David M. Straub, Roman Zwicky +// https://arxiv.org/abs/1503.05534 - void decay( EvtParticle* p ) override; - - private: - double m_poleSize{ 0.00000002 }; +class EvtbTosllBSZFF : public EvtbTosllFF { + public: + EvtbTosllBSZFF() {} - // Following are rho mass and width, but in order to keep consistency - // with what was done before do not use data from particle table. - const double m_m0Sq{ 0.768 * 0.768 }; - const double m_m0SqG0Sq{ m_m0Sq * 0.151 * 0.151 }; + void getScalarFF( EvtId, EvtId, double, double, double&, double&, + double& ) override; + void getVectorFF( EvtId parent, EvtId vmeson, double q2, double vmesonmass, + double& a1, double& a2, double& a0, double& v, double& t1, + double& t2, double& t3 ) override; }; #endif diff --git a/EvtGenModels/EvtPi0Dalitz.hh b/EvtGenModels/EvtbTosllNP.hh copy from EvtGenModels/EvtPi0Dalitz.hh copy to EvtGenModels/EvtbTosllNP.hh --- a/EvtGenModels/EvtPi0Dalitz.hh +++ b/EvtGenModels/EvtbTosllNP.hh @@ -18,30 +18,32 @@ * along with EvtGen. If not, see . * ***********************************************************************/ -#ifndef EVTPI0DALITZ_HH -#define EVTPI0DALITZ_HH +#ifndef EVTBTOSLLNP_HH +#define EVTBTOSLLNP_HH -#include "EvtGenBase/EvtDecayProb.hh" +#include "EvtGenBase/EvtDecayAmp.hh" +#include + +class EvtbTosllFF; +class EvtbTosllAmp; class EvtParticle; -class EvtPi0Dalitz : public EvtDecayProb { +// Description:Implementation of the b->sll decays according to Rusa and Rahul + +class EvtbTosllNP : public EvtDecayAmp { public: std::string getName() override; EvtDecayBase* clone() override; + void decay( EvtParticle* p ) override; void init() override; void initProbMax() override; - void decay( EvtParticle* p ) override; - private: - double m_poleSize{ 0.00000002 }; - - // Following are rho mass and width, but in order to keep consistency - // with what was done before do not use data from particle table. - const double m_m0Sq{ 0.768 * 0.768 }; - const double m_m0SqG0Sq{ m_m0Sq * 0.151 * 0.151 }; + std::unique_ptr _calcamp; + std::unique_ptr _ffmodel; + double _poleSize; }; #endif diff --git a/EvtGenModels/EvtPi0Dalitz.hh b/EvtGenModels/EvtbTosllNPR.hh copy from EvtGenModels/EvtPi0Dalitz.hh copy to EvtGenModels/EvtbTosllNPR.hh --- a/EvtGenModels/EvtPi0Dalitz.hh +++ b/EvtGenModels/EvtbTosllNPR.hh @@ -18,30 +18,34 @@ * along with EvtGen. If not, see . * ***********************************************************************/ -#ifndef EVTPI0DALITZ_HH -#define EVTPI0DALITZ_HH +#ifndef EVTBTOSLLNPR_HH +#define EVTBTOSLLNPR_HH -#include "EvtGenBase/EvtDecayProb.hh" +#include "EvtGenBase/EvtDecayAmp.hh" +#include "EvtGenBase/EvtLinSample.hh" +#include + +class EvtbTosllFF; +class EvtbTosllAmp; class EvtParticle; -class EvtPi0Dalitz : public EvtDecayProb { +// Description: Implementation of the b->sll decays with resonances according to Rusa and Rahul + +class EvtbTosllNPR : public EvtDecayAmp { public: std::string getName() override; EvtDecayBase* clone() override; + void decay( EvtParticle* p ) override; void init() override; void initProbMax() override; - void decay( EvtParticle* p ) override; - private: - double m_poleSize{ 0.00000002 }; - - // Following are rho mass and width, but in order to keep consistency - // with what was done before do not use data from particle table. - const double m_m0Sq{ 0.768 * 0.768 }; - const double m_m0SqG0Sq{ m_m0Sq * 0.151 * 0.151 }; + std::unique_ptr _calcamp; + std::unique_ptr _ffmodel; + double _poleSize; + EvtLinSample _ls; }; #endif diff --git a/EvtGenModels/EvtbTosllVectorAmpNP.hh b/EvtGenModels/EvtbTosllVectorAmpNP.hh new file mode 100644 --- /dev/null +++ b/EvtGenModels/EvtbTosllVectorAmpNP.hh @@ -0,0 +1,64 @@ + +/*********************************************************************** +* 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 EVTBTOSLLVECTORAMPNP_HH +#define EVTBTOSLLVECTORAMPNP_HH + +#include "EvtGenBase/EvtComplex.hh" + +#include "EvtGenModels/EvtbTosllAmp.hh" + +#include +#include + +class EvtAmp; +class EvtParticle; +class EvtbTosllFF; + +class EvtbTosllVectorAmpNP : public EvtbTosllAmp { + public: + EvtbTosllVectorAmpNP() : + m_dc7( 0 ), + m_dc9( 0 ), + m_dc10( 0 ), + m_c7p( 0 ), + m_c9p( 0 ), + m_c10p( 0 ), + m_cS( 0 ), + m_cP( 0 ) + { + } + void CalcAmp( EvtParticle* parent, EvtAmp& amp, + EvtbTosllFF* formFactors ) override; + void CalcSAmp( EvtParticle* parent, EvtAmp& amp, EvtbTosllFF* formFactors ); + void CalcVAmp( EvtParticle* parent, EvtAmp& amp, EvtbTosllFF* formFactors ); + + EvtComplex m_dc7; // delta C_7eff -- addition to NNLO SM value + EvtComplex m_dc9; // delta C_9eff -- addition to NNLO SM value + EvtComplex m_dc10; // delta C_10eff -- addition to NNLO SM value + EvtComplex m_c7p; // C'_7eff -- right hand polarizations + EvtComplex m_c9p; // C'_9eff -- right hand polarizations + EvtComplex m_c10p; // c'_10eff -- right hand polarizations + EvtComplex m_cS; // (C_S - C'_S) -- scalar right and left polarizations + EvtComplex m_cP; // (C_P - C'_P) -- pseudo-scalar right and left polarizations + std::vector>> m_reso; // resonance contribution +}; + +#endif diff --git a/src/EvtGenBase/EvtAmp.cpp b/src/EvtGenBase/EvtAmp.cpp --- a/src/EvtGenBase/EvtAmp.cpp +++ b/src/EvtGenBase/EvtAmp.cpp @@ -141,64 +141,31 @@ EvtSpinDensity EvtAmp::getSpinDensity() { + int np = _pstates; EvtSpinDensity rho; - rho.setDim( _pstates ); - - EvtComplex temp; - - int i, j, n; - - if ( _pstates == 1 ) { - if ( _nontrivial == 0 ) { - rho.set( 0, 0, _amp[0] * conj( _amp[0] ) ); - return rho; - } - - n = 1; - - temp = EvtComplex( 0.0 ); - - for ( i = 0; i < _nontrivial; i++ ) { + rho.setDim( np ); + if ( np == 1 ) { + double a = abs2( _amp[0] ); + int n = 1; + for ( int i = 0; i < _nontrivial; i++ ) n *= _nstate[i]; - } - - for ( i = 0; i < n; i++ ) { - temp += _amp[i] * conj( _amp[i] ); - } - - rho.set( 0, 0, temp ); - ; - - return rho; - - } - - else { - for ( i = 0; i < _pstates; i++ ) { - for ( j = 0; j < _pstates; j++ ) { - temp = EvtComplex( 0.0 ); - - int kk; - - int allloop = 1; - for ( kk = 0; kk < _ndaug; kk++ ) { - allloop *= dstates[kk]; - } - - for ( kk = 0; kk < allloop; kk++ ) { - temp += _amp[_pstates * kk + i] * - conj( _amp[_pstates * kk + j] ); - } - - // if (_nontrivial>3){ - //EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Can't handle so many states in EvtAmp!"< 0 ) { - _spinorRest[0].set( EvtComplex( sqrt( 2.0 * mass() ), 0.0 ), - EvtComplex( 0.0, 0.0 ), EvtComplex( 0.0, 0.0 ), - EvtComplex( 0.0, 0.0 ) ); - _spinorRest[1].set( EvtComplex( 0.0, 0.0 ), - EvtComplex( sqrt( 2.0 * mass() ), 0.0 ), - EvtComplex( 0.0, 0.0 ), EvtComplex( 0.0, 0.0 ) ); - - _spinorParent[0] = boostTo( _spinorRest[0], p4 ); - _spinorParent[1] = boostTo( _spinorRest[1], p4 ); - + _spinorRest[0].set( l, o, o, o ); + _spinorRest[1].set( o, l, o, o ); + _spinorParent[0].set( t, 0, uz, uxy ); + _spinorParent[1].set( 0, t, cxy, -uz ); } else { - _spinorRest[0].set( EvtComplex( 0.0, 0.0 ), EvtComplex( 0.0, 0.0 ), - EvtComplex( sqrt( 2.0 * mass() ), 0.0 ), - EvtComplex( 0.0, 0.0 ) ); - _spinorRest[1].set( EvtComplex( 0.0, 0.0 ), EvtComplex( 0.0, 0.0 ), - EvtComplex( 0.0, 0.0 ), - EvtComplex( sqrt( 2.0 * mass() ), 0.0 ) ); - - _spinorParent[0] = boostTo( _spinorRest[0], p4 ); - _spinorParent[1] = boostTo( _spinorRest[1], p4 ); + _spinorRest[0].set( o, o, l, o ); + _spinorRest[1].set( o, o, o, l ); + _spinorParent[0].set( uz, uxy, t, 0 ); + _spinorParent[1].set( cxy, -uz, 0, t ); } + // EvtDiracSpinor s0 = boostTo( _spinorRest[0], p4 ); + // EvtDiracSpinor s1 = boostTo( _spinorRest[1], p4 ); + // std::cout<< EvtPDL::getStdHep( part_n )<<" "<<_spinorParent[0] - s0<<" "<<_spinorParent[1] - s1 < using std::ostream; -EvtDiracSpinor::EvtDiracSpinor( const EvtComplex& sp0, const EvtComplex& sp1, - const EvtComplex& sp2, const EvtComplex& sp3 ) -{ - set( sp0, sp1, sp2, sp3 ); -} - -void EvtDiracSpinor::set( const EvtComplex& sp0, const EvtComplex& sp1, - const EvtComplex& sp2, const EvtComplex& sp3 ) -{ - spinor[0] = sp0; - spinor[1] = sp1; - spinor[2] = sp2; - spinor[3] = sp3; -} - -void EvtDiracSpinor::set_spinor( int i, const EvtComplex& sp ) -{ - spinor[i] = sp; -} - ostream& operator<<( ostream& s, const EvtDiracSpinor& sp ) { s << "[" << sp.spinor[0] << "," << sp.spinor[1] << "," << sp.spinor[2] @@ -58,11 +38,6 @@ return s; } -const EvtComplex& EvtDiracSpinor::get_spinor( int i ) const -{ - return spinor[i]; -} - EvtDiracSpinor rotateEuler( const EvtDiracSpinor& sp, double alpha, double beta, double gamma ) { @@ -87,6 +62,7 @@ void EvtDiracSpinor::applyBoostTo( const EvtVector4R& p4 ) { +#if 0 double e = p4.get( 0 ); EvtVector3R boost( p4.get( 1 ) / e, p4.get( 2 ) / e, p4.get( 3 ) / e ); @@ -94,10 +70,25 @@ applyBoostTo( boost ); return; +#else + double px = p4.get( 1 ), py = p4.get( 2 ), pz = p4.get( 3 ), e = p4.get( 0 ); + double p2 = px * px + py * py + pz * pz, m2 = e * e - p2, m = sqrt( m2 ), + em = e + m; + double f2 = sqrt( 0.5 / ( em * m ) ), f1 = f2 * em; + double bx = f2 * px, by = f2 * py, bz = f2 * pz; + EvtComplex t = spinor[0], x = spinor[1], y = spinor[2], z = spinor[3], + bxy( bx, by ); + + spinor[0] = f1 * t + ( bxy & z ) + bz * y; + spinor[1] = f1 * x + ( bxy * y ) - bz * z; + spinor[2] = f1 * y + ( bxy & x ) + bz * t; + spinor[3] = f1 * z + ( bxy * t ) - bz * x; +#endif } void EvtDiracSpinor::applyBoostTo( const EvtVector3R& boost ) { +#if 0 double bx, by, bz, gamma, b2, f1, f2; EvtComplex spinorp[4]; @@ -135,6 +126,22 @@ spinor[3] = spinorp[3]; return; +#else + double bx = boost.get( 0 ), by = boost.get( 1 ), bz = boost.get( 2 ); + double b2 = bx * bx + by * by + bz * bz, c2 = 1 - b2, c = sqrt( c2 ); + double f2 = sqrt( 0.5 / ( c2 + c ) ); + bx *= f2; + by *= f2; + bz *= f2; + double f1 = f2 + f2 * c; + EvtComplex t = spinor[0], x = spinor[1], y = spinor[2], z = spinor[3], + bxy( bx, by ); + + spinor[0] = f1 * t + ( bxy & z ) + bz * y; + spinor[1] = f1 * x + ( bxy * y ) - bz * z; + spinor[2] = f1 * y + ( bxy & x ) + bz * t; + spinor[3] = f1 * z + ( bxy * t ) - bz * x; +#endif } void EvtDiracSpinor::applyRotateEuler( double alpha, double beta, double gamma ) @@ -178,175 +185,385 @@ EvtVector4C EvtLeptonVACurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) { - //Old code; below is a new specialized code that does it more efficiently. - //EvtGammaMatrix mat; - //EvtVector4C temp; - //mat.va0(); - //temp.set(0,d*(mat*dp)); - //mat.va1(); - //temp.set(1,d*(mat*dp)); - //mat.va2(); - //temp.set(2,d*(mat*dp)); - //mat.va3(); - //temp.set(3,d*(mat*dp)); - //return temp; - - EvtComplex u02 = ::conj( d.spinor[0] - d.spinor[2] ); - EvtComplex u13 = ::conj( d.spinor[1] - d.spinor[3] ); - - EvtComplex v02 = dp.spinor[0] - dp.spinor[2]; - EvtComplex v13 = dp.spinor[1] - dp.spinor[3]; + // Old code; below is a new specialized code that does it more efficiently. + // EvtVector4C tmp; + // tmp.set(0,d*(EvtGammaMatrix::va0()*dp)); + // tmp.set(1,d*(EvtGammaMatrix::va1()*dp)); + // tmp.set(2,d*(EvtGammaMatrix::va2()*dp)); + // tmp.set(3,d*(EvtGammaMatrix::va3()*dp)); + // return tmp; + + EvtComplex u02( real( d[0] ) - real( d[2] ), imag( d[2] ) - imag( d[0] ) ); + EvtComplex u13( real( d[1] ) - real( d[3] ), imag( d[3] ) - imag( d[1] ) ); + + EvtComplex v02 = dp[2] - dp[0]; + EvtComplex v13 = dp[1] - dp[3]; EvtComplex a = u02 * v02; EvtComplex b = u13 * v13; EvtComplex c = u02 * v13; EvtComplex e = u13 * v02; - - return EvtVector4C( a + b, -( c + e ), EvtComplex( 0, 1 ) * ( c - e ), b - a ); + EvtComplex ce = c + e; + EvtVector4C res( b - a, e - c, EvtComplex( -imag( ce ), real( ce ) ), b + a ); + // std::cout<< res - tmp <conj(); // first conjugate, then multiply with gamma0 - EvtGammaMatrix g0 = EvtGammaMatrix::g0(); - EvtDiracSpinor result; // automatically initialized to 0 - - for ( int i = 0; i < 4; ++i ) - for ( int j = 0; j < 4; ++j ) - result.spinor[i] += d.spinor[j] * g0._gamma[i][j]; - - return result; + // EvtDiracSpinor d = this->conj(); // first conjugate, then multiply with gamma0 + // EvtGammaMatrix g0 = EvtGammaMatrix::g0(); + // EvtDiracSpinor result; // automatically initialized to 0 + + // for ( int i = 0; i < 4; ++i ) + // for ( int j = 0; j < 4; ++j ) + // result.spinor[i] += d.spinor[j] * g0._gamma[i][j]; + + // return result; + const EvtDiracSpinor& d = *this; + double a0r = real( d[0] ), a0i = imag( d[0] ); + double a1r = real( d[1] ), a1i = imag( d[1] ); + double a2r = real( d[2] ), a2i = imag( d[2] ); + double a3r = real( d[3] ), a3i = imag( d[3] ); + return EvtDiracSpinor( EvtComplex( a0r, -a0i ), EvtComplex( a1r, -a1i ), + EvtComplex( -a2r, a2i ), EvtComplex( -a3r, a3i ) ); } EvtComplex operator*( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) { - int i; - EvtComplex temp; + // int i; + // EvtComplex temp; - temp = EvtComplex( 0.0, 0.0 ); + // temp = EvtComplex( 0.0, 0.0 ); - for ( i = 0; i < 4; i++ ) { - temp += conj( d.get_spinor( i ) ) * dp.get_spinor( i ); - } - return temp; + // for ( i = 0; i < 4; i++ ) { + // temp += conj( d.get_spinor( i ) ) * dp.get_spinor( i ); + // } + // return temp; + + double a0r = real( d[0] ), a0i = imag( d[0] ); + double a1r = real( d[1] ), a1i = imag( d[1] ); + double a2r = real( d[2] ), a2i = imag( d[2] ); + double a3r = real( d[3] ), a3i = imag( d[3] ); + + double b0r = real( dp[0] ), b0i = imag( dp[0] ); + double b1r = real( dp[1] ), b1i = imag( dp[1] ); + double b2r = real( dp[2] ), b2i = imag( dp[2] ); + double b3r = real( dp[3] ), b3i = imag( dp[3] ); + + double w0r = a0r * b0r + a0i * b0i, w0i = a0r * b0i - a0i * b0r; + double w1r = a1r * b1r + a1i * b1i, w1i = a1r * b1i - a1i * b1r; + double w2r = a2r * b2r + a2i * b2i, w2i = a2r * b2i - a2i * b2r; + double w3r = a3r * b3r + a3i * b3i, w3i = a3r * b3i - a3i * b3r; + + return EvtComplex( w0r + w1r + ( w2r + w3r ), w0i + w1i + ( w2i + w3i ) ); } diff --git a/src/EvtGenBase/EvtGenKine.cpp b/src/EvtGenBase/EvtGenKine.cpp --- a/src/EvtGenBase/EvtGenKine.cpp +++ b/src/EvtGenBase/EvtGenKine.cpp @@ -21,324 +21,912 @@ #include "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtConst.hh" +#include "EvtGenBase/EvtLinSample.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtVector4R.hh" +#include #include -#include using std::endl; -double EvtPawt( double a, double b, double c ) +inline void orderdouble( double* x, int i, int j ) { - double temp = ( a * a - ( b + c ) * ( b + c ) ) * - ( a * a - ( b - c ) * ( b - c ) ); + auto t = std::min( x[i], x[j] ); + x[j] = std::max( x[i], x[j] ); + x[i] = t; +} - if ( temp <= 0 ) { - return 0.0; - } +#define Od( i, j ) orderdouble( x, i, j ) - return sqrt( temp ) / ( 2.0 * a ); +void sort15( double* x ) +{ + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); + Od( 6, 7 ); + Od( 8, 9 ); + Od( 10, 11 ); + Od( 12, 13 ); + Od( 0, 2 ); + Od( 4, 6 ); + Od( 8, 10 ); + Od( 12, 14 ); + Od( 1, 3 ); + Od( 5, 7 ); + Od( 9, 11 ); + Od( 0, 4 ); + Od( 8, 12 ); + Od( 1, 5 ); + Od( 9, 13 ); + Od( 2, 6 ); + Od( 10, 14 ); + Od( 3, 7 ); + Od( 0, 8 ); + Od( 1, 9 ); + Od( 2, 10 ); + Od( 3, 11 ); + Od( 4, 12 ); + Od( 5, 13 ); + Od( 6, 14 ); + Od( 5, 10 ); + Od( 6, 9 ); + Od( 3, 12 ); + Od( 13, 14 ); + Od( 7, 11 ); + Od( 1, 2 ); + Od( 4, 8 ); + Od( 1, 4 ); + Od( 7, 13 ); + Od( 2, 8 ); + Od( 11, 14 ); + Od( 2, 4 ); + Od( 5, 6 ); + Od( 9, 10 ); + Od( 11, 13 ); + Od( 3, 8 ); + Od( 7, 12 ); + Od( 6, 8 ); + Od( 10, 12 ); + Od( 3, 5 ); + Od( 7, 9 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 7, 8 ); + Od( 9, 10 ); + Od( 11, 12 ); + Od( 6, 7 ); + Od( 8, 9 ); } -double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30], - double mp ) +void sort14( double* x ) +{ + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); + Od( 6, 7 ); + Od( 8, 9 ); + Od( 10, 11 ); + Od( 12, 13 ); + Od( 0, 2 ); + Od( 4, 6 ); + Od( 8, 10 ); + Od( 1, 3 ); + Od( 5, 7 ); + Od( 9, 11 ); + Od( 0, 4 ); + Od( 8, 12 ); + Od( 1, 5 ); + Od( 9, 13 ); + Od( 2, 6 ); + Od( 3, 7 ); + Od( 0, 8 ); + Od( 1, 9 ); + Od( 2, 10 ); + Od( 3, 11 ); + Od( 4, 12 ); + Od( 5, 13 ); + Od( 5, 10 ); + Od( 6, 9 ); + Od( 3, 12 ); + Od( 7, 11 ); + Od( 1, 2 ); + Od( 4, 8 ); + Od( 1, 4 ); + Od( 7, 13 ); + Od( 2, 8 ); + Od( 2, 4 ); + Od( 5, 6 ); + Od( 9, 10 ); + Od( 11, 13 ); + Od( 3, 8 ); + Od( 7, 12 ); + Od( 6, 8 ); + Od( 10, 12 ); + Od( 3, 5 ); + Od( 7, 9 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 7, 8 ); + Od( 9, 10 ); + Od( 11, 12 ); + Od( 6, 7 ); + Od( 8, 9 ); +} -// N body phase space routine. Send parent with -// daughters already defined ( Number and masses ) -// Returns four vectors in parent frame. +void sort13( double* x ) +{ + Od( 1, 3 ); + Od( 5, 11 ); + Od( 4, 7 ); + Od( 8, 9 ); + Od( 0, 12 ); + Od( 6, 10 ); + Od( 1, 4 ); + Od( 6, 8 ); + Od( 9, 10 ); + Od( 11, 12 ); + Od( 0, 5 ); + Od( 3, 7 ); + Od( 1, 2 ); + Od( 5, 9 ); + Od( 10, 12 ); + Od( 0, 6 ); + Od( 8, 11 ); + Od( 0, 1 ); + Od( 4, 5 ); + Od( 3, 8 ); + Od( 2, 6 ); + Od( 1, 2 ); + Od( 3, 4 ); + Od( 6, 7 ); + Od( 8, 9 ); + Od( 5, 11 ); + Od( 1, 3 ); + Od( 6, 10 ); + Od( 2, 4 ); + Od( 5, 8 ); + Od( 9, 11 ); + Od( 7, 12 ); + Od( 2, 3 ); + Od( 4, 6 ); + Od( 7, 10 ); + Od( 4, 5 ); + Od( 7, 9 ); + Od( 10, 11 ); + Od( 6, 8 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 7, 8 ); + Od( 9, 10 ); + Od( 6, 7 ); + Od( 8, 9 ); +} +void sort12( double* x ) { - double energy, p3, alpha, beta; + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); + Od( 6, 7 ); + Od( 8, 9 ); + Od( 10, 11 ); + Od( 1, 3 ); + Od( 5, 7 ); + Od( 9, 11 ); + Od( 0, 2 ); + Od( 4, 6 ); + Od( 8, 10 ); + Od( 5, 9 ); + Od( 6, 10 ); + Od( 1, 2 ); + Od( 7, 11 ); + Od( 0, 4 ); + Od( 1, 5 ); + Od( 9, 10 ); + Od( 2, 6 ); + Od( 3, 7 ); + Od( 4, 8 ); + Od( 5, 9 ); + Od( 1, 2 ); + Od( 6, 10 ); + Od( 0, 4 ); + Od( 7, 11 ); + Od( 3, 8 ); + Od( 1, 4 ); + Od( 7, 10 ); + Od( 2, 3 ); + Od( 5, 6 ); + Od( 8, 9 ); + Od( 3, 5 ); + Od( 6, 8 ); + Od( 2, 4 ); + Od( 7, 9 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 7, 8 ); +} - if ( ndaug == 1 ) { - p4[0].set( mass[0], 0.0, 0.0, 0.0 ); - return 1.0; - } +void sort11( double* x ) +{ + Od( 0, 9 ); + Od( 1, 8 ); + Od( 2, 7 ); + Od( 3, 6 ); + Od( 4, 5 ); + Od( 0, 3 ); + Od( 4, 10 ); + Od( 1, 2 ); + Od( 6, 9 ); + Od( 7, 8 ); + Od( 0, 1 ); + Od( 2, 3 ); + Od( 5, 8 ); + Od( 9, 10 ); + Od( 6, 7 ); + Od( 1, 2 ); + Od( 4, 6 ); + Od( 8, 10 ); + Od( 5, 9 ); + Od( 0, 4 ); + Od( 7, 8 ); + Od( 1, 5 ); + Od( 2, 9 ); + Od( 3, 6 ); + Od( 1, 4 ); + Od( 5, 7 ); + Od( 2, 3 ); + Od( 6, 9 ); + Od( 2, 4 ); + Od( 6, 7 ); + Od( 8, 9 ); + Od( 3, 5 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 7, 8 ); +} - if ( ndaug == 2 ) { - //Two body phase space +void sort10( double* x ) +{ + Od( 1, 8 ); + Od( 0, 4 ); + Od( 5, 9 ); + Od( 3, 7 ); + Od( 2, 6 ); + Od( 0, 3 ); + Od( 6, 9 ); + Od( 4, 7 ); + Od( 0, 1 ); + Od( 3, 6 ); + Od( 8, 9 ); + Od( 2, 5 ); + Od( 1, 5 ); + Od( 7, 9 ); + Od( 0, 2 ); + Od( 4, 8 ); + Od( 1, 2 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 7, 8 ); + Od( 1, 3 ); + Od( 6, 8 ); + Od( 2, 4 ); + Od( 5, 7 ); + Od( 2, 3 ); + Od( 6, 7 ); + Od( 4, 6 ); + Od( 3, 5 ); + Od( 4, 5 ); +} - energy = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) / - ( 2.0 * mp ); +void sort9( double* x ) +{ + Od( 1, 2 ); + Od( 4, 5 ); + Od( 7, 8 ); + Od( 0, 1 ); + Od( 3, 4 ); + Od( 6, 7 ); + Od( 1, 2 ); + Od( 4, 5 ); + Od( 7, 8 ); + Od( 3, 6 ); + Od( 0, 3 ); + Od( 5, 8 ); + Od( 4, 7 ); + Od( 3, 6 ); + Od( 2, 5 ); + Od( 1, 4 ); + Od( 1, 3 ); + Od( 5, 8 ); + Od( 4, 7 ); + Od( 2, 6 ); + Od( 5, 7 ); + Od( 2, 3 ); + Od( 4, 6 ); + Od( 3, 4 ); + Od( 5, 6 ); +} - p3 = 0.0; - if ( energy > mass[0] ) { - p3 = sqrt( energy * energy - mass[0] * mass[0] ); - } +void sort8( double* x ) +{ + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); + Od( 6, 7 ); + Od( 1, 3 ); + Od( 5, 7 ); + Od( 0, 2 ); + Od( 4, 6 ); + Od( 1, 2 ); + Od( 5, 6 ); + Od( 3, 7 ); + Od( 2, 6 ); + Od( 1, 5 ); + Od( 0, 4 ); + Od( 3, 5 ); + Od( 2, 4 ); + Od( 1, 2 ); + Od( 3, 4 ); + Od( 5, 6 ); +} - p4[0].set( energy, 0.0, 0.0, p3 ); +void sort7( double* x ) +{ + Od( 1, 2 ); + Od( 3, 4 ); + Od( 5, 6 ); + Od( 0, 2 ); + Od( 4, 6 ); + Od( 3, 5 ); + Od( 2, 6 ); + Od( 1, 5 ); + Od( 0, 4 ); + Od( 2, 5 ); + Od( 0, 3 ); + Od( 2, 4 ); + Od( 1, 3 ); + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); +} - energy = mp - energy; - p3 = -1.0 * p3; - p4[1].set( energy, 0.0, 0.0, p3 ); +void sort6( double* x ) +{ + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); + Od( 0, 2 ); + Od( 3, 5 ); + Od( 1, 4 ); + Od( 0, 1 ); + Od( 2, 3 ); + Od( 4, 5 ); + Od( 1, 2 ); + Od( 3, 4 ); + Od( 2, 3 ); +} - //Now rotate four vectors. +void sort5( double* x ) +{ + Od( 1, 2 ); + Od( 3, 4 ); + Od( 1, 3 ); + Od( 0, 2 ); + Od( 2, 4 ); + Od( 0, 3 ); + Od( 0, 1 ); + Od( 2, 3 ); + Od( 1, 2 ); +} - alpha = EvtRandom::Flat( EvtConst::twoPi ); - beta = acos( EvtRandom::Flat( -1.0, 1.0 ) ); +void sort4( double* x ) +{ + Od( 0, 1 ); + Od( 2, 3 ); + Od( 0, 2 ); + Od( 1, 3 ); + Od( 1, 2 ); +} - p4[0].applyRotateEuler( alpha, beta, -alpha ); - p4[1].applyRotateEuler( alpha, beta, -alpha ); +void sort3( double* x ) +{ + Od( 0, 1 ); + Od( 1, 2 ); + Od( 0, 1 ); +} - return 1.0; - } +void sort2( double* x ) +{ + Od( 0, 1 ); +} - if ( ndaug != 2 ) { - double wtmax = 0.0; - double pm[5][30], pmin, pmax, psum, rnd[30]; - double ran, wt, pa, costh, sinth, phi, p[4][30], be[4], bep, temp; - int i, il, ilr, i1, il1u, il1, il2r, ilu; - int il2 = 0; +void sort1( double* ) +{ +} - for ( i = 0; i < ndaug; i++ ) { - pm[4][i] = 0.0; - rnd[i] = 0.0; - } +typedef void( fun_t )( double* ); +fun_t* sortfuns[] = { sort1, sort2, sort3, sort4, sort5, + sort6, sort7, sort8, sort9, sort10, + sort11, sort12, sort13, sort14, sort15 }; - pm[0][0] = mp; - pm[1][0] = 0.0; - pm[2][0] = 0.0; - pm[3][0] = 0.0; - pm[4][0] = mp; +inline double sqr( double x ) +{ + return x * x; +} - psum = 0.0; - for ( i = 1; i < ndaug + 1; i++ ) { - psum = psum + mass[i - 1]; - } +// return s = sin and c = cos of phi = k/2^32*2*M_PI +// tested that |s|<=1 and |c|<=1 +// max |s^2 + c^2 - 1| <= 4.440892e-16 +// max relative difference |s-s_true|/s_true < 1.716228e-15 +void __attribute__( ( noinline ) ) +usincos( unsigned long kw, double& s, double& c ) +{ + const static double st[] = { 0, + 4.90676743274180142550e-2, + 9.80171403295606019942e-2, + 1.46730474455361751659e-1, + 1.95090322016128267848e-1, + 2.42980179903263889948e-1, + 2.90284677254462367636e-1, + 3.36889853392220050689e-1, + 3.82683432365089771728e-1, + 4.27555093430282094321e-1, + 4.71396736825997648556e-1, + 5.14102744193221726594e-1, + 5.55570233019602224743e-1, + 5.95699304492433343467e-1, + 6.34393284163645498215e-1, + 6.71558954847018400625e-1, + 7.07106781186547524401e-1, + 7.40951125354959091176e-1, + 7.73010453362736960811e-1, + 8.03207531480644909807e-1, + 8.31469612302545237079e-1, + 8.57728610000272069902e-1, + 8.81921264348355029713e-1, + 9.03989293123443331586e-1, + 9.23879532511286756128e-1, + 9.41544065183020778413e-1, + 9.56940335732208864936e-1, + 9.70031253194543992604e-1, + 9.80785280403230449126e-1, + 9.89176509964780973452e-1, + 9.95184726672196886245e-1, + 9.98795456205172392715e-1, + 1 }; + double x = (long)( kw << ( 5 + 2 ) ), x2 = x * x; + static const double as[3] = { 2.661032484442617284e-21, + -3.140503474026838861e-63, + 1.111886075860967104e-105 }, + ac[3] = { -3.540546941629423600e-42, + 2.089245440416171469e-84, + -4.931274944723895814e-127 }; + unsigned k = (unsigned long)kw >> 32; + int jk = ( k + ( 1 << 24 ) ) << 1, mask = jk >> 31, + absj = ( ( jk >> 26 ) ^ mask ) - mask; + double s0 = st[absj], c0 = st[32 - absj]; + static const double sign[2] = { 1, -1 }; + s0 *= sign[( k + 0 ) >> 31]; + c0 *= sign[( k + ( 1 << 30 ) ) >> 31]; + double sn = ( as[0] + x2 * ( as[1] + x2 * ( as[2] ) ) ); + double dcs = x * ( ac[0] + x2 * ( ac[1] + x2 * ( ac[2] ) ) ); + s = s0 + x * ( sn * c0 + dcs * s0 ); + c = c0 + x * ( dcs * c0 - sn * s0 ); +} - pm[4][ndaug - 1] = mass[ndaug - 1]; +// return s = sin and c = cos of phi = r*2*M_PI +// 0 <= r < 1 +// minimal nonzero r is supposed to be 2^-53 +void rsincos( double r, double& s, double& c ) +{ + long kw = r * 9007199254740992ul; + usincos( kw << 11, s, c ); +} - switch ( ndaug ) { - case 1: - wtmax = 1.0 / 16.0; - break; - case 2: - wtmax = 1.0 / 150.0; - break; - case 3: - wtmax = 1.0 / 2.0; - break; - case 4: - wtmax = 1.0 / 5.0; - break; - case 5: - wtmax = 1.0 / 15.0; - break; - case 6: - wtmax = 1.0 / 15.0; - break; - case 7: - wtmax = 1.0 / 15.0; - break; - case 8: - wtmax = 1.0 / 15.0; - break; - case 9: - wtmax = 1.0 / 15.0; - break; - case 10: - wtmax = 1.0 / 15.0; - break; - case 11: - wtmax = 1.0 / 15.0; - break; - case 12: - wtmax = 1.0 / 15.0; - break; - case 13: - wtmax = 1.0 / 15.0; - break; - case 14: - wtmax = 1.0 / 15.0; - break; - case 15: - wtmax = 1.0 / 15.0; - break; - default: - EvtGenReport( EVTGEN_ERROR, "EvtGen" ) - << "too many daughters for phase space..." << ndaug << " " - << mp << endl; - ; - break; - } +double EvtGenKine::PhaseSpace( int ndaug, const double mass[30], + EvtVector4R p4[30], double mp ) - pmax = mp - psum + mass[ndaug - 1]; +// N body phase space routine. Send parent with +// daughters already defined ( Number and masses ) +// Returns four vectors in parent frame. - pmin = 0.0; +{ + // Maximal weight scale -- based on numerical experiments + const static double wtscale[] = { + 1.000000e+00, 1.000000e+00, 5.000707e-01, 1.924501e-01, 6.249988e-02, + 1.789312e-02, 4.632982e-03, 1.104825e-03, 2.450442e-04, 5.095439e-05, + 1.022054e-05, 1.834890e-06, 3.216923e-07, 5.540728e-08, 9.468573e-09 }; - for ( ilr = 2; ilr < ndaug + 1; ilr++ ) { - il = ndaug + 1 - ilr; - pmax = pmax + mass[il - 1]; - pmin = pmin + mass[il + 1 - 1]; - wtmax = wtmax * EvtPawt( pmax, pmin, mass[il - 1] ); + if ( ndaug == 1 ) { + p4[0].set( mass[0], 0.0, 0.0, 0.0 ); + } else if ( ndaug == 2 ) { + //Two body phase space + double en = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) / + ( 2.0 * mp ); + double p3 = ( en > mass[0] ) ? sqrt( en * en - mass[0] * mass[0] ) : 0; + + //Now uniformly distribute over sphere + double s, c; + usincos( EvtRandom::urandom(), s, c ); + double z = EvtRandom::random() * 2 - 1, r = sqrt( 1 - z * z ); + double pt = p3 * r, px = pt * c, py = pt * s, pz = p3 * z; + p4[0].set( en, px, py, pz ); + p4[1].set( mp - en, -px, -py, -pz ); + } else if ( ndaug == 3 ) { +#define ipck( i, j, k ) ( i | ( j << 2 ) | ( k << 4 ) ) + const static unsigned char indx[9] = { ipck( 2, 1, 0 ), + ipck( 2, 0, 1 ), + 0, + ipck( 1, 2, 0 ), + 0, + ipck( 0, 2, 1 ), + 0, + ipck( 1, 0, 2 ), + ipck( 0, 1, 2 ) }; +#undef ipck + double m0 = mass[0], m1 = mass[1], m2 = mass[2]; + unsigned i0 = ( m0 > m1 ) + ( m0 > m2 ); + unsigned i1 = ( m1 >= m0 ) + ( m1 > m2 ); + unsigned i2 = ( m2 >= m0 ) + ( m2 >= m1 ); + unsigned I = i0 + 3 * ( i1 + 3 * i2 ), J = indx[( I - 5 ) >> 1], + j0 = J & 3, j1 = ( J >> 2 ) & 3, j2 = ( J >> 4 ) & 3; + double wtmax = wtscale[ndaug - 1] * wtscale[ndaug - 1]; + auto order = []( double& a, double& b ) -> void { + auto t = std::max( a, b ); + a = std::min( a, b ); + b = t; + }; + order( m0, m1 ); + order( m1, m2 ); + double u0 = m0 + m1, v0 = m0 - m1, M1max = mp - m2, M02 = sqr( M1max ); + double u02 = u0 * u0, v02 = v0 * v0; + order( m0, m1 ); + double u1 = u0 + m2, v1 = u0 - m2, M12 = sqr( mp ), dE = mp - u1; + if ( dE <= 0 ) { + printf( "Not enough energy: Mtot = %f Etot = %f\n", u1, mp ); + exit( -1 ); } + wtmax /= M12 * M02; + wtmax *= ( M02 - u02 ) * ( M02 - v02 ) * ( M12 - u1 * u1 ) * + ( M12 - v1 * v1 ); + double wt, wtd, R, M1, p0, p1; do { - rnd[0] = 1.0; - il1u = ndaug - 1; - - for ( il1 = 2; il1 < il1u + 1; il1++ ) { - ran = EvtRandom::Flat(); - for ( il2r = 2; il2r < il1 + 1; il2r++ ) { - il2 = il1 + 1 - il2r; - if ( ran <= rnd[il2 - 1] ) - goto two39; - rnd[il2 + 1 - 1] = rnd[il2 - 1]; - } - two39: - rnd[il2 + 1 - 1] = ran; - } + M1 = u0 + dE * EvtRandom::random(); + R = wtmax * sqr( EvtRandom::random() ); + double M02 = M1 * M1; + p0 = ( M02 - u02 ) * ( M02 - v02 ); + double u1 = M1 + m2, v1 = M1 - m2; + p1 = ( M12 - u1 * u1 ) * ( M12 - v1 * v1 ); + wt = p0 * p1; + wtd = M02 * M12; + } while ( wt < wtd * R ); + double iM1 = 1 / M1; + p0 = sqrt( p0 ) * 0.5 * iM1; + p1 = sqrt( p1 ) / ( 2 * mp ); + + // for tree particle decay all momenta lie in a plane so let us + // generate momenta in the x-y plane with the third particle + // momentum along x-axis and then randomly rotate them + double E3 = sqrt( p1 * p1 + m2 * m2 ); + double E12 = mp - E3; + + // energies of particles 1 and 2 in their rest frame where they are back-to-back + double E1 = sqrt( p0 * p0 + m0 * m0 ), E2 = sqrt( p0 * p0 + m1 * m1 ); + + // projection on the x-y plane of uniform rotations in 3D + double cth = EvtRandom::random() * 2 - 1, sth = sqrt( 1 - cth * cth ); + double px = p0 * cth, py = p0 * sth; + + // boost 1 and 2 particle into the 1-2-3 rest frame + double g = E12 * iM1, gb = p1 * iM1; + double gpx = g * px, gbpx = gb * px; + double v0x = gpx - gb * E1, v0e = g * E1 - gbpx; + double v1x = -gpx - gb * E2, v1e = g * E2 + gbpx; + double v2x = p1; + + double x, y; + usincos( EvtRandom::urandom(), x, y ); + double c, s; + usincos( EvtRandom::urandom(), c, s ); + double u = EvtRandom::random() * 2, z = u - 1, r = sqrt( 1 - z * z ); + double Rx = s * y - c * x, Ry = c * y + s * x, ux = u * x, uy = u * y; + double R00 = ux * Rx + c, R01 = s - ux * Ry, R10 = uy * Rx - s, + R11 = c - uy * Ry, R20 = -r * Rx, R21 = r * Ry; + double pyx = R01 * py, pyy = R11 * py, pyz = R21 * py; + + p4[j0].set( v0e, R00 * v0x + pyx, R10 * v0x + pyy, R20 * v0x + pyz ); + p4[j1].set( v1e, R00 * v1x - pyx, R10 * v1x - pyy, R20 * v1x - pyz ); + p4[j2].set( E3, R00 * v2x, R10 * v2x, R20 * v2x ); + } else if ( ndaug < 16 ) { + const int nmax = 15; + double M[nmax], rndf[nmax], p[nmax - 1]; + const double* m = mass; + + double E0 = 0.0; + for ( int i = 0; i < ndaug; i++ ) { + E0 += m[i]; + M[i] = m[i]; + } + double dE = mp - E0; + if ( dE <= 0 ) { + printf( "Not enough energy: Mtot = %f Etot = %f\n", E0, mp ); + exit( -1 ); + } + rndf[0] = 0.0; + rndf[ndaug - 1] = 1; + + double wtmax = wtscale[ndaug - 1] * wtscale[ndaug - 1]; + // std::sort(M, M + ndaug, std::less()); + sortfuns[ndaug - 1]( M ); + m = M; + double Mmin = 0.0, Mmax = dE + m[0], wtmaxd = 1.0; + for ( int i = 0; i < ndaug - 1; i++ ) { + Mmin += m[i]; + Mmax += m[i + 1]; + double u = Mmin + m[i + 1], v = Mmin - m[i + 1], M2 = Mmax * Mmax; + wtmax *= ( M2 - u * u ) * ( M2 - v * v ); + wtmaxd *= M2; + } + wtmax /= wtmaxd; + m = mass; + fun_t* sortfun = sortfuns[ndaug - 3]; + double wt, wtd, R; + do { + for ( int i = 1; i < ndaug - 1; i++ ) + rndf[i] = EvtRandom::random(); + sortfun( rndf + 1 ); - rnd[ndaug - 1] = 0.0; wt = 1.0; - for ( ilr = 2; ilr < ndaug + 1; ilr++ ) { - il = ndaug + 1 - ilr; - pm[4][il - 1] = pm[4][il + 1 - 1] + mass[il - 1] + - ( rnd[il - 1] - rnd[il + 1 - 1] ) * ( mp - psum ); - wt = wt * - EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] ); - } - if ( wt > wtmax ) { - EvtGenReport( EVTGEN_ERROR, "EvtGen" ) - << "wtmax to small in EvtPhaseSpace with " << ndaug - << " daughters" << endl; - ; - } - } while ( wt < EvtRandom::Flat( wtmax ) ); - - ilu = ndaug - 1; - - for ( il = 1; il < ilu + 1; il++ ) { - pa = EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] ); - costh = EvtRandom::Flat( -1.0, 1.0 ); - sinth = sqrt( 1.0 - costh * costh ); - phi = EvtRandom::Flat( EvtConst::twoPi ); - p[1][il - 1] = pa * sinth * cos( phi ); - p[2][il - 1] = pa * sinth * sin( phi ); - p[3][il - 1] = pa * costh; - pm[1][il + 1 - 1] = -p[1][il - 1]; - pm[2][il + 1 - 1] = -p[2][il - 1]; - pm[3][il + 1 - 1] = -p[3][il - 1]; - p[0][il - 1] = sqrt( pa * pa + mass[il - 1] * mass[il - 1] ); - pm[0][il + 1 - 1] = sqrt( pa * pa + - pm[4][il + 1 - 1] * pm[4][il + 1 - 1] ); + wtd = 1.0; + R = wtmax * sqr( EvtRandom::random() ); + double M0 = m[0]; + int i = 1; + do { + double f = rndf[i] - rndf[i - 1]; + double m1 = m[i], M1 = m1 + M0 + f * dE; + double ma = M0 + m1, ms = M0 - m1, M2 = M1 * M1; + double t = ( M2 - ma * ma ) * ( M2 - ms * ms ); + wt *= t; + wtd *= M2; + p[i - 1] = t; + M[i] = M1; + M0 = M1; + } while ( ++i < ndaug ); + } while ( wt < wtd * R ); + + if ( wt > wtd * wtmax ) { + printf( "Warning: current weight is higher than supposed maximum: %e > %e\n", + sqrt( wt / wtd ), sqrt( wtmax ) ); } - p[0][ndaug - 1] = pm[0][ndaug - 1]; - p[1][ndaug - 1] = pm[1][ndaug - 1]; - p[2][ndaug - 1] = pm[2][ndaug - 1]; - p[3][ndaug - 1] = pm[3][ndaug - 1]; - - for ( ilr = 2; ilr < ndaug + 1; ilr++ ) { - il = ndaug + 1 - ilr; - be[0] = pm[0][il - 1] / pm[4][il - 1]; - be[1] = pm[1][il - 1] / pm[4][il - 1]; - be[2] = pm[2][il - 1] / pm[4][il - 1]; - be[3] = pm[3][il - 1] / pm[4][il - 1]; - - for ( i1 = il; i1 < ndaug + 1; i1++ ) { - bep = be[1] * p[1][i1 - 1] + be[2] * p[2][i1 - 1] + - be[3] * p[3][i1 - 1] + be[0] * p[0][i1 - 1]; - temp = ( p[0][i1 - 1] + bep ) / ( be[0] + 1.0 ); - p[1][i1 - 1] = p[1][i1 - 1] + temp * be[1]; - p[2][i1 - 1] = p[2][i1 - 1] + temp * be[2]; - p[3][i1 - 1] = p[3][i1 - 1] + temp * be[3]; - p[0][i1 - 1] = bep; - } + double* iM = rndf; + for ( int i = 0; i < ndaug - 1; i++ ) { + iM[i + 1] = 1 / M[i + 1]; + p[i] = sqrt( p[i] ) * ( 0.5 * iM[i + 1] ); } - for ( ilr = 0; ilr < ndaug; ilr++ ) { - p4[ilr].set( p[0][ilr], p[1][ilr], p[2][ilr], p[3][ilr] ); - } + p4[0].set( sqrt( p[0] * p[0] + m[0] * m[0] ), 0, p[0], 0 ); + int i = 1; + while ( 1 ) { + p4[i].set( sqrt( p[i - 1] * p[i - 1] + m[i] * m[i] ), 0, -p[i - 1], + 0 ); + + double cz = EvtRandom::random() * 2 - 1, sz = sqrt( 1 - cz * cz ); + double sy, cy; + usincos( EvtRandom::urandom(), sy, cy ); + for ( int j = 0; j <= i; j++ ) { + double x = p4[j].get( 1 ), y = p4[j].get( 2 ), z = p4[j].get( 3 ); + double xp = cz * x - sz * y, + yp = sz * x + cz * y; // rotation around z + double zp = sy * xp + cy * z; + xp = cy * xp - sy * z; // rotation around y + p4[j].set( 1, xp ); + p4[j].set( 2, yp ); + p4[j].set( 3, zp ); + } + + if ( i == ( ndaug - 1 ) ) + break; - return 1.0; + double E = sqrt( p[i] * p[i] + M[i] * M[i] ), gamma = E * iM[i], + betagamma = p[i] * iM[i]; + for ( int j = 0; j <= i; j++ ) { + double e = p4[j].get( 0 ), py = p4[j].get( 2 ); + p4[j].set( 0, gamma * e + betagamma * py ); + p4[j].set( 2, gamma * py + betagamma * e ); + } + i++; + } + } else { + printf( "No more than 15 particles! Ndaughter = %d", ndaug ); + exit( -1 ); } return 1.0; } -double EvtGenKine::PhaseSpacePole( double M, double m1, double m2, double m3, - double a, EvtVector4R p4[10] ) +double PhaseSpacePole1( double M, double m1, double m2, double m3, double a, + EvtVector4R p4[10] ) // generate kinematics for 3 body decays, pole for the m1,m2 mass. - { - //f1 = 1 (phasespace) - //f2 = a*(1/m12sq)^2 - - double m12sqmax = ( M - m3 ) * ( M - m3 ); - double m12sqmin = ( m1 + m2 ) * ( m1 + m2 ); - - double m13sqmax = ( M - m2 ) * ( M - m2 ); - double m13sqmin = ( m1 + m3 ) * ( m1 + m3 ); - - double v1 = ( m12sqmax - m12sqmin ) * ( m13sqmax - m13sqmin ); - double v2 = a * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) * ( m13sqmax - m13sqmin ); - + //f1 = 1 (phasespace) + //f2 = a*(1/(p4[0]+p4[1])^2) + + double m12sqmin = ( m1 + m2 ) * ( m1 + m2 ), + m12sqmax = ( M - m3 ) * ( M - m3 ); + double m13sqmin = ( m1 + m3 ) * ( m1 + m3 ), + m13sqmax = ( M - m2 ) * ( M - m2 ); + double d12 = m12sqmax - m12sqmin, d13 = m13sqmax - m13sqmin; + double M2 = M * M, m32 = m3 * m3, m12 = m1 * m1, m22 = m2 * m2; + double c0 = M2 - m32, c1 = m12 - m22, c2 = m32 + m12; + double ab12 = log( m12sqmax / m12sqmin ); // \int_m12sqmin^m12sqmax dx/x + double r = d12 / ( d12 + a * ab12 ); double m12sq, m13sq; + do { + double z0 = EvtRandom::random(), z1 = EvtRandom::random(), + z2 = EvtRandom::random(); + m13sq = m13sqmin + z0 * d13; + m12sq = ( r > z1 ) ? m12sqmin + z2 * d12 : m12sqmin * exp( z2 * ab12 ); + double E3s = c0 - m12sq, E1s = m12sq + c1; + double w = 2 * m12sq, e = 4 * m12sq; + double A = ( m13sq - c2 ) * w - E3s * E1s; + double B = ( E3s * E3s - e * m32 ) * ( E1s * E1s - e * m12 ); + if ( A * A < B ) + break; + } while ( true ); + double iM = 0.5 / M; + double E2 = ( M2 + m22 - m13sq ) * iM; + double E3 = ( M2 + m32 - m12sq ) * iM; + double E1 = M - E2 - E3; + double p1 = sqrt( E1 * E1 - m12 ); + double p3 = sqrt( E3 * E3 - m32 ); + double cost13 = ( 2.0 * E1 * E3 + ( m12 + m32 - m13sq ) ) / ( 2.0 * p1 * p3 ); + + double px = p1 * cost13; + double v0x = px; + double v1x = -px - p3; + double py = p1 * sqrt( 1.0 - cost13 * cost13 ); + double v2x = p3; + + double x, y; + usincos( EvtRandom::urandom(), x, y ); + double c, s; + usincos( EvtRandom::urandom(), c, s ); + double u = EvtRandom::random() * 2, z = u - 1, t = sqrt( 1 - z * z ); + double Rx = s * y - c * x, Ry = c * y + s * x, ux = u * x, uy = u * y; + double R00 = ux * Rx + c, R01 = s - ux * Ry, R10 = uy * Rx - s, + R11 = c - uy * Ry, R20 = -t * Rx, R21 = t * Ry; + double pyx = R01 * py, pyy = R11 * py, pyz = R21 * py; + + p4[0].set( E1, R00 * v0x + pyx, R10 * v0x + pyy, R20 * v0x + pyz ); + p4[1].set( E2, R00 * v1x - pyx, R10 * v1x - pyy, R20 * v1x - pyz ); + p4[2].set( E3, R00 * v2x, R10 * v2x, R20 * v2x ); + + return 1.0 + a / m12sq; +} - double r = v1 / ( v1 + v2 ); - - double m13min, m13max; - +double EvtGenKine::PhaseSpacePole2( double M, double m1, double m2, double m3, + EvtVector4R* p4, const EvtLinSample& g ) +{ + // static linsample g("b2sllprob_wide.dat"); + // static linsample g("b2sllprob_amp.dat"); + + double m13sqmin = ( m1 + m3 ) * ( m1 + m3 ), + m13sqmax = ( M - m2 ) * ( M - m2 ); + double d13 = m13sqmax - m13sqmin; + double M2 = M * M, m32 = m3 * m3, m12 = m1 * m1, m22 = m2 * m2; + double c0 = M2 - m32, c1 = m12 - m22, c2 = m32 + m12; + double m12sq, m13sq, w; do { - m13sq = EvtRandom::Flat( m13sqmin, m13sqmax ); - - if ( r > EvtRandom::Flat() ) { - m12sq = EvtRandom::Flat( m12sqmin, m12sqmax ); - } else { - m12sq = 1.0 / - ( 1.0 / m12sqmin - - EvtRandom::Flat() * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) ); - } + double z0 = EvtRandom::random(), z1 = EvtRandom::random(); + m13sq = m13sqmin + z0 * d13; + auto t = g( z1 ); + m12sq = t.first; + w = t.second; + double E3s = c0 - m12sq, E1s = m12sq + c1; + double w = 2 * m12sq, e = 4 * m12sq; + double A = ( m13sq - c2 ) * w - E3s * E1s; + double B = ( E3s * E3s - e * m32 ) * ( E1s * E1s - e * m12 ); + if ( A * A < B ) + break; + } while ( true ); + double iM = 0.5 / M; + double E2 = ( M2 + m22 - m13sq ) * iM; + double E3 = ( M2 + m32 - m12sq ) * iM; + double E1 = M - E2 - E3; + double p1 = sqrt( E1 * E1 - m12 ); + double p3 = sqrt( E3 * E3 - m32 ); + double cost13 = ( 2.0 * E1 * E3 + ( m12 + m32 - m13sq ) ) / ( 2.0 * p1 * p3 ); + + double px = p1 * cost13; + double v0x = px; + double v1x = -px - p3; + double py = p1 * sqrt( 1.0 - cost13 * cost13 ); + double v2x = p3; + + double x, y; + usincos( EvtRandom::urandom(), x, y ); + double c, s; + usincos( EvtRandom::urandom(), c, s ); + double u = EvtRandom::random() * 2, z = u - 1, t = sqrt( 1 - z * z ); + double Rx = s * y - c * x, Ry = c * y + s * x, ux = u * x, uy = u * y; + double R00 = ux * Rx + c, R01 = s - ux * Ry, R10 = uy * Rx - s, + R11 = c - uy * Ry, R20 = -t * Rx, R21 = t * Ry; + double pyx = R01 * py, pyy = R11 * py, pyz = R21 * py; + + p4[0].set( E1, R00 * v0x + pyx, R10 * v0x + pyy, R20 * v0x + pyz ); + p4[1].set( E2, R00 * v1x - pyx, R10 * v1x - pyy, R20 * v1x - pyz ); + p4[2].set( E3, R00 * v2x, R10 * v2x, R20 * v2x ); + + return w; +} - //kinematically allowed? - double E3star = ( M * M - m12sq - m3 * m3 ) / sqrt( 4 * m12sq ); - double E1star = ( m12sq + m1 * m1 - m2 * m2 ) / sqrt( 4 * m12sq ); - double p3star = sqrt( E3star * E3star - m3 * m3 ); - double p1star = sqrt( E1star * E1star - m1 * m1 ); - m13max = ( E3star + E1star ) * ( E3star + E1star ) - - ( p3star - p1star ) * ( p3star - p1star ); - m13min = ( E3star + E1star ) * ( E3star + E1star ) - - ( p3star + p1star ) * ( p3star + p1star ); +double EvtGenKine::PhaseSpacePole( double M, double m1, double m2, double m3, + double a, EvtVector4R p4[10] ) - } while ( m13sq < m13min || m13sq > m13max ); +// generate kinematics for 3 body decays, pole for the m1,m2 mass. - double E2 = ( M * M + m2 * m2 - m13sq ) / ( 2.0 * M ); - double E3 = ( M * M + m3 * m3 - m12sq ) / ( 2.0 * M ); +{ + // return PhaseSpacePole2(M, m1, m2, m3, p4); + return PhaseSpacePole1( M, m1, m2, m3, a, p4 ); + //f1 = 1 (phasespace) + //f2 = a*(1/m12sq)^2 = a*(1/(p4[0]+p4[1])^4) + + double m12sqmin = ( m1 + m2 ) * ( m1 + m2 ), + m12sqmax = ( M - m3 ) * ( M - m3 ); + double m13sqmin = ( m1 + m3 ) * ( m1 + m3 ), + m13sqmax = ( M - m2 ) * ( M - m2 ); + double d12 = m12sqmax - m12sqmin, d13 = m13sqmax - m13sqmin; + double M2 = M * M, m32 = m3 * m3, m12 = m1 * m1, m22 = m2 * m2; + double c0 = M2 - m32, c1 = m12 - m22, c2 = m32 + m12; + double a12 = 1.0 / m12sqmin, ab12 = a12 - 1.0 / m12sqmax; + double r = d12 / ( d12 + a * ab12 ); + double m12sq, m13sq; + do { + double z0 = EvtRandom::random(), z1 = EvtRandom::random(), + z2 = EvtRandom::random(); + m13sq = m13sqmin + z0 * d13; + m12sq = ( r > z1 ) ? m12sqmin + z2 * d12 : 1.0 / ( a12 - z2 * ab12 ); + double E3s = c0 - m12sq, E1s = m12sq + c1; + double w = 2 * m12sq, e = 4 * m12sq; + double A = ( m13sq - c2 ) * w - E3s * E1s; + double B = ( E3s * E3s - e * m32 ) * ( E1s * E1s - e * m12 ); + if ( A * A < B ) + break; + } while ( true ); + double iM = 0.5 / M; + double E2 = ( M2 + m22 - m13sq ) * iM; + double E3 = ( M2 + m32 - m12sq ) * iM; double E1 = M - E2 - E3; - double p1mom = sqrt( E1 * E1 - m1 * m1 ); - double p3mom = sqrt( E3 * E3 - m3 * m3 ); - double cost13 = ( 2.0 * E1 * E3 + m1 * m1 + m3 * m3 - m13sq ) / - ( 2.0 * p1mom * p3mom ); - - //EvtGenReport(EVTGEN_INFO,"EvtGen") << m13sq << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << m12sq << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << E1 << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << E2 << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << E3 << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << p1mom << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << p3mom << endl; - //EvtGenReport(EVTGEN_INFO,"EvtGen") << cost13 << endl; - - p4[2].set( E3, 0.0, 0.0, p3mom ); - p4[0].set( E1, p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, p1mom * cost13 ); - p4[1].set( E2, -p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, - -p1mom * cost13 - p3mom ); - - //EvtGenReport(EVTGEN_INFO,"EvtGen") << "p4:"< +#include +#include +#include + +EvtLinSample::EvtLinSample( const char* fname ) +{ + std::ifstream IN( fname ); + float x, y; + while ( IN >> x >> y ) { + v.push_back( { x, y } ); + } + std::cout << "Point read -- " << v.size() << std::endl; + I.push_back( 0 ); + for ( unsigned i = 0; i < v.size() - 1; i++ ) { + I.push_back( ( v[i + 1].y + v[i + 0].y ) * ( v[i + 1].x - v[i + 0].x ) + + I.back() ); + } + double norm = 1 / I.back(); + for ( unsigned i = 0; i < v.size(); i++ ) + I[i] *= norm; +} + +void EvtLinSample::init( const std::vector& _v ) +{ + v = _v; + std::cout << "Point provided -- " << v.size() << std::endl; + I.push_back( 0 ); + for ( unsigned i = 0; i < v.size() - 1; i++ ) { + I.push_back( ( v[i + 1].y + v[i + 0].y ) * ( v[i + 1].x - v[i + 0].x ) + + I.back() ); + } + double norm = 1 / I.back(); + for ( unsigned i = 0; i < v.size(); i++ ) + I[i] *= norm; +} + +std::pair EvtLinSample::operator()( double r ) const +{ + int j = upper_bound( I.begin(), I.end(), r ) - I.begin(); + double dI = I[j] - I[j - 1]; + r = ( r - I[j - 1] ) / dI; + double f0 = v[j - 1].y, f1 = v[j].y, x0 = v[j - 1].x, x1 = v[j].x, + df = f1 - f0, dx = x1 - x0, z; + if ( fabs( df ) > f0 * 1e-3 ) { + z = ( f1 * x0 - x1 * f0 + dx * sqrt( df * ( f0 + f1 ) * r + f0 * f0 ) ) / + df; + } else { + if ( f0 > f1 ) { + z = x0 + ( r * dx ) * ( f1 - r * df ) / f0; + } else { + z = x1 - ( r * dx ) * ( f0 + r * df ) / f1; + } + } + return std::make_pair( z, f0 + ( df / dx ) * ( z - x0 ) ); +} diff --git a/src/EvtGenBase/EvtMTRandomEngine.cpp b/src/EvtGenBase/EvtMTRandomEngine.cpp --- a/src/EvtGenBase/EvtMTRandomEngine.cpp +++ b/src/EvtGenBase/EvtMTRandomEngine.cpp @@ -25,8 +25,7 @@ #include -EvtMTRandomEngine::EvtMTRandomEngine( unsigned int seed ) : - engine_( seed ), distribution_( URDist( 0.0, 1.0 ) ) +EvtMTRandomEngine::EvtMTRandomEngine( unsigned int seed ) : engine_( seed ) { EvtGenReport( EVTGEN_INFO, "EvtMTRandomEngine" ) << "Mersenne-Twister random number generator with seed = " << seed @@ -35,5 +34,11 @@ double EvtMTRandomEngine::random() { - return distribution_( engine_ ); + return (long)( engine_() >> 1 ) * ( 1.0 / ( 1ul << 63 ) ); + // return distribution_( engine_ ); +} + +unsigned long EvtMTRandomEngine::urandom() +{ + return engine_(); } diff --git a/src/EvtGenBase/EvtModel.cpp b/src/EvtGenBase/EvtModel.cpp --- a/src/EvtGenBase/EvtModel.cpp +++ b/src/EvtGenBase/EvtModel.cpp @@ -96,7 +96,7 @@ model = _commandNameHash[cmd]; } - assert( model != 0 ); + assert( model != nullptr ); model->command( cnfgstr ); } diff --git a/src/EvtGenBase/EvtParticle.cpp b/src/EvtGenBase/EvtParticle.cpp --- a/src/EvtGenBase/EvtParticle.cpp +++ b/src/EvtGenBase/EvtParticle.cpp @@ -101,7 +101,9 @@ void EvtParticle::setLifetime() { if ( _genlifetime ) { - _t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( getId() ); + _t = EvtPDL::getctau( getId() ); + if ( _t > 0 ) + _t *= -log( EvtRandom::Flat() ); } } diff --git a/src/EvtGenBase/EvtRandom.cpp b/src/EvtGenBase/EvtRandom.cpp --- a/src/EvtGenBase/EvtRandom.cpp +++ b/src/EvtGenBase/EvtRandom.cpp @@ -51,6 +51,11 @@ return _randomEngine->random(); } +unsigned long EvtRandom::urandom() +{ + return _randomEngine->urandom(); +} + // Random number routine to generate numbers between // min and max. By djl on July 27, 1995. double EvtRandom::Flat( double min, double max ) diff --git a/src/EvtGenBase/EvtSemiLeptonicVectorAmp.cpp b/src/EvtGenBase/EvtSemiLeptonicVectorAmp.cpp --- a/src/EvtGenBase/EvtSemiLeptonicVectorAmp.cpp +++ b/src/EvtGenBase/EvtSemiLeptonicVectorAmp.cpp @@ -33,6 +33,15 @@ #include "EvtGenBase/EvtVector4C.hh" using std::endl; +double gmn( int m, int n ) +{ + if ( m != n ) + return 0; + if ( m > 0 ) + return -1; + return 1; +} + void EvtSemiLeptonicVectorAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp, EvtSemiLeptonicFF* FormFactors ) { @@ -118,9 +127,40 @@ } } - EvtVector4C et0 = tds.cont1( parent->getDaug( 0 )->epsParent( 0 ).conj() ); - EvtVector4C et1 = tds.cont1( parent->getDaug( 0 )->epsParent( 1 ).conj() ); - EvtVector4C et2 = tds.cont1( parent->getDaug( 0 )->epsParent( 2 ).conj() ); + EvtParticle* D = parent->getDaug( 0 ); + EvtVector4C t0 = D->epsParent( 0 ), t1 = D->epsParent( 1 ), + t2 = D->epsParent( 2 ); + /* + EvtVector4R k = D->getP4(); + std::cout< +#include #include #include #include @@ -35,46 +35,25 @@ { dim = 0; rho = nullptr; - - int i, j; setDim( density.dim ); - - for ( i = 0; i < dim; i++ ) { - for ( j = 0; j < dim; j++ ) { - rho[i][j] = density.rho[i][j]; - } - } + // memmove(rho,density.rho,dim*dim*sizeof(EvtComplex)); + for ( int i = 0; i < dim * dim; i++ ) + rho[i] = density.rho[i]; } EvtSpinDensity& EvtSpinDensity::operator=( const EvtSpinDensity& density ) { - int i, j; setDim( density.dim ); - - for ( i = 0; i < dim; i++ ) { - for ( j = 0; j < dim; j++ ) { - rho[i][j] = density.rho[i][j]; - } - } - + // memmove(rho,density.rho,dim*dim*sizeof(EvtComplex)); + for ( int i = 0; i < dim * dim; i++ ) + rho[i] = density.rho[i]; return *this; } EvtSpinDensity::~EvtSpinDensity() { - if ( dim != 0 ) { - int i; - for ( i = 0; i < dim; i++ ) - delete[] rho[i]; - } - - delete[] rho; -} - -EvtSpinDensity::EvtSpinDensity() -{ - dim = 0; - rho = nullptr; + if ( dim != 0 ) + delete[] rho; } void EvtSpinDensity::setDim( int n ) @@ -82,9 +61,6 @@ if ( dim == n ) return; if ( dim != 0 ) { - int i; - for ( i = 0; i < dim; i++ ) - delete[] rho[i]; delete[] rho; rho = nullptr; dim = 0; @@ -92,61 +68,33 @@ if ( n == 0 ) return; dim = n; - rho = new EvtComplexPtr[n]; - int i; - for ( i = 0; i < n; i++ ) { - rho[i] = new EvtComplex[n]; - } -} - -int EvtSpinDensity::getDim() const -{ - return dim; -} - -void EvtSpinDensity::set( int i, int j, const EvtComplex& rhoij ) -{ - assert( i < dim && j < dim ); - rho[i][j] = rhoij; -} - -const EvtComplex& EvtSpinDensity::get( int i, int j ) const -{ - assert( i < dim && j < dim ); - return rho[i][j]; + rho = new EvtComplex[n * n]; } void EvtSpinDensity::setDiag( int n ) { setDim( n ); - int i, j; - - for ( i = 0; i < n; i++ ) { - for ( j = 0; j < n; j++ ) { - rho[i][j] = EvtComplex( 0.0 ); - } - rho[i][i] = EvtComplex( 1.0 ); - } + for ( int i = 0; i < n * n; i++ ) + rho[i] = 0; + for ( int i = 0; i < n; i++ ) + rho[i * dim + i] = 1.0; } -double EvtSpinDensity::normalizedProb( const EvtSpinDensity& d ) +double EvtSpinDensity::normalizedProb( const EvtSpinDensity& d ) const { - int i, j; - EvtComplex prob( 0.0, 0.0 ); - double norm = 0.0; - if ( dim != d.dim ) { EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Not matching dimensions in NormalizedProb" << endl; ::abort(); } - for ( i = 0; i < dim; i++ ) { - norm += real( rho[i][i] ); - for ( j = 0; j < dim; j++ ) { - prob += rho[i][j] * d.rho[i][j]; - } - } + double norm = 0.0; + for ( int i = 0; i < dim; i++ ) + norm += real( rho[i * dim + i] ); + + EvtComplex prob; + for ( int i = 0, imax = dim * dim; i < imax; i++ ) + prob += rho[i] * d.rho[i]; if ( imag( prob ) > 0.00000001 * real( prob ) ) { EvtGenReport( EVTGEN_ERROR, "EvtGen" ) @@ -172,13 +120,13 @@ double trace( 0.0 ); for ( i = 0; i < dim; i++ ) { - trace += abs( rho[i][i] ); + trace += abs( rho[i * dim + i] ); } for ( i = 0; i < dim; i++ ) { - if ( real( rho[i][i] ) < 0.0 ) + if ( real( rho[i * dim + i] ) < 0.0 ) return 0; - if ( imag( rho[i][i] ) * 1000000.0 > trace ) { + if ( imag( rho[i * dim + i] ) * 1000000.0 > trace ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << *this << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << trace << endl; EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 1" << endl; @@ -188,13 +136,15 @@ for ( i = 0; i < dim; i++ ) { for ( j = i + 1; j < dim; j++ ) { - if ( fabs( real( rho[i][j] - rho[j][i] ) ) > - 0.00000001 * ( abs( rho[i][i] ) + abs( rho[j][j] ) ) ) { + if ( fabs( real( rho[i * dim + j] - rho[j * dim + i] ) ) > + 0.00000001 * + ( abs( rho[i * dim + i] ) + abs( rho[j * dim + j] ) ) ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 2" << endl; return 0; } - if ( fabs( imag( rho[i][j] + rho[j][i] ) ) > - 0.00000001 * ( abs( rho[i][i] ) + abs( rho[j][j] ) ) ) { + if ( fabs( imag( rho[i * dim + j] + rho[j * dim + i] ) ) > + 0.00000001 * + ( abs( rho[i * dim + i] ) + abs( rho[j * dim + j] ) ) ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Failing 3" << endl; return 0; } @@ -206,14 +156,12 @@ ostream& operator<<( ostream& s, const EvtSpinDensity& d ) { - int i, j; - s << endl; s << "Dimension:" << d.dim << endl; - for ( i = 0; i < d.dim; i++ ) { - for ( j = 0; j < d.dim; j++ ) { - s << d.rho[i][j] << " "; + for ( int i = 0; i < d.dim; i++ ) { + for ( int j = 0; j < d.dim; j++ ) { + s << d.get( i, j ) << " "; } s << endl; } diff --git a/src/EvtGenBase/EvtTensor4C.cpp b/src/EvtGenBase/EvtTensor4C.cpp --- a/src/EvtGenBase/EvtTensor4C.cpp +++ b/src/EvtGenBase/EvtTensor4C.cpp @@ -272,7 +272,7 @@ for ( i = 0; i < 4; i++ ) { for ( j = 0; j < 4; j++ ) { - t[i][j] *= EvtComplex( d, 0.0 ); + t[i][j] *= d; } } return *this; @@ -280,12 +280,12 @@ EvtTensor4C operator*( const EvtTensor4C& t1, double d ) { - return EvtTensor4C( t1 ) *= EvtComplex( d, 0.0 ); + return EvtTensor4C( t1 ) *= d; } EvtTensor4C operator*( double d, const EvtTensor4C& t1 ) { - return EvtTensor4C( t1 ) *= EvtComplex( d, 0.0 ); + return EvtTensor4C( t1 ) *= d; } EvtComplex cont( const EvtTensor4C& t1, const EvtTensor4C& t2 ) @@ -553,3 +553,120 @@ } } } + +EvtTensor4C& EvtTensor4C::addScaled( EvtComplex s, const EvtTensor4C& a ) +{ + EvtComplex* u0 = reinterpret_cast( t ); + const EvtComplex* u1 = reinterpret_cast( a.t ); + for ( int i = 0; i < 16; i++ ) + u0[i] += s * u1[i]; + return *this; +} + +EvtTensor4C EvtGenFunctions::asymProd( const EvtVector4R& a, const EvtVector4R& b ) +{ + // t_{ij} = eps_{ijkl} a^{k} b^{l} + // eps_{ijkl} is the Levi-Civita symbol -- antisymmetric tensor + + EvtTensor4C res; + EvtComplex( &t )[4][4] = res.t; + double t01 = a.get( 2 ) * b.get( 3 ) - a.get( 3 ) * b.get( 2 ); + t[0][1] = t01; + t[1][0] = -t01; + double t02 = a.get( 3 ) * b.get( 1 ) - a.get( 1 ) * b.get( 3 ); + t[0][2] = t02; + t[2][0] = -t02; + double t03 = a.get( 1 ) * b.get( 2 ) - a.get( 2 ) * b.get( 1 ); + t[0][3] = t03; + t[3][0] = -t03; + double t12 = a.get( 0 ) * b.get( 3 ) - a.get( 3 ) * b.get( 0 ); + t[1][2] = t12; + t[2][1] = -t12; + double t13 = a.get( 2 ) * b.get( 0 ) - a.get( 0 ) * b.get( 2 ); + t[1][3] = t13; + t[3][1] = -t13; + double t23 = a.get( 0 ) * b.get( 1 ) - a.get( 1 ) * b.get( 0 ); + t[2][3] = t23; + t[3][2] = -t23; + return res; +} + +EvtVector4R EvtGenFunctions::asymProd( const EvtVector4R& _a, + const EvtVector4R& _b, + const EvtVector4R& _c ) +{ + // t_{ij} = eps_{ijkl} a^{j} b^{k} c^{l} + // eps_{ijkl} is the Levi-Civita symbol -- antisymmetric tensor + + double a[4] = { _a.get( 0 ), _a.get( 1 ), _a.get( 2 ), _a.get( 3 ) }; + /* + double b[4] = {_b.get(0), _b.get(1), _b.get(2), _b.get(3)}; + double c[4] = {_c.get(0), _c.get(1), _c.get(2), _c.get(3)}; + double t[4] = {0}; + + t[0] += a[1]*b[2]*c[3]; + t[0] -= a[1]*b[3]*c[2]; + t[0] -= a[2]*b[1]*c[3]; + t[0] += a[2]*b[3]*c[1]; + t[0] += a[3]*b[1]*c[2]; + t[0] -= a[3]*b[2]*c[1]; + + t[1] -= a[0]*b[2]*c[3]; + t[1] += a[0]*b[3]*c[2]; + t[1] += a[2]*b[0]*c[3]; + t[1] -= a[2]*b[3]*c[0]; + t[1] -= a[3]*b[0]*c[2]; + t[1] += a[3]*b[2]*c[0]; + + t[2] += a[0]*b[1]*c[3]; + t[2] -= a[0]*b[3]*c[1]; + t[2] -= a[1]*b[0]*c[3]; + t[2] += a[1]*b[3]*c[0]; + t[2] += a[3]*b[0]*c[1]; + t[2] -= a[3]*b[1]*c[0]; + + t[3] -= a[0]*b[1]*c[2]; + t[3] += a[0]*b[2]*c[1]; + t[3] += a[1]*b[0]*c[2]; + t[3] -= a[1]*b[2]*c[0]; + t[3] -= a[2]*b[0]*c[1]; + t[3] += a[2]*b[1]*c[0]; + + EvtVector4R res(t[0],t[1],t[2],t[3]); + // std::cout< using std::ostream; -EvtVector4C::EvtVector4C() -{ - v[0] = EvtComplex( 0.0 ); - v[1] = EvtComplex( 0.0 ); - v[2] = EvtComplex( 0.0 ); - v[3] = EvtComplex( 0.0 ); -} - -EvtVector4C::EvtVector4C( const EvtComplex& e0, const EvtComplex& e1, - const EvtComplex& e2, const EvtComplex& e3 ) -{ - v[0] = e0; - v[1] = e1; - v[2] = e2; - v[3] = e3; -} - EvtVector4C rotateEuler( const EvtVector4C& rs, double alpha, double beta, double gamma ) { @@ -80,50 +63,24 @@ void EvtVector4C::applyBoostTo( const EvtVector3R& boost ) { - double bx, by, bz, gamma, b2; - - bx = boost.get( 0 ); - by = boost.get( 1 ); - bz = boost.get( 2 ); - - double bxx = bx * bx; - double byy = by * by; - double bzz = bz * bz; - - b2 = bxx + byy + bzz; - - if ( b2 == 0.0 ) { - return; - } + double bx = boost.get( 0 ); + double by = boost.get( 1 ); + double bz = boost.get( 2 ); + double b2 = bx * bx + by * by + bz * bz; assert( b2 < 1.0 ); - gamma = 1.0 / sqrt( 1 - b2 ); - - double gb2 = ( gamma - 1.0 ) / b2; - - double gb2xy = gb2 * bx * by; - double gb2xz = gb2 * bx * bz; - double gb2yz = gb2 * by * bz; - - double gbx = gamma * bx; - double gby = gamma * by; - double gbz = gamma * bz; - - EvtComplex e2 = v[0]; - EvtComplex px2 = v[1]; - EvtComplex py2 = v[2]; - EvtComplex pz2 = v[3]; + double gamma2 = 1 / ( 1 - b2 ); + double gamma = sqrt( gamma2 ); + double F = gamma2 / ( gamma + 1 ); - v[0] = gamma * e2 + gbx * px2 + gby * py2 + gbz * pz2; + EvtComplex bp = bx * v[1] + by * v[2] + bz * v[3]; + EvtComplex c = F * bp + gamma * v[0]; - v[1] = gbx * e2 + gb2 * bxx * px2 + px2 + gb2xy * py2 + gb2xz * pz2; - - v[2] = gby * e2 + gb2 * byy * py2 + py2 + gb2xy * px2 + gb2yz * pz2; - - v[3] = gbz * e2 + gb2 * bzz * pz2 + pz2 + gb2yz * py2 + gb2xz * px2; - - return; + v[1] += c * bx; + v[2] += c * by; + v[3] += c * bz; + v[0] = gamma * ( v[0] + bp ); } void EvtVector4C::applyRotateEuler( double phi, double theta, double ksi ) diff --git a/src/EvtGenBase/EvtVectorParticle.cpp b/src/EvtGenBase/EvtVectorParticle.cpp --- a/src/EvtGenBase/EvtVectorParticle.cpp +++ b/src/EvtGenBase/EvtVectorParticle.cpp @@ -56,21 +56,21 @@ setLifetime(); } -void EvtVectorParticle::init( EvtId part_n, const EvtVector4R& p4, - const EvtVector4C& epsin1, - const EvtVector4C& epsin2, - const EvtVector4C& epsin3 ) -{ - _validP4 = true; - setp( p4 ); - setpart_num( part_n ); - - _eps[0] = epsin1; - _eps[1] = epsin2; - _eps[2] = epsin3; - - setLifetime(); -} +// void EvtVectorParticle::init( EvtId part_n, const EvtVector4R& p4, +// const EvtVector4C& epsin1, +// const EvtVector4C& epsin2, +// const EvtVector4C& epsin3 ) +// { +// _validP4 = true; +// setp( p4 ); +// setpart_num( part_n ); + +// _eps[0] = epsin1; +// _eps[1] = epsin2; +// _eps[2] = epsin3; + +// setLifetime(); +// } EvtSpinDensity EvtVectorParticle::rotateToHelicityBasis() const { @@ -121,3 +121,30 @@ return R; } + +EvtVector4C EvtVectorParticle::epsParent( int i ) const +{ + /* + assuming that polarization vectors have this structure: + _eps[0] = ( 0.0, 1.0, 0.0, 0.0 ); + _eps[1] = ( 0.0, 0.0, 1.0, 0.0 ); + _eps[2] = ( 0.0, 0.0, 0.0, 1.0 ); + boost them into the parent frame + */ + const EvtVector4R& p = getP4(); + int j = i + 1; + double e = p.get( 0 ), px = p.get( 1 ), py = p.get( 2 ), pz = p.get( 3 ); + double m2 = ( e * e - px * px ) - ( py * py + pz * pz ), m = sqrt( m2 ); + double pj = p.get( j ), c = pj / ( m * e + m2 ); + double v0 = c * ( e + m ); + double v1 = c * px; + double v2 = c * py; + double v3 = c * pz; + EvtVector4C t; + t.set( v0, v1, v2, v3 ); + t.set( j, t.get( j ) + 1 ); + return t; + // EvtVector4C o = boostTo( _eps[i], p ); + // std::cout<<"epsParent "<* extraModels ) { EvtModel& modelist = EvtModel::instance(); - if ( extraModels ) { for ( std::list::const_iterator it = extraModels->begin(); it != extraModels->end(); ++it ) { @@ -239,6 +240,8 @@ modelist.registerModel( new EvtTauScalarnu ); modelist.registerModel( new EvtKstarnunu ); modelist.registerModel( new EvtbTosllBall ); + modelist.registerModel( new EvtbTosllNP ); + modelist.registerModel( new EvtbTosllNPR ); modelist.registerModel( new EvtBto2piCPiso ); modelist.registerModel( new EvtBtoKpiCPiso ); modelist.registerModel( new EvtSVSCPiso ); diff --git a/src/EvtGenModels/EvtPi0Dalitz.cpp b/src/EvtGenModels/EvtPi0Dalitz.cpp --- a/src/EvtGenModels/EvtPi0Dalitz.cpp +++ b/src/EvtGenModels/EvtPi0Dalitz.cpp @@ -52,30 +52,24 @@ // contractions, which is: // 1/((m_R^2-q^2)^2+m_R^2 Gamma_R^2) 1/(2q^2) (M^2-q^2)^2 beta_l // (1+cos(theta)^2) where we set cos(theta)=1 - auto daughter1 = getDaug( 0 ); - auto daughter2 = getDaug( 1 ); - double q2Min = EvtPDL::getMass( daughter1 ) + EvtPDL::getMass( daughter2 ); + double q2Min = EvtPDL::getMass( getDaug( 0 ) ) + + EvtPDL::getMass( getDaug( 1 ) ); q2Min *= q2Min; double q2Max = EvtPDL::getMass( getParentId() ); q2Max *= q2Max; - const int steps = 20000; - const double step = ( q2Max - q2Min ) / steps; - double maxProb = 0; - for ( int ii = 0; ii < steps; ++ii ) { - double q2 = q2Min + ii * step; - const double mSqDiff = m_m0Sq - q2; - const double q2Sq = q2 * q2; - double prob = ( q2Max - q2 ) * ( q2Max - q2 ) * ( q2 - q2Min ) / ( q2Sq ); - prob *= ( 1.0 / ( mSqDiff * mSqDiff + m_m0SqG0Sq ) ); - // When generating events, we do not start from phase-space, but - // add some pole to it, weight of which is taken into account - // elsewhere - prob /= 1.0 + m_poleSize / ( q2Sq ); - if ( prob > maxProb ) { + const int steps = 100; + const double step = expm1( log( q2Max / q2Min ) / steps ); + double maxProb = 0, q2 = q2Min; + do { + double prob = ( q2Max - q2 ) * ( q2Max - q2 ) * ( q2 - q2Min ); + double m2d = m_m0Sq - q2; + prob /= q2 * ( m2d * m2d + m_m0SqG0Sq ) * ( q2 + m_poleSize ); + if ( prob > maxProb ) maxProb = prob; - } - } - setProbMax( maxProb * 1.05 ); + } while ( ( q2 += q2 * step ) < q2Max ); + setProbMax( maxProb * 1.01 ); + EvtGenReport( EVTGEN_INFO, "EvtGenProbMax" ) + << " " << sqrt( q2 ) << " " << maxProb << std::endl; } void EvtPi0Dalitz::init() @@ -99,41 +93,33 @@ void EvtPi0Dalitz::decay( EvtParticle* p ) { - EvtParticle *ep, *em, *gamma; - setWeight( p->initializePhaseSpace( getNDaug(), getDaugs(), false, - m_poleSize, 0, 1 ) ); - ep = p->getDaug( 0 ); - em = p->getDaug( 1 ); - gamma = p->getDaug( 2 ); - - // the next four lines generates events with a weight such that - // the efficiency for selecting them is good. The parameter below of - // 0.1 is the size of the peak at low q^2 (in arbitrary units). - // The value of 0.1 is appropriate for muons. - // when you use this remember to remove the cut on q^2! - - //ep em invariant mass^2 - double m2 = ( ep->getP4() + em->getP4() ).mass2(); - EvtVector4R q = ep->getP4() + em->getP4(); - //Just use the prob summed over spins... + using EvtGenFunctions::directProd; + double ew = p->initializePhaseSpace( getNDaug(), getDaugs(), false, + m_poleSize, 0, 1 ); + setWeight( ew ); - EvtTensor4C w, v; + const EvtVector4R& pp = p->getDaug( 0 )->getP4(); + const EvtVector4R& pm = p->getDaug( 1 )->getP4(); + const EvtVector4R& pg = p->getDaug( 2 )->getP4(); + EvtVector4R q = pp + pm; + double q2 = q.mass2(), t = pg * q; - v = 2.0 * ( gamma->getP4() * q ) * - EvtGenFunctions::directProd( q, gamma->getP4() ) - - ( gamma->getP4() * q ) * ( gamma->getP4() * q ) * EvtTensor4C::g() - - m2 * EvtGenFunctions::directProd( gamma->getP4(), gamma->getP4() ); + //Just use the prob summed over spins... + //EvtTensor4C v = (2.0 * t) * directProd( q, pg ) - t2 * EvtTensor4C::g() - q2 * directProd( pg, pg ); + //EvtTensor4C w = 4.0 * ( directProd( pp, pm ) + directProd( pm, pp ) - EvtTensor4C::g() * ( pp*pm - pp.mass2() ) ); + //double prob = real( cont( v, w ) ) / ( q2 * q2 ); + EvtTensor4C v = directProd( q, pg ) + .addScaled( -0.5 * t, EvtTensor4C::g() ) + .addScaled( -0.5 * q2 / t, directProd( pg, pg ) ); + EvtTensor4C w = directProd( pp, pm ).addDirProd( pm, pp ).addScaled( + pp.mass2() - pp * pm, EvtTensor4C::g() ); - w = 4.0 * ( EvtGenFunctions::directProd( ep->getP4(), em->getP4() ) + - EvtGenFunctions::directProd( em->getP4(), ep->getP4() ) - - EvtTensor4C::g() * - ( ep->getP4() * em->getP4() - ep->getP4().mass2() ) ); + double prob = 8 * t * real( cont( v, w ) ) / ( q2 * q2 ); + double m2d = m_m0Sq - q2; + prob /= m2d * m2d + m_m0SqG0Sq; - double prob = ( real( cont( v, w ) ) ) / ( m2 * m2 ); - const double m2Diff = m_m0Sq - m2; - prob *= ( 1.0 / ( m2Diff * m2Diff + m_m0SqG0Sq ) ); + // EvtGenReport(EVTGEN_INFO,"EvtGenDDD") <<" "<getDaug( 0 )->getP4(); double norm = 1.0 / pDaug.d3mag(); - + pDaug *= norm; + // std::cout<<"VSS "<eps( i ) ) ); + vertex( i, pDaug * ( p->eps( i ) ) ); return; } diff --git a/src/EvtGenModels/EvtbTosllAmp.cpp b/src/EvtGenModels/EvtbTosllAmp.cpp --- a/src/EvtGenModels/EvtbTosllAmp.cpp +++ b/src/EvtGenModels/EvtbTosllAmp.cpp @@ -34,6 +34,45 @@ #include "EvtGenBase/EvtVector4C.hh" #include "EvtGenBase/EvtVectorParticle.hh" +#include + +double getKineWeight( EvtParticle* p ) +{ + int nd = p->getNDaug(); + double prod = 1; + for ( int i = 0; i < nd - 1; i++ ) { + EvtVector4R M0; + for ( int j = i + 1; j < nd; j++ ) + M0 += p->getDaug( j )->getP4(); + EvtVector4R M1 = M0 + p->getDaug( i )->getP4(); + double Mi12 = M1.mass2(), Mi = M0.mass(), + mi = p->getDaug( i )->getP4().mass(), Mipm = Mi + mi, + Mimm = Mi - mi; + // std::cout << Mi12 << std::endl; + prod *= 0.5 * + sqrt( ( Mi12 - Mipm * Mipm ) * ( Mi12 - Mimm * Mimm ) / Mi12 ); + } + return prod; +} + +double getKineWeight2( EvtParticle* p ) +{ + int nd = p->getNDaug(); + double prod = 1; + for ( int i = nd - 1; i > 0; i-- ) { + EvtVector4R M0; + for ( int j = i - 1; j >= 0; j-- ) + M0 += p->getDaug( j )->getP4(); + EvtVector4R M1 = M0 + p->getDaug( i )->getP4(); + double Mi12 = M1.mass2(), Mi = M0.mass(), + mi = p->getDaug( i )->getP4().mass(), Mipm = Mi + mi, + Mimm = Mi - mi; + prod *= 0.5 * + sqrt( ( Mi12 - Mipm * Mipm ) * ( Mi12 - Mimm * Mimm ) / Mi12 ); + } + return prod; +} + double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson, EvtId lepton1, EvtId lepton2, EvtbTosllFF* FormFactors, double& poleSize ) @@ -85,8 +124,7 @@ lep1->noLifeTime(); lep2->noLifeTime(); - //Initial particle is unpolarized, well it is a scalar so it is - //trivial + //Initial particle is unpolarized, well it is a scalar so it is trivial EvtSpinDensity rho; rho.setDiag( root_part->getSpinStates() ); @@ -95,120 +133,88 @@ double m = root_part->mass(); EvtVector4R p4meson, p4lepton1, p4lepton2, p4w; - double q2max; - - double q2, elepton, plepton; - int i, j; - double erho, prho, costl; double maxfoundprob = 0.0; double prob = -10.0; - int massiter; - double maxpole = 0; - for ( massiter = 0; massiter < 3; massiter++ ) { - mass[0] = EvtPDL::getMeanMass( meson ); - mass[1] = EvtPDL::getMeanMass( lepton1 ); - mass[2] = EvtPDL::getMeanMass( lepton2 ); - if ( massiter == 1 ) { - mass[0] = EvtPDL::getMinMass( meson ); - } - if ( massiter == 2 ) { - mass[0] = EvtPDL::getMaxMass( meson ); - if ( ( mass[0] + mass[1] + mass[2] ) > m ) - mass[0] = m - mass[1] - mass[2] - 0.00001; - } - - q2max = ( m - mass[0] ) * ( m - mass[0] ); + mass[1] = EvtPDL::getMeanMass( lepton1 ); + mass[2] = EvtPDL::getMeanMass( lepton2 ); + std::vector mH; + mH.push_back( EvtPDL::getMeanMass( meson ) ); + //if the particle is narrow dont bother with changing the mass. + double g = EvtPDL::getWidth( meson ); + if ( g > 0 ) { + mH.push_back( EvtPDL::getMinMass( meson ) ); + mH.push_back( + std::min( EvtPDL::getMaxMass( meson ), m - mass[1] - mass[2] ) ); + mH.push_back( mH.front() - g ); + mH.push_back( mH.front() + g ); + } + double q2min = ( mass[1] + mass[2] ) * ( mass[1] + mass[2] ); + for ( double m0 : mH ) { + double q2max = ( m - m0 ) * ( m - m0 ); //loop over q2 - //cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl; - for ( i = 0; i < 25; i++ ) { - //want to avoid picking up the tail of the photon propagator - q2 = ( ( i + 1.5 ) * q2max ) / 26.0; - + for ( int i = 0, n = 25; i < n; i++ ) { + // want to avoid picking up the tail of the photon propagator + double q2 = q2min + ( i + 0.5 ) / n * ( q2max - q2min ); + // std::cout<<"m0 = "< prob ) - prob = probctl[1]; - if ( probctl[2] > prob ) - prob = probctl[2]; - - if ( fabs( c ) > 1e-20 ) { - double ctlx = -0.5 * b / c; - if ( fabs( ctlx ) < 1.0 ) { - double probtmp = a + b * ctlx + c * ctlx * ctlx; - if ( probtmp > prob ) - prob = probtmp; + for ( int j = 0; j < nj - 2; j++ ) { + // pm,p0,pp contain probabilities at x=-1,0,1. + // prob = a + b*x + c*x^2 + double pm = p[j + 0], p0 = p[j + 1], pp = p[j + 2]; + double a = p0, b = 0.5 * ( pp - pm ), c = 0.5 * ( pp + pm ) - p0; + prob = std::max( p0, std::max( pp, pm ) ); + if ( fabs( c ) > 1e-20 ) { + double x = -0.5 * b / c; + if ( fabs( x ) < 1.0 ) + prob = std::max( prob, a + b * x + c * x * x ); } } - //EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "prob,probctl:"< maxfoundprob ) { maxfoundprob = prob; } - - //cout << "q2,maxfoundprob:"<deleteTree(); + if ( 0 ) { + std::vector v; + float x, y; + std::ifstream INR( "eRe_wide.dat" ); + while ( INR >> x >> y ) + v.push_back( x ); + std::cout << "q2 points read -- " << v.size() << std::endl; + + std::ofstream OUT( "b2sllprob_amp.dat" ); + double m0 = mH[0]; + m0 = EvtPDL::getMinMass( meson ); + double q2max = ( m - m0 ) * ( m - m0 ); + //loop over q2 + for ( double q2 : v ) { + // want to avoid picking up the tail of the photon propagator + if ( !( q2min <= q2 && q2 < q2max ) ) + continue; + double Erho = ( m * m + m0 * m0 - q2 ) / ( 2.0 * m ), + Prho = sqrt( ( Erho - m0 ) * ( Erho + m0 ) ); + + p4meson.set( Erho, 0.0, 0.0, -Prho ); + p4w.set( m - Erho, 0.0, 0.0, Prho ); + + //This is in the W rest frame + double Elep = sqrt( q2 ) * 0.5, + Plep = sqrt( Elep * Elep - mass[1] * mass[1] ); + + const int nj = 3 + 2 + 4 + 8 + 32; + double cmin = -1, cmax = 1, dc = ( cmax - cmin ) / ( nj - 1 ); + double maxprob = 0; + for ( int j = 0; j < nj; j++ ) { + double c = cmin + dc * j; - poleSize = 0.04 * ( maxpole / maxfoundprob ) * 4 * ( mass[1] * mass[1] ); + //These are in the W rest frame. Need to boost out into the B frame. + double Py = Plep * sqrt( 1.0 - c * c ), Pz = Plep * c; + p4lepton1.set( Elep, 0.0, Py, Pz ); + p4lepton2.set( Elep, 0.0, -Py, -Pz ); - //poleSize=0.002; + p4lepton1 = boostTo( p4lepton1, p4w ); + p4lepton2 = boostTo( p4lepton2, p4w ); - //cout <<"maxfoundprob,maxpole,poleSize:"<init( meson, p4meson ); + lep1->init( lepton1, p4lepton1 ); + lep2->init( lepton2, p4lepton2 ); + CalcAmp( root_part, amp, FormFactors ); + double prob = rho.normalizedProb( amp.getSpinDensity() ); + maxprob = std::max( maxprob, prob ); + } + OUT << q2 << " " << maxprob << "\n"; + } + } + + root_part->deleteTree(); - maxfoundprob *= 1.15; + // poleSize = 0.04 * ( maxpole / maxfoundprob ) * q2min; + poleSize = ( maxpole / maxfoundprob ) * q2min; + std::cout << "max " << maxfoundprob << " " << maxpole << " " << poleSize + << std::endl; + // maxfoundprob *= 1.5e-3; return maxfoundprob; // with resonances + // maxfoundprob *= 4e-7; return maxfoundprob; + // maxfoundprob *= 3.7e-4; return maxfoundprob; // with resonances - return maxfoundprob; + return maxfoundprob * 3; // without resonances + return maxfoundprob / 1000 * 3.5; } EvtComplex EvtbTosllAmp::GetC7Eff( double q2, bool nnlo ) diff --git a/src/EvtGenModels/EvtbTosllBSZFF.cpp b/src/EvtGenModels/EvtbTosllBSZFF.cpp new file mode 100644 --- /dev/null +++ b/src/EvtGenModels/EvtbTosllBSZFF.cpp @@ -0,0 +1,130 @@ +#include "EvtGenModels/EvtbTosllBSZFF.hh" + +#include "EvtGenBase/EvtPDL.hh" + +#include + +inline double poly( double x, int n, const double* c ) +{ + double t = c[--n]; + while ( n-- ) + t = c[n] + x * t; + return t; +} + +void EvtbTosllBSZFF::getScalarFF( EvtId parent, EvtId daught, double t, + double /*mass*/, double& fp, double& f0, + double& ft ) +{ + double m = EvtPDL::getMeanMass( parent ); + t /= m * m; + + // B --> K form factors + //this is Ali-Ball '01 (or really Ali-Ball'99 minimum allowed) + fp = t * ( 1.568 + t * ( 0.470 + t * ( 0.885 ) ) ); + f0 = t * ( 0.740 + t * ( 0.080 + t * ( 0.425 ) ) ); + ft = t * ( 1.600 + t * ( 0.501 + t * ( 0.796 ) ) ); + + fp = 0.278 * exp( fp ); + f0 = 0.278 * exp( f0 ); + ft = 0.300 * exp( ft ); +} + +void EvtbTosllBSZFF::getVectorFF( EvtId pId, EvtId, double q2, double mV, + double& a1, double& a2, double& a0, double& v, + double& t1, double& t2, double& t3 ) +{ + const static double alfa[7][4] = { + // coefficients are from https://arxiv.org/src/1503.05534v3/anc/BKstar_LCSR-Lattice.json + // m_res, c0, c1, c2 + { 5.415000, 0.376313, -1.165970, 2.424430 }, // V + { 5.366000, 0.369196, -1.365840, 0.128191 }, // A0 + { 5.829000, 0.297250, 0.392378, 1.189160 }, // A1 + { 5.829000, 0.265375, 0.533638, 0.483166 }, // A12 + { 5.415000, 0.312055, -1.008930, 1.527200 }, // T1 + { 5.829000, 0.312055, 0.496846, 1.614310 }, // T2 + { 5.829000, 0.667412, 1.318120, 3.823340 } // T12 + // ffv1 + // {5.415000, 0.376313+0.033259, -1.165970+0.123496, 2.424430-0.038413}, + // {5.366000, 0.369196-0.000870, -1.365840+0.158886, 0.128191+0.705349}, + // {5.829000, 0.297250-0.023744, 0.392378-0.128900, 1.189160-0.207408}, + // {5.829000, 0.265375-0.008545, 0.533638+0.108205, 0.483166+0.435457}, + // {5.415000, 0.312055-0.001346, -1.008930-0.198232, 1.527200-2.387283}, + // {5.829000, 0.312055-0.038173, 0.496846-0.067891, 1.614310+0.349193}, + // {5.829000, 0.667412+0.036716, 1.318120-0.038069, 3.823340+0.606166}, + + // ffv2 + // {5.415000, 0.376313-0.041113, -1.165970-0.045535, 2.424430+1.137686}, + // {5.366000, 0.369196-0.003901, -1.365840-0.121212, 0.128191-0.729504}, + // {5.829000, 0.297250+0.048609, 0.392378+0.210965, 1.189160+0.906521}, + // {5.829000, 0.265375+0.000100, 0.533638-0.042370, 0.483166+0.284835}, + // {5.415000, 0.312055+0.002605, -1.008930+0.215231, 1.527200+2.604514}, + // {5.829000, 0.312055-0.003748, 0.496846-0.039576, 1.614310-0.241782}, + // {5.829000, 0.667412-0.007496, 1.318120-0.252147, 3.823340-1.210450}, + + // ffv3 + // {5.415000, 0.376313+0.000735, -1.165970-0.007790, 2.424430-1.608712}, + // {5.366000, 0.369196+0.050484, -1.365840+0.217869, 0.128191-0.479841}, + // {5.829000, 0.297250-0.031106, 0.392378-0.314969, 1.189160-2.128420}, + // {5.829000, 0.265375+0.027843, 0.533638+0.155898, 0.483166+0.413373}, + // {5.415000, 0.312055+0.006707, -1.008930-0.112003, 1.527200-1.779484}, + // {5.829000, 0.312055-0.021665, 0.496846-0.400001, 1.614310-2.018040}, + // {5.829000, 0.667412-0.078741, 1.318120-0.140305, 3.823340+3.390407}, + + // ffv4 + // {5.415000, 0.376313-0.015611, -1.165970-0.474344, 2.424430-1.821098}, + // {5.366000, 0.369196+0.051546, -1.365840+0.113835, 0.128191-0.334548}, + // {5.829000, 0.297250-0.009712, 0.392378-0.144326, 1.189160-0.880261}, + // {5.829000, 0.265375-0.006262, 0.533638-0.036838, 0.483166+0.208676}, + // {5.415000, 0.312055+0.012515, -1.008930+0.057467, 1.527200+1.067426}, + // {5.829000, 0.312055+0.009879, 0.496846-0.025578, 1.614310+0.063207}, + // {5.829000, 0.667412+0.088066, 1.318120+0.306504, 3.823340-1.852577}, + + // ffv5 + // {5.415000, 0.376313+0.009252, -1.165970-0.073415, 2.424430-1.479015}, + // {5.366000, 0.369196-0.037582, -1.365840-0.015347, 0.128191+1.217427}, + // {5.829000, 0.297250-0.039144, 0.392378-0.016291, 1.189160+1.097737}, + // {5.829000, 0.265375-0.003710, 0.533638-0.047800, 0.483166-0.066831}, + // {5.415000, 0.312055-0.012360, -1.008930+0.139269, 1.527200+1.496325}, + // {5.829000, 0.312055+0.002256, 0.496846+0.142894, 1.614310+1.115639}, + // {5.829000, 0.667412-0.089004, 1.318120-0.037696, 3.823340+2.278156}, + }; + + double mB = EvtPDL::getMeanMass( pId ); + double mBaV = mB + mV, mBsV = mB - mV; + double tp = mBaV * mBaV; // t_{+} = (m_B + m_V)^2 + double s0 = sqrt( + 2 * mBaV * + sqrt( mB * mV ) ); // sqrt(t_{+} - t_{+}*(1 - sqrt(1 - t_{-}/t_{+}))) + double z0 = ( mBaV - s0 ) / + ( mBaV + s0 ); // (sqrt(t_{+}) - s0)/(sqrt(t_{+}) + s0) + + double s = sqrt( tp - q2 ), z = ( s - s0 ) / ( s + s0 ), dz = z - z0, ff[7]; + for ( int j = 0; j < 7; j++ ) { + double mR = alfa[j][0], mR2 = mR * mR; + ff[j] = ( mR2 / ( mR2 - q2 ) ) * poly( dz, 3, alfa[j] + 1 ); + } + + // Källén-function + // arXiv:1503.05534 Eq. D.3 + double lambda = ( mBaV * mBaV - q2 ) * ( mBsV * mBsV - q2 ); + + v = ff[0]; + a0 = ff[1]; + a1 = ff[2]; + // Eq. D.5 arXiv:1503.05534 + // A12 = (mBaV*mBaV*(mBaV*mBsV - q2)*A1 - lambda(q2)*A2)/(16*mB*mV*mV*mBaV); + double a12 = ff[3]; + a2 = mBaV * + ( ( mBaV * ( mBaV * mBsV - q2 ) ) * a1 - ( 16 * mB * mV * mV ) * a12 ) / + lambda; + t1 = ff[4]; + t2 = ff[5]; + // Eq. D.5 arXiv:1503.05534 + // T23 = mBaV*mBsV*(mB*mB + 3*mV*mV - q2)*T2 - lambda(q2)*T3)/(8*mB*mV*mV*mBsV); + double t23 = ff[6]; + t3 = mBsV * + ( mBaV * ( mB * mB + 3 * mV * mV - q2 ) * t2 - + ( 8 * mB * mV * mV ) * t23 ) / + lambda; +} diff --git a/src/EvtGenModels/EvtbTosllBallFF.cpp b/src/EvtGenModels/EvtbTosllBallFF.cpp --- a/src/EvtGenModels/EvtbTosllBallFF.cpp +++ b/src/EvtGenModels/EvtbTosllBallFF.cpp @@ -30,10 +30,27 @@ _theFFModel = ffmodel; } +static EvtId IdKp, IdKm, IdKS, IdK0, IdaK0, IdKL, IdPip, IdPim, IdPi0, IdEta, + IdEtap; +static bool sfirst = true; void EvtbTosllBallFF::getScalarFF( EvtId parent, EvtId daught, double t, double /*mass*/, double& fp, double& f0, double& ft ) { + if ( sfirst ) { + sfirst = false; + IdKp = EvtPDL::getId( "K+" ); + IdKm = EvtPDL::getId( "K-" ); + IdKS = EvtPDL::getId( "K_S0" ); + IdK0 = EvtPDL::getId( "K0" ); + IdaK0 = EvtPDL::getId( "anti-K0" ); + IdKL = EvtPDL::getId( "K_L0" ); + IdPip = EvtPDL::getId( "pi+" ); + IdPim = EvtPDL::getId( "pi-" ); + IdPi0 = EvtPDL::getId( "pi0" ); + IdEta = EvtPDL::getId( "eta" ); + IdEtap = EvtPDL::getId( "eta'" ); + } int model = _theFFModel; double m = EvtPDL::getMeanMass( parent ); @@ -43,12 +60,8 @@ double shat2 = shat * shat; double shat3 = shat2 * shat; - if ( daught == EvtPDL::getId( std::string( "K+" ) ) || - daught == EvtPDL::getId( std::string( "K-" ) ) || - daught == EvtPDL::getId( std::string( "K_S0" ) ) || - daught == EvtPDL::getId( std::string( "K0" ) ) || - daught == EvtPDL::getId( std::string( "anti-K0" ) ) || - daught == EvtPDL::getId( std::string( "K_L0" ) ) ) { + if ( daught == IdKp || daught == IdKm || daught == IdKS || daught == IdK0 || + daught == IdaK0 || daught == IdKL ) { // B --> K form factors if ( model == 1 ) { //this is Ali-Ball '01 (or really Ali-Ball'99 minimum allowed) @@ -97,9 +110,7 @@ ft = ( 0.1851 / ( 1. - ( t / 29.30 ) ) ) + ( 0.1905 / ( 1. - ( t / 29.30 ) ) / ( 1. - ( t / 29.30 ) ) ); } - } else if ( daught == EvtPDL::getId( std::string( "pi+" ) ) || - daught == EvtPDL::getId( std::string( "pi-" ) ) || - daught == EvtPDL::getId( std::string( "pi0" ) ) ) { + } else if ( daught == IdPip || daught == IdPim || daught == IdPi0 ) { if ( model == 1 ) { // B --> pi form factors from Ball-Zwicky'01 (tabulated in hep-ph/0306251) fp = 0.261 / ( 1. - 2.03 * shat + 1.293 * shat * shat ); @@ -155,8 +166,7 @@ ft = ( 0.152 / ( 1. - ( t / ( 5.32 * 5.32 ) ) ) ) + ( 0.122 / ( 1. - ( t / 28.40 ) ) / ( 1. - ( t / 28.40 ) ) ); } - } else if ( daught == EvtPDL::getId( std::string( "eta" ) ) || - daught == EvtPDL::getId( std::string( "eta'" ) ) ) { + } else if ( daught == IdEta || daught == IdEtap ) { if ( model == 1 ) { // B --> eta form factors fp = 0.261 / ( 1. - 2.03 * shat + 1.293 * shat * shat ); @@ -185,11 +195,30 @@ // <<"ft "< K* form factors - daught == EvtPDL::getId( std::string( "K*0" ) ) || - daught == EvtPDL::getId( std::string( "anti-K*0" ) ) ) { + if ( parent == IdBs || parent == IdaBs ) { // B_s decay form factors + if ( daught == IdKst0 || daught == IdaKst0 ) { // B_s --> K* form factors if ( model == 6 ) { // Ball-Zwicky LCSR '05 (mb = 4.8) a1 = 0.231 / ( 1. - ( t / 32.94 ) ); @@ -228,9 +251,7 @@ } } - else if ( - // B_s --> phi form factors - daught == EvtPDL::getId( std::string( "phi" ) ) ) { + else if ( daught == IdPhi ) { // B_s --> phi form factors if ( model == 6 ) { // Ball-Zwicky LCSR '05 (mb = 4.8) a1 = 0.308 / ( 1. - ( t / 36.54 ) ); @@ -253,14 +274,8 @@ } } } - } else if ( daught == EvtPDL::getId( std::string( "K*+" ) ) || - daught == EvtPDL::getId( std::string( "K*-" ) ) || - daught == EvtPDL::getId( std::string( "K*0" ) ) || - daught == EvtPDL::getId( std::string( "anti-K*0" ) ) || - daught == EvtPDL::getId( std::string( "K_1+" ) ) || - daught == EvtPDL::getId( std::string( "K_1-" ) ) - - ) { + } else if ( daught == IdKstp || daught == IdKstm || daught == IdKst0 || + daught == IdaKst0 || daught == IdK1p || daught == IdK1m ) { if ( model == 1 ) { //this is Ali-Ball '01 a1 = 0.294 * exp( 0.656 * shat + 0.456 * shat2 ); @@ -372,9 +387,7 @@ ( 1 + ( -10.3 ) * ( z - z0 + 0.5 * ( ( z * z ) - ( z0 * z0 ) ) ) ); } - } else if ( daught == EvtPDL::getId( std::string( "rho+" ) ) || - daught == EvtPDL::getId( std::string( "rho-" ) ) || - daught == EvtPDL::getId( std::string( "rho0" ) ) ) { + } else if ( daught == IdRhop || daught == IdRhom || daught == IdRho0 ) { if ( model == 1 ) { // B --> rho form factors a1 = 0.261 / ( 1. - 0.29 * shat - 0.415 * shat * shat ); @@ -406,7 +419,7 @@ t3 = ( m * m - md * md ) * ( t3tilde - t2 ) / t; } } - } else if ( daught == EvtPDL::getId( std::string( "omega" ) ) ) { + } else if ( daught == IdOmega ) { if ( model == 1 ) { // B --> rho form factors a1 = 0.261 / ( 1. - 0.29 * shat - 0.415 * shat * shat ); diff --git a/src/EvtGenModels/EvtbTosllNP.cpp b/src/EvtGenModels/EvtbTosllNP.cpp new file mode 100644 --- /dev/null +++ b/src/EvtGenModels/EvtbTosllNP.cpp @@ -0,0 +1,142 @@ + +/*********************************************************************** +* 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 "EvtGenModels/EvtbTosllNP.hh" + +#include "EvtGenBase/EvtGenKine.hh" +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtParticle.hh" + +#include "EvtGenModels/EvtbTosllAmp.hh" +#include "EvtGenModels/EvtbTosllBSZFF.hh" +#include "EvtGenModels/EvtbTosllVectorAmpNP.hh" + +using std::endl; + +std::string EvtbTosllNP::getName() +{ + return "BTOSLLNP"; +} + +EvtDecayBase* EvtbTosllNP::clone() +{ + return new EvtbTosllNP; +} + +void EvtbTosllNP::decay( EvtParticle* p ) +{ + setWeight( p->initializePhaseSpace( getNDaug(), getDaugs(), false, + _poleSize, 1, 2 ) ); + _calcamp->CalcAmp( p, _amp2, _ffmodel.get() ); +} + +void EvtbTosllNP::initProbMax() +{ + EvtId pId = getParentId(), mId = getDaug( 0 ), l1Id = getDaug( 1 ), + l2Id = getDaug( 2 ); + //This routine sets the _poleSize. + double mymaxprob = _calcamp->CalcMaxProb( pId, mId, l1Id, l2Id, + _ffmodel.get(), _poleSize ); + setProbMax( mymaxprob ); +} + +void EvtbTosllNP::init() +{ + // First choose format of NP coefficients from the .DEC file + // Cartesian(x,y)(0) or Polar(R,phase)(1) + int n = getNArg(); + if ( !( n == 0 || ( n - 1 ) % 3 == 0 ) ) { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Error in parameters in the BTOSLLNP decay model." << endl; + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Will terminate execution!" << endl; + ::abort(); + } + + checkNDaug( 3 ); + + //We expect the parent to be a scalar + //and the daughters to be K* lepton+ lepton- or K lepton+ lepton- + + checkSpinParent( EvtSpinType::SCALAR ); + checkSpinDaughter( 1, EvtSpinType::DIRAC ); + checkSpinDaughter( 2, EvtSpinType::DIRAC ); + + EvtId mId = getDaug( 0 ); + EvtId IdKst0 = EvtPDL::getId( "K*0" ), IdaKst0 = EvtPDL::getId( "anti-K*0" ), + IdKstp = EvtPDL::getId( "K*+" ), IdKstm = EvtPDL::getId( "K*-" ); + EvtId IdK0 = EvtPDL::getId( "K0" ), IdaK0 = EvtPDL::getId( "anti-K0" ), + IdKp = EvtPDL::getId( "K+" ), IdKm = EvtPDL::getId( "K-" ); + + if ( mId != IdKst0 && mId != IdaKst0 && mId != IdKstp && mId != IdKstm && + mId != IdK0 && mId != IdaK0 && mId != IdKp && mId != IdKm ) { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "EvtbTosllNP generator expected a K(*) 1st daughter, found: " + << EvtPDL::name( mId ) << endl; + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Will terminate execution!" << endl; + ::abort(); + } + + _ffmodel = std::make_unique(); + _calcamp = std::make_unique(); + + auto getInteger = [this]( int i ) -> int { + double a = getArg( i ); + if ( a - static_cast( a ) != 0 ) { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Parameters is not integer in the BTOSLLNP decay model: " << i + << " " << a << endl; + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Will terminate execution!" << endl; + ::abort(); + } + return static_cast( a ); + }; + EvtbTosllVectorAmpNP* amp = static_cast( + _calcamp.get() ); + if ( n > 0 ) { // parse arguments + int i = 0, coordsyst = getInteger( i++ ); + while ( i < n ) { + int id = getInteger( i++ ); // New Physics cooeficient Id + double a0 = getArg( i++ ); + double a1 = getArg( i++ ); + EvtComplex c = ( coordsyst ) + ? EvtComplex( a0 * cos( a1 ), a0 * sin( a1 ) ) + : EvtComplex( a0, a1 ); + if ( id == 0 ) + amp->m_dc7 = c; // delta C_7eff + if ( id == 1 ) + amp->m_dc9 = c; // delta C_9eff + if ( id == 2 ) + amp->m_dc10 = c; // delta C_10eff + if ( id == 3 ) + amp->m_c7p = c; // C'_7eff -- right hand polarizations + if ( id == 4 ) + amp->m_c9p = c; // C'_9eff -- right hand polarizations + if ( id == 5 ) + amp->m_c10p = c; // c'_10eff -- right hand polarizations + if ( id == 6 ) + amp->m_cS = c; // (C_S - C'_S) -- scalar right and left polarizations + if ( id == 7 ) + amp->m_cP = c; // (C_P - C'_P) -- pseudo-scalar right and left polarizations + } + } +} diff --git a/src/EvtGenModels/EvtbTosllNPR.cpp b/src/EvtGenModels/EvtbTosllNPR.cpp new file mode 100644 --- /dev/null +++ b/src/EvtGenModels/EvtbTosllNPR.cpp @@ -0,0 +1,405 @@ + +/*********************************************************************** +* 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 "EvtGenModels/EvtbTosllNPR.hh" + +#include "EvtGenBase/EvtGenKine.hh" +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtScalarParticle.hh" + +#include "EvtGenModels/EvtbTosllAmp.hh" +#include "EvtGenModels/EvtbTosllBSZFF.hh" +#include "EvtGenModels/EvtbTosllVectorAmpNP.hh" + +#include +#include +using std::endl; + +std::string EvtbTosllNPR::getName() +{ + return "BTOSLLNPR"; +} + +EvtDecayBase* EvtbTosllNPR::clone() +{ + return new EvtbTosllNPR; +} + +void EvtbTosllNPR::decay( EvtParticle* p ) +{ + static EvtVector4R p4[3]; + static double mass[3]; + double m_b = p->mass(); + for ( int i = 0; i < 3; i++ ) + mass[i] = p->getDaug( i )->mass(); + EvtId* daughters = getDaugs(); + double weight = EvtGenKine::PhaseSpacePole2( m_b, mass[2], mass[1], mass[0], + p4, _ls ); + p->getDaug( 0 )->init( daughters[0], p4[2] ); + p->getDaug( 1 )->init( daughters[1], p4[1] ); + p->getDaug( 2 )->init( daughters[2], p4[0] ); + setWeight( weight ); + _calcamp->CalcAmp( p, _amp2, _ffmodel.get() ); +} + +struct ABW_t { + double Mv2, MGv, A; + + ABW_t( double m, double gtot, double, double bll ) + { + const double alpha = 1 / 137.035999084; // (+-21) PDG 2021 + Mv2 = m * m; + MGv = m * gtot; + A = 3 * M_PI / ( alpha * alpha ) * bll * gtot; + } + + std::complex R( double s ) const + { + return A / std::complex( s - Mv2, MGv ); + } + + void addknots( std::vector& en ) const + { + double Mv = sqrt( Mv2 ), Gv = MGv / Mv; + double w = 0.6; + double smin = ( Mv - w * Gv ) * ( Mv - w * Gv ); + double smax = ( Mv + w * Gv ) * ( Mv + w * Gv ); + en.push_back( Mv * Mv ); + for ( int i = 0, n = 30; i < n; i++ ) { + double s = smin + ( smax - smin ) * ( i + 0.5 ) / n; + en.push_back( s ); + } + w = 2; + double smin1 = ( Mv - w * Gv ) * ( Mv - w * Gv ); + double smax1 = ( Mv + w * Gv ) * ( Mv + w * Gv ); + for ( int i = 0, n = 30; i < n; i++ ) { + double s = smin1 + ( smax1 - smin1 ) * ( i + 0.5 ) / n; + if ( s >= smin && s <= smax ) + continue; + en.push_back( s ); + } + w = 8; + double smin2 = ( Mv - w * Gv ) * ( Mv - w * Gv ); + double smax2 = ( Mv + w * Gv ) * ( Mv + w * Gv ); + for ( int i = 0, n = 30; i < n; i++ ) { + double s = smin2 + ( smax2 - smin2 ) * ( i + 0.5 ) / n; + if ( s >= smin1 && s <= smax1 ) + continue; + en.push_back( s ); + } + w = 30; + double smin3 = ( Mv - w * Gv ) * ( Mv - w * Gv ); + double smax3 = ( Mv + w * Gv ) * ( Mv + w * Gv ); + for ( int i = 0, n = 30; i < n; i++ ) { + double s = smin3 + ( smax3 - smin3 ) * ( i + 0.5 ) / n; + if ( s >= smin2 && s <= smax2 ) + continue; + en.push_back( s ); + } + w = 100; + double smin4 = ( Mv - w * Gv ) * ( Mv - w * Gv ); + double smax4 = ( Mv + w * Gv ) * ( Mv + w * Gv ); + for ( int i = 0, n = 20; i < n; i++ ) { + double s = smin4 + ( smax4 - smin4 ) * ( i + 0.5 ) / n; + if ( s >= smin3 && s <= smax3 ) + continue; + en.push_back( s ); + } + } +}; + +ABW_t ajpsi( 3.0969, 0.0926 / 1000, 0.0926 / 1000 * 0.877, 0.05971 ); +ABW_t apsi2s( 3.6861, 0.294 / 1000, 0.294 / 1000 * 0.9785, 0.00793 ); +ABW_t apsi3770( 3773.7 / 1000, 27.2 / 1000, 27.2 / 1000, 1e-5 ); +ABW_t apsi4040( 4039 / 1000., 80 / 1000., 80 / 1000., 1.07e-5 ); +ABW_t apsi4160( 4191 / 1000., 70 / 1000., 70 / 1000., 6.9e-6 ); +ABW_t apsi4230( 4220 / 1000., 60 / 1000., 60 / 1000., 1e-5 ); + +std::complex aR( double s ) +{ + return ajpsi.R( s ) + apsi2s.R( s ) + apsi3770.R( s ) + apsi4040.R( s ) + + apsi4160.R( s ) + apsi4230.R( s ); +} + +std::vector getgrid() +{ + const double m_c = 1.3, m_s = 0.2, m_e = 0.511e-3, MD0 = 1864.84 / 1000; + + std::vector v = { ajpsi, apsi2s, apsi3770, + apsi4040, apsi4160, apsi4230 }; + std::vector en; + for ( const auto& t : v ) + t.addknots( en ); + + double smax = 25; //4.37*4.37; + for ( unsigned i = 0, n = 200; i < n; i++ ) { + double s = smax / n * ( i + 0.5 ); + en.push_back( s ); + } + en.push_back( ajpsi.Mv2 + 0.1 ); + en.push_back( ajpsi.Mv2 - 0.1 ); + en.push_back( ajpsi.Mv2 + 0.08 ); + en.push_back( ajpsi.Mv2 - 0.08 ); + en.push_back( ajpsi.Mv2 + 0.065 ); + en.push_back( ajpsi.Mv2 - 0.065 ); + + double t0; + for ( t0 = 4 * m_e * m_e; t0 < 0.15; t0 *= 1.5 ) { + en.push_back( t0 ); + } + + for ( t0 = 4 * m_e * m_e; t0 < 0.5; t0 *= 1.5 ) { + en.push_back( 4 * MD0 * MD0 + t0 ); + en.push_back( 4 * MD0 * MD0 - t0 ); + } + + en.push_back( 4 * m_c * m_c ); + for ( t0 = 0.00125; t0 < 0.05; t0 *= 1.5 ) { + en.push_back( 4 * m_c * m_c + t0 ); + en.push_back( 4 * m_c * m_c - t0 ); + } + + en.push_back( 4 * m_s * m_s ); + for ( t0 = 0.00125; t0 < 0.1; t0 *= 1.5 ) { + en.push_back( 4 * m_s * m_s + t0 ); + en.push_back( 4 * m_s * m_s - t0 ); + } + + std::sort( en.begin(), en.end() ); + + for ( std::vector::iterator it = en.begin(); it != en.end(); it++ ) { + if ( *it > smax ) { + en.erase( it, en.end() ); + break; + } + } + return en; +} + +void EvtbTosllNPR::initProbMax() +{ + EvtId pId = getParentId(), mId = getDaug( 0 ), l1Id = getDaug( 1 ), + l2Id = getDaug( 2 ); + // //This routine sets the _poleSize. + // double mymaxprob = _calcamp->CalcMaxProb( pId, mId, l1Id, l2Id, _ffmodel.get(), _poleSize ); + // setProbMax( mymaxprob ); + + std::vector v = getgrid(); + std::vector pp; + EvtScalarParticle* scalar_part; + EvtParticle* root_part; + + scalar_part = new EvtScalarParticle; + + //cludge to avoid generating random numbers! + scalar_part->noLifeTime(); + + EvtVector4R p_init; + + p_init.set( EvtPDL::getMass( pId ), 0.0, 0.0, 0.0 ); + scalar_part->init( pId, p_init ); + root_part = (EvtParticle*)scalar_part; + root_part->setDiagonalSpinDensity(); + + EvtParticle *daughter, *lep1, *lep2; + + EvtAmp amp; + + EvtId listdaug[3]; + listdaug[0] = mId; + listdaug[1] = l1Id; + listdaug[2] = l2Id; + + amp.init( pId, 3, listdaug ); + + root_part->makeDaughters( 3, listdaug ); + daughter = root_part->getDaug( 0 ); + lep1 = root_part->getDaug( 1 ); + lep2 = root_part->getDaug( 2 ); + + //cludge to avoid generating random numbers! + daughter->noLifeTime(); + lep1->noLifeTime(); + lep2->noLifeTime(); + + //Initial particle is unpolarized, well it is a scalar so it is trivial + EvtSpinDensity rho; + rho.setDiag( root_part->getSpinStates() ); + + double mass[3]; + + double m = root_part->mass(); + + EvtVector4R p4meson, p4lepton1, p4lepton2, p4w; + + double maxprobfound = 0.0; + + mass[1] = EvtPDL::getMeanMass( l1Id ); + mass[2] = EvtPDL::getMeanMass( l2Id ); + std::vector mH; + mH.push_back( EvtPDL::getMeanMass( mId ) ); + //if the particle is narrow dont bother with changing the mass. + double g = EvtPDL::getWidth( mId ); + if ( g > 0 ) { + mH.push_back( EvtPDL::getMinMass( mId ) ); + mH.push_back( + std::min( EvtPDL::getMaxMass( mId ), m - mass[1] - mass[2] ) ); + mH.push_back( mH.front() - g ); + mH.push_back( mH.front() + g ); + } + + double q2min = ( mass[1] + mass[2] ) * ( mass[1] + mass[2] ); + + double m0 = EvtPDL::getMinMass( mId ); + double q2max = ( m - m0 ) * ( m - m0 ); + m0 = mH[0]; + //loop over q2 + for ( double q2 : v ) { + // want to avoid picking up the tail of the photon propagator + if ( !( q2min <= q2 && q2 < q2max ) ) + continue; + double Erho = ( m * m + m0 * m0 - q2 ) / ( 2.0 * m ), + Prho = sqrt( ( Erho - m0 ) * ( Erho + m0 ) ); + + p4meson.set( Erho, 0.0, 0.0, -Prho ); + p4w.set( m - Erho, 0.0, 0.0, Prho ); + + //This is in the W rest frame + double Elep = sqrt( q2 ) * 0.5, + Plep = sqrt( Elep * Elep - mass[1] * mass[1] ); + + const int nj = 3 + 2 + 4 + 8 + 32; + double cmin = -1, cmax = 1, dc = ( cmax - cmin ) / ( nj - 1 ); + double maxprob = 0; + for ( int j = 0; j < nj; j++ ) { + double c = cmin + dc * j; + + //These are in the W rest frame. Need to boost out into the B frame. + double Py = Plep * sqrt( 1.0 - c * c ), Pz = Plep * c; + p4lepton1.set( Elep, 0.0, Py, Pz ); + p4lepton2.set( Elep, 0.0, -Py, -Pz ); + + p4lepton1 = boostTo( p4lepton1, p4w ); + p4lepton2 = boostTo( p4lepton2, p4w ); + + //Now initialize the daughters... + daughter->init( mId, p4meson ); + lep1->init( l1Id, p4lepton1 ); + lep2->init( l2Id, p4lepton2 ); + _calcamp->CalcAmp( root_part, amp, _ffmodel.get() ); + double prob = rho.normalizedProb( amp.getSpinDensity() ); + maxprob = std::max( maxprob, prob ); + } + // printf("%f %f\n",q2,maxprob); + pp.push_back( { (float)q2, (float)maxprob } ); + maxprobfound = std::max( maxprobfound, maxprob ); + } + root_part->deleteTree(); + + _ls.init( pp ); + maxprobfound *= 8e-8 * 1.1; + setProbMax( maxprobfound ); +} + +void EvtbTosllNPR::init() +{ + // First choose format of NP coefficients from the .DEC file + // Cartesian(x,y)(0) or Polar(R,phase)(1) + int n = getNArg(); + checkNDaug( 3 ); + + //We expect the parent to be a scalar + //and the daughters to be K* lepton+ lepton- + + checkSpinParent( EvtSpinType::SCALAR ); + checkSpinDaughter( 1, EvtSpinType::DIRAC ); + checkSpinDaughter( 2, EvtSpinType::DIRAC ); + + EvtId mId = getDaug( 0 ); + EvtId IdKst0 = EvtPDL::getId( "K*0" ), IdaKst0 = EvtPDL::getId( "anti-K*0" ), + IdKstp = EvtPDL::getId( "K*+" ), IdKstm = EvtPDL::getId( "K*-" ); + EvtId IdK0 = EvtPDL::getId( "K0" ), IdaK0 = EvtPDL::getId( "anti-K0" ), + IdKp = EvtPDL::getId( "K+" ), IdKm = EvtPDL::getId( "K-" ); + if ( mId != IdKst0 && mId != IdaKst0 && mId != IdKstp && mId != IdKstm && + mId != IdK0 && mId != IdaK0 && mId != IdKp && mId != IdKm ) { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "EvtbTosllNPR generator expected a K(*) 1st daughter, found: " + << EvtPDL::name( getDaug( 0 ) ) << endl; + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Will terminate execution!" << endl; + ::abort(); + } + + _ffmodel = std::make_unique(); + _calcamp = std::make_unique(); + + auto getInteger = [this]( int i ) -> int { + double a = getArg( i ); + if ( a - static_cast( a ) != 0 ) { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Parameters is not integer in the BTOSLLNPR decay model: " << i + << " " << a << endl; + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Will terminate execution!" << endl; + ::abort(); + } + return static_cast( a ); + }; + EvtbTosllVectorAmpNP* amp = static_cast( + _calcamp.get() ); + double phi = 0; + if ( n > 0 ) { // parse arguments + int i = 0, coordsyst = getInteger( i++ ); + auto getcomplex = [this, &i, coordsyst]() -> EvtComplex { + double a0 = getArg( i++ ); + double a1 = getArg( i++ ); + return ( coordsyst ) ? EvtComplex( a0 * cos( a1 ), a0 * sin( a1 ) ) + : EvtComplex( a0, a1 ); + }; + auto getreal = [this, &i]() -> double { return getArg( i++ ); }; + while ( i < n ) { + int id = getInteger( i++ ); // New Physics cooeficient Id + if ( id == 0 ) + amp->m_dc7 = getcomplex(); // delta C_7eff + if ( id == 1 ) + amp->m_dc9 = getcomplex(); // delta C_9eff + if ( id == 2 ) + amp->m_dc10 = getcomplex(); // delta C_10eff + if ( id == 3 ) + amp->m_c7p = getcomplex(); // C'_7eff -- right hand polarizations + if ( id == 4 ) + amp->m_c9p = getcomplex(); // C'_9eff -- right hand polarizations + if ( id == 5 ) + amp->m_c10p = getcomplex(); // c'_10eff -- right hand polarizations + if ( id == 6 ) + amp->m_cS = getcomplex(); // (C_S - C'_S) -- scalar right and left polarizations + if ( id == 7 ) + amp->m_cP = getcomplex(); // (C_P - C'_P) -- pseudo-scalar right and left polarizations + if ( id == 10 ) + phi = getreal(); + } + } + + std::vector v = getgrid(); + std::complex eiphi = std::complex( cos( phi ), sin( phi ) ); + for ( auto t : v ) + amp->m_reso.push_back( std::make_pair( t, -aR( t ) * eiphi ) ); +} diff --git a/src/EvtGenModels/EvtbTosllVectorAmpNP.cpp b/src/EvtGenModels/EvtbTosllVectorAmpNP.cpp new file mode 100644 --- /dev/null +++ b/src/EvtGenModels/EvtbTosllVectorAmpNP.cpp @@ -0,0 +1,396 @@ + +/*********************************************************************** +* 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 "EvtGenModels/EvtbTosllVectorAmpNP.hh" + +#include "EvtGenBase/EvtAmp.hh" +#include "EvtGenBase/EvtDiracSpinor.hh" +#include "EvtGenBase/EvtId.hh" +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtTensor4C.hh" +#include "EvtGenBase/EvtVector4C.hh" + +#include "EvtGenModels/EvtbTosllAmp.hh" +#include "EvtGenModels/EvtbTosllFF.hh" + +#include +#include + +static std::complex h( double q2, double mq ) +{ + const double mu = 4.8; // mu = m_b = 4.8 GeV + if ( mq == 0.0 ) + return std::complex( 8. / 27 + 4. / 9 * log( mu * mu / q2 ), + 4. / 9 * M_PI ); + double z = 4 * mq * mq / q2; + std::complex t; + if ( z > 1 ) { + t = atan( 1 / sqrt( z - 1 ) ); + } else { + t = std::complex( log( ( 1 + sqrt( 1 - z ) ) / sqrt( z ) ), + -M_PI / 2 ); + } + return -4 / 9. * ( 2 * log( mq / mu ) - 2 / 3. - z ) - + ( 4 / 9. * ( 2 + z ) * sqrt( fabs( z - 1 ) ) ) * t; +} + +static EvtComplex C9( double q2, + std::complex hReso = std::complex( 0, 0 ) ) +{ + double m_c = 1.3, m_b = 4.8; // quark masses + + double C1 = -0.257, C2 = 1.009, C3 = -0.005, C4 = -0.078, C5 = 0, C6 = 0.001; + + std::complex Y = ( h( q2, m_c ) + hReso ) * + ( 4 / 3. * C1 + C2 + 6 * C3 + 60 * C5 ); + Y -= 0.5 * h( q2, m_b ) * ( 7 * C3 + 4 / 3. * C4 + 76 * C5 + 64 / 3. * C6 ); + Y -= 0.5 * h( q2, 0.0 ) * ( C3 + 4 / 3. * C4 + 16 * C5 + 64 / 3. * C6 ); + Y += 4 / 3. * C3 + 64 / 9. * C5 + 64 / 27. * C6; + return EvtComplex( 4.211 + real( Y ), imag( Y ) ); +} + +static std::complex interpol( + const std::vector>>& v, double s ) +{ + if ( !v.size() ) + return 0; + int j = std::lower_bound( v.begin(), v.end(), s, + []( const std::pair>& a, + float b ) { return a.first < b; } ) - + v.begin(); + + double dx = v[j].first - v[j - 1].first; + auto dy = v[j].second - v[j - 1].second; + return v[j - 1].second + ( s - v[j - 1].first ) * ( dy / dx ); +} + +static EvtId IdB0, IdaB0, IdBp, IdBm, IdBs, IdaBs, IdRhop, IdRhom, IdRho0, + IdOmega, IdKst0, IdaKst0, IdKstp, IdKstm, IdK0, IdaK0, IdKp, IdKm; +static bool cafirst = true; +void EvtbTosllVectorAmpNP::CalcAmp( EvtParticle* parent, EvtAmp& amp, + EvtbTosllFF* formFactors ) +{ + if ( cafirst ) { + cafirst = false; + IdB0 = EvtPDL::getId( "B0" ); + IdaB0 = EvtPDL::getId( "anti-B0" ); + IdBp = EvtPDL::getId( "B+" ); + IdBm = EvtPDL::getId( "B-" ); + IdBs = EvtPDL::getId( "B_s0" ); + IdaBs = EvtPDL::getId( "anti-B_s0" ); + IdRhop = EvtPDL::getId( "rho+" ); + IdRhom = EvtPDL::getId( "rho-" ); + IdRho0 = EvtPDL::getId( "rho0" ); + IdOmega = EvtPDL::getId( "omega" ); + IdKst0 = EvtPDL::getId( "K*0" ); + IdaKst0 = EvtPDL::getId( "anti-K*0" ); + IdKstp = EvtPDL::getId( "K*+" ); + IdKstm = EvtPDL::getId( "K*-" ); + IdK0 = EvtPDL::getId( "K0" ); + IdaK0 = EvtPDL::getId( "anti-K0" ); + IdKp = EvtPDL::getId( "K+" ); + IdKm = EvtPDL::getId( "K-" ); + } + + EvtId dId = parent->getDaug( 0 )->getId(); + if ( dId == IdKst0 || dId == IdaKst0 || dId == IdKstp || dId == IdKstm ) { + CalcVAmp( parent, amp, formFactors ); + } + if ( dId == IdK0 || dId == IdaK0 || dId == IdKp || dId == IdKm ) { + CalcSAmp( parent, amp, formFactors ); + } +} + +void EvtbTosllVectorAmpNP::CalcVAmp( EvtParticle* parent, EvtAmp& amp, + EvtbTosllFF* formFactors ) +{ + // Add the lepton and neutrino 4 momenta to find q2 + EvtId pId = parent->getId(); + EvtId dId = parent->getDaug( 0 )->getId(); + + EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4(); + double q2 = q.mass2(); + double mesonmass = parent->getDaug( 0 )->mass(); + + double a1, a2, a0, v, t1, t2, t3; // form factors + formFactors->getVectorFF( pId, dId, q2, mesonmass, a1, a2, a0, v, t1, t2, t3 ); + + const EvtVector4R& p4meson = parent->getDaug( 0 )->getP4(); + double pmass = parent->mass(), ipmass = 1 / pmass; + const EvtVector4R p4b( pmass, 0.0, 0.0, 0.0 ); + EvtVector4R pbhat = p4b * ipmass; + EvtVector4R qhat = q * ipmass; + EvtVector4R pkstarhat = p4meson * ipmass; + EvtVector4R phat = pbhat + pkstarhat; + + EvtComplex c7 = -0.304; + EvtComplex c9 = C9( q2, interpol( m_reso, q2 ) ); + EvtComplex c10 = -4.103; + c7 += m_dc7; + c9 += m_dc9; + c10 += m_dc10; + + EvtComplex I( 0.0, 1.0 ); + + double mb = 4.8 /*GeV/c^2*/ * ipmass, ms = 0.093 /*GeV/c^2*/ * ipmass; + double mH = mesonmass * ipmass, oamH = 1 + mH, osmH = 1 - mH, + osmH2 = oamH * osmH, iosmH2 = 1 / osmH2; // mhatkstar + double is = pmass * pmass / q2; // 1/shat + a1 *= oamH; + a2 *= osmH; + a0 *= 2 * mH; + double cs0 = ( a1 - a2 - a0 ) * is; + a2 *= iosmH2; + v *= 2 / oamH; + + EvtComplex a = ( c9 + m_c9p ) * v + ( c7 + m_c7p ) * ( 4 * mb * is * t1 ); + EvtComplex b = ( c9 - m_c9p ) * a1 + + ( c7 - m_c7p ) * ( 2 * mb * is * osmH2 * t2 ); + EvtComplex c = ( c9 - m_c9p ) * a2 + + ( c7 - m_c7p ) * ( 2 * mb * ( t3 * iosmH2 + t2 * is ) ); + EvtComplex d = ( c9 - m_c9p ) * cs0 - ( c7 - m_c7p ) * ( 2 * mb * is * t3 ); + EvtComplex e = ( c10 + m_c10p ) * v; + EvtComplex f = ( c10 - m_c10p ) * a1; + EvtComplex g = ( c10 - m_c10p ) * a2; + EvtComplex h = ( c10 - m_c10p ) * cs0; + double sscale = a0 / ( mb + ms ); + EvtComplex hs = m_cS * sscale, hp = m_cP * sscale; + + int charge1 = EvtPDL::chg3( parent->getDaug( 1 )->getId() ); + + EvtParticle* lepPos = ( charge1 > 0 ) ? parent->getDaug( 1 ) + : parent->getDaug( 2 ); + EvtParticle* lepNeg = ( charge1 < 0 ) ? parent->getDaug( 1 ) + : parent->getDaug( 2 ); + + EvtDiracSpinor lp0( lepPos->spParent( 0 ) ), lp1( lepPos->spParent( 1 ) ); + EvtDiracSpinor lm0( lepNeg->spParent( 0 ) ), lm1( lepNeg->spParent( 1 ) ); + + EvtVector4C l11, l12, l21, l22, a11, a12, a21, a22; + EvtComplex s11, s12, s21, s22, p11, p12, p21, p22; + EvtTensor4C tt0( EvtGenFunctions::asymProd( pbhat, pkstarhat ) ); + + EvtTensor4C T1( tt0 ), T2( tt0 ); + const EvtTensor4C& gmn = EvtTensor4C::g(); + EvtTensor4C tt1( EvtGenFunctions::directProd( pbhat, phat ) ); + EvtTensor4C tt2( EvtGenFunctions::directProd( pbhat, qhat ) ); + + b *= I; + c *= I; + d *= I; + f *= I; + g *= I; + h *= I; + if ( pId == IdBm || pId == IdaB0 || pId == IdaBs ) { + T1 *= a; + T1.addScaled( -b, gmn ); + T1.addScaled( c, tt1 ); + T1.addScaled( d, tt2 ); + T2 *= e; + T2.addScaled( -f, gmn ); + T2.addScaled( g, tt1 ); + T2.addScaled( h, tt2 ); + + EvtLeptonVandACurrents( l11, a11, lp0, lm0 ); + EvtLeptonVandACurrents( l21, a21, lp1, lm0 ); + EvtLeptonVandACurrents( l12, a12, lp0, lm1 ); + EvtLeptonVandACurrents( l22, a22, lp1, lm1 ); + + s11 = EvtLeptonSCurrent( lp0, lm0 ); + p11 = EvtLeptonPCurrent( lp0, lm0 ); + s21 = EvtLeptonSCurrent( lp1, lm0 ); + p21 = EvtLeptonPCurrent( lp1, lm0 ); + s12 = EvtLeptonSCurrent( lp0, lm1 ); + p12 = EvtLeptonPCurrent( lp0, lm1 ); + s22 = EvtLeptonSCurrent( lp1, lm1 ); + p22 = EvtLeptonPCurrent( lp1, lm1 ); + } else if ( pId == IdBp || pId == IdB0 || pId == IdBs ) { + T1 *= -a; + T1.addScaled( -b, gmn ); + T1.addScaled( c, tt1 ); + T1.addScaled( d, tt2 ); + T2 *= -e; + T2.addScaled( -f, gmn ); + T2.addScaled( g, tt1 ); + T2.addScaled( h, tt2 ); + + EvtLeptonVandACurrents( l11, a11, lp1, lm1 ); + EvtLeptonVandACurrents( l21, a21, lp0, lm1 ); + EvtLeptonVandACurrents( l12, a12, lp1, lm0 ); + EvtLeptonVandACurrents( l22, a22, lp0, lm0 ); + + s11 = EvtLeptonSCurrent( lp1, lm1 ); + p11 = EvtLeptonPCurrent( lp1, lm1 ); + s21 = EvtLeptonSCurrent( lp0, lm1 ); + p21 = EvtLeptonPCurrent( lp0, lm1 ); + s12 = EvtLeptonSCurrent( lp1, lm0 ); + p12 = EvtLeptonPCurrent( lp1, lm0 ); + s22 = EvtLeptonSCurrent( lp0, lm0 ); + p22 = EvtLeptonPCurrent( lp0, lm0 ); + } else { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Wrong lepton number\n"; + T1.zero(); + T2.zero(); // Set all tensor terms to zero. + } + + for ( int i = 0; i < 3; i++ ) { + EvtVector4C eps = parent->getDaug( 0 )->epsParent( i ).conj(); + + EvtVector4C E1 = T1.cont1( eps ); + EvtVector4C E2 = T2.cont1( eps ); + + EvtComplex epsq = I * ( eps * q ), epsqs = epsq * hs, epsqp = epsq * hp; + + amp.vertex( i, 0, 0, l11 * E1 + a11 * E2 - s11 * epsqs - p11 * epsqp ); + amp.vertex( i, 0, 1, l12 * E1 + a12 * E2 - s12 * epsqs - p12 * epsqp ); + amp.vertex( i, 1, 0, l21 * E1 + a21 * E2 - s21 * epsqs - p21 * epsqp ); + amp.vertex( i, 1, 1, l22 * E1 + a22 * E2 - s22 * epsqs - p22 * epsqp ); + } +} + +void getScalarFF( double q2, double& fp, double& f0, double& fT ) +{ + // The form factors are taken from: + // B -> K and D -> K form factors from fully relativistic lattice QCD + // W. G. Parrott, C. Bouchard, C. T. H. Davies + // https://arxiv.org/abs/2207.12468 + static const double MK = 0.495644, MB = 5.279495108051249, + MBs0 = 5.729495108051249, MBsstar = 5.415766151925566, + logsB = 1.3036556717286227; + static const int N = 3; + static const double a0B[] = { 0.2546162729845652, 0.211016713841977, + 0.02690776720598433, 0.2546162729845652, + -0.7084710010870424, 0.3096901516968882, + 0.2549078414069112, -0.6549384905373116, + 0.36265904141973127 }, + *apB = a0B + N, *aTB = apB + N; + + double pole0 = 1 / ( 1 - q2 / ( MBs0 * MBs0 ) ); + double polep = 1 / ( 1 - q2 / ( MBsstar * MBsstar ) ); + double B = MB + MK, A = sqrt( B * B - q2 ), z = ( A - B ) / ( A + B ), + zN = z * z * z / N, zn = z; + f0 = a0B[0]; + fp = apB[0]; + fT = aTB[0]; + for ( int i = 1; i < N; i++ ) { + f0 += a0B[i] * zn; + fp += apB[i] * zn; + fT += aTB[i] * zn; + double izN = i * zN; + if ( ( i - N ) & 1 ) { + fp += apB[i] * izN; + fT += aTB[i] * izN; + } else { + fp -= apB[i] * izN; + fT -= aTB[i] * izN; + } + zn *= z; + } + f0 *= pole0 * logsB; + fp *= polep * logsB; + fT *= polep * logsB; +} + +void EvtbTosllVectorAmpNP::CalcSAmp( EvtParticle* parent, EvtAmp& amp, + EvtbTosllFF* formFactors ) +{ + // Add the lepton and neutrino 4 momenta to find q2 + EvtId pId = parent->getId(); + EvtId dId = parent->getDaug( 0 )->getId(); + + EvtVector4R q = parent->getDaug( 1 )->getP4() + parent->getDaug( 2 )->getP4(); + double q2 = q.mass2(); + double dmass = parent->getDaug( 0 )->mass(); + + double fp, f0, ft; // form factors + getScalarFF( q2, fp, f0, ft ); + const EvtVector4R& k = parent->getDaug( 0 )->getP4(); + double pmass = parent->mass(); + const EvtVector4R p( pmass, 0.0, 0.0, 0.0 ); + EvtVector4R pk = p + k; + + EvtComplex c7 = -0.304; + EvtComplex c9 = C9( q2, interpol( m_reso, q2 ) ); + EvtComplex c10 = -4.103; + c7 += m_dc7; + c9 += m_dc9; + c10 += m_dc10; + + double mb = 4.8 /*GeV/c^2*/, ms = 0.093 /*GeV/c^2*/; + + int charge1 = EvtPDL::chg3( parent->getDaug( 1 )->getId() ); + + EvtParticle* lepPos = ( charge1 > 0 ) ? parent->getDaug( 1 ) + : parent->getDaug( 2 ); + EvtParticle* lepNeg = ( charge1 < 0 ) ? parent->getDaug( 1 ) + : parent->getDaug( 2 ); + + EvtDiracSpinor lp0( lepPos->spParent( 0 ) ), lp1( lepPos->spParent( 1 ) ); + EvtDiracSpinor lm0( lepNeg->spParent( 0 ) ), lm1( lepNeg->spParent( 1 ) ); + + EvtVector4C l11, l12, l21, l22, a11, a12, a21, a22; + EvtComplex s11, s12, s21, s22, p11, p12, p21, p22; + + if ( pId == IdBm || pId == IdaB0 || pId == IdaBs ) { + EvtLeptonVandACurrents( l11, a11, lp0, lm0 ); + EvtLeptonVandACurrents( l21, a21, lp1, lm0 ); + EvtLeptonVandACurrents( l12, a12, lp0, lm1 ); + EvtLeptonVandACurrents( l22, a22, lp1, lm1 ); + + s11 = EvtLeptonSCurrent( lp0, lm0 ); + p11 = EvtLeptonPCurrent( lp0, lm0 ); + s21 = EvtLeptonSCurrent( lp1, lm0 ); + p21 = EvtLeptonPCurrent( lp1, lm0 ); + s12 = EvtLeptonSCurrent( lp0, lm1 ); + p12 = EvtLeptonPCurrent( lp0, lm1 ); + s22 = EvtLeptonSCurrent( lp1, lm1 ); + p22 = EvtLeptonPCurrent( lp1, lm1 ); + } else if ( pId == IdBp || pId == IdB0 || pId == IdBs ) { + EvtLeptonVandACurrents( l11, a11, lp1, lm1 ); + EvtLeptonVandACurrents( l21, a21, lp0, lm1 ); + EvtLeptonVandACurrents( l12, a12, lp1, lm0 ); + EvtLeptonVandACurrents( l22, a22, lp0, lm0 ); + + s11 = EvtLeptonSCurrent( lp1, lm1 ); + p11 = EvtLeptonPCurrent( lp1, lm1 ); + s21 = EvtLeptonSCurrent( lp0, lm1 ); + p21 = EvtLeptonPCurrent( lp0, lm1 ); + s12 = EvtLeptonSCurrent( lp1, lm0 ); + p12 = EvtLeptonPCurrent( lp1, lm0 ); + s22 = EvtLeptonSCurrent( lp0, lm0 ); + p22 = EvtLeptonPCurrent( lp0, lm0 ); + } else { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Wrong lepton number\n"; + } + double dm2 = pmass * pmass - dmass * dmass, t0 = dm2 / q2, + t1 = 2 * mb * ft / ( pmass + dmass ); + c7 += m_c7p; + c9 += m_c9p; + c10 += m_c10p; + EvtVector4C E1 = ( c9 * fp + c7 * t1 ) * pk + + ( t0 * ( c9 * ( f0 - fp ) - c7 * t1 ) ) * q; + EvtVector4C E2 = ( c10 * fp ) * pk + ( t0 * ( f0 - fp ) ) * q; + double s = dm2 / ( mb - ms ) * f0; + amp.vertex( 0, 0, l11 * E1 + a11 * E2 + ( m_cS * s11 + m_cP * p11 ) * s ); + amp.vertex( 0, 1, l12 * E1 + a12 * E2 + ( m_cS * s12 + m_cP * p12 ) * s ); + amp.vertex( 1, 0, l21 * E1 + a21 * E2 + ( m_cS * s21 + m_cP * p21 ) * s ); + amp.vertex( 1, 1, l22 * E1 + a22 * E2 + ( m_cS * s22 + m_cP * p22 ) * s ); +} diff --git a/test/evtgenlhc_test1.cc b/test/evtgenlhc_test1.cc --- a/test/evtgenlhc_test1.cc +++ b/test/evtgenlhc_test1.cc @@ -1609,14 +1609,14 @@ directName = strcpy( directName, tkname.c_str() ); directName = strcat( directName, "Direct" ); histo2[ik] = new TH1F( directName, directName, 30, 0.0, 3.0 ); - delete directName; + delete[] directName; char* massName; massName = new char[( strlen( tkname.c_str() ) + 4 )]; massName = strcpy( massName, tkname.c_str() ); massName = strcat( massName, "Mass" ); massHisto[ik] = new TH1F( massName, massName, 3000, 0.0, 5.0 ); - delete massName; + delete[] massName; } count = 1; diff --git a/test/exampleFiles/PI0DALITZ.DEC b/test/exampleFiles/PI0DALITZ.DEC --- a/test/exampleFiles/PI0DALITZ.DEC +++ b/test/exampleFiles/PI0DALITZ.DEC @@ -3,5 +3,8 @@ Decay pi0 1.000 e+ e- gamma PI0_DALITZ; Enddecay +Decay eta +1.000 mu+ mu- gamma PI0_DALITZ; +Enddecay # End