diff --git a/EvtGenBase/EvtAmplitude.hh b/EvtGenBase/EvtAmplitude.hh
index b0b598c..3ecbfe5 100644
--- a/EvtGenBase/EvtAmplitude.hh
+++ b/EvtGenBase/EvtAmplitude.hh
@@ -1,47 +1,51 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtAmplitude.hh,v 1.2 2009-03-16 16:43:40 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// Complex-valued amplitude
#ifndef EVT_AMPLITUDE_HH
#define EVT_AMPLITUDE_HH
#include "EvtGenBase/EvtComplex.hh"
template <class T>
class EvtAmplitude {
public:
- EvtAmplitude() {}
- EvtAmplitude(const EvtAmplitude&) {}
- virtual ~EvtAmplitude() {}
+ EvtAmplitude() = default;
+ EvtAmplitude(const EvtAmplitude&) = default;
+ EvtAmplitude(EvtAmplitude&&) = default;
+ EvtAmplitude& operator=(const EvtAmplitude&) = default;
+ EvtAmplitude& operator=(EvtAmplitude&&) = default;
+ virtual ~EvtAmplitude() = default;
+
virtual EvtAmplitude<T>* clone() const = 0;
EvtComplex evaluate(const T& p) const
{
EvtComplex ret(0.,0.);
if(p.isValid()) ret = amplitude(p);
return ret;
}
protected:
// Derive in subclasses to define amplitude computation
// for a fully constructed amplitude object.
virtual EvtComplex amplitude(const T&) const = 0;
};
#endif
diff --git a/EvtGenBase/EvtDalitzPlot.hh b/EvtGenBase/EvtDalitzPlot.hh
index 3fdf968..c05adfe 100644
--- a/EvtGenBase/EvtDalitzPlot.hh
+++ b/EvtGenBase/EvtDalitzPlot.hh
@@ -1,121 +1,119 @@
//-----------------------------------------------------------------------
// File and Version Information:
// $Id: EvtDalitzPlot.hh,v 1.2 2009-03-16 16:44:53 robbep Exp $
//
// Environment:
// This software is part of the EvtGen package developed jointly
// for the BaBar and CLEO collaborations. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information:
// Copyright (C) 1998 Caltech, UCSB
//
// Module creator:
// Alexei Dvoretskii, Caltech, 2001-2002.
//-----------------------------------------------------------------------
#ifndef EVT_DALITZ_PLOT_HH
#define EVT_DALITZ_PLOT_HH
#include <assert.h>
#include "EvtGenBase/EvtCyclic3.hh"
#include "EvtGenBase/EvtTwoBodyVertex.hh"
#include "EvtGenBase/EvtDecayMode.hh"
class EvtDalitzPlot {
public:
EvtDalitzPlot();
EvtDalitzPlot(double mA, double mB, double mC, double bigM, double ldel = 0., double rdel = 0.);
EvtDalitzPlot(const EvtDecayMode& mode, double ldel = 0., double rdel = 0.);
- EvtDalitzPlot(const EvtDalitzPlot& other);
- ~EvtDalitzPlot();
bool operator==(const EvtDalitzPlot& other) const;
const EvtDalitzPlot* clone() const;
// Absolute limits for masses squared in the Dalitz plot
// e.g. qAbsMin(0) is the lowest possible value
// for m2 of particles {12}
double qAbsMin(EvtCyclic3::Pair i) const;
double qAbsMax(EvtCyclic3::Pair i) const;
double mAbsMin(EvtCyclic3::Pair i) const;
double mAbsMax(EvtCyclic3::Pair i) const;
// Absolute limits for Zemach coordinate qres and qhel (approximate)
// qHelAbsMin(BC,CA) means absolute minimum for (qCA-qAB)/2.
double qResAbsMin(EvtCyclic3::Pair i) const;
double qResAbsMax(EvtCyclic3::Pair i) const;
double qHelAbsMin(EvtCyclic3::Pair i) const;
double qHelAbsMax(EvtCyclic3::Pair i) const;
inline double qSumMin() const { return sum() + _ldel; }
inline double qSumMax() const { return sum() + _rdel; }
inline bool fuzzy() const { return (_rdel - _ldel != 0.); }
// Find the area of the Dalitz plot by numeric integration. (N bins for variable q(i) are used).
// Very large numbers of N can result in a very long calculation. It should not
// matter which two pairs f variables are used. The integral should eventually
// converge to the same number
double getArea(int N = 1000, EvtCyclic3::Pair i = EvtCyclic3::AB, EvtCyclic3::Pair j = EvtCyclic3::BC) const;
// Limits for masses squared when one mass squared is known
double qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j, double q) const;
double qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j, double q) const;
// Coordinate transformations
double cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const;
double e(EvtCyclic3::Index i, EvtCyclic3::Pair j, double q) const;
double p(EvtCyclic3::Index i, EvtCyclic3::Pair j, double q) const;
double q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const;
// |J| of transformation of qi to cosTh in the rest-frame of j
double jacobian(EvtCyclic3::Pair i, double q) const;
// Given resonance index and mass returns decay
// and birth vertices
EvtTwoBodyVertex vD(EvtCyclic3::Pair iRes, double m0, int L) const;
EvtTwoBodyVertex vB(EvtCyclic3::Pair iRes, double m0, int L) const;
// Accessors
double sum() const;
inline double bigM() const { return _bigM; }
inline double mA() const { return _mA; }
inline double mB() const { return _mB; }
inline double mC() const { return _mC; }
double m(EvtCyclic3::Index i) const;
void print() const;
void sanityCheck() const;
protected:
// Defines two dimensional dalitz plot
double _mA;
double _mB;
double _mC;
double _bigM;
// Defines third dimension, or fuzziness. M^2 + ldel < M^2 < M^2 + rdel
double _ldel;
double _rdel;
};
#endif
diff --git a/EvtGenBase/EvtDalitzPoint.hh b/EvtGenBase/EvtDalitzPoint.hh
index 2a8ac56..2a07b7d 100644
--- a/EvtGenBase/EvtDalitzPoint.hh
+++ b/EvtGenBase/EvtDalitzPoint.hh
@@ -1,75 +1,74 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtDalitzPoint.hh,v 1.2 2009-03-16 16:44:53 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// This class describes the complete kinematics of the Dalitz decay.
// It holds all the six invariant momentum products, three daughter
// particle masses and three invariant masses of pairs of particles.
// This description is completely symmetric with respect to particle
// permutations.
//
// Another way to slice the six coordinate is to make a transformation
// to the mass of the decaying particle. The four masses make up a
// Dalitz plot. The other two are coordinates of a point in the plot.
#ifndef EVT_DALITZ_POINT_HH
#define EVT_DALITZ_POINT_HH
#include "EvtGenBase/EvtCyclic3.hh"
#include "EvtGenBase/EvtDalitzCoord.hh"
#include "EvtGenBase/EvtDalitzPlot.hh"
class EvtDalitzPoint final {
public:
EvtDalitzPoint();
EvtDalitzPoint(double mA, double mB, double mC,
double qAB, double qBC, double qCA);
EvtDalitzPoint(double mA, double mB, double mC,
EvtCyclic3::Pair i, double qres, double qhel, double qsum);
EvtDalitzPoint(const EvtDalitzPlot&, const EvtDalitzCoord&);
- EvtDalitzPoint(const EvtDalitzPoint& other);
EvtDalitzCoord getDalitzPoint(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const;
EvtDalitzPlot getDalitzPlot() const;
double q(EvtCyclic3::Pair) const;
double bigM() const;
double m(EvtCyclic3::Index) const;
// Zemach variables
double qres(EvtCyclic3::Pair i) const;
double qhel(EvtCyclic3::Pair i) const;
double qsum() const;
// Kinematic quantities
//
// pp - 4 momentum product
// e,p,cosTh - energy/moementum in rest-frame of j
double qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const;
double qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const;
double pp(EvtCyclic3::Index i, EvtCyclic3::Index j) const;
double e(EvtCyclic3::Index i, EvtCyclic3::Pair j) const;
double p(EvtCyclic3::Index i, EvtCyclic3::Pair j) const;
double cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const;
bool isValid() const;
void print() const;
private:
double _mA, _mB, _mC; // masses
double _qAB, _qBC, _qCA; // masses squared
};
#endif
diff --git a/EvtGenBase/EvtPdfMax.hh b/EvtGenBase/EvtPdfMax.hh
index 1fee11d..32e7d94 100644
--- a/EvtGenBase/EvtPdfMax.hh
+++ b/EvtGenBase/EvtPdfMax.hh
@@ -1,52 +1,48 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtPdfMax.hh,v 1.2 2009-03-16 16:40:15 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// Pdf maximum and its location
#ifndef EVT_PDF_MAX_HH
#define EVT_PDF_MAX_HH
#include "EvtGenBase/EvtMacros.hh"
// PDF maximum - helper class
template <class Point>
class EvtPdfMax {
public:
EvtPdfMax()
: _value(-1),_valueKnown(false), _locKnown(false)
{}
EvtPdfMax(double value)
: _value(value),_valueKnown(true), _locKnown(false)
{}
EvtPdfMax(Point p, double value)
: _value(value), _valueKnown(true), _locKnown(true), _loc(p)
{}
- EvtPdfMax(const EvtPdfMax& other)
- : COPY_MEM(_value), COPY_MEM(_valueKnown), COPY_MEM(_locKnown), COPY_MEM(_loc)
- {}
- ~EvtPdfMax() {}
-
+
bool valueKnown() const { return _valueKnown; }
double value() const { assert(_valueKnown); return _value; }
bool locKnown() const { return _locKnown; }
Point loc() const { assert(_locKnown); return _loc; }
private:
double _value;
bool _valueKnown;
bool _locKnown;
Point _loc;
};
#endif
diff --git a/EvtGenBase/EvtPropBreitWigner.hh b/EvtGenBase/EvtPropBreitWigner.hh
index 83b2686..757a2dd 100644
--- a/EvtGenBase/EvtPropBreitWigner.hh
+++ b/EvtGenBase/EvtPropBreitWigner.hh
@@ -1,34 +1,33 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtPropBreitWigner.hh,v 1.2 2009-03-16 16:40:16 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// Non-relativistic Breit-Wigner propagator
#ifndef EVT_PROP_BREIT_WIGNER_HH
#define EVT_PROP_BREIT_WIGNER_HH
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtPropagator.hh"
class EvtPropBreitWigner : public EvtPropagator {
public:
EvtPropBreitWigner(double m0, double g0);
- EvtPropBreitWigner(const EvtPropBreitWigner& other);
EvtAmplitude<EvtPoint1D>* clone() const override;
protected:
EvtComplex amplitude(const EvtPoint1D& m) const override;
};
#endif
diff --git a/EvtGenBase/EvtPropBreitWignerRel.hh b/EvtGenBase/EvtPropBreitWignerRel.hh
index 110ada9..b9aecde 100644
--- a/EvtGenBase/EvtPropBreitWignerRel.hh
+++ b/EvtGenBase/EvtPropBreitWignerRel.hh
@@ -1,32 +1,31 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtPropBreitWignerRel.hh,v 1.2 2009-03-16 16:42:03 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// Relativistic Breit-Wigner Propagator
#ifndef EVT_PROP_BREIT_WIGNER_REL_HH
#define EVT_PROP_BREIT_WIGNER_REL_HH
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtPropagator.hh"
class EvtPropBreitWignerRel : public EvtPropagator {
public:
EvtPropBreitWignerRel(double m0, double g0);
- EvtPropBreitWignerRel(const EvtPropBreitWignerRel& other);
EvtAmplitude<EvtPoint1D>* clone() const override;
protected:
EvtComplex amplitude(const EvtPoint1D& x) const override;
};
#endif
diff --git a/EvtGenBase/EvtPropFlatte.hh b/EvtGenBase/EvtPropFlatte.hh
index 666e18c..2f0cdb0 100644
--- a/EvtGenBase/EvtPropFlatte.hh
+++ b/EvtGenBase/EvtPropFlatte.hh
@@ -1,42 +1,41 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* Author: Denis Dujmic, SLAC
*
* Copyright (C) 2005 SLAC
*******************************************************************************/
// Flatte propagator: S.M.Flatte, Phys. Lett. B63, 224 (1976)
#ifndef EVT_PROP_FLATTE_HH
#define EVT_PROP_FLATTE_HH
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtPropagator.hh"
class EvtPropFlatte : public EvtPropagator {
public:
EvtPropFlatte(double m0,
double g0, double m0a, double m0b,
double g1, double m1a, double m1b);
- EvtPropFlatte(const EvtPropFlatte& other);
EvtAmplitude<EvtPoint1D>* clone() const override;
protected:
EvtComplex amplitude(const EvtPoint1D& x) const override;
double _m0a;
double _m0b;
double _g1;
double _m1a;
double _m1b;
};
#endif
diff --git a/EvtGenBase/EvtPropagator.hh b/EvtGenBase/EvtPropagator.hh
index d4dc47b..97fa4b0 100644
--- a/EvtGenBase/EvtPropagator.hh
+++ b/EvtGenBase/EvtPropagator.hh
@@ -1,53 +1,48 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtPropagator.hh,v 1.2 2009-03-16 16:42:03 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// Defines propagator as a function of mass and width
#ifndef EVT_PROPAGATOR_HH
#define EVT_PROPAGATOR_HH
#include <assert.h>
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtAmplitude.hh"
#include "EvtGenBase/EvtPoint1D.hh"
class EvtPropagator : public EvtAmplitude<EvtPoint1D> {
public:
EvtPropagator(double m0, double g0)
: _m0(m0), _g0(g0)
{
assert(m0 > 0);
assert(g0 >= 0);
}
- EvtPropagator(const EvtPropagator& other)
- : EvtAmplitude<EvtPoint1D>(other), _m0(other._m0), _g0(other._g0)
- {}
- virtual ~EvtPropagator()
- {}
-
+
// Accessors
inline double m0() const { return _m0; }
inline double g0() const { return _g0; }
// Modifiers (can be useful e.g. for fitting!)
inline void set_m0(double m0) { assert(m0>0); _m0 = m0; }
inline void set_g0(double g0) { assert(g0>=0); _g0 = g0; }
protected:
double _m0;
double _g0;
};
#endif
diff --git a/EvtGenBase/EvtTwoBodyKine.hh b/EvtGenBase/EvtTwoBodyKine.hh
index 3bbf4b6..55fae63 100644
--- a/EvtGenBase/EvtTwoBodyKine.hh
+++ b/EvtGenBase/EvtTwoBodyKine.hh
@@ -1,55 +1,53 @@
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtTwoBodyKine.hh,v 1.2 2009-03-16 16:34:38 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
// Descriptions of the kinematics of a two-body decay.
#ifndef EVT_TWO_BODY_KINE_HH
#define EVT_TWO_BODY_KINE_HH
#include <iostream>
class EvtTwoBodyKine {
public:
enum Index {A,B,AB};
EvtTwoBodyKine();
EvtTwoBodyKine(double mA, double mB, double mAB);
- EvtTwoBodyKine(const EvtTwoBodyKine& other);
- ~EvtTwoBodyKine();
// Accessors
inline double mA() const { return _mA; }
inline double mB() const { return _mB; }
inline double mAB() const { return _mAB; }
double m(Index i) const;
// Momentum of the other two particles in the
// rest-frame of particle i.
double p(Index i = AB) const;
// Energy of particle i in the rest frame of particle j
double e(Index i, Index j) const;
void print(std::ostream& os) const;
private:
double _mA;
double _mB;
double _mAB;
};
std::ostream& operator<<(std::ostream& os, const EvtTwoBodyKine& p);
#endif
diff --git a/EvtGenBase/EvtVector4C.hh b/EvtGenBase/EvtVector4C.hh
index d8e75a2..061884e 100644
--- a/EvtGenBase/EvtVector4C.hh
+++ b/EvtGenBase/EvtVector4C.hh
@@ -1,211 +1,201 @@
//--------------------------------------------------------------------------
//
// Environment:
// This software is part of the EvtGen package developed jointly
// for the BaBar and CLEO collaborations. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
// Copyright (C) 1998 Caltech, UCSB
//
// Module: EvtGen/EvtVector4C.hh
//
// Description: Class for complex 4 vectors
//
// Modification history:
//
// DJL/RYD September 25, 1996 Module created
//
//------------------------------------------------------------------------
#ifndef EVTVECTOR4C_HH
#define EVTVECTOR4C_HH
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtVector3C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include <iosfwd>
class EvtVector4C final {
- friend EvtVector4C rotateEuler(const EvtVector4C& e,
- double alpha,double beta,double gamma);
- friend EvtVector4C boostTo(const EvtVector4C& e,
- const EvtVector4R p4);
- friend EvtVector4C boostTo(const EvtVector4C& e,
- const EvtVector3R boost);
inline friend EvtVector4C operator*(double d,const EvtVector4C& v2);
inline friend EvtVector4C operator*(const EvtComplex& c,const EvtVector4C& v2);
inline friend EvtVector4C operator*(const EvtVector4C& v2,const EvtComplex& c);
inline friend EvtVector4C operator*(const EvtComplex& c,const EvtVector4R& v2);
inline friend EvtComplex operator*(const EvtVector4R& v1,const EvtVector4C& v2);
inline friend EvtComplex operator*(const EvtVector4C& v1,const EvtVector4R& v2);
inline friend EvtComplex operator*(const EvtVector4C& v1,const EvtVector4C& v2);
friend EvtVector4C operator+(const EvtVector4C& v1,const EvtVector4C& v2);
friend EvtVector4C operator-(const EvtVector4C& v1,const EvtVector4C& v2);
public:
EvtVector4C();
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&);
inline void set(double,double,double,double);
inline EvtVector4C(const EvtVector4R& v1);
inline const EvtComplex& get(int) const;
inline EvtComplex cont(const EvtVector4C& v4) const;
inline EvtVector4C conj() const;
EvtVector3C vec() const;
- inline EvtVector4C& operator=(const EvtVector4C& v2);
inline EvtVector4C& operator-=(const EvtVector4C& v2);
inline EvtVector4C& operator+=(const EvtVector4C& v2);
inline EvtVector4C& operator*=(const EvtComplex& c);
void applyRotateEuler(double alpha,double beta,double gamma);
void applyBoostTo(const EvtVector4R& p4);
void applyBoostTo(const EvtVector3R& boost);
friend std::ostream& operator<<(std::ostream& s, const EvtVector4C& v);
double dot( const EvtVector4C& p2 );
private:
EvtComplex v[4];
};
-inline EvtVector4C& EvtVector4C::operator=(const EvtVector4C& v2){
-
- v[0]=v2.v[0];
- v[1]=v2.v[1];
- v[2]=v2.v[2];
- v[3]=v2.v[3];
-
- return *this;
-}
+EvtVector4C rotateEuler(const EvtVector4C& e,
+ double alpha,double beta,double gamma);
+EvtVector4C boostTo(const EvtVector4C& e,
+ const EvtVector4R p4);
+EvtVector4C boostTo(const EvtVector4C& e,
+ const EvtVector3R boost);
inline EvtVector4C& EvtVector4C::operator+=(const EvtVector4C& v2){
v[0]+=v2.v[0];
v[1]+=v2.v[1];
v[2]+=v2.v[2];
v[3]+=v2.v[3];
return *this;
}
inline EvtVector4C& EvtVector4C::operator-=(const EvtVector4C& v2){
v[0]-=v2.v[0];
v[1]-=v2.v[1];
v[2]-=v2.v[2];
v[3]-=v2.v[3];
return *this;
}
inline void EvtVector4C::set(int i,const EvtComplex& c){
v[i]=c;
}
inline EvtVector3C EvtVector4C::vec() const {
return EvtVector3C(v[1],v[2],v[3]);
}
inline void EvtVector4C::set(const EvtComplex& e,const EvtComplex& p1,
const EvtComplex& p2,const EvtComplex& p3){
v[0]=e; v[1]=p1; v[2]=p2; v[3]=p3;
}
inline void EvtVector4C::set(double e,double p1,
double p2,double p3){
v[0]=EvtComplex(e); v[1]=EvtComplex(p1); v[2]=EvtComplex(p2); v[3]=EvtComplex(p3);
}
inline const EvtComplex& EvtVector4C::get(int i) const {
return v[i];
}
inline EvtVector4C operator+(const EvtVector4C& v1,const EvtVector4C& v2) {
return EvtVector4C(v1)+=v2;
}
inline EvtVector4C operator-(const EvtVector4C& v1,const EvtVector4C& v2) {
return EvtVector4C(v1)-=v2;
}
inline EvtComplex EvtVector4C::cont(const EvtVector4C& v4) const {
return v[0]*v4.v[0]-v[1]*v4.v[1]-
v[2]*v4.v[2]-v[3]*v4.v[3];
}
inline EvtVector4C& EvtVector4C::operator*=(const EvtComplex& c) {
v[0]*=c;
v[1]*=c;
v[2]*=c;
v[3]*=c;
return *this;
}
inline EvtVector4C operator*(double d,const EvtVector4C& v2){
return EvtVector4C(v2.v[0]*d,v2.v[1]*d,v2.v[2]*d,v2.v[3]*d);
}
inline EvtVector4C operator*(const EvtComplex& c,const EvtVector4C& v2){
return EvtVector4C(v2)*=c;
}
inline EvtVector4C operator*(const EvtVector4C& v2,const EvtComplex& c){
return EvtVector4C(v2)*=c;
}
inline EvtVector4C operator*(const EvtComplex& c,const EvtVector4R& v2){
return EvtVector4C(c*v2.get(0),c*v2.get(1),c*v2.get(2),c*v2.get(3));
}
inline EvtVector4C::EvtVector4C(const EvtVector4R& v1){
v[0]=EvtComplex(v1.get(0)); v[1]=EvtComplex(v1.get(1));
v[2]=EvtComplex(v1.get(2)); v[3]=EvtComplex(v1.get(3));
}
inline EvtComplex operator*(const EvtVector4R& v1,const EvtVector4C& v2){
return v1.get(0)*v2.v[0]-v1.get(1)*v2.v[1]-
v1.get(2)*v2.v[2]-v1.get(3)*v2.v[3];
}
inline EvtComplex operator*(const EvtVector4C& v1,const EvtVector4R& v2){
return v1.v[0]*v2.get(0)-v1.v[1]*v2.get(1)-
v1.v[2]*v2.get(2)-v1.v[3]*v2.get(3);
}
inline EvtComplex operator*(const EvtVector4C& v1,const EvtVector4C& v2){
return v1.v[0]*v2.v[0]-v1.v[1]*v2.v[1]-
v1.v[2]*v2.v[2]-v1.v[3]*v2.v[3];
}
inline EvtVector4C EvtVector4C::conj() const {
return EvtVector4C(::conj(v[0]),::conj(v[1]),
::conj(v[2]),::conj(v[3]));
}
#endif
diff --git a/EvtGenBase/EvtVector4R.hh b/EvtGenBase/EvtVector4R.hh
index 73fa602..0fcd696 100644
--- a/EvtGenBase/EvtVector4R.hh
+++ b/EvtGenBase/EvtVector4R.hh
@@ -1,204 +1,183 @@
//--------------------------------------------------------------------------
//
// Environment:
// This software is part of the EvtGen package developed jointly
// for the BaBar and CLEO collaborations. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
// Copyright (C) 1998 Caltech, UCSB
//
// Module: EvtGen/EvtVector4R.hh
//
// Description: Class to describe real 4 vectors
//
// Modification history:
//
// DJL/RYD September 25, 1996 Module created
//
//------------------------------------------------------------------------
#ifndef EVTVECTOR4R_HH
#define EVTVECTOR4R_HH
#include <iostream>
#include <math.h>
class EvtVector3R;
-class EvtVector4R;
-
-EvtVector4R rotateEuler(const EvtVector4R& rs,
- double alpha,double beta,double gamma);
-EvtVector4R boostTo(const EvtVector4R& rs,
- const EvtVector4R& p4, bool inverse = false);
-EvtVector4R boostTo(const EvtVector4R& rs,
- const EvtVector3R& boost, bool inverse = false);
class EvtVector4R {
- friend EvtVector4R rotateEuler(const EvtVector4R& rs,
- double alpha,double beta,double gamma);
- friend EvtVector4R boostTo(const EvtVector4R& rs,
- const EvtVector4R& p4, bool inverse);
- friend EvtVector4R boostTo(const EvtVector4R& rs,
- const EvtVector3R& boost, bool inverse);
-
-
inline friend EvtVector4R operator*(double d,const EvtVector4R& v2);
inline friend EvtVector4R operator*(const EvtVector4R& v2,double d);
inline friend EvtVector4R operator/(const EvtVector4R& v2,double d);
inline friend double operator*(const EvtVector4R& v1,const EvtVector4R& v2);
inline friend EvtVector4R operator+(const EvtVector4R& v1,const EvtVector4R& v2);
inline friend EvtVector4R operator-(const EvtVector4R& v1,const EvtVector4R& v2);
public:
EvtVector4R();
EvtVector4R(double e,double px,double py ,double pz);
inline void set(int i,double d);
inline void set(double e,double px,double py ,double pz);
inline EvtVector4R& operator*=(double c);
inline EvtVector4R& operator/=(double c);
- inline EvtVector4R& operator=(const EvtVector4R& v2);
inline EvtVector4R& operator+=(const EvtVector4R& v2);
inline EvtVector4R& operator-=(const EvtVector4R& v2);
inline double get(int i) const;
inline double cont(const EvtVector4R& v4) const;
friend std::ostream& operator<<(std::ostream& s, const EvtVector4R& v);
double mass2() const;
double mass() const;
void applyRotateEuler(double alpha,double beta,double gamma);
void applyBoostTo(const EvtVector4R& p4, bool inverse = false);
void applyBoostTo(const EvtVector3R& boost, bool inverse = false);
EvtVector4R cross(const EvtVector4R& v2);
double dot(const EvtVector4R& v2) const;
double d3mag() const;
// Added by AJB - calculate scalars in the rest frame of the current object
double scalartripler3( const EvtVector4R& p1, const EvtVector4R& p2,
const EvtVector4R& p3 ) const;
double dotr3( const EvtVector4R& p1, const EvtVector4R& p2 ) const;
double mag2r3( const EvtVector4R& p1 ) const;
double magr3( const EvtVector4R& p1 ) const;
private:
double v[4];
inline double Square( double x ) const { return x*x; }
};
-
-inline EvtVector4R& EvtVector4R::operator=(const EvtVector4R& v2){
-
- v[0]=v2.v[0];
- v[1]=v2.v[1];
- v[2]=v2.v[2];
- v[3]=v2.v[3];
-
- return *this;
-}
+EvtVector4R rotateEuler(const EvtVector4R& rs,
+ double alpha,double beta,double gamma);
+EvtVector4R boostTo(const EvtVector4R& rs,
+ const EvtVector4R& p4, bool inverse = false);
+EvtVector4R boostTo(const EvtVector4R& rs,
+ const EvtVector3R& boost, bool inverse = false);
inline EvtVector4R& EvtVector4R::operator+=(const EvtVector4R& v2){
v[0]+=v2.v[0];
v[1]+=v2.v[1];
v[2]+=v2.v[2];
v[3]+=v2.v[3];
return *this;
}
inline EvtVector4R& EvtVector4R::operator-=(const EvtVector4R& v2){
v[0]-=v2.v[0];
v[1]-=v2.v[1];
v[2]-=v2.v[2];
v[3]-=v2.v[3];
return *this;
}
inline double EvtVector4R::mass2() const{
return v[0]*v[0]-v[1]*v[1]-v[2]*v[2]-v[3]*v[3];
}
inline EvtVector4R operator*(double c,const EvtVector4R& v2){
return EvtVector4R(v2)*=c;
}
inline EvtVector4R operator*(const EvtVector4R& v2,double c){
return EvtVector4R(v2)*=c;
}
inline EvtVector4R operator/(const EvtVector4R& v2,double c){
return EvtVector4R(v2)/=c;
}
inline EvtVector4R& EvtVector4R::operator*=(double c){
v[0]*=c;
v[1]*=c;
v[2]*=c;
v[3]*=c;
return *this;
}
inline EvtVector4R& EvtVector4R::operator/=(double c){
double cinv=1.0/c;
v[0]*=cinv;
v[1]*=cinv;
v[2]*=cinv;
v[3]*=cinv;
return *this;
}
inline double operator*(const EvtVector4R& v1,const EvtVector4R& v2){
return v1.v[0]*v2.v[0]-v1.v[1]*v2.v[1]-
v1.v[2]*v2.v[2]-v1.v[3]*v2.v[3];
}
inline double EvtVector4R::cont(const EvtVector4R& v4) const {
return v[0]*v4.v[0]-v[1]*v4.v[1]-
v[2]*v4.v[2]-v[3]*v4.v[3];
}
inline EvtVector4R operator-(const EvtVector4R& v1,const EvtVector4R& v2){
return EvtVector4R(v1)-=v2;
}
inline EvtVector4R operator+(const EvtVector4R& v1,const EvtVector4R& v2){
return EvtVector4R(v1)+=v2;
}
inline double EvtVector4R::get(int i) const {
return v[i];
}
inline void EvtVector4R::set(int i,double d){
v[i]=d;
}
inline void EvtVector4R::set(double e,double p1,double p2, double p3){
v[0]=e;
v[1]=p1;
v[2]=p2;
v[3]=p3;
}
#endif
diff --git a/src/EvtGenBase/EvtDalitzPlot.cpp b/src/EvtGenBase/EvtDalitzPlot.cpp
index fe540c1..b95b254 100644
--- a/src/EvtGenBase/EvtDalitzPlot.cpp
+++ b/src/EvtGenBase/EvtDalitzPlot.cpp
@@ -1,353 +1,342 @@
//-----------------------------------------------------------------------
// File and Version Information:
// $Id: EvtDalitzPlot.cpp,v 1.3 2009-03-16 15:53:27 robbep Exp $
//
// Environment:
// This software is part of the EvtGen package developed jointly
// for the BaBar and CLEO collaborations. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information:
// Copyright (C) 1998 Caltech, UCSB
//
// Module creator:
// Alexei Dvoretskii, Caltech, 2001-2002.
//-----------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"
// Global 3-body Dalitz decay kinematics as defined by the mass
// of the mother and the daughters. Spins are not considered.
#include <math.h>
#include <assert.h>
#include <stdio.h>
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtDalitzPlot.hh"
#include "EvtGenBase/EvtTwoBodyVertex.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtDecayMode.hh"
using namespace EvtCyclic3;
EvtDalitzPlot::EvtDalitzPlot()
: _mA(0.), _mB(0.), _mC(0.), _bigM(0.),
_ldel(0.), _rdel(0.)
{}
EvtDalitzPlot::EvtDalitzPlot(double mA, double mB, double mC, double bigM,
double ldel, double rdel)
: _mA(mA), _mB(mB), _mC(mC), _bigM(bigM),
_ldel(ldel), _rdel(rdel)
{
sanityCheck();
}
EvtDalitzPlot::EvtDalitzPlot(const EvtDecayMode& mode, double ldel, double rdel )
{
_mA = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(A)));
_mB = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(B)));
_mC = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(C)));
_bigM = EvtPDL::getMeanMass(EvtPDL::getId(mode.mother()));
_ldel = ldel;
_rdel = rdel;
sanityCheck();
}
-
-EvtDalitzPlot::EvtDalitzPlot(const EvtDalitzPlot& other)
- : _mA(other._mA), _mB(other._mB), _mC(other._mC), _bigM(other._bigM),
- _ldel(other._ldel), _rdel(other._rdel)
-{}
-
-
-EvtDalitzPlot::~EvtDalitzPlot()
-{}
-
-
bool EvtDalitzPlot::operator==(const EvtDalitzPlot& other) const
{
bool ret = false;
if(_mA == other._mA &&
_mB == other._mB &&
_mC == other._mC &&
_bigM == other._bigM) ret = true;
return ret;
}
const EvtDalitzPlot* EvtDalitzPlot::clone() const
{
return new EvtDalitzPlot(*this);
}
void EvtDalitzPlot::sanityCheck() const
{
if(_mA < 0 || _mB < 0 || _mC < 0 || _bigM <= 0 || _bigM - _mA - _mB - _mC < 0.) {
printf("Invalid Dalitz plot %f %f %f %f\n",_mA,_mB,_mC,_bigM);
assert(0);
}
assert(_ldel <= 0.);
assert(_rdel >= 0.);
}
double EvtDalitzPlot::m(Index i) const {
double m = _mA;
if(i == B) m = _mB;
else
if(i == C) m = _mC;
return m;
}
double EvtDalitzPlot::sum() const
{
return _mA*_mA + _mB*_mB + _mC*_mC + _bigM*_bigM;
}
double EvtDalitzPlot::qAbsMin(Pair i) const
{
Index j = first(i);
Index k = second(i);
return (m(j) + m(k))*(m(j) + m(k));
}
double EvtDalitzPlot::qAbsMax(Pair i) const
{
Index j = other(i);
return (_bigM-m(j))*(_bigM-m(j));
}
double EvtDalitzPlot::qResAbsMin(EvtCyclic3::Pair i) const
{
return qAbsMin(i) - sum()/3.;
}
double EvtDalitzPlot::qResAbsMax(EvtCyclic3::Pair i) const
{
return qAbsMax(i) - sum()/3.;
}
double EvtDalitzPlot::qHelAbsMin(EvtCyclic3::Pair i) const
{
Pair j = next(i);
Pair k = prev(i);
return (qAbsMin(j) - qAbsMax(k))/2.;
}
double EvtDalitzPlot::qHelAbsMax(EvtCyclic3::Pair i) const
{
Pair j = next(i);
Pair k = prev(i);
return (qAbsMax(j) - qAbsMin(k))/2.;
}
double EvtDalitzPlot::mAbsMin(Pair i) const
{
return sqrt(qAbsMin(i));
}
double EvtDalitzPlot::mAbsMax(Pair i) const
{
return sqrt(qAbsMax(i));
}
// parallel
double EvtDalitzPlot::qMin(Pair i, Pair j, double q) const
{
if(i == j) return q;
else {
// Particle pair j defines the rest-frame
// 0 - particle common to r.f. and angle calculations
// 1 - particle belonging to r.f. but not angle
// 2 - particle not belonging to r.f.
Index k0 = common(i,j);
Index k2 = other(j);
Index k1 = other(k0,k2);
// Energy, momentum of particle common to rest-frame and angle
EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
double pk = jpair.p();
double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
// Energy and momentum of the other particle
EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
double pj = mother.p(EvtTwoBodyKine::A);
// See PDG 34.4.3.1
return (ek+ej)*(ek+ej) - (pk+pj)*(pk+pj);
}
}
// antiparallel
double EvtDalitzPlot::qMax(Pair i, Pair j, double q) const
{
if(i == j) return q;
else {
// Particle pair j defines the rest-frame
// 0 - particle common to r.f. and angle calculations
// 1 - particle belonging to r.f. but not angle
// 2 - particle not belonging to r.f.
Index k0 = common(i,j);
Index k2 = other(j);
Index k1 = other(k0,k2);
// Energy, momentum of particle common to rest-frame and angle
EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
double pk = jpair.p();
// Energy and momentum of the other particle
EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
double pj = mother.p(EvtTwoBodyKine::A);
// See PDG 34.4.3.1
return (ek+ej)*(ek+ej) - (pk-pj)*(pk-pj);
}
}
double EvtDalitzPlot::getArea(int N, Pair i, Pair j) const
{
// Trapezoidal integral over qi. qj can be calculated.
// The first and the last point are zero, so they are not counted
double dh = (qAbsMax(i) - qAbsMin(i))/((double) N);
double sum = 0;
int ii;
for(ii=1;ii<N;ii++) {
double x = qAbsMin(i) + ii*dh;
double dy = qMax(j,i,x) - qMin(j,i,x);
sum += dy;
}
return sum * dh;
}
double EvtDalitzPlot::cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
{
if(i1 == i2) return 1.;
double qmax = qMax(i1,i2,q2);
double qmin = qMin(i1,i2,q2);
double cos = (qmax + qmin - 2*q1)/(qmax - qmin);
return cos;
}
double EvtDalitzPlot::e(Index i, Pair j, double q) const
{
if(i == other(j)) {
// i does not belong to pair j
return (bigM()*bigM()-q-m(i)*m(i))/2/sqrt(q);
}
else {
// i and k make pair j
Index k;
if(first(j) == i) k = second(j);
else k = first(j);
double e = (q + m(i)*m(i) - m(k)*m(k))/2/sqrt(q);
return e;
}
}
double EvtDalitzPlot::p(Index i, Pair j, double q) const
{
double en = e(i,j,q);
double p2 = en*en - m(i)*m(i);
if(p2 < 0) {
printf("Bad value of p2 %f %d %d %f %f\n",p2,i,j,en,m(i));
assert(0);
}
return sqrt(p2);
}
double EvtDalitzPlot::q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const
{
if(i1 == i2) return q2;
EvtCyclic3::Index f = first(i1);
EvtCyclic3::Index s = second(i1);
return m(f)*m(f) + m(s)*m(s) + 2*e(f,i2,q2)*e(s,i2,q2) - 2*p(f,i2,q2)*p(s,i2,q2)*cosTh;
}
double EvtDalitzPlot::jacobian(EvtCyclic3::Pair i, double q) const
{
return 2*p(first(i),i,q)*p(other(i),i,q); // J(BC) = 2pA*pB = 2pA*pC
}
EvtTwoBodyVertex EvtDalitzPlot::vD(Pair iRes, double m0, int L) const
{
return EvtTwoBodyVertex(m(first(iRes)),
m(second(iRes)),m0,L);
}
EvtTwoBodyVertex EvtDalitzPlot::vB(Pair iRes, double m0, int L) const
{
return EvtTwoBodyVertex(m0,m(other(iRes)),bigM(),L);
}
void EvtDalitzPlot::print() const
{
// masses
printf("Mass M %f\n",bigM());
printf("Mass mA %f\n",_mA);
printf("Mass mB %f\n",_mB);
printf("Mass mC %f\n",_mC);
// limits
printf("Limits qAB %f : %f\n",qAbsMin(AB),qAbsMax(AB));
printf("Limits qBC %f : %f\n",qAbsMin(BC),qAbsMax(BC));
printf("Limits qCA %f : %f\n",qAbsMin(CA),qAbsMax(CA));
printf("Sum q %f\n",sum());
printf("Limit qsum %f : %f\n",qSumMin(),qSumMax());
}
diff --git a/src/EvtGenBase/EvtDalitzPoint.cpp b/src/EvtGenBase/EvtDalitzPoint.cpp
index be5312a..8c3c227 100644
--- a/src/EvtGenBase/EvtDalitzPoint.cpp
+++ b/src/EvtGenBase/EvtDalitzPoint.cpp
@@ -1,185 +1,180 @@
#include "EvtGenBase/EvtPatches.hh"
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtDalitzPoint.cpp,v 1.3 2009-03-16 15:53:27 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
#include "EvtGenBase/EvtPatches.hh"
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "EvtGenBase/EvtDalitzPoint.hh"
using namespace EvtCyclic3;
EvtDalitzPoint::EvtDalitzPoint()
: _mA(-1.), _mB(-1.), _mC(-1.), _qAB(-1.), _qBC(-1.), _qCA(-1.)
{}
EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC, double qAB, double qBC, double qCA)
: _mA(mA), _mB(mB), _mC(mC), _qAB(qAB), _qBC(qBC), _qCA(qCA)
{}
// Constructor from Zemach coordinates
EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC,
EvtCyclic3::Pair i,
double qres, double qhel, double qsum)
: _mA(mA), _mB(mB), _mC(mC)
{
double qi = qres + qsum/3.;
double qj = -qres/2. + qhel + qsum/3.;
double qk = -qres/2. - qhel + qsum/3.;
if(i == AB) { _qAB = qi; _qBC = qj; _qCA = qk; }
else if(i == BC) { _qAB = qk; _qBC = qi; _qCA = qj; }
else if(i == CA) { _qAB = qj; _qBC = qk; _qCA = qi; }
}
EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPlot& dp, const EvtDalitzCoord& x)
: _mA(dp.m(A)), _mB(dp.m(B)), _mC(dp.m(C))
{
if(x.pair1() == AB) _qAB = x.q1();
else
if(x.pair2() == AB) _qAB = x.q2();
else _qAB = dp.sum() - x.q1() - x.q2();
if(x.pair1() == BC) _qBC = x.q1();
else
if(x.pair2() == BC) _qBC = x.q2();
else _qBC = dp.sum() - x.q1() - x.q2();
if(x.pair1() == CA) _qCA = x.q1();
else
if(x.pair2() == CA) _qCA = x.q2();
else _qCA = dp.sum() - x.q1() - x.q2();
}
-EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPoint& other)
- : _mA(other._mA), _mB(other._mB), _mC(other._mC),
- _qAB(other._qAB), _qBC(other._qBC), _qCA(other._qCA)
-{}
-
double EvtDalitzPoint::q(EvtCyclic3::Pair i) const
{
double ret = _qAB;
if(BC == i) ret = _qBC;
else
if(CA == i) ret = _qCA;
return ret;
}
double EvtDalitzPoint::m(EvtCyclic3::Index i) const
{
double ret = _mA;
if(B == i) ret = _mB;
else
if(C == i) ret = _mC;
return ret;
}
// Zemach variables
double EvtDalitzPoint::qres(EvtCyclic3::Pair i) const
{
return (2.*q(i) - q(EvtCyclic3::prev(i)) - q(EvtCyclic3::next(i)))/3.;
}
double EvtDalitzPoint::qhel(EvtCyclic3::Pair i) const
{
Pair j = next(i);
Pair k = prev(i);
return (q(j) - q(k))/2.;
}
double EvtDalitzPoint::qsum() const
{
return _qAB + _qBC + _qCA;
}
double EvtDalitzPoint::qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
{
EvtDalitzPlot dp = getDalitzPlot();
return dp.qMin(i,j,q(j));
}
double EvtDalitzPoint::qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
{
EvtDalitzPlot dp = getDalitzPlot();
return dp.qMax(i,j,q(j));
}
double EvtDalitzPoint::pp(EvtCyclic3::Index i, EvtCyclic3::Index j) const
{
if(i == j) return m(i)*m(i);
else return (q(combine(i,j)) - m(i)*m(i) - m(j)*m(j))/2.;
}
double EvtDalitzPoint::e(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
{
EvtDalitzPlot dp = getDalitzPlot();
return dp.e(i,j,q(j));
}
double EvtDalitzPoint::p(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
{
EvtDalitzPlot dp = getDalitzPlot();
return dp.p(i,j,q(j));
}
double EvtDalitzPoint::cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const
{
EvtDalitzPlot dp = getDalitzPlot();
return dp.cosTh(pairAng,q(pairAng),pairRes,q(pairRes));
}
EvtDalitzCoord EvtDalitzPoint::getDalitzPoint(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
{
return EvtDalitzCoord(i,q(i),j,q(j));
}
EvtDalitzPlot EvtDalitzPoint::getDalitzPlot() const
{
return EvtDalitzPlot(_mA,_mB,_mC,bigM());
}
bool EvtDalitzPoint::isValid() const
{
// Check masses
double M = bigM();
if(_mA < 0 || _mB < 0 || _mC < 0 || M <= 0) return false;
if(M < _mA + _mB + _mC) return false;
// Check that first coordinate is within absolute limits
bool inside = false;
EvtDalitzPlot dp = getDalitzPlot();
if(dp.qAbsMin(AB) <= _qAB && _qAB <= dp.qAbsMax(AB))
if(qMin(BC,AB) <= _qBC && _qBC <= qMax(BC,AB))
inside = true;
return inside;
}
double EvtDalitzPoint::bigM() const
{
return sqrt(_qAB+_qBC+_qCA - _mA*_mA - _mB*_mB - _mC*_mC);
}
void EvtDalitzPoint::print() const
{
getDalitzPlot().print();
printf("%f %f %f\n",_qAB,_qBC,_qCA);
}
diff --git a/src/EvtGenBase/EvtPropBreitWigner.cpp b/src/EvtGenBase/EvtPropBreitWigner.cpp
index 5c1ade9..8a98547 100644
--- a/src/EvtGenBase/EvtPropBreitWigner.cpp
+++ b/src/EvtGenBase/EvtPropBreitWigner.cpp
@@ -1,37 +1,32 @@
#include "EvtGenBase/EvtPatches.hh"
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtPropBreitWigner.cpp,v 1.3 2009-03-16 15:44:41 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
#include <math.h>
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtPropBreitWigner.hh"
EvtPropBreitWigner::EvtPropBreitWigner(double m0, double g0)
: EvtPropagator(m0,g0)
{}
-EvtPropBreitWigner::EvtPropBreitWigner(const EvtPropBreitWigner& other)
- : EvtPropagator(other)
-{}
-
-
EvtAmplitude<EvtPoint1D>* EvtPropBreitWigner::clone() const
{
return new EvtPropBreitWigner(*this);
}
EvtComplex EvtPropBreitWigner::amplitude(const EvtPoint1D& x) const
{
double m = x.value();
EvtComplex value = sqrt(_g0/EvtConst::twoPi)/(m-_m0-EvtComplex(0.0,_g0/2.));
return value;
}
diff --git a/src/EvtGenBase/EvtPropBreitWignerRel.cpp b/src/EvtGenBase/EvtPropBreitWignerRel.cpp
index aa75d68..6c2e4af 100644
--- a/src/EvtGenBase/EvtPropBreitWignerRel.cpp
+++ b/src/EvtGenBase/EvtPropBreitWignerRel.cpp
@@ -1,35 +1,30 @@
#include "EvtGenBase/EvtPatches.hh"
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtPropBreitWignerRel.cpp,v 1.3 2009-03-16 15:44:41 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
#include <math.h>
#include "EvtGenBase/EvtPropBreitWignerRel.hh"
EvtPropBreitWignerRel::EvtPropBreitWignerRel(double m0, double g0)
: EvtPropagator(m0,g0)
{}
-EvtPropBreitWignerRel::EvtPropBreitWignerRel(const EvtPropBreitWignerRel& other)
- : EvtPropagator(other)
-{}
-
-
EvtAmplitude<EvtPoint1D>* EvtPropBreitWignerRel::clone() const
{
return new EvtPropBreitWignerRel(*this);
}
EvtComplex EvtPropBreitWignerRel::amplitude(const EvtPoint1D& x) const
{
double m = x.value();
return 1./(_m0*_m0-m*m-EvtComplex(0.,_m0*_g0));
}
diff --git a/src/EvtGenBase/EvtPropFlatte.cpp b/src/EvtGenBase/EvtPropFlatte.cpp
index 2d51631..fa084b0 100644
--- a/src/EvtGenBase/EvtPropFlatte.cpp
+++ b/src/EvtGenBase/EvtPropFlatte.cpp
@@ -1,93 +1,83 @@
#include "EvtGenBase/EvtPatches.hh"
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* Author : D. Dujmic, J. Thompson
*
* Copyright (C) 2005 SLAC
*******************************************************************************/
#include <math.h>
#include "EvtGenBase/EvtPropFlatte.hh"
#include <iostream>
using std::cout;
using std::endl;
EvtPropFlatte::EvtPropFlatte(double m0,
double g0, double m0a, double m0b,
double g1, double m1a, double m1b) :
EvtPropagator( m0, g0),
_m0a(m0a),
_m0b(m0b),
_g1 (g1),
_m1a(m1a),
_m1b(m1b)
{}
-EvtPropFlatte::EvtPropFlatte(const EvtPropFlatte& other) :
- EvtPropagator(other),
- _m0a (other._m0a),
- _m0b (other._m0b),
- _g1 (other._g1),
- _m1a (other._m1a),
- _m1b (other._m1b)
-{}
-
-
EvtAmplitude<EvtPoint1D>* EvtPropFlatte::clone() const
{
return new EvtPropFlatte(*this);
}
EvtComplex EvtPropFlatte::amplitude(const EvtPoint1D& x) const
{
/*
Use BES parameterization:
1.
-----------------------------------------
m0^2 - m^2 - i*m0*( g1*rho1 + g2*rho2 )
Resonance mass: m0
Channel1: m0a, m0b, g0
Channel2: m1a, m1b, g1
where breakup momenta q's are:
E0a = (m^2 + m0a^2 - m0b^2) / 2m
q0 = sqrt( E0a^2 - m0a^2 )
E1a = (m^2 + m1a^2 - m1b^2) / 2m
q1 = sqrt( E1a^2 - m1a^2 )
*/
double s = x.value()*x.value();
double m = x.value();
double E0a = 0.5 * (s + _m0a*_m0a - _m0b*_m0b) / m;
double qSq0 = E0a*E0a - _m0a*_m0a;
double E1a = 0.5 * (s + _m1a*_m1a - _m1b*_m1b) / m;
double qSq1 = E1a*E1a - _m1a*_m1a;
EvtComplex gamma0 = qSq0 >= 0 ? EvtComplex( _g0 * sqrt(qSq0), 0) : EvtComplex( 0, _g0 * sqrt(-qSq0) );
EvtComplex gamma1 = qSq1 >= 0 ? EvtComplex( _g1 * sqrt(qSq1), 0) : EvtComplex( 0, _g1 * sqrt(-qSq1) );
EvtComplex gamma = gamma0 + gamma1;
EvtComplex a = 1.0/( _m0*_m0 - s - EvtComplex(0.0,2*_m0/m)*gamma );
return a;
}
diff --git a/src/EvtGenBase/EvtTwoBodyKine.cpp b/src/EvtGenBase/EvtTwoBodyKine.cpp
index 5be2837..b59644c 100644
--- a/src/EvtGenBase/EvtTwoBodyKine.cpp
+++ b/src/EvtGenBase/EvtTwoBodyKine.cpp
@@ -1,105 +1,98 @@
#include "EvtGenBase/EvtPatches.hh"
/*******************************************************************************
* Project: BaBar detector at the SLAC PEP-II B-factory
* Package: EvtGenBase
* File: $Id: EvtTwoBodyKine.cpp,v 1.3 2009-03-16 15:37:54 robbep Exp $
* Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
*
* Copyright (C) 2002 Caltech
*******************************************************************************/
#include <iostream>
#include <assert.h>
#include <math.h>
#include "EvtGenBase/EvtTwoBodyKine.hh"
#include "EvtGenBase/EvtReport.hh"
using std::endl;
using std::ostream;
EvtTwoBodyKine::EvtTwoBodyKine()
: _mA(0.), _mB(0.), _mAB(0.)
{}
EvtTwoBodyKine::EvtTwoBodyKine(double mA, double mB, double mAB)
: _mA(mA), _mB(mB), _mAB(mAB)
{
if(mAB < mA + mB) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << mAB << " < " << mA << " + " << mB << endl;
assert(0);
}
}
-EvtTwoBodyKine::EvtTwoBodyKine(const EvtTwoBodyKine& other)
- : _mA(other._mA), _mB(other._mB), _mAB(other._mAB)
-{}
-
-EvtTwoBodyKine::~EvtTwoBodyKine()
-{}
-
double EvtTwoBodyKine::m(Index i) const
{
double ret = _mAB;
if(A == i) ret = _mA;
else
if(B == i) ret = _mB;
return ret;
}
double EvtTwoBodyKine::p(Index i) const
{
double p0 = 0.;
if(i == AB) {
double x = _mAB*_mAB - _mA*_mA - _mB*_mB;
double y = 2*_mA*_mB;
p0 = sqrt(x*x - y*y)/2./_mAB;
}
else
if(i == A) {
double x = _mA*_mA - _mAB*_mAB - _mB*_mB;
double y = 2*_mAB*_mB;
p0 = sqrt(x*x - y*y)/2./_mA;
}
else {
double x = _mB*_mB - _mAB*_mAB - _mA*_mA;
double y = 2*_mAB*_mA;
p0 = sqrt(x*x - y*y)/2./_mB;
}
return p0;
}
double EvtTwoBodyKine::e(Index i, Index j) const
{
double ret = m(i);
if(i != j) {
double pD = p(j);
ret = sqrt(ret*ret + pD*pD);
}
return ret;
}
void EvtTwoBodyKine::print(ostream& os) const
{
os << " mA = " << _mA << endl;
os << " mB = " << _mB << endl;
os << "mAB = " << _mAB << endl;
}
ostream& operator<<(ostream& os, const EvtTwoBodyKine& p)
{
p.print(os);
return os;
}
diff --git a/src/EvtGenModels/EvtBTo4piCP.cpp b/src/EvtGenModels/EvtBTo4piCP.cpp
index bac780c..6831cb0 100644
--- a/src/EvtGenModels/EvtBTo4piCP.cpp
+++ b/src/EvtGenModels/EvtBTo4piCP.cpp
@@ -1,257 +1,258 @@
//--------------------------------------------------------------------------
//
// Environment:
// This software is part of the EvtGen package developed jointly
// for the BaBar and CLEO collaborations. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
// Copyright (C) 1998 Caltech, UCSB
//
// Module: EvtBTo4piCP.cc
//
// Description: Routine to decay B->pi+ pi- pi+ pi-.
//
// Modification history:
//
// RYD March 2, 1997 Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
+#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtCPUtil.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtBTo4piCP.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtConst.hh"
#include <string>
EvtComplex EvtAmpA2(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
//added by Lange Jan4,2000
static EvtId A2M=EvtPDL::getId("a_2-");
static EvtId RHO0=EvtPDL::getId("rho0");
EvtVector4R p4a2,p4rho,p4b;
p4rho=p4pi1+p4pi2;
p4a2=p4rho+p4pi3;
p4b=p4a2+p4pi4;
EvtVector4R p4b_a2,p4rho_a2,p4pi1_a2,p4a2_a2;
p4b_a2=boostTo(p4b,p4a2);
p4rho_a2=boostTo(p4rho,p4a2);
p4pi1_a2=boostTo(p4pi1,p4a2);
p4a2_a2=boostTo(p4a2,p4a2);
EvtVector4R p4pi1_rho;
p4pi1_rho=boostTo(p4pi1_a2,p4rho_a2);
EvtVector4R vb,vrho,vpi,t;
vb=p4b_a2/p4b_a2.d3mag();
vrho=p4rho_a2/p4rho_a2.d3mag();
vpi=p4pi1_rho/p4pi1_rho.d3mag();
t.set(1.0,0.0,0.0,0.0);
// EvtComplex amp_a1,amp_a2;
EvtComplex amp_a2;
// double bwm_a1=EvtPDL::getMeanMass(A1M);
// double gamma_a1=EvtPDL::getWidth(A1M);
double bwm_a2=EvtPDL::getMeanMass(A2M);
double gamma_a2=EvtPDL::getWidth(A2M);
double bwm_rho=EvtPDL::getMeanMass(RHO0);
double gamma_rho=EvtPDL::getWidth(RHO0);
amp_a2=(sqrt(gamma_a2/EvtConst::twoPi)/
((p4a2).mass()-bwm_a2-EvtComplex(0.0,0.5*gamma_a2)))*
(sqrt(gamma_rho/EvtConst::twoPi)/
((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
return amp_a2*
(vb.get(1)*vrho.get(1)+vb.get(2)*vrho.get(2)+vb.get(3)*vrho.get(3))*
(
vpi.get(1)*(vb.get(2)*vrho.get(3)-vb.get(3)*vrho.get(2))+
vpi.get(2)*(vb.get(3)*vrho.get(1)-vb.get(1)*vrho.get(3))+
vpi.get(3)*(vb.get(1)*vrho.get(2)-vb.get(2)*vrho.get(1))
);
}
EvtComplex EvtAmpA1(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
//added by Lange Jan4,2000
static EvtId A1M=EvtPDL::getId("a_1-");
static EvtId RHO0=EvtPDL::getId("rho0");
EvtVector4R p4a1,p4rho,p4b;
p4rho=p4pi1+p4pi2;
p4a1=p4rho+p4pi3;
p4b=p4a1+p4pi4;
EvtVector4R p4b_a1,p4rho_a1,p4pi1_a1,p4a1_a1;
p4b_a1=boostTo(p4b,p4a1);
p4rho_a1=boostTo(p4rho,p4a1);
p4pi1_a1=boostTo(p4pi1,p4a1);
p4a1_a1=boostTo(p4a1,p4a1);
EvtVector4R p4pi1_rho;
p4pi1_rho=boostTo(p4pi1_a1,p4rho_a1);
EvtVector4R vb,vrho,vpi,t;
vb=p4b_a1/p4b_a1.d3mag();
vrho=p4rho_a1/p4rho_a1.d3mag();
vpi=p4pi1_rho/p4pi1_rho.d3mag();
t.set(1.0,0.0,0.0,0.0);
EvtComplex amp_a1;
double bwm_a1=EvtPDL::getMeanMass(A1M);
double gamma_a1=EvtPDL::getWidth(A1M);
// double bwm_a2=EvtPDL::getMeanMass(A2M);
// double gamma_a2=EvtPDL::getWidth(A2M);
double bwm_rho=EvtPDL::getMeanMass(RHO0);
double gamma_rho=EvtPDL::getWidth(RHO0);
amp_a1=(sqrt(gamma_a1/EvtConst::twoPi)/
((p4a1).mass()-bwm_a1-EvtComplex(0.0,0.5*gamma_a1)))*
(sqrt(gamma_rho/EvtConst::twoPi)/
((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
return amp_a1*
(vb.get(1)*vpi.get(1)+vb.get(2)*vpi.get(2)+vb.get(3)*vpi.get(3));
}
std::string EvtBTo4piCP::getName(){
return "BTO4PI_CP";
}
EvtBTo4piCP* EvtBTo4piCP::clone(){
return new EvtBTo4piCP;
}
void EvtBTo4piCP::init(){
// check that there are 18 arguments
checkNArg(18);
checkNDaug(4);
checkSpinParent(EvtSpinType::SCALAR);
checkSpinDaughter(0,EvtSpinType::SCALAR);
checkSpinDaughter(1,EvtSpinType::SCALAR);
checkSpinDaughter(2,EvtSpinType::SCALAR);
checkSpinDaughter(3,EvtSpinType::SCALAR);
}
void EvtBTo4piCP::decay( EvtParticle *p){
//added by Lange Jan4,2000
static EvtId B0=EvtPDL::getId("B0");
static EvtId B0B=EvtPDL::getId("anti-B0");
double t;
EvtId other_b;
EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5);
p->initializePhaseSpace(getNDaug(),getDaugs());
EvtVector4R mom1 = p->getDaug(0)->getP4();
EvtVector4R mom2 = p->getDaug(1)->getP4();
EvtVector4R mom3 = p->getDaug(2)->getP4();
EvtVector4R mom4 = p->getDaug(3)->getP4();
// double alpha=getArg(0);
//double dm=getArg(1);
EvtComplex amp;
EvtComplex A,Abar;
EvtComplex A_a1p,Abar_a1p,A_a2p,Abar_a2p;
EvtComplex A_a1m,Abar_a1m,A_a2m,Abar_a2m;
A_a1p=EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
Abar_a1p=EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
A_a2p=EvtComplex(getArg(6)*cos(getArg(7)),getArg(6)*sin(getArg(7)));
Abar_a2p=EvtComplex(getArg(8)*cos(getArg(9)),getArg(8)*sin(getArg(9)));
A_a1m=EvtComplex(getArg(10)*cos(getArg(11)),getArg(10)*sin(getArg(11)));
Abar_a1m=EvtComplex(getArg(12)*cos(getArg(13)),getArg(12)*sin(getArg(13)));
A_a2m=EvtComplex(getArg(14)*cos(getArg(15)),getArg(14)*sin(getArg(15)));
Abar_a2m=EvtComplex(getArg(16)*cos(getArg(17)),getArg(16)*sin(getArg(17)));
EvtComplex a2p_amp=EvtAmpA2(mom1,mom2,mom3,mom4)+
EvtAmpA2(mom1,mom4,mom3,mom2)+
EvtAmpA2(mom3,mom2,mom1,mom4)+
EvtAmpA2(mom3,mom4,mom1,mom2);
EvtComplex a2m_amp=EvtAmpA2(mom2,mom3,mom4,mom1)+
EvtAmpA2(mom2,mom1,mom4,mom3)+
EvtAmpA2(mom4,mom3,mom2,mom1)+
EvtAmpA2(mom4,mom1,mom2,mom3);
EvtComplex a1p_amp=EvtAmpA1(mom1,mom2,mom3,mom4)+
EvtAmpA1(mom1,mom4,mom3,mom2)+
EvtAmpA1(mom3,mom2,mom1,mom4)+
EvtAmpA1(mom3,mom4,mom1,mom2);
EvtComplex a1m_amp=EvtAmpA1(mom2,mom3,mom4,mom1)+
EvtAmpA1(mom2,mom1,mom4,mom3)+
EvtAmpA1(mom4,mom3,mom2,mom1)+
EvtAmpA1(mom4,mom1,mom2,mom3);
A=A_a2p*a2p_amp+A_a1p*a1p_amp+
A_a2m*a2m_amp+A_a1m*a1m_amp;
Abar=Abar_a2p*a2p_amp+Abar_a1p*a1p_amp+
Abar_a2m*a2m_amp+Abar_a1m*a1m_amp;
if (other_b==B0B){
amp=A*cos(getArg(1)*t/(2*EvtConst::c))+
EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
getArg(2)*EvtComplex(0.0,1.0)*Abar*sin(getArg(1)*t/(2*EvtConst::c));
}
if (other_b==B0){
amp=A*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
getArg(2)*Abar*cos(getArg(1)*t/(2*EvtConst::c));
}
vertex(amp);
return ;
}