Belle II Software development
generators

Modules

 generators data objects
 
 generators modules
 

Classes

class  AAFHInterface
 Class to inferface AAFH/DIAG36 Generator written in Fortran. More...
 
class  BabayagaNLO
 C++ Interface for the Fortran Bhabha and exclusive two photon generator BABAYAGA.NLO. More...
 
class  BBBrem
 Generator for low scattering angle radiative Bhabha events (Beam-Beam Bremsstrahlung). More...
 
class  BHWide
 C++ Interface for the Fortran Bhabha scattering generator BHWide. More...
 
class  SGCosmic
 Class to generate tracks in the cosmics generator and store them in a MCParticle graph. More...
 
class  CRY
 C++ Interface for the generator CRY. More...
 
class  EvtBSemiTauonicDecayRateCalculator
 Class to calculate the differential decay rate, R(D), R(D*), polarizations of tau and D* using BSTD model based on [M. More...
 
class  EvtGenInterface
 Class to interface EvtGen. More...
 
class  EvtGenModelRegister
 Class to keep a register of all Belle2 EvtDecayBases. More...
 
class  EvtB0toKsKK
 Register Decay model EvtB0toKsKK. More...
 
class  EvtBCL
 The class provides the form factors for orbitally excited semileptonic decays. More...
 
class  EvtBCLFF
 BCL Form Factors. More...
 
class  EvtBSemiTauonic
 The EvtGen model of semi-tauonic B decays including new physics effects based on [M. More...
 
class  EvtBSemiTauonic2HDMType2
 The EvtGen model of semi-tauonic B decays including the charged higgs effect of 2HDM Type2 based on [M. More...
 
class  EvtBSemiTauonicAmplitude
 The base class for the calculation of the spin dependent amplitudes for the BSemiTauonic model based on [M. More...
 
class  EvtBSemiTauonicHelicityAmplitudeCalculator
 The class calculates the helicity amplitude of semi-tauonic B decays including new physics effects based on [M. More...
 
class  EvtBSemiTauonicScalarMesonAmplitude
 The class calculates the spin dependent amplitudes of B->Dtaunu decays for the BSemiTauonic model based on [M. More...
 
class  EvtBSemiTauonicVectorMesonAmplitude
 The class calculates the spin dependent amplitudes of B->D*taunu decays for the BSemiTauonic model based on [M. More...
 
class  EvtBtoXsnunu_FERMI
 The evtgen model to produce non-resonant B-> Xs nu nubar decay sample. More...
 
class  EvtEtaFullDalitz
 Class for the simulation of the eta -> pi+pi-pi0 decay with an improved dalitz description. More...
 
class  EvtEtaPi0Dalitz
 Class for the simulation of the eta -> 3pi0 and eta'->3pi0 decays. More...
 
class  EvtEtaPrimeDalitz
 Class for the simulation of the eta' -> pi+ pi- eta and pi0 pi0 eta decays. More...
 
class  EvtFlatDaughter
 The evtgen model to produce flat invariant mass distribution for 2nd and 3rd daughters. More...
 
class  EvtKnunu
 The evtgen model to produce B-> K nu nubar decay sample. More...
 
class  EvtKstarnunu_REV
 The evtgen model to produce B-> Kstar nu nubar decay sample. More...
 
class  EvtPhspCP
 Register Decay model EvtPhspCP. More...
 
class  HepevtReader
 Class to read Hepevt files and store the content in a MCParticle graph. More...
 
class  HepMCReader
 Class to read HepMC files and store the content in a MCParticle graph. More...
 
class  KKGenInterface
 Interface class for using the KKMC generator. More...
 
class  KoralW
 C++ interface for the FORTRAN 4-fermion final state generator KoralW. More...
 
class  LHEReader
 Class to read LHE files and store the content in a MCParticle graph. More...
 
class  ParticleGun
 Class to generate tracks in the particle gun and store them in a MCParticle graph. More...
 
class  Phokhara
 C++ Interface for the Fortran generator PHOKHARA. More...
 
class  ReaderSAD
 Class to read files that have been created by SAD and store their content in a MCParticle graph. More...
 
class  Teegg
 C++ Interface for the Fortran generator TEEGG. More...
 
class  EvtBCLFFTest
 The fixture for testing the EvtBCLFF. More...
 
class  ParticleGunTest
 The fixture for testing the Particlegun. More...
 
class  TouschekReaderTURTLE
 Class to read Touschek files and store their content in a MCParticle graph. More...
 
struct  DualNumber
 Simple structure implementing dual numbers which are used for exact evaluation of the derivatives, see https://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentiation_using_dual_numbers. More...
 
struct  GeneralVector< T >
 3-vector with members of arbitrary type, especially members can be dual numbers More...
 
class  InitialParticleGeneration
 Generate Collision. More...
 

Macros

#define B2_EVTGEN_REGISTER_MODEL(classname)
 Class to register B2_EVTGEN_REGISTER_MODEL.
 

Functions

 B2_EVTGEN_REGISTER_MODEL (EvtB0toK0K0K0)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtB0toKsKK)
 register the model in EvtGen
 
double s13_min (const double &s12, const double &M, const double &m1, const double &m2, const double &m3)
 minimum s13
 
double s13_max (const double &s12, const double &M, const double &m1, const double &m2, const double &m3)
 maximum s13
 
 B2_EVTGEN_REGISTER_MODEL (EvtBCL)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtBSemiTauonic)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtBSemiTauonic2HDMType2)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtBtoKK0K0)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtBtoXsnunu_FERMI)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtEtaFullDalitz)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtEtaPi0Dalitz)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtEtaPrimeDalitz)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtFlatDaughter)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtKnunu)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtKstarnunu_REV)
 register the model in EvtGen
 
 B2_EVTGEN_REGISTER_MODEL (EvtPhspCP)
 register the model in EvtGen
 
 TEST_F (EvtBCLFFTest, Scalar)
 Test calculation of scalar BCL FF.
 
 TEST_F (EvtBCLFFTest, Vector)
 Test calculation of vector BCL FF.
 
 TEST_F (ParticleGunTest, MomentumDistributions)
 Tests momentum generation parameters.
 
 TEST_F (ParticleGunTest, VertexDistributions)
 Tests vertex generation parameters.
 
 TEST_F (ParticleGunTest, AngularDistributions)
 Tests angular generation parameters.
 
double sqrt (double a)
 sqrt for double
 
double tan (double a)
 tan for double
 
double atan (double a)
 atan for double
 
DualNumber operator+ (DualNumber a, DualNumber b)
 addition of dual numbers
 
DualNumber operator+ (DualNumber a, double b)
 addition of dual number and double
 
DualNumber operator- (DualNumber a, double b)
 subtraction of dual number and double
 
DualNumber operator- (double a, DualNumber b)
 subtraction of double and dual number
 
DualNumber operator/ (double a, DualNumber b)
 division of double and dual number
 
DualNumber operator/ (DualNumber a, DualNumber b)
 division of two dual numbers
 
DualNumber operator- (DualNumber a, DualNumber b)
 subtraction of two dual numbers
 
DualNumber operator* (DualNumber a, DualNumber b)
 multiplication of two dual numbers
 
DualNumber operator* (double a, DualNumber b)
 multiplication of double and dual number
 
DualNumber sqrt (DualNumber a)
 sqrt of dual number
 
DualNumber atan (DualNumber a)
 atan of dual number
 
DualNumber tan (DualNumber a)
 tan of dual number
 
template<typename T >
GeneralVector< T > operator+ (GeneralVector< T > a, GeneralVector< T > b)
 addition of two general vectors
 
template<typename T >
dot (GeneralVector< T > a, GeneralVector< T > b)
 dot product of two general vectors
 
template<typename T >
GeneralVector< T > operator* (T a, GeneralVector< T > b)
 product of scalar and general vector
 
template<typename T >
getEcms (GeneralVector< T > p1, GeneralVector< T > p2)
 get CMS energy from HER & LER momentum 3-vectors

 
template<typename T >
GeneralVector< T > getBoost (GeneralVector< T > p1, GeneralVector< T > p2)
 get boost vector from HER & LER momentum 3-vectors

 
template<typename T >
std::pair< T, T > getAnglesCMS (GeneralVector< T > p1, GeneralVector< T > p2)
 get collision axis angles in CM frame from HER & LER momentum 3-vectors

 
template<typename T >
GeneralVector< T > getFourVector (T energy, T angleX, T angleY, bool isHER)
 get 4-momentum from energy and angles of beam
 
Eigen::MatrixXd getGradientMatrix (double eH, double thXH, double thYH, double eL, double thXL, double thYL)
 get gradient matrix of (eH, thXH, thYH, eL, thXL, thYL) -> (eCMS, boostX, boostY, boostY, angleXZ, angleYZ) transformation
 
Eigen::VectorXd getCentralValues (double eH, double thXH, double thYH, double eL, double thXL, double thYL)
 get vector including (Ecms, boostX, boostY, boostZ, angleXZ, angleYZ) from beam energies and angles
 
Eigen::MatrixXd transformCov (const Eigen::MatrixXd &covHER, const Eigen::MatrixXd &covLER, const Eigen::MatrixXd &grad)
 transform covariance matrix to the variables (Ecms, boostX, boostY, boostZ, angleXZ, angleYZ)
 
double getTotalEnergy (const std::vector< ROOT::Math::PxPyPzMVector > &ps)
 Main function scaleParticleEnergies in this header scales momenta in CMS of all particles in the final state by a constant factor such that the overall collision energy is slightly changed.
 
double getScaleFactor (const std::vector< ROOT::Math::PxPyPzMVector > &ps, double EcmsTarget)
 Find approximative value of the momentum scaling factor which makes the total energy of the particles provided in the vector ps equal to EcmsTarget.
 
std::vector< ROOT::Math::PxPyPzMVector > scaleMomenta (const std::vector< ROOT::Math::PxPyPzMVector > &ps, double C)
 Scale momenta of the particles in the vector ps with a factor C.
 
void scaleParticleEnergies (MCParticleGraph &mpg, double EcmsTarget)
 Scale momenta of the particles by a constant factor such that total energy is changed to EcmsTarget.
 
static Eigen::MatrixXd toEigen (TMatrixDSym m)
 transform matrix from ROOT to Eigen format
 
void init ()
 Initializes the generator.
 
void generateEvent (MCParticleGraph &mcGraph)
 Generates one single event.
 
void term ()
 Terminates the generator.
 
std::string getName ()
 Get function Name

 
EvtDecayBase * clone ()
 Clone the decay of B0toKsKK.
 
void init ()
 Initialize standard stream objects

 
void initProbMax ()
 Initialize standard stream objects for probability function

 
void decay (EvtParticle *p)
 Member of particle in EvtGen.
 
EvtVector4R umu (const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
 Function 4Vector umu.
 
EvtVector4R Smu (const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
 Function 4Vector Smu.
 
EvtVector4R Lmu (const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
 Function 4Vector Lmu.
 
EvtTensor4C gmunu_tilde (const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
 Function Tensor gmunu

 
EvtTensor4C Tmunu (const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
 Function Tensor Tmunu

 
EvtTensor4C Multiply (const EvtTensor4C &t1, const EvtTensor4C &t2)
 Function Tensor Multiply

 
EvtTensor4C RaiseIndices (const EvtTensor4C &t)
 Function RaiseIndices

 
void RaiseIndex (EvtVector4R &vector)
 Member function RaiseIndices.
 
EvtTensor4C Mmunu (const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
 Function Tensor Mmunu.
 
double BWBF (const double &q, const unsigned int &L)
 Meson radius

 
double BWBF (const double &q, const double &q0, const unsigned int &L)
 Meson radius

 
EvtComplex BreitWigner (const double &m, const double &m0, const double &Gamma0, const double &q, const double &q0, const unsigned int &L)
 BreitWigner Shape.
 
EvtVector4R Boost (const EvtVector4R &p4, const EvtVector4R &boost)
 Parameter for boost frame

 
double p (const double &mab, const double &M, const double &mc)
 Constants p

 
double q (const double &mab, const double &ma, const double &mb)
 Constants q.
 
EvtComplex Flatte_k (const double &s, const double &m_h)
 Constant Flatte_k.
 
EvtComplex Flatte (const double &m, const double &m0)
 Constant Flatte.
 
EvtComplex A_f0ks (const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
 A_f0ks is amplitude of f0.
 
EvtComplex A_phiks (const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
 A_phiks is amplitude of phi.
 
EvtComplex A_fxks (const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
 A_fxks is amplitude of fxks.
 
EvtComplex A_chic0ks (const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
 A_chic0ks is amplitude of chic0ks.
 
EvtComplex A_kknr (const EvtVector4R &p4k1, const EvtVector4R &p4k2, const double &alpha_kk)
 A_kknr is amplitude of kknr.
 
 EvtBCL ()
 Default constructor.
 
virtual ~EvtBCL ()
 virtual destructor
 
std::string getName ()
 Returns name of module.
 
EvtDecayBase * clone ()
 Clones module.
 
void decay (EvtParticle *p)
 Creates a decay.
 
void initProbMax ()
 Sets maximal probab.
 
void init ()
 Initializes module.
 
 EvtBCLFF (int numarg, double *arglist)
 constructor
 
void getscalarff (EvtId parent, EvtId daughter, double t, double, double *fpf, double *f0f)
 Scalar FF's.
 
void getvectorff (EvtId parent, EvtId daughter, double t, double, double *a1f, double *a2f, double *vf, double *a0f)
 Vector FF's.
 
void gettensorff (EvtId parent, EvtId daughter, double t, double, double *hf, double *kf, double *bp, double *bm)
 Not Implemented.
 
void getbaryonff (EvtId, EvtId, double, double, double *, double *, double *, double *)
 Not Implemented.
 
void getdiracff (EvtId, EvtId, double, double, double *, double *, double *, double *, double *, double *)
 Not Implemented.
 
void getraritaff (EvtId, EvtId, double, double, double *, double *, double *, double *, double *, double *, double *, double *)
 Not Implemented.
 
 EvtBSemiTauonic ()
 The default constructor

 
virtual ~EvtBSemiTauonic ()
 The destructor

 
std::string getName ()
 The function returns the model name.
 
EvtDecayBase * clone ()
 The function makes a copy of an EvtBSTD object.
 
void decay (EvtParticle *p)
 The function evaluates the decay amplitude of the parent particle.
 
void initProbMax ()
 The function sets the maximum value of the probability.
 
void init ()
 The function initializes the decay.
 
 EvtBSemiTauonic2HDMType2 ()
 The default constructor

 
virtual ~EvtBSemiTauonic2HDMType2 ()
 The destructor

 
std::string getName ()
 The function returns the model name.
 
EvtDecayBase * clone ()
 The function makes a copy of an EvtBSemiTauonic2HDMType2 object.
 
void decay (EvtParticle *p)
 The function evaluates the decay amplitude of the parent particle.
 
void initProbMax ()
 The function sets the maximum value of the probability.
 
void init ()
 The function initializes the decay.
 
EvtSpinDensity RotateToHelicityBasisInBoostedFrame (const EvtParticle *p, EvtVector4R p4boost)
 The function calculates the rotation matrix to convert the spin basis to the helicity basis in the boosted frame.
 
double CalcMaxProb (EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator)
 The function calculates the maximum probability.
 
 EvtBSemiTauonicHelicityAmplitudeCalculator ()
 The default constructor.
 
 EvtBSemiTauonicHelicityAmplitudeCalculator (const double rho12, const double rhoA12, const double ffR11, const double ffR21, const double AS1, const double AR3, const double bottomMass, const double charmMass, const EvtComplex &CV1, const EvtComplex &CV2, const EvtComplex &CS1, const EvtComplex &CS2, const EvtComplex &CT, const double parentMass, const double DMass, const double DstarMass)
 The constructor with HQET form factor parameters, Wilson coefficients of new physics contributions and parent B, daughter D(*) meson masses.
 
EvtComplex helAmp (double mtau, int tauhel, int Dhel, double w, double costau) const
 The function calculates the helicity amplitude.
 
EvtComplex helAmp (const EvtComplex &CV1, const EvtComplex &CV2, const EvtComplex &CS1, const EvtComplex &CS2, const EvtComplex &CT, double mtau, int tauhel, int Dhel, double w, double costau) const
 The function calculates helicity amplitudes with given Wilson coefficients.
 
double Lep (const double mtau, int tauhel, int whel, double q2, double costau) const
 The function to calculate the Leptonic Amplitudes for B->D*taunu decay of the vector type contribution.
 
double Lep (const double mtau, int tauhel, double q2, double costau) const
 The function to calculate the Leptonic Amplitudes for B->Dtaunu decay of the scalar type contribution.
 
double Lep (const double mtau, int tauhel, int whel1, int whel2, double q2, double costau) const
 The function to calculate the Leptonic Amplitudes for B->D*taunu decay of the tensor type contribution.
 
double HadV1 (int Dhel, int whel, double w) const
 The function to calculate the Hadronic Amplitudes of left handed (V-A) type contribution.
 
double HadV2 (int Dhel, int whel, double w) const
 The function to calculate the Hadronic Amplitudes of right handed (V+A) type contribution.
 
double HadS1 (int Dhel, double w) const
 The function to calculate the Hadronic Amplitudes of scalar (S+P) type contribution.
 
double HadS2 (int Dhel, double w) const
 The function to calculate the Hadronic Amplitudes of scalar (S-P) type contribution.
 
double HadT (int Dhel, int whel1, int whel2, double w) const
 The function to calculate the Hadronic Amplitudes of tensor type contribution.
 
double helampSM (double mtau, int tauhel, int Dhel, double w, double costau) const
 Helicity Amplitudes of SM (left handed) contribution.
 
double helampV1 (double mtau, int tauhel, int Dhel, double w, double costau) const
 Helicity Amplitudes of left handed (V-A) contribution.
 
double helampV2 (double mtau, int tauhel, int Dhel, double w, double costau) const
 Helicity Amplitudes of right handed (V+A) contribution.
 
double helampS1 (double mtau, int tauhel, int Dhel, double w, double costau) const
 Helicity Amplitudes of scalar (S+P) type contribution.
 
double helampS2 (double mtau, int tauhel, int Dhel, double w, double costau) const
 Helicity Amplitudes of scalar (S-P) type contribution.
 
double helampT (double mtau, int tauhel, int Dhel, double w, double costau) const
 Helicity Amplitudes of tensor type contribution.
 
double hp (double w) const
 HQET D vector form factor h_+(w).
 
double hm (double w) const
 HQET D vector form factor h_-(w).
 
double hA1 (double w) const
 HQET D* axial vector form factor h_{A1}(w).
 
double hV (double w) const
 HQET D* axial vector form factor h_V(w).
 
double hA2 (double w) const
 HQET D* axial vector form factor h_{A2}(w).
 
double hA3 (double w) const
 HQET D* axial vector form factor h_{A3}(w).
 
double hS (double w) const
 D scalar form factor h_S(w) in terms of vector form factors.
 
double hP (double w) const
 D* pseudo scalar form factor h_P(w) in terms of axial vector form factors.
 
double hT (double w) const
 D tensor form factor h_T(w) in terms of vector form factors.
 
double hT1 (double w) const
 D* tensor form factor h_{T1}(w) in terms of axial vector form factors.
 
double hT2 (double w) const
 D* tensor form factor h_{T2}(w).
 
double hT3 (double w) const
 D* tensor form factor h_{T3}(w).
 
double z (double w) const
 CLN form factor z.
 
double ffV1 (double w) const
 CLN form factor V1.
 
double ffS1 (double w) const
 CLN form factor S1.
 
double ffA1 (double w) const
 CLN form factor A1.
 
double ffR1 (double w) const
 CLN form factor R1.
 
double ffR2 (double w) const
 CLN form factor R2.
 
double ffR3 (double w) const
 CLN form factor R3.
 
double dS1 (double w) const
 HQET correction factor for the scalar form factor for B->Dtaunu.
 
double dR3 (double w) const
 HQET correction factor for the scalar form factor for B->D*taunu.
 
double mD (int Dhel) const
 Daughter D(*) meson mass.
 
double v (double mtau, double q2) const
 Function to calculate the tau velocity.
 
double qh2 (int Dhel, double w) const
 Function to calculate the q^2 divided by the square of parent mass (m_B^2).
 
double q2 (int Dhel, double w) const
 Function to calculate the q^2 of the decay (square of l+nu invariant mass).
 
bool chkDhel (int Dhel) const
 sanity checkers
 
bool chkwhel (int whel) const
 Function to check if whel is in the valid range.
 
bool chktauhel (int tauhel) const
 Function to check if tauhel is in the valid range.
 
void CalcAmp (EvtParticle *parent, EvtAmp &amp, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator) override
 The function calculates the spin dependent amplitude.
 
void CalcAmp (EvtParticle *parent, EvtAmp &amp, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator) override
 The function calculates the spin dependent amplitude.
 
virtual ~EvtBtoXsnunu_FERMI ()
 Destructor.
 
std::string getName ()
 The function which returns the name of the model.
 
EvtDecayBase * clone ()
 The function which makes a copy of the model.
 
void decay (EvtParticle *p)
 The function to determine kinematics of daughter particles based on dGdsb distribution.
 
void initProbMax ()
 The function to sets a maximum probability.
 
void init ()
 The function for an initialization.
 
double dGdsbProb (double _sb)
 The function returns the probability density value depending on sb.
 
double FermiMomentum (double pf)
 The function returns a momentum of b quark.
 
double FermiMomentumProb (double pb, double pf)
 The function returns a probability based on the Fermi motion model.
 
virtual ~EvtEtaFullDalitz ()
 Default Destructor.
 
std::string getName ()
 Returns the model name: ETA_FULLDALITZ.
 
EvtDecayBase * clone ()
 Makes a copy of the pointer to the class.
 
void init ()
 Checks that the number of input parameters are correct:
 
void initProbMax ()
 Sets the Maximum probability for the PHSP reweight.
 
void decay (EvtParticle *p)
 Function that implements the energy-dependent Dalitz.
 
virtual ~EvtEtaPi0Dalitz ()
 Default destructor.
 
std::string getName ()
 Returns the name of the model: ETA_PI0DALITZ.
 
EvtDecayBase * clone ()
 Makes a copy of the class object.
 
void init ()
 Checks that the number of input parameters are correct:
 
void initProbMax ()
 Sets the Maximum probability for the PHSP reweight.
 
void decay (EvtParticle *p)
 Function that implements the energy-dependent Dalitz.
 
virtual ~EvtEtaPrimeDalitz ()
 Default destructor.
 
std::string getName ()
 Returns the model name: ETAPRIME_DALITZ.
 
EvtDecayBase * clone ()
 Returns a copy of the class object.
 
void init ()
 Checks that the number of input parameters are correct:
 
void initProbMax ()
 Sets the Maximum probability for the PHSP reweight.
 
void decay (EvtParticle *p)
 Function that implements the energy-dependent Dalitz.
 
virtual ~EvtFlatDaughter ()
 Destructor.
 
std::string getName ()
 The function which returns the name of the model.
 
EvtDecayBase * clone ()
 The function which makes a copy of the model.
 
void decay (EvtParticle *p)
 The function to determine kinematics of daughter particles.
 
void initProbMax ()
 The function to sets a maximum probability.
 
void init ()
 The function for an initialization.
 
virtual ~EvtKnunu ()
 Destructor.
 
std::string getName ()
 The function which returns the name of the model.
 
EvtDecayBase * clone ()
 The function which makes a copy of the model.
 
void decay (EvtParticle *p)
 The function to calculate a quark decay amplitude.
 
void init ()
 The function for an initialization.
 
void initProbMax ()
 The function to sets a maximum probability.
 
virtual ~EvtKstarnunu_REV ()
 Destructor.
 
std::string getName ()
 The function which returns the name of the model.
 
EvtDecayBase * clone ()
 The function which makes a copy of the model.
 
void decay (EvtParticle *p)
 The function to calculate a quark decay amplitude.
 
void init ()
 The function for an initialization.
 
void initProbMax ()
 The function to sets a maximum probability.
 
std::string getName ()
 Get function Name

 
EvtDecayBase * clone ()
 Clone the decay.
 
void init ()
 Initialize standard stream objects

 
void initProbMax ()
 Initialize standard stream objects for probability function

 
void decay (EvtParticle *p)
 Member of particle in EvtGen.
 
double dGammadwdcostau (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double w, double costau)
 Function calculates the differential decay rate dGamma/dw/dcostau.
 
double dGammadw (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double w)
 Function calculates the differential decay rate dGamma/dw, integrated for costau.
 
double dGammadcostau (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double costau)
 Function calculates the differential decay rate dGamma/dcostau, integrated for w.
 
double Gamma (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel)
 Function calculates the helicity dependent decay rate Gamma, integrated for w and costau.
 
double GammaD (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau)
 Function calculates the decay rate Gamma for Dtaunu decay, integrated for w and costau and summed for helicities.
 
double GammaDstar (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau)
 Function calculates the differential decay rate Gamma for D*taunu decay, integrated for w and costau and summed for helicities.
 
double GammaSMD (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mlep=0.0005110)
 Function calculates the SM decay rate Gamma for Dlnu decay, integrated for w and costau and summed for helicities.
 
double GammaSMDstar (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mlep=0.0005110)
 Function calculates the SM decay rate Gamma for D*lnu decay, integrated for w and costau and summed for helicities.
 
double RGammaD (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau, const double mlep=0.0005110)
 Function calculates the ratio of Br(B->Dtaunu)/Br(B->Dlnu), R(D).
 
double RGammaDstar (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau, const double mlep=0.0005110)
 Function calculates the ratio of Br(B->Dtaunu)/Br(B->Dlnu), R(D*).
 
double PtauD (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau)
 Function calculates the polarization of tau, (RH - LH)/(LH + RH), in B->Dtaunu decay.
 
double PtauDstar (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau)
 Function calculates the polarization of tau, (RH - LH)/(LH + RH), in B->D*taunu decay.
 
double PDstar (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mtau)
 Function calculates the polarization of D*, longitudinal/(longitudinal + transverse), in B->D*taunu decay.
 
double pf (const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int Dhel, double w)
 Phase space factor, which is multiplied to the helicity amplitude to calculate the decay rate.
 
double EvaluateByW (double *x, double *param)
 Function used internally for numerical integration.
 
double EvaluateByCostau (double *x, double *param)
 Function used internally for numerical integration.
 
double EvaluateBy2D (double *x, double *param)
 Function used internally for numerical integration.
 
static EvtGenModelRegistergetInstance ()
 Return reference to the instance.
 
static std::list< EvtDecayBase * > getModels ()
 Return a list of models.
 
void checkVariable (const std::string &name, const std::map< int, bool > &allowed, ParticleGun::EDistribution &dist, std::vector< double > &pars)
 check one of the variables given a list of allowed and excluded distributions
 
void initialize ()
 function to be executed on initialize()
 
ROOT::Math::XYZVector generateVertex (const ROOT::Math::XYZVector &initial, const TMatrixDSym &cov, MultivariateNormalGenerator &gen) const
 generate the vertex
 
TMatrixDSym adjustCovMatrix (TMatrixDSym cov) const
 adjust smearing covariance matrix based on the generation flags
 
ConditionalGaussGenerator initConditionalGenerator (const ROOT::Math::PxPyPzEVector &pHER, const ROOT::Math::PxPyPzEVector &pLER, const TMatrixDSym &covHER, const TMatrixDSym &covLER)
 Initialize the conditional generator using HER & LER 4-vectors and HER & LER covariance matrices describing spread.
 
ROOT::Math::PxPyPzEVector generateBeam (const ROOT::Math::PxPyPzEVector &initial, const TMatrixDSym &cov, MultivariateNormalGenerator &gen) const
 generate 4 vector for one beam
 
MCInitialParticlesgenerate ()
 Generate a new event.
 
MCInitialParticlesgenerate (int allowedFlags)
 Generate a new event wit a particular set of allowed flags.
 
ROOT::Math::XYZVector getVertexConditional ()
 Generate vertex position and possibly update the generator of Lorentz transformation.
 
ROOT::Math::XYZVector updateVertex (bool force=false)
 Update the vertex position:
 

Variables

static const double me = Const::electron.getMass()
 electron mass
 
static EvtGen * s_evtgen = nullptr
 pointer to the evtgen instance
 

Detailed Description

Macro Definition Documentation

◆ B2_EVTGEN_REGISTER_MODEL

#define B2_EVTGEN_REGISTER_MODEL (   classname)
Value:
namespace {\
Belle2::EvtGenModelRegister::Factory<classname> EvtGenModelFactory_##classname; \
}
Helper Class to easily register new EvtDecayBases for Belle2.

Class to register B2_EVTGEN_REGISTER_MODEL.

Definition at line 75 of file EvtGenModelRegister.h.

Function Documentation

◆ A_chic0ks()

EvtComplex A_chic0ks ( const EvtVector4R &  p4ks,
const EvtVector4R &  p4kp,
const EvtVector4R &  p4km 
)

A_chic0ks is amplitude of chic0ks.

Definition at line 800 of file EvtB0toKsKK.cc.

803 {
804 static EvtId chic0 = EvtPDL::getId("chi_c0");
805
806 const double chic0_m = EvtPDL::getMeanMass(chic0);
807 const double chic0_w = EvtPDL::getWidth(chic0);
808
809 const EvtVector4R p4kpkm = p4kp + p4km;
810 const double mkpkm = p4kpkm.mass();
811
812 //Angular distribution
813 const EvtComplex H_chic0ks = 1.0;
814
815 //Barrier factors
816 const EvtComplex D_chic0ks = 1.0;
817
818 //Line shape
819 const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
820 const double q0_chic0 = q(chic0_m, p4kp.mass(), p4km.mass());
821 const EvtComplex L_chic0ks =
822 BreitWigner(mkpkm, chic0_m, chic0_w, p4kp_kpkm.d3mag(), q0_chic0, 0);
823
824 //Amplitude
825 const EvtComplex A_chic0ks = D_chic0ks * H_chic0ks * L_chic0ks;
826
827 return A_chic0ks;
828 }
double q(const double &mab, const double &ma, const double &mb)
Constants q.
Definition: EvtB0toKsKK.cc:594
EvtComplex BreitWigner(const double &m, const double &m0, const double &Gamma0, const double &q, const double &q0, const unsigned int &L)
BreitWigner Shape.
Definition: EvtB0toKsKK.cc:521
EvtVector4R Boost(const EvtVector4R &p4, const EvtVector4R &boost)
Parameter for boost frame
Definition: EvtB0toKsKK.cc:538
EvtComplex A_chic0ks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_chic0ks is amplitude of chic0ks.
Definition: EvtB0toKsKK.cc:800

◆ A_f0ks()

EvtComplex A_f0ks ( const EvtVector4R &  p4ks,
const EvtVector4R &  p4kp,
const EvtVector4R &  p4km 
)

A_f0ks is amplitude of f0.

Definition at line 639 of file EvtB0toKsKK.cc.

642 {
643 static EvtId f0 = EvtPDL::getId("f_0");
644
645 const double f0_m = EvtPDL::getMeanMass(f0);
646
647 const EvtVector4R p4kpkm = p4kp + p4km;
648 const double mkpkm = p4kpkm.mass();
649
650 //Angular distribution
651 const EvtComplex H_f0ks = 1.0;
652
653 //Barrier factors
654 const EvtComplex D_f0ks = 1.0;
655
656 //Line shape
657 const EvtComplex L_f0ks = Flatte(mkpkm, f0_m);
658
659 //Amplitude
660 const EvtComplex A_f0ks = D_f0ks * H_f0ks * L_f0ks;
661
662 return A_f0ks;
663 }
EvtComplex A_f0ks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_f0ks is amplitude of f0.
Definition: EvtB0toKsKK.cc:639
EvtComplex Flatte(const double &m, const double &m0)
Constant Flatte.
Definition: EvtB0toKsKK.cc:617

◆ A_fxks()

EvtComplex A_fxks ( const EvtVector4R &  p4ks,
const EvtVector4R &  p4kp,
const EvtVector4R &  p4km 
)

A_fxks is amplitude of fxks.

Definition at line 770 of file EvtB0toKsKK.cc.

773 {
774 static EvtId fx = EvtPDL::getId("f_0(1500)");
775
776 const double fx_m = EvtPDL::getMeanMass(fx);
777 const double fx_w = EvtPDL::getWidth(fx);
778
779 const EvtVector4R p4kpkm = p4kp + p4km;
780 const double mkpkm = p4kpkm.mass();
781
782 //Angular distribution
783 const EvtComplex H_fxks = 1.0;
784
785 //Barrier factors
786 const EvtComplex D_fxks = 1.0;
787
788 //Line shape
789 const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
790 const double q0_fx = q(fx_m, p4kp.mass(), p4km.mass());
791 const EvtComplex L_fxks =
792 BreitWigner(mkpkm, fx_m, fx_w, p4kp_kpkm.d3mag(), q0_fx, 0);
793
794 //Amplitude
795 const EvtComplex A_fxks = D_fxks * H_fxks * L_fxks;
796
797 return A_fxks;
798 }
EvtComplex A_fxks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_fxks is amplitude of fxks.
Definition: EvtB0toKsKK.cc:770

◆ A_kknr()

EvtComplex A_kknr ( const EvtVector4R &  p4k1,
const EvtVector4R &  p4k2,
const double &  alpha_kk 
)

A_kknr is amplitude of kknr.

Definition at line 830 of file EvtB0toKsKK.cc.

833 {
834 const EvtVector4R p4kk = p4k1 + p4k2;
835 const double m2kk = p4kk.mass2();
836
837 const EvtComplex A_kknr = exp(-alpha_kk * m2kk);
838
839 return A_kknr;
840 }
EvtComplex A_kknr(const EvtVector4R &p4k1, const EvtVector4R &p4k2, const double &alpha_kk)
A_kknr is amplitude of kknr.
Definition: EvtB0toKsKK.cc:830

◆ A_phiks()

EvtComplex A_phiks ( const EvtVector4R &  p4ks,
const EvtVector4R &  p4kp,
const EvtVector4R &  p4km 
)

A_phiks is amplitude of phi.

Definition at line 695 of file EvtB0toKsKK.cc.

698 {
699 static EvtId phi = EvtPDL::getId("phi");
700
701 const double phi_m = EvtPDL::getMeanMass(phi);
702 const double phi_w = EvtPDL::getWidth(phi);
703
704 const EvtVector4R p4kpkm = p4kp + p4km;
705 const double mkpkm = p4kpkm.mass();
706
707 //Angular distribution
708 const EvtVector4R S_mu = Smu(p4kp, p4km, p4ks);
709 const EvtVector4R L_mu = Lmu(p4kp, p4km, p4ks);
710
711 const EvtTensor4C g_munu_tilde = gmunu_tilde(p4kp, p4km, p4ks);
712 const EvtTensor4C g__munu_tilde = RaiseIndices(g_munu_tilde);
713
714 EvtComplex H_phiks = 0.0;
715 for (unsigned int i = 0; i < 4; i++)
716 for (unsigned int j = 0; j < 4; j++)
717 H_phiks += g__munu_tilde.get(i, j) * S_mu.get(i) * L_mu.get(j);
718
719 //Barrier factors
720 const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
721
722 const EvtComplex D_phiks = BWBF(p4kp_kpkm.d3mag(), 1);
723
724 //Line shape
725 const double q0_phi = q(phi_m, p4kp.mass(), p4km.mass());
726 const EvtComplex L_phiks =
727 BreitWigner(mkpkm, phi_m, phi_w, p4kp_kpkm.d3mag(), q0_phi, 1);
728
729 //Amplitude
730 const EvtComplex A_phiks = D_phiks * H_phiks * L_phiks;
731
733 /*const EvtVector4R p4ks_kpkm = Boost(p4ks, p4kpkm);
734 const double p4ks_kpkm_alt = p(mkpkm, (p4ks+p4kp+p4km).mass(), p4ks.mass());
735 const double p4kp_kpkm_alt = q(mkpkm, p4kp.mass(), p4km.mass());
736 //std::cout << p4ks_kpkm.d3mag() << std::endl;
737 //std::cout << p4ks_kpkm_alt << std::endl;
738 //std::cout << p4kp_kpkm.d3mag() << std::endl;
739 //std::cout << p4kp_kpkm_alt << std::endl;
740 const double costheta =
741 p4ks_kpkm.dot(p4kp_kpkm)/p4ks_kpkm.d3mag()/p4kp_kpkm.d3mag();
742 //const double z = p4ks_kpkm.d3mag()/(p4ks+p4kp+p4km).mass();
743 //const double H_phiks_alt = -sqrt(1.0+(z*z))*costheta;
744 const double sp_min =
745 s13_min( mkpkm*mkpkm, (p4ks+p4kp+p4km).mass(),
746 p4kp.mass(), p4km.mass(), p4ks.mass() );
747 const double sp_max =
748 s13_max( mkpkm*mkpkm, (p4ks+p4kp+p4km).mass(),
749 p4kp.mass(), p4km.mass(), p4ks.mass() );
750 const double sp = (p4ks+p4kp).mass2();
751 const double costheta_alt = (sp_max + sp_min - (2.0*sp))/(sp_max - sp_min);
752 //std::cout << costheta << std::endl;
753 //std::cout << costheta_alt << std::endl;
754 //exit(1);
755 const double z_alt = p4ks_kpkm_alt/(p4ks+p4kp+p4km).mass();
756 const double H_phiks_alt = -sqrt(1.0+(z_alt*z_alt))*costheta_alt;
757 //std::cout << H_phiks << std::endl;
758 //std::cout << H_phiks_alt << std::endl;
759 if( real(H_phiks) != H_phiks_alt )
760 {
761 std::cout << H_phiks << std::endl;
762 std::cout << H_phiks_alt << std::endl;
763 }
764 const EvtComplex A_phiks = D_phiks*H_phiks_alt*L_phiks;*/
766
767 return A_phiks;
768 }
EvtVector4R Smu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector Smu.
Definition: EvtB0toKsKK.cc:370
EvtComplex A_phiks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_phiks is amplitude of phi.
Definition: EvtB0toKsKK.cc:695
EvtTensor4C gmunu_tilde(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor gmunu
Definition: EvtB0toKsKK.cc:409
EvtTensor4C RaiseIndices(const EvtTensor4C &t)
Function RaiseIndices
Definition: EvtB0toKsKK.cc:453
EvtVector4R Lmu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector Lmu.
Definition: EvtB0toKsKK.cc:390
double BWBF(const double &q, const unsigned int &L)
Meson radius
Definition: EvtB0toKsKK.cc:479

◆ adjustCovMatrix()

TMatrixDSym adjustCovMatrix ( TMatrixDSym  cov) const
private

adjust smearing covariance matrix based on the generation flags

Definition at line 38 of file InitialParticleGeneration.cc.

39 {
40 // no smearing
41 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearBeam))
42 return TMatrixDSym(3);
43
44 // don't smear energy
45 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearBeamEnergy)) {
46 for (int i = 0; i < 3; ++i)
47 cov(i, 0) = cov(0, i) = 0;
48 }
49
50 // don't smear angles
51 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearBeamDirection)) {
52 for (int i = 1; i < 3; ++i)
53 for (int j = 1; j < 3; ++j)
54 cov(i, j) = cov(j, i) = 0;
55 }
56
57 return cov;
58
59 }
DBObjPtr< BeamParameters > m_beamParams
Datastore object containing the nominal beam parameters.
@ c_smearBeam
smear the full beam momentum (energy and direction)
@ c_smearBeamEnergy
smear energy of HER and LER (but not direction)
@ c_smearBeamDirection
smear direction of HER and LER (but not energy)

◆ atan() [1/2]

double atan ( double  a)
inline

atan for double

Definition at line 34 of file beamHelpers.h.

34{ return std::atan(a); }

◆ atan() [2/2]

DualNumber atan ( DualNumber  a)
inline

atan of dual number

Definition at line 121 of file beamHelpers.h.

122 {
123 return DualNumber(std::atan(a.x), 1. / (1 + a.x * a.x) * a.dx);
124 }

◆ Boost()

EvtVector4R Boost ( const EvtVector4R &  p4,
const EvtVector4R &  boost 
)

Parameter for boost frame

Definition at line 538 of file EvtB0toKsKK.cc.

540 {
541 const double bx = boost.get(1) / boost.get(0);
542 const double by = boost.get(2) / boost.get(0);
543 const double bz = boost.get(3) / boost.get(0);
544
545 const double bx2 = bx * bx;
546 const double by2 = by * by;
547 const double bz2 = bz * bz;
548
549 const double b2 = bx2 + by2 + bz2;
550 if (b2 == 0.0)
551 return p4;
552 assert(b2 < 1.0);
553
554 const double gamma = 1.0 / sqrt(1 - b2);
555
556 const double gb2 = (gamma - 1.0) / b2;
557
558 const double gb2xy = gb2 * bx * by;
559 const double gb2xz = gb2 * bx * bz;
560 const double gb2yz = gb2 * by * bz;
561
562 const double gbx = gamma * bx;
563 const double gby = gamma * by;
564 const double gbz = gamma * bz;
565
566 const double e_b =
567 (gamma * p4.get(0)) - (gbx * p4.get(1)) - (gby * p4.get(2)) - (gbz * p4.get(3));
568 const double x_b =
569 (-gbx * p4.get(0)) + ((1.0 + (gb2 * bx2)) * p4.get(1)) +
570 (gb2xy * p4.get(2)) + (gb2xz * p4.get(3));
571 const double y_b =
572 (-gby * p4.get(0)) + (gb2xy * p4.get(1)) +
573 ((1.0 + (gb2 * by2)) * p4.get(2)) + (gb2yz * p4.get(3));
574 const double z_b =
575 (-gbz * p4.get(0)) + (gb2xz * p4.get(1)) +
576 (gb2yz * p4.get(2)) + ((1.0 + (gb2 * bz2)) * p4.get(3));
577
578 const EvtVector4R p4_b(e_b, x_b, y_b, z_b);
579
580 return p4_b;
581 }

◆ BreitWigner()

EvtComplex BreitWigner ( const double &  m,
const double &  m0,
const double &  Gamma0,
const double &  q,
const double &  q0,
const unsigned int &  L 
)

BreitWigner Shape.

Definition at line 521 of file EvtB0toKsKK.cc.

525 {
526 const double s = m * m;
527 const double s0 = m0 * m0;
528
529 const double Gamma =
530 Gamma0 * m0 / m * pow(q / q0, (2 * L) + 1) *
531 BWBF(q, q0, L) * BWBF(q, q0, L);
532
533 const EvtComplex BreitWigner = 1.0 / EvtComplex(s0 - s, -m0 * Gamma);
534
535 return BreitWigner;
536 }

◆ BWBF() [1/2]

double BWBF ( const double &  q,
const double &  q0,
const unsigned int &  L 
)

Meson radius

Definition at line 499 of file EvtB0toKsKK.cc.

501 {
502 //Meson radius
503 const double d = 1.0;
504
505 const double z = q * q * d * d;
506 const double z0 = q0 * q0 * d * d;
507
508 double bwbf = 1.0;
509
510 if (L == 1)
511 bwbf = sqrt((1.0 + z0) / (1.0 + z));
512 if (L == 2)
513 bwbf = sqrt((((z0 - 3.0) * (z0 - 3.0)) + (9.0 * z0)) / (((z - 3.0) * (z - 3.0)) + (9.0 * z)));
514 if (L > 2) {
515 std::cout << "ERROR: BWBF not implemented for L>2" << std::endl;
516 exit(1);
517 }
518 return bwbf;
519 }

◆ BWBF() [2/2]

double BWBF ( const double &  q,
const unsigned int &  L 
)

Meson radius

Definition at line 479 of file EvtB0toKsKK.cc.

480 {
481 //Meson radius
482 const double d = 1.0;
483
484 const double z = q * q * d * d;
485
486 double bwbf = 1.0;
487
488 if (L == 1)
489 bwbf = sqrt(2.0 * z / (1.0 + z));
490 if (L == 2)
491 bwbf = sqrt(13.0 * z * z / (((z - 3.0) * (z - 3.0)) + (9.0 * z)));
492 if (L > 2) {
493 std::cout << "ERROR: BWBF not implemented for L>2" << std::endl;
494 exit(1);
495 }
496 return bwbf;
497 }

◆ CalcAmp() [1/2]

void CalcAmp ( EvtParticle *  parent,
EvtAmp &  amp,
EvtBSemiTauonicHelicityAmplitudeCalculator HelicityAmplitudeCalculator 
)
overridevirtual

The function calculates the spin dependent amplitude.

Parameters
parenta pointer to the parent particle.
ampa pointer to fill the calculated spin dependent amplitude.
HelicityAmplitudeCalculatora pointer to the calculator of the helicity dependent amplitude. The function calculate the spin dependent amplitude of the semi-tauonic decay to a scalar meson (D meson).

Implements EvtBSemiTauonicAmplitude.

Definition at line 25 of file EvtBSemiTauonicScalarMesonAmplitude.cc.

28 {
29
30 static EvtId EM = EvtPDL::getId("e-");
31 static EvtId MUM = EvtPDL::getId("mu-");
32 static EvtId TAUM = EvtPDL::getId("tau-");
33// static EvtId EP = EvtPDL::getId("e+");
34// static EvtId MUP = EvtPDL::getId("mu+");
35// static EvtId TAUP = EvtPDL::getId("tau+");
36
37 // calculate w and costau
38
39 EvtVector4R p4d = p->getDaug(0)->getP4();
40 EvtVector4R p4l = p->getDaug(1)->getP4();
41 EvtVector4R p4n = p->getDaug(2)->getP4();
42 EvtVector4R p4ln(p4l + p4n);
43
44 EvtVector4R p4dln = boostTo(p4d, p4ln, true);
45 EvtVector4R p4lln = boostTo(p4l, p4ln, true);
46
47 const double gmB = p->getP4().mass();
48 const double gmd = p4d.mass();
49 const double gr = gmd / gmB;
50
51 const double q2 = (p4l + p4n).mass2();
52 const double w = (1. + gr * gr - q2 / gmB / gmB) / 2. / gr;
53 // const double w=CalcHelAmp->wfunc(2,q2); avoid possible w<0 caused by the decay width
54 const double costau = p4dln.dot(p4lln) / p4dln.d3mag() / p4lln.d3mag();
55 const double ml = p4l.mass();
56
57 // obtain helicity amplitudes
58 EvtComplex helamp[2]; // tauhel={1,-1}
59 helamp[0] = CalcHelAmp->helAmp(ml, 1, 2, w, costau); // note the parameter order is tauhel, Dhel
60 helamp[1] = CalcHelAmp->helAmp(ml, -1, 2, w, costau); // Dhel=2 ==> D meson
61
62 // lepton theta and phi in l+nu rest frame
63 //const double l_theta=acos(p4lln.get(3)/p4lln.d3mag());
64 //const double l_phi=atan2(p4lln.get(2),p4lln.get(1));
65
66 // spin (in l rest frame) -> helicity (in l+nu rest frame) rotation matrix
67 // ( sp0->hel0 , sp1->hel0 )
68 // ( sp0->hel1 , sp1->hel1 )
69 EvtSpinDensity l_HelFromSp = RotateToHelicityBasisInBoostedFrame(p->getDaug(1),
70 p4ln);
71// l_phi,
72// l_theta,
73// -l_phi);
74
75 // helicity (in l+nu rest frame) -> spin (in l rest frame) rotation matrix
76 // ( hel0->sp0 , hel1->sp0 )
77 // ( hel0->sp1 , hel1->sp1 )
78 EvtComplex l_SpFromHel[2][2]; // {0,1} from {1,-1}
79 EvtId l_num = p->getDaug(1)->getId();
80// if (l_num == EM || l_num == MUM || l_num == TAUM) {
81 l_SpFromHel[0][0] = conj(l_HelFromSp.get(0, 0));
82 l_SpFromHel[0][1] = conj(l_HelFromSp.get(1, 0));
83 l_SpFromHel[1][0] = conj(l_HelFromSp.get(0, 1));
84 l_SpFromHel[1][1] = conj(l_HelFromSp.get(1, 1));
85// } else {
86// l_SpFromHel[0][1] = conj(l_HelFromSp.get(0, 0));
87// l_SpFromHel[0][0] = conj(l_HelFromSp.get(1, 0));
88// l_SpFromHel[1][1] = conj(l_HelFromSp.get(0, 1));
89// l_SpFromHel[1][0] = conj(l_HelFromSp.get(1, 1));
90// }
91
92 // calculate spin amplitudes
93 EvtComplex spinamp[2];
94
95 for (int lsp = 0; lsp < 2; lsp++) {
96 for (int lhel = 0; lhel < 2; lhel++) {
97 // b -> l
98 if (l_num == EM || l_num == MUM || l_num == TAUM) {
99 spinamp[lsp] += l_SpFromHel[lsp][lhel] * helamp[lhel];
100 }
101 // b-bar -> anti-l
102 else {
103 spinamp[lsp] += l_SpFromHel[lsp][lhel] * (lhel == 0 ? +1 : -1) * conj(helamp[1 - lhel]);
104 }
105 }
106 }
107
108 amp.vertex(0, spinamp[0]);
109 amp.vertex(1, spinamp[1]);
110
111 // consistency check
112 double helprob = abs2(helamp[0]) + abs2(helamp[1]);
113 double spinprob = abs2(spinamp[0]) + abs2(spinamp[1]);
114 if (fabs(helprob - spinprob) / helprob > 1E-6 || !finite(helprob) || !finite(spinprob)) {
115 B2ERROR("EvtBSemiTauonicScalarMesonAmplitude total helicity prob does not match with total spin prob.");
116 B2ERROR("helprob: " << helprob << " spinprob: " << spinprob);
117 B2ERROR("w: " << w << " costau: " << costau << " hel probs: " << abs2(helamp[0])
118 << "\t" << abs2(helamp[1])
119 << "\t" << abs2(helamp[0]) + abs2(helamp[1]));
120
121 B2ERROR("w: " << w << " costau: " << costau << " spin probs: " << abs2(spinamp[0])
122 << "\t" << abs2(spinamp[1])
123 << "\t" << abs2(spinamp[0]) + abs2(spinamp[1]));
124
125// EvtGenReport(EVTGEN_ERROR, "EvtGen") <<
126// "EvtBSemiTauonicScalarMesonAmplitude total helicity prob does not match with total spin prob."
127// << std::endl;
128// EvtGenReport(EVTGEN_ERROR, "EvtGen") << "helprob: "<<helprob<<" spinprob: "<<spinprob<<std::endl;
129// EvtGenReport(EVTGEN_ERROR, "EvtGen") << "w: "<<w<<" costau: "<<costau
130// <<" hel probs: "<<abs2(helamp[0])
131// <<"\t"<<abs2(helamp[1])
132// <<"\t"<<abs2(helamp[0])+abs2(helamp[1])<<std::endl;
133//
134// EvtGenReport(EVTGEN_ERROR, "EvtGen") << "w: "<<w<<" costau: "<<costau
135// <<" spin probs: "<<abs2(spinamp[0])
136// <<"\t"<<abs2(spinamp[1])
137// <<"\t"<<abs2(spinamp[0])+abs2(spinamp[1])<<std::endl;
138 // abort();
139 }
140
141 return;
142 }
R E
internal precision of FFTW codelets
EvtSpinDensity RotateToHelicityBasisInBoostedFrame(const EvtParticle *p, EvtVector4R p4boost)
The function calculates the rotation matrix to convert the spin basis to the helicity basis in the bo...

◆ CalcAmp() [2/2]

void CalcAmp ( EvtParticle *  parent,
EvtAmp &  amp,
EvtBSemiTauonicHelicityAmplitudeCalculator HelicityAmplitudeCalculator 
)
overridevirtual

The function calculates the spin dependent amplitude.

Parameters
parenta pointer to the parent particle.
ampa pointer to fill the calculated spin dependent amplitude.
HelicityAmplitudeCalculatora pointer to the calculator of the helicity dependent amplitude. The function calculate the spin dependent amplitude of the semi-tauonic decay to a vector meson (D* meson).

Implements EvtBSemiTauonicAmplitude.

Definition at line 24 of file EvtBSemiTauonicVectorMesonAmplitude.cc.

27 {
28 static EvtId EM = EvtPDL::getId("e-");
29 static EvtId MUM = EvtPDL::getId("mu-");
30 static EvtId TAUM = EvtPDL::getId("tau-");
31// static EvtId EP = EvtPDL::getId("e+");
32// static EvtId MUP = EvtPDL::getId("mu+");
33// static EvtId TAUP = EvtPDL::getId("tau+");
34
35 // calculate w and costau
36
37 EvtVector4R p4d = p->getDaug(0)->getP4();
38 EvtVector4R p4l = p->getDaug(1)->getP4();
39 EvtVector4R p4n = p->getDaug(2)->getP4();
40 EvtVector4R p4ln(p4l + p4n);
41
42 EvtVector4R p4dln = boostTo(p4d, p4ln, true);
43 EvtVector4R p4lln = boostTo(p4l, p4ln, true);
44
45 const double gmB = p->getP4().mass();
46 const double gmd = p4d.mass();
47 const double gr = gmd / gmB;
48
49 const double q2 = (p4l + p4n).mass2();
50 const double w = (1. + gr * gr - q2 / gmB / gmB) / 2. / gr;
51 // const double w=CalcHelAmp->wfunc(1,q2); avoid w<1 caused by the D* width
52 const double costau = p4dln.dot(p4lln) / p4dln.d3mag() / p4lln.d3mag();
53 const double ml = p4l.mass();
54
55 // Set D* mass from its momentum (to take into account fluctuation due to its width)
56 const double orig_mD = CalcHelAmp->getMDst();
57 CalcHelAmp->setMDst(gmd);
58
59 // obtain helicity amplitudes
60 EvtComplex helamp[3][2]; // Dhel={1,0,-1}, tauhel={1,-1}
61 helamp[0][0] = CalcHelAmp->helAmp(ml, 1, 1, w, costau); // note the parameter order is tauhel, Dhel
62 helamp[0][1] = CalcHelAmp->helAmp(ml, -1, 1, w, costau);
63 helamp[1][0] = CalcHelAmp->helAmp(ml, 1, 0, w, costau);
64 helamp[1][1] = CalcHelAmp->helAmp(ml, -1, 0, w, costau);
65 helamp[2][0] = CalcHelAmp->helAmp(ml, 1, -1, w, costau);
66 helamp[2][1] = CalcHelAmp->helAmp(ml, -1, -1, w, costau);
67
68 // lepton theta and phi in l+nu rest frame
69 // const double l_theta=acos(p4lln.get(3)/p4lln.d3mag());
70 // const double l_phi=atan2(p4lln.get(2),p4lln.get(1));
71
72 // spin (in l rest frame) -> helicity (in l+nu rest frame) rotation matrix
73 // ( sp0->hel0 , sp1->hel0 )
74 // ( sp0->hel1 , sp1->hel1 )
75 EvtSpinDensity l_HelFromSp = RotateToHelicityBasisInBoostedFrame(p->getDaug(1),
76 p4ln);
77// l_phi,
78// l_theta,
79// -l_phi);
80
81 // helicity (in l+nu rest frame) -> spin (in l rest frame) rotation matrix
82 // ( hel0->sp0 , hel1->sp0 )
83 // ( hel0->sp1 , hel1->sp1 )
84 EvtComplex l_SpFromHel[2][2]; // {0,1} from {1,-1}
85 EvtId l_num = p->getDaug(1)->getId();
86// if (l_num == EM || l_num == MUM || l_num == TAUM) {
87 l_SpFromHel[0][0] = conj(l_HelFromSp.get(0, 0));
88 l_SpFromHel[0][1] = conj(l_HelFromSp.get(1, 0));
89 l_SpFromHel[1][0] = conj(l_HelFromSp.get(0, 1));
90 l_SpFromHel[1][1] = conj(l_HelFromSp.get(1, 1));
91// } else {
92// l_SpFromHel[0][1] = conj(l_HelFromSp.get(0, 0));
93// l_SpFromHel[0][0] = conj(l_HelFromSp.get(1, 0));
94// l_SpFromHel[1][1] = conj(l_HelFromSp.get(0, 1));
95// l_SpFromHel[1][0] = conj(l_HelFromSp.get(1, 1));
96// }
97
98 // meson spin state to helicity state
99 const double D_theta = acos(p4d.get(3) / p4d.d3mag());
100 const double D_phi = atan2(p4d.get(2), p4d.get(1));
101 EvtSpinDensity D_HelFromSp = p->getDaug(0)->rotateToHelicityBasis(D_phi, D_theta, -D_phi);
102
103 EvtComplex D_SpFromHel[3][3]; // {0,1,2} from {1,0,-1}
104 D_SpFromHel[0][0] = conj(D_HelFromSp.get(0, 0));
105 D_SpFromHel[0][1] = conj(D_HelFromSp.get(1, 0));
106 D_SpFromHel[0][2] = conj(D_HelFromSp.get(2, 0));
107 D_SpFromHel[1][0] = conj(D_HelFromSp.get(0, 1));
108 D_SpFromHel[1][1] = conj(D_HelFromSp.get(1, 1));
109 D_SpFromHel[1][2] = conj(D_HelFromSp.get(2, 1));
110 D_SpFromHel[2][0] = conj(D_HelFromSp.get(0, 2));
111 D_SpFromHel[2][1] = conj(D_HelFromSp.get(1, 2));
112 D_SpFromHel[2][2] = conj(D_HelFromSp.get(2, 2));
113
114 // calculate spin amplitudes
115 EvtComplex spinamp[3][2];
116
117 for (int dsp = 0; dsp < 3; dsp++) {
118 for (int lsp = 0; lsp < 2; lsp++) {
119 for (int dhel = 0; dhel < 3; dhel++) {
120 for (int lhel = 0; lhel < 2; lhel++) {
121 // b -> l, D*+
122 if (l_num == EM || l_num == MUM || l_num == TAUM) {
123 spinamp[dsp][lsp] += l_SpFromHel[lsp][lhel] * D_SpFromHel[dsp][dhel] * helamp[dhel][lhel];
124 }
125 // b-bar -> anti-l, D*-
126 else {
127 spinamp[dsp][lsp] += l_SpFromHel[lsp][lhel] * D_SpFromHel[dsp][dhel] * (lhel == 0 ? +1 : -1) * (dhel == 1 ? +1 : -1)
128 * conj(helamp[2 - dhel][1 - lhel]);
129 }
130 }
131 }
132 }
133 }
134
135 amp.vertex(0, 0, spinamp[0][0]);
136 amp.vertex(0, 1, spinamp[0][1]);
137 amp.vertex(1, 0, spinamp[1][0]);
138 amp.vertex(1, 1, spinamp[1][1]);
139 amp.vertex(2, 0, spinamp[2][0]);
140 amp.vertex(2, 1, spinamp[2][1]);
141
142 // Set D* mass to its original value
143 CalcHelAmp->setMDst(orig_mD);
144
145 // consistency check
146 double helprob = abs2(helamp[0][0]) + abs2(helamp[0][1]) + abs2(helamp[1][0]) + abs2(helamp[1][1]) + abs2(helamp[2][0])
147 + abs2(helamp[2][1]);
148 double spinprob = abs2(spinamp[0][0]) + abs2(spinamp[0][1]) + abs2(spinamp[1][0]) + abs2(spinamp[1][1])
149 + abs2(spinamp[2][0]) + abs2(spinamp[2][1]);
150 if (fabs(helprob - spinprob) / helprob > 1E-6 || !finite(helprob) || !finite(spinprob)) {
151 B2ERROR("EvtBSemiTauonicVectorMesonAmplitude total helicity prob does not match with total spin prob, or nan.");
152 B2ERROR("helprob: " << helprob << " spinprob: " << spinprob);
153 B2ERROR("w: " << w << " costau: " << costau << " hel probs: "
154 << abs2(helamp[0][0]) << "\t" << abs2(helamp[0][1]) << "\t"
155 << abs2(helamp[1][0]) << "\t" << abs2(helamp[1][1]) << "\t"
156 << abs2(helamp[2][0]) << "\t" << abs2(helamp[2][1]) << "\t"
157 << abs2(helamp[0][0]) + abs2(helamp[0][1]) + abs2(helamp[1][0]) + abs2(helamp[1][1]) + abs2(helamp[2][0]) + abs2(helamp[2][1])
158 );
159
160 B2ERROR("w: " << w << " costau: " << costau << " spin probs: "
161 << abs2(spinamp[0][0]) << "\t" << abs2(spinamp[0][1]) << "\t"
162 << abs2(spinamp[1][0]) << "\t" << abs2(spinamp[1][1]) << "\t"
163 << abs2(spinamp[2][0]) << "\t" << abs2(spinamp[2][1]) << "\t"
164 << abs2(spinamp[0][0]) + abs2(spinamp[0][1]) + abs2(spinamp[1][0]) + abs2(spinamp[1][1]) + abs2(spinamp[2][0]) + abs2(spinamp[2][1])
165 );
166
167// EvtGenReport(EVTGEN_ERROR, "EvtGen") <<
168// "EvtBSemiTauonicVectorMesonAmplitude total helicity prob does not match with total spin prob, or nan."
169// << std::endl;
170// EvtGenReport(EVTGEN_ERROR, "EvtGen") "helprob: "<<helprob<<" spinprob: "<<spinprob<< std::endl;
171// EvtGenReport(EVTGEN_ERROR, "EvtGen") "w: "<<w<<" costau: "<<costau<<" hel probs: "
172// <<abs2(helamp[0][0])<<"\t"<<abs2(helamp[0][1])<<"\t"
173// <<abs2(helamp[1][0])<<"\t"<<abs2(helamp[1][1])<<"\t"
174// <<abs2(helamp[2][0])<<"\t"<<abs2(helamp[2][1])<<"\t"
175// <<abs2(helamp[0][0]) + abs2(helamp[0][1]) + abs2(helamp[1][0]) + abs2(helamp[1][1]) + abs2(helamp[2][0]) + abs2(helamp[2][1])
176// <<std::endl;
177//
178// EvtGenReport(EVTGEN_ERROR, "EvtGen") "w: "<<w<<" costau: "<<costau<<" spin probs: "
179// <<abs2(spinamp[0][0])<<"\t"<<abs2(spinamp[0][1])<<"\t"
180// <<abs2(spinamp[1][0])<<"\t"<<abs2(spinamp[1][1])<<"\t"
181// <<abs2(spinamp[2][0])<<"\t"<<abs2(spinamp[2][1])<<"\t"
182// <<abs2(spinamp[0][0]) + abs2(spinamp[0][1]) + abs2(spinamp[1][0]) + abs2(spinamp[1][1]) + abs2(spinamp[2][0]) + abs2(spinamp[2][1])
183// <<std::endl;
184
185 // Debugging information
186// fprintf(stderr, "q2 by Helamp: %g\n", CalcHelAmp->q2(1, w));
187// fprintf(stderr, "helampSM by Helamp: %g\n", CalcHelAmp->helampSM(ml, 1, 1, w, costau));
188// fprintf(stderr, "Lep by Helamp: %g\n", CalcHelAmp->Lep(ml, 1, 1, CalcHelAmp->q2(1, w), costau));
189// fprintf(stderr, "HadV2 by Helamp: %g\n", CalcHelAmp->HadV2(1, 1, w));
190// fprintf(stderr, "v by Helamp: %g\n", CalcHelAmp->v(ml, CalcHelAmp->v(ml, CalcHelAmp->q2(1, w))));
191//
192// fprintf(stderr, "B mass : %g \t nominal %g\n", gmB, EvtPDL::getMeanMass(p->getId()));
193// fprintf(stderr, "D mass : %g \t nominal %g\n", gmd, EvtPDL::getMeanMass(p->getDaug(0)->getId()));
194// fprintf(stderr, "lepton mass : %g \t nominal %g\n", ml, EvtPDL::getMeanMass(p->getDaug(1)->getId()));
195
196 // abort();
197 }
198
199 return;
200 }

◆ CalcMaxProb()

double CalcMaxProb ( EvtId  parent,
EvtId  meson,
EvtId  lepton,
EvtId  nudaug,
EvtBSemiTauonicHelicityAmplitudeCalculator HelicityAmplitudeCalculator 
)

The function calculates the maximum probability.

Parameters
parenta ID of the parent meson.
mesona ID of the daughter meson.
leptona ID of the daughter lepton.
nudauga ID of the daughter neutrino.
HelicityAmplitudeCalculatora pointer to the calculator of the helicity dependent amplitude.
Returns
the maxmum probability multiplied by 1.1. The function scan the q^2 and angle between the daughter meson and lepton and search for the the maximum probability value.

Definition at line 77 of file EvtBSemiTauonicAmplitude.cc.

80 {
81 //This routine takes the arguements parent, meson, and lepton
82 //number, and a form factor model, and returns a maximum
83 //probability for this semileptonic form factor model. A
84 //brute force method is used. The 2D cos theta lepton and
85 //q2 phase space is probed.
86
87 //Start by declaring a particle at rest.
88
89 //It only makes sense to have a scalar parent. For now.
90 //This should be generalized later.
91
92 EvtScalarParticle* scalar_part;
93 EvtParticle* root_part;
94
95 scalar_part = new EvtScalarParticle;
96
97 //cludge to avoid generating random numbers!
98 scalar_part->noLifeTime();
99
100 EvtVector4R p_init;
101
102 p_init.set(EvtPDL::getMass(parent), 0.0, 0.0, 0.0);
103 scalar_part->init(parent, p_init);
104 //root_part = (EvtParticle*)scalar_part;
105 root_part = static_cast<EvtParticle*>(scalar_part);
106
107 root_part->setDiagonalSpinDensity();
108
109 EvtParticle* daughter, *lep, *trino;
110
111 EvtAmp amp;
112
113 EvtId listdaug[3];
114 listdaug[0] = meson;
115 listdaug[1] = lepton;
116 listdaug[2] = nudaug;
117
118 amp.init(parent, 3, listdaug);
119
120 root_part->makeDaughters(3, listdaug);
121 daughter = root_part->getDaug(0);
122 lep = root_part->getDaug(1);
123 trino = root_part->getDaug(2);
124
125 //cludge to avoid generating random numbers!
126 daughter->noLifeTime();
127 lep->noLifeTime();
128 trino->noLifeTime();
129
130
131 //Initial particle is unpolarized, well it is a scalar so it is
132 //trivial
133 EvtSpinDensity rho;
134 rho.setDiag(root_part->getSpinStates());
135
136 double mass[3];
137
138 double m = root_part->mass();
139
140 EvtVector4R p4meson, p4lepton, p4nu, p4w;
141
142 double q2, elepton, plepton;
143 int i, j;
144 double erho, prho, costl;
145
146 double maxfoundprob = 0.0;
147 double prob;
148 int massiter;
149
150 for (massiter = 0; massiter < 3; massiter++) {
151
152 mass[0] = EvtPDL::getMeanMass(meson);
153 mass[1] = EvtPDL::getMeanMass(lepton);
154 mass[2] = EvtPDL::getMeanMass(nudaug);
155 if (massiter == 1) {
156 mass[0] = EvtPDL::getMinMass(meson);
157 }
158 if (massiter == 2) {
159 mass[0] = EvtPDL::getMaxMass(meson);
160 if ((mass[0] + mass[1] + mass[2]) > m) mass[0] = m - mass[1] - mass[2] - 0.00001;
161 }
162
163 double q2min = mass[1] * mass[1]; // limit to minimum=lepton mass square
164 double q2max = (m - mass[0]) * (m - mass[0]);
165 double dq2 = (q2max - q2min) / 25;
166// std::cout<<"m: "<<m<<" mass[0]: "<<mass[0]<<" q2min: "<<q2min<<" q2max: "<<q2max<<std::endl;
167
168 //loop over q2
169
170 for (i = 0; i < 25; i++) {
171 q2 = q2min + (i + 0.5) * dq2; // <-- !! not start from unphysical q2=0 !!
172
173 erho = (m * m + mass[0] * mass[0] - q2) / (2.0 * m);
174
175 prho = sqrt(erho * erho - mass[0] * mass[0]);
176
177 p4meson.set(erho, 0.0, 0.0, -1.0 * prho);
178 p4w.set(m - erho, 0.0, 0.0, prho);
179
180// std::cout<<"q2: "<<q2<<std::endl;
181// std::cout<<"p4meson: "<<p4meson<<std::endl;
182
183 //This is in the W rest frame
184 elepton = (q2 + mass[1] * mass[1]) / (2.0 * sqrt(q2));
185 plepton = sqrt(elepton * elepton - mass[1] * mass[1]);
186// std::cout<<"elepton: "<<elepton<<" plepton: "<<plepton<<std::endl;
187
188 double probctl[3];
189
190 for (j = 0; j < 3; j++) {
191
192 costl = 0.99 * (j - 1.0);
193
194 //These are in the W rest frame. Need to boost out into
195 //the B frame.
196 p4lepton.set(elepton, 0.0,
197 plepton * sqrt(1.0 - costl * costl), plepton * costl);
198 p4nu.set(plepton, 0.0,
199 -1.0 * plepton * sqrt(1.0 - costl * costl), -1.0 * plepton * costl);
200
201 EvtVector4R boost((m - erho), 0.0, 0.0, 1.0 * prho);
202 p4lepton = boostTo(p4lepton, boost);
203 p4nu = boostTo(p4nu, boost);
204
205 //Now initialize the daughters...
206
207 daughter->init(meson, p4meson);
208 lep->init(lepton, p4lepton);
209 trino->init(nudaug, p4nu);
210
211 CalcAmp(root_part, amp, CalcHelAmp);
212
213 //Now find the probability at this q2 and cos theta lepton point
214 //and compare to maxfoundprob.
215
216 //Do a little magic to get the probability!!
217 prob = rho.normalizedProb(amp.getSpinDensity());
218
219 probctl[j] = prob;
220 }
221
222 //probclt contains prob at ctl=-1,0,1.
223 //prob=a+b*ctl+c*ctl^2
224
225 double a = probctl[1];
226 double b = 0.5 * (probctl[2] - probctl[0]);
227 double c = 0.5 * (probctl[2] + probctl[0]) - probctl[1];
228
229 prob = probctl[0];
230 if (probctl[1] > prob) prob = probctl[1];
231 if (probctl[2] > prob) prob = probctl[2];
232
233 if (fabs(c) > 1e-20) {
234 double ctlx = -0.5 * b / c;
235 if (fabs(ctlx) < 1.0) {
236 double probtmp = a + b * ctlx + c * ctlx * ctlx;
237 if (probtmp > prob) prob = probtmp;
238 }
239
240 }
241
242 if (prob > maxfoundprob) {
243 maxfoundprob = prob;
244 }
245
246 }
247 if (EvtPDL::getWidth(meson) <= 0.0) {
248 //if the particle is narrow dont bother with changing the mass.
249 massiter = 4;
250 }
251
252 }
253 root_part->deleteTree();
254
255 maxfoundprob *= 1.1;
256 return maxfoundprob;
257
258 }
virtual void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator)=0
The function calculates the spin dependent amplitude.
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28

◆ checkVariable()

void checkVariable ( const std::string &  name,
const std::map< int, bool > &  allowed,
ParticleGun::EDistribution dist,
std::vector< double > &  pars 
)
protected

check one of the variables given a list of allowed and excluded distributions

Definition at line 41 of file particlegun.cc.

43 {
44 parameters = ParticleGun::Parameters();
45 parameters.pdgCodes = {11};
46 for (auto item : allowed) {
47 dist = (ParticleGun::EDistribution) item.first;
48 //We need at least for parameters for some distributions.
49 pars = {1, 1, 1, 1};
50 //Check if distribution is accepted
51 ASSERT_EQ(item.second, pgun.setParameters(parameters)) << name << ": " << item.first << ", " << item.second;
52 if (item.second) {
53 //And if it is accepted, check number of parameters
54 for (int i = 0; i < 10; ++i) {
55 pars.resize(i, 1);
56 int minpar = (item.first == ParticleGun::c_fixedValue) ? 1 : 2;
57 //Polyline needs at least two points, so four parameters
59 || dist == ParticleGun::c_polylineCosDistribution) minpar = 4;
60 bool valid = i >= minpar;
61 //discrete and polyline need even number of parameters
64 valid &= i % 2 == 0;
65 }
66 //Check if numberof parameters is accepted
67 ASSERT_EQ(valid, pgun.setParameters(parameters)) << name << " " << pars.size();
68 }
69 }
70 //Check that zeros are forbidden for inversePt spectrum
71 if (item.first == ParticleGun::c_inversePtDistribution && item.second) {
72 pars = {1, 1};
73 ASSERT_TRUE(pgun.setParameters(parameters));
74 pars = {1, 0};
75 ASSERT_FALSE(pgun.setParameters(parameters)) << name << ": " << parameters.momentumParams[0] << ", " <<
77 pars = {0, 1};
78 ASSERT_FALSE(pgun.setParameters(parameters)) << name << ": " << parameters.momentumParams[0] << ", " <<
80 }
81 //Check polyline stuff
83 || dist == ParticleGun::c_polylineCosDistribution) && item.second) {
84 pars = {0, 1, 0, 1};
85 ASSERT_TRUE(pgun.setParameters(parameters));
86 //x needs to be ascending
87 pars = {1, 0, 0, 1};
88 ASSERT_FALSE(pgun.setParameters(parameters));
89 //at least one positive y value
90 pars = {0, 1, 0, 0};
91 ASSERT_FALSE(pgun.setParameters(parameters));
92 //no negative y values
93 pars = {0, 1, 1, -1};
94 ASSERT_FALSE(pgun.setParameters(parameters));
95 }
96 //Check discrete spectrum has non-negative weights
98 pars = {0, 1, 0, 1};
99 ASSERT_TRUE(pgun.setParameters(parameters));
100 //at least one weight must be positive
101 pars = {0, 1, 0, 0};
102 ASSERT_FALSE(pgun.setParameters(parameters));
103 //no negative weights
104 pars = {0, 1, 1, -1};
105 ASSERT_FALSE(pgun.setParameters(parameters));
106 }
107 }
108 std::cout << std::endl;
109 }
ParticleGun::Parameters parameters
Variable parameters.
Definition: particlegun.cc:34
ParticleGun pgun
Variable of pgun.
Definition: particlegun.cc:32
EDistribution
enum containing all known distributions available for generation of values
Definition: ParticleGun.h:29
@ c_polylineDistribution
Distribution given as list of (x,y) points.
Definition: ParticleGun.h:57
@ c_polylinePtDistribution
Same as polylineDistribution but for the transverse momentum.
Definition: ParticleGun.h:59
@ c_discreteSpectrum
Discrete spectrum, parameters are first the values and then the weights (non-negative) for each value...
Definition: ParticleGun.h:49
@ c_fixedValue
Fixed value, no random generation at all, 1 parameter.
Definition: ParticleGun.h:31
@ c_inversePtDistribution
Distribution uniform in the inverse pt distribution, that is uniform in track curvature.
Definition: ParticleGun.h:53
@ c_polylineCosDistribution
Same as polylineDistribution but for the cosine of the angle.
Definition: ParticleGun.h:61
bool setParameters(const Parameters &parameters)
Set the parameters for generating the Particles.
std::vector< int > pdgCodes
List of PDG particle codes to pick from when generating particles.
Definition: ParticleGun.h:83
std::vector< double > momentumParams
Parameters for the momentum generation, meaning depends on chosen distribution.
Definition: ParticleGun.h:85

◆ chkDhel()

bool chkDhel ( int  Dhel) const
private

sanity checkers

Function to check if Dhel is in the valid range.

Parameters
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.

Definition at line 528 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

529 {
530 if (Dhel == -1 || Dhel == 0 || Dhel == 1 || Dhel == 2)return true;
531 else return false;
532 }

◆ chktauhel()

bool chktauhel ( int  tauhel) const
private

Function to check if tauhel is in the valid range.

Parameters
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.

Definition at line 540 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

541 {
542 if (tauhel == -1 || tauhel == 1)return true;
543 else return false;
544 }

◆ chkwhel()

bool chkwhel ( int  whel) const
private

Function to check if whel is in the valid range.

Parameters
whelhelicity of the virtual vector boson {+1,0,1,2}.

Definition at line 534 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

535 {
536 if (whel == -1 || whel == 0 || whel == 1 || whel == 2)return true;
537 else return false;
538 }

◆ clone() [1/12]

EvtDecayBase * clone ( )

Clone the decay of B0toKsKK.

Definition at line 36 of file EvtB0toKsKK.cc.

37 {
38 return new EvtB0toKsKK;
39 }

◆ clone() [2/12]

EvtDecayBase * clone ( )

Clones module.

Definition at line 44 of file EvtBCL.cc.

45 {
46 return new EvtBCL;
47 }
EvtBCL()
Default constructor.
Definition: EvtBCL.cc:29

◆ clone() [3/12]

EvtDecayBase * clone ( )

The function makes a copy of an EvtBSTD object.

Returns
a pointer to a new uninitialized EvtBSemiTauonic object

Definition at line 47 of file EvtBSemiTauonic.cc.

48 {
49 return new EvtBSemiTauonic;
50 }
EvtBSemiTauonic()
The default constructor

◆ clone() [4/12]

EvtDecayBase * clone ( )

The function makes a copy of an EvtBSemiTauonic2HDMType2 object.

Returns
a pointer to a new uninitialized EvtBSemiTauonic2HDMType2 object

Definition at line 47 of file EvtBSemiTauonic2HDMType2.cc.

48 {
49 return new EvtBSemiTauonic2HDMType2;
50 }
EvtBSemiTauonic2HDMType2()
The default constructor

◆ clone() [5/12]

EvtDecayBase * clone ( )

The function which makes a copy of the model.

Definition at line 41 of file EvtBtoXsnunu_FERMI.cc.

42 {
43 return new EvtBtoXsnunu_FERMI;
44 }

◆ clone() [6/12]

EvtDecayBase * clone ( )

Makes a copy of the pointer to the class.

Definition at line 34 of file EvtEtaFullDalitz.cc.

35 {
36
37 return new EvtEtaFullDalitz;
38
39 }
EvtEtaFullDalitz()
Default Constructor.

◆ clone() [7/12]

EvtDecayBase * clone ( )

Makes a copy of the class object.

Definition at line 49 of file EvtEtaPi0Dalitz.cc.

50 {
51
52 return new EvtEtaPi0Dalitz;
53
54 }
EvtEtaPi0Dalitz()
Default constructor.

◆ clone() [8/12]

EvtDecayBase * clone ( )

Returns a copy of the class object.

Definition at line 51 of file EvtEtaPrimeDalitz.cc.

52 {
53
54 return new EvtEtaPrimeDalitz;
55
56 }
EvtEtaPrimeDalitz()
Default constructor.

◆ clone() [9/12]

EvtDecayBase * clone ( )

The function which makes a copy of the model.

Definition at line 41 of file EvtFlatDaughter.cc.

42 {
43 return new EvtFlatDaughter;
44 }
EvtFlatDaughter()
Constructor.

◆ clone() [10/12]

EvtDecayBase * clone ( )

The function which makes a copy of the model.

Definition at line 45 of file EvtKnunu.cc.

46 {
47 return new EvtKnunu;
48 }
EvtKnunu()
Constructor.
Definition: EvtKnunu.h:33

◆ clone() [11/12]

EvtDecayBase * clone ( )

The function which makes a copy of the model.

Definition at line 44 of file EvtKstarnunu_REV.cc.

45 {
46 return new EvtKstarnunu_REV;
47 }
EvtKstarnunu_REV()
Constructor.

◆ clone() [12/12]

EvtDecayBase * clone ( )

Clone the decay.

Definition at line 35 of file EvtPhspCP.cc.

36 {
37 return new EvtPhspCP;
38 }

◆ decay() [1/12]

void decay ( EvtParticle *  p)

Member of particle in EvtGen.

Definition at line 146 of file EvtB0toKsKK.cc.

147 {
148 // Btag
149 static EvtId B0 = EvtPDL::getId("B0");
150 static EvtId B0B = EvtPDL::getId("anti-B0");
151
152 double t;
153 EvtId other_b;
154
155 //std::cout << EvtCPUtil::getInstance()->getMixingType() << std::endl;
156 //EvtCPUtil::getInstance()->setMixingType(0);
157 EvtCPUtil::getInstance()->OtherB(p, t, other_b, 0.5);
158 //EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.4);
159 //EvtCPUtil::getInstance()->OtherB(p,t,other_b);
160
161 // Brec
162 p->initializePhaseSpace(getNDaug(), getDaugs());
163 EvtVector4R p4ks, p4kp, p4km;
164 if (p->getId() == B0) {
165 p4ks = p->getDaug(0)->getP4();
166 p4kp = p->getDaug(1)->getP4();
167 p4km = p->getDaug(2)->getP4();
170 /*p4ks = p->getDaug(0)->getP4();
171 p4kp = p->getDaug(2)->getP4();
172 p4km = p->getDaug(1)->getP4();*/
173 } else {
174 p4ks = p->getDaug(0)->getP4();
175 p4kp = p->getDaug(2)->getP4();
176 p4km = p->getDaug(1)->getP4();
178 /*p4ks = p->getDaug(0)->getP4();
179 p4kp = p->getDaug(1)->getP4();
180 p4km = p->getDaug(2)->getP4();*/
181 }
182
183 /*std::cout << (p4ks + p4kp).mass() << " " << (p4ks + p4km).mass() << " "
184 << (p4kp + p4km).mass() << std::endl;
185 if( p->getId() == other_b )
186 std::cout << "same flavour" << std::endl;
187 std::cout << p->getId() << " --> " << p->getDaug(0)->getId() << " "
188 << p->getDaug(1)->getId() << " " << p->getDaug(2)->getId()
189 << std::endl;
190 std::cout << B0 << std::endl;
191 std::cout << B0B << std::endl;
192 std::cout << other_b << std::endl;*/
193 /*std::cout << p->getLifetime() << std::endl;
194 std::cout << p->get4Pos() << std::endl;
195 std::cout << p->getP4() << std::endl;
196 std::cout << p->getP4Lab() << std::endl;
197 EvtParticle *parent=p->getParent();
198 if (parent->getDaug(0)!=p){
199 std::cout << parent->getDaug(0)->getLifetime() << std::endl;
200 std::cout << parent->getDaug(0)->get4Pos() << std::endl;
201 std::cout << parent->getDaug(0)->getP4() << std::endl;
202 std::cout << parent->getDaug(0)->getP4Lab() << std::endl;
203 std::cout << parent->getDaug(0)->getP4Lab()+p->getP4Lab() << std::endl;
204 }
205 else{
206 std::cout << parent->getDaug(1)->getLifetime() << std::endl;
207 std::cout << parent->getDaug(1)->get4Pos() << std::endl;
208 std::cout << parent->getDaug(1)->getP4() << std::endl;
209 std::cout << parent->getDaug(1)->getP4Lab() << std::endl;
210 std::cout << parent->getDaug(1)->getP4Lab()+p->getP4Lab() << std::endl;
211 }
212 std::cout << parent->getP4() << std::endl;
213 std::cout << t << std::endl;*/
214
215 // Relative amplides and phases with direct CP violation
216 a_f0ks_ =
217 (1.0 + getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
218 getArg(1) * sin(getArg(2) * M_PI / 180.0));
219 a_phiks_ =
220 (1.0 + getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
221 getArg(5) * sin(getArg(6) * M_PI / 180.0));
222 a_fxks_ =
223 (1.0 + getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
224 getArg(9) * sin(getArg(10) * M_PI / 180.0));
225 a_chic0ks_ =
226 (1.0 + getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
227 getArg(13) * sin(getArg(14) * M_PI / 180.0));
228 a_kpkmnr_ =
229 (1.0 + getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
230 getArg(17) * sin(getArg(18) * M_PI / 180.0));
231 a_kskpnr_ =
232 (1.0 + getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
233 getArg(22) * sin(getArg(23) * M_PI / 180.0));
234 a_kskmnr_ =
235 (1.0 + getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
236 getArg(27) * sin(getArg(28) * M_PI / 180.0));
237
238 abar_f0ks_ =
239 (1.0 - getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
240 getArg(1) * sin(getArg(2) * M_PI / 180.0));
242 (1.0 - getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
243 getArg(5) * sin(getArg(6) * M_PI / 180.0));
244 abar_fxks_ =
245 (1.0 - getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
246 getArg(9) * sin(getArg(10) * M_PI / 180.0));
248 (1.0 - getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
249 getArg(13) * sin(getArg(14) * M_PI / 180.0));
251 (1.0 - getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
252 getArg(17) * sin(getArg(18) * M_PI / 180.0));
254 (1.0 - getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
255 getArg(22) * sin(getArg(23) * M_PI / 180.0));
257 (1.0 - getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
258 getArg(27) * sin(getArg(28) * M_PI / 180.0));
259
260 // Mixing-induced CP asymmetry
261 const double pCP_f0ks = getArg(4) * M_PI / 180.0;
262 const double pCP_phiks = getArg(8) * M_PI / 180.0;
263 const double pCP_fxks = getArg(12) * M_PI / 180.0;
264 const double pCP_chic0ks = getArg(16) * M_PI / 180.0;
265 const double pCP_kpkmnr = getArg(20) * M_PI / 180.0;
266 const double pCP_kskpnr = getArg(25) * M_PI / 180.0;
267 const double pCP_kskmnr = getArg(30) * M_PI / 180.0;
268
269 if (other_b == B0) {
270 a_f0ks_ *=
271 EvtComplex(cos(+2.0 * pCP_f0ks), sin(+2.0 * pCP_f0ks));
272 a_phiks_ *=
273 EvtComplex(cos(+2.0 * pCP_phiks), sin(+2.0 * pCP_phiks));
274 a_fxks_ *=
275 EvtComplex(cos(+2.0 * pCP_fxks), sin(+2.0 * pCP_fxks));
276 a_chic0ks_ *=
277 EvtComplex(cos(+2.0 * pCP_chic0ks), sin(+2.0 * pCP_chic0ks));
278 a_kpkmnr_ *=
279 EvtComplex(cos(+2.0 * pCP_kpkmnr), sin(+2.0 * pCP_kpkmnr));
280 a_kskpnr_ *=
281 EvtComplex(cos(+2.0 * pCP_kskpnr), sin(+2.0 * pCP_kskpnr));
282 a_kskmnr_ *=
283 EvtComplex(cos(+2.0 * pCP_kskmnr), sin(+2.0 * pCP_kskmnr));
284 }
285 if (other_b == B0B) {
286 abar_f0ks_ *=
287 EvtComplex(cos(-2.0 * pCP_f0ks), sin(-2.0 * pCP_f0ks));
288 abar_phiks_ *=
289 EvtComplex(cos(-2.0 * pCP_phiks), sin(-2.0 * pCP_phiks));
290 abar_fxks_ *=
291 EvtComplex(cos(-2.0 * pCP_fxks), sin(-2.0 * pCP_fxks));
293 EvtComplex(cos(-2.0 * pCP_chic0ks), sin(-2.0 * pCP_chic0ks));
294 abar_kpkmnr_ *=
295 EvtComplex(cos(-2.0 * pCP_kpkmnr), sin(-2.0 * pCP_kpkmnr));
296 abar_kskpnr_ *=
297 EvtComplex(cos(-2.0 * pCP_kskpnr), sin(-2.0 * pCP_kskpnr));
298 abar_kskmnr_ *=
299 EvtComplex(cos(-2.0 * pCP_kskmnr), sin(-2.0 * pCP_kskmnr));
300 }
301
302 //Form Factors
303 EvtComplex Amp_f0ks = A_f0ks(p4ks, p4kp, p4km);
304 EvtComplex Amp_phiks = A_phiks(p4ks, p4kp, p4km);
305 EvtComplex Amp_fxks = A_fxks(p4ks, p4kp, p4km);
306 EvtComplex Amp_chic0ks = A_chic0ks(p4ks, p4kp, p4km);
307 EvtComplex Amp_kpkmnr = A_kknr(p4kp, p4km, alpha_kpkmnr);
308 EvtComplex Amp_kskpnr = A_kknr(p4ks, p4kp, alpha_kskpnr);
309 EvtComplex Amp_kskmnr = A_kknr(p4ks, p4km, alpha_kskmnr);
310
311 EvtComplex Ampbar_f0ks = A_f0ks(p4ks, p4km, p4kp);
312 EvtComplex Ampbar_phiks = A_phiks(p4ks, p4km, p4kp);
313 EvtComplex Ampbar_fxks = A_fxks(p4ks, p4km, p4kp);
314 EvtComplex Ampbar_chic0ks = A_chic0ks(p4ks, p4km, p4kp);
315 EvtComplex Ampbar_kpkmnr = A_kknr(p4km, p4kp, alpha_kpkmnr);
316 EvtComplex Ampbar_kskpnr = A_kknr(p4ks, p4km, alpha_kskpnr);
317 EvtComplex Ampbar_kskmnr = A_kknr(p4ks, p4kp, alpha_kskmnr);
318
319 const EvtComplex A_B0toKsKK =
320 (a_f0ks_ * Amp_f0ks) +
321 (a_phiks_ * Amp_phiks) +
322 (a_fxks_ * Amp_fxks) +
323 (a_chic0ks_ * Amp_chic0ks) +
324 (a_kpkmnr_ * Amp_kpkmnr) +
325 (a_kskpnr_ * Amp_kskpnr) +
326 (a_kskmnr_ * Amp_kskmnr);
327
328 const EvtComplex Abar_B0toKsKK =
329 (abar_f0ks_ * Ampbar_f0ks) +
330 (abar_phiks_ * Ampbar_phiks) +
331 (abar_fxks_ * Ampbar_fxks) +
332 (abar_chic0ks_ * Ampbar_chic0ks) +
333 (abar_kpkmnr_ * Ampbar_kpkmnr) +
334 (abar_kskpnr_ * Ampbar_kskpnr) +
335 (abar_kskmnr_ * Ampbar_kskmnr);
336
337 // CP asymmetry
338 const double dm = getArg(0);
339
340 EvtComplex amp;
341 if (other_b == B0B) {
342 amp =
343 A_B0toKsKK * cos(dm * t / (2.0 * EvtConst::c)) +
344 EvtComplex(0.0, 1.0) * Abar_B0toKsKK * sin(dm * t / (2.0 * EvtConst::c));
345 }
346 if (other_b == B0) {
347 amp =
348 A_B0toKsKK *
349 EvtComplex(0.0, 1.0) * sin(dm * t / (2.0 * EvtConst::c)) +
350 Abar_B0toKsKK * cos(dm * t / (2.0 * EvtConst::c));
351 }
352
353 if (abs2(amp) > 50000.0)
354 debugfile_ << abs2(amp) << std::endl;
355 vertex(amp);
356
357 return;
358 }
EvtComplex abar_kskpnr_
Variable member abar_kskpnr_
Definition: EvtB0toKsKK.h:100
double alpha_kskmnr
Variable member alpha_kskmnr.
Definition: EvtB0toKsKK.h:105
EvtComplex a_chic0ks_
Variable member a_chic0ks_.
Definition: EvtB0toKsKK.h:90
EvtComplex a_fxks_
Variable member a_fxks_
Definition: EvtB0toKsKK.h:89
std::ofstream debugfile_
debuging stream
Definition: EvtB0toKsKK.h:107
EvtComplex a_f0ks_
<Variable names for form factors
Definition: EvtB0toKsKK.h:87
double alpha_kpkmnr
Variable member alpha_kpkmnr.
Definition: EvtB0toKsKK.h:103
EvtComplex abar_fxks_
Variable member abar_fxks_
Definition: EvtB0toKsKK.h:97
EvtComplex abar_f0ks_
Variable member abar_f0ks_
Definition: EvtB0toKsKK.h:95
double alpha_kskpnr
Variable member alpha_kskpnr.
Definition: EvtB0toKsKK.h:104
EvtComplex abar_chic0ks_
Variable member abar_chic0ks_.
Definition: EvtB0toKsKK.h:98
EvtComplex a_kskmnr_
Variable member a_kskmnr_.
Definition: EvtB0toKsKK.h:93
EvtComplex abar_kskmnr_
Variable member abar_kskmnr_
Definition: EvtB0toKsKK.h:101
EvtComplex a_kpkmnr_
Variable member a_kpkmnr_.
Definition: EvtB0toKsKK.h:91
EvtComplex abar_phiks_
Variable member abar_phiks_.
Definition: EvtB0toKsKK.h:96
EvtComplex a_kskpnr_
Variable member a_kskpnr_.
Definition: EvtB0toKsKK.h:92
EvtComplex a_phiks_
Variable member a_phiks_
Definition: EvtB0toKsKK.h:88
EvtComplex abar_kpkmnr_
Variable member abar_kpkmnr_
Definition: EvtB0toKsKK.h:99
double p(const double &mab, const double &M, const double &mc)
Constants p
Definition: EvtB0toKsKK.cc:583

◆ decay() [2/12]

void decay ( EvtParticle *  p)

Creates a decay.

Definition at line 49 of file EvtBCL.cc.

50 {
51 p->initializePhaseSpace(getNDaug(), getDaugs());
52 calcamp->CalcAmp(p, _amp2, bclmodel);
53 }
EvtSemiLeptonicAmp * calcamp
Pointers needed to calculate amplitude.
Definition: EvtBCL.h:55
EvtSemiLeptonicFF * bclmodel
Pointers needed for FFs.
Definition: EvtBCL.h:52

◆ decay() [3/12]

void decay ( EvtParticle *  p)

The function evaluates the decay amplitude of the parent particle.

Parameters
pa pointer to the parent particle

Definition at line 53 of file EvtBSemiTauonic.cc.

54 {
55 p->initializePhaseSpace(getNDaug(), getDaugs());
56 m_CalcAmp->CalcAmp(p, _amp2, m_CalcHelAmp);
57 }
EvtBSemiTauonicAmplitude * m_CalcAmp
A pointer to the spin-dependent amplitude calculator specific to the spin type of the daughter meson.
EvtBSemiTauonicHelicityAmplitudeCalculator * m_CalcHelAmp
A pointer to the helicity amplitude calculator.

◆ decay() [4/12]

void decay ( EvtParticle *  p)

The function evaluates the decay amplitude of the parent particle.

a pointer to the parent particle

Definition at line 53 of file EvtBSemiTauonic2HDMType2.cc.

54 {
55 p->initializePhaseSpace(getNDaug(), getDaugs());
56 m_CalcAmp->CalcAmp(p, _amp2, m_CalcHelAmp);
57 }
EvtBSemiTauonicAmplitude * m_CalcAmp
A pointer to the spin-dependent amplitude calculator specific to the spin type of the daughter meson.
EvtBSemiTauonicHelicityAmplitudeCalculator * m_CalcHelAmp
A pointer to the helicity amplitude calculator.

◆ decay() [5/12]

void decay ( EvtParticle *  p)

The function to determine kinematics of daughter particles based on dGdsb distribution.

Definition at line 46 of file EvtBtoXsnunu_FERMI.cc.

47 {
48 p->makeDaughters(getNDaug(), getDaugs());
49
50 EvtParticle* xhadron = p->getDaug(0);
51 EvtParticle* leptonp = p->getDaug(1);
52 EvtParticle* leptonn = p->getDaug(2);
53
54 double mass[3];
55
56 findMasses(p, getNDaug(), getDaugs(), mass);
57
58 double mB = p->mass();
59 double ml = mass[1];
60 double pb(0.); // fermi momentum of b-quark
61
62 double xhadronMass = -999.0;
63
64 EvtVector4R p4xhadron;
65 EvtVector4R p4leptonp;
66 EvtVector4R p4leptonn;
67
68 while (xhadronMass < m_mxmin) {
69
70 // Apply Fermi motion and determine effective b-quark mass
71
72 double mb = 0.0;
73
74 bool FailToSetmb = true; // true when an appropriate mb cannot be found
75 while (FailToSetmb) {
76 pb = FermiMomentum(m_pf);
77
78 // effective b-quark mass
79 mb = mB * mB + m_mq * m_mq - 2.0 * mB * sqrt(pb * pb + m_mq * m_mq);
80 if (mb > 0. && sqrt(mb) - m_ms < 2.0 * ml) FailToSetmb = true;
81 else if (mb <= 0.0) FailToSetmb = true;
82 else FailToSetmb = false;
83 }
84 mb = sqrt(mb);
85
86 double mb_prob = m_mb_prob; // b-quark mass for probability density
87 double ms_prob = m_ms_prob; // s-quark mass for probability density
88 double mstilda = ms_prob / mb_prob;
89 double sb = 0.0;
90 double sbmin = 0;
91 double sbmax = (1 - mstilda) * (1 - mstilda);
92 while (sb == 0.0) {
93 double xbox = EvtRandom::Flat(sbmin, sbmax);
94 double ybox = EvtRandom::Flat(m_dGdsbProbMax);
95 double prob = dGdsbProb(xbox);
96 if (ybox < prob) sb = xbox;
97 }
98
99 // b->s (nu nubar)
100 EvtVector4R p4sdilep[2];
101
102 double msdilep[2];
103 msdilep[0] = m_ms;
104 msdilep[1] = sqrt(sb * mb_prob * mb_prob);
105
106 EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
107
108 // we do not care about neutrino
109 // just (nu nubar) -> nu nubar by phase space
110 EvtVector4R p4ll[2];
111
112 double mll[2];
113 mll[0] = ml;
114 mll[1] = ml;
115
116 EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
117
118 // boost to b-quark rest frame
119
120 p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
121 p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
122
123 // assign 4-momenta to valence quarks inside B meson in B rest frame
124 double phi = EvtRandom::Flat(EvtConst::twoPi);
125 double costh = EvtRandom::Flat(-1.0, 1.0);
126 double sinth = sqrt(1.0 - costh * costh);
127
128 // b-quark four-momentum in B meson rest frame
129
130 EvtVector4R p4b(sqrt(mb * mb + pb * pb),
131 pb * sinth * sin(phi),
132 pb * sinth * cos(phi),
133 pb * costh);
134
135 EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
136 p4leptonp = boostTo(p4ll[0], p4b);
137 p4leptonn = boostTo(p4ll[1], p4b);
138
139 // spectator quark in B meson rest frame
140 EvtVector4R p4q(sqrt(pb * pb + m_mq * m_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3));
141
142 // hadron system in B meson rest frame
143 p4xhadron = p4s + p4q;
144 xhadronMass = p4xhadron.mass();
145
146 }
147
148 // initialize the decay products
149 xhadron->init(getDaug(0), p4xhadron);
150
151 // assign the momentum of neutrino
152 // it is out of our interest
153 leptonp->init(getDaug(1), p4leptonp);
154 leptonn->init(getDaug(2), p4leptonn);
155
156 return;
157
158 }
double m_mq
mass of spectator quark for the Fermi motion model.
double m_dGdsbProbMax
The maximum value of dGdsb.
double m_mxmin
Minimum mass of Xs.
double m_mb_prob
b-quark mass to calculate dGdsb.
double m_ms
mass of s-quark for the Fermi motion model and EvtGenKine::PhaseSpace.
double m_pf
Parameter for the Fermi motioin model.
double m_ms_prob
s-quark mass to calculate dGdsb.
double dGdsbProb(double _sb)
The function returns the probability density value depending on sb.
double FermiMomentum(double pf)
The function returns a momentum of b quark.

◆ decay() [6/12]

void decay ( EvtParticle *  p)

Function that implements the energy-dependent Dalitz.

Definition at line 64 of file EvtEtaFullDalitz.cc.

65 {
66
67 const double a = getArg(0);
68 const double b = getArg(1);
69 const double c = getArg(2);
70 const double d = getArg(3);
71 const double e = getArg(4);
72 const double f = getArg(5);
73
74 p->initializePhaseSpace(getNDaug(), getDaugs());
75
76 EvtVector4R mompip = p->getDaug(0)->getP4();
77 EvtVector4R mompim = p->getDaug(1)->getP4();
78 EvtVector4R mompi0 = p->getDaug(2)->getP4();
79
80 double m_eta = p->mass();
81 double m_pip = p->getDaug(0)->mass();
82 double m_pi0 = p->getDaug(2)->mass();
83
84 //Q value
85 double deltaM = m_eta - 2 * m_pip - m_pi0;
86
87 //The decay amplitude comes from BESIII collab, Phys.Rev. D92 (2015) 012014
88
89 //Kinetic energies T
90 double Tpip = (mompip.get(0) - m_pip);
91 double Tpim = (mompim.get(0) - m_pip);
92 double Tpi0 = (mompi0.get(0) - m_pip);
93
94 //Dalitz variables
95 double X = sqrt(3.) * (Tpip - Tpim) / deltaM;
96 double Y = 3.*Tpi0 / deltaM - 1. ;
97
98 double amp2 = 1. + a * Y + b * Y * Y + c * X + d * X * X + e * X * Y + f * Y * Y * Y;
99
100 EvtComplex amp(sqrt(amp2), 0.0);
101
102 vertex(amp);
103
104 return ;
105
106 }

◆ decay() [7/12]

void decay ( EvtParticle *  p)

Function that implements the energy-dependent Dalitz.

Definition at line 79 of file EvtEtaPi0Dalitz.cc.

80 {
81
82 const double alpha = getArg(0);
83
84 p->initializePhaseSpace(getNDaug(), getDaugs());
85
86 EvtVector4R mompi0_0 = p->getDaug(0)->getP4();
87 EvtVector4R mompi0_1 = p->getDaug(1)->getP4();
88 EvtVector4R mompi0_2 = p->getDaug(2)->getP4();
89
90 double m_eta = p->mass();
91 double m_pi0 = p->getDaug(0)->mass();
92 double deltaM = m_eta - 3 * m_pi0;
93
94 //The decay amplitude comes from KLOE collab., Phys.Lett. B694 (2011) 16-21
95 double z0 = (3 * mompi0_0.get(0) - m_eta) / deltaM;
96 double z1 = (3 * mompi0_1.get(0) - m_eta) / deltaM;
97 double z2 = (3 * mompi0_2.get(0) - m_eta) / deltaM;
98
99 double z = (2. / 3.) * (z0 * z0 + z1 * z1 + z2 * z2) ;
100
101 double Amp2 = 1.0 + alpha * 2.0 * z;
102 if (Amp2 < 0)
103 Amp2 = 0;
104
105 EvtComplex amp(sqrt(Amp2), 0.);
106
107 vertex(amp);
108
109 return ;
110
111 }

◆ decay() [8/12]

void decay ( EvtParticle *  p)

Function that implements the energy-dependent Dalitz.

Definition at line 81 of file EvtEtaPrimeDalitz.cc.

82 {
83
84 const double a = getArg(0);
85 const double b = getArg(1);
86 const double c = getArg(2);
87 const double d = getArg(3);
88
89 p->initializePhaseSpace(getNDaug(), getDaugs());
90
91 EvtVector4R mompip = p->getDaug(0)->getP4();
92 EvtVector4R mompim = p->getDaug(1)->getP4();
93 EvtVector4R mometa = p->getDaug(2)->getP4();
94
95 double m_etaprime = p->mass();
96 double m_eta = p->getDaug(2)->mass();
97 double m_pi = p->getDaug(1)->mass();
98
99 //Q value
100 double deltaM = m_etaprime - 2 * m_pi - m_eta;
101
102 //The decay amplitude comes from BESIII collab, Phys.Rev. D92 (2015) 012014
103
104 //Kinetic energies T
105 double Tpip = (mompip.get(0) - m_pi);
106 double Tpim = (mompim.get(0) - m_pi);
107 double Teta = (mometa.get(0) - m_eta);
108
109 //Dalitz variables
110 double X = sqrt(3.) * (Tpip - Tpim) / deltaM;
111 double Y = ((m_eta + 2 * m_pi) / m_pi) * Teta / deltaM - 1. ;
112
113 double amp2 = 1. + a * Y + b * Y * Y + c * X + d * X * X;
114
115 EvtComplex amp(sqrt(amp2), 0.0);
116
117 vertex(amp);
118
119 return ;
120
121 }

◆ decay() [9/12]

void decay ( EvtParticle *  p)

The function to determine kinematics of daughter particles.

Definition at line 46 of file EvtFlatDaughter.cc.

47 {
48 // variables to find.
49 EvtVector4R p4KorKstar;
50 EvtVector4R p4norantin;
51 EvtVector4R p4antinorn;
52
53 // start!
54 p->makeDaughters(getNDaug(), getDaugs());
55
56 EvtParticle* KorKstar = p->getDaug(0);
57 EvtParticle* norantin = p->getDaug(1);
58 EvtParticle* antinorn = p->getDaug(2);
59
60 double mass[3];
61
62 findMasses(p, getNDaug(), getDaugs(), mass);
63
64 double m_B = p->mass();
65 double m_K = mass[0];
66 double m_norantin = mass[1];
67 double m_antinorn = mass[2];
68
69 double M_max = m_B - m_K;
70 double M_min = m_norantin + m_antinorn;
71
72 // determined M_nn
73 double Mnn = EvtRandom::Flat(M_min, M_max);
74
75 // B -> K(*) (nnbar) decay. Momentums are determined in B meson's rest frame
76 EvtVector4R p4Kdin[2];
77
78 double mKdin[2];
79 mKdin[0] = m_K;
80 mKdin[1] = Mnn;
81
82 EvtGenKine::PhaseSpace(2, mKdin, p4Kdin, m_B);
83
84 // (nnbar) -> n nbar decay. Momentums are determined in (nnbar) rest frame
85 EvtVector4R p4nn[2];
86
87 double mnn[2];
88 mnn[0] = m_norantin;
89 mnn[1] = m_antinorn;
90
91 EvtGenKine::PhaseSpace(2, mnn, p4nn, Mnn);
92
93 // boost to B meson rest frame
94 p4norantin = boostTo(p4nn[0], p4Kdin[1]);
95 p4antinorn = boostTo(p4nn[1], p4Kdin[1]);
96
97 // initialize the decay products
98 KorKstar->init(getDaug(0), p4Kdin[0]);
99 norantin->init(getDaug(1), p4norantin);
100 antinorn->init(getDaug(2), p4antinorn);
101
102
103 return;
104
105 }

◆ decay() [10/12]

void decay ( EvtParticle *  p)

The function to calculate a quark decay amplitude.

Definition at line 50 of file EvtKnunu.cc.

51 {
52 static EvtId NUE = EvtPDL::getId("nu_e");
53 static EvtId NUM = EvtPDL::getId("nu_mu");
54 static EvtId NUT = EvtPDL::getId("nu_tau");
55 static EvtId NUEB = EvtPDL::getId("anti-nu_e");
56 static EvtId NUMB = EvtPDL::getId("anti-nu_mu");
57 static EvtId NUTB = EvtPDL::getId("anti-nu_tau");
58
59 p->initializePhaseSpace(getNDaug(), getDaugs());
60
61 double m_b = p->mass();
62
63 EvtParticle* meson, * neutrino1, * neutrino2;
64 meson = p->getDaug(0);
65 neutrino1 = p->getDaug(1);
66 neutrino2 = p->getDaug(2);
67 EvtVector4R momnu1 = neutrino1->getP4();
68 EvtVector4R momnu2 = neutrino2->getP4();
69 EvtVector4R momk = meson->getP4();
70
71 EvtVector4R q = momnu1 + momnu2;
72 double q2 = q.mass2();
73
74 EvtVector4R p4b; p4b.set(m_b, 0., 0., 0.); // Do calcs in mother rest frame
75
76 double m_k = meson->mass();
77 double mkhat = m_k / (m_b);
78 double shat = q2 / (m_b * m_b);
79
80 EvtVector4R phat = (p4b + momk) / m_b;
81 EvtVector4R qhat = q / m_b;
82
83 // calculate form factor from [arXiv:1409.4557v2]
84 // see eq 104, eq 105, eq 106, eq 107
85 double alpha0 = 0.432;
86 double alpha1 = -0.664;
87 double alpha2 = -1.2;
88 double mp = m_b + 0.046;
89 double tp = (m_b + m_k) * (m_b + m_k);
90 double tm = (m_b - m_k) * (m_b - m_k);
91 double t0 = tp * (1 - sqrt(1 - tm / tp));
92 double z = (sqrt(tp - q2) - sqrt(tp - t0)) / (sqrt(tp - q2) + sqrt(tp - t0));
93 double fp = (1 / (1 - q2 / (mp * mp))) * (alpha0 + alpha1 * z + alpha2 * z * z + (-alpha1 + 2 * alpha2) * z * z * z / 3);
94 double fm = -fp * (1 - mkhat * mkhat) / shat;
95
96 // calculate quark decay amplitude from [arXiv:hep-ph/9910221v2]
97 // see eq 3.1, eq 3.2, eq 4.1, eq 4.2, and eq 4.3
98 // but in B->K nu nubar, fT and f0 terms are ignored
99
100 EvtVector4C T1;
101
102 EvtComplex aprime;
103 aprime = fp;
104 EvtComplex bprime;
105 bprime = fm;
106
107 T1 = aprime * phat + bprime * qhat;
108
109 EvtVector4C l;
110
111 if (getDaug(1) == NUE || getDaug(1) == NUM || getDaug(1) == NUT) {
112 l = EvtLeptonVACurrent(neutrino1->spParentNeutrino(),
113 neutrino2->spParentNeutrino());
114 }
115 if (getDaug(1) == NUEB || getDaug(1) == NUMB || getDaug(1) == NUTB) {
116 l = EvtLeptonVACurrent(neutrino2->spParentNeutrino(),
117 neutrino1->spParentNeutrino());
118 }
119
120 vertex(l * T1);
121
122 }

◆ decay() [11/12]

void decay ( EvtParticle *  p)

The function to calculate a quark decay amplitude.

Definition at line 49 of file EvtKstarnunu_REV.cc.

50 {
51
52 static EvtId NUE = EvtPDL::getId("nu_e");
53 static EvtId NUM = EvtPDL::getId("nu_mu");
54 static EvtId NUT = EvtPDL::getId("nu_tau");
55 static EvtId NUEB = EvtPDL::getId("anti-nu_e");
56 static EvtId NUMB = EvtPDL::getId("anti-nu_mu");
57 static EvtId NUTB = EvtPDL::getId("anti-nu_tau");
58
59 p->initializePhaseSpace(getNDaug(), getDaugs());
60
61 double m_b = p->mass();
62
63 EvtParticle* meson, * neutrino1, * neutrino2;
64 meson = p->getDaug(0);
65 neutrino1 = p->getDaug(1);
66 neutrino2 = p->getDaug(2);
67 EvtVector4R momnu1 = neutrino1->getP4();
68 EvtVector4R momnu2 = neutrino2->getP4();
69 EvtVector4R momkstar = meson->getP4();
70
71 EvtVector4R q = momnu1 + momnu2;
72 double q2 = q.mass2();
73
74 double m_k = meson->mass();
75
76 EvtVector4R p4b; p4b.set(m_b, 0., 0., 0.); // Do calcs in mother rest frame
77
78 // calculate form factors [arXiv:1503.05534v3]
79 // see Table 15, eq 15, eq 16, Table 3
80
81 // FFs
82 double tp = (m_b + m_k) * (m_b + m_k);
83 double tm = (m_b - m_k) * (m_b - m_k);
84 double t0 = tp * (1 - sqrt(1 - tm / tp));
85 double z = (sqrt(tp - q2) - sqrt(tp - t0)) / (sqrt(tp - q2) + sqrt(tp - t0));
86 double z0 = (sqrt(tp) - sqrt(tp - t0)) / (sqrt(tp) + sqrt(tp - t0));
87
88 // v0
89 double alpha0_v0 = 0.38;
90 double alpha1_v0 = -1.17;
91 double alpha2_v0 = 2.42;
92 double mR_v0 = 5.415;
93 double v0 = (1 / (1 - q2 / (mR_v0 * mR_v0))) * (alpha0_v0 + alpha1_v0 * (z - z0) + alpha2_v0 * (z - z0) * (z - z0));
94
95 // A1
96 double alpha0_A1 = 0.3;
97 double alpha1_A1 = 0.39;
98 double alpha2_A1 = 1.19;
99 double mR_A1 = 5.829;
100 double A1 = (1 / (1 - q2 / (mR_A1 * mR_A1))) * (alpha0_A1 + alpha1_A1 * (z - z0) + alpha2_A1 * (z - z0) * (z - z0));
101
102 // A12
103 double alpha0_A12 = 0.27;
104 double alpha1_A12 = 0.53;
105 double alpha2_A12 = 0.48;
106 double mR_A12 = 5.829;
107 double A12 = (1 / (1 - q2 / (mR_A12 * mR_A12))) * (alpha0_A12 + alpha1_A12 * (z - z0) + alpha2_A12 * (z - z0) * (z - z0));
108
109 // lambda
110 double lambda = (tp - q2) * (tm - q2);
111
112 // A2
113 double A2 = ((m_b + m_k) * (m_b + m_k) * (m_b * m_b - m_k * m_k - q2) * A1 - A12 * 16 * m_b * m_k * m_k * (m_b + m_k)) / lambda;
114
115 // calculate quark decay amplitude from [arXiv:hep-ph/9910221v2]
116 // see eq 3.3, eq 3.4, eq 3.5, eq 4.1, eq 4.4, and eq 4.5
117 // but in B->Kstar nu nubar, A3, A0, and all T terms are ignored
118
119 // definition of A12 can be found from [arXiv:1303.5794v2]
120
121 // definition of A3 can be found from [arXiv:1408.5614v1]
122
123 EvtTensor4C tds = (-2 * v0 / (m_b + m_k)) * dual(EvtGenFunctions::directProd(p4b, momkstar))
124 - EvtComplex(0.0, 1.0) *
125 ((m_b + m_k) * A1 * EvtTensor4C::g()
126 - (A2 / (m_b + m_k)) * EvtGenFunctions::directProd(p4b - momkstar, p4b + momkstar));
127
128 EvtVector4C l;
129
130 if (getDaug(1) == NUE || getDaug(1) == NUM || getDaug(1) == NUT) {
131 l = EvtLeptonVACurrent(neutrino1->spParentNeutrino(),
132 neutrino2->spParentNeutrino());
133 }
134 if (getDaug(1) == NUEB || getDaug(1) == NUMB || getDaug(1) == NUTB) {
135 l = EvtLeptonVACurrent(neutrino2->spParentNeutrino(),
136 neutrino1->spParentNeutrino());
137 }
138
139 EvtVector4C et0, et1, et2;
140 et0 = tds.cont1(meson->epsParent(0).conj());
141 et1 = tds.cont1(meson->epsParent(1).conj());
142 et2 = tds.cont1(meson->epsParent(2).conj());
143
144 vertex(0, l * et0);
145 vertex(1, l * et1);
146 vertex(2, l * et2);
147
148 return;
149
150 }

◆ decay() [12/12]

void decay ( EvtParticle *  p)

Member of particle in EvtGen.

Definition at line 51 of file EvtPhspCP.cc.

52 {
53 static EvtId B0 = EvtPDL::getId("B0");
54 static EvtId B0B = EvtPDL::getId("anti-B0");
55
56 double t;
57 EvtId other_b;
58
59 EvtCPUtil::getInstance()->OtherB(p, t, other_b, 0.5);
60
61 p->initializePhaseSpace(getNDaug(), getDaugs());
62
63 EvtComplex amp;
64
65 EvtComplex A, Abar;
66
67 A = EvtComplex(getArg(3) * cos(getArg(4)), getArg(3) * sin(getArg(4)));
68 Abar = EvtComplex(getArg(5) * cos(getArg(6)), getArg(5) * sin(getArg(6)));
69
70 // CP asymmetry
71 const double angle = getArg(0); // q/p = exp(-2i*angle). |q/p| = 1
72 const double dm = getArg(1); // Delta m_d
73 const double etaCP = getArg(2); // CP eigenvalue
74
75 if (other_b == B0B) {
76 amp = A * cos(dm * t / (2 * EvtConst::c)) +
77 EvtComplex(cos(-2.0 * angle), sin(-2.0 * angle)) *
78 etaCP * EvtComplex(0.0, 1.0) * Abar * sin(dm * t / (2 * EvtConst::c));
79
80 }
81 if (other_b == B0) {
82 amp = A * EvtComplex(cos(2.0 * angle), sin(2.0 * angle)) *
83 EvtComplex(0.0, 1.0) * sin(dm * t / (2 * EvtConst::c)) +
84 etaCP * Abar * cos(dm * t / (2 * EvtConst::c));
85 }
86
87 setProb(abs2(amp));
88
89 return;
90 }

◆ dGammadcostau()

double dGammadcostau ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau,
int  tauhel,
int  Dhel,
double  costau 
)

Function calculates the differential decay rate dGamma/dcostau, integrated for w.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame. Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 59 of file EvtBSemiTauonicDecayRateCalculator.cc.

61 {
62 m_BSTD = &BSTD;
63 TF1 f1("f1", this, &EvtBSemiTauonicDecayRateCalculator::EvaluateByW, BSTD.wmin(), BSTD.wmax(mtau, Dhel), 4,
64 "EvtBSemiTauonicDecayRateCalculator", "EvaluateByW");
65 f1.SetParameter(0, mtau);
66 f1.SetParameter(1, tauhel);
67 f1.SetParameter(2, Dhel);
68 f1.SetParameter(3, costau);
69 ROOT::Math::WrappedTF1 wf1(f1);
70// ROOT::Math::GSLIntegrator ig;
71 ROOT::Math::GaussIntegrator ig;
72 ig.SetFunction(wf1);
73 return ig.Integral(BSTD.wmin(), BSTD.wmax(mtau, Dhel));
74 }
const EvtBSemiTauonicHelicityAmplitudeCalculator * m_BSTD
temporal pointer to the helicity amplitude calculator for EvaluateBy* functions
double EvaluateByW(double *x, double *param)
Function used internally for numerical integration.

◆ dGammadw()

double dGammadw ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau,
int  tauhel,
int  Dhel,
double  w 
)

Function calculates the differential decay rate dGamma/dw, integrated for costau.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable. Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 41 of file EvtBSemiTauonicDecayRateCalculator.cc.

43 {
44 m_BSTD = &BSTD;
45 TF1 f1("f1", this, &EvtBSemiTauonicDecayRateCalculator::EvaluateByCostau, -1, 1, 4,
46 "EvtBSemiTauonicDecayRateCalculator", "EvaluateByCostau");
47 f1.SetParameter(0, mtau);
48 f1.SetParameter(1, tauhel);
49 f1.SetParameter(2, Dhel);
50 f1.SetParameter(3, w);
51 ROOT::Math::WrappedTF1 wf1(f1);
52// ROOT::Math::GSLIntegrator ig;
53 ROOT::Math::GaussIntegrator ig;
54 ig.SetFunction(wf1);
55 return ig.Integral(-1, 1);
56 }
double EvaluateByCostau(double *x, double *param)
Function used internally for numerical integration.

◆ dGammadwdcostau()

double dGammadwdcostau ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
)

Function calculates the differential decay rate dGamma/dw/dcostau.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame. Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 31 of file EvtBSemiTauonicDecayRateCalculator.cc.

33 {
34 EvtComplex temp = BSTD.helAmp(mtau, tauhel, Dhel, w, costau);
35 std::complex<double> amp(real(temp), imag(temp));
36
37 return std::norm(amp) * pf(BSTD, mtau, Dhel, w);
38 }
double pf(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int Dhel, double w)
Phase space factor, which is multiplied to the helicity amplitude to calculate the decay rate.

◆ dGdsbProb()

double dGdsbProb ( double  _sb)

The function returns the probability density value depending on sb.

Definition at line 291 of file EvtBtoXsnunu_FERMI.cc.

292 {
293 // dGdsb: arXiv:1509.06248v2
294 // see (eq.41)
295
296 double mb = m_mb_prob; // b-quark mass for probability density
297 double ms = m_ms_prob; // s-quark mass for probability density
298 double mstilda = ms / mb;
299
300 double sb = _sb;
301 double mstilda2 = mstilda * mstilda;
302
303 double lambda = 1 + mstilda2 * mstilda2 + sb * sb - 2 * (mstilda2 + sb + mstilda2 * sb);
304
305 double prob = sqrt(lambda) * (3 * sb * (1 + mstilda2 - sb) + lambda);
306
307 return prob;
308 }

◆ dot()

T dot ( GeneralVector< T >  a,
GeneralVector< T >  b 
)

dot product of two general vectors

Definition at line 163 of file beamHelpers.h.

164 {
165 return (a.el[0] * b.el[0] + a.el[1] * b.el[1] + a.el[2] * b.el[2]);
166 }

◆ dR3()

double dR3 ( double  w) const

HQET correction factor for the scalar form factor for B->D*taunu.

Parameters
wvelocity transfer variable.
Returns
calculated factor.

Definition at line 494 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

495 {
496 return 0.22 - 0.052 * (w - 1.) + 0.026 * (w - 1.) * (w - 1.);
497 }

◆ dS1()

double dS1 ( double  w) const

HQET correction factor for the scalar form factor for B->Dtaunu.

Parameters
wvelocity transfer variable.
Returns
calculated factor.

Definition at line 489 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

490 {
491 return -0.019 + 0.041 * (w - 1.) - 0.015 * (w - 1.) * (w - 1.);
492 }

◆ EvaluateBy2D()

double EvaluateBy2D ( double *  x,
double *  param 
)
private

Function used internally for numerical integration.

Definition at line 219 of file EvtBSemiTauonicDecayRateCalculator.cc.

220 {
221 // x={w,costau}
222 return dGammadwdcostau(*m_BSTD, param[0], (int)param[1], (int)param[2], x[0], x[1]);
223 }
double dGammadwdcostau(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel, double w, double costau)
Function calculates the differential decay rate dGamma/dw/dcostau.

◆ EvaluateByCostau()

double EvaluateByCostau ( double *  x,
double *  param 
)
private

Function used internally for numerical integration.

Definition at line 213 of file EvtBSemiTauonicDecayRateCalculator.cc.

214 {
215 // x=costau
216 return dGammadwdcostau(*m_BSTD, param[0], (int)param[1], (int)param[2], param[3], x[0]);
217 }

◆ EvaluateByW()

double EvaluateByW ( double *  x,
double *  param 
)
private

Function used internally for numerical integration.

Definition at line 207 of file EvtBSemiTauonicDecayRateCalculator.cc.

208 {
209 // x=w
210 return dGammadwdcostau(*m_BSTD, param[0], (int)param[1], (int)param[2], x[0], param[3]);
211 }

◆ EvtBCL()

EvtBCL ( )

Default constructor.

Definition at line 29 of file EvtBCL.cc.

29: bclmodel(nullptr), calcamp(nullptr) {}

◆ EvtBCLFF()

EvtBCLFF ( int  numarg,
double *  arglist 
)

constructor

Definition at line 23 of file EvtBCLFF.cc.

24 {
25 for (int i = 0; i < m_numBCLFFCoefficients; ++i) {
26 m_BCLFFCoefficients[i] = arglist[i];
27 }
28 }
double m_BCLFFCoefficients[19]
Parameters passed to the model; BCL expansion coefficients.
Definition: EvtBCLFF.h:78
int m_numBCLFFCoefficients
Total number of parameters passed to the model.
Definition: EvtBCLFF.h:76

◆ EvtBSemiTauonic()

The default constructor

Definition at line 32 of file EvtBSemiTauonic.cc.

32: m_CalcHelAmp(0), m_CalcAmp(0) {}

◆ EvtBSemiTauonic2HDMType2()

The default constructor

Definition at line 32 of file EvtBSemiTauonic2HDMType2.cc.

32: m_CalcHelAmp(0), m_CalcAmp(0) {}

◆ EvtBSemiTauonicHelicityAmplitudeCalculator() [1/2]

The default constructor.

Initializes with the default parameter values used by the aurhos of PRD87,034028.

Definition at line 22 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

22 :
23 m_CV1(0), m_CV2(0), m_CS1(0), m_CS2(0), m_CT(0)
24 {
25 m_mB = 5.279; // value used in PRD87,034028
26 m_mD = 1.866; // value used in PRD87,034028
27 m_mDst = 2.009; // value used in PRD87,034028
28
29 m_rho12 = 1.186; // HFAG end of year 2011
30 m_rhoA12 = 1.207; // HFAG end of year 2011
31
32 m_ffR11 = 1.403; // HFAG end of year 2011
33 m_ffR21 = 0.854; // HFAG end of year 2011
34
35 m_aS1 = 1.;
36 m_aR3 = 1.;
37
38 m_mBottom = 4.20; // PRD77,113016
39 m_mCharm = 0.901; // PRD77,113016 running mass at m_b scale
40
41 B2INFO("EvtBSemiTauonicHelicityAmplitudeCalculator initialized with the default values.");
42 B2INFO("rho_1^2 : " << m_rho12);
43 B2INFO("rho_A1^22 : " << m_rhoA12);
44 B2INFO("R_1(1) : " << m_ffR11);
45 B2INFO("R_2(1) : " << m_ffR21);
46 B2INFO("a_S1: " << m_aS1);
47 B2INFO("a_R3: " << m_aR3);
48 B2INFO("bottom quark mass: " << m_mBottom);
49 B2INFO("charm quark mass: " << m_mCharm);
50 B2INFO("CV1 : " << m_CV1);
51 B2INFO("CV2 : " << m_CV2);
52 B2INFO("CS1 : " << m_CS1);
53 B2INFO("CS2 : " << m_CS2);
54 B2INFO("CT : " << m_CT);
55 B2INFO("B meson mass: " << m_mB);
56 B2INFO("D meson mass: " << m_mD);
57 B2INFO("D* meson mass: " << m_mDst);
58 }
double m_mCharm
c quark mass (running mass at m_b scale), used for scalar form factor term )
double m_mBottom
b quark mass (running mass at m_b scale), used for scalar form factor term

◆ EvtBSemiTauonicHelicityAmplitudeCalculator() [2/2]

EvtBSemiTauonicHelicityAmplitudeCalculator ( const double  rho12,
const double  rhoA12,
const double  ffR11,
const double  ffR21,
const double  AS1,
const double  AR3,
const double  bottomMass,
const double  charmMass,
const EvtComplex &  CV1,
const EvtComplex &  CV2,
const EvtComplex &  CS1,
const EvtComplex &  CS2,
const EvtComplex &  CT,
const double  parentMass,
const double  DMass,
const double  DstarMass 
)

The constructor with HQET form factor parameters, Wilson coefficients of new physics contributions and parent B, daughter D(*) meson masses.

Parameters
rho12HQET form factor parameter rho_1^2 obtained by Dlnu decay data.
rhoA12HQET form factor parameter rho_A1^2 obtained by D*lnu decay data.
ffR11HQET form factor parameter R_1(1) obtained by D*lnu decay data.
ffR21HQET form factor parameter R_2(1) obtained by D*lnu decay data.
AS1a parameter to take into account the theoretical error of the scalar form factor for Dtaunu decay.
AR3a parameter to take into account the theoretical error of the scalar form factor for D*taunu decay.
CV1Wilson coeffcient of the left handed vector type NP contribution.
CV2Wilson coeffcient of the right handed vector type NP contribution.
CS1Wilson coeffcient of the scalar (S+P) type NP contribution.
CS2Wilson coeffcient of the scalar (S-P) type NP contribution.
CTWilson coeffcient of the tensor type NP contribution.
parentMassmass of the parent (B) meson.
DMassmass of the scalar type daughter (D) meson.
DstarMassmass of the vector type daughter (D*) meson.
bottomMassmass of the bottom quark mass (running mass at the energy of bottom quark mass)
charmMassmass of the charm quark mass (running mass at the energy of bottom quark mass) The constructor initializes the parameters of the decay. The recommended values of AS1 and AR3 by authors of PRD87,034028 are 1+/-1.

Definition at line 60 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

64 :
65 m_CV1(CV1), m_CV2(CV2), m_CS1(CS1), m_CS2(CS2), m_CT(CT)
66 {
67 m_mB = parentMass;
68 m_mD = DMass;
69 m_mDst = DstarMass;
70 m_mBottom = bottomMass;
71 m_mCharm = charmMass;
72
73 m_rho12 = rho12;
74 m_rhoA12 = rhoA12;
75
76 m_ffR11 = ffR11;
77 m_ffR21 = ffR21;
78
79 m_aS1 = aS1;
80 m_aR3 = aR3;
81
82 B2INFO("EvtBSemiTauonicHelicityAmplitudeCalculator initialized.");
83 B2INFO("rho_1^2 : " << rho12);
84 B2INFO("rho_A1^2 : " << rhoA12);
85 B2INFO("R_1(1) : " << ffR11);
86 B2INFO("R_2(1) : " << ffR21);
87 B2INFO("a_S1: " << m_aS1);
88 B2INFO("a_R3: " << m_aR3);
89 B2INFO("bottom quark mass: " << m_mBottom);
90 B2INFO("charm quark mass: " << m_mCharm);
91 B2INFO("CV1 : " << m_CV1);
92 B2INFO("CV2 : " << m_CV2);
93 B2INFO("CS1 : " << m_CS1);
94 B2INFO("CS2 : " << m_CS2);
95 B2INFO("CT : " << m_CT);
96 B2INFO("Parent meson mass: " << m_mB);
97 B2INFO("D meson mass: " << m_mD);
98 B2INFO("D* meson mass: " << m_mDst);
99 }
double aR3() const
HQET correction factor for the uncertainty of 1/m_Q correction.
double aS1() const
HQET correction factor for the uncertainty of 1/m_Q correction.

◆ FermiMomentum()

double FermiMomentum ( double  pf)

The function returns a momentum of b quark.

The distribution of the momentum is based on the Fermi motion model.

Definition at line 310 of file EvtBtoXsnunu_FERMI.cc.

311 {
312 // Pick a value for the b-quark Fermi motion momentum
313 // according to Ali's Gaussian model
314
315 // reference: Ali, Ahmed, et al. "Power corrections in the decay rate and distributions in B->Xs l+l- 2 in the standard model"
316 // see (eq.57)
317
318 double pb, pbmax;
319 pb = 0.0;
320 pbmax = 5.0 * pf;
321
322 while (pb == 0.0) {
323 double xbox = EvtRandom::Flat(pbmax);
324 double ybox = EvtRandom::Flat();
325 if (ybox < FermiMomentumProb(xbox, pf)) { pb = xbox; }
326 }
327
328 return pb;
329 }
double FermiMomentumProb(double pb, double pf)
The function returns a probability based on the Fermi motion model.

◆ FermiMomentumProb()

double FermiMomentumProb ( double  pb,
double  pf 
)

The function returns a probability based on the Fermi motion model.

reference: Ali, Ahmed, et al. "Power corrections in the decay rate and distributions in B->Xs l+l- 2 in the standard model" see (eq.57)

Parameters
pbmomentum of the b quark
pfmomentum width parameter in the Fermi motion model

Definition at line 331 of file EvtBtoXsnunu_FERMI.cc.

332 {
333 // Compute probability according to Ali's Gaussian model
334 // the function chosen has a convenient maximum value of 1 for pb = pf
335
336 // reference: Ali, Ahmed, et al. "Power corrections in the decay rate and distributions in B->Xs l+l- 2 in the standard model"
337 // see (eq.57)
338
339 double prsq = (pb * pb) / (pf * pf);
340 double prob = prsq * exp(1.0 - prsq);
341
342 return prob;
343 }

◆ ffA1()

double ffA1 ( double  w) const

CLN form factor A1.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 467 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

468 {
469 return ffA11() * (1 - 8 * m_rhoA12 * z(w)
470 + (53 * m_rhoA12 - 15) * z(w) * z(w)
471 - (231 * m_rhoA12 - 91) * z(w) * z(w) * z(w));
472 }
double ffA11() const
Form factor normalization factor for B->D*lnu.

◆ ffR1()

double ffR1 ( double  w) const

CLN form factor R1.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 474 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

475 {
476 return m_ffR11 - 0.12 * (w - 1) + 0.05 * (w - 1) * (w - 1);
477 }

◆ ffR2()

double ffR2 ( double  w) const

CLN form factor R2.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 479 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

480 {
481 return m_ffR21 + 0.11 * (w - 1) - 0.06 * (w - 1) * (w - 1);
482 }

◆ ffR3()

double ffR3 ( double  w) const

CLN form factor R3.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 484 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

485 {
486 return 1 + aR3() * dR3(w);
487 }
double dR3(double w) const
HQET correction factor for the scalar form factor for B->D*taunu.

◆ ffS1()

double ffS1 ( double  w) const

CLN form factor S1.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 462 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

463 {
464 return (1 + aS1() * dS1(w)) * ffV1(w);
465 }
double dS1(double w) const
HQET correction factor for the scalar form factor for B->Dtaunu.

◆ ffV1()

double ffV1 ( double  w) const

CLN form factor V1.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 455 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

456 {
457 return ffV11() * (1 - 8 * m_rho12 * z(w)
458 + (51 * m_rho12 - 10) * z(w) * z(w)
459 - (252 * m_rho12 - 84) * z(w) * z(w) * z(w));
460 }
double ffV11() const
Form factor normalization factor for B->Dlnu.

◆ Flatte()

EvtComplex Flatte ( const double &  m,
const double &  m0 
)

Constant Flatte.

Definition at line 617 of file EvtB0toKsKK.cc.

618 {
619 const double g_pipi = 0.165;
620 const double g_kk = 4.21 * g_pipi;
621
622 static EvtId pip = EvtPDL::getId("pi+");
623 static EvtId kp = EvtPDL::getId("K+");
624
625 const double s = m * m;
626 const double s0 = m0 * m0;
627
628 const EvtComplex rho_pipi = 2.0 * Flatte_k(s, EvtPDL::getMeanMass(pip)) / m;
629 const EvtComplex rho_kk = 2.0 * Flatte_k(s, EvtPDL::getMeanMass(kp)) / m;
630
631 const EvtComplex Gamma =
632 EvtComplex(0.0, 1.0) * ((g_pipi * rho_pipi) + (g_kk * rho_kk));
633
634 const EvtComplex Flatte = 1.0 / (s0 - s - Gamma);
635
636 return Flatte;
637 }
EvtComplex Flatte_k(const double &s, const double &m_h)
Constant Flatte_k.
Definition: EvtB0toKsKK.cc:605

◆ Flatte_k()

EvtComplex Flatte_k ( const double &  s,
const double &  m_h 
)

Constant Flatte_k.

Definition at line 605 of file EvtB0toKsKK.cc.

606 {
607 const double k2 = 1.0 - (4.0 * m_h * m_h / s);
608 EvtComplex k;
609 if (k2 < 0.0)
610 k = EvtComplex(0.0, sqrt(fabs(k2)));
611 else
612 k = sqrt(k2);
613
614 return k;
615 }

◆ Gamma()

double Gamma ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau,
int  tauhel,
int  Dhel 
)

Function calculates the helicity dependent decay rate Gamma, integrated for w and costau.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D. Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 77 of file EvtBSemiTauonicDecayRateCalculator.cc.

79 {
80 m_BSTD = &BSTD;
82 BSTD.wmin(), BSTD.wmax(mtau, Dhel), -1, 1, 3,
83 "EvtBSemiTauonicDecayRateCalculator", "EvaluateBy2D");
84 f1.SetParameter(0, mtau);
85 f1.SetParameter(1, tauhel);
86 f1.SetParameter(2, Dhel);
87 ROOT::Math::WrappedMultiTF1 wf1(f1);
88
89 // Create the Integrator
90 ROOT::Math::AdaptiveIntegratorMultiDim ig;
91 ig.SetFunction(wf1);
92 double xmin[] = {BSTD.wmin(), -1};
93 double xmax[] = {BSTD.wmax(mtau, Dhel), 1};
94 return ig.Integral(xmin, xmax);
95 }
double EvaluateBy2D(double *x, double *param)
Function used internally for numerical integration.

◆ GammaD()

double GammaD ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau 
)

Function calculates the decay rate Gamma for Dtaunu decay, integrated for w and costau and summed for helicities.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter lepton mass. Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 98 of file EvtBSemiTauonicDecayRateCalculator.cc.

99 {
100 double sum(0);
101 sum += Gamma(BSTD, mtau, -1, 2);
102 sum += Gamma(BSTD, mtau, +1, 2);
103 return sum;
104 }
double Gamma(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, double mtau, int tauhel, int Dhel)
Function calculates the helicity dependent decay rate Gamma, integrated for w and costau.

◆ GammaDstar()

double GammaDstar ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau 
)

Function calculates the differential decay rate Gamma for D*taunu decay, integrated for w and costau and summed for helicities.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter lepton mass. Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 106 of file EvtBSemiTauonicDecayRateCalculator.cc.

107 {
108 double sum(0);
109 for (int Dhel = -1; Dhel <= 1; Dhel++) {
110 sum += Gamma(BSTD, mtau, -1, Dhel);
111 sum += Gamma(BSTD, mtau, +1, Dhel);
112 }
113 return sum;
114 }

◆ GammaSMD()

double GammaSMD ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mlep = 0.0005110 
)

Function calculates the SM decay rate Gamma for Dlnu decay, integrated for w and costau and summed for helicities.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mlepdaughter lepton mass (default is the electron mass). Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 117 of file EvtBSemiTauonicDecayRateCalculator.cc.

118 {
119
120 EvtBSemiTauonicHelicityAmplitudeCalculator SM(BSTD);
121 SM.setCV1(0);
122 SM.setCV2(0);
123 SM.setCS1(0);
124 SM.setCS2(0);
125 SM.setCT(0);
126 double sum(0);
127 sum += Gamma(SM, mlep, -1, 2);
128 sum += Gamma(SM, mlep, +1, 2);
129 return sum;
130 }

◆ GammaSMDstar()

double GammaSMDstar ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mlep = 0.0005110 
)

Function calculates the SM decay rate Gamma for D*lnu decay, integrated for w and costau and summed for helicities.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mlepdaughter lepton mass (default is the electron mass). Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 132 of file EvtBSemiTauonicDecayRateCalculator.cc.

133 {
134 EvtBSemiTauonicHelicityAmplitudeCalculator SM(BSTD);
135 SM.setCV1(0);
136 SM.setCV2(0);
137 SM.setCS1(0);
138 SM.setCS2(0);
139 SM.setCT(0);
140 double sum(0);
141 for (int Dhel = -1; Dhel <= 1; Dhel++) {
142 sum += Gamma(SM, mlep, -1, Dhel);
143 sum += Gamma(SM, mlep, +1, Dhel);
144 }
145 return sum;
146 }

◆ generate() [1/2]

MCInitialParticles & generate ( )

Generate a new event.

Definition at line 140 of file InitialParticleGeneration.cc.

141 {
142 return generate(m_allowedFlags);
143 }
int m_allowedFlags
Allowed generation flags.
MCInitialParticles & generate()
Generate a new event.

◆ generate() [2/2]

MCInitialParticles & generate ( int  allowedFlags)
private

Generate a new event wit a particular set of allowed flags.

Parameters
[in]allowedFlagsAllowed flags.

Definition at line 145 of file InitialParticleGeneration.cc.

146 {
147 if (!m_event) {
148 m_event.create();
149 }
150 if (!m_beamParams.isValid()) {
151 B2FATAL("Cannot generate beam without valid BeamParameters");
152 }
153 if (m_beamParams.hasChanged()) {
157 }
158 m_event->setGenerationFlags(m_beamParams->getGenerationFlags() & allowedFlags);
159 ROOT::Math::PxPyPzEVector her = generateBeam(m_beamParams->getHER(), m_beamParams->getCovHER(), m_generateHER);
160 ROOT::Math::PxPyPzEVector ler = generateBeam(m_beamParams->getLER(), m_beamParams->getCovLER(), m_generateLER);
161 ROOT::Math::XYZVector vtx = generateVertex(m_beamParams->getVertex(), m_beamParams->getCovVertex(), m_generateVertex);
162 m_event->set(her, ler, vtx);
163 //Check if we want to go to CMS, if so boost both
164 if (m_beamParams->hasGenerationFlags(BeamParameters::c_generateCMS)) {
165 m_event->setGenerationFlags(0);
166 her = m_event->getLabToCMS() * her;
167 ler = m_event->getLabToCMS() * ler;
168 m_event->set(her, ler, vtx);
169 m_event->setGenerationFlags(m_beamParams->getGenerationFlags() & allowedFlags);
170 }
171 return *m_event;
172 }
MultivariateNormalGenerator m_generateHER
Generator for HER.
StoreObjPtr< MCInitialParticles > m_event
Datastore object containing the generated event.
MultivariateNormalGenerator m_generateLER
Generator for LER.
MultivariateNormalGenerator m_generateVertex
Generator for Vertex.
@ c_generateCMS
generate initial event in CMS instead of lab
void reset()
reset the generator setting the size to 0.
ROOT::Math::XYZVector generateVertex(const ROOT::Math::XYZVector &initial, const TMatrixDSym &cov, MultivariateNormalGenerator &gen) const
generate the vertex
ROOT::Math::PxPyPzEVector generateBeam(const ROOT::Math::PxPyPzEVector &initial, const TMatrixDSym &cov, MultivariateNormalGenerator &gen) const
generate 4 vector for one beam

◆ generateBeam()

ROOT::Math::PxPyPzEVector generateBeam ( const ROOT::Math::PxPyPzEVector &  initial,
const TMatrixDSym &  cov,
MultivariateNormalGenerator gen 
) const
private

generate 4 vector for one beam

Parameters
initialbeam
covcovariance of the beam momentum (E, theta_x, theta_y)
genmultivariate normal generator to be used

Definition at line 101 of file InitialParticleGeneration.cc.

103 {
104 //no smearing, nothing to do
105 if (!(m_event->getGenerationFlags() & BeamParameters::c_smearBeam))
106 return initial;
107
108
109 //Calculate initial parameters of the beam before smearing
110 double E0 = initial.E();
111 double thX0 = std::atan(initial.Px() / initial.Pz());
112 double thY0 = std::atan(initial.Py() / initial.Pz());
113
114 //init random generator if there is a change in beam parameters or at the beginning
115 if (gen.size() != 3) gen.setMeanCov(ROOT::Math::XYZVector(E0, thX0, thY0), cov);
116
117 //generate the actual smeared vector (E, thX, thY)
118 Eigen::VectorXd p = gen.generate();
119
120 //if we don't smear beam energy set the E to initial value
121 if (!m_event->hasGenerationFlags(BeamParameters::c_smearBeamEnergy))
122 p[0] = E0;
123
124 //if we don't smear beam direction set thX and thY to initial values
125 if (!m_event->hasGenerationFlags(BeamParameters::c_smearBeamDirection)) {
126 p[1] = thX0;
127 p[2] = thY0;
128 }
129
130 //store smeared values
131 double E = p[0];
132 double thX = p[1];
133 double thY = p[2];
134
135 //Convert values back to 4-vector
136 bool isHER = initial.Z() > 0;
137 return BeamParameters::getFourVector(E, thX, thY, isHER);
138 }
static ROOT::Math::PxPyPzEVector getFourVector(double energy, double angleX, double angleY, bool isHER)
Calculate FourVector of a beam from energy and angles in xz and yz planes.

◆ generateEvent()

void generateEvent ( MCParticleGraph mcGraph)

Generates one single event.

Parameters
mcGraphReference to the MonteCarlo graph into which the generated particles will be stored.

Definition at line 105 of file CRY.cc.

106 {
107 bool eventInAcceptance = 0;
108 static std::vector<CRYParticle*> ev;
109
110 //force at least particle in acceptance box
111 for (int iTrial = 0; iTrial < m_maxTrials; ++iTrial) {
113 // Generate one event
114 ev.clear();
115 m_cryGenerator->genEvent(&ev);
116 // check all particles
117 for (auto* p : ev) {
118 const int pdg = p->PDGid();
119 const double kineticEnergy = p->ke() * Unit::MeV;
120 if (kineticEnergy < m_kineticEnergyThreshold) continue;
121
122 const double mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
123 const double etot = kineticEnergy + mass;
124 // since etot is at least mass this cannot be negative
125 const double ptot = sqrt(etot * etot - mass * mass);
126
127 // Momentum
128 // We have x horizontal, y up and z along the beam. So uvw -> zxy, xc yc zc -> zxy
129 const double px = ptot * p->v();
130 const double py = ptot * p->w();
131 const double pz = ptot * p->u();
132 // Vertex
133 const double vx = p->y() * Unit::m;
134 const double vy = p->z() * Unit::m;
135 const double vz = p->x() * Unit::m;
136
137 // Time
138 /* In basf2, it is assumed that t = 0 when an event was produced,
139 For the cosmic case, we set t = 0 when particle cross y=0 plane;
140 The output time from CRY (p->t()) is too large (order of second) and also
141 increase as simulated time, so it is impossible to handle in basf2.
142 For event which has more then one particle, the difference between their production
143 times is also too large (> micro-second) to keep in basf2, so the time relation
144 between particles in each event is also reset. Production of each particle in event is set t=0 at y=0.
145 if one need production from CRY for a special study, you have to find a way to handle it...*/
146 double time = 0;
147
148 vecgeom::Vector3D<double> pos(vx, vy, vz);
149 vecgeom::LorentzVector<double> mom(px, py, pz, etot);
150 const double speed = (mass == 0 ? 1 : mom.Beta()) * Const::speedOfLight;
151
152 // Project on the boundary of the world box
153 auto inside = m_world->Inside(pos);
154 if (inside == vecgeom::kInside) {
155 // inside the world volume, go backwards in time to the world box
156 const auto dir = -mom.vect().Unit();
157 double dist = m_world->DistanceToOut(pos, dir);
158 pos += dist * dir;
159 time -= dist / speed;
160 } else if (inside == vecgeom::kOutside) {
161 // outside the world volume, go forwards in time to the world box
162 // this should not happen but better safe then sorry ...
163 const auto dir = mom.vect().Unit();
164 double dist = m_world->DistanceToIn(pos, dir);
165 if (dist == vecgeom::InfinityLength<double>()) continue;
166 pos += dist * dir;
167 time += dist / speed;
168 }
169 // Intersect with the acceptance box
170 double dist = m_acceptance->DistanceToIn(pos, mom.vect().Unit());
171 if (dist == vecgeom::InfinityLength<double>()) continue;
172
173 // We want to keep this one
174 auto& particle = mcGraph.addParticle();
175 // all particle of a generator are primary
176 particle.addStatus(MCParticle::c_PrimaryParticle);
177 // all particle of CRY are stable
178 particle.addStatus(MCParticle::c_StableInGenerator);
179 particle.setPDG(pdg);
180 particle.setFirstDaughter(0);
181 particle.setLastDaughter(0);
182 particle.setMomentum(ROOT::Math::XYZVector(mom.x(), mom.y(), mom.z()));
183 particle.setMass(mass);
184 particle.setEnergy(mom.e());
185 particle.setProductionVertex(ROOT::Math::XYZVector(pos.x(), pos.y(), pos.z()));
186 particle.setProductionTime(time);
187 eventInAcceptance = true;
188 }
189 // clean up CRY event
190 for (auto* p : ev) delete p;
191 // and if we have something in the acceptance then we're done
192 if (eventInAcceptance) {
193 return;
194 }
195
196 }
197 B2ERROR("Number of trials exceeds maxTrials increase number of maxTrials" << LogVar("maxTrials", m_maxTrials));
198 }
std::unique_ptr< vecgeom::VUnplacedVolume > m_acceptance
acceptance shape
Definition: CRY.h:160
std::unique_ptr< vecgeom::VUnplacedVolume > m_world
world box shape
Definition: CRY.h:159
double m_kineticEnergyThreshold
kinetic energy threshold.
Definition: CRY.h:145
std::unique_ptr< CRYGenerator > m_cryGenerator
The CRY generator.
Definition: CRY.h:158
int m_totalTrials
total number of thrown events.
Definition: CRY.h:148
int m_maxTrials
number of trials per event.
Definition: CRY.h:147
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition: MCParticle.h:49
static const double m
[meters]
Definition: Unit.h:69
static const double MeV
[megaelectronvolt]
Definition: Unit.h:114
Class to store variables with their name which were sent to the logging service.

◆ generateVertex()

ROOT::Math::XYZVector generateVertex ( const ROOT::Math::XYZVector &  initial,
const TMatrixDSym &  cov,
MultivariateNormalGenerator gen 
) const
private

generate the vertex

Parameters
initialnominal vertex position
covcovariance of the vertex position
genmultivariate normal generator to be used

Definition at line 27 of file InitialParticleGeneration.cc.

29 {
30 if (m_event->hasGenerationFlags(BeamParameters::c_smearVertex)) {
31 if (!gen.size()) gen.setMeanCov(initial, cov);
32 return gen.generateVec3();
33 }
34 return initial;
35 }

◆ getAnglesCMS()

std::pair< T, T > getAnglesCMS ( GeneralVector< T >  p1,
GeneralVector< T >  p2 
)

get collision axis angles in CM frame from HER & LER momentum 3-vectors

Definition at line 204 of file beamHelpers.h.

205 {
206 GeneralVector<T> bv = getBoost(p1, p2);
207
208 T gamma = 1.0 / sqrt(1 - bv.norm2());
209
210 T E1 = sqrt(p1.norm2() + me * me);
211
212 GeneralVector<T> pCMS = p1 + ((gamma - 1) / bv.norm2() * dot(p1, bv) - gamma * E1) * bv;
213
214 return std::make_pair(pCMS.angleX(), pCMS.angleY());
215 }

◆ getbaryonff()

void getbaryonff ( EvtId  ,
EvtId  ,
double  ,
double  ,
double *  ,
double *  ,
double *  ,
double *   
)

Not Implemented.

Definition at line 139 of file EvtBCLFF.cc.

140 {
141 B2FATAL("Not implemented :getbaryonff in EvtBCLFF.");
142 }

◆ getBoost()

GeneralVector< T > getBoost ( GeneralVector< T >  p1,
GeneralVector< T >  p2 
)

get boost vector from HER & LER momentum 3-vectors

Definition at line 193 of file beamHelpers.h.

194 {
195 T E1 = sqrt(p1.norm2() + me * me);
196 T E2 = sqrt(p2.norm2() + me * me);
197
198 return 1. / (E1 + E2) * (p1 + p2);
199 }

◆ getCentralValues()

Eigen::VectorXd getCentralValues ( double  eH,
double  thXH,
double  thYH,
double  eL,
double  thXL,
double  thYL 
)
inline

get vector including (Ecms, boostX, boostY, boostZ, angleXZ, angleYZ) from beam energies and angles

Definition at line 267 of file beamHelpers.h.

269 {
270 Eigen::VectorXd mu(6);
271
272 GeneralVector<double> pH = getFourVector<double>(eH, thXH, thYH, true);
273 GeneralVector<double> pL = getFourVector<double>(eL, thXL, thYL, false);
274
275 double Ecms = getEcms<double>(pH, pL);
276 GeneralVector<double> boost = getBoost<double>(pH, pL);
277
278 double angleX, angleY;
279 std::tie(angleX, angleY) = getAnglesCMS<double>(pH, pL);
280
281 mu(0) = Ecms;
282 mu(1) = boost.el[0];
283 mu(2) = boost.el[1];
284 mu(3) = boost.el[2];
285 mu(4) = angleX;
286 mu(5) = angleY;
287
288 return mu;
289 }

◆ getdiracff()

void getdiracff ( EvtId  ,
EvtId  ,
double  ,
double  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *   
)

Not Implemented.

Definition at line 144 of file EvtBCLFF.cc.

145 {
146 B2FATAL("Not implemented :getbaryonff in EvtBCLFF.");
147 }

◆ getEcms()

T getEcms ( GeneralVector< T >  p1,
GeneralVector< T >  p2 
)

get CMS energy from HER & LER momentum 3-vectors

Definition at line 182 of file beamHelpers.h.

183 {
184 T E1 = sqrt(p1.norm2() + me * me);
185 T E2 = sqrt(p2.norm2() + me * me);
186
187 return sqrt((E1 + E2) * (E1 + E2) - (p1 + p2).norm2());
188 }

◆ getFourVector()

GeneralVector< T > getFourVector ( energy,
angleX,
angleY,
bool  isHER 
)

get 4-momentum from energy and angles of beam

Definition at line 220 of file beamHelpers.h.

221 {
222 T p = sqrt(energy * energy - me * me);
223
224 double dir = isHER ? 1 : -1;
225
226 T pz = dir * p / sqrt(tan(angleX) * tan(angleX) + tan(angleY) * tan(angleY) + 1.0);
227
228 return GeneralVector<T>(pz * tan(angleX), pz * tan(angleY), pz);
229 }

◆ getGradientMatrix()

Eigen::MatrixXd getGradientMatrix ( double  eH,
double  thXH,
double  thYH,
double  eL,
double  thXL,
double  thYL 
)
inline

get gradient matrix of (eH, thXH, thYH, eL, thXL, thYL) -> (eCMS, boostX, boostY, boostY, angleXZ, angleYZ) transformation

Definition at line 233 of file beamHelpers.h.

235 {
236 Eigen::MatrixXd grad(6, 6);
237
238 //calculate derivatives wrt all 6 input variables
239 for (int j = 0; j < 6; ++j) {
240
241 std::vector<double> eps(6, 0.0);
242 eps[j] = 1.0;
243
244 GeneralVector<DualNumber> pH = getFourVector(DualNumber(eH, eps[0]), DualNumber(thXH, eps[1]), DualNumber(thYH, eps[2]), true);
245 GeneralVector<DualNumber> pL = getFourVector(DualNumber(eL, eps[3]), DualNumber(thXL, eps[4]), DualNumber(thYL, eps[5]), false);
246
247
248 DualNumber Ecms = getEcms(pH, pL);
249 GeneralVector<DualNumber> boost = getBoost(pH, pL);
250
251 DualNumber angleX, angleY;
252 std::tie(angleX, angleY) = getAnglesCMS(pH, pL);
253
254 grad(0, j) = Ecms.dx;
255 grad(1, j) = boost.el[0].dx;
256 grad(2, j) = boost.el[1].dx;
257 grad(3, j) = boost.el[2].dx;
258 grad(4, j) = angleX.dx;
259 grad(5, j) = angleY.dx;
260 }
261
262 return grad;
263 }

◆ getInstance()

EvtGenModelRegister & getInstance ( )
staticprivate

Return reference to the instance.

This class behaves like a purely static class but we need a singleton pattern to avoid initialisation hell.

Definition at line 21 of file EvtGenModelRegister.cc.

22 {
23 static unique_ptr<EvtGenModelRegister> instance(new EvtGenModelRegister());
24 return *instance;
25 }
EvtGenModelRegister()
Singleton: private constructor.

◆ getModels()

list< EvtDecayBase * > getModels ( )
static

Return a list of models.

Caller takes responsibility to free the instances when no longer needed

Returns
list with pointers to instances of all registered models.

Definition at line 27 of file EvtGenModelRegister.cc.

28 {
29 list<EvtDecayBase*> modelList;
30 for (auto factory : getInstance().m_models) {
31 modelList.push_back(factory());
32 }
33 return modelList;
34 }
std::vector< ModelFactory * > m_models
List of all registered EvtGenModels.
static EvtGenModelRegister & getInstance()
Return reference to the instance.

◆ getName() [1/12]

std::string getName ( )

Get function Name

Definition at line 31 of file EvtB0toKsKK.cc.

32 {
33 return "B0toKsKK";
34 }

◆ getName() [2/12]

std::string getName ( )

Returns name of module.

Definition at line 39 of file EvtBCL.cc.

40 {
41 return "BCL";
42 }

◆ getName() [3/12]

std::string getName ( )

The function returns the model name.

Returns
name of the model "BSTD"

Definition at line 40 of file EvtBSemiTauonic.cc.

41 {
42 return "BSTD";
43 }

◆ getName() [4/12]

std::string getName ( )

The function returns the model name.

Returns
name of the model "BSTD_2HDMTYPE2"

Definition at line 40 of file EvtBSemiTauonic2HDMType2.cc.

41 {
42 return "BSTD_2HDMTYPE2";
43 }

◆ getName() [5/12]

std::string getName ( )

The function which returns the name of the model.

Definition at line 36 of file EvtBtoXsnunu_FERMI.cc.

37 {
38 return "BTOXSNUNU_FERMI";
39 }

◆ getName() [6/12]

std::string getName ( )

Returns the model name: ETA_FULLDALITZ.

Definition at line 26 of file EvtEtaFullDalitz.cc.

27 {
28
29 return "ETA_FULLDALITZ";
30
31 }

◆ getName() [7/12]

std::string getName ( )

Returns the name of the model: ETA_PI0DALITZ.

Definition at line 41 of file EvtEtaPi0Dalitz.cc.

42 {
43
44 return "ETA_PI0DALITZ";
45
46 }

◆ getName() [8/12]

std::string getName ( )

Returns the model name: ETAPRIME_DALITZ.

Definition at line 43 of file EvtEtaPrimeDalitz.cc.

44 {
45
46 return "ETAPRIME_DALITZ";
47
48 }

◆ getName() [9/12]

std::string getName ( )

The function which returns the name of the model.

Definition at line 36 of file EvtFlatDaughter.cc.

37 {
38 return "FLAT_DAUGHTER";
39 }

◆ getName() [10/12]

std::string getName ( )

The function which returns the name of the model.

Definition at line 40 of file EvtKnunu.cc.

41 {
42 return "KNUNU";
43 }

◆ getName() [11/12]

std::string getName ( )

The function which returns the name of the model.

Definition at line 39 of file EvtKstarnunu_REV.cc.

40 {
41 return "KSTARNUNU_REV";
42 }

◆ getName() [12/12]

std::string getName ( )

Get function Name

Definition at line 30 of file EvtPhspCP.cc.

31 {
32 return "PHSP_CP";
33 }

◆ getraritaff()

void getraritaff ( EvtId  ,
EvtId  ,
double  ,
double  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *   
)

Not Implemented.

Definition at line 149 of file EvtBCLFF.cc.

150 {
151 B2FATAL("Not implemented :getbaryonff in EvtBCLFF.");
152 }

◆ getscalarff()

void getscalarff ( EvtId  parent,
EvtId  daughter,
double  t,
double  ,
double *  fpf,
double *  f0f 
)

Scalar FF's.

Implementation follows arXiv:1509.06938v3.

For the scalar FF, the arglist in the constructor should contain 8 expansion parameters: b+_0, b+_1, b+_2, b+_3, b0_0, b0_1, b0_2, b0_3

Parameters
parent
daughter
tMomentum transfer, also called q2. q2 = (p_B - p_M)^2
fpff_+(q2).
f0ff_0(q2).

Definition at line 30 of file EvtBCLFF.cc.

31 {
32
33 if (m_numBCLFFCoefficients != 8) {
34 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Wrong number of arguments for EvtBCLFF::getscalarff!\n";
35 }
36
37 const auto mB = EvtPDL::getMeanMass(parent);
38 const auto mM = EvtPDL::getMeanMass(daughter);
39
40 const auto tplus = (mB + mM) * (mB + mM);
41 const auto tzero = (mB + mM) * (std::sqrt(mB) - std::sqrt(mM)) * (std::sqrt(mB) - std::sqrt(mM));
42
43 const auto mR2 = m_resonance1Minus * m_resonance1Minus;
44 const auto pole = 1 / (1 - t / mR2);
45
46 const std::vector<double> bplus = {m_BCLFFCoefficients[0], m_BCLFFCoefficients[1], m_BCLFFCoefficients[2], m_BCLFFCoefficients[3]};
47 const std::vector<double> bzero = {m_BCLFFCoefficients[4], m_BCLFFCoefficients[5], m_BCLFFCoefficients[6], m_BCLFFCoefficients[7]};
48
49 const auto N_fpf = bplus.size();
50 const auto N_f0f = bzero.size();
51
52 auto z = [tplus, tzero](decltype(t) q2) {
53 const auto term1 = std::sqrt(tplus - q2);
54 const auto term2 = std::sqrt(tplus - tzero);
55 return (term1 - term2) / (term1 + term2);
56 };
57
58 double sum_fpf = 0;
59 for (unsigned int n = 0; n < N_fpf; ++n) {
60 sum_fpf += bplus[n] * (std::pow(z(t), n) - std::pow(-1, n - N_fpf) * n / N_fpf * std::pow(z(t), N_fpf));
61 }
62 *fpf = pole * sum_fpf;
63
64 double sum_f0f = 0;
65 for (unsigned int n = 0; n < N_f0f; ++n) {
66 sum_f0f += bzero[n] * std::pow(z(t), n);
67 }
68 *f0f = sum_f0f;
69 }
static constexpr double m_resonance1Minus
Mass of the 1- resonance for the parametrization of the vector FF.
Definition: EvtBCLFF.h:83

◆ getScaleFactor()

double getScaleFactor ( const std::vector< ROOT::Math::PxPyPzMVector > &  ps,
double  EcmsTarget 
)
inline

Find approximative value of the momentum scaling factor which makes the total energy of the particles provided in the vector ps equal to EcmsTarget.

Definition at line 44 of file scaleParticleEnergies.h.

45 {
46 double s = 0;
47 for (auto p : ps)
48 s += pow(p.P(), 2) / (2 * p.E());
49
50 double dE = getTotalEnergy(ps) - EcmsTarget;
51 double L = dE / s;
52
53 return sqrt(1 - L);
54 }

◆ gettensorff()

void gettensorff ( EvtId  parent,
EvtId  daughter,
double  t,
double  ,
double *  hf,
double *  kf,
double *  bp,
double *  bm 
)

Not Implemented.

Definition at line 134 of file EvtBCLFF.cc.

135 {
136 B2FATAL("Not implemented :gettensorff in EvtBCLFF.");
137 }

◆ getTotalEnergy()

double getTotalEnergy ( const std::vector< ROOT::Math::PxPyPzMVector > &  ps)
inline

Main function scaleParticleEnergies in this header scales momenta in CMS of all particles in the final state by a constant factor such that the overall collision energy is slightly changed.

It is used for MC generators which do not provide option to generate events with smeared energy of initial particles. There is an assumption that matrix element of the process does not change much when Ecms is varied by 5MeV/10580MeV = 0.05% (typical Ecms smearing at Belle II). This is typically the case if the cross section does not have the resonance peak around collision energy. Get total energy of all particles in the provided vector

Definition at line 34 of file scaleParticleEnergies.h.

35 {
36 double eTot = 0;
37 for (auto p : ps)
38 eTot += p.E();
39 return eTot;
40 }

◆ getvectorff()

void getvectorff ( EvtId  parent,
EvtId  daughter,
double  t,
double  ,
double *  a1f,
double *  a2f,
double *  vf,
double *  a0f 
)

Vector FF's.

Implementation follows arXiv:1503.05534v3. It is assumed that each expansion has three terms (hardcoded). However, this can be easily expanded or generalized. It is not done, because this way we can check if the number of arguments in the decay file is the correct one.

For the vector FF, the arglist in the constructor should contain 11 expansion parameters: A0_1, A0_2, A1_0, A1_1, A1_2, A12_0, A12_1, A12_2, V_0, V_1, V_2 Nota bene: A0_0 is correlated to A12_0.

Parameters
parent
daughter
tMomentum transfer, also called q2. q2 = (p_B - p_M)^2
a1fA1(q2)
a2fA2(q2)
vfV(q2)
a0fA0(q2)

Definition at line 71 of file EvtBCLFF.cc.

72 {
73
74 if (m_numBCLFFCoefficients != 11) {
75 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Wrong number of arguments for EvtBCLFF::getvectorff!\n";
76 }
77
78 const auto mB = EvtPDL::getMeanMass(parent);
79 const auto mB2 = mB * mB;
80 const auto mM = EvtPDL::getMeanMass(daughter);
81 const auto mM2 = mM * mM;
82
83 const auto tplus = (mB + mM) * (mB + mM);
84 const auto tminus = (mB - mM) * (mB - mM);
85 const auto tzero = tplus * (1 - std::sqrt(1 - tminus / tplus));
86
87 const auto mR2A0 = m_resonance0Minus * m_resonance0Minus;
88 const auto mR2A1 = m_resonance1Plus * m_resonance1Plus;
89 const auto mR2A12 = m_resonance1Plus * m_resonance1Plus;
90 const auto mR2V = m_resonance1Minus * m_resonance1Minus;
91
92 const auto poleA0 = 1. / (1 - t / mR2A0);
93 const auto poleA1 = 1. / (1 - t / mR2A1);
94 const auto poleA12 = 1. / (1 - t / mR2A12);
95 const auto poleV = 1. / (1 - t / mR2V);
96
97 const std::vector<double> A0 = {8 * mB* mM / (mB2 - mM2)* m_BCLFFCoefficients[5], m_BCLFFCoefficients[0], m_BCLFFCoefficients[1]};
98 const std::vector<double> A1 = {m_BCLFFCoefficients[2], m_BCLFFCoefficients[3], m_BCLFFCoefficients[4]};
99 const std::vector<double> A12 = {m_BCLFFCoefficients[5], m_BCLFFCoefficients[6], m_BCLFFCoefficients[7]};
100 const std::vector<double> V = {m_BCLFFCoefficients[8], m_BCLFFCoefficients[9], m_BCLFFCoefficients[10]};
101
102 auto z = [tplus, tzero](decltype(t) q2) {
103 const auto term1 = std::sqrt(tplus - q2);
104 const auto term2 = std::sqrt(tplus - tzero);
105 return (term1 - term2) / (term1 + term2);
106 };
107
108 auto sum = [&z](decltype(t) q2, std::vector<double> par) {
109 double tot = 0;
110 for (unsigned int n = 0; n < par.size(); ++n) {
111 tot += par[n] * std::pow(z(q2) - z(0), n);
112 }
113 return tot;
114 };
115
116 auto kaellen = [mB, mM](decltype(t) q2) {
117 return ((mB + mM) * (mB + mM) - q2) * ((mB - mM) * (mB - mM) - q2);
118 };
119
120 const auto ffA0 = poleA0 * sum(t, A0);
121 const auto ffA1 = poleA1 * sum(t, A1);
122 const auto ffA12 = poleA12 * sum(t, A12);
123 const auto ffV = poleV * sum(t, V);
124
125 const auto ffA2 = ((mB + mM) * (mB + mM) * (mB2 - mM2 - t) * ffA1 - (16 * mB * mM2 * (mB + mM)) * ffA12) / kaellen(t);
126
127 *a0f = ffA0;
128 *a1f = ffA1;
129 *a2f = ffA2;
130 *vf = ffV;
131
132 }
static constexpr double m_resonance0Minus
Mass of the 0- resonance for the parametrization of the vector FF.
Definition: EvtBCLFF.h:81
static constexpr double m_resonance1Plus
Mass of the 1+ resonance for the parametrization of the vector FF.
Definition: EvtBCLFF.h:85

◆ getVertexConditional()

ROOT::Math::XYZVector getVertexConditional ( )

Generate vertex position and possibly update the generator of Lorentz transformation.

Definition at line 176 of file InitialParticleGeneration.cc.

177 {
178 if (!m_event) {
179 m_event.create();
180 }
181 if (!m_beamParams.isValid()) {
182 B2FATAL("Cannot generate beam without valid BeamParameters");
183 }
184 if (m_beamParams.hasChanged()) {
186 m_beamParams->getCovHER(), m_beamParams->getCovLER());
188 }
189 m_event->setGenerationFlags(m_beamParams->getGenerationFlags() & m_allowedFlags);
190 ROOT::Math::XYZVector vtx = generateVertex(m_beamParams->getVertex(), m_beamParams->getCovVertex(), m_generateVertex);
191
192 //TODO is m_event of type MCInitialParticles important?
193 //Eigen::VectorXd collSys = m_generateLorentzTransformation.generate(Ecms);
194 //m_event->set(her, ler, vtx);
195 //m_event->setByLorentzTransformation(collSys[0], collSys[1], collSys[2], collSys[3], collSys[4], collSys[5], vtx);
196
197 return vtx;
198 }
ConditionalGaussGenerator m_generateLorentzTransformation
Generator of the Lorentz transformation.
ConditionalGaussGenerator initConditionalGenerator(const ROOT::Math::PxPyPzEVector &pHER, const ROOT::Math::PxPyPzEVector &pLER, const TMatrixDSym &covHER, const TMatrixDSym &covLER)
Initialize the conditional generator using HER & LER 4-vectors and HER & LER covariance matrices desc...

◆ gmunu_tilde()

EvtTensor4C gmunu_tilde ( const EvtVector4R &  p4a,
const EvtVector4R &  p4b,
const EvtVector4R &  p4c 
)

Function Tensor gmunu

Definition at line 409 of file EvtB0toKsKK.cc.

412 {
413 const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
414
415 const EvtVector4R umu = (p4a + p4b + p4c) / sqrt(s);
416
417 const EvtTensor4C gmunu_tilde =
418 EvtTensor4C::g() - EvtGenFunctions::directProd(umu, umu);
419
420 return gmunu_tilde;
421 }
EvtVector4R umu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector umu.
Definition: EvtB0toKsKK.cc:360

◆ hA1()

double hA1 ( double  w) const

HQET D* axial vector form factor h_{A1}(w).


Definition at line 368 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

369 {
370 return ffA1(w);
371 }

◆ hA2()

double hA2 ( double  w) const

HQET D* axial vector form factor h_{A2}(w).


Definition at line 380 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

381 {
382 return (ffR2(w) - ffR3(w)) * ffA1(w) / (2 * r(1));
383 }
double r(int Dhel) const
Ratio of the daughter meson mass to the parent meson.

◆ hA3()

double hA3 ( double  w) const

HQET D* axial vector form factor h_{A3}(w).


Definition at line 386 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

387 {
388 return (ffR2(w) + ffR3(w)) * ffA1(w) / 2;
389 }

◆ HadS1()

double HadS1 ( int  Dhel,
double  w 
) const

The function to calculate the Hadronic Amplitudes of scalar (S+P) type contribution.

Parameters
Dhelhelicity of the daughter D(*) meson {+1,0,1} for D* and 2 for D.
wvelocity transfer variable.
Returns
calculated amplitude value.

Definition at line 243 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

244 {
245 // sanity check
246 assert(chkDhel(Dhel));
247
248 if (Dhel == 2) return m_mB * sqrt(r(2)) * (w + 1) * hS(w);
249 if (Dhel == 0) return -m_mB * sqrt(r(0)) * sqrt(w * w - 1) * hP(w);
250 return 0.;
251
252 }
double hS(double w) const
D scalar form factor h_S(w) in terms of vector form factors.
double hP(double w) const
D* pseudo scalar form factor h_P(w) in terms of axial vector form factors.

◆ HadS2()

double HadS2 ( int  Dhel,
double  w 
) const

The function to calculate the Hadronic Amplitudes of scalar (S-P) type contribution.

Parameters
Dhelhelicity of the daughter D(*) meson {+1,0,1} for D* and 2 for D.
wvelocity transfer variable.
Returns
calculated amplitude value.

Definition at line 255 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

256 {
257 // sanity check
258 assert(chkDhel(Dhel));
259
260 if (Dhel == 2) return HadS1(2, w);
261 else return -HadS1(Dhel, w);
262 }
double HadS1(int Dhel, double w) const
The function to calculate the Hadronic Amplitudes of scalar (S+P) type contribution.

◆ HadT()

double HadT ( int  Dhel,
int  whel1,
int  whel2,
double  w 
) const

The function to calculate the Hadronic Amplitudes of tensor type contribution.

Parameters
Dhelhelicity of the daughter D(*) meson {+1,0,1} for D* and 2 for D.
whel1helicity of the one virtual vector boson {+1,0,1,2}.
whel2helicity of the another virtual vector boson {+1,0,1,2}.
wvelocity transfer variable.
Returns
calculated amplitude value. The tensor contribution is described by two vector contributions.

Definition at line 265 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

266 {
267 // sanity check
268 assert(chkDhel(Dhel) && chkwhel(whel1) && chkwhel(whel2));
269
270 if (whel1 == whel2)return 0.;
271 if (whel1 > whel2)return -HadT(Dhel, whel2, whel1, w);
272
273 const double r0 = r(Dhel);
274 if (Dhel == 2 && whel1 == -1 && whel2 == +1)return m_mB * sqrt(r(2)) * sqrt(w * w - 1) * hT(w);
275 if (Dhel == 2 && whel1 == 0 && whel2 == 2) return -HadT(2, -1, +1, w);
276 if (Dhel == -1 && whel1 == -1 && whel2 == 0)
277 return -m_mB * sqrt(r0 / qh2(-1, w)) * (1 - r0 * (w + sqrt(w * w - 1))) *
278 (hT1(w) + hT2(w) + (w - sqrt(w * w - 1)) * (hT1(w) - hT2(w)));
279 if (Dhel == -1 && whel1 == -1 && whel2 == 2) return -HadT(-1, -1, 0, w);
280 if (Dhel == 0 && whel1 == -1 && whel2 == +1)
281 return m_mB * sqrt(r0) * ((w + 1) * hT1(w) + (w - 1) * hT2(w) + 2 * (w * w - 1) * hT3(w));
282 if (Dhel == 0 && whel1 == 0 && whel2 == 2) return -HadT(0, -1, +1, w);
283 if (Dhel == +1 && whel1 == 0 && whel2 == +1)
284 return -m_mB * sqrt(r0 / qh2(+1, w)) * (1 - r0 * (w - sqrt(w * w - 1))) *
285 (hT1(w) + hT2(w) + (w + sqrt(w * w - 1)) * (hT1(w) - hT2(w)));
286
287 if (Dhel == +1 && whel1 == +1 && whel2 == 2) return -HadT(+1, 0, +1, w);
288
289 // other cases
290 return 0.;
291 }
double qh2(int Dhel, double w) const
Function to calculate the q^2 divided by the square of parent mass (m_B^2).
double hT1(double w) const
D* tensor form factor h_{T1}(w) in terms of axial vector form factors.
bool chkwhel(int whel) const
Function to check if whel is in the valid range.
double HadT(int Dhel, int whel1, int whel2, double w) const
The function to calculate the Hadronic Amplitudes of tensor type contribution.
double hT3(double w) const
D* tensor form factor h_{T3}(w).
double hT(double w) const
D tensor form factor h_T(w) in terms of vector form factors.
double hT2(double w) const
D* tensor form factor h_{T2}(w).

◆ HadV1()

double HadV1 ( int  Dhel,
int  whel,
double  w 
) const

The function to calculate the Hadronic Amplitudes of left handed (V-A) type contribution.

Parameters
Dhelhelicity of the daughter D(*) meson {+1,0,1} for D* and 2 for D.
whelhelicity of the virtual vector boson {+1,0,1,2}.
wvelocity transfer variable.
Returns
calculated amplitude value.

Definition at line 198 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

199 {
200 // sanity check
201 assert(chkDhel(Dhel) && chkwhel(whel));
202
203 const double r0 = r(Dhel);
204 if (Dhel == 2 && whel == 0) {
205 return m_mB * sqrt(r0 * (w * w - 1.) / qh2(2, w)) * ((1 + r0) * hp(w) - (1 - r0) * hm(w));
206 }
207 if (Dhel == 2 && whel == 2) {
208 return m_mB * sqrt(r0 / qh2(2, w)) * ((1 - r0) * (w + 1) * hp(w) - (1 + r0) * (w - 1) * hm(w));
209 }
210 if (Dhel == +1 && whel == +1) {
211 return m_mB * sqrt(r0) * ((w + 1) * hA1(w) - sqrt(w * w - 1) * hV(w));
212 }
213 if (Dhel == -1 && whel == -1) {
214 return m_mB * sqrt(r0) * ((w + 1) * hA1(w) + sqrt(w * w - 1) * hV(w));
215 }
216 if (Dhel == 0 && whel == 0) {
217 return m_mB * sqrt(r0 / qh2(0, w)) * (w + 1) *
218 (-(w - r0) * hA1(w) + (w - 1) * (r0 * hA2(w) + hA3(w)));
219 }
220 if (Dhel == 0 && whel == 2) {
221 return m_mB * sqrt(r0 * (w * w - 1) / qh2(0, w)) * (-(w + 1) * hA1(w) + (1 - r0 * w) * hA2(w) + (w - r0) * hA3(w));
222 }
223
224 // other cases
225 return 0.;
226
227 }
double hV(double w) const
HQET D* axial vector form factor h_V(w).
double hm(double w) const
HQET D vector form factor h_-(w).
double hA1(double w) const
HQET D* axial vector form factor h_{A1}(w).
double hA2(double w) const
HQET D* axial vector form factor h_{A2}(w).
double hp(double w) const
HQET D vector form factor h_+(w).
double hA3(double w) const
HQET D* axial vector form factor h_{A3}(w).

◆ HadV2()

double HadV2 ( int  Dhel,
int  whel,
double  w 
) const

The function to calculate the Hadronic Amplitudes of right handed (V+A) type contribution.

Parameters
Dhelhelicity of the daughter D(*) meson {+1,0,1} for D* and 2 for D.
whelhelicity of the virtual vector boson {+1,0,1,2}.
wvelocity transfer variable.
Returns
calculated amplitude value.

Definition at line 230 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

231 {
232 // sanity check
233 assert(chkDhel(Dhel) && chkwhel(whel));
234
235 if (Dhel == 2) return HadV1(2, whel, w);
236 if (Dhel == +1 && whel == +1)return -HadV1(-1, -1, w);
237 if (Dhel == -1 && whel == -1)return -HadV1(+1, +1, w);
238 if (Dhel == 0) return -HadV1(0, whel, w);
239 return 0;
240 }
double HadV1(int Dhel, int whel, double w) const
The function to calculate the Hadronic Amplitudes of left handed (V-A) type contribution.

◆ helAmp() [1/2]

EvtComplex helAmp ( const EvtComplex &  CV1,
const EvtComplex &  CV2,
const EvtComplex &  CS1,
const EvtComplex &  CS2,
const EvtComplex &  CT,
double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

The function calculates helicity amplitudes with given Wilson coefficients.

Parameters
CV1Wilson coeffcient of the left handed vector type NP contribution.
CV2Wilson coeffcient of the right handed vector type NP contribution.
CS1Wilson coeffcient of the scalar (S+P) type NP contribution.
CS2Wilson coeffcient of the scalar (S-P) type NP contribution.
CTWilson coeffcient of the tensor type NP contribution.
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame. The overall factor GF/sqrt(2) Vcb omitted because it does not change the distribution.

Definition at line 109 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

112 {
113 // sanity check
114 assert(chktauhel(tauhel) && chkDhel(Dhel));
115
116 return //(GF/sqrt(2))*Vcb* // <-- constants which does not affect the distribution omitted
117 (1.*helampSM(mtau, tauhel, Dhel, w, costau)
118 + CV1 * helampV1(mtau, tauhel, Dhel, w, costau)
119 + CV2 * helampV2(mtau, tauhel, Dhel, w, costau)
120 + CS1 * helampS1(mtau, tauhel, Dhel, w, costau)
121 + CS2 * helampS2(mtau, tauhel, Dhel, w, costau)
122 + CT * helampT(mtau, tauhel, Dhel, w, costau));
123 }
double helampT(double mtau, int tauhel, int Dhel, double w, double costau) const
Helicity Amplitudes of tensor type contribution.
double helampV1(double mtau, int tauhel, int Dhel, double w, double costau) const
Helicity Amplitudes of left handed (V-A) contribution.
bool chktauhel(int tauhel) const
Function to check if tauhel is in the valid range.
double helampS1(double mtau, int tauhel, int Dhel, double w, double costau) const
Helicity Amplitudes of scalar (S+P) type contribution.
double helampSM(double mtau, int tauhel, int Dhel, double w, double costau) const
Helicity Amplitudes of SM (left handed) contribution.
double helampV2(double mtau, int tauhel, int Dhel, double w, double costau) const
Helicity Amplitudes of right handed (V+A) contribution.
double helampS2(double mtau, int tauhel, int Dhel, double w, double costau) const
Helicity Amplitudes of scalar (S-P) type contribution.

◆ helAmp() [2/2]

EvtComplex helAmp ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

The function calculates the helicity amplitude.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame. The overall factor GF/sqrt(2) Vcb omitted because it does not change the distribution.

Definition at line 102 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

103 {
104 return helAmp(m_CV1, m_CV2, m_CS1, m_CS2, m_CT,
105 mtau, tauhel, Dhel, w, costau);
106 }
EvtComplex helAmp(double mtau, int tauhel, int Dhel, double w, double costau) const
The function calculates the helicity amplitude.

◆ helampS1()

double helampS1 ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

Helicity Amplitudes of scalar (S+P) type contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value. Overall factor GF/sqrt(2) Vcb omitted. Wilson coefficients CXX ommited.

Definition at line 326 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

327 {
328 return -Lep(mtau, tauhel, q2(Dhel, w), costau) * HadS1(Dhel, w);
329 }
double Lep(const double mtau, int tauhel, int whel, double q2, double costau) const
The function to calculate the Leptonic Amplitudes for B->D*taunu decay of the vector type contributio...

◆ helampS2()

double helampS2 ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

Helicity Amplitudes of scalar (S-P) type contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value. Overall factor GF/sqrt(2) Vcb omitted. Wilson coefficients CXX ommited.

Definition at line 332 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

333 {
334 return -Lep(mtau, tauhel, q2(Dhel, w), costau) * HadS2(Dhel, w);
335 }
double HadS2(int Dhel, double w) const
The function to calculate the Hadronic Amplitudes of scalar (S-P) type contribution.

◆ helampSM()

double helampSM ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

Helicity Amplitudes of SM (left handed) contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame. Overall factor GF/sqrt(2) Vcb omitted. Wilson coefficients CXX ommited.

Definition at line 298 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

299 {
300 double amp(0.);
301 for (int whel = -1; whel <= 2; whel++) {
302 amp += eta(whel) * Lep(mtau, tauhel, whel, q2(Dhel, w), costau)
303 * HadV1(Dhel, whel, w);
304 }
305 return amp;
306 }

◆ helampT()

double helampT ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

Helicity Amplitudes of tensor type contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value. Overall factor GF/sqrt(2) Vcb omitted. Wilson coefficients CXX ommited.

Definition at line 338 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

339 {
340 double amp(0.);
341 for (int whel1 = -1; whel1 <= 2; whel1++) {
342 for (int whel2 = -1; whel2 <= 2; whel2++) {
343 amp += -1 * eta(whel1) * eta(whel2) * Lep(mtau, tauhel, whel1, whel2, q2(Dhel, w), costau)
344 * HadT(Dhel, whel1, whel2, w);
345 }
346 }
347 return amp;
348 }

◆ helampV1()

double helampV1 ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

Helicity Amplitudes of left handed (V-A) contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value. Overall factor GF/sqrt(2) Vcb omitted. Wilson coefficients CXX ommited.

Definition at line 309 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

310 {
311 return helampSM(mtau, tauhel, Dhel, w, costau);
312 }

◆ helampV2()

double helampV2 ( double  mtau,
int  tauhel,
int  Dhel,
double  w,
double  costau 
) const

Helicity Amplitudes of right handed (V+A) contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value. Overall factor GF/sqrt(2) Vcb omitted. Wilson coefficients CXX ommited.

Definition at line 315 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

316 {
317 double amp(0.);
318 for (int whel = -1; whel <= 2; whel++) {
319 amp += eta(whel) * Lep(mtau, tauhel, whel, q2(Dhel, w), costau)
320 * HadV2(Dhel, whel, w);
321 }
322 return amp;
323 }
double HadV2(int Dhel, int whel, double w) const
The function to calculate the Hadronic Amplitudes of right handed (V+A) type contribution.

◆ hm()

double hm ( double  w) const

HQET D vector form factor h_-(w).


Definition at line 362 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

363 {
364 return - (1 - r(2) * r(2)) * (w + 1) * (ffV1(w) - ffS1(w)) / (2 * qh2(2, w));
365 }

◆ hp()

double hp ( double  w) const

HQET D vector form factor h_+(w).


Definition at line 355 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

356 {
357 const double r0 = r(2);
358 return -((1 + r0) * (1 + r0) * (w - 1) * ffV1(w)
359 - (1 - r0) * (1 - r0) * (w + 1) * ffS1(w)) / (2 * qh2(2, w));
360 }

◆ hP()

double hP ( double  w) const

D* pseudo scalar form factor h_P(w) in terms of axial vector form factors.


Definition at line 405 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

406 {
407 // M. Tanaka, private communication, without approximating quark masses to the meson mass
408 const double r0 = r(1);
409 return m_mB / m_mBottom * ((w + 1) * hA1(w) - (1 - r0 * w) * hA2(w) - (w - r0) * hA3(w)) / (1 + rq());
410
411 // when quark masses are approximated to the meson masses
412 //return ( (w+1) * hA1(w) - (1-r0*w) * hA2(w) - (w-r0) * hA3(w) ) / (1+r0); // no Lambda/mQ
413 }
double rq() const
Ratio of the charm quark mass to the charm quark mass.

◆ hS()

double hS ( double  w) const

D scalar form factor h_S(w) in terms of vector form factors.


Definition at line 393 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

394 {
395 // M. Tanaka, private communication, without approximating quark masses to the meson mass
396 const double r0 = r(2);
397 return m_mB / m_mBottom * ((1 - r0) / (1 - rq()) * hp(w)
398 - (1 + r0) / (1 - rq()) * (w - 1) / (w + 1) * hm(w));
399
400 // when quark masses are approximated to the meson masses
401 //return ffS1(w);
402 }

◆ hT()

double hT ( double  w) const

D tensor form factor h_T(w) in terms of vector form factors.


Definition at line 417 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

418 {
419 const double r0 = r(2);
420 return hp(w) - (1 + r0) * hm(w) / (1 - r0);
421 }

◆ hT1()

double hT1 ( double  w) const

D* tensor form factor h_{T1}(w) in terms of axial vector form factors.


Definition at line 424 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

425 {
426 const double r0 = r(1);
427 return -((1 + r0) * (1 + r0) * (w - 1) * hV(w)
428 - (1 - r0) * (1 - r0) * (w + 1) * hA1(w))
429 / (2 * qh2(1, w));
430 }

◆ hT2()

double hT2 ( double  w) const

D* tensor form factor h_{T2}(w).


Definition at line 433 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

434 {
435 const double r0 = r(1);
436 return -(1 - r0 * r0) * (w + 1) * (hV(w) - hA1(w))
437 / (2 * qh2(1, w));
438 }

◆ hT3()

double hT3 ( double  w) const

D* tensor form factor h_{T3}(w).


Definition at line 441 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

442 {
443 const double r0 = r(1);
444 return ((1 + r0) * (1 + r0) * hV(w) - 2 * r0 * (w + 1) * hA1(w)
445 + qh2(1, w) * (r0 * hA2(w) - hA3(w)))
446 / (2 * qh2(1, w) * (1 + r0));
447 }

◆ hV()

double hV ( double  w) const

HQET D* axial vector form factor h_V(w).


Definition at line 374 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

375 {
376 return ffR1(w) * ffA1(w);
377 }

◆ init() [1/13]

void init ( )

Initializes the generator.

Definition at line 37 of file CRY.cc.

38 {
39 std::stringstream setupString;
40 setupString << " returnGammas " << m_returnGammas << std::endl;
41 setupString << " returnKaons " << m_returnKaons << std::endl;
42 setupString << " returnPions " << m_returnPions << std::endl;
43 setupString << " returnProtons " << m_returnProtons << std::endl;
44 setupString << " returnNeutrons " << m_returnNeutrons << std::endl;
45 setupString << " returnElectrons " << m_returnElectrons << std::endl;
46 setupString << " returnMuons " << m_returnMuons << std::endl;
47 setupString << " date " << m_date << std::endl;
48 setupString << " latitude " << 36.0 << std::endl;
49 setupString << " altitude " << 0 << std::endl;
50 setupString << " subboxLength " << m_subboxLength << std::endl;
51
52 IOIntercept::OutputToLogMessages initLogCapture("CRY", LogConfig::c_Debug, LogConfig::c_Warning, 100, 100);
53 initLogCapture.start();
54
55 // CRY setup
56 m_crySetup.reset(new CRYSetup(setupString.str(), m_cosmicDataDir));
57 // Set the random generator to the setup
58 m_crySetup->setRandomFunction([]()->double { return gRandom->Rndm(); });
59 // Set up the generator
60 m_cryGenerator.reset(new CRYGenerator(m_crySetup.get()));
61
62 initLogCapture.finish();
63
64 // set up the acceptance box
65 const auto coords = m_acceptSize.size();
66 if (coords == 1) {
67 m_acceptance.reset(new vecgeom::UnplacedOrb(m_acceptSize[0]));
68 } else if (coords == 2) {
69 m_acceptance.reset(new vecgeom::SUnplacedTube<vecgeom::TubeTypes::NonHollowTube>(0, m_acceptSize[0], m_acceptSize[1], 0, 2 * M_PI));
70 } else if (coords == 3) {
71 m_acceptance.reset(new vecgeom::UnplacedBox(m_acceptSize[0], m_acceptSize[1], m_acceptSize[2]));
72 } else {
73 B2FATAL("Acceptance volume needs to have one, two or three values for sphere, cylinder and box respectively");
74 }
75 //get information from geometry and create the world box
76 G4VPhysicalVolume* volume = geometry::GeometryManager::getInstance().getTopVolume();
77 if (!volume) {
78 B2FATAL("No geometry found -> Add the Geometry module to the path before the CRY module.");
79 }
80 G4Box* topbox = (G4Box*) volume->GetLogicalVolume()->GetSolid();
81 if (!topbox) {
82 B2FATAL("No G4Box found -> Check the logical volume of the geometry.");
83 }
84
85 // Wrap the world coordinates (G4 coordinates are mm, Belle 2 units are currently cm)
86 double halfLength_x = topbox->GetXHalfLength() * Belle2::Unit::mm;
87 double halfLength_y = topbox->GetYHalfLength() * Belle2::Unit::mm;
88 double halfLength_z = topbox->GetZHalfLength() * Belle2::Unit::mm;
89
90 // m_world.reset(new vecgeom::UnplacedBox(m_subboxLength / 2. * Unit::m, m_subboxLength / 2. * Unit::m,
91 // m_subboxLength / 2. * Unit::m));
92 m_world.reset(new vecgeom::UnplacedBox(halfLength_x / Unit::cm, halfLength_y / Unit::cm,
93 halfLength_z / Unit::cm));
94 B2INFO("World size [" << halfLength_x / Unit::cm << ", " << halfLength_x / Unit::cm << ", " << halfLength_x / Unit::cm << "]");
95
96 const double maxSize = *std::max_element(m_acceptSize.begin(), m_acceptSize.end());
97 const double minWorld = std::min({halfLength_x, halfLength_y, halfLength_z, });
98 if (maxSize > (minWorld / Unit::cm)) {
99 B2FATAL("Acceptance bigger than world volume " << LogVar("acceptance", maxSize)
100 << LogVar("World size", minWorld / Unit::cm));
101 }
102
103 }
bool m_returnElectrons
Whether or not CRY should return electrons.
Definition: CRY.h:154
std::string m_cosmicDataDir
directory that holds cosmic data files.
Definition: CRY.h:141
bool m_returnKaons
Whether or not CRY should return kaons.
Definition: CRY.h:150
bool m_returnMuons
Whether or not CRY should return muons.
Definition: CRY.h:155
int m_subboxLength
length of the square n-n plane in Cry in meters
Definition: CRY.h:142
bool m_returnGammas
Whether or not CRY should return gammas.
Definition: CRY.h:149
std::vector< double > m_acceptSize
Shape parameters for the acceptance shape.
Definition: CRY.h:143
std::string m_date
date used for generation (month-day-year).
Definition: CRY.h:146
std::unique_ptr< CRYSetup > m_crySetup
The CRY generator setup.
Definition: CRY.h:157
bool m_returnProtons
Whether or not CRY should return protons.
Definition: CRY.h:152
bool m_returnNeutrons
Whether or not CRY should return neutrons.
Definition: CRY.h:153
bool m_returnPions
Whether or not CRY should return pions.
Definition: CRY.h:151
@ c_Debug
Debug: for code development.
Definition: LogConfig.h:26
@ c_Warning
Warning: for potential problems that the user should pay attention to.
Definition: LogConfig.h:29
static const double mm
[millimeters]
Definition: Unit.h:70
static const double cm
Standard units with the value = 1.
Definition: Unit.h:47
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
static GeometryManager & getInstance()
Return a reference to the instance.

◆ init() [2/13]

void init ( )

Initialize standard stream objects

Definition at line 41 of file EvtB0toKsKK.cc.

42 {
43 // Check number of arguments
44 checkNArg(32);
45
46 // Check number of daughters
47 checkNDaug(3);
48
49 // Check decay chain
50 if ((getParentId() != EvtPDL::getId("B0") &&
51 getParentId() != EvtPDL::getId("anti-B0")) ||
52 (getDaug(0) != EvtPDL::getId("K_S0")) ||
53 (getDaug(1) != EvtPDL::getId("K+") &&
54 getDaug(1) != EvtPDL::getId("K-")) ||
55 (getDaug(2) != EvtPDL::getId("K+") &&
56 getDaug(2) != EvtPDL::getId("K-"))) {
57 std::cout << "ERROR: Invalid decay" << std::endl;
58 std::cout << "USAGE: K_S0 K+ K-" << std::endl;
59 exit(1);
60 }
61
62 a_f0ks_ =
63 (1.0 + getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
64 getArg(1) * sin(getArg(2) * M_PI / 180.0));
65 a_phiks_ =
66 (1.0 + getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
67 getArg(5) * sin(getArg(6) * M_PI / 180.0));
68 a_fxks_ =
69 (1.0 + getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
70 getArg(9) * sin(getArg(10) * M_PI / 180.0));
72 (1.0 + getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
73 getArg(13) * sin(getArg(14) * M_PI / 180.0));
74 a_kpkmnr_ =
75 (1.0 + getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
76 getArg(17) * sin(getArg(18) * M_PI / 180.0));
77 a_kskpnr_ =
78 (1.0 + getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
79 getArg(22) * sin(getArg(23) * M_PI / 180.0));
80 a_kskmnr_ =
81 (1.0 + getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
82 getArg(27) * sin(getArg(28) * M_PI / 180.0));
83
85 (1.0 - getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
86 getArg(1) * sin(getArg(2) * M_PI / 180.0));
88 (1.0 - getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
89 getArg(5) * sin(getArg(6) * M_PI / 180.0));
91 (1.0 - getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
92 getArg(9) * sin(getArg(10) * M_PI / 180.0));
94 (1.0 - getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
95 getArg(13) * sin(getArg(14) * M_PI / 180.0));
97 (1.0 - getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
98 getArg(17) * sin(getArg(18) * M_PI / 180.0));
100 (1.0 - getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
101 getArg(22) * sin(getArg(23) * M_PI / 180.0));
103 (1.0 - getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
104 getArg(27) * sin(getArg(28) * M_PI / 180.0));
105
106 alpha_kpkmnr = getArg(21);
107 alpha_kskpnr = getArg(26);
108 alpha_kskmnr = getArg(31);
109
110 std::cout << setiosflags(std::ios::left) << std::endl
111 << "B0 Channel " << std::setw(20) << "Relative amplitude" << std::setw(20) << "Relative phase" << std::endl
112 << "f0ks " << std::setw(20) << real(a_f0ks_) << std::setw(20) << imag(a_f0ks_) << std::endl
113 << "phiks " << std::setw(20) << real(a_phiks_) << std::setw(20) << imag(a_phiks_) << std::endl
114 << "fxks " << std::setw(20) << real(a_fxks_) << std::setw(20) << imag(a_fxks_) << std::endl
115 << "chic0ks " << std::setw(20) << real(a_chic0ks_) << std::setw(20) << imag(a_chic0ks_) << std::endl
116 << "kpkmnr " << std::setw(20) << real(a_kpkmnr_) << std::setw(20) << imag(a_kpkmnr_) << std::endl
117 << "kskpnr " << std::setw(20) << real(a_kskpnr_) << std::setw(20) << imag(a_kskpnr_) << std::endl
118 << "kskmnr " << std::setw(20) << real(a_kskmnr_) << std::setw(20) << imag(a_kskmnr_) << std::endl
119 << std::endl;
120
121 std::cout << setiosflags(std::ios::left) << std::endl
122 << "B0B Channel " << std::setw(20) << "Relative amplitude" << std::setw(20) << "Relative phase" << std::endl
123 << "f0ks " << std::setw(20) << real(abar_f0ks_) << std::setw(20) << imag(abar_f0ks_) << std::endl
124 << "phiks " << std::setw(20) << real(abar_phiks_) << std::setw(20) << imag(abar_phiks_) << std::endl
125 << "fxks " << std::setw(20) << real(abar_fxks_) << std::setw(20) << imag(abar_fxks_) << std::endl
126 << "chic0ks " << std::setw(20) << real(abar_chic0ks_) << std::setw(20) << imag(abar_chic0ks_) << std::endl
127 << "kpkmnr " << std::setw(20) << real(abar_kpkmnr_) << std::setw(20) << imag(abar_kpkmnr_) << std::endl
128 << "kskpnr " << std::setw(20) << real(abar_kskpnr_) << std::setw(20) << imag(abar_kskpnr_) << std::endl
129 << "kskmnr " << std::setw(20) << real(abar_kskmnr_) << std::setw(20) << imag(abar_kskmnr_) << std::endl
130 << std::endl;
131
132 // Open debugging file
133 debugfile_.open("debug.dat", std::ios::out);
134 }

◆ init() [3/13]

void init ( )

Initializes module.

Definition at line 72 of file EvtBCL.cc.

73 {
74
75 checkNDaug(3);
76
77 //We expect the parent to be a scalar
78 //and the daughters to be X lepton neutrino
79
80 checkSpinParent(EvtSpinType::SCALAR);
81 checkSpinDaughter(1, EvtSpinType::DIRAC);
82 checkSpinDaughter(2, EvtSpinType::NEUTRINO);
83
84 EvtSpinType::spintype mesontype = EvtPDL::getSpinType(getDaug(0));
85
86 bclmodel = new EvtBCLFF(getNArg(), getArgs());
87
88 if (mesontype == EvtSpinType::SCALAR) {
89 calcamp = new EvtSemiLeptonicScalarAmp;
90 }
91 if (mesontype == EvtSpinType::VECTOR) {
92 calcamp = new EvtSemiLeptonicVectorAmp;
93 }
94 // Tensor Meson implementation is possible here.
95
96 }

◆ init() [4/13]

void init ( )

The function initializes the decay.

The function reads paramters from the decay file and initializes the parameters for the decay amplitude calculator.

Definition at line 78 of file EvtBSemiTauonic.cc.

79 {
80 checkNDaug(3);
81
82 //We expect the parent to be a scalar
83 //and the daughters to be X lepton neutrino
84 checkSpinParent(EvtSpinType::SCALAR);
85
86 checkSpinDaughter(1, EvtSpinType::DIRAC);
87 checkSpinDaughter(2, EvtSpinType::NEUTRINO);
88
89 EvtSpinType::spintype d1type = EvtPDL::getSpinType(getDaug(0));
90
91 if (d1type == EvtSpinType::VECTOR) {
92 B2INFO("EvtBSemiTauonic::init()>> Initializing for decays to a vector type meson ");
93 checkNArg(16);
94 const double rhoa12 = getArg(0);
95 const double R11 = getArg(1);
96 const double R21 = getArg(2);
97 const double aR3 = getArg(3);
98
99 const double rho12 = 0; // not used for D*taunu decays
100 const double aS1 = 0; // not used for D*taunu decays
101
102 const double m_b = getArg(4);
103 const double m_c = getArg(5);
104
105 EvtComplex Coeffs[5];
106 Coeffs[0] = EvtComplex(getArg(6) * cos(getArg(7)), getArg(6) * sin(getArg(7)));
107 Coeffs[1] = EvtComplex(getArg(8) * cos(getArg(9)), getArg(8) * sin(getArg(9)));
108 Coeffs[2] = EvtComplex(getArg(10) * cos(getArg(11)), getArg(10) * sin(getArg(11)));
109 Coeffs[3] = EvtComplex(getArg(12) * cos(getArg(13)), getArg(12) * sin(getArg(13)));
110 Coeffs[4] = EvtComplex(getArg(14) * cos(getArg(15)), getArg(14) * sin(getArg(15)));
111
112 m_CalcHelAmp = new EvtBSemiTauonicHelicityAmplitudeCalculator(rho12, rhoa12, R11, R21, aS1, aR3, m_b, m_c,
113 Coeffs[0], Coeffs[1], Coeffs[2], Coeffs[3], Coeffs[4],
114 EvtPDL::getMeanMass(getParentId()),
115 -1, /*dummy for Dmass*/
116 EvtPDL::getMeanMass(getDaug(0)));
117 m_CalcAmp = new EvtBSemiTauonicVectorMesonAmplitude();
118 } else if (d1type == EvtSpinType::SCALAR) {
119 B2INFO("EvtBSemiTauonic::init()>> Initializing for decays to a scalar type meson ");
120 checkNArg(14);
121 const double rho12 = getArg(0);
122 const double aS1 = getArg(1);
123
124 const double rhoa12 = 0; // not used for Dtaunu decays
125 const double R11 = 0; // not used for Dtaunu decays
126 const double R21 = 0; // not used for Dtaunu decays
127 const double aR3 = 0; // not used for Dtaunu decays
128
129 const double m_b = getArg(2);
130 const double m_c = getArg(3);
131
132 EvtComplex Coeffs[5];
133 Coeffs[0] = EvtComplex(getArg(4) * cos(getArg(5)), getArg(4) * sin(getArg(5)));
134 Coeffs[1] = EvtComplex(getArg(6) * cos(getArg(7)), getArg(6) * sin(getArg(7)));
135 Coeffs[2] = EvtComplex(getArg(8) * cos(getArg(9)), getArg(8) * sin(getArg(9)));
136 Coeffs[3] = EvtComplex(getArg(10) * cos(getArg(11)), getArg(10) * sin(getArg(11)));
137 Coeffs[4] = EvtComplex(getArg(12) * cos(getArg(13)), getArg(12) * sin(getArg(13)));
138
139 m_CalcHelAmp = new EvtBSemiTauonicHelicityAmplitudeCalculator(rho12, rhoa12, R11, R21, aS1, aR3, m_b, m_c,
140 Coeffs[0], Coeffs[1], Coeffs[2], Coeffs[3], Coeffs[4],
141 EvtPDL::getMeanMass(getParentId()),
142 EvtPDL::getMeanMass(getDaug(0)),
143 -1 /*dummy for D*mass*/);
144 m_CalcAmp = new EvtBSemiTauonicScalarMesonAmplitude();
145 } else {
146 B2ERROR("BSemiTauonic model handles only scalar and vector meson daughters.");
147// EvtGenReport(EVTGEN_ERROR, "EvtGen") << "BSemiTauonic model handles only scalar and vector meson daughters." << std::endl;
148 ::abort();
149 }
150 }

◆ init() [5/13]

void init ( )

The function initializes the decay.

The function reads paramters from the decay file and initializes the parameters for the decay amplitude calculator.

Definition at line 78 of file EvtBSemiTauonic2HDMType2.cc.

79 {
80 checkNDaug(3);
81
82 //We expect the parent to be a scalar
83 //and the daughters to be X lepton neutrino
84 checkSpinParent(EvtSpinType::SCALAR);
85
86 checkSpinDaughter(1, EvtSpinType::DIRAC);
87 checkSpinDaughter(2, EvtSpinType::NEUTRINO);
88
89 const double m_tau = EvtPDL::getMeanMass(getDaug(1));
90
91 EvtSpinType::spintype d1type = EvtPDL::getSpinType(getDaug(0));
92
93 if (d1type == EvtSpinType::VECTOR) {
94 B2INFO("EvtBSemiTauonic2HDMType2::init()>> Initializing for decays to a vector type meson ");
95 checkNArg(7, 8);
96
97 const double rhoa12 = getArg(0);
98 const double R11 = getArg(1);
99 const double R21 = getArg(2);
100 const double aR3 = getArg(3);
101
102 const double rho12 = 0; // not used for D*taunu decays
103 const double aS1 = 0; // not used for D*taunu decays
104
105 EvtComplex Coeffs[5];
106
107 const double m_b = getArg(4);
108 const double m_c = getArg(5);
109 const double tanBetaOverMH = getArg(6);
110 B2INFO("tan(beta)/m_H+ = " << tanBetaOverMH);
111
112 Coeffs[0] = 0; // CV1
113 Coeffs[1] = 0; // CV2
114 Coeffs[2] = -m_b * m_tau * tanBetaOverMH * tanBetaOverMH; // CS1
115 Coeffs[3] = 0; // CS2, neglected by default
116 Coeffs[4] = 0; // CT
117
118 if (getNArg() == 8) {
119 // if m_H is explicitly given
120 const double m_H = getArg(7);
121 Coeffs[3] = -m_c * m_tau / m_H / m_H; // CS2
122 }
123
124 m_CalcHelAmp = new EvtBSemiTauonicHelicityAmplitudeCalculator(rho12, rhoa12, R11, R21, aS1, aR3, m_b, m_c,
125 Coeffs[0], Coeffs[1], Coeffs[2], Coeffs[3], Coeffs[4],
126 EvtPDL::getMeanMass(getParentId()),
127 -1, /*dummy for Dmass*/
128 EvtPDL::getMeanMass(getDaug(0)));
129 m_CalcAmp = new EvtBSemiTauonicVectorMesonAmplitude();
130 } else if (d1type == EvtSpinType::SCALAR) {
131 B2INFO("EvtBSemiTauonic2HDMType2::init()>> Initializing for decays to a scalar type meson ");
132 checkNArg(5, 6);
133
134 const double rho12 = getArg(0);
135 const double aS1 = getArg(1);
136
137 const double rhoa12 = 0; // not used for Dtaunu decays
138 const double R11 = 0; // not used for Dtaunu decays
139 const double R21 = 0; // not used for Dtaunu decays
140 const double aR3 = 0; // not used for Dtaunu decays
141
142 EvtComplex Coeffs[5];
143
144 const double m_b = getArg(2);
145 const double m_c = getArg(3);
146 const double tanBetaOverMH = getArg(4);
147 B2INFO("tan(beta)/m_H+ = " << tanBetaOverMH);
148
149 Coeffs[0] = 0; // CV1
150 Coeffs[1] = 0; // CV2
151 Coeffs[2] = -m_b * m_tau * tanBetaOverMH * tanBetaOverMH; // CS1
152 Coeffs[3] = 0; // CS2, neglected by default
153 Coeffs[4] = 0; // CT
154
155 if (getNArg() == 6) {
156 // if m_H is explicitly given
157 const double m_H = getArg(5);
158 Coeffs[3] = -m_c * m_tau / m_H / m_H; // CS2
159 }
160
161 m_CalcHelAmp = new EvtBSemiTauonicHelicityAmplitudeCalculator(rho12, rhoa12, R11, R21, aS1, aR3, m_b, m_c,
162 Coeffs[0], Coeffs[1], Coeffs[2], Coeffs[3], Coeffs[4],
163 EvtPDL::getMeanMass(getParentId()),
164 EvtPDL::getMeanMass(getDaug(0)),
165 -1 /*dummy for D*mass*/);
166 m_CalcAmp = new EvtBSemiTauonicScalarMesonAmplitude();
167 } else {
168 B2ERROR("BSemiTauonic2HDMType2 model handles only scalar and vector meson daughters.");
169// EvtGenReport(EVTGEN_ERROR, "EvtGen") << "BSemiTauonic2HDMType2 model handles only scalar and vector meson daughters."<<std::endl;
170 ::abort();
171 }
172 }

◆ init() [6/13]

void init ( )

The function for an initialization.

Definition at line 167 of file EvtBtoXsnunu_FERMI.cc.

168 {
169
170 // check that there are no arguments
171
172 checkNArg(0, 3, 4, 6);
173
174 checkNDaug(3);
175
176 // Check that the two leptons are the same type
177
178 EvtId lepton1type = getDaug(1);
179 EvtId lepton2type = getDaug(2);
180
181 int etyp = 0;
182 int mutyp = 0;
183 int tautyp = 0;
184 if (lepton1type == EvtPDL::getId("anti-nu_e") ||
185 lepton1type == EvtPDL::getId("nu_e")) {
186 etyp++;
187 }
188 if (lepton2type == EvtPDL::getId("anti-nu_e") ||
189 lepton2type == EvtPDL::getId("nu_e")) {
190 etyp++;
191 }
192 if (lepton1type == EvtPDL::getId("anti-nu_mu") ||
193 lepton1type == EvtPDL::getId("nu_mu")) {
194 mutyp++;
195 }
196 if (lepton2type == EvtPDL::getId("anti-nu_mu") ||
197 lepton2type == EvtPDL::getId("nu_mu")) {
198 mutyp++;
199 }
200 if (lepton1type == EvtPDL::getId("anti-nu_tau") ||
201 lepton1type == EvtPDL::getId("nu_tau")) {
202 tautyp++;
203 }
204 if (lepton2type == EvtPDL::getId("anti-nu_tau") ||
205 lepton2type == EvtPDL::getId("nu_tau")) {
206 tautyp++;
207 }
208
209 if (etyp != 2 && mutyp != 2 && tautyp != 2) {
210
211 std::cout << "Expect two neutrinos of the same type in EvtBtoXsnunu.cc\n";
212 ::abort();
213 }
214
215 // Check that the second and third entries are leptons with positive
216 // and negative charge, respectively
217
218 int lpos = 0;
219 int lneg = 0;
220 if (lepton1type == EvtPDL::getId("anti-nu_e") ||
221 lepton1type == EvtPDL::getId("anti-nu_mu") ||
222 lepton1type == EvtPDL::getId("anti-nu_tau")) {
223 lpos++;
224 } else if (lepton1type == EvtPDL::getId("nu_e") ||
225 lepton1type == EvtPDL::getId("nu_mu") ||
226 lepton1type == EvtPDL::getId("nu_tau")) {
227 lneg++;
228 }
229 if (lepton2type == EvtPDL::getId("nu_e") ||
230 lepton2type == EvtPDL::getId("nu_mu") ||
231 lepton2type == EvtPDL::getId("nu_tau")) {
232 lneg++;
233 } else if (lepton2type == EvtPDL::getId("anti-nu_e") ||
234 lepton2type == EvtPDL::getId("anti-nu_mu") ||
235 lepton2type == EvtPDL::getId("anti-nu_tau")) {
236 lpos++;
237 }
238
239 if (lpos != 1 || lneg != 1) {
240
241 std::cout << "There should be one neutrino and one anti-neutrino in EvtBtoXsnunu.cc\n";
242 ::abort();
243 }
244
245 if (getNArg() >= 3) {
246 // s-quark mass for fermi motion
247 m_ms = getArg(0);
248 // spectator quark mass for fermi motion
249 m_mq = getArg(1);
250 // Fermi motion parameter for fermi motion
251 m_pf = getArg(2);
252 }
253 if (getNArg() >= 4) {
254 m_mxmin = getArg(3);
255 }
256 if (getNArg() == 6) {
257 m_mb_prob = getArg(4);
258 m_ms_prob = getArg(5);
259 }
260
261 // get a maximum probability
262 double mb = m_mb_prob;
263 double ms = m_ms_prob;
264 double mstilda = ms / mb;
265 double mstilda2 = mstilda * mstilda;
266
267 int nsteps = 100;
268 double sbmin = 0;
269 double sbmax = (1 - mstilda) * (1 - mstilda);
270 double probMax = -10000.0;
271 double sbProbMax = -10.0;
272
273 for (int i = 0; i < nsteps; i++) {
274 double sb = sbmin + (i + 0.0005) * (sbmax - sbmin) / (double)nsteps;
275 double lambda = 1 + mstilda2 * mstilda2 + sb * sb - 2 * (mstilda2 + sb + mstilda2 * sb);
276 double prob = sqrt(lambda) * (3 * sb * (1 + mstilda2 - sb) + lambda);
277 if (prob > probMax) {
278 sbProbMax = sb;
279 probMax = prob;
280 }
281 }
282
283 if (verbose()) {
284 std::cout << "dGdsbProbMax = " << probMax << " for sb = " << sbProbMax << std::endl;
285 }
286
287 m_dGdsbProbMax = probMax;
288
289 }

◆ init() [7/13]

void init ( )

Checks that the number of input parameters are correct:

  • 3 scalar particles (check on spin but not on charge)
  • 6 real parameters (no check on the actual values)

Definition at line 41 of file EvtEtaFullDalitz.cc.

42 {
43
44 // check that there is are arguments
45 checkNArg(6);
46 checkNDaug(3);
47
48
49 checkSpinParent(EvtSpinType::SCALAR);
50
51 checkSpinDaughter(0, EvtSpinType::SCALAR);
52 checkSpinDaughter(1, EvtSpinType::SCALAR);
53 checkSpinDaughter(2, EvtSpinType::SCALAR);
54 }

◆ init() [8/13]

void init ( )

Checks that the number of input parameters are correct:

  • 3 scalar particles (check on spin but not on charge)
  • 1 real parameters (no check on the actual values)

Definition at line 56 of file EvtEtaPi0Dalitz.cc.

57 {
58
59 // check that there is 1 argument
60 checkNArg(1);
61 checkNDaug(3);
62
63
64 checkSpinParent(EvtSpinType::SCALAR);
65
66 checkSpinDaughter(0, EvtSpinType::SCALAR);
67 checkSpinDaughter(1, EvtSpinType::SCALAR);
68 checkSpinDaughter(2, EvtSpinType::SCALAR);
69 }

◆ init() [9/13]

void init ( )

Checks that the number of input parameters are correct:

  • 3 scalar particles (check on spin but not on charge)
  • 4 real parameters (no check on the actual values)

Definition at line 58 of file EvtEtaPrimeDalitz.cc.

59 {
60
61 // check that there is are arguments
62 checkNArg(4);
63 checkNDaug(3);
64
65
66 checkSpinParent(EvtSpinType::SCALAR);
67
68 checkSpinDaughter(0, EvtSpinType::SCALAR);
69 checkSpinDaughter(1, EvtSpinType::SCALAR);
70 checkSpinDaughter(2, EvtSpinType::SCALAR);
71 }

◆ init() [10/13]

void init ( )

The function for an initialization.

Definition at line 114 of file EvtFlatDaughter.cc.

115 {
116
117 // check that there are two arguments - Mmin Mmax
118 checkNArg(0);
119
120 // check the number of daughters. It should be 3
121 checkNDaug(3);
122
123 }

◆ init() [11/13]

void init ( )

The function for an initialization.

Definition at line 124 of file EvtKnunu.cc.

125 {
126 // check that there are 0 arguments
127 checkNArg(0);
128 checkNDaug(3);
129
130 //We expect the parent to be a scalar
131 //and the daughters to be K neutrino netrino
132
133 checkSpinParent(EvtSpinType::SCALAR);
134
135 checkSpinDaughter(0, EvtSpinType::SCALAR);
136 checkSpinDaughter(1, EvtSpinType::NEUTRINO);
137 checkSpinDaughter(2, EvtSpinType::NEUTRINO);
138 }

◆ init() [12/13]

void init ( )

The function for an initialization.

Definition at line 152 of file EvtKstarnunu_REV.cc.

153 {
154
155 // check that there are 0 arguments
156 checkNArg(0);
157 checkNDaug(3);
158
159 //We expect the parent to be a scalar
160 //and the daughters to be Kstar neutrino netrino
161
162 checkSpinParent(EvtSpinType::SCALAR);
163
164 checkSpinDaughter(0, EvtSpinType::VECTOR);
165 checkSpinDaughter(1, EvtSpinType::NEUTRINO);
166 checkSpinDaughter(2, EvtSpinType::NEUTRINO);
167 }

◆ init() [13/13]

void init ( )

Initialize standard stream objects

Definition at line 40 of file EvtPhspCP.cc.

41 {
42 // Check number of arguments
43 checkNArg(7);
44 }

◆ initConditionalGenerator()

ConditionalGaussGenerator initConditionalGenerator ( const ROOT::Math::PxPyPzEVector &  pHER,
const ROOT::Math::PxPyPzEVector &  pLER,
const TMatrixDSym &  covHER,
const TMatrixDSym &  covLER 
)

Initialize the conditional generator using HER & LER 4-vectors and HER & LER covariance matrices describing spread.

Definition at line 75 of file InitialParticleGeneration.cc.

77 {
78 //Calculate initial parameters of the HER beam before smearing
79 double E0her = pHER.E();
80 double thX0her = std::atan(pHER.Px() / pHER.Pz());
81 double thY0her = std::atan(pHER.Py() / pHER.Pz());
82
83 //Calculate initial parameters of the LER beam before smearing
84 double E0ler = pLER.E();
85 double thX0ler = std::atan(pLER.Px() / pLER.Pz());
86 double thY0ler = std::atan(pLER.Py() / pLER.Pz());
87
88 // calculate gradient of transformation (EH, thXH, thYH, EL, thXL, thYL) -> (Ecms, bX, bY, bZ, angleCmsX, angleCmsY)
89 Eigen::MatrixXd grad = getGradientMatrix(E0her, thX0her, thY0her, E0ler, thX0ler, thY0ler);
90 Eigen::VectorXd mu = getCentralValues(E0her, thX0her, thY0her, E0ler, thX0ler, thY0ler);
91
92 // calculate covariance smearing matrix in the (Ecms, bX, bY, bZ, angleCmsX, angleCmsY) coordinate system
93 Eigen::MatrixXd covT = transformCov(toEigen(adjustCovMatrix(covHER)), toEigen(adjustCovMatrix(covLER)), grad);
94
95 // return random number generator which allows to generate Lorentz transformation parameters based on the cmsEnergy
96 return ConditionalGaussGenerator(mu, covT);
97
98 }
static Eigen::MatrixXd toEigen(TMatrixDSym m)
transform matrix from ROOT to Eigen format
Eigen::MatrixXd transformCov(const Eigen::MatrixXd &covHER, const Eigen::MatrixXd &covLER, const Eigen::MatrixXd &grad)
transform covariance matrix to the variables (Ecms, boostX, boostY, boostZ, angleXZ,...
Definition: beamHelpers.h:293
Eigen::VectorXd getCentralValues(double eH, double thXH, double thYH, double eL, double thXL, double thYL)
get vector including (Ecms, boostX, boostY, boostZ, angleXZ, angleYZ) from beam energies and angles
Definition: beamHelpers.h:267
TMatrixDSym adjustCovMatrix(TMatrixDSym cov) const
adjust smearing covariance matrix based on the generation flags
Eigen::MatrixXd getGradientMatrix(double eH, double thXH, double thYH, double eL, double thXL, double thYL)
get gradient matrix of (eH, thXH, thYH, eL, thXL, thYL) -> (eCMS, boostX, boostY, boostY,...
Definition: beamHelpers.h:233

◆ initialize()

void initialize ( )

function to be executed on initialize()

Definition at line 22 of file InitialParticleGeneration.cc.

23 {
24 m_event.registerInDataStore();
25 }

◆ initProbMax() [1/12]

void initProbMax ( )

Initialize standard stream objects for probability function

Definition at line 136 of file EvtB0toKsKK.cc.

137 {
138 //setProbMax(50000.0);
139 //setProbMax(100000.0);
140 //setProbMax(200000.0);
141 //setProbMax(400000.0);
142 //setProbMax(600000.0);
143 setProbMax(1100000.0);
144 }

◆ initProbMax() [2/12]

void initProbMax ( )

Sets maximal probab.

Definition at line 56 of file EvtBCL.cc.

57 {
58
59 EvtId parnum, mesnum, lnum, nunum;
60
61 parnum = getParentId();
62 mesnum = getDaug(0);
63 lnum = getDaug(1);
64 nunum = getDaug(2);
65
66 double mymaxprob = calcamp->CalcMaxProb(parnum, mesnum, lnum, nunum, bclmodel);
67
68 setProbMax(mymaxprob);
69 }

◆ initProbMax() [3/12]

void initProbMax ( )

The function sets the maximum value of the probability.

Definition at line 59 of file EvtBSemiTauonic.cc.

60 {
61 EvtId parnum, mesnum, lnum, nunum;
62
63 parnum = getParentId();
64 mesnum = getDaug(0);
65 lnum = getDaug(1);
66 nunum = getDaug(2);
67
68 double maxprob = m_CalcAmp->CalcMaxProb(parnum, mesnum,
69 lnum, nunum,
71
72 B2INFO("EvtBSemiTauonic::initProbMax()>> maxprob: " << maxprob);
73
74 setProbMax(maxprob);
75 }
double CalcMaxProb(EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator)
The function calculates the maximum probability.

◆ initProbMax() [4/12]

void initProbMax ( )

The function sets the maximum value of the probability.


Definition at line 59 of file EvtBSemiTauonic2HDMType2.cc.

60 {
61 EvtId parnum, mesnum, lnum, nunum;
62
63 parnum = getParentId();
64 mesnum = getDaug(0);
65 lnum = getDaug(1);
66 nunum = getDaug(2);
67
68 double maxprob = m_CalcAmp->CalcMaxProb(parnum, mesnum,
69 lnum, nunum,
71
72 B2INFO("EvtBSemiTauonic2HDMType2::initProbMax()>> maxprob: " << maxprob);
73
74 setProbMax(maxprob);
75 }

◆ initProbMax() [5/12]

void initProbMax ( )

The function to sets a maximum probability.

In this model, noProbMax() is called because maximum probability is determined by m_dGdsbProbMax

Definition at line 161 of file EvtBtoXsnunu_FERMI.cc.

162 {
163 noProbMax();
164 }

◆ initProbMax() [6/12]

void initProbMax ( )

Sets the Maximum probability for the PHSP reweight.

Maximum value is hardcoded and inherited from the EvtEtaDalitz class.

Definition at line 57 of file EvtEtaFullDalitz.cc.

58 {
59
60 setProbMax(2.5);
61
62 }

◆ initProbMax() [7/12]

void initProbMax ( )

Sets the Maximum probability for the PHSP reweight.

Maximum value is hardcoded and inherited from the EvtEtaDalitz class.

Definition at line 72 of file EvtEtaPi0Dalitz.cc.

73 {
74
75 setProbMax(2.5);
76
77 }

◆ initProbMax() [8/12]

void initProbMax ( )

Sets the Maximum probability for the PHSP reweight.

Maximum value is hardcoded and inherited from the EvtEtaDalitz class.

Definition at line 74 of file EvtEtaPrimeDalitz.cc.

75 {
76
77 setProbMax(2.5);
78
79 }

◆ initProbMax() [9/12]

void initProbMax ( )

The function to sets a maximum probability.

At this model, noProbMax() is used.

Definition at line 108 of file EvtFlatDaughter.cc.

109 {
110 noProbMax();
111 }

◆ initProbMax() [10/12]

void initProbMax ( )

The function to sets a maximum probability.

Definition at line 140 of file EvtKnunu.cc.

141 {
142 // Maximum probability was obtained by 10^6 MC samples. It was 71.177
143 // About 10% higher value is used.
144 setProbMax(80.0);
145 }

◆ initProbMax() [11/12]

void initProbMax ( )

The function to sets a maximum probability.

Definition at line 169 of file EvtKstarnunu_REV.cc.

170 {
171 // Maximum probability was obtained by 10^6 MC samples. It was 20000.169
172 // About 10% higher value is used.
173 setProbMax(22000.0);
174 }

◆ initProbMax() [12/12]

void initProbMax ( )

Initialize standard stream objects for probability function

Definition at line 46 of file EvtPhspCP.cc.

47 {
48 setProbMax(2 * (getArg(3)*getArg(3) + getArg(5)*getArg(5)));
49 }

◆ Lep() [1/3]

double Lep ( const double  mtau,
int  tauhel,
double  q2,
double  costau 
) const

The function to calculate the Leptonic Amplitudes for B->Dtaunu decay of the scalar type contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
q2q^2 of the decay (square of l+nu invariant mass).
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value.

Definition at line 151 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

153 {
154 // sanity check
155 assert(chktauhel(tauhel));
156
157 if (tauhel == -1)return 0.;
158 if (tauhel == +1)return -2.*sqrt(q2) * v(mtau, q2);
159
160 // should never reach here ...
161 assert(0);
162 return 0.;
163 }
double v(double mtau, double q2) const
Function to calculate the tau velocity.

◆ Lep() [2/3]

double Lep ( const double  mtau,
int  tauhel,
int  whel,
double  q2,
double  costau 
) const

The function to calculate the Leptonic Amplitudes for B->D*taunu decay of the vector type contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
whelhelicity of the virtual vector boson {+1,0,1,2}.
q2q^2 of the decay (square of l+nu invariant mass).
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value.

Definition at line 128 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

131 {
132 // sanity check
133 assert(chktauhel(tauhel) && chkwhel(whel));
134
135 if (tauhel == -1 && whel == -1)return sqrt(2.*q2) * v(mtau, q2) * (1. - costau);
136 if (tauhel == -1 && whel == 0) return -2.*sqrt(q2) * v(mtau, q2) * sqrt(1 - costau * costau);
137 if (tauhel == -1 && whel == +1)return sqrt(2.*q2) * v(mtau, q2) * (1. + costau);
138 if (tauhel == -1 && whel == 2) return 0.;
139
140 if (tauhel == +1 && whel == -1)return -sqrt(2.) * mtau * v(mtau, q2) * sqrt(1. - costau * costau);
141 if (tauhel == +1 && whel == 0) return 2.*mtau * v(mtau, q2) * costau;
142 if (tauhel == +1 && whel == +1)return sqrt(2.) * mtau * v(mtau, q2) * sqrt(1 - costau * costau);
143 if (tauhel == +1 && whel == 2) return -2.*mtau * v(mtau, q2);
144
145 // should never reach here ...
146 assert(0);
147 return 0.;
148 }

◆ Lep() [3/3]

double Lep ( const double  mtau,
int  tauhel,
int  whel1,
int  whel2,
double  q2,
double  costau 
) const

The function to calculate the Leptonic Amplitudes for B->D*taunu decay of the tensor type contribution.

Parameters
mtaudaughter lepton mass.
tauhelhelicity of the lepton in the (l+nu) rest frame {+1,-1}.
whel1helicity of the one virtual vector boson {+1,0,1,2}.
whel2helicity of the another virtual vector boson {+1,0,1,2}.
q2q^2 of the decay (square of l+nu invariant mass).
costaucosine of the angle between D(*) meson and the lepton in the (l+nu) rest frame.
Returns
calculated amplitude value. The tensor contribution is described by two vector contributions.

Definition at line 166 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

169 {
170 // sanity check
171 assert(chktauhel(tauhel) && chkwhel(whel1) && chkwhel(whel2));
172
173 if (whel1 == whel2)return 0;
174 if (whel1 > whel2)return -Lep(mtau, tauhel, whel2, whel1, q2, costau);
175
176 if (tauhel == -1 && whel1 == -1 && whel2 == 0) return -sqrt(2.) * mtau * v(mtau, q2) * (1. - costau);
177 if (tauhel == -1 && whel1 == -1 && whel2 == +1) return 2 * mtau * v(mtau, q2) * sqrt(1. - costau * costau);
178 if (tauhel == -1 && whel1 == -1 && whel2 == 2) return Lep(mtau, -1, -1, 0, q2, costau);
179 if (tauhel == -1 && whel1 == 0 && whel2 == +1) return -sqrt(2.) * mtau * v(mtau, q2) * (1. + costau);
180 if (tauhel == -1 && whel1 == 0 && whel2 == 2) return Lep(mtau, -1, -1, +1, q2, costau);
181 if (tauhel == -1 && whel1 == +1 && whel2 == 2) return Lep(mtau, -1, 0, +1, q2, costau);
182
183 if (tauhel == +1 && whel1 == -1 && whel2 == 0) return sqrt(2.*q2) * v(mtau, q2) * sqrt(1. - costau * costau);
184 if (tauhel == +1 && whel1 == -1 && whel2 == +1) return -2.*sqrt(q2) * v(mtau, q2) * costau;
185 if (tauhel == +1 && whel1 == -1 && whel2 == 2) return Lep(mtau, +1, -1, 0, q2, costau);
186 if (tauhel == +1 && whel1 == 0 && whel2 == +1) return -Lep(mtau, +1, -1, 0, q2, costau);
187 if (tauhel == +1 && whel1 == 0 && whel2 == 2) return Lep(mtau, +1, -1, +1, q2, costau);
188 if (tauhel == +1 && whel1 == +1 && whel2 == 2) return -Lep(mtau, +1, -1, 0, q2, costau);
189
190 // should never reach here ...
191 assert(0);
192 return 0.;
193 }

◆ Lmu()

EvtVector4R Lmu ( const EvtVector4R &  p4a,
const EvtVector4R &  p4b,
const EvtVector4R &  p4c 
)

Function 4Vector Lmu.

Definition at line 390 of file EvtB0toKsKK.cc.

392 {
393 const double mc = p4c.mass();
394 const double mres = (p4a + p4b).mass();
395 const double mc2 = mc * mc;
396 const double mres2 = mres * mres;
397 const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
398
399 const double N2 =
400 sqrt(s) /
401 (sqrt(s - ((mres + mc) * (mres + mc))) * sqrt(s - ((mres - mc) * (mres - mc))));
402
403 const EvtVector4R Lmu =
404 N2 * ((p4c - (p4a + p4b)) - ((p4c + (p4a + p4b)) * (mc2 - mres2) / s));
405
406 return Lmu;
407 }

◆ mD()

double mD ( int  Dhel) const

Daughter D(*) meson mass.

kinematics

Parameters
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
Returns
daughter D(*) meson mass corresponding to the Dhel.

Definition at line 500 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

501 {
502 assert(chkDhel(Dhel));
503 double mesonMass = m_mDst;
504 if (Dhel == 2) mesonMass = m_mD;
505 assert(mesonMass >= 0.);
506 return mesonMass;
507 }

◆ Mmunu()

EvtTensor4C Mmunu ( const EvtVector4R &  p4a,
const EvtVector4R &  p4b,
const EvtVector4R &  p4c 
)

Function Tensor Mmunu.

Definition at line 468 of file EvtB0toKsKK.cc.

470 {
471 const EvtTensor4C Mmunu =
472 sqrt(3.0 / 2.0) *
473 (EvtGenFunctions::directProd(Lmu(p4a, p4b, p4c), Lmu(p4a, p4b, p4c)) +
474 (1.0 / 3.0 * gmunu_tilde(p4a, p4b, p4c)));
475
476 return Mmunu;
477 }
EvtTensor4C Mmunu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor Mmunu.
Definition: EvtB0toKsKK.cc:468

◆ Multiply()

EvtTensor4C Multiply ( const EvtTensor4C &  t1,
const EvtTensor4C &  t2 
)

Function Tensor Multiply

Definition at line 437 of file EvtB0toKsKK.cc.

439 {
440 EvtTensor4C t;
441 for (unsigned int i = 0; i < 4; i++)
442 for (unsigned int j = 0; j < 4; j++) {
443 const EvtComplex element =
444 t1.get(i, 0) * t2.get(0, j) +
445 t1.get(i, 1) * t2.get(1, j) +
446 t1.get(i, 2) * t2.get(2, j) +
447 t1.get(i, 3) * t2.get(3, j);
448 t.set(i, j, element);
449 }
450 return t;
451 }

◆ operator*() [1/3]

DualNumber operator* ( double  a,
DualNumber  b 
)
inline

multiplication of double and dual number

Definition at line 107 of file beamHelpers.h.

108 {
109 return DualNumber(a * b.x, a * b.dx);
110 }

◆ operator*() [2/3]

DualNumber operator* ( DualNumber  a,
DualNumber  b 
)
inline

multiplication of two dual numbers

Definition at line 100 of file beamHelpers.h.

101 {
102 return DualNumber(a.x * b.x, a.x * b.dx + b.x * a.dx);
103 }

◆ operator*() [3/3]

GeneralVector< T > operator* ( a,
GeneralVector< T >  b 
)

product of scalar and general vector

Definition at line 171 of file beamHelpers.h.

172 {
173 return GeneralVector<T>(a * b.el[0], a * b.el[1], a * b.el[2]);
174 }

◆ operator+() [1/3]

DualNumber operator+ ( DualNumber  a,
double  b 
)
inline

addition of dual number and double

Definition at line 58 of file beamHelpers.h.

59 {
60 return DualNumber(a.x + b, a.dx);
61 }

◆ operator+() [2/3]

DualNumber operator+ ( DualNumber  a,
DualNumber  b 
)
inline

addition of dual numbers

Definition at line 51 of file beamHelpers.h.

52 {
53 return DualNumber(a.x + b.x, a.dx + b.dx);
54 }

◆ operator+() [3/3]

GeneralVector< T > operator+ ( GeneralVector< T >  a,
GeneralVector< T >  b 
)

addition of two general vectors

Definition at line 155 of file beamHelpers.h.

156 {
157 return GeneralVector<T>(a.el[0] + b.el[0], a.el[1] + b.el[1], a.el[2] + b.el[2]);
158 }

◆ operator-() [1/3]

DualNumber operator- ( double  a,
DualNumber  b 
)
inline

subtraction of double and dual number

Definition at line 72 of file beamHelpers.h.

73 {
74 return DualNumber(a - b.x, -b.dx);
75 }

◆ operator-() [2/3]

DualNumber operator- ( DualNumber  a,
double  b 
)
inline

subtraction of dual number and double

Definition at line 65 of file beamHelpers.h.

66 {
67 return DualNumber(a.x - b, a.dx);
68 }

◆ operator-() [3/3]

DualNumber operator- ( DualNumber  a,
DualNumber  b 
)
inline

subtraction of two dual numbers

Definition at line 93 of file beamHelpers.h.

94 {
95 return DualNumber(a.x - b.x, a.dx - b.dx);
96 }

◆ operator/() [1/2]

DualNumber operator/ ( double  a,
DualNumber  b 
)
inline

division of double and dual number

Definition at line 79 of file beamHelpers.h.

80 {
81 return DualNumber(a / b.x, - a / (b.x * b.x) * b.dx);
82 }

◆ operator/() [2/2]

DualNumber operator/ ( DualNumber  a,
DualNumber  b 
)
inline

division of two dual numbers

Definition at line 86 of file beamHelpers.h.

87 {
88 return DualNumber(a.x / b.x, (a.dx * b.x - a.x * b.dx) / (b.x * b.x));
89 }

◆ p()

double p ( const double &  mab,
const double &  M,
const double &  mc 
)

Constants p

Definition at line 583 of file EvtB0toKsKK.cc.

584 {
585 const double sab = mab * mab;
586 const double M_p_mc2 = (M + mc) * (M + mc);
587 const double M_m_mc2 = (M - mc) * (M - mc);
588
589 const double p = sqrt((sab - M_p_mc2) * (sab - M_m_mc2)) / (2.0 * mab);
590
591 return p;
592 }

◆ PDstar()

double PDstar ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mtau 
)

Function calculates the polarization of D*, longitudinal/(longitudinal + transverse), in B->D*taunu decay.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter tau lepton mass for B->Dtaunu decay.

Definition at line 187 of file EvtBSemiTauonicDecayRateCalculator.cc.

188 {
189 double transverse(0), longitudinal(0);
190 transverse += Gamma(BSTD, mtau, -1, -1);
191 transverse += Gamma(BSTD, mtau, +1, -1);
192 transverse += Gamma(BSTD, mtau, -1, +1);
193 transverse += Gamma(BSTD, mtau, +1, +1);
194 longitudinal += Gamma(BSTD, mtau, -1, 0);
195 longitudinal += Gamma(BSTD, mtau, +1, 0);
196 return longitudinal / (longitudinal + transverse);
197 }

◆ pf()

double pf ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
double  mtau,
int  Dhel,
double  w 
)

Phase space factor, which is multiplied to the helicity amplitude to calculate the decay rate.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter tau lepton mass for B->Dtaunu decay.
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.

Definition at line 199 of file EvtBSemiTauonicDecayRateCalculator.cc.

201 {
202 return 1. / (2 * BSTD.getMB()) * BSTD.getMB() * BSTD.getMB() * BSTD.r(Dhel) * BSTD.r(Dhel)
203 * BSTD.v(mtau, BSTD.q2(Dhel, w)) * BSTD.v(mtau, BSTD.q2(Dhel, w))
204 * sqrt(w * w - 1) / (64 * M_PI * M_PI * M_PI);
205 }

◆ PtauD()

double PtauD ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mtau 
)

Function calculates the polarization of tau, (RH - LH)/(LH + RH), in B->Dtaunu decay.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter tau lepton mass for B->Dtaunu decay.

Definition at line 170 of file EvtBSemiTauonicDecayRateCalculator.cc.

171 {
172 double left(0), right(0);
173 left += Gamma(BSTD, mtau, -1, 2);
174 right += Gamma(BSTD, mtau, +1, 2);
175 return (right - left) / (right + left);
176 }

◆ PtauDstar()

double PtauDstar ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mtau 
)

Function calculates the polarization of tau, (RH - LH)/(LH + RH), in B->D*taunu decay.

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter tau lepton mass for B->Dtaunu decay.

Definition at line 177 of file EvtBSemiTauonicDecayRateCalculator.cc.

178 {
179 double left(0), right(0);
180 for (int Dhel = -1; Dhel <= 1; Dhel++) {
181 left += Gamma(BSTD, mtau, -1, Dhel);
182 right += Gamma(BSTD, mtau, +1, Dhel);
183 }
184 return (right - left) / (right + left);
185 }

◆ q()

double q ( const double &  mab,
const double &  ma,
const double &  mb 
)

Constants q.

Definition at line 594 of file EvtB0toKsKK.cc.

595 {
596 const double mab2 = mab * mab;
597 const double ma_p_mb2 = (ma + mb) * (ma + mb);
598 const double ma_m_mb2 = (ma - mb) * (ma - mb);
599
600 const double q = sqrt((mab2 - ma_p_mb2) * (mab2 - ma_m_mb2)) / (2.0 * mab);
601
602 return q;
603 }

◆ q2()

double q2 ( int  Dhel,
double  w 
) const

Function to calculate the q^2 of the decay (square of l+nu invariant mass).

Parameters
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
Returns
calculated q^2.

Definition at line 522 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

523 {
524 return m_mB * m_mB * qh2(Dhel, w);
525 }

◆ qh2()

double qh2 ( int  Dhel,
double  w 
) const

Function to calculate the q^2 divided by the square of parent mass (m_B^2).

Parameters
Dhelhelicity of the D(*) meson in the rest frame of the parent meson {+1,0,-1} for D* and 2 for D.
wvelocity transfer variable.
Returns
calculated q^2/m_B^2.

Definition at line 516 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

517 {
518 return 1 - 2 * r(Dhel) * w + r(Dhel) * r(Dhel);
519 }

◆ RaiseIndex()

void RaiseIndex ( EvtVector4R &  vector)

Member function RaiseIndices.

Definition at line 461 of file EvtB0toKsKK.cc.

462 {
463 vector.set(1, -vector.get(1));
464 vector.set(2, -vector.get(2));
465 vector.set(3, -vector.get(3));
466 }

◆ RaiseIndices()

EvtTensor4C RaiseIndices ( const EvtTensor4C &  t)

Function RaiseIndices

Definition at line 453 of file EvtB0toKsKK.cc.

454 {
455 const EvtTensor4C tmp = Multiply(t, EvtTensor4C::g());
456 const EvtTensor4C t_raised = Multiply(EvtTensor4C::g(), tmp);
457
458 return t_raised;
459 }
EvtTensor4C Multiply(const EvtTensor4C &t1, const EvtTensor4C &t2)
Function Tensor Multiply
Definition: EvtB0toKsKK.cc:437

◆ RGammaD()

double RGammaD ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mtau,
const double  mlep = 0.0005110 
)

Function calculates the ratio of Br(B->Dtaunu)/Br(B->Dlnu), R(D).

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter tau lepton mass for B->Dtaunu decay.
mlepdaughter lepton mass for B->Dlnu decay (default is the electron mass). Note that overall factor (GF/sqrt(2) Vcb)^2 is omitted.

Definition at line 149 of file EvtBSemiTauonicDecayRateCalculator.cc.

151 {
152 double sum(0);
153 sum += Gamma(BSTD, mtau, -1, 2);
154 sum += Gamma(BSTD, mtau, +1, 2);
155 return sum / GammaSMD(BSTD, mlep);
156 }
double GammaSMD(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mlep=0.0005110)
Function calculates the SM decay rate Gamma for Dlnu decay, integrated for w and costau and summed fo...

◆ RGammaDstar()

double RGammaDstar ( const EvtBSemiTauonicHelicityAmplitudeCalculator BSTD,
const double  mtau,
const double  mlep = 0.0005110 
)

Function calculates the ratio of Br(B->Dtaunu)/Br(B->Dlnu), R(D*).

Parameters
BSTDhelicity ampliutude calculator initialized properly by user.
mtaudaughter tau lepton mass for B->Dtaunu decay.
mlepdaughter lepton mass for B->Dlnu decay (default is the electron mass).

Definition at line 158 of file EvtBSemiTauonicDecayRateCalculator.cc.

160 {
161 double sum(0);
162 for (int Dhel = -1; Dhel <= 1; Dhel++) {
163 sum += Gamma(BSTD, mtau, -1, Dhel);
164 sum += Gamma(BSTD, mtau, +1, Dhel);
165 }
166 return sum / GammaSMDstar(BSTD, mlep);
167 }
double GammaSMDstar(const EvtBSemiTauonicHelicityAmplitudeCalculator &BSTD, const double mlep=0.0005110)
Function calculates the SM decay rate Gamma for D*lnu decay, integrated for w and costau and summed f...

◆ RotateToHelicityBasisInBoostedFrame()

EvtSpinDensity RotateToHelicityBasisInBoostedFrame ( const EvtParticle *  p,
EvtVector4R  p4boost 
)

The function calculates the rotation matrix to convert the spin basis to the helicity basis in the boosted frame.

Parameters
pa pointer to the particle
p4boosta 4 vector to specify the rest frame to boost to. The function calculate the rotation matrix from the spin basis defined in the p rest frame to the helicity basis in the rest frame of p4boost.

Definition at line 27 of file EvtBSemiTauonicAmplitude.cc.

29 {
30 // theta and phi of p momentum in p4boost rest frame
31 EvtVector4R p4Boosted = boostTo(p->getP4(), p4boost, true);
32 const double theta = acos(p4Boosted.get(3) / p4Boosted.d3mag());
33 const double phi = atan2(p4Boosted.get(2), p4Boosted.get(1));
34
35 // here p must be EvtDiracParticle (if not EvtGen will abort)
36 EvtDiracSpinor spRest[2] = {p->sp(0), p->sp(1)};
37 EvtDiracSpinor sp[2];
38 sp[0] = boostTo(spRest[0], p4boost);
39 sp[1] = boostTo(spRest[1], p4boost);
40
41 EvtDiracSpinor spplus;
42 EvtDiracSpinor spminus;
43
44 double norm;
45
46 if (EvtPDL::getStdHep(p->getId()) > 0) {
47 spplus.set(1.0, 0.0, 0.0, 0.0);
48 spminus.set(0.0, 1.0, 0.0, 0.0);
49 norm = sqrt(real(sp[0].get_spinor(0) * sp[0].get_spinor(0) + sp[0].get_spinor(1) * sp[0].get_spinor(1)));
50 } else {
51 spplus.set(0.0, 0.0, 0.0, 1.0);
52 spminus.set(0.0, 0.0, 1.0, 0.0);
53 norm = sqrt(real(sp[0].get_spinor(2) * sp[0].get_spinor(2) + sp[0].get_spinor(3) * sp[0].get_spinor(3)));
54 }
55
56 spplus.applyRotateEuler(phi, theta, -phi);
57 spminus.applyRotateEuler(phi, theta, -phi);
58
59 EvtSpinDensity R;
60 R.setDim(2);
61
62 for (int i = 0; i < 2; i++) {
63 if (EvtPDL::getStdHep(p->getId()) > 0) {
64 R.set(0, i, (spplus * sp[i]) / norm);
65 R.set(1, i, (spminus * sp[i]) / norm);
66 } else {
67 R.set(0, i, (sp[i]*spplus) / norm);
68 R.set(1, i, (sp[i]*spminus) / norm);
69 }
70 }
71
72 return R;
73
74 }
double R
typedef autogenerated by FFTW

◆ s13_max()

double s13_max ( const double &  s12,
const double &  M,
const double &  m1,
const double &  m2,
const double &  m3 
)

maximum s13

Definition at line 681 of file EvtB0toKsKK.cc.

683 {
684 const double E2star = (s12 - (m2 * m2) + (m1 * m1)) / 2.0 / sqrt(s12);
685 const double E3star = (M * M - s12 - (m3 * m3)) / 2.0 / sqrt(s12);
686
687 const double s23_max =
688 ((E2star + E3star) * (E2star + E3star)) -
689 ((sqrt((E2star * E2star) - (m1 * m1)) - sqrt((E3star * E3star) - (m3 * m3))) *
690 (sqrt((E2star * E2star) - (m1 * m1)) - sqrt((E3star * E3star) - (m3 * m3))));
691
692 return s23_max;
693 }

◆ s13_min()

double s13_min ( const double &  s12,
const double &  M,
const double &  m1,
const double &  m2,
const double &  m3 
)

minimum s13

Definition at line 666 of file EvtB0toKsKK.cc.

668 {
669 const double E2star = (s12 - (m2 * m2) + (m1 * m1)) / 2.0 / sqrt(s12);
670 const double E3star = (M * M - s12 - (m3 * m3)) / 2.0 / sqrt(s12);
671
672 const double s23_min =
673 ((E2star + E3star) * (E2star + E3star)) -
674 ((sqrt((E2star * E2star) - (m1 * m1)) + sqrt((E3star * E3star) - (m3 * m3))) *
675 (sqrt((E2star * E2star) - (m1 * m1)) + sqrt((E3star * E3star) - (m3 * m3))));
676
677 return s23_min;
678 }

◆ scaleMomenta()

std::vector< ROOT::Math::PxPyPzMVector > scaleMomenta ( const std::vector< ROOT::Math::PxPyPzMVector > &  ps,
double  C 
)
inline

Scale momenta of the particles in the vector ps with a factor C.

Result is returned as returned as a vector with the scaled momenta

Definition at line 58 of file scaleParticleEnergies.h.

59 {
60 std::vector<ROOT::Math::PxPyPzMVector> psS(ps.size());
61 for (unsigned i = 0; i < ps.size(); ++i)
62 psS[i].SetCoordinates(C * ps[i].Px(), C * ps[i].Py(), C * ps[i].Pz(), ps[i].M());
63
64 return psS;
65 }

◆ scaleParticleEnergies()

void scaleParticleEnergies ( MCParticleGraph mpg,
double  EcmsTarget 
)
inline

Scale momenta of the particles by a constant factor such that total energy is changed to EcmsTarget.

It also changes momenta of the incoming particles such that energy is conserved. Should be called in CM reference frame.

Definition at line 71 of file scaleParticleEnergies.h.

72 {
73 // scale energy of incoming particles
74 double eIn = EcmsTarget / 2;
75 double pIn = sqrt(pow(eIn, 2) - pow(mpg[0].getMass(), 2));
76
77 mpg[0].setMomentum(0, 0, pIn);
78 mpg[1].setMomentum(0, 0, -pIn);
79 mpg[0].setEnergy(eIn);
80 mpg[1].setEnergy(eIn);
81
82
83 // extract momenta of final state particles
84 std::vector<ROOT::Math::PxPyPzMVector> cmsMomenta(mpg.size() - 2);
85 for (size_t i = 2; i < mpg.size(); ++i) {
86 ROOT::Math::XYZVector p = mpg[i].getMomentum();
87 cmsMomenta[i - 2].SetCoordinates(p.X(), p.Y(), p.Z(), mpg[i].getMass());
88 }
89
90 // calculate the scaling factor in an iterative way
91 double C = 1;
92 for (int i = 0; i < 10; ++i) {
93 auto cmsMomentaScaled = scaleMomenta(cmsMomenta, C);
94 double dC = getScaleFactor(cmsMomentaScaled, EcmsTarget);
95 C *= dC;
96 if (abs(dC - 1) < 3 * std::numeric_limits<double>::epsilon()) break;
97 }
98
99 // apply scaling factor on final-state particles
100 auto cmsMomentaScaled = scaleMomenta(cmsMomenta, C);
101
102 // change momenta in the MCParticleGraph
103 for (unsigned i = 2; i < mpg.size(); ++i) {
104 mpg[i].setMomentum(cmsMomentaScaled[i - 2].Px(), cmsMomentaScaled[i - 2].Py(), cmsMomentaScaled[i - 2].Pz());
105 mpg[i].setEnergy(cmsMomentaScaled[i - 2].E());
106 }
107 }
std::vector< ROOT::Math::PxPyPzMVector > scaleMomenta(const std::vector< ROOT::Math::PxPyPzMVector > &ps, double C)
Scale momenta of the particles in the vector ps with a factor C.

◆ Smu()

EvtVector4R Smu ( const EvtVector4R &  p4a,
const EvtVector4R &  p4b,
const EvtVector4R &  p4c 
)

Function 4Vector Smu.

Definition at line 370 of file EvtB0toKsKK.cc.

372 {
373 const double ma = p4a.mass();
374 const double mb = p4b.mass();
375 const double mres = (p4a + p4b).mass();
376 const double ma2 = ma * ma;
377 const double mb2 = mb * mb;
378 const double mres2 = mres * mres;
379
380 const double N1 =
381 mres /
382 (sqrt(mres2 - ((ma + mb) * (ma + mb))) * sqrt(mres2 - ((ma - mb) * (ma - mb))));
383
384 const EvtVector4R Smu =
385 N1 * ((p4a - p4b) - ((p4a + p4b) * (ma2 - mb2) / mres2));
386
387 return Smu;
388 }

◆ sqrt() [1/2]

double sqrt ( double  a)
inline

sqrt for double

Definition at line 28 of file beamHelpers.h.

28{ return std::sqrt(a); }

◆ sqrt() [2/2]

DualNumber sqrt ( DualNumber  a)
inline

sqrt of dual number

Definition at line 114 of file beamHelpers.h.

115 {
116 return DualNumber(std::sqrt(a.x), 1. / (2 * std::sqrt(a.x)) * a.dx);
117 }

◆ tan() [1/2]

double tan ( double  a)
inline

tan for double

Definition at line 31 of file beamHelpers.h.

31{ return std::tan(a); }

◆ tan() [2/2]

DualNumber tan ( DualNumber  a)
inline

tan of dual number

Definition at line 128 of file beamHelpers.h.

129 {
130 return DualNumber(std::tan(a.x), (1.0 + std::tan(a.x) * std::tan(a.x)) * a.dx);
131 }

◆ term()

void term ( )

Terminates the generator.

Closes the internal generator.

Definition at line 200 of file CRY.cc.

201 {
202 B2RESULT("Total time simulated: " << m_cryGenerator->timeSimulated() << " seconds");
203 B2RESULT("Total number of events simulated: " << m_totalTrials);
204 }

◆ TEST_F() [1/5]

TEST_F ( EvtBCLFFTest  ,
Scalar   
)

Test calculation of scalar BCL FF.

Definition at line 50 of file evtBCLFF.cc.

51 {
52 EvtId B0 = EvtPDL::getId("B0");
53 EvtId M = EvtPDL::getId("pi+");
54 const auto mB = EvtPDL::getMeanMass(B0);
55 const auto mB2 = mB * mB;
56 const auto mM = EvtPDL::getMeanMass(M);
57 const auto mM2 = mM * mM;
58 const auto q2min = 0.0;
59 const auto q2max = mB2 + mM2 - 2 * mB * mM;
60
61 int nArguments = 8;
62 double arguments[] = {0.419, -0.495, -0.43, 0.22, 0.510, -1.700, 1.53, 4.52};
63
64 EvtBCLFF bclff(nArguments, arguments);
65
66 double fplus = 0;
67 double fzero = 0;
68
69 bclff.getscalarff(B0, M, q2min, 0, &fplus, &fzero);
70
71 ASSERT_NEAR(0.253, fplus, 0.003);
72 ASSERT_NEAR(0.253, fzero, 0.003);
73
74 bclff.getscalarff(B0, M, q2max, 0, &fplus, &fzero);
75
76 ASSERT_NEAR(7.620, fplus, 0.003);
77 ASSERT_NEAR(1.006, fzero, 0.003);
78
79
80 } // Scalar Tests

◆ TEST_F() [2/5]

TEST_F ( EvtBCLFFTest  ,
Vector   
)

Test calculation of vector BCL FF.

Definition at line 83 of file evtBCLFF.cc.

84 {
85 EvtId B0 = EvtPDL::getId("B0");
86 EvtId M = EvtPDL::getId("rho+");
87 const auto mB = EvtPDL::getMeanMass(B0);
88 const auto mB2 = mB * mB;
89 const auto mM = EvtPDL::getMeanMass(M);
90 const auto mM2 = mM * mM;
91 const auto q2min = 0.0;
92 const auto q2max = mB2 + mM2 - 2 * mB * mM;
93
94 int nArguments = 11;
95 double arguments[] = { -0.833, 1.331, 0.262, 0.394, 0.163, 0.297, 0.759, 0.465, 0.327, -0.86 , 1.802};
96
97 EvtBCLFF bclFF(nArguments, arguments);
98
99 double A0 = 0;
100 double A1 = 0;
101 double A2 = 0;
102 double V = 0;
103
104 bclFF.getvectorff(B0, M, q2min, 0, &A1, &A2, &V, &A0);
105
106 ASSERT_NEAR(0.356, A0, 0.003);
107 ASSERT_NEAR(0.262, A1, 0.003);
108 ASSERT_NEAR(0.327, V, 0.003);
109
110 bclFF.getvectorff(B0, M, q2max, 0, &A1, &A2, &V, &A0);
111
112 ASSERT_NEAR(2.123, A0, 0.003);
113 ASSERT_NEAR(0.497, A1, 0.003);
114 ASSERT_NEAR(2.014, V, 0.003);
115
116 } // Vector Tests

◆ TEST_F() [3/5]

TEST_F ( ParticleGunTest  ,
AngularDistributions   
)

Tests angular generation parameters.

Definition at line 154 of file particlegun.cc.

155 {
156 std::map<int, bool> distributions = {
157 {ParticleGun::c_fixedValue, true},
158 {ParticleGun::c_uniformDistribution, true},
159 {ParticleGun::c_uniformPtDistribution, false},
160 {ParticleGun::c_uniformCosDistribution, true},
161 {ParticleGun::c_normalDistribution, true},
162 {ParticleGun::c_normalPtDistribution, false},
163 {ParticleGun::c_normalCosDistribution, true},
164 {ParticleGun::c_discreteSpectrum, true},
165 {ParticleGun::c_inversePtDistribution, false},
166 {ParticleGun::c_polylineDistribution, true},
167 {ParticleGun::c_polylinePtDistribution, false},
168 {ParticleGun::c_polylineCosDistribution, true}
169 };
170 checkVariable("theta", distributions, parameters.thetaDist, parameters.thetaParams);
171 checkVariable("phi", distributions, parameters.phiDist, parameters.phiParams);
172 }

◆ TEST_F() [4/5]

TEST_F ( ParticleGunTest  ,
MomentumDistributions   
)

Tests momentum generation parameters.

Definition at line 112 of file particlegun.cc.

113 {
114 std::map<int, bool> distributions = {
115 {ParticleGun::c_fixedValue, true},
116 {ParticleGun::c_uniformDistribution, true},
117 {ParticleGun::c_uniformPtDistribution, true},
118 {ParticleGun::c_uniformCosDistribution, false},
119 {ParticleGun::c_normalDistribution, true},
120 {ParticleGun::c_normalPtDistribution, true},
121 {ParticleGun::c_normalCosDistribution, false},
122 {ParticleGun::c_discreteSpectrum, true},
123 {ParticleGun::c_inversePtDistribution, true},
124 {ParticleGun::c_polylineDistribution, true},
125 {ParticleGun::c_polylinePtDistribution, true},
126 {ParticleGun::c_polylineCosDistribution, false}
127 };
128 checkVariable("momentum", distributions, parameters.momentumDist, parameters.momentumParams);
129 }

◆ TEST_F() [5/5]

TEST_F ( ParticleGunTest  ,
VertexDistributions   
)

Tests vertex generation parameters.

Definition at line 132 of file particlegun.cc.

133 {
134 std::map<int, bool> distributions = {
135 {ParticleGun::c_fixedValue, true},
136 {ParticleGun::c_uniformDistribution, true},
137 {ParticleGun::c_uniformPtDistribution, false},
138 {ParticleGun::c_uniformCosDistribution, false},
139 {ParticleGun::c_normalDistribution, true},
140 {ParticleGun::c_normalPtDistribution, false},
141 {ParticleGun::c_normalCosDistribution, false},
142 {ParticleGun::c_discreteSpectrum, true},
143 {ParticleGun::c_inversePtDistribution, false},
144 {ParticleGun::c_polylineDistribution, true},
145 {ParticleGun::c_polylinePtDistribution, false},
146 {ParticleGun::c_polylineCosDistribution, false}
147 };
148 checkVariable("xvertex", distributions, parameters.xVertexDist, parameters.xVertexParams);
149 checkVariable("yvertex", distributions, parameters.yVertexDist, parameters.yVertexParams);
150 checkVariable("zvertex", distributions, parameters.zVertexDist, parameters.zVertexParams);
151 }

◆ Tmunu()

EvtTensor4C Tmunu ( const EvtVector4R &  p4a,
const EvtVector4R &  p4b,
const EvtVector4R &  p4c 
)

Function Tensor Tmunu

Definition at line 423 of file EvtB0toKsKK.cc.

425 {
426 const double mres = (p4a + p4b).mass();
427 const EvtVector4R wmu = (p4a + p4b) / mres;
428
429 const EvtTensor4C Tmunu =
430 sqrt(3.0 / 2.0) *
431 (EvtGenFunctions::directProd(Smu(p4a, p4b, p4c), Smu(p4a, p4b, p4c)) +
432 (1.0 / 3.0 * (EvtTensor4C::g() - EvtGenFunctions::directProd(wmu, wmu))));
433
434 return Tmunu;
435 }
EvtTensor4C Tmunu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor Tmunu
Definition: EvtB0toKsKK.cc:423

◆ toEigen()

static Eigen::MatrixXd toEigen ( TMatrixDSym  m)
static

transform matrix from ROOT to Eigen format

Definition at line 63 of file InitialParticleGeneration.cc.

64 {
65 Eigen::MatrixXd mEig(m.GetNrows(), m.GetNcols());
66 for (int i = 0; i < m.GetNrows(); ++i)
67 for (int j = 0; j < m.GetNcols(); ++j)
68 mEig(i, j) = m(i, j);
69
70 return mEig;
71 }

◆ transformCov()

Eigen::MatrixXd transformCov ( const Eigen::MatrixXd &  covHER,
const Eigen::MatrixXd &  covLER,
const Eigen::MatrixXd &  grad 
)
inline

transform covariance matrix to the variables (Ecms, boostX, boostY, boostZ, angleXZ, angleYZ)

Definition at line 293 of file beamHelpers.h.

294 {
295
296 Eigen::MatrixXd cov = Eigen::MatrixXd::Zero(6, 6);
297 cov.block<3, 3>(0, 0) = covHER;
298 cov.block<3, 3>(3, 3) = covLER;
299
300 Eigen::MatrixXd covNew = grad * cov * grad.transpose();
301
302 return covNew;
303
304 }

◆ umu()

EvtVector4R umu ( const EvtVector4R &  p4a,
const EvtVector4R &  p4b,
const EvtVector4R &  p4c 
)

Function 4Vector umu.

Definition at line 360 of file EvtB0toKsKK.cc.

362 {
363 const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
364
365 const EvtVector4R umu = (p4a + p4b + p4c) / sqrt(s);
366
367 return umu;
368 }

◆ updateVertex()

ROOT::Math::XYZVector updateVertex ( bool  force = false)

Update the vertex position:

  1. If there is no initial particles object generate a new one with nominal values without smearing
  2. If the BeamParameters disallow smearing of the vertex it does nothing
  3. If initial particles already exist check if the vertex smearing has already been applied. If not, apply vertex smearing if allowed.
  4. Return the shift in vertex to the previous value (or the origin if there was no previous value).

This function does not update the energies as this would possibly introduce inconsistency between values used by the generator and the values contained in the initial particles. But it is useful to smear the vertex after generation.

Parameters
forceif true the vertex will be regenerated even if vertex smearing was already applied.

Definition at line 202 of file InitialParticleGeneration.cc.

203 {
204 if (!m_beamParams.isValid()) {
205 B2FATAL("Cannot generate beam without valid BeamParameters");
206 }
207 if (!m_event) {
208 // generate a new mc initial particle without smearing except for the vertex
210 return m_event->getVertex();
211 }
212 if (!m_beamParams->hasGenerationFlags(BeamParameters::c_smearVertex) or
213 (m_event->hasGenerationFlags(BeamParameters::c_smearVertex) and not force)) {
214 // already has a smeared vertex or smearing disallowed. nothing to do
215 return {0, 0, 0};
216 }
217 // Ok, now we need to just update the vertex, make sure the parameters are up to date ...
218 if (m_beamParams.hasChanged()) {
222 }
223 // Add the smear vertex flag
224 m_event->setGenerationFlags(m_event->getGenerationFlags() | BeamParameters::c_smearVertex);
225 // And generate a vertex ...
226 auto previous = m_event->getVertex();
227 auto vtx = generateVertex(m_beamParams->getVertex(), m_beamParams->getCovVertex(), m_generateVertex);
228 m_event->setVertex(vtx);
229 // Everything done
230 return vtx - previous;
231 }

◆ v()

double v ( double  mtau,
double  q2 
) const

Function to calculate the tau velocity.

Parameters
mtaudaughter lepton mass.
q2q^2 of the decay (square of l+nu invariant mass).
Returns
calculated tau velocity.

Definition at line 510 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

511 {
512 return sqrt(1 - mtau * mtau / q2);
513 }

◆ z()

double z ( double  w) const

CLN form factor z.

Parameters
wvelocity transfer variable.
Returns
calculated form factor value.

Definition at line 450 of file EvtBSemiTauonicHelicityAmplitudeCalculator.cc.

451 {
452 return (sqrt(w + 1) - sqrt(2)) / (sqrt(w + 1) + sqrt(2));
453 }

◆ ~EvtB0toKsKK()

~EvtB0toKsKK ( )
virtual

Definition at line 29 of file EvtB0toKsKK.cc.

29{}

◆ ~EvtBCL()

~EvtBCL ( )
virtual

virtual destructor

Definition at line 31 of file EvtBCL.cc.

32 {
33 delete bclmodel;
34 bclmodel = nullptr;
35 delete calcamp;
36 calcamp = nullptr;
37 }

◆ ~EvtBSemiTauonic()

~EvtBSemiTauonic ( )
virtual

The destructor

Definition at line 34 of file EvtBSemiTauonic.cc.

35 {
36 if (m_CalcHelAmp)delete m_CalcHelAmp;
37 if (m_CalcAmp)delete m_CalcAmp;
38 }

◆ ~EvtBSemiTauonic2HDMType2()

The destructor

Definition at line 34 of file EvtBSemiTauonic2HDMType2.cc.

35 {
36 if (m_CalcHelAmp)delete m_CalcHelAmp;
37 if (m_CalcAmp)delete m_CalcAmp;
38 }

◆ ~EvtBtoXsnunu_FERMI()

~EvtBtoXsnunu_FERMI ( )
virtual

Destructor.

Definition at line 34 of file EvtBtoXsnunu_FERMI.cc.

34{}

◆ ~EvtEtaFullDalitz()

~EvtEtaFullDalitz ( )
virtual

Default Destructor.

Definition at line 24 of file EvtEtaFullDalitz.cc.

24{}

◆ ~EvtEtaPi0Dalitz()

~EvtEtaPi0Dalitz ( )
virtual

Default destructor.

Definition at line 39 of file EvtEtaPi0Dalitz.cc.

39{}

◆ ~EvtEtaPrimeDalitz()

~EvtEtaPrimeDalitz ( )
virtual

Default destructor.

Definition at line 41 of file EvtEtaPrimeDalitz.cc.

41{}

◆ ~EvtFlatDaughter()

~EvtFlatDaughter ( )
virtual

Destructor.

Definition at line 34 of file EvtFlatDaughter.cc.

34{}

◆ ~EvtKnunu()

~EvtKnunu ( )
virtual

Destructor.

Definition at line 38 of file EvtKnunu.cc.

38{ }

◆ ~EvtKstarnunu_REV()

~EvtKstarnunu_REV ( )
virtual

Destructor.

Definition at line 37 of file EvtKstarnunu_REV.cc.

37{ }

◆ ~EvtPhspCP()

~EvtPhspCP ( )
virtual

Definition at line 28 of file EvtPhspCP.cc.

28{}

Variable Documentation

◆ me

const double me = Const::electron.getMass()
static

electron mass

Definition at line 177 of file beamHelpers.h.

◆ s_evtgen

EvtGen * s_evtgen = nullptr
staticprotected

pointer to the evtgen instance

Definition at line 43 of file evtBCLFF.cc.