Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F19264916
D87.1759132347.diff
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Award Token
Flag For Later
Size
186 KB
Referenced Files
None
Subscribers
None
D87.1759132347.diff
View Options
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 <vector>
+struct EvtPointf {
+ float x, y;
+};
+
+struct EvtLinSample {
+ EvtLinSample() {}
+ EvtLinSample( const char* );
+ EvtLinSample( const std::vector<EvtPointf>& _v ) { init( _v ); }
+ void init( const std::vector<EvtPointf>& _v );
+ std::pair<double, double> operator()( double ) const;
+
+ std::vector<EvtPointf> v;
+ std::vector<double> 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<double> 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<T>& 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 <cassert>
+
// 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 <https://www.gnu.org/licenses/>. *
***********************************************************************/
-#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 <https://www.gnu.org/licenses/>. *
***********************************************************************/
-#ifndef EVTPI0DALITZ_HH
-#define EVTPI0DALITZ_HH
+#ifndef EVTBTOSLLNP_HH
+#define EVTBTOSLLNP_HH
-#include "EvtGenBase/EvtDecayProb.hh"
+#include "EvtGenBase/EvtDecayAmp.hh"
+#include <memory>
+
+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<EvtbTosllAmp> _calcamp;
+ std::unique_ptr<EvtbTosllFF> _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 <https://www.gnu.org/licenses/>. *
***********************************************************************/
-#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 <memory>
+
+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<EvtbTosllAmp> _calcamp;
+ std::unique_ptr<EvtbTosllFF> _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 <https://www.gnu.org/licenses/>. *
+***********************************************************************/
+
+#ifndef EVTBTOSLLVECTORAMPNP_HH
+#define EVTBTOSLLVECTORAMPNP_HH
+
+#include "EvtGenBase/EvtComplex.hh"
+
+#include "EvtGenModels/EvtbTosllAmp.hh"
+
+#include <ccomplex>
+#include <vector>
+
+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<std::pair<double, std::complex<double>>> 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!"<<endl;
- //}
-
- rho.set( i, j, temp );
+ for ( int i = 1; i < n; i++ )
+ a += abs2( _amp[i] );
+ rho.set( 0, 0, EvtComplex( a ) );
+ } else {
+ int n = 1;
+ for ( int i = 0; i < _ndaug; i++ )
+ n *= dstates[i];
+ for ( int i = 0; i < np; i++ ) {
+ for ( int j = 0; j < np; j++ ) {
+ EvtComplex a;
+ for ( int k = 0; k < n; k++ )
+ a += _amp[np * k + i] * conj( _amp[np * k + j] );
+ rho.set( i, j, a );
}
}
- return rho;
}
+ return rho;
}
EvtSpinDensity EvtAmp::getBackwardSpinDensity( EvtSpinDensity* rho_list )
@@ -258,76 +225,109 @@
return ampprime.contract( _dnontrivial[i], ( *this ) );
}
-EvtAmp EvtAmp::contract( int k, const EvtSpinDensity& rho )
+EvtAmp EvtAmp::contract( int istate, const EvtSpinDensity& rho )
{
+ // contract amplitude tensor with spin density tensor
+ // istate -- contraction index
+ // dimension of rho has to be _nstate[istate]
EvtAmp temp;
-
- int i, j;
temp._ndaug = _ndaug;
temp._pstates = _pstates;
temp._nontrivial = _nontrivial;
-
- for ( i = 0; i < _ndaug; i++ ) {
+ for ( int i = 0; i < _ndaug; i++ ) {
temp.dstates[i] = dstates[i];
temp._dnontrivial[i] = _dnontrivial[i];
}
-
- if ( _nontrivial == 0 ) {
- EvtGenReport( EVTGEN_ERROR, "EvtGen" )
- << "Should not be here EvtAmp!" << endl;
- }
-
- for ( i = 0; i < _nontrivial; i++ ) {
+ for ( int i = 0; i < _nontrivial; i++ )
temp._nstate[i] = _nstate[i];
+ int nup = 1, stride = 1;
+ for ( int i = istate + 1; i < _nontrivial; i++ )
+ nup *= _nstate[i];
+ for ( int i = 0; i < istate; i++ )
+ stride *= _nstate[i];
+
+ int ns = _nstate[istate], offstride = stride * ns, namp = offstride * nup;
+#if 1
+ for ( int i = 0; i < stride; i++ ) {
+ for ( int off = i; off < namp; off += offstride ) {
+ for ( int k = 0; k < ns; k++ ) {
+ EvtComplex c;
+ for ( int j = 0; j < ns; j++ )
+ c += rho( j, k ) * _amp[off + j * stride];
+ temp._amp[off + k * stride] = c;
+ }
+ }
}
-
+ // temp.dump();
+#endif
+#if 0
+ int index[10] = {0};
+ for (int i = 0; i < namp; i++ ) {
EvtComplex c;
+ int tempint = index[istate];
+ for (int j = 0; j < _nstate[istate]; j++ ) {
+ // std::cout<<i<<" "<<j<<" "<<tempint<<std::endl;
+ index[istate] = j;
+ c += rho.get( j, tempint ) * getAmp( index );
+ }
+ index[istate] = tempint;
+ temp.setAmp( index, c );
+ int indflag = 0;
+ for (int j = 0; j < _nontrivial; j++ ) {
+ if ( indflag == 0 ) {
+ if ( index[j] == ( _nstate[j] - 1 ) ) {
+ index[j] = 0;
+ } else {
+ indflag = 1;
+ index[j] += 1;
+ }
+ }
+ }
+ }
+ // temp.dump();
+#endif
+ return temp;
+}
- int index[10];
- for ( i = 0; i < 10; i++ ) {
- index[i] = 0;
- }
-
- int allloop = 1;
- int indflag, ii;
- for ( i = 0; i < _nontrivial; i++ ) {
- allloop *= _nstate[i];
- }
-
- for ( i = 0; i < allloop; i++ ) {
- c = EvtComplex( 0.0 );
- int tempint = index[k];
- for ( j = 0; j < _nstate[k]; j++ ) {
- index[k] = j;
- c += rho.get( j, tempint ) * getAmp( index );
- }
- index[k] = tempint;
-
- temp.setAmp( index, c );
-
- indflag = 0;
- for ( ii = 0; ii < _nontrivial; ii++ ) {
- if ( indflag == 0 ) {
- if ( index[ii] == ( _nstate[ii] - 1 ) ) {
- index[ii] = 0;
- } else {
- indflag = 1;
- index[ii] += 1;
+EvtSpinDensity EvtAmp::contract( int istate, const EvtAmp& amp2 )
+{
+#if 1
+ int nup = 1, stride = 1;
+ for ( int i = istate + 1; i < _nontrivial; i++ )
+ nup *= _nstate[i];
+ for ( int i = 0; i < istate; i++ )
+ stride *= _nstate[i];
+
+ int ns = _nstate[istate], offstride = stride * ns, namp = offstride * nup;
+ EvtSpinDensity mrho( ns );
+ for ( int i = 0; i < stride; i++ ) {
+ for ( int off = i; off < namp; off += offstride ) {
+ for ( int k = 0; k < ns; k++ ) {
+ EvtComplex ak = conj( amp2._amp[off + k * stride] );
+ for ( int j = k; j < ns; j++ ) {
+ EvtComplex aj = _amp[off + j * stride];
+ mrho( j, k ) += aj * ak;
+ // if(k==j)std::cout<<" "<<i<<" "<<off<<" "<<k<<" "<<ak<<" "<<aj<<" "<<aj*ak<<std::endl;
}
}
}
}
- return temp;
-}
+ for ( int k = 0; k < ns; k++ ) {
+ for ( int j = 0; j < k; j++ ) {
+ mrho( j, k ) = conj( mrho( k, j ) );
+ }
+ }
+ // std::cout<<mrho;
+ return mrho;
+#endif
-EvtSpinDensity EvtAmp::contract( int k, const EvtAmp& amp2 )
-{
+#if 0
+ int k = istate;
int i, j, l;
EvtComplex temp;
- EvtSpinDensity rho;
-
- rho.setDim( _nstate[k] );
+ EvtSpinDensity rho( _nstate[k] );
+ // rho.setDim( _nstate[k] );
int allloop = 1;
int indflag, ii;
@@ -377,8 +377,9 @@
rho.set( i, j, temp );
}
}
-
+ // std::cout<<rho-mrho;
return rho;
+#endif
}
EvtAmp EvtAmp::contract( int, const EvtAmp&, const EvtAmp& )
@@ -450,7 +451,7 @@
EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
<< "-----------------------------------" << endl;
}
-
+/*
void EvtAmp::vertex( const EvtComplex& c )
{
int list[1];
@@ -486,6 +487,7 @@
{
setAmp( i1, c );
}
+*/
EvtAmp& EvtAmp::operator=( const EvtAmp& amp )
{
diff --git a/src/EvtGenBase/EvtCPUtil.cpp b/src/EvtGenBase/EvtCPUtil.cpp
--- a/src/EvtGenBase/EvtCPUtil.cpp
+++ b/src/EvtGenBase/EvtCPUtil.cpp
@@ -198,7 +198,7 @@
<< "p=" << EvtPDL::name( p->getId() )
<< " parent=" << EvtPDL::name( parent->getId() ) << endl;
}
- assert( parent == 0 );
+ assert( parent == nullptr );
p->setLifetime();
t = p->getLifetime();
bool needToChargeConj = false;
diff --git a/src/EvtGenBase/EvtComplex.cpp b/src/EvtGenBase/EvtComplex.cpp
--- a/src/EvtGenBase/EvtComplex.cpp
+++ b/src/EvtGenBase/EvtComplex.cpp
@@ -31,27 +31,3 @@
s << "(" << c._rpart << "," << c._ipart << ")";
return s;
}
-
-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;
-}
-
-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;
-}
diff --git a/src/EvtGenBase/EvtDiracParticle.cpp b/src/EvtGenBase/EvtDiracParticle.cpp
--- a/src/EvtGenBase/EvtDiracParticle.cpp
+++ b/src/EvtGenBase/EvtDiracParticle.cpp
@@ -46,29 +46,30 @@
::abort();
}
+ 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 t = sqrt( em ), q = 1 / t;
+ double ux = px * q, uy = py * q, uz = pz * q;
+ double lm = sqrt( 2.0 * m );
+ EvtComplex o( 0.0, 0.0 ), l( lm, 0.0 );
+ EvtComplex uxy( ux, uy ), cxy = conj( uxy );
if ( EvtPDL::getStdHep( part_n ) > 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 <<std::endl;
+
setLifetime();
}
diff --git a/src/EvtGenBase/EvtDiracSpinor.cpp b/src/EvtGenBase/EvtDiracSpinor.cpp
--- a/src/EvtGenBase/EvtDiracSpinor.cpp
+++ b/src/EvtGenBase/EvtDiracSpinor.cpp
@@ -31,26 +31,6 @@
#include <math.h>
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 <<std::endl; // check that everything is correct
+ return res;
}
EvtVector4C EvtLeptonVCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp )
{
- EvtVector4C temp;
+ // EvtVector4C temp;
// no conjugate here; done in the multiplication
// yes this is stupid and fooled me to for a long time (ryd)
- temp.set( 0, d * ( EvtGammaMatrix::v0() * dp ) );
- temp.set( 1, d * ( EvtGammaMatrix::v1() * dp ) );
- temp.set( 2, d * ( EvtGammaMatrix::v2() * dp ) );
- temp.set( 3, d * ( EvtGammaMatrix::v3() * dp ) );
-
- return temp;
+ // temp.set( 0, d * ( EvtGammaMatrix::v0() * dp ) );
+ // temp.set( 1, d * ( EvtGammaMatrix::v1() * dp ) );
+ // temp.set( 2, d * ( EvtGammaMatrix::v2() * dp ) );
+ // temp.set( 3, d * ( EvtGammaMatrix::v3() * dp ) );
+
+ // 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;
+ double t0r = a0r * b3r + a0i * b3i, t0i = a0r * b3i - a0i * b3r;
+ double t1r = a1r * b2r + a1i * b2i, t1i = a1r * b2i - a1i * b2r;
+ double t2r = a2r * b1r + a2i * b1i, t2i = a2r * b1i - a2i * b1r;
+ double t3r = a3r * b0r + a3i * b0i, t3i = a3r * b0i - a3i * b0r;
+ double q0r = a0r * b2r + a0i * b2i, q0i = a0r * b2i - a0i * b2r;
+ double q1r = a1r * b3r + a1i * b3i, q1i = a1r * b3i - a1i * b3r;
+ double q2r = a2r * b0r + a2i * b0i, q2i = a2r * b0i - a2i * b0r;
+ double q3r = a3r * b1r + a3i * b1i, q3i = a3r * b1i - a3i * b1r;
+
+ EvtComplex c0( w0r + w1r + ( w2r + w3r ), w0i + w1i + ( w2i + w3i ) );
+ EvtComplex c1( t0r + t1r + ( t2r + t3r ), t0i + t1i + ( t2i + t3i ) );
+ EvtComplex c2( t0i + t2i - ( t1i + t3i ), t1r + t3r - ( t0r + t2r ) );
+ EvtComplex c3( q0r + q2r - ( q1r + q3r ), q0i + q2i - ( q1i + q3i ) );
+
+ return EvtVector4C( c0, c1, c2, c3 );
+
+ // EvtComplex t0 = temp.get(0), t1 = temp.get(1), t2 = temp.get(2), t3 = temp.get(3);
+ // std::cout<<c0-t0<<" "<<c1-t1<<" "<<c2-t2<<" "<<c3-t3<<std::endl;
}
EvtVector4C EvtLeptonACurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp )
{
- EvtVector4C temp;
+ // EvtVector4C temp;
- EvtGammaMatrix mat;
+ // EvtGammaMatrix mat;
// no conjugate here; done in the multiplication
// yes this is stupid and fooled me to for a long time (ryd)
- mat = EvtGammaMatrix::v0() - EvtGammaMatrix::va0();
- temp.set( 0, d * ( mat * dp ) );
+ // mat = EvtGammaMatrix::v0() - EvtGammaMatrix::va0();
+ // temp.set( 0, d * ( mat * dp ) );
- mat = EvtGammaMatrix::v1() - EvtGammaMatrix::va1();
- temp.set( 1, d * ( mat * dp ) );
+ // mat = EvtGammaMatrix::v1() - EvtGammaMatrix::va1();
+ // temp.set( 1, d * ( mat * dp ) );
- mat = EvtGammaMatrix::v2() - EvtGammaMatrix::va2();
- temp.set( 2, d * ( mat * dp ) );
+ // mat = EvtGammaMatrix::v2() - EvtGammaMatrix::va2();
+ // temp.set( 2, d * ( mat * dp ) );
- mat = EvtGammaMatrix::v3() - EvtGammaMatrix::va3();
- temp.set( 3, d * ( mat * dp ) );
+ // mat = EvtGammaMatrix::v3() - EvtGammaMatrix::va3();
+ // temp.set( 3, d * ( mat * dp ) );
- 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;
+
+ double q0r = a0r * b2r + a0i * b2i, q0i = a0r * b2i - a0i * b2r;
+ double q1r = a1r * b3r + a1i * b3i, q1i = a1r * b3i - a1i * b3r;
+ double q2r = a2r * b0r + a2i * b0i, q2i = a2r * b0i - a2i * b0r;
+ double q3r = a3r * b1r + a3i * b1i, q3i = a3r * b1i - a3i * b1r;
+
+ double t0r = a0r * b1r + a0i * b1i, t0i = a0r * b1i - a0i * b1r;
+ double t1r = a1r * b0r + a1i * b0i, t1i = a1r * b0i - a1i * b0r;
+ double t2r = a2r * b3r + a2i * b3i, t2i = a2r * b3i - a2i * b3r;
+ double t3r = a3r * b2r + a3i * b2i, t3i = a3r * b2i - a3i * b2r;
+
+ EvtComplex c0( q0r + q2r + ( q1r + q3r ), q0i + q2i + ( q1i + q3i ) );
+ EvtComplex c1( t0r + t1r + ( t2r + t3r ), t0i + t1i + ( t2i + t3i ) );
+ EvtComplex c2( t0i + t2i - ( t1i + t3i ), t1r + t3r - ( t0r + t2r ) );
+ EvtComplex c3( w0r + w2r - ( w1r + w3r ), w0i + w2i - ( w1i + w3i ) );
+
+ return EvtVector4C( c0, c1, c2, c3 );
+
+ // std::cout<<EvtVector4C(c0,c1,c2,c3)-temp<<std::endl;
+ // return temp;
+}
+
+void EvtLeptonVandACurrents( EvtVector4C& v, EvtVector4C& a,
+ const EvtDiracSpinor& x, const EvtDiracSpinor& y )
+{
+ EvtComplex w0 = x[0] & y[0], w1 = x[1] & y[1], w2 = x[2] & y[2],
+ w3 = x[3] & y[3];
+ EvtComplex w02 = w0 + w2, w13 = w1 + w3;
+ EvtComplex W1 = w02 + w13, W2 = w02 - w13;
+ EvtComplex q0 = x[0] & y[2], q1 = x[1] & y[3], q2 = x[2] & y[0],
+ q3 = x[3] & y[1];
+ EvtComplex q02 = q0 + q2, q13 = q1 + q3;
+ EvtComplex Q1 = q02 + q13, Q2 = q02 - q13;
+ EvtComplex e0 = x[0] & y[3], e1 = x[1] & y[2], e2 = x[2] & y[1],
+ e3 = x[3] & y[0];
+ EvtComplex e20 = e0 + e2, e13 = e1 + e3;
+ EvtComplex E1 = e13 + e20, E2 = e13 - e20;
+ v.set( W1, E1, EvtComplex( -imag( E2 ), real( E2 ) ), Q2 );
+ EvtComplex t0 = x[0] & y[1], t1 = x[1] & y[0], t2 = x[2] & y[3],
+ t3 = x[3] & y[2];
+ EvtComplex t20 = t0 + t2, t13 = t1 + t3;
+ EvtComplex T1 = t13 + t20, T2 = t13 - t20;
+ a.set( Q1, T1, EvtComplex( -imag( T2 ), real( T2 ) ), W2 );
}
EvtComplex EvtLeptonSCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp )
{
- EvtComplex temp;
+ // EvtComplex temp;
// no conjugate here; done in the multiplication
// yes this is stupid and fooled me to for a long time (ryd)
- temp = d * ( EvtGammaMatrix::g0() * dp );
+ // temp = d * ( EvtGammaMatrix::g0() * dp );
+
+ // 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] );
- return temp;
+ 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] );
+
+ EvtComplex res( ( a0r * b0r + a0i * b0i ) + ( a1r * b1r + a1i * b1i ) -
+ ( a2r * b2r + a2i * b2i ) + ( a3r * b3r + a3i * b3i ),
+ ( a0r * b0i - a0i * b0r ) + ( a1r * b1i - a1i * b1r ) -
+ ( a2r * b2i - a2i * b2r ) + ( a3r * b3i - a3i * b3r ) );
+ // std::cout<<res-temp<<std::endl;
+ return res;
}
EvtComplex EvtLeptonPCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp )
{
- EvtComplex temp;
+ // EvtComplex temp;
// no conjugate here; done in the multiplication
// yes this is stupid and fooled me to for a long time (ryd)
- static EvtGammaMatrix m = EvtGammaMatrix::g0() * EvtGammaMatrix::g5();
- temp = d * ( m * dp );
-
- return temp;
+ // static EvtGammaMatrix m = EvtGammaMatrix::g0() * EvtGammaMatrix::g5();
+ // temp = d * ( m * dp );
+
+ // 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] );
+
+ EvtComplex res( ( a0r * b2r + a0i * b2i ) + ( a1r * b3r + a1i * b3i ) -
+ ( a2r * b0r + a2i * b0i ) - ( a3r * b1r + a3i * b1i ),
+ ( a0r * b2i - a0i * b2r ) + ( a1r * b3i - a1i * b3r ) -
+ ( a2r * b0i - a2i * b0r ) - ( a3r * b1i - a3i * b1r ) );
+ // std::cout<<res-temp<<std::endl;
+ return res;
}
EvtTensor4C EvtLeptonTCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp )
{
- EvtTensor4C temp;
- temp.zero();
- EvtComplex i2( 0, 0.5 );
-
- static EvtGammaMatrix mat01 =
- EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g1() -
- EvtGammaMatrix::g1() * EvtGammaMatrix::g0() );
- static EvtGammaMatrix mat02 =
- EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g2() -
- EvtGammaMatrix::g2() * EvtGammaMatrix::g0() );
- static EvtGammaMatrix mat03 =
- EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g3() -
- EvtGammaMatrix::g3() * EvtGammaMatrix::g0() );
- static EvtGammaMatrix mat12 =
- EvtGammaMatrix::g0() * ( EvtGammaMatrix::g1() * EvtGammaMatrix::g2() -
- EvtGammaMatrix::g2() * EvtGammaMatrix::g1() );
- static EvtGammaMatrix mat13 =
- EvtGammaMatrix::g0() * ( EvtGammaMatrix::g1() * EvtGammaMatrix::g3() -
- EvtGammaMatrix::g3() * EvtGammaMatrix::g1() );
- static EvtGammaMatrix mat23 =
- EvtGammaMatrix::g0() * ( EvtGammaMatrix::g2() * EvtGammaMatrix::g3() -
- EvtGammaMatrix::g3() * EvtGammaMatrix::g2() );
-
- temp.set( 0, 1, i2 * ( d * ( mat01 * dp ) ) );
- temp.set( 1, 0, -temp.get( 0, 1 ) );
-
- temp.set( 0, 2, i2 * ( d * ( mat02 * dp ) ) );
- temp.set( 2, 0, -temp.get( 0, 2 ) );
-
- temp.set( 0, 3, i2 * ( d * ( mat03 * dp ) ) );
- temp.set( 3, 0, -temp.get( 0, 3 ) );
-
- temp.set( 1, 2, i2 * ( d * ( mat12 * dp ) ) );
- temp.set( 2, 1, -temp.get( 1, 2 ) );
-
- temp.set( 1, 3, i2 * ( d * ( mat13 * dp ) ) );
- temp.set( 3, 1, -temp.get( 1, 3 ) );
-
- temp.set( 2, 3, i2 * ( d * ( mat23 * dp ) ) );
- temp.set( 3, 2, -temp.get( 2, 3 ) );
-
- return temp;
+ // EvtTensor4C temp;
+ // temp.zero();
+ // EvtComplex i2( 0, 0.5 );
+
+ // static EvtGammaMatrix mat01 =
+ // EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g1() -
+ // EvtGammaMatrix::g1() * EvtGammaMatrix::g0() );
+ // static EvtGammaMatrix mat02 =
+ // EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g2() -
+ // EvtGammaMatrix::g2() * EvtGammaMatrix::g0() );
+ // static EvtGammaMatrix mat03 =
+ // EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g3() -
+ // EvtGammaMatrix::g3() * EvtGammaMatrix::g0() );
+ // static EvtGammaMatrix mat12 =
+ // EvtGammaMatrix::g0() * ( EvtGammaMatrix::g1() * EvtGammaMatrix::g2() -
+ // EvtGammaMatrix::g2() * EvtGammaMatrix::g1() );
+ // static EvtGammaMatrix mat13 =
+ // EvtGammaMatrix::g0() * ( EvtGammaMatrix::g1() * EvtGammaMatrix::g3() -
+ // EvtGammaMatrix::g3() * EvtGammaMatrix::g1() );
+ // static EvtGammaMatrix mat23 =
+ // EvtGammaMatrix::g0() * ( EvtGammaMatrix::g2() * EvtGammaMatrix::g3() -
+ // EvtGammaMatrix::g3() * EvtGammaMatrix::g2() );
+
+ // temp.set( 0, 1, i2 * ( d * ( mat01 * dp ) ) );
+ // temp.set( 1, 0, -temp.get( 0, 1 ) );
+
+ // temp.set( 0, 2, i2 * ( d * ( mat02 * dp ) ) );
+ // temp.set( 2, 0, -temp.get( 0, 2 ) );
+
+ // temp.set( 0, 3, i2 * ( d * ( mat03 * dp ) ) );
+ // temp.set( 3, 0, -temp.get( 0, 3 ) );
+
+ // temp.set( 1, 2, i2 * ( d * ( mat12 * dp ) ) );
+ // temp.set( 2, 1, -temp.get( 1, 2 ) );
+
+ // temp.set( 1, 3, i2 * ( d * ( mat13 * dp ) ) );
+ // temp.set( 3, 1, -temp.get( 1, 3 ) );
+
+ // temp.set( 2, 3, i2 * ( d * ( mat23 * dp ) ) );
+ // temp.set( 3, 2, -temp.get( 2, 3 ) );
+
+ // 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 t0r = a0r * b3r + a0i * b3i, t0i = a0r * b3i - a0i * b3r;
+ double t1r = a1r * b2r + a1i * b2i, t1i = a1r * b2i - a1i * b2r;
+ double t2r = a2r * b1r + a2i * b1i, t2i = a2r * b1i - a2i * b1r;
+ double t3r = a3r * b0r + a3i * b0i, t3i = a3r * b0i - a3i * b0r;
+
+ double t01r = t2i + t3i - ( t0i + t1i );
+ double t01i = t0r + t1r - ( t2r + t3r );
+
+ double t02r = t0r + t3r - ( t1r + t2r );
+ double t02i = t0i + t3i - ( t1i + t2i );
+
+ double t03r = ( a1r * b3i - a1i * b3r ) + ( a2r * b0i - a2i * b0r ) -
+ ( a0r * b2i - a0i * b2r ) - ( a3r * b1i - a3i * b1r );
+ double t03i = ( a0r * b2r + a0i * b2i ) - ( a1r * b3r + a1i * b3i ) -
+ ( a2r * b0r + a2i * b0i ) + ( a3r * b1r + a3i * b1i );
+
+ double t12r = ( a0r * b0r + a0i * b0i ) - ( a1r * b1r + a1i * b1i ) -
+ ( a2r * b2r + a2i * b2i ) + ( a3r * b3r + a3i * b3i );
+ double t12i = ( a0r * b0i - a0i * b0r ) - ( a1r * b1i - a1i * b1r ) -
+ ( a2r * b2i - a2i * b2r ) + ( a3r * b3i - a3i * b3r );
+
+ double q0r = a0r * b1r + a0i * b1i, q0i = a0r * b1i - a0i * b1r;
+ double q1r = a1r * b0r + a1i * b0i, q1i = a1r * b0i - a1i * b0r;
+ double q2r = a2r * b3r + a2i * b3i, q2i = a2r * b3i - a2i * b3r;
+ double q3r = a3r * b2r + a3i * b2i, q3i = a3r * b2i - a3i * b2r;
+
+ double t13r = q1i + q2i - ( q0i + q3i );
+ double t13i = q0r + q3r - ( q1r + q2r );
+
+ double t23r = q0r + q1r - ( q2r + q3r );
+ double t23i = q0i + q1i - ( q2i + q3i );
+
+ EvtTensor4C res;
+ res.set( 0, 1, EvtComplex( t01r, t01i ) );
+ res.set( 1, 0, EvtComplex( -t01r, -t01i ) );
+
+ res.set( 0, 2, EvtComplex( t02r, t02i ) );
+ res.set( 2, 0, EvtComplex( -t02r, -t02i ) );
+
+ res.set( 0, 3, EvtComplex( t03r, t03i ) );
+ res.set( 3, 0, EvtComplex( -t03r, -t03i ) );
+
+ res.set( 1, 2, EvtComplex( t12r, t12i ) );
+ res.set( 2, 1, EvtComplex( -t12r, -t12i ) );
+
+ res.set( 1, 3, EvtComplex( t13r, t13i ) );
+ res.set( 3, 1, EvtComplex( -t13r, -t13i ) );
+
+ res.set( 2, 3, EvtComplex( t23r, t23i ) );
+ res.set( 3, 2, EvtComplex( -t23r, -t23i ) );
+
+ // std::cout<<res - temp<<std::endl;
+
+ return res;
}
EvtDiracSpinor operator*( const EvtComplex& c, const EvtDiracSpinor& d )
{
- EvtDiracSpinor result;
- result.spinor[0] = c * d.spinor[0];
- result.spinor[1] = c * d.spinor[1];
- result.spinor[2] = c * d.spinor[2];
- result.spinor[3] = c * d.spinor[3];
- return result;
+ // EvtDiracSpinor result;
+ // result.spinor[0] = c * d.spinor[0];
+ // result.spinor[1] = c * d.spinor[1];
+ // result.spinor[2] = c * d.spinor[2];
+ // result.spinor[3] = c * d.spinor[3];
+ // return result;
+ return EvtDiracSpinor( c * d[0], c * d[1], c * d[2], c * d[3] );
}
EvtDiracSpinor EvtDiracSpinor::adjoint() const
{
- 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;
+ // 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 <algorithm>
#include <cmath>
-#include <iostream>
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<double>());
+ 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:"<<p4[0]<<p4[1]<<p4[2]<<endl;
-
- double alpha = EvtRandom::Flat( EvtConst::twoPi );
- double beta = acos( EvtRandom::Flat( -1.0, 1.0 ) );
- double gamma = EvtRandom::Flat( EvtConst::twoPi );
-
- p4[0].applyRotateEuler( alpha, beta, gamma );
- p4[1].applyRotateEuler( alpha, beta, gamma );
- p4[2].applyRotateEuler( alpha, beta, gamma );
+ 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 * m12sq );
}
diff --git a/src/EvtGenBase/EvtHepMCEvent.cpp b/src/EvtGenBase/EvtHepMCEvent.cpp
--- a/src/EvtGenBase/EvtHepMCEvent.cpp
+++ b/src/EvtGenBase/EvtHepMCEvent.cpp
@@ -165,7 +165,7 @@
} // Have daughter products
- } // hepMCDaughter != 0
+ } // hepMCDaughter != nullptr
} // Loop over daughters
}
diff --git a/src/EvtGenBase/EvtIdSet.cpp b/src/EvtGenBase/EvtIdSet.cpp
--- a/src/EvtGenBase/EvtIdSet.cpp
+++ b/src/EvtGenBase/EvtIdSet.cpp
@@ -445,9 +445,7 @@
{
int combLen = _numInList + set1.sizeOfSet();
int uniqueLen = 0;
- EvtId* combSet;
-
- combSet = new EvtId[combLen];
+ EvtId* combSet = new EvtId[combLen];
int i;
for ( i = 0; i < combLen; i++ ) {
@@ -468,17 +466,17 @@
combSet[uniqueLen] = _list[i];
uniqueLen += 1;
}
+ }
- delete _list;
- _list = new EvtId[uniqueLen];
-
- _numInList = uniqueLen;
- for ( i = 0; i < _numInList; i++ ) {
- _list[i] = combSet[i];
- }
+ delete[] _list;
+ _list = new EvtId[uniqueLen];
- delete combSet;
+ _numInList = uniqueLen;
+ for ( i = 0; i < _numInList; i++ ) {
+ _list[i] = combSet[i];
}
+
+ delete[] combSet;
}
int EvtIdSet::sizeOfSet() const
diff --git a/src/EvtGenBase/EvtLinSample.cpp b/src/EvtGenBase/EvtLinSample.cpp
new file mode 100644
--- /dev/null
+++ b/src/EvtGenBase/EvtLinSample.cpp
@@ -0,0 +1,58 @@
+#include "EvtGenBase/EvtLinSample.hh"
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+
+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<EvtPointf>& _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<double, double> 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 <iostream>
-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<<t0<<endl;
+ std::cout<<t1<<endl;
+ std::cout<<t2<<endl;
+
+ std::cout<<t0*t0<<endl;
+ std::cout<<t1*t1<<endl;
+ std::cout<<t2*t2<<endl;
+ std::cout<<t0*t1<<endl;
+ std::cout<<t0*t2<<endl;
+ std::cout<<t1*t2<<endl;
+
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ std::cout<<t0.get(i)*t0.get(j) + t1.get(i)*t1.get(j) + t2.get(i)*t2.get(j)<<" ";
+ }
+ std::cout<<endl;
+ }
+
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ std::cout<<-gmn(i,j)+k.get(i)*k.get(j)/k.mass2()<<" ";
+ }
+ std::cout<<endl;
+ }
+ */
+
+ EvtVector4C et0 = tds.cont1( t0.conj() );
+ EvtVector4C et1 = tds.cont1( t1.conj() );
+ EvtVector4C et2 = tds.cont1( t2.conj() );
amp.vertex( 0, 0, l1.cont( et0 ) );
amp.vertex( 0, 1, l2.cont( et0 ) );
diff --git a/src/EvtGenBase/EvtSimpleRandomEngine.cpp b/src/EvtGenBase/EvtSimpleRandomEngine.cpp
--- a/src/EvtGenBase/EvtSimpleRandomEngine.cpp
+++ b/src/EvtGenBase/EvtSimpleRandomEngine.cpp
@@ -33,3 +33,11 @@
return ( temp + 1.0 ) / 32769.0;
}
+
+unsigned long EvtSimpleRandomEngine::urandom()
+{
+ _next = _next * 1103515245 + 123345;
+ unsigned temp = (unsigned)( _next / 65536 ) % 32768;
+
+ return (long)temp << 49;
+}
diff --git a/src/EvtGenBase/EvtSpinDensity.cpp b/src/EvtGenBase/EvtSpinDensity.cpp
--- a/src/EvtGenBase/EvtSpinDensity.cpp
+++ b/src/EvtGenBase/EvtSpinDensity.cpp
@@ -24,7 +24,7 @@
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh"
-#include <assert.h>
+#include <cstring>
#include <iostream>
#include <math.h>
#include <stdlib.h>
@@ -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<EvtComplex*>( t );
+ const EvtComplex* u1 = reinterpret_cast<const EvtComplex*>( 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<<res<<std::endl;
+ return res;
+ */
+ double t01 = _b.get( 2 ) * _c.get( 3 ) - _b.get( 3 ) * _c.get( 2 );
+ double t02 = _b.get( 3 ) * _c.get( 1 ) - _b.get( 1 ) * _c.get( 3 );
+ double t03 = _b.get( 1 ) * _c.get( 2 ) - _b.get( 2 ) * _c.get( 1 );
+ double t12 = _b.get( 0 ) * _c.get( 3 ) - _b.get( 3 ) * _c.get( 0 );
+ double t13 = _b.get( 2 ) * _c.get( 0 ) - _b.get( 0 ) * _c.get( 2 );
+ double t23 = _b.get( 0 ) * _c.get( 1 ) - _b.get( 1 ) * _c.get( 0 );
+ double u[4] = { 0 };
+ u[0] = t01 * a[1] + t02 * a[2] + t03 * a[3];
+ u[1] = t12 * a[2] + t13 * a[3] - t01 * a[0];
+ u[2] = t23 * a[3] - ( t02 * a[0] + t12 * a[1] );
+ u[3] = -( t03 * a[0] + t13 * a[1] + t23 * a[2] );
+ // std::cout<<u[0]<<" "<<u[1]<<" "<<u[2]<<" "<<u[3]<<"\n";
+ return EvtVector4R( u[0], u[1], u[2], u[3] );
+}
+
+EvtVector4C EvtGenFunctions::asymProd( const EvtVector4C& _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
+
+ EvtComplex a[4] = { _a.get( 0 ), _a.get( 1 ), _a.get( 2 ), _a.get( 3 ) };
+ double t01 = _b.get( 2 ) * _c.get( 3 ) - _b.get( 3 ) * _c.get( 2 );
+ double t02 = _b.get( 3 ) * _c.get( 1 ) - _b.get( 1 ) * _c.get( 3 );
+ double t03 = _b.get( 1 ) * _c.get( 2 ) - _b.get( 2 ) * _c.get( 1 );
+ double t12 = _b.get( 0 ) * _c.get( 3 ) - _b.get( 3 ) * _c.get( 0 );
+ double t13 = _b.get( 2 ) * _c.get( 0 ) - _b.get( 0 ) * _c.get( 2 );
+ double t23 = _b.get( 0 ) * _c.get( 1 ) - _b.get( 1 ) * _c.get( 0 );
+ return EvtVector4C( t01 * a[1] + t02 * a[2] + t03 * a[3],
+ t12 * a[2] + t13 * a[3] - t01 * a[0],
+ t23 * a[3] - ( t02 * a[0] + t12 * a[1] ),
+ -( t03 * a[0] + t13 * a[1] + t23 * a[2] ) );
+}
diff --git a/src/EvtGenBase/EvtVector4C.cpp b/src/EvtGenBase/EvtVector4C.cpp
--- a/src/EvtGenBase/EvtVector4C.cpp
+++ b/src/EvtGenBase/EvtVector4C.cpp
@@ -28,23 +28,6 @@
#include <math.h>
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 "<<t-o<<std::endl;
+ // return o;
+}
diff --git a/src/EvtGenModels/EvtModelReg.cpp b/src/EvtGenModels/EvtModelReg.cpp
--- a/src/EvtGenModels/EvtModelReg.cpp
+++ b/src/EvtGenModels/EvtModelReg.cpp
@@ -158,6 +158,8 @@
#include "EvtGenModels/EvtbTosllBall.hh"
#include "EvtGenModels/EvtbTosllMS.hh"
#include "EvtGenModels/EvtbTosllMSExt.hh"
+#include "EvtGenModels/EvtbTosllNP.hh"
+#include "EvtGenModels/EvtbTosllNPR.hh"
#include "EvtGenModels/Evtbs2llGammaISRFSR.hh"
#include "EvtGenModels/Evtbs2llGammaMNT.hh"
#include "EvtGenModels/EvtbsToLLLL.hh"
@@ -177,7 +179,6 @@
EvtModelReg::EvtModelReg( const std::list<EvtDecayBase*>* extraModels )
{
EvtModel& modelist = EvtModel::instance();
-
if ( extraModels ) {
for ( std::list<EvtDecayBase*>::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") <<" "<<q2<<" "<<prob<<" "<<ew<<std::endl;
- // EvtGenReport(EVTGEN_INFO,"EvtGen") << "prob is "<<prob<<endl;
setProb( prob );
return;
diff --git a/src/EvtGenModels/EvtVSS.cpp b/src/EvtGenModels/EvtVSS.cpp
--- a/src/EvtGenModels/EvtVSS.cpp
+++ b/src/EvtGenModels/EvtVSS.cpp
@@ -67,9 +67,10 @@
EvtVector4R pDaug = p->getDaug( 0 )->getP4();
double norm = 1.0 / pDaug.d3mag();
-
+ pDaug *= norm;
+ // std::cout<<"VSS "<<norm<<std::endl;
for ( int i = 0; i < 3; i++ )
- vertex( i, norm * pDaug * ( p->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 <fstream>
+
+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<double> 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 = "<<m0<<" q2 = "<<q2<<" prob = "<<prob<<std::endl;
if ( i == 0 )
- q2 = 4 * ( mass[1] * mass[1] );
-
- erho = ( m * m + mass[0] * mass[0] - q2 ) / ( 2.0 * m );
+ q2 = q2min; // pole
+ double Erho = ( m * m + m0 * m0 - q2 ) / ( 2.0 * m ),
+ Prho = sqrt( ( Erho - m0 ) * ( Erho + m0 ) );
- prho = sqrt( erho * erho - mass[0] * mass[0] );
-
- p4meson.set( erho, 0.0, 0.0, -1.0 * prho );
- p4w.set( m - erho, 0.0, 0.0, prho );
+ p4meson.set( Erho, 0.0, 0.0, -Prho );
+ p4w.set( m - Erho, 0.0, 0.0, Prho );
//This is in the W rest frame
- elepton = ( q2 + mass[1] * mass[1] ) / ( 2.0 * sqrt( q2 ) );
- plepton = sqrt( elepton * elepton - mass[1] * mass[1] );
-
- double probctl[3];
+ double Elep = sqrt( q2 ) * 0.5,
+ Plep = sqrt( Elep * Elep - mass[1] * mass[1] );
- for ( j = 0; j < 3; j++ ) {
- costl = 0.99 * ( j - 1.0 );
+ const int nj = 3 + 2 + 4 + 8;
+ double p[nj], cmin = -1, cmax = 1, dc = ( cmax - cmin ) / ( nj - 1 );
+ 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.
- p4lepton1.set( elepton, 0.0,
- plepton * sqrt( 1.0 - costl * costl ),
- plepton * costl );
- p4lepton2.set( elepton, 0.0,
- -1.0 * plepton * sqrt( 1.0 - costl * costl ),
- -1.0 * plepton * costl );
+ //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 );
- EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho );
- p4lepton1 = boostTo( p4lepton1, boost );
- p4lepton2 = boostTo( p4lepton2, boost );
+ p4lepton1 = boostTo( p4lepton1, p4w );
+ p4lepton2 = boostTo( p4lepton2, p4w );
//Now initialize the daughters...
daughter->init( meson, p4meson );
lep1->init( lepton1, p4lepton1 );
lep2->init( lepton2, p4lepton2 );
-
+ // std::cout<<j<<" kw "<<q2<<" "<<getKineWeight(root_part)<<" "<<getKineWeight2(root_part)<<std::endl;
CalcAmp( root_part, amp, FormFactors );
- //Now find the probability at this q2 and cos theta lepton point
+ //Now find the probability at this q2 and cos(theta) lepton point
//and compare to maxfoundprob.
//Do a little magic to get the probability!!
- //cout <<"amp:"<<amp.getSpinDensity()<<endl;
-
- prob = rho.normalizedProb( amp.getSpinDensity() );
-
- //cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl;
-
- probctl[j] = prob;
+ p[j] = rho.normalizedProb( amp.getSpinDensity() );
+ // std::cout<<j<<" "<<p[j]<<std::endl;
}
- //probclt contains prob at ctl=-1,0,1.
- //prob=a+b*ctl+c*ctl^2
-
- double a = probctl[1];
- double b = 0.5 * ( probctl[2] - probctl[0] );
- double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1];
-
- prob = probctl[0];
- if ( probctl[1] > 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:"<<prob<<" "
- // << probctl[0]<<" "
- // << probctl[1]<<" "
- // << probctl[2]<<endl;
-
if ( i == 0 ) {
maxpole = prob;
continue;
@@ -217,27 +223,74 @@
if ( prob > maxfoundprob ) {
maxfoundprob = prob;
}
-
- //cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl;
- }
- if ( EvtPDL::getWidth( meson ) <= 0.0 ) {
- //if the particle is narrow dont bother with changing the mass.
- massiter = 4;
}
}
- root_part->deleteTree();
+ if ( 0 ) {
+ std::vector<float> 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:"<<maxfoundprob<<" "
- // <<maxpole<<" "<<poleSize<<endl;
+ //Now initialize the daughters...
+ daughter->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 <cmath>
+
+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 "<<ft<<endl;
}
+static EvtId IdBs, IdaBs, IdRhop, IdRhom, IdRho0, IdOmega, IdKst0, IdaKst0,
+ IdKstp, IdKstm, IdPhi, IdK1p, IdK1m;
+static bool vfirst = true;
void EvtbTosllBallFF::getVectorFF( EvtId parent, EvtId daught, double t,
double /*mass*/, double& a1, double& a2,
double& a0, double& v, double& t1,
double& t2, double& t3 )
{
+ if ( vfirst ) {
+ vfirst = false;
+ 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*-" );
+ IdPhi = EvtPDL::getId( "phi" );
+ IdK1p = EvtPDL::getId( "K_1+" );
+ IdK1m = EvtPDL::getId( "K_1-" );
+ }
int model = _theFFModel;
double m = EvtPDL::getMeanMass( parent );
@@ -198,14 +227,8 @@
double shat = t / ( m * m );
double shat2 = shat * shat;
- if (
- // B_s decay form factors
- parent == EvtPDL::getId( std::string( "B_s0" ) ) ||
- parent == EvtPDL::getId( std::string( "anti-B_s0" ) ) ) {
- if (
- // B_s --> 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 <https://www.gnu.org/licenses/>. *
+***********************************************************************/
+
+#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<EvtbTosllBSZFF>();
+ _calcamp = std::make_unique<EvtbTosllVectorAmpNP>();
+
+ auto getInteger = [this]( int i ) -> int {
+ double a = getArg( i );
+ if ( a - static_cast<int>( 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<int>( a );
+ };
+ EvtbTosllVectorAmpNP* amp = static_cast<EvtbTosllVectorAmpNP*>(
+ _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 <https://www.gnu.org/licenses/>. *
+***********************************************************************/
+
+#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 <algorithm>
+#include <ccomplex>
+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<double> R( double s ) const
+ {
+ return A / std::complex<double>( s - Mv2, MGv );
+ }
+
+ void addknots( std::vector<double>& 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<double> aR( double s )
+{
+ return ajpsi.R( s ) + apsi2s.R( s ) + apsi3770.R( s ) + apsi4040.R( s ) +
+ apsi4160.R( s ) + apsi4230.R( s );
+}
+
+std::vector<double> getgrid()
+{
+ const double m_c = 1.3, m_s = 0.2, m_e = 0.511e-3, MD0 = 1864.84 / 1000;
+
+ std::vector<ABW_t> v = { ajpsi, apsi2s, apsi3770,
+ apsi4040, apsi4160, apsi4230 };
+ std::vector<double> 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<double>::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<double> v = getgrid();
+ std::vector<EvtPointf> 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<double> 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<EvtbTosllBSZFF>();
+ _calcamp = std::make_unique<EvtbTosllVectorAmpNP>();
+
+ auto getInteger = [this]( int i ) -> int {
+ double a = getArg( i );
+ if ( a - static_cast<int>( 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<int>( a );
+ };
+ EvtbTosllVectorAmpNP* amp = static_cast<EvtbTosllVectorAmpNP*>(
+ _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<double> v = getgrid();
+ std::complex<double> eiphi = std::complex<double>( 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 <https://www.gnu.org/licenses/>. *
+***********************************************************************/
+
+#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 <algorithm>
+#include <ccomplex>
+
+static std::complex<double> h( double q2, double mq )
+{
+ const double mu = 4.8; // mu = m_b = 4.8 GeV
+ if ( mq == 0.0 )
+ return std::complex<double>( 8. / 27 + 4. / 9 * log( mu * mu / q2 ),
+ 4. / 9 * M_PI );
+ double z = 4 * mq * mq / q2;
+ std::complex<double> t;
+ if ( z > 1 ) {
+ t = atan( 1 / sqrt( z - 1 ) );
+ } else {
+ t = std::complex<double>( 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<double> hReso = std::complex<double>( 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<double> 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<double> interpol(
+ const std::vector<std::pair<double, std::complex<double>>>& v, double s )
+{
+ if ( !v.size() )
+ return 0;
+ int j = std::lower_bound( v.begin(), v.end(), s,
+ []( const std::pair<double, std::complex<double>>& 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
File Metadata
Details
Attached
Mime Type
text/plain
Expires
Mon, Sep 29, 8:52 AM (9 h, 17 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6560951
Default Alt Text
D87.1759132347.diff (186 KB)
Attached To
Mode
D87: Summary: New Physics generator B->K*ll and performance changes associated with it.
Attached
Detach File
Event Timeline
Log In to Comment