diff --git a/EvtGenModels/EvtBToDiBaryonlnupQCD.hh b/EvtGenModels/EvtBToDiBaryonlnupQCD.hh index 99d1376..48932ba 100644 --- a/EvtGenModels/EvtBToDiBaryonlnupQCD.hh +++ b/EvtGenModels/EvtBToDiBaryonlnupQCD.hh @@ -1,56 +1,57 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: // Copyright (C) 1998 Caltech, UCSB // // Module: EvtGen/EvtBToDiBaryonlnupQCD.hh // // Description: Class to handle semileptonic B -> Baryon Anti-baryon l nu decays // using the using form factor predictions from pQCD counting rules. Taken // from arXiv:1107.0801 // // // Modification history: // // Mark Smith July 17, 2017 Module created +// John B Oct 2018 Code simplification // //------------------------------------------------------------------------ #ifndef EVTBTODIBARYONLNUPQCD_HH #define EVTBTODIBARYONLNUPQCD_HH #include "EvtGenBase/EvtDecayAmp.hh" #include class EvtParticle; class EvtBToDiBaryonlnupQCDFF; class EvtSLDiBaryonAmp; class EvtBToDiBaryonlnupQCD: public EvtDecayAmp { public: EvtBToDiBaryonlnupQCD(); virtual ~EvtBToDiBaryonlnupQCD(); std::string getName(); EvtDecayBase* clone(); void decay(EvtParticle *p); void initProbMax(); void init(); private: - EvtBToDiBaryonlnupQCDFF* ffModel; - EvtSLDiBaryonAmp* calcAmp; + EvtBToDiBaryonlnupQCDFF* ffModel_; + EvtSLDiBaryonAmp* calcAmp_; }; #endif diff --git a/EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh b/EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh index c27b903..687f19b 100644 --- a/EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh +++ b/EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh @@ -1,65 +1,66 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: // Copyright (C) 1998 Caltech, UCSB // // Module: EvtGen/EvtBToDiBaryonlnupQCDFF.hh.hh // // Description: Module for computation of B->ppbar form factors according // to pQCD counting rules, see arXiv:1107.0801. // // Modification history: // // Mark Smith July 17, 2017 Module created -// +// John B Oct 2018 Added FormFactors class +// //------------------------------------------------------------------------ #ifndef EVTBTODIBARYONLNUPQCDFF_HH #define EVTBTODIBARYONLNUPQCDFF_HH class EvtParticle; #include class EvtBToDiBaryonlnupQCDFF { public: class FormFactors { public: FormFactors() {;} virtual ~FormFactors() {;} double F1, F2, F3, F4, F5; double G1, G2, G3, G4, G5; }; EvtBToDiBaryonlnupQCDFF(); EvtBToDiBaryonlnupQCDFF(std::vector& DParameters); void getDiracFF(EvtParticle* parent, double dibaryonMass, - EvtBToDiBaryonlnupQCDFF::FormFactors& FF); + EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const; void getRaritaFF(EvtParticle* parent, double dibaryonMass, - EvtBToDiBaryonlnupQCDFF::FormFactors& FF); + EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const; void getFF(EvtParticle* parent, double dibaryonMass, - EvtBToDiBaryonlnupQCDFF::FormFactors& FF); + EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const; private: std::vector DPars; int nDPars; }; #endif diff --git a/EvtGenModels/EvtBcVHad.hh b/EvtGenModels/EvtBcVHad.hh index dd89c21..c21fe0a 100644 --- a/EvtGenModels/EvtBcVHad.hh +++ b/EvtGenModels/EvtBcVHad.hh @@ -1,73 +1,74 @@ +//-------------------------------------------------------------------------- +// +// Environment: +// This software is part of the EvtGen package. If you use all or part +// of it, please give an appropriate acknowledgement. +// +// Copyright Information: See EvtGen/COPYRIGHT +// +// Module: EvtBcVHad.hh +// +// Description: Module to implement Bc -> psi + (n pi) + (m K) decays +// +// Modification history: +// +// A V Luchinsky Jan 29, 2013 Module created +// A V Luchinsky Apr 30, 2019 psi K_S K node added +// +//------------------------------------------------------------------------ + #ifndef EvtBcVHad_HH #define EvtBcVHad_HH -#include -#include "EvtGenBase/EvtPDL.hh" -#include "EvtGenBase/EvtParticle.hh" -#include "EvtGenBase/EvtScalarParticle.hh" - -#include "EvtGenBase/EvtPatches.hh" -#include -#include -#include -#include -#include -#include "EvtGenBase/EvtGenKine.hh" -#include "EvtGenModels/EvtTauHadnu.hh" -#include "EvtGenBase/EvtDiracSpinor.hh" -#include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtVector4C.hh" -#include "EvtGenBase/EvtTensor4C.hh" -#include "EvtGenBase/EvtIdSet.hh" -#include "EvtGenBase/EvtParser.hh" - -#include "EvtGenModels/EvtWHad.hh" -#include "EvtGenModels/EvtBCVFF2.hh" +#include "EvtGenBase/EvtDecayAmp.hh" +#include -using std::endl; -using std::fstream; -using std::ifstream; -using std::cout; -using std::string; +class EvtBCVFF2; +class EvtParticle; +class EvtWHad; -class EvtBcVHad:public EvtDecayAmp { +class EvtBcVHad : public EvtDecayAmp { public: - EvtBcVHad() : nCall( 0 ) , whichfit( 0 ) , idVector( 0 ) , out_code( 0 ) , - ffmodel( 0 ) , wcurr( 0 ) { }; + + EvtBcVHad(); virtual ~EvtBcVHad(); + std::string getName(); EvtDecayBase* clone(); void initProbMax(); void init(); - void decay(EvtParticle *p); + void decay(EvtParticle *p); + protected: - int nCall; - -// whichfit --- code of the Bc -> VW formfactor set: -// 1 - SR -// 2 - PM - int whichfit; - -// idVector --- final vector particle code - int idVector; + // Hadronic current function + EvtVector4C hardCurr(EvtParticle *root_particle) const; + +private: + // whichfit --- code of the Bc -> VW formfactor set: + // 1 - SR + // 2 - PM + int whichfit; + + // idVector --- final vector particle code + int idVector; -// out_code: code of the hadronic final state -// 1 - pi+ -// 2 - pi+ pi0 -// 3 - pi+ pi+ pi- -// 4 - 4pi -// 5 - pi+ pi+ pi- pi- pi+ -// 6 - K+ K- pi+ -// 7 - K+ pi+ pi- - int out_code; - - EvtBCVFF2 *ffmodel; - EvtWHad *wcurr; - EvtVector4C _hardCurr(EvtParticle *root_particle); + // out_code: code of the hadronic final state + // 1 - pi+ + // 2 - pi+ pi0 + // 3 - pi+ pi+ pi- + // 4 - 4pi + // 5 - pi+ pi+ pi- pi- pi+ + // 6 - K+ K- pi+ + // 7 - K+ pi+ pi- + // 8 - K_S0 K+ + int out_code; + EvtBCVFF2 *ffmodel; + EvtWHad *wcurr; }; -#endif +#endif diff --git a/EvtGenModels/EvtBsMuMuKK.hh b/EvtGenModels/EvtBsMuMuKK.hh new file mode 100644 index 0000000..5ea2c2f --- /dev/null +++ b/EvtGenModels/EvtBsMuMuKK.hh @@ -0,0 +1,90 @@ +//////////////////////////////////////////////////////////////////////////////// +// // +// Environment: // +// This software is part of the EvtGen package. // +// // +// Copyright Information: See EvtGen/COPYRIGHT // +// // +// Module: EvtBsMUMUKK.hh // +// // +// Description: Routine to implement Bs -> J/psi KK // +// // +// Modification history: // +// Veronika Chobanova, Jeremy Dalseno, Diego Martinez Santos // +// April 21, 2016 Module created LHCb collaboration // +// Marcos Romero Lamas // +// February 8, 2019 Module updated LHCb collaboration // +// // +//////////////////////////////////////////////////////////////////////////////// + +#ifndef EvtBsMuMuKK_HH +#define EvtBsMuMuKK_HH + +#include +#include "EvtGenBase/EvtDecayAmp.hh" +#include "EvtGenBase/EvtComplex.hh" + +class EvtParticle; + +class EvtBsMuMuKK : public EvtDecayAmp +{ + +public: + + EvtBsMuMuKK() {;} + + virtual std::string getName() { return "BS_MUMUKK"; } + + virtual EvtDecayBase* clone(); + + virtual void init(); + + virtual void initProbMax(); + + virtual void decay(EvtParticle *p); + +protected: + + EvtComplex Flatte(const double m0, const double m) const; + + EvtComplex GetRho(const double m0, const double m) const; + + EvtComplex Breit_Wigner(const double Gamma0, const double m0, const double m, + const int J, const double q0, const double q) const; + + double Integral(const double Gamma0, const double m0, const int JR, const int JB, + const double q0, const double M_KK_ll, const double M_KK_ul, const int fcntype) const; + + double X_J(const int J, const double q, const int isB) const; + + double Wignerd(int J, int l, int alpha, double theta) const; + + EvtComplex AngularDist(int J, int l, int alpha, double cK, double cL, double chi) const; + + EvtComplex AmpTime(const int q, const EvtComplex& gplus, const EvtComplex& gminus, + const double delta, const double lambda_abs, const double Amp, + const double phis, const int eta) const; + +private: + + double MBs, MJpsi, Mf0, Mphi, Mf2p, MKp, MKm, MK0, Mpip, Mpi0, Mmu; + double Gamma0phi, Gamma0f2p; + double kin_lower_limit, kin_upper_limit, kin_middle; + double p30Kp_mid_CMS, p30Kp_ll_CMS, p30Kp_phi_CMS, p30Kp_f2p_CMS; + double p30Jpsi_mid_CMS, p30Jpsi_ll_CMS, p30Jpsi_phi_CMS, p30Jpsi_f2p_CMS; + double int_const_NR, int_Flatte_f0, int_BW_phi, int_BW_f2p; + double f_S_NR, f_f0, f_phi, f_f2p, f_phi_0, f_phi_perp, f_f2p_0, f_f2p_perp; + double A_S_NR, A_f0, A_phi_0, A_phi_perp, A_phi_par, A_f2p_0, A_f2p_perp; + double A_f2p_par; + double delta_S_NR, delta_f0, delta_phi_0, delta_phi_perp, delta_phi_par; + double delta_f2p_0, delta_f2p_perp, delta_f2p_par; + double phis_S_NR, phis_f0, phis_phi_0, phis_phi_perp, phis_phi_par; + double phis_f2p_0, phis_f2p_perp, phis_f2p_par; + double lambda_S_NR_abs, lambda_f0_abs, lambda_phi_0_abs, lambda_phi_perp_abs; + double lambda_phi_par_abs, lambda_f2p_0_abs, lambda_f2p_perp_abs; + double lambda_f2p_par_abs; + double Gamma, deltaGamma, ctau, deltaMs; + +}; + +#endif diff --git a/EvtGenModels/EvtSLDiBaryonAmp.hh b/EvtGenModels/EvtSLDiBaryonAmp.hh index bc13bf2..4fd575a 100644 --- a/EvtGenModels/EvtSLDiBaryonAmp.hh +++ b/EvtGenModels/EvtSLDiBaryonAmp.hh @@ -1,48 +1,63 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 1998 Caltech, UCSB // // Module: EvtGen/EvtSLDiBaryonAmp.hh // // Description: // // Modification history: // // Mark Smith July 18, 2017 Created +// John B Oct 2018 Code simplification // //------------------------------------------------------------------------ #ifndef EVTSLDIBARYONAMP_HH #define EVTSLDIBARYONAMP_HH #include "EvtGenBase/EvtAmp.hh" +#include "EvtGenBase/EvtDiracSpinor.hh" +#include "EvtGenBase/EvtId.hh" +#include "EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh" +#include "EvtGenBase/EvtVector4R.hh" +#include "EvtGenBase/EvtVector4C.hh" + #include class EvtParticle; -class EvtBToDiBaryonlnupQCDFF; class EvtSLDiBaryonAmp { public: - EvtSLDiBaryonAmp(EvtBToDiBaryonlnupQCDFF*); + EvtSLDiBaryonAmp(const EvtBToDiBaryonlnupQCDFF&); - virtual ~EvtSLDiBaryonAmp() {;} - void CalcAmp(EvtParticle *parent, EvtAmp& amp) const; -private: +protected: + + int checkDibaryonParity(const EvtId& id1, const EvtId& id2, + const int J1, const int J2) const; - EvtBToDiBaryonlnupQCDFF* ffModel; + int getBaryonParity(const EvtId& id) const; + std::vector getHadronicCurrents(const EvtDiracSpinor& u, const EvtDiracSpinor& v, + const EvtVector4R& p, const EvtVector4R& gMtmTerms, + const EvtVector4R& fMtmTerms) const; + +private: + + EvtBToDiBaryonlnupQCDFF ffModel_; + }; #endif diff --git a/EvtGenModels/EvtWHad.hh b/EvtGenModels/EvtWHad.hh index a2413fd..5cc069b 100644 --- a/EvtGenModels/EvtWHad.hh +++ b/EvtGenModels/EvtWHad.hh @@ -1,63 +1,74 @@ //-------------------------------------------------------------------------- // // Environment: -// This software is part of the EvtGen package developed jointly -// for the BaBar and CLEO collaborations. If you use all or part +// This software is part of the EvtGen package. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT -// Copyright (C) 1998 Caltech, UCSB // -// Module: EvtWHad.cc +// Module: EvtWHad.hh // // Description: Routine to calculate W -> (n pi) + (m K) current // according to [Kuhn, Was, Acta.Phys.Polon B39 (2008) 147] // // Modification history: -// AVL 20 Jan, 2013 Module created +// A V Luchinsky 20 Jan, 2013 Module created // //------------------------------------------------------------------------ -// - #ifndef EvtWHad_HH #define EvtWHad_HH -#include "EvtGenBase/EvtPatches.hh" -#include "EvtGenBase/EvtVector4R.hh" +#include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtVector4C.hh" -#include "EvtGenBase/EvtTensor4C.hh" - -#include -#include -#include -#include -#include -#include -#include "EvtGenBase/EvtParticle.hh" -#include "EvtGenBase/EvtPDL.hh" -#include "EvtGenBase/EvtGenKine.hh" -#include "EvtGenModels/EvtTauHadnu.hh" -#include "EvtGenBase/EvtDiracSpinor.hh" -#include "EvtGenBase/EvtReport.hh" -#include "EvtGenBase/EvtIdSet.hh" -#include "EvtGenBase/EvtParser.hh" +#include "EvtGenBase/EvtVector4R.hh" +#include class EvtWHad { + public: - EvtVector4C WCurrent(EvtVector4R q1); - EvtVector4C WCurrent(EvtVector4R q1, EvtVector4R q2); - EvtVector4C WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3); - EvtVector4C WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3, EvtVector4R q4, EvtVector4R q5); - EvtVector4C WCurrent_KKP(EvtVector4R pKplus, EvtVector4R pKminus, EvtVector4R pPiPlus); - EvtVector4C WCurrent_KPP(EvtVector4R pKplus, EvtVector4R pPiPlus, EvtVector4R pPiMinus); + + EvtWHad(); + + EvtVector4C WCurrent(const EvtVector4R& q1) const; + + EvtVector4C WCurrent(const EvtVector4R& q1, const EvtVector4R& q2) const; + + EvtVector4C WCurrent(const EvtVector4R& q1, const EvtVector4R& q2, const EvtVector4R& q3) const; + + EvtVector4C WCurrent(const EvtVector4R& q1, const EvtVector4R& q2, const EvtVector4R& q3, + const EvtVector4R& q4, const EvtVector4R& q5) const; + + EvtVector4C WCurrent_KKP(const EvtVector4R& pKplus, const EvtVector4R& pKminus, + const EvtVector4R& pPiPlus) const; + + EvtVector4C WCurrent_KPP(const EvtVector4R& pKplus, const EvtVector4R& pPiPlus, + const EvtVector4R& pPiMinus) const; + + EvtVector4C WCurrent_KSK(const EvtVector4R& pKS, const EvtVector4R& pKplus) const; + protected: - EvtVector4C JB(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3, EvtVector4R q4, EvtVector4R q5); - EvtComplex BWa( EvtVector4R q); - EvtComplex BWf( EvtVector4R q); - EvtComplex BWr( EvtVector4R q); - double pi3G(double Q2); + + EvtVector4C JB(const EvtVector4R& q1, const EvtVector4R& q2, const EvtVector4R& q3, + const EvtVector4R& q4, const EvtVector4R& q5) const; + + EvtComplex BWa(const EvtVector4R& q) const; + + EvtComplex BWf(const EvtVector4R& q) const; + + EvtComplex BWr(const EvtVector4R& q) const; + + EvtComplex BWKK(double s, int i) const; + + double pi3G(double Q2) const; + + EvtComplex pcm(double s) const; + +private: + + std::vector mRho_, gamma0_, cK_; + }; #endif diff --git a/History.txt b/History.txt index e7fa142..ccf42ec 100644 --- a/History.txt +++ b/History.txt @@ -1,555 +1,565 @@ //========================================================================== // // History file for EvtGen // //=========================================================================== -21 & 24 Jan 2019 John Back - Correct B 4-momentum in EvtSLDiBaryonAmp for the BToDiBaryonlnupQCD model, - as well as including a parity factor for the vector amplitude terms. +3rd May 2019 John Back + Corrected EvtSLDiBaryonAmp bugs/issues in the BToDiBaryonlnupQCD model: + parity, amplitude terms and B momentum reference frame variables. + Corrected treament of spinor indices in EvtRareLb2Lll, courtesy of + Tom Blake and Michal Kreps (LHCb). + Updated the EvtBcVHad model to also handle Bc -> psi Ks K decays, + courtesy of Aleksei Luchinsky (LHCb). + Add new decay model EvtBsMuMuKK (BS_MUMUKK) for Bs to J/psi (mu+mu-) K+K-, + courtesy of Veronika Chobanova, Jeremy Dalseno, Diego Martinez Santos and + Marcos Romero Lamas (LHCb). + +15th March 2019 Tom Latham + Implement cmake build system, replacing the old config method. 7th December 2018 John Back Added the EvtBLLNuL (BLLNUL) model that generates rare B -> ell ell nu ell decays, where ell = e or mu, courtesy of Anna Danilina and Nikolai Nikitin (LHCb). Removed the EvtB2MuMuMuNu (BUTOMMMN) model, since its now replaced by the more general BLLNuL one. 5th November 2018 John Back Added the BToDiBaryonlnupQCD model for generating B to p N* l nu decays, where N can be any (exited) charged baryon (spin 1/2 or 3/2), courtesy of Mark Smith and Ryan Newcombe (LHCb), with added code optimisations. 17th October 2018 John Back Added various decay models from LHCb EvtGenExtras package: EvtBcVHad ("BC_VHAD"), Evtbs2llGammaMNT ("BSTOGLLMNT"), Evtbs2llGammaISRFSR ("BSTOGLLISRFSR"), EvtbTosllMS ("BTOSLLMS"), EvtbTosllMSExt ("BTOSLLMSEXT"), EvtLb2Baryonlnu ("Lb2Baryonlnu"), EvtLb2plnuLCSR ("Lb2plnuLCSR"), EvtLb2plnuLQCD ("Lb2plnuLQCD"), EvtFlatSqDalitz ("FLATSQDALITZ") and EvtPhspFlatLifetime ("PHSPFLATLIFETIME"). 5th October 2018 John Back Updated setupEvtGen.sh to work with the new HepForge Phabricator site. 13th March 2018 John Back Updated EvtPythiaEngine to correctly handle updates of various particle properties so that Pythia uses the same information as EvtGen (evt.pdl) for the generic and alias PYTHIA decay model. 12th March 2018 John Back Updated EvtBcXMuNu models (X = Scalar, Vector, Tensor) to generate Bc to D0(star) mu nu decays, with associated form factors in EvtBCXFF, - courtesy of Aleksei Luchinskii (LHCb). Also generalised the calculation + courtesy of Aleksei Luchinsky (LHCb). Also generalised the calculation of their maximum probabilities by reusing the CalcMaxProb method in EvtSemiLeptonicAmp, which now allows for different Q^2 binning (default remains at 25 bins). R01-07-00------------------------------------------------------------------- 13th December 2017 John Back New tag incorporating all changes below. Recommended external packages are HepMC 2.06.09, pythia 8.230, Photos++ 3.61 and Tauola++ 1.1.6c, as used in the setupEvtGen.sh script. 12th December 2017 John Back Changed Pythia validation example decay files to use Pythia8 codes. 6th December 2017 John Back Modified the examples to use DECAY.DEC (see 25th April 2016) instead of DECAY_2010.DEC. Changed EvtExternalGenList to assume Pythia8 codes are used in decay files by default, which is the case for DECAY.DEC. Also updated the setupEvtGen.sh script to work with Pythia 8.2x versions. 29th November 2017 John Back Modified EvtSVP, EvtVVP and EvtTVP models to handle both radiative and - two-lepton decays, courtesy of Aleksei Luchinskii (LHCb). + two-lepton decays, courtesy of Aleksei Luchinsky (LHCb). 14th July 2017 John Back Only create external generator objects if they don't already exist in EvtExternalGenFactory. Modified configure script to work with Pythia 8.2x 5th July 2017 Michal Kreps Register the VTOSLL model. 14th June 2017 John Back Add isNeutralKaon() boolean function and corrected comments in EvtDDalitz. 8th May 2017 Michal Kreps Fix bug in EvtbTosllVectorAmp to recognise Bs --> K*bar mu mu decay as b --> d ll transition. 8th May 2017 Michal Kreps Significantly simplify way how we decide on decay mode and daughters ordering in DDalitz model. With new code by definition all orderings of daughters in the decay file will yield same output. 4th May 2017 John Back Further fixes to DDalitz particle ordering (including charge-conjugates): Mode 5: D0 -> K- K0bar K+ and K+ K- K0bar Mode 12: D0 -> pi0 pi- pi+ and pi+ pi0 pi- Removed unneeded index ordering checks for mode 10 (D+ -> pi- pi+ pi+) and mode 11 (Ds+ -> pi- pi+ pi+) 27th April 2017 John Back Fixed DDalitz particle ordering for mode 7: D+ -> pi+ K- K+ and K+ pi+ K- and their charge-conjugates 7th April 2017 John Back Modified EvtGenExternal/EvtPythiaEngine to ensure that the EvtGen-based instances of Pythia8 (for generic and alias decays) use the same particle properties as defined by EvtGen, courtesy Patrick Robbe (LHCb). 5th April 2017 Michal Kreps Fixed indexing in copy constructor of Evt3Rank3C, which would otherwise produce an infinite loop; bug report from David Grellscheid. 3rd November 2016 John Back Modified EvtFlatQ2 model to work for all B -> X lepton lepton modes, as well as adding an extra phase space factor to correct for the dip at low q^2, courtesy of Marcin Chrzaszcz & Thomas Blake (LHCb), which is enabled by using "FLATQ2 1" instead of just "FLATQ2" in the decay file. 13th October 2016 John Back Added the TauolaCurrentOption decfile keyword to select the hadronic current in Tauola; default is the BaBar-tuned current option (int = 1). EvtParticles can store double attributes using the functions setAttributeDouble(name, double) and getAttributeDouble(name), which can be useful for storing and retrieving amplitude weights, for example. The analogous EvtParticle integer attribute interface remains unchanged: setAttribute(name, int) and getAttribute(name). 13th September 2016 John Back Modified EvtTauolaEngine to use internal Tauola spin matrices for tau pair events by temporarily setting the PDG id of the mother as a boson, keeping the same 4-momentum. The BaBar hadronic currents are now used by default. Also added the ability to change some Tauola parameters using the "Define" keyword in decay files. Added an example decay file illustrating the new features: validation/TauolaFiles/Btautau.dec 9th September 2016 Michal Kreps Reimplement code in EvtBTo3pi.F, EvtBTo3piMPP.F, EvtBTo3piP00.F and EvtBToKpipi.F in C++ in order to remove dependence on Fortran compiler. With this, there is no internal Fortran code in EvtGen. R01-06-00------------------------------------------------------------------- 1st June 2016 John Back New tag incorporating all changes below. Recommended external packages are HepMC 2.06.09, pythia 8.186, Photos++ 3.61 and Tauola++ 1.1.5. 28th April 2016 Michal Kreps For Ds+ --> 2pi+ pi- there was double counting of branching fraction resulting in total branching fraction being 1.5 times larger than measured one. Fix by revisiting submodes, which now fill total Ds --> 3pi. 25th April 2016 Michal Kreps Added DECAY.DEC/XML, which contain updated semileptonic charm and beauty branching fractions using the 2014 PDG, tuned to minimize disagreements between measurements and EvtGen for both inclusive and exclusive decays. Updated the evt.pdl particle properties file to the PDG 2014 edition. Implemented new LQCD form factors for Lb --> L mu mu from arXiv paper 1602.01399 (EvtRareLbToLllFFlQCD); old LQCD form factors are removed. 18th March 2016 John Back Fixed incorrect spinor algebra used in S1 -> 1/2 S2, 1/2 -> S3 S4 decays in EvtDiracParticle::rotateToHelicityBasis() functions, courtesy of Luis Miguel Garcia Martin and the IFIC Valencia LHCb group. 19th Feburary 2016 John Back Fixed bug in the definition of the initial spinor term Sinit in EvtRareLbToLll::HadronicAmpRS(), from Tom Blake (LHCb). 12th February 2016 John Back From LHCb, added extensions to the EvtHQET2(FF) model for semitauonic decays from Brian Hamilton, which needs a patch to EvtSemiLeptonicAmp from Jack Wimberley to ensure that the q^2 range is physical when finding the maximum amplitude probability. 2nd December 2015 John Back From LHCb, added EvtKStopizmumu model for KS -> pi0 mu mu decays based on JHEP08(1998)004, courtesy of Veronika Chobanova, Diego Martinez Santos and Jeremy Dalseno. Added EvtConst::Fermi for Fermi coupling constant. R01-05-00------------------------------------------------------------------- 21st October 2015 John Back New tag incorporating all changes below. Recommended external packages are HepMC 2.06.09, pythia 8.186, Photos++ 3.61 and Tauola++ 1.1.5. Added the EvtB2MuMuMuNu model for simulating the very rare four-leptonic decays B- -> mu+ mu- anti-nu_mu mu-, courtesy Nikolai Nikitin. 16th October 2015 John Back Updated the configure script to automatically select the library names for PHOTOS++; version 3.56 and below uses Fortran, version 3.61 and above uses C++ only (default). Avoid using v3.60, since it does not work. This needs the PHOTOS libraries built before EvtGen is configured. Modified setupEvtGen.sh to use Photos++ v3.61. 7th October 2015 John Back Updated EvtGenExternal/EvtPhotosEngine to check that additional particles from the outgoing vertex are indeed (FSR) photons, since later versions of PHOTOS introduce pair emission, where particles may not always be photons. Added the genRootDecayChain.cc validation program to create ROOT files containing information about the complete decay tree. Two example test decay files BKstarGamma.dec and BuDst0rhop.dec can be used with this; the first tests PHOTOS, the second looks at sequential decay chain storage. The plotBKstarGamma.C ROOT macro can be used for B -> K* gamma plots. 2nd October 2015 John Back Modified EvtSVPHelAmp and added a new EvtSVPHelCPMix model, implementing the complete mixing phenomenology of Bs to vector gamma decays, courtesy of Clara Remon (LHCb). EvtD0mixDalitz code: cleanup, inverted q/p for decays of D0bar (simplifies user decay files) and fixed y parameter bug, courtesy of Jordi Tico (LHCb). Changed the initialisation order of the infrared cut-off in EvtPhotosEngine. This actually has no effect, since the exponentiation function sets it to the same 1e-7 value, but it's now in the correct order if we need to update it. Removed all remaining obsolete pragma (Win32) warnings from some classes. 23rd September 2015 Michal Kreps Reimplement the real Spence function in C++ and removed its fortran implementation. 15th September 2015 Michal Kreps Fixed accessed uninitialised memory in EvtPDL.cpp, line 213. Modified the configure and setupEvtGen.sh scripts to work on Mac; needed Mac compilation patch files added to the new "platform" subdirectory. 10th September 2015 John Back Updated setupEvtGen.sh to use the recommended external packages: HepMC 2.06.09, pythia 8.186, Photos++ 3.56 and Tauola++ 1.1.5. Fixed form-factor calculations for the BTOSLLBALL model 6 used to generate b -> sll decays, courtesy of Christoph Langenbruch and David Loh (LHCb). Affects B->K*ll, B->rholl and B->omegall, particularly the electron modes. In the validation directory, added runPhotosTest.sh for testing FSR in Upsilon(4S) -> e+ e- decays, and changed the plot comparison scripts to use the 2nd directory "oldRootFiles" (which could be a soft-link) for including ROOT histograms made from a previous version of EvtGen. 27th August 2015 John Back Added Mersenne-Twister random number generator (RNG) EvtMTRandomEngine. It requires c++11 compiler features (>= gcc 4.7), which should automatically be enabled by the configure script. Introduced the preprocessor environment variable EVTGEN_CPP11 for c++11 features. EvtMTRandomEngine is the default RNG for the validation and test examples if c++11 features are enabled. Added a phase-space test validation/genPHSP.sh and PhaseSpacePlots.C to visually check the flatness of Dalitz plots in order to ensure that the RNG is not producing biased results that depend on particle ordering. Added the models EvtbsToLLLLAmp and EvtbsToLLLLHyperCP for B0_q -> l+ l- l+ l- decays (SM and one supersymmetric scenario), courtesy of Nikolai Nikitin and Konstantin Toms. Documentation provided in doc/evt_BQTOLLLL_model.pdf and doc/evt_BQTOLLLLHYPERCP_model.pdf. Changed the installation and set-up script name to be just setupEvtGen.sh; it uses the VERSION variable to specify the required tag. List of tags are available using either "svn ls -v http://svn.cern.ch/guest/evtgen/tags" or by going to http://svn.cern.ch/guest/evtgen/tags in a web browser. 12th June 2015 John Back Changed the width of chi_b1 in evt.pdl from 9.8928 GeV (!) to zero. 1st May 2015 John Back Added Bc -> scalar ell nu (EvtBcSMuNu) and Bc -> tensor ell nu (EvtBcTMuNu) decays, courtesy of Jack Wimberley, LHCb. Also included the chi_c1 mode for EvtBcVMuNu. R01-04-00------------------------------------------------------------------- 2nd April 2015 John Back Removed the EvtStdlibRandomEngine class since this can produce biases to kinematic distributions when one or more of the daughters is a resonance, such as B0 -> K pi psi (thanks to Antonio Augusto Alves Jr who discovered this issue). EvtSimpleRandomEngine is now the default random number generator in the validation and test examples. Incorporated several additions and modifications from LHCb: a) From Michal Kreps, Tom Blake & Christoph Langenbruch, added rare Lb --> Lambda^(*) ell ell models described in arXiv:1108.6129, with various form factors from Gutsche et al. (arXiv:1301.3737) and lattice QCD (arXiv:1212.4827) b) From Andrew Crocombe, implemented Bs --> K* form factors from Ball-Zwicky and z-parametrization form factors from arXiv:1006.4945 for EvtbTosllBallFF c) Christoph Langenbruch fixed the Bs -> phi ll form factors in EvtbTosllBallFF; T3 showed a non-physical pole at very low q2 which significantly affected the electron mode d) From Michal Kreps, removed semicolons from wrong places to clear warnings when compiled with the "-pedantic" option. 9th October 2014 John Back Change svnweb.cern.ch to svn.cern.ch in the setup script. 1st April 2014 John Back In EvtReport, modified the logging output severity status flags to have the "EVTGEN_" prefix, e.g. INFO becomes EVTGEN_INFO. The global report() function has been renamed to EvtGenReport(). 31st March 2014 John Back Added the ability to store named attributes for EvtParticles in the form of a map. The setAttribute(name, value) stores the required value, while getAttribute(name) retrieves the integer value. This is used in EvtPhotosEngine to specify the final-state radiation "FSR" attribute to 1 for any additional photons (EvtPhotonParticles) created by Photos++. It also stores the "ISR" attribute, but this is always set to zero, since only FSR photons are created. If the named attribute doesn't exist, then getAttribute() returns zero. 29th January 2014 Daniel Craik Removed mass assertion on GS shape in EvtDalitzReso to allow it to also be used for charged rho resonances. 27th January 2014 John Back Minor corrections to Vub models to remove further gcc 4.8 warnings. Updated configure script to work for MacOS clang (from Genser team). R01-03-00------------------------------------------------------------------- 9th January 2014 John Back New tag version "1.3.0", incorporating all changes below. Replaced auto-install script to work with this version as well as the latest versions of all external generator packages. Updated README to mention the new CERN-based web pages for Photos++ and Tauola++. 8th January 2014 John Back Fix gcc 4.6 and 4.8 compilation warnings, courtesy of Patrick Robbe (LHCb); main changes are removal of unused variables. Changed the EvtPythiaEngine class and configure script to use new Pythia 8 header locations; Pythia 8.180 or above is now required. 7th January 2014 John Back Modified EvtBCVFF to correct the Kiselev form factors from Jack Wimberley (LHCb). 9th October 2013 Daniel Craik Added Gounaris-Sakurai and Gaussian shapes to EvtGenericDalitz and set sensible defaults for LASS parameters. 19th September 2013 John Back Modified EvtGenExternal/EvtPythiaEngine to keep track of any new particles that are added to the default Pythia database to avoid duplicating particle/anti-particle entries that could override previously defined Pythia decay chains. 18th September 2013 John Back Added Mac OS flags for the configure script and src/Makefile. 15th July 2013 Daniel Craik Added flag to turn on scaling of LASS amplitude by M/q in EvtDalitzReso 15th July 2013 Daniel Craik EvtParserXML now accepts file names containing environment variables, exponential non-resonant shape in EvtDalitzReso now defined as exp(-alpha*m^2), LASS shape in EvtDalitzReso now takes a cutoff parameter 4th July 2013 Daniel Craik Added LASS, exponential non-resonant and linear non-resonant shapes to EvtGenericDalitz. 3rd July 2013 Daniel Craik Fixed auto-install script for R01-02-00. 1st July 2013 Daniel Craik Added auto-install script for R01-02-00. R01-02-00------------------------------------------------------------------- 15th May 2013 John Back New tag, version "1.2.0", incorporating all changes below. 14th May 2013 Michal Kreps Added Blatt-Weisskopf barrier factors up to L=5 in EvtGenBase/EvtBlattWeisskopf::compute(). 14th May 2013 John Back Added additional entries (appended at the end) to the evt.pdl particle data file courtesy of Romulus Godang and Belle II colleagues. 14th March 2013 John Back Added the method EvtParticle::getPDGId() to get the PDG integer for a particle directly (which just calls EvtPDL::getStdHep()). Added a check in EvtPhotosEngine::doDecay to skip Photos if a given particle has too many daughters (>= 10) to avoid a problem with a hard coded upper limit in the Photos PHOENE subroutine. 2nd February 2013 Daniel Craik Updated EvtDalitzTable to estimate probMax when it is missing from a Dalitz model. 1st February 2013 John Back Added the ability to read in Pythia 6 commands in ascii decay files in EvtDecayTable::readDecayFile (this was already possible in xml files). Modified the Photos++ engine default settings to be more suited to B decays (from LHCb defaults). 31st January 2013 John Back Added the ability to read in Pythia 8 commands in ascii decay files in EvtDecayTable::readDecayFile. They can be set using the syntax: "PythiaTypeParam module:variable=value", where Type = Generic, Alias or Both for specifying whether the parameter is for the generic or alias Pythia decay engines (or both). The 2nd argument must not contain any spaces. Fixed the list of commands strings being used in the EvtPythiaEngine class (i.e. Pythia parameters that can be set via decay files). 31st January 2013 Daniel Craik Added named parameters to various decay models. 30th January 2013 John Back Fixed some of the parameter arguments used in the EvtSVSCPiso model. 24th January 2013 John Back Set the Photos++ and Tauola++ models to use the EvtGen random number engine when useEvtGenRandom is set to true in the EvtExternalGenList constructor. 23rd January 2013 John Back Added EvtGenExternal/EvtPythiaRandom to allow the use of the EvtGen random number engine to also be used for the random engine for Pythia 8. Added a boolean (useEvtGenRandom, default = true) within the EvtExternalGenList constructor to use this feature. 18th December 2012 John Back Corrected some wrong daughter ordering assignments for decay modes 5 and 12 in EvtDDalitz. Updated validation/DalitzDecays.xml to also contain D decay mode 12, as well as various modes that may use K_S0 and K_L0. Added validation/genDDalitzModes.sh and updated validation/compareDalitz.C to do a complete comparison of the D Dalitz modes with the xml versions. 11th December 2012 Daniel Craik Updated the Xml parser to support named model parameters. Updated the generic Dalitz model to use named model parameters as an example. 15th October 2012 John Back Make EvtSimpleRandomEngine inherit from EvtRandomEngine to avoid crash in EvtGen.cpp when no random engine is defined (from Bjoern Spruck). R01-01-00------------------------------------------------------------------- 4th October 2012 John Back New tag, version "1.1.0", incorporating all changes below. Provide proper default constructors for EvtVector4R and EvtPhotonParticle. Modified the validation and test code to also compile/link in the case of no external generators being included. 3rd October 2012 John Back Corrected the t3 vector form factor values for the Ball-Zwicky '05 model (modelId = 6) in EvtbTosllBallFF::getVectorFF(), which were set to t3tilde instead. 18th September 2012 John Back Moved the external generator engines to a new sub-directory EvtGenExternal. Building the code now creates 2 libraries: libEvtGen.so (Base+Models) and libEvtGenExternal.so. This allows anyone to ignore using the new external generators if required (by not creating/loading the 2nd library). Added prefix option to the configure script/Makefile to allow the user to specify an installation directory for the include files, libraries, DECAY.DEC and evt.pdl files (for Genser). 14th September 2012 Michal Kreps Fixed the calculation of the angle between decay planes in the function EvtKine::EvtDecayAngleChi. Fixed typo in EvtLb2Lll decay model. Only some NP scenarious could be affected, SM one is definitely unaffected. 13th September 2012 John Back Added the use of the environment variables EVTGEN_PHOTOS, EVTGEN_PYTHIA and EVTGEN_TAUOLA to specify if the Photos, Pythia and/or Tauola engine classes are used or not. These variables are set by the configure script, depending if the library paths are specified for these generators. R01-00-01-------------------------------------------------------------------- 12th September 2012 John Back New tag incorporating all changes below, since R01-00-00. 11th September 2012 John Back Modified the Photos and Tauola engine classes to use the new Photospp and Tauolapp namespaces that are present in the latest versions of Photos++(3.5) and Tauola++(1.0.7). Updated the configure file to get the correct location of the Tauola++ include files. Added the D0->pi+pi-pi0 decay mode in EvtDDalitz from Marco Gersabeck and Frederic Dreyer (LHCb). Added new decay models/classes from Alexey Luchinsky (LHCb): EvtBcVMuNu, EvtTVP, EvtWnPi, EvtSVP, EvtXPsiGamma, EvtBcVNpi 29th June 2012 John Back Corrected mass(squared) variables filled in the Dalitz TTree in validation/genExampleRootFiles. 15th May 2012 Daniel Craik Updated EvtD0gammaDalitz to deal with D mesons from neutral B->DK Added save function to validation/compareDalitz.C. 11th May 2012 Daniel Craik Replaced BaBar specific configuration for BlattWeisskopf birth factors. Updated XML conversion script to handle new configuration. Fixed some bugs in the XML conversion script related to particle modifications. 9th May 2012 Daniel Craik Added latex documentation for xml decay files. 2nd May 2012 Daniel Craik Added resDaughters attribute to the Dalitz resonance xml tag to simplify defining symmetric resonances. Updated validation xml files to use the new functionality. 27th April 2012 Daniel Craik Upgraded EvtGenericDalitz to use EvtDalitzReso for resonances. Added validation to compare EvtGenericDalitz to all 11 EvtDDalitz modes. Added a root macro to quickly compare two Dalitz decays for validation. 24th April 2012 John Back Solved two bugs in the EvtD0gammaDalitz model (from Jordi Tico, LHCb): configuration of the conjugated model, and using only the B charge to determine the model used, not the D flavour. 17th April 2012 Daniel Craik Updated the GenericDalitz validation code to use the same probMax values as DDalitz. Added XML decay file parsing to EvtGen::readUDecay. Dec files are still the default. 30th March 2012 John Back Update maximum probability values in EvtDDalitz::initProbMax() for all DDalitz modes. 23rd March 2012 John Back Added the EvtEta2MuMuGamma decay model from LHCb. 21st March 2012 John Back Added EvtD0gammaDalitz decay model from LHCb. 20th March 2012 Daniel Craik Added backwards compatibility for Pythia 6 commands in the XML configuration. Updated decay file conversion tool to convert JetSetPar lines to pythia6Param tags. 19th March 2012 Daniel Craik Added infrastructure to pass commands to external generators. XML config now takes Pythia8 configuration commands. 16th March 2012 Daniel Craik Added the ability to define particles from the PDL for Dalitz decay resonances instead of defining mass, width and spin seperately. Renamed the lifetime attribute of Dalitz decay resonaces to width to avoid confusion. Added further validation code for the generic Dalitz model. 15th March 2012 Daniel Craik Added validation code for xml decay files and the generic Dalitz model. R01-00-00 ------------------------------------------------------------------ 6th March 2012 John Back First official version for Genser (evtgen 1.0.0) that includes support for external generators: Pythia8, Photos++ and Tauola++. This also includes a preliminary version of creating Dalitz plot decay models using EvtGenericDalitz. diff --git a/src/EvtGenModels/EvtBToDiBaryonlnupQCD.cpp b/src/EvtGenModels/EvtBToDiBaryonlnupQCD.cpp index 174516e..3dec9a4 100644 --- a/src/EvtGenModels/EvtBToDiBaryonlnupQCD.cpp +++ b/src/EvtGenModels/EvtBToDiBaryonlnupQCD.cpp @@ -1,253 +1,249 @@ - - // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: // Copyright (C) 1998 Caltech, UCSB // // Module: EvtBToDiBaryonlnupQCD.cc // // Description: Routine to implement B -> Baryon Anti-baryon l nu decays. // The form factors are from arXiv.1107.0801 (2011) // // // Modification history: // // Mark Smith 17/07/2017 Module created // Ryan Newcombe May 2018 Additional baryons and // Rarita-Schwinger-type particles // John B Oct 2018 Optimisations // //------------------------------------------------------------------------ #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtSpinType.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtIdSet.hh" #include "EvtGenBase/EvtScalarParticle.hh" #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenModels/EvtBToDiBaryonlnupQCD.hh" #include "EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh" #include "EvtGenModels/EvtSLDiBaryonAmp.hh" using std::endl; EvtBToDiBaryonlnupQCD::EvtBToDiBaryonlnupQCD() : - ffModel(0), - calcAmp(0) + ffModel_(0), + calcAmp_(0) { } EvtBToDiBaryonlnupQCD::~EvtBToDiBaryonlnupQCD() { - delete ffModel; - ffModel = 0; - delete calcAmp; - calcAmp = 0; + delete ffModel_; + ffModel_ = 0; + delete calcAmp_; + calcAmp_ = 0; } std::string EvtBToDiBaryonlnupQCD::getName() { return "BToDiBaryonlnupQCD"; } EvtDecayBase* EvtBToDiBaryonlnupQCD::clone() { return new EvtBToDiBaryonlnupQCD; } void EvtBToDiBaryonlnupQCD::decay(EvtParticle *p) { p->initializePhaseSpace(getNDaug(), getDaugs(), true); - calcAmp->CalcAmp(p, _amp2); + calcAmp_->CalcAmp(p, _amp2); } void EvtBToDiBaryonlnupQCD::init() { if ( !(getNArg() == 6 || getNArg() == 7) ) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected " << " 6 or 7 arguments but found:" << getNArg() << endl; EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate execution!" << endl; ::abort(); } if (getNDaug() != 4) { - EvtGenReport(EVTGEN_ERROR,"EvtGen") - << "Wrong number of daughters in EvtBToDiBaryonlnupQCD model: " - << "4 daughters expected but found: " << getNDaug() << endl; + EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Wrong number of daughters in EvtBToDiBaryonlnupQCD model: " + << "4 daughters expected but found: " << getNDaug() << endl; EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate execution!" << endl; ::abort(); } // We expect B -> baryon baryon lepton neutrino EvtSpinType::spintype parentType = EvtPDL::getSpinType(getParentId()); EvtSpinType::spintype leptonType = EvtPDL::getSpinType(getDaug(2)); EvtSpinType::spintype neutrinoType = EvtPDL::getSpinType(getDaug(3)); if (parentType != EvtSpinType::SCALAR) { - EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected " + EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected " << " a SCALAR parent, found:" << EvtPDL::name(getParentId()) << endl; EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate execution!" << endl; ::abort(); } if (leptonType != EvtSpinType::DIRAC) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected " << " a DIRAC 3rd daughter, found:" << EvtPDL::name(getDaug(2)) << endl; EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate execution!" << endl; ::abort(); } if (neutrinoType != EvtSpinType::NEUTRINO) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtBToDiBaryonlnupQCD model expected " << " a NEUTRINO 4th daughter, found:" << EvtPDL::name(getDaug(3)) << endl; EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate execution!" << endl; ::abort(); } // Get the 6 form factor D parameters from model arguments in the decay file std::vector DPars(6); for (int i = 0; i < 6; i++) { DPars[i] = getArg(i); } // Form factor model - ffModel = new EvtBToDiBaryonlnupQCDFF(DPars); + ffModel_ = new EvtBToDiBaryonlnupQCDFF(DPars); // Set amplitude calculation pointer. // Accomodate for spin 1/2 (DIRAC) or 3/2 (RARITASCHWINGER) baryons EvtSpinType::spintype baryon1Type = EvtPDL::getSpinType(getDaug(0)); EvtSpinType::spintype baryon2Type = EvtPDL::getSpinType(getDaug(1)); if ( (baryon1Type == EvtSpinType::DIRAC && baryon2Type == EvtSpinType::RARITASCHWINGER) || (baryon1Type == EvtSpinType::RARITASCHWINGER && baryon2Type == EvtSpinType::DIRAC) || (baryon1Type == EvtSpinType::DIRAC && baryon2Type == EvtSpinType::DIRAC) ) { - calcAmp = new EvtSLDiBaryonAmp(ffModel); + calcAmp_ = new EvtSLDiBaryonAmp(*ffModel_); } else { - EvtGenReport(EVTGEN_ERROR,"EvtGen") - << "Wrong baryon spin type in EvtBToDiBaryonlnupQCD model. " - << "Expected spin type " << EvtSpinType::DIRAC - << " or " << EvtSpinType::RARITASCHWINGER - << ", found spin types " << baryon1Type << " and " << baryon2Type << endl; + EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Wrong baryon spin type in EvtBToDiBaryonlnupQCD model. " + << "Expected spin type " << EvtSpinType::DIRAC + << " or " << EvtSpinType::RARITASCHWINGER + << ", found spin types " << baryon1Type << " and " << baryon2Type << endl; EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate execution!" << endl; ::abort(); } } void EvtBToDiBaryonlnupQCD::initProbMax() { // Set maximum prob using dec file parameter if present if (getNArg() == 7) { setProbMax(getArg(6)); } else { // Default probability for the B -> p p l nu mode, where l = e, mu or tau setProbMax(3.0e6); // Specific decay modes, where we have one proton plus a second // baryon that can be any (excited) state. They all have lower // maximum probabilities compared to the default pp mode in order // to improve accept/reject generation efficiency static EvtIdSet BMesons("B-", "B+"); static EvtIdSet Delta("Delta+", "anti-Delta-"); static EvtIdSet LambdaC("Lambda_c+", "anti-Lambda_c-"); static EvtIdSet LambdaC1("Lambda_c(2593)+", "anti-Lambda_c(2593)-"); static EvtIdSet LambdaC2("Lambda_c(2625)+", "anti-Lambda_c(2625)-"); static EvtIdSet N1440("N(1440)+", "anti-N(1440)-"); static EvtIdSet N1520("N(1520)+", "anti-N(1520)-"); static EvtIdSet N1535("N(1535)+", "anti-N(1535)-"); static EvtIdSet N1650("N(1650)+", "anti-N(1650)-"); static EvtIdSet N1700("N(1700)+", "anti-N(1700)-"); static EvtIdSet N1710("N(1710)+", "anti-N(1710)-"); static EvtIdSet N1720("N(1720)+", "anti-N(1720)-"); EvtId parId = getParentId(); EvtId bar1Id = getDaug(0); EvtId bar2Id = getDaug(1); // These probabilties are sensitive to the sub-decay modes of the excited baryon states, // which limit the available phase space and allows for events to be generated within the // 10,000 event trial limit. Otherwise the amplitude varies too much (by more than a factor // of a million) and events fail to be generated correctly. In case of problems, specify // the maximum probability by passing an extra 7th model parameter if (BMesons.contains(parId)) { if (Delta.contains(bar1Id) || Delta.contains(bar2Id)) { // Delta - setProbMax(2.2e6); + setProbMax(1e7); } else if (LambdaC.contains(bar1Id) || LambdaC.contains(bar2Id)) { // Lambda_c+ setProbMax(1000.0); } else if (LambdaC1.contains(bar1Id) || LambdaC1.contains(bar2Id)) { // Lambda_c+(2593) setProbMax(200.0); } else if (LambdaC2.contains(bar1Id) || LambdaC2.contains(bar2Id)) { // Lambda_c+(2625) setProbMax(500.0); } else if (N1440.contains(bar1Id) || N1440.contains(bar2Id)) { // N(1440) - setProbMax(1.25e6); + setProbMax(8e5); } else if (N1520.contains(bar1Id) || N1520.contains(bar2Id)) { // N(1520) - setProbMax(1.25e6); + setProbMax(8e6); } else if (N1535.contains(bar1Id) || N1535.contains(bar2Id)) { // N(1535) - setProbMax(1.25e6); + setProbMax(8e5); } else if (N1650.contains(bar1Id) || N1650.contains(bar2Id)) { // N(1650) - setProbMax(1.25e6); + setProbMax(8e5); } else if (N1700.contains(bar1Id) || N1700.contains(bar2Id)) { // N(1700) - setProbMax(1.25e6); + setProbMax(4e6); } else if (N1710.contains(bar1Id) || N1710.contains(bar2Id)) { // N(1710) - setProbMax(1.25e6); + setProbMax(5e5); } else if (N1720.contains(bar1Id) || N1720.contains(bar2Id)) { // N(1720) - setProbMax(1.25e6); + setProbMax(4e6); } // Baryon combinations } // B parent } // Specific modes } diff --git a/src/EvtGenModels/EvtBToDiBaryonlnupQCDFF.cpp b/src/EvtGenModels/EvtBToDiBaryonlnupQCDFF.cpp index 2e73c5e..45036fc 100644 --- a/src/EvtGenModels/EvtBToDiBaryonlnupQCDFF.cpp +++ b/src/EvtGenModels/EvtBToDiBaryonlnupQCDFF.cpp @@ -1,89 +1,89 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: // Copyright (C) 1998 Caltech, UCSB // // Module: EvtBTopplnupQCDFF.cc // // Description: Routine to implement form factor calculation for // B->Baryon Anti-baryon l nu from pQCD counting rules. // Taken from arXiv:1107.0801 // // // Modification history: // // Mark Smith 17/07/2017 Module created // Ryan Newcombe May 2018 Added function to get form factors for // Rarita-Schwinger daughters -// John B Oct 2018 Simplified code +// John B Oct 2018 Added FormFactors class // //-------------------------------------------------------------------------- #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh" EvtBToDiBaryonlnupQCDFF::EvtBToDiBaryonlnupQCDFF() : DPars(), nDPars(0) { } EvtBToDiBaryonlnupQCDFF::EvtBToDiBaryonlnupQCDFF(std::vector& DParameters): DPars(DParameters), nDPars(DParameters.size()) { } void EvtBToDiBaryonlnupQCDFF::getFF(EvtParticle*, double dibaryonMass, - EvtBToDiBaryonlnupQCDFF::FormFactors& FF) + EvtBToDiBaryonlnupQCDFF::FormFactors& FF) const { if (nDPars == 6 && dibaryonMass > 0.0) { // 5/3*[1/M^2]^3 double t = 5.0/(3.0*pow(dibaryonMass, 6.0)); double Dp = DPars[0]; double Dpb = DPars[1]; double D2 = DPars[2]; double D3 = DPars[3]; double D4 = DPars[4]; double D5 = DPars[5]; FF.F1 = (Dp + 0.2*Dpb)*t; FF.F2 = -D2*t; FF.F3 = -D3*t; FF.F4 = -D4*t; FF.F5 = -D5*t; FF.G1 = (Dp - 0.2*Dpb)*t; FF.G2 = -FF.F2; FF.G3 = -FF.F3; FF.G4 = -FF.F4; FF.G5 = -FF.F5; } } void EvtBToDiBaryonlnupQCDFF::getDiracFF(EvtParticle* parent, double dibaryonMass, - EvtBToDiBaryonlnupQCDFF::FormFactors &FF) + EvtBToDiBaryonlnupQCDFF::FormFactors &FF) const { this->getFF(parent, dibaryonMass, FF); } void EvtBToDiBaryonlnupQCDFF::getRaritaFF(EvtParticle* parent, double dibaryonMass, - EvtBToDiBaryonlnupQCDFF::FormFactors &FF) + EvtBToDiBaryonlnupQCDFF::FormFactors &FF) const { this->getFF(parent, dibaryonMass, FF); } diff --git a/src/EvtGenModels/EvtBcVHad.cpp b/src/EvtGenModels/EvtBcVHad.cpp index 3801980..3351730 100644 --- a/src/EvtGenModels/EvtBcVHad.cpp +++ b/src/EvtGenModels/EvtBcVHad.cpp @@ -1,215 +1,276 @@ //-------------------------------------------------------------------------- // // Environment: -// This software is part of the EvtGen package developed jointly -// for the BaBar and CLEO collaborations. If you use all or part +// This software is part of the EvtGen package. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT -// Copyright (C) 1998 Caltech, UCSB // // Module: EvtBcVHad.cc // // Description: Module to implement Bc -> psi + (n pi) + (m K) decays // // Modification history: // -// AVL Jan 29, 2013 Module created +// A V Luchinsky Jan 29, 2013 Module created +// A V Luchinsky Apr 30, 2019 psi K_S K mode added // //------------------------------------------------------------------------ -// + #include "EvtGenModels/EvtBcVHad.hh" +#include "EvtGenBase/EvtIdSet.hh" +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtReport.hh" +#include "EvtGenBase/EvtSpinType.hh" +#include "EvtGenBase/EvtTensor4C.hh" +#include "EvtGenBase/EvtVector4R.hh" +#include "EvtGenModels/EvtBCVFF2.hh" +#include "EvtGenModels/EvtWHad.hh" -EvtBcVHad::~EvtBcVHad() {} +#include + +EvtBcVHad::EvtBcVHad() : + whichfit(0), + idVector(0), + out_code(0), + ffmodel(0), + wcurr(0) +{ +} -std::string EvtBcVHad::getName(){ return "BC_VHAD";} +EvtBcVHad::~EvtBcVHad() +{ -EvtDecayBase* EvtBcVHad::clone() { return new EvtBcVHad;} + if (ffmodel) {delete ffmodel;} + if (wcurr) {delete wcurr;} + +} + +std::string EvtBcVHad::getName() { + return "BC_VHAD"; +} + +EvtDecayBase* EvtBcVHad::clone() { + return new EvtBcVHad; +} //====================================================== -void EvtBcVHad::init(){ - + +void EvtBcVHad::init() { + checkNArg(1); checkSpinParent(EvtSpinType::SCALAR); - checkSpinDaughter(0,EvtSpinType::VECTOR); - for (int i=1; i<=(getNDaug()-1);i++) { - checkSpinDaughter(i,EvtSpinType::SCALAR); - }; - - if(getNDaug()<2 || getNDaug()>6) { - EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Have not yet implemented this final state in BcVNpi model" << endl; - EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Ndaug="< 6) { + EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Have not yet implemented this final state in BC_VHAD model" << std::endl; + EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Ndaug=" << getNDaug() << std::endl; + for (int id = 0; id < (getNDaug() - 1); id++) { + EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + } + return; + } + + idVector = getDaug(0).getId(); + whichfit = int(getArg(0) + 0.1); + ffmodel = new EvtBCVFF2(idVector, whichfit); - idVector = getDaug(0).getId(); - whichfit = int(getArg(0)+0.1); - ffmodel = new EvtBCVFF2(idVector,whichfit); - wcurr = new EvtWHad(); - -// determine the code of final hadronic state - EvtIdSet thePis("pi+","pi-","pi0"); - EvtIdSet theK("K+","K-"); - if( getNDaug()==2 && thePis.contains(getDaug(1)) ) - out_code=1; - else if( getNDaug()==3 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) ) - out_code=2; - else if( getNDaug()==4 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) && thePis.contains(getDaug(3)) ) - out_code=3; - else if( getNDaug()==5 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) && thePis.contains(getDaug(3)) && thePis.contains(getDaug(4)) ) - out_code=4; - else if( getNDaug()==6 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) - && thePis.contains(getDaug(3)) && thePis.contains(getDaug(4)) && thePis.contains(getDaug(5)) ) - out_code=5; - else if( getNDaug()==4 && theK.contains(getDaug(1)) && theK.contains(getDaug(2)) && thePis.contains(getDaug(3)) ) - out_code=6; - else if( getNDaug()==4 && theK.contains(getDaug(1)) && thePis.contains(getDaug(2)) && thePis.contains(getDaug(3)) ) - out_code=7; - -// cout<<"[AVL] BcVHad: whichfit ="<WCurrent( root_particle->getDaug(1)->getP4() ); - } - else if( out_code==2) { // pi+ pi0 - hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() , - root_particle->getDaug(2)->getP4() - ); - } - else if( out_code==3) {// pi+ pi+ pi- - hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() , - root_particle->getDaug(2)->getP4(), - root_particle->getDaug(3)->getP4() - ); - } - else if( out_code == 5) // Bc -> psi pi+ pi+ pi- pi- pi+ from [Kuhn, Was, hep-ph/0602162 - { - hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(), - root_particle->getDaug(2)->getP4(), - root_particle->getDaug(3)->getP4(), - root_particle->getDaug(4)->getP4(), - root_particle->getDaug(5)->getP4() - ); - } - else if( out_code==6) { // K+ K- pi+ - hardCur =wcurr->WCurrent_KKP( root_particle->getDaug(1)->getP4() , - root_particle->getDaug(2)->getP4(), - root_particle->getDaug(3)->getP4() - ); - } - else if( out_code==7) {//K+ pi+ pi- - hardCur = wcurr->WCurrent_KPP( root_particle->getDaug(1)->getP4() , // K+ - root_particle->getDaug(2)->getP4(), // pi+ - root_particle->getDaug(3)->getP4() // pi- - ); + + if (out_code == 1) { + + // pi+ + hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4()); + + } else if (out_code == 2) { + + // pi+ pi0 + hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(), + root_particle->getDaug(2)->getP4()); + + } else if (out_code == 3) { + + // pi+ pi+ pi- + hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(), + root_particle->getDaug(2)->getP4(), + root_particle->getDaug(3)->getP4()); + + } else if (out_code == 5) { + + // Bc -> psi pi+ pi+ pi- pi- pi+ from Kuhn, Was, hep-ph/0602162 + hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(), + root_particle->getDaug(2)->getP4(), + root_particle->getDaug(3)->getP4(), + root_particle->getDaug(4)->getP4(), + root_particle->getDaug(5)->getP4()); + + } else if (out_code == 6) { + + // K+ K- pi+ + hardCur = wcurr->WCurrent_KKP(root_particle->getDaug(1)->getP4(), + root_particle->getDaug(2)->getP4(), + root_particle->getDaug(3)->getP4()); + + } else if (out_code == 7) { + + // K+ pi+ pi- + hardCur = wcurr->WCurrent_KPP(root_particle->getDaug(1)->getP4(), + root_particle->getDaug(2)->getP4(), + root_particle->getDaug(3)->getP4()); + + } else if (out_code == 8) { + + // K_S0 K+ + hardCur = wcurr->WCurrent_KSK(root_particle->getDaug(1)->getP4(), + root_particle->getDaug(2)->getP4()); + + } else { + + EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Have not yet implemented this final state in BC_VHAD model" << std::endl; + EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Ndaug=" << getNDaug() << std::endl; + for (int id = 0; id < (getNDaug() - 1); id++) { + EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + } + ::abort(); + } - else { - EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl; - EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Ndaug="<initializePhaseSpace(getNDaug(),getDaugs()); +void EvtBcVHad::decay(EvtParticle *root_particle) { - // Calculate hadronic current - EvtVector4C hardCur=_hardCurr(root_particle); + root_particle->initializePhaseSpace(getNDaug(), getDaugs()); + // Calculate hadronic current + EvtVector4C hardCur = hardCurr(root_particle); + + EvtParticle* Jpsi = root_particle->getDaug(0); + EvtVector4R - p4b(root_particle->mass(), 0., 0., 0.), // Bc momentum - p4meson=root_particle->getDaug(0)->getP4(), // J/psi momenta - Q=p4b-p4meson; - double Q2=Q.mass2(); - - - -// calculate Bc -> V W form-factors - double a1f, a2f, vf, a0f; - double m_meson = root_particle->getDaug(0)->mass(); - double m_b = root_particle->mass(); - ffmodel->getvectorff(root_particle->getId(), - root_particle->getDaug(0)->getId(), - Q2, - m_meson, - &a1f, - &a2f, - &vf, - &a0f); - double a3f = ((m_b+m_meson)/(2.0*m_meson))*a1f - - ((m_b-m_meson)/(2.0*m_meson))*a2f; - -// calculate Bc -> V W current - EvtTensor4C H; - H = a1f*(m_b+m_meson)*EvtTensor4C::g(); - H.addDirProd((-a2f/(m_b+m_meson))*p4b,p4b+p4meson); - H+=EvtComplex(0.0,vf/(m_b+m_meson))*dual(EvtGenFunctions::directProd(p4meson+p4b,p4b-p4meson)); - H.addDirProd((a0f-a3f)*2.0*(m_meson/Q2)*p4b,p4b-p4meson); - EvtVector4C Heps=H.cont2(hardCur); - - for(int i=0; i<4; i++) { - EvtVector4C eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector - EvtComplex amp=eps*Heps; - vertex(i,amp); - }; + p4b(root_particle->mass(), 0., 0., 0.), // Bc momentum + p4meson = Jpsi->getP4(), // J/psi momenta + Q = p4b - p4meson, + p4Sum = p4meson + p4b; + double Q2 = Q.mass2(); + + // Calculate Bc -> V W form-factors + double a1f(0.0), a2f(0.0), vf(0.0), a0f(0.0); + + double m_meson = Jpsi->mass(); + double m_b = root_particle->mass(); + double mVar = m_b + m_meson; + + ffmodel->getvectorff(root_particle->getId(), + Jpsi->getId(), + Q2, m_meson, &a1f, &a2f, &vf, &a0f); + + double a3f = (mVar/(2.0*m_meson))*a1f - ((m_b - m_meson)/(2.0*m_meson))*a2f; + + // Calculate Bc -> V W current + EvtTensor4C H = a1f*mVar*EvtTensor4C::g(); + H.addDirProd((-a2f/mVar)*p4b, p4Sum); + H += EvtComplex(0.0, vf/mVar)*dual(EvtGenFunctions::directProd(p4Sum, Q)); + H.addDirProd((a0f - a3f)*2.0*(m_meson/Q2)*p4b, Q); + EvtVector4C Heps = H.cont2(hardCur); + + for (int i = 0; i < 4; i++) { + EvtVector4C eps = Jpsi->epsParent(i).conj(); // psi-meson polarization vector + EvtComplex amp = eps*Heps; + vertex(i, amp); + } } diff --git a/src/EvtGenModels/EvtBsMuMuKK.cpp b/src/EvtGenModels/EvtBsMuMuKK.cpp new file mode 100644 index 0000000..14f6c29 --- /dev/null +++ b/src/EvtGenModels/EvtBsMuMuKK.cpp @@ -0,0 +1,643 @@ +//////////////////////////////////////////////////////////////////////////////// +// // +// Environment: // +// This software is part of the EvtGen package. // +// // +// Copyright Information: See EvtGen/COPYRIGHT // +// // +// Module: EvtBsMUMUKK.cc // +// // +// Description: Routine to implement Bs -> J/psi KK // +// // +// Modification history: // +// Veronika Chobanova, Jeremy Dalseno, Diego Martinez Santos // +// April 21, 2016 Module created LHCb collaboration // +// Marcos Romero Lamas // +// February 8, 2019 Module updated LHCb collaboration // +// // +//////////////////////////////////////////////////////////////////////////////// + +#include "EvtGenBase/EvtPatches.hh" + +#include "EvtGenBase/EvtGenKine.hh" +#include "EvtGenBase/EvtConst.hh" +#include "EvtGenBase/EvtCPUtil.hh" +#include "EvtGenBase/EvtdFunction.hh" +#include "EvtGenBase/EvtKine.hh" +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtRandom.hh" +#include "EvtGenBase/EvtReport.hh" +#include "EvtGenBase/EvtVector3R.hh" +#include "EvtGenBase/EvtVector4R.hh" +#include "EvtGenBase/EvtVector4C.hh" + +#include "EvtGenModels/EvtBsMuMuKK.hh" + +const double pi = EvtConst::pi; +const EvtComplex I = EvtComplex(0.0,1.0); +const double sq2 = sqrt(2.0); + +EvtDecayBase* EvtBsMuMuKK::clone() +{ + return new EvtBsMuMuKK; +} + +void EvtBsMuMuKK::init() +{ + + // DecFile parameters + checkNArg(37); + + // Non-resonant S wave + f_S_NR = getArg(0); + delta_S_NR = getArg(1); + phis_S_NR = getArg(2); + lambda_S_NR_abs = getArg(3); + + // f0 (S wave) + f_f0 = getArg(4); + delta_f0 = getArg(5); + phis_f0 = getArg(6); + lambda_f0_abs = getArg(7); + + // phi (P wave) + f_phi = getArg(8); + f_phi_0 = getArg(9); + delta_phi_0 = getArg(10); + phis_phi_0 = getArg(11); + lambda_phi_0_abs = getArg(12); + f_phi_perp = getArg(13); + delta_phi_perp = pi - getArg(14); + phis_phi_perp = getArg(15); + lambda_phi_perp_abs = getArg(16); + delta_phi_par = pi - getArg(17); + phis_phi_par = getArg(18); + lambda_phi_par_abs = getArg(19); + + // f2' (D wave) + f_f2p_0 = getArg(20); + delta_f2p_0 = getArg(21); + phis_f2p_0 = getArg(22); + lambda_f2p_0_abs = getArg(23); + f_f2p_perp = getArg(24); + delta_f2p_perp = pi - getArg(25); + phis_f2p_perp = getArg(26); + lambda_f2p_perp_abs = getArg(27); + delta_f2p_par = pi - getArg(28); + phis_f2p_par = getArg(29); + lambda_f2p_par_abs = getArg(30); + + // Time dependence + Gamma = getArg(31); + deltaGamma = getArg(32); + deltaMs = getArg(33); + + // mKK window + Mf0 = getArg(34); + kin_lower_limit = getArg(35); // the minimum is approx 2.03*MKp + kin_upper_limit = getArg(36); + + // PDG masses + MBs = EvtPDL::getMass(EvtPDL::getId("B_s0")); + MJpsi = EvtPDL::getMeanMass(EvtPDL::getId("J/psi")); + Mphi = EvtPDL::getMeanMass(EvtPDL::getId("phi")); + Mf2p = EvtPDL::getMeanMass(EvtPDL::getId("f'_2")); + MKp = EvtPDL::getMass(EvtPDL::getId("K+")); + MKm = EvtPDL::getMass(EvtPDL::getId("K-")); + MK0 = EvtPDL::getMass(EvtPDL::getId("K0")); + Mpip = EvtPDL::getMass(EvtPDL::getId("pi+")); + Mpi0 = EvtPDL::getMass(EvtPDL::getId("pi0")); + Mmu = EvtPDL::getMass(EvtPDL::getId("mu+")); + + double MBsSq = MBs*MBs; + + // Amplitudes and other time parameters + A_S_NR = sqrt(f_S_NR); + A_f0 = sqrt(f_f0); + + A_phi_0 = sqrt(f_phi_0*f_phi); + A_phi_perp = sqrt(f_phi_perp*f_phi); + // Use fabs to make sure subtractions are >= 0, since subtracting 0 from 0 can give -0 + A_phi_par = sqrt(fabs(f_phi-A_phi_perp*A_phi_perp-A_phi_0*A_phi_0)); + + f_f2p = fabs(1.0 - f_S_NR - f_f0 - f_phi); + A_f2p_0 = sqrt(f_f2p_0*f_f2p); + A_f2p_perp = sqrt(f_f2p_perp*f_f2p); + A_f2p_par = sqrt(fabs(f_f2p-A_f2p_perp*A_f2p_perp-A_f2p_0*A_f2p_0)); + + ctau = 1.0/Gamma; + Gamma0phi = EvtPDL::getWidth(EvtPDL::getId("phi")); + Gamma0f2p = EvtPDL::getWidth(EvtPDL::getId("f'_2")); + + kin_middle = 0.5*(kin_upper_limit+kin_lower_limit); + + int_const_NR = sqrt(Integral(1.0, 1.0, 0, 1, 1.0, kin_lower_limit, kin_upper_limit, 0)); + + int_Flatte_f0 = sqrt(Integral(1.0, Mf0, 0, 1, 1.0, kin_lower_limit, kin_upper_limit, 1)); + + p30Kp_mid_CMS = sqrt((pow(kin_middle,2) - pow(MKp+MKm,2)) * (pow(kin_middle,2) - pow(MKp-MKm,2)))/(2.0*kin_middle); + + p30Kp_ll_CMS = sqrt((pow(kin_lower_limit,2) - pow(MKp+MKm,2)) * (pow(kin_lower_limit,2) - pow(MKp-MKm,2)))/(2.0*kin_lower_limit); + + p30Kp_phi_CMS = sqrt((Mphi*Mphi - pow(MKp+MKm,2)) * (Mphi*Mphi - pow(MKp-MKm,2)))/(2.0*Mphi); + + p30Kp_f2p_CMS = sqrt((Mf2p*Mf2p - pow(MKp+MKm,2)) * (Mf2p*Mf2p - pow(MKp-MKm,2)))/(2.0*Mf2p); + + p30Jpsi_mid_CMS = sqrt((MBsSq - pow(kin_middle+MJpsi,2)) * (MBsSq - pow(kin_middle-MJpsi,2)))/(2.0*MBs); + + p30Jpsi_ll_CMS = sqrt((MBsSq - pow(kin_lower_limit+MJpsi,2)) * (MBsSq - pow(kin_lower_limit-MJpsi,2)))/(2.0*MBs); + + p30Jpsi_phi_CMS = sqrt((MBsSq - pow(Mphi+MJpsi,2)) * (MBsSq - pow(Mphi-MJpsi,2)))/(2.0*MBs); + + p30Jpsi_f2p_CMS = sqrt((MBsSq - pow(Mf2p+MJpsi,2)) * (MBsSq - pow(Mf2p-MJpsi,2)))/(2.0*MBs); + + int_BW_phi = sqrt(Integral(Gamma0phi, Mphi, 1, 0, p30Kp_phi_CMS, kin_lower_limit, kin_upper_limit, 2)); + + int_BW_f2p = sqrt(Integral(Gamma0f2p, Mf2p, 2, 1, p30Kp_f2p_CMS, kin_lower_limit, kin_upper_limit, 2)); + + // 4 daughters + checkNDaug(4); + + // Spin-0 parent + checkSpinParent(EvtSpinType::SCALAR); // B_s0 (anti-B_s0) + + // Daughters + checkSpinDaughter(0,EvtSpinType::DIRAC); // mu+ (mu-) + checkSpinDaughter(1,EvtSpinType::DIRAC); // mu- (mu+) + checkSpinDaughter(2,EvtSpinType::SCALAR); // K+ (K-) + checkSpinDaughter(3,EvtSpinType::SCALAR); // K- (K+) + + // B_s0 parent (Parent must be B_s0 or anti-B_s0) + const EvtId p = getParentId(); + if (p != EvtPDL::getId("B_s0") && p != EvtPDL::getId("anti-B_s0")) { + assert(0); + } + + // Daughter types and ordering (should be mu+-, mu-+, K+-, K-+) + const EvtId d1 = getDaug(0); + const EvtId d2 = getDaug(1); + const EvtId d3 = getDaug(2); + const EvtId d4 = getDaug(3); + if (!((d1 == EvtPDL::getId("mu+") || d1 == EvtPDL::getId("mu-")) && + (d2 == EvtPDL::getId("mu-") || d2 == EvtPDL::getId("mu+")) && + (d3 == EvtPDL::getId("K+") || d3 == EvtPDL::getId("K-")) && + (d4 == EvtPDL::getId("K-") || d4 == EvtPDL::getId("K+")) )) { + assert(0); + } + +} + +// Get ProbMax +void EvtBsMuMuKK::initProbMax() +{ + const EvtComplex term11 = sqrt(p30Jpsi_f2p_CMS*p30Kp_f2p_CMS); + + const EvtComplex term12 = X_J(2, p30Kp_f2p_CMS, 0)*X_J(1, p30Jpsi_f2p_CMS, 1)* + p30Kp_f2p_CMS*p30Kp_f2p_CMS*p30Jpsi_f2p_CMS* + (A_f2p_0 + 0.3*A_f2p_perp + 0.3*A_f2p_par); + + const EvtComplex term13 = f_f2p*Breit_Wigner(Gamma0f2p, Mf2p, Mf2p, 2, + p30Kp_f2p_CMS, p30Kp_f2p_CMS)/int_BW_f2p; + + const EvtComplex term21 = sqrt(p30Jpsi_phi_CMS*p30Kp_phi_CMS); + + const EvtComplex term22 = X_J(1, p30Kp_phi_CMS, 0)* + p30Kp_phi_CMS* + (0.65*A_phi_0 + 0.6*A_phi_perp + 0.6*A_phi_par); + + const EvtComplex term23 = f_phi*Breit_Wigner(Gamma0phi, Mphi, Mphi, 1, + p30Kp_phi_CMS, p30Kp_phi_CMS)/int_BW_phi; + + const EvtComplex term31 = sqrt(p30Jpsi_ll_CMS*p30Kp_ll_CMS); + + const EvtComplex term32 = X_J(1, p30Jpsi_ll_CMS, 1)*p30Jpsi_ll_CMS; + + const EvtComplex term33 = f_f0*Flatte(Mf0, kin_lower_limit)/int_Flatte_f0; + + const EvtComplex term41 = sqrt(p30Jpsi_mid_CMS*p30Kp_mid_CMS); + + const EvtComplex term42 = X_J(1, p30Jpsi_mid_CMS, 1)*p30Jpsi_mid_CMS; + + const EvtComplex term43 = 1.2*f_S_NR/int_const_NR; + + const EvtComplex hm = term11*term12*term13 + term21*term22*term23 + + term31*term32*term33 + term41*term42*term43; + + // Increase by 10% + setProbMax(0.5*abs2(hm)*1.1); + +} + +// Decay function +void EvtBsMuMuKK::decay(EvtParticle *p) +{ + + EvtId other_b; + double time(0.0); + EvtCPUtil::getInstance()->OtherB(p,time,other_b); + time = -log(EvtRandom::Flat())*ctau; // This overrules the ctau made in OtherB + + if (EvtCPUtil::getInstance()->isBsMixed(p)) { + p->getParent()->setLifetime(time*EvtConst::c/1e12); // units: mm + } else { + p->setLifetime(time*EvtConst::c/1e12); // units: mm + } + + double DGtime = 0.25*deltaGamma*time; + double DMtime = 0.5*deltaMs*time; + double mt = exp(-DGtime); + double pt = exp(+DGtime); + double cDMt = cos(DMtime); + double sDMt = sin(DMtime); + EvtComplex termplus = EvtComplex(cDMt, sDMt); + EvtComplex terminus = EvtComplex(cDMt, -sDMt); + + EvtComplex gplus = 0.5*(mt*termplus + pt*terminus); + EvtComplex gminus = 0.5*(mt*termplus - pt*terminus); + + EvtId BSB = EvtPDL::getId("anti-B_s0"); + + // Flavour: first assume B_s0, otherwise choose anti-B_s0 + int q(1); + if (other_b == BSB) {q = -1;} + p->setAttribute("q",q); + + // Amplitudes + EvtComplex a_S_NR = AmpTime(q, gplus, gminus, delta_S_NR, + lambda_S_NR_abs, A_S_NR, phis_S_NR, -1); + + EvtComplex a_f0 = AmpTime(q, gplus, gminus, delta_f0, + lambda_f0_abs, A_f0, phis_f0, -1); + + EvtComplex a0_phi = AmpTime(q, gplus, gminus, delta_phi_0, + lambda_phi_0_abs, A_phi_0, phis_phi_0, 1); + + EvtComplex aperp_phi = AmpTime(q, gplus, gminus, delta_phi_perp, + lambda_phi_perp_abs, A_phi_perp, phis_phi_perp, -1); + + EvtComplex apar_phi = AmpTime(q, gplus, gminus, delta_phi_par, + lambda_phi_par_abs, A_phi_par, phis_phi_par, 1); + + EvtComplex a0_f2p = AmpTime(q, gplus, gminus, delta_f2p_0, + lambda_f2p_0_abs, A_f2p_0, phis_f2p_0, -1); + + EvtComplex aperp_f2p = AmpTime(q, gplus, gminus, delta_f2p_perp, + lambda_f2p_perp_abs, A_f2p_perp, phis_f2p_perp, 1); + + EvtComplex apar_f2p = AmpTime(q, gplus, gminus, delta_f2p_par, + lambda_f2p_par_abs, A_f2p_par, phis_f2p_par, -1); + + // Generate 4-momenta + double mKK = EvtRandom::Flat(kin_lower_limit, kin_upper_limit); + double mass[10] = {MJpsi, mKK, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + double Kmass[10] = {MKp, MKm, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + double muMass[10] = {Mmu, Mmu, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + + EvtVector4R mypV[2], mypK[2], mypmu[2]; + EvtGenKine::PhaseSpace(2, mass, mypV, MBs); + EvtGenKine::PhaseSpace(2, Kmass, mypK, mKK); + EvtGenKine::PhaseSpace(2, muMass, mypmu, MJpsi); + + EvtVector4R p4mup = boostTo(mypmu[0], mypV[0]); + EvtVector4R p4mum = boostTo(mypmu[1], mypV[0]); + EvtVector4R p4Kp = boostTo(mypK[0], mypV[1]); + EvtVector4R p4Km = boostTo(mypK[1], mypV[1]); + + p->makeDaughters(getNDaug(), getDaugs()); + + EvtParticle *thisparticle; + EvtParticle *muplus, *muminus, *Kplus, *Kminus; + + // Check particle ID + for (int k = 0; k <= 3; k++) { + + thisparticle = p->getDaug(k); + EvtId pId = thisparticle->getId(); + + if (pId == EvtPDL::getId("mu+")) { + + muplus = thisparticle; + muplus->init(getDaug(k), p4mup); + + } else if (pId == EvtPDL::getId("mu-")) { + + muminus = thisparticle; + muminus->init(getDaug(k), p4mum); + + } else if (pId == EvtPDL::getId("K+")) { + + Kplus = thisparticle; + Kplus->init(getDaug(k), p4Kp); + + } else if (pId == EvtPDL::getId("K-")) { + + Kminus = thisparticle; + Kminus->init(getDaug(k), p4Km); + } + + } + + EvtVector4R p4KK = p4Kp + p4Km; + EvtVector4R p4mumu = p4mup + p4mum; + EvtVector4R p4Bs = p4mumu + p4KK; + + double p4KK_mass2 = p4KK.mass2(); + double p4KK_mass = p4KK.mass(); + double p4Bs_mass2 = p4Bs.mass2(); + double p4Bs_mass = p4Bs.mass(); + + // Kp momentum in the KK CMS + double p3Kp_KK_CMS = sqrt( (p4KK_mass2 - pow(MKp+MKm,2)) * + (p4KK_mass2 - pow(MKp-MKm,2)) )/(2.0*p4KK_mass); + + // J/psi momentum in the KK CMS + double p3Jpsi_KK_CMS = sqrt( (p4Bs_mass2 - pow(p4KK_mass+MJpsi,2)) * + (p4Bs_mass2 - pow(p4KK_mass-MJpsi,2)) )/(2.0*p4Bs_mass); + + // Mass lineshapes + + // Non-resonant S wave + EvtComplex P_NR = 1.0/int_const_NR; + + // f0 Flatte + EvtComplex F_f0 = Flatte(Mf0, p4KK_mass)/int_Flatte_f0; + + // phi Breit Wigner + EvtComplex BW_phi = Breit_Wigner(Gamma0phi, Mphi, p4KK_mass, 1, + p30Kp_phi_CMS, p3Kp_KK_CMS)/int_BW_phi; + + // f2' Breit Wigner + EvtComplex BW_f2p = Breit_Wigner(Gamma0f2p, Mf2p, p4KK_mass, 1, + p30Kp_f2p_CMS, p3Kp_KK_CMS )/int_BW_f2p; + + // Barrier factors: Always taking the lowest Bs L + double X_KK_0 = 1.0; + double X_KK_1 = X_J(1, p3Kp_KK_CMS, 0); + double X_KK_2 = X_J(2, p3Kp_KK_CMS, 0); + double X_NR_Jpsi_1 = X_J(1, p3Jpsi_KK_CMS, 1); + double X_f0_Jpsi_1 = X_J(1, p3Jpsi_KK_CMS, 1); + double X_phi_Jpsi_0 = 1.0; + double X_f2p_Jpsi_1 = X_J(1, p3Jpsi_KK_CMS, 1); + + // Birth momentum factors: pow(p3(K+),LR)* pow(p3(J/psi),LB) + double f_PHSP = sqrt(p3Jpsi_KK_CMS*p3Kp_KK_CMS); + double f_BMF_NR = p3Jpsi_KK_CMS; + double f_BMF_f0 = p3Jpsi_KK_CMS; + double f_BMF_phi = p3Kp_KK_CMS; + double f_BMF_f2p = p3Kp_KK_CMS*p3Kp_KK_CMS*p3Jpsi_KK_CMS; + + // Angular distribution and sum over KK states + double CosK = EvtDecayAngle(p4Bs, p4KK, p4Kp); + double CosMu = EvtDecayAngle(p4Bs, p4mumu, p4mup); + double chi = EvtDecayAngleChi(p4Bs, p4mup, p4mum, p4Kp, p4Km); + + // Build helicity amplitudes + + // phi + EvtComplex H0_phi = a0_phi; + EvtComplex Hp_phi = (apar_phi + aperp_phi)/sq2; + EvtComplex Hm_phi = (apar_phi - aperp_phi)/sq2; + + // f2p + EvtComplex H0_f2p = a0_f2p; + EvtComplex Hp_f2p = (apar_f2p + aperp_f2p)/sq2; + EvtComplex Hm_f2p = (apar_f2p - aperp_f2p)/sq2; + + // muon polarization +1 + EvtComplex swaveangdist1 = AngularDist(0, 0, 1, CosK,CosMu,chi); + + // KK Spin-0 NR + EvtComplex mp_hS_NR = a_S_NR*swaveangdist1; + EvtComplex Amp_p_NR = P_NR*X_KK_0*X_NR_Jpsi_1*f_BMF_NR * mp_hS_NR; + + // KK Spin-0 f0 + EvtComplex mp_h_f0 = a_f0*swaveangdist1; + EvtComplex Amp_p_f0 = F_f0*X_KK_0*X_f0_Jpsi_1*f_BMF_f0 * mp_h_f0; + + // KK Spin-1 + EvtComplex mp_h0_phi = H0_phi*AngularDist(1, 0, 1, CosK,CosMu,chi); + EvtComplex mp_hp_phi = Hp_phi*AngularDist(1, 1, 1, CosK,CosMu,chi); + EvtComplex mp_hm_phi = Hm_phi*AngularDist(1,-1, 1, CosK,CosMu,chi); + EvtComplex Amp_p_phi = BW_phi*X_KK_1*X_phi_Jpsi_0*f_BMF_phi * + (mp_h0_phi + mp_hp_phi + mp_hm_phi); + + // KK Spin-2 + EvtComplex mp_h0_f2p = H0_f2p*AngularDist(2, 0, 1, CosK,CosMu,chi); + EvtComplex mp_hp_f2p = Hp_f2p*AngularDist(2, 1, 1, CosK,CosMu,chi); + EvtComplex mp_hm_f2p = Hm_f2p*AngularDist(2,-1, 1, CosK,CosMu,chi); + EvtComplex Amp_p_f2p = BW_f2p*X_KK_2*X_f2p_Jpsi_1*f_BMF_f2p * + (mp_h0_f2p + mp_hp_f2p + mp_hm_f2p); + + // muon polarization -1 + EvtComplex swaveangdist2 = AngularDist(0, 0,-1, CosK,CosMu,chi); + + // KK Spin-0 NR + EvtComplex mm_hS_NR = a_S_NR*swaveangdist2; + EvtComplex Amp_m_NR = P_NR*X_KK_0*X_NR_Jpsi_1*f_BMF_NR * mm_hS_NR; + + // KK Spin-0 + EvtComplex mm_h_f0 = a_f0*swaveangdist2; + EvtComplex Amp_m_f0 = F_f0*X_KK_0*X_f0_Jpsi_1*f_BMF_f0 * mm_h_f0; + + // KK Spin-1 + EvtComplex mm_h0_phi = H0_phi*AngularDist(1, 0,-1, CosK,CosMu,chi); + EvtComplex mm_hp_phi = Hp_phi*AngularDist(1,+1,-1, CosK,CosMu,chi); + EvtComplex mm_hm_phi = Hm_phi*AngularDist(1,-1,-1, CosK,CosMu,chi); + EvtComplex Amp_m_phi = BW_phi*X_KK_1*X_phi_Jpsi_0*f_BMF_phi * + (mm_h0_phi + mm_hp_phi + mm_hm_phi); + + // KK Spin-2 + EvtComplex mm_h0_f2p = H0_f2p*AngularDist(2, 0,-1, CosK,CosMu,chi); + EvtComplex mm_hp_f2p = Hp_f2p*AngularDist(2, 1,-1, CosK,CosMu,chi); + EvtComplex mm_hm_f2p = Hm_f2p*AngularDist(2,-1,-1, CosK,CosMu,chi); + EvtComplex Amp_m_f2p = BW_f2p*X_KK_2*X_f2p_Jpsi_1*f_BMF_f2p * + (mm_h0_f2p + mm_hp_f2p + mm_hm_f2p); + + // Total amplitudes + EvtComplex Amp_tot_plus = f_PHSP*(Amp_p_NR + Amp_p_f0 + Amp_p_phi + Amp_p_f2p); + EvtComplex Amp_tot_minus = f_PHSP*(Amp_m_NR + Amp_m_f0 + Amp_m_phi + Amp_m_f2p); + + vertex(0, 0, 0.0); + vertex(0, 1, Amp_tot_plus); + vertex(1, 0, Amp_tot_minus); + vertex(1, 1, 0.0); + +} + +// Rho function +EvtComplex EvtBsMuMuKK::GetRho(const double m0, const double m) const +{ + + double rho_sq = 1.0 - (4.0*m0*m0/(m*m)); + EvtComplex rho; + + if (rho_sq > 0.0) { + rho = EvtComplex(sqrt(rho_sq), 0.0); + } else { + rho = EvtComplex(0.0, sqrt(-rho_sq)); + } + + return rho; + +} + +// Flatte function +EvtComplex EvtBsMuMuKK::Flatte(const double m0, const double m) const +{ + double gpipi = 0.167; + double gKK = 3.05*gpipi; + + EvtComplex term1 = (2.0*GetRho(Mpip, m) + GetRho(Mpi0, m))/3.0; + EvtComplex term2 = (GetRho(MKp, m) + GetRho(MK0, m))/2.0; + + EvtComplex w = gpipi*term1 + gKK*term2; + + EvtComplex Flatte_0 = 1.0/(m0*m0 - m*m - I*m0*w); + + return Flatte_0; + +} + +// Breit-Wigner function +EvtComplex EvtBsMuMuKK::Breit_Wigner(const double Gamma0, const double m0, const double m, + const int J, const double q0, const double q) const +{ + double X_J_q0_sq = pow(X_J(J, q0, 0), 2); + double X_J_q_sq = pow(X_J(J, q, 0), 2); + + double Gamma = Gamma0*pow(q/q0,2*J+1)*(m0/m)*(X_J_q_sq/X_J_q0_sq); + + return 1.0/(m0*m0-m*m - I*m0*Gamma); +} + +// Integral +double EvtBsMuMuKK::Integral(const double Gamma0, const double m0, const int JR, const int JB, const double q0, + const double M_KK_ll, const double M_KK_ul, const int fcntype) const +{ + int bins = 1000; + double bin_width = (M_KK_ul-M_KK_ll)/static_cast(bins); + EvtComplex integral(0.0, 0.0); + double MKpiKm2 = pow(MKp+MKm,2); + double MBs2 = pow(MBs,2); + + for (int i = 0; i < bins; i++) { + + double M_KK_i = M_KK_ll + static_cast(i)*bin_width; + double M_KK_f = M_KK_ll + static_cast(i+1)*bin_width; + + double p3Kp_KK_CMS_i = (pow(M_KK_i,2) - MKpiKm2)/(2.0*M_KK_i); + double p3Kp_KK_CMS_f = (pow(M_KK_f,2) - MKpiKm2)/(2.0*M_KK_f); + + double p3Jpsi_Bs_CMS_i = sqrt((MBs2 - pow(M_KK_i+MJpsi,2)) * + (MBs2 - pow(M_KK_i-MJpsi,2)) )/(2.0*MBs); + double p3Jpsi_Bs_CMS_f = sqrt((MBs2 - pow(M_KK_f+MJpsi,2)) * + (MBs2 - pow(M_KK_f-MJpsi,2)) )/(2.0*MBs); + + double f_PHSP_i = sqrt(p3Kp_KK_CMS_i*p3Jpsi_Bs_CMS_i); + double f_PHSP_f = sqrt(p3Kp_KK_CMS_f*p3Jpsi_Bs_CMS_f); + + double f_MBF_KK_i = pow(p3Kp_KK_CMS_i, JR); + double f_MBF_KK_f = pow(p3Kp_KK_CMS_f, JR); + + double f_MBF_Bs_i = pow(p3Jpsi_Bs_CMS_i, JB); + double f_MBF_Bs_f = pow(p3Jpsi_Bs_CMS_f, JB); + + double X_JR_i = X_J(JR, p3Kp_KK_CMS_i, 0); + double X_JR_f = X_J(JR, p3Kp_KK_CMS_f, 0); + + double X_JB_i = X_J(JB, p3Jpsi_Bs_CMS_i, 1); + double X_JB_f = X_J(JB, p3Jpsi_Bs_CMS_f, 1); + + EvtComplex fcn_i(1.0, 0.0), fcn_f(1.0, 0.0); + + if (fcntype == 1) { + + fcn_i = Flatte(m0, M_KK_i); + fcn_f = Flatte(m0, M_KK_f); + + } else if (fcntype == 2) { + + fcn_i = Breit_Wigner(Gamma0, m0, M_KK_i, JR, q0, p3Kp_KK_CMS_i); + fcn_f = Breit_Wigner(Gamma0, m0, M_KK_f, JR, q0, p3Kp_KK_CMS_f); + + } + + EvtComplex a_i = f_PHSP_i*f_MBF_KK_i*f_MBF_Bs_i*X_JR_i*X_JB_i*fcn_i; + EvtComplex a_st_i = conj(a_i); + EvtComplex a_f = f_PHSP_f*f_MBF_KK_f*f_MBF_Bs_f*X_JR_f*X_JB_f*fcn_f; + EvtComplex a_st_f = conj(a_f); + + integral += 0.5*bin_width*(a_i*a_st_i + a_f*a_st_f); + + } + + return sqrt(abs2(integral)); + +} + +// Blatt-Weisskopf barrier factors +double EvtBsMuMuKK::X_J(const int J, const double q, const int isB) const +{ + + double r_BW = 1.0; + + if (isB == 0) { + r_BW = 1.5; + } else if (isB == 1) { + r_BW = 5.0; + } + + double zsq = pow(r_BW*q,2); + + double X_J(1.0); + + if (J == 1) { + X_J = sqrt(1.0/(1.0+zsq)); + } else if (J == 2) { + X_J = sqrt(1.0/(zsq*zsq+3.0*zsq+9.0)); + } + + return X_J; +} + + +// EvtGen d matrix: Input is 2J instead of J etc +double EvtBsMuMuKK::Wignerd(const int J, const int l, const int alpha, const double theta) const +{ + return EvtdFunction::d(2*J, 2*l, 2*alpha, theta); +} + + + +// J spin of KK, l spin proj of J/psi, alpha dimuon spin +EvtComplex EvtBsMuMuKK::AngularDist(const int J, const int l, const int alpha, + const double cK, const double cL, const double chi) const +{ + + double thetaL = acos(cL); + double thetaK = acos(cK); + + EvtComplex out = 0.5*sqrt((2*J+1)/pi)*exp(EvtComplex(0,-l*chi)); + + out *= Wignerd(1, l, alpha, thetaL)*Wignerd(J, -l, 0, thetaK); + + return out; + +} + +// Time-dependent amplitude calculation +EvtComplex EvtBsMuMuKK::AmpTime(const int q, const EvtComplex& gplus, const EvtComplex& gminus, + const double delta, const double lambda_abs, const double Amp, + const double phis, const int eta) const +{ + + EvtComplex amp_time = Amp * EvtComplex(cos(-delta), sin(-delta)); + double qphis = q*phis; + amp_time *= (gplus + eta*pow(lambda_abs,-1.0*q) * EvtComplex(cos(qphis), sin(qphis)) * gminus); + + if (q == 1) {amp_time *= eta;} + + return amp_time; + +} diff --git a/src/EvtGenModels/EvtModelReg.cpp b/src/EvtGenModels/EvtModelReg.cpp index b476060..ac72a59 100644 --- a/src/EvtGenModels/EvtModelReg.cpp +++ b/src/EvtGenModels/EvtModelReg.cpp @@ -1,353 +1,354 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 1998 Caltech, UCSB // // Module: EvtModel.cc // // Description: // // Modification history: // // RYD September 25, 1996 Module created // //------------------------------------------------------------------------ // #include "EvtGenBase/EvtPatches.hh" #include #include #include #include #include #include #include "EvtGenBase/EvtModel.hh" #include "EvtGenModels/EvtTauVectornu.hh" #include "EvtGenModels/EvtVVP.hh" #include "EvtGenModels/EvtSLN.hh" #include "EvtGenModels/EvtISGW2.hh" #include "EvtGenModels/EvtMelikhov.hh" #include "EvtGenModels/EvtSLPole.hh" #include "EvtGenModels/EvtPropSLPole.hh" #include "EvtGenModels/EvtSLBKPole.hh" #include "EvtGenModels/EvtISGW.hh" #include "EvtGenModels/EvtBHadronic.hh" #include "EvtGenModels/EvtVSS.hh" #include "EvtGenModels/EvtVSSMix.hh" #include "EvtGenModels/EvtVSSBMixCPT.hh" #include "EvtGenModels/EvtVSPPwave.hh" #include "EvtGenModels/EvtGoityRoberts.hh" #include "EvtGenModels/EvtSVS.hh" #include "EvtGenModels/EvtTSS.hh" #include "EvtGenModels/EvtTVSPwave.hh" #include "EvtGenModels/EvtSVVHelAmp.hh" #include "EvtGenModels/EvtSVPHelAmp.hh" #include "EvtGenModels/EvtSVPCP.hh" #include "EvtGenModels/EvtVVSPwave.hh" #include "EvtGenModels/EvtDDalitz.hh" #include "EvtGenModels/EvtOmegaDalitz.hh" #include "EvtGenModels/EvtPi0Dalitz.hh" #include "EvtGenModels/EvtEtaDalitz.hh" #include "EvtGenModels/EvtPhsp.hh" #include "EvtGenModels/EvtBtoXsgamma.hh" #include "EvtGenModels/EvtBtoXsll.hh" #include "EvtGenModels/EvtBtoXsEtap.hh" #include "EvtGenModels/EvtSSSCP.hh" #include "EvtGenModels/EvtSSSCPpng.hh" #include "EvtGenModels/EvtSTSCP.hh" #include "EvtGenModels/EvtSTS.hh" #include "EvtGenModels/EvtSSSCPT.hh" #include "EvtGenModels/EvtSVSCP.hh" #include "EvtGenModels/EvtSSDCP.hh" #include "EvtGenModels/EvtSVSNONCPEIGEN.hh" #include "EvtGenModels/EvtSVVNONCPEIGEN.hh" #include "EvtGenModels/EvtSVVCP.hh" #include "EvtGenModels/EvtSVVCPLH.hh" #include "EvtGenModels/EvtSVSCPLH.hh" #include "EvtGenModels/EvtSll.hh" #include "EvtGenModels/EvtVll.hh" #include "EvtGenModels/EvtTaulnunu.hh" #include "EvtGenModels/EvtTauHadnu.hh" #include "EvtGenModels/EvtTauScalarnu.hh" #include "EvtGenModels/EvtKstarnunu.hh" #include "EvtGenModels/EvtbTosllBall.hh" #include "EvtGenModels/EvtSingleParticle.hh" #include "EvtGenModels/EvtVectorIsr.hh" #include "EvtGenModels/EvtBToPlnuBK.hh" #include "EvtGenModels/EvtBToVlnuBall.hh" #include "EvtGenModels/EvtSVVHelCPMix.hh" #include "EvtGenModels/EvtSVPHelCPMix.hh" #include "EvtGenModels/EvtLb2Lll.hh" #include "EvtGenModels/EvtRareLbToLll.hh" #include "EvtGenModels/EvtHypNonLepton.hh" #include "EvtGenModels/EvtbTosllAli.hh" #include "EvtGenModels/EvtBToDDalitzCPK.hh" #include "EvtGenModels/EvtPVVCPLH.hh" #include "EvtGenModels/EvtLambdaB2LambdaV.hh" #include "EvtGenModels/EvtSSD_DirectCP.hh" #include "EvtGenModels/EvtHQET.hh" #include "EvtGenModels/EvtHQET2.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenModels/EvtHelAmp.hh" #include "EvtGenModels/EvtPartWave.hh" #include "EvtGenModels/EvtBto2piCPiso.hh" #include "EvtGenModels/EvtBtoKpiCPiso.hh" #include "EvtGenModels/EvtSVSCPiso.hh" #include "EvtGenModels/EvtVVpipi.hh" #include "EvtGenModels/EvtY3SToY1SpipiMoxhay.hh" #include "EvtGenModels/EvtYmSToYnSpipiCLEO.hh" #include "EvtGenModels/EvtVVPIPI_WEIGHTED.hh" #include "EvtGenModels/EvtVPHOtoVISRHi.hh" #include "EvtGenModels/EvtBTo4piCP.hh" #include "EvtGenModels/EvtBTo3piCP.hh" #include "EvtGenModels/EvtCBTo3piP00.hh" #include "EvtGenModels/EvtCBTo3piMPP.hh" #include "EvtGenModels/EvtBToKpipiCP.hh" #include "EvtGenModels/EvtBsquark.hh" #include "EvtGenModels/EvtPhiDalitz.hh" #include "EvtGenModels/EvtLNuGamma.hh" #include "EvtGenModels/EvtVub.hh" #include "EvtGenModels/EvtVubHybrid.hh" #include "EvtGenModels/EvtVubNLO.hh" #include "EvtGenModels/EvtVubBLNP.hh" #include "EvtGenModels/EvtVubBLNPHybrid.hh" #include "EvtGenModels/EvtPto3P.hh" #include "EvtGenModels/EvtBtoKD3P.hh" #include "EvtGenModels/EvtKstarstargamma.hh" #include "EvtGenModels/EvtFlatQ2.hh" #include "EvtGenModels/EvtLambdaP_BarGamma.hh" #include "EvtGenModels/EvtBBScalar.hh" #include "EvtGenModels/EvtKKLambdaC.hh" #include "EvtGenModels/EvtMultibody.hh" #include "EvtGenModels/EvtBaryonPCR.hh" #include "EvtGenModels/EvtDMix.hh" #include "EvtGenModels/EvtD0mixDalitz.hh" #include "EvtGenModels/EvtD0gammaDalitz.hh" #include "EvtGenModels/EvtEta2MuMuGamma.hh" #include "EvtGenModels/EvtBcToNPi.hh" #include "EvtGenModels/EvtBcPsiNPi.hh" #include "EvtGenModels/EvtBcBsNPi.hh" #include "EvtGenModels/EvtBcBsStarNPi.hh" #include "EvtGenModels/EvtBcSMuNu.hh" #include "EvtGenModels/EvtBcVMuNu.hh" #include "EvtGenModels/EvtBcTMuNu.hh" #include "EvtGenModels/EvtBcVNpi.hh" #include "EvtGenModels/EvtSVP.hh" #include "EvtGenModels/EvtTVP.hh" #include "EvtGenModels/EvtXPsiGamma.hh" #include "EvtGenModels/EvtbsToLLLL.hh" #include "EvtGenModels/EvtbsToLLLLHyperCP.hh" #include "EvtGenModels/EvtBLLNuL.hh" #include "EvtGenModels/EvtKStopizmumu.hh" #include "EvtGenModels/EvtVtoSll.hh" +#include "EvtGenModels/EvtBsMuMuKK.hh" #include "EvtGenModels/EvtGenericDalitz.hh" #include "EvtGenModels/EvtBcVHad.hh" #include "EvtGenModels/Evtbs2llGammaMNT.hh" #include "EvtGenModels/Evtbs2llGammaISRFSR.hh" #include "EvtGenModels/EvtbTosllMS.hh" #include "EvtGenModels/EvtbTosllMSExt.hh" #include "EvtGenModels/EvtbsToLLLL.hh" #include "EvtGenModels/EvtbsToLLLLHyperCP.hh" #include "EvtGenModels/EvtLb2plnuLCSR.hh" #include "EvtGenModels/EvtLb2plnuLQCD.hh" #include "EvtGenModels/EvtLb2Baryonlnu.hh" #include "EvtGenModels/EvtBToDiBaryonlnupQCD.hh" #include "EvtGenModels/EvtFlatSqDalitz.hh" #include "EvtGenModels/EvtPhspFlatLifetime.hh" #include "EvtGenModels/EvtModelReg.hh" using std::fstream; using std::cout; using std::endl; EvtModelReg::EvtModelReg(const std::list* extraModels) { EvtModel &modelist=EvtModel::instance(); if(extraModels){ for(std::list::const_iterator it = extraModels->begin(); it != extraModels->end(); ++it){ modelist.registerModel(*it); } } modelist.registerModel(new EvtBBScalar); modelist.registerModel(new EvtLambdaP_BarGamma); modelist.registerModel(new EvtFlatQ2); modelist.registerModel(new EvtTauHadnu); modelist.registerModel(new EvtTauVectornu); modelist.registerModel(new EvtVVP); modelist.registerModel(new EvtSLN); modelist.registerModel(new EvtISGW2); modelist.registerModel(new EvtMelikhov); modelist.registerModel(new EvtSLPole); modelist.registerModel(new EvtPropSLPole); modelist.registerModel(new EvtSLBKPole); modelist.registerModel(new EvtHQET); modelist.registerModel(new EvtHQET2); modelist.registerModel(new EvtISGW); modelist.registerModel(new EvtBHadronic); modelist.registerModel(new EvtVSS); modelist.registerModel(new EvtVSSMix); modelist.registerModel(new EvtVSSBMixCPT); modelist.registerModel(new EvtVSPPwave); modelist.registerModel(new EvtGoityRoberts); modelist.registerModel(new EvtSVS); modelist.registerModel(new EvtTSS); modelist.registerModel(new EvtTVSPwave); modelist.registerModel(new EvtSVVHelAmp); modelist.registerModel(new EvtSVPHelAmp); modelist.registerModel(new EvtSVPCP); modelist.registerModel(new EvtVVSPwave); modelist.registerModel(new EvtDDalitz); modelist.registerModel(new EvtOmegaDalitz); modelist.registerModel(new EvtEtaDalitz); modelist.registerModel(new EvtPhsp); modelist.registerModel(new EvtBtoXsgamma); modelist.registerModel(new EvtBtoXsll); modelist.registerModel(new EvtBtoXsEtap); modelist.registerModel(new EvtSSSCP); modelist.registerModel(new EvtSSSCPpng); modelist.registerModel(new EvtSTSCP); modelist.registerModel(new EvtSTS); modelist.registerModel(new EvtSSSCPT); modelist.registerModel(new EvtSVSCP); modelist.registerModel(new EvtSSDCP); modelist.registerModel(new EvtSVSNONCPEIGEN); modelist.registerModel(new EvtSVVNONCPEIGEN); modelist.registerModel(new EvtSVVCP); modelist.registerModel(new EvtSVVCPLH); modelist.registerModel(new EvtSVSCPLH); modelist.registerModel(new EvtSll); modelist.registerModel(new EvtVll); modelist.registerModel(new EvtTaulnunu); modelist.registerModel(new EvtTauScalarnu); modelist.registerModel(new EvtKstarnunu); modelist.registerModel(new EvtbTosllBall); modelist.registerModel(new EvtBto2piCPiso); modelist.registerModel(new EvtBtoKpiCPiso); modelist.registerModel(new EvtSVSCPiso); modelist.registerModel(new EvtSingleParticle); modelist.registerModel(new EvtVectorIsr); modelist.registerModel(new EvtPi0Dalitz); modelist.registerModel(new EvtHelAmp); modelist.registerModel(new EvtPartWave); modelist.registerModel(new EvtVVpipi); modelist.registerModel(new EvtY3SToY1SpipiMoxhay); modelist.registerModel(new EvtYmSToYnSpipiCLEO); modelist.registerModel(new EvtBsquark); modelist.registerModel(new EvtPhiDalitz); modelist.registerModel(new EvtBToPlnuBK); modelist.registerModel(new EvtBToVlnuBall); modelist.registerModel(new EvtVVPIPI_WEIGHTED); modelist.registerModel(new EvtVPHOtoVISRHi); modelist.registerModel(new EvtBTo4piCP); modelist.registerModel(new EvtBTo3piCP); modelist.registerModel(new EvtCBTo3piP00); modelist.registerModel(new EvtCBTo3piMPP); modelist.registerModel(new EvtBToKpipiCP); modelist.registerModel(new EvtLb2Lll); modelist.registerModel(new EvtRareLbToLll); modelist.registerModel(new EvtHypNonLepton); modelist.registerModel(new EvtSVVHelCPMix); modelist.registerModel(new EvtSVPHelCPMix); modelist.registerModel(new EvtLNuGamma); modelist.registerModel(new EvtKstarstargamma); modelist.registerModel(new EvtVub); modelist.registerModel(new EvtVubHybrid); modelist.registerModel(new EvtVubNLO); modelist.registerModel(new EvtVubBLNP); modelist.registerModel(new EvtVubBLNPHybrid); modelist.registerModel(new EvtPto3P); modelist.registerModel(new EvtBtoKD3P); modelist.registerModel(new EvtKKLambdaC); modelist.registerModel(new EvtMultibody); modelist.registerModel(new EvtDMix); modelist.registerModel(new EvtD0mixDalitz); modelist.registerModel(new EvtD0gammaDalitz); modelist.registerModel(new EvtEta2MuMuGamma); modelist.registerModel(new EvtbTosllAli); modelist.registerModel(new EvtBaryonPCR); modelist.registerModel(new EvtBToDDalitzCPK); modelist.registerModel(new EvtLambdaB2LambdaV); modelist.registerModel(new EvtLambda2PPiForLambdaB2LambdaV); modelist.registerModel(new EvtV2VpVmForLambdaB2LambdaV); modelist.registerModel(new EvtPVVCPLH); modelist.registerModel(new EvtSSD_DirectCP); modelist.registerModel(new EvtBcToNPi(true)); // true = print author info modelist.registerModel(new EvtBcPsiNPi); modelist.registerModel(new EvtBcBsNPi); modelist.registerModel(new EvtBcBsStarNPi); modelist.registerModel(new EvtBcSMuNu); modelist.registerModel(new EvtBcVMuNu); modelist.registerModel(new EvtBcTMuNu); modelist.registerModel(new EvtBcVNpi); modelist.registerModel(new EvtSVP); modelist.registerModel(new EvtTVP); modelist.registerModel(new EvtXPsiGamma); modelist.registerModel(new EvtbsToLLLL); modelist.registerModel(new EvtbsToLLLLHyperCP); modelist.registerModel(new EvtBLLNuL); modelist.registerModel(new EvtKStopizmumu); modelist.registerModel(new EvtVtoSll); modelist.registerModel(new EvtGenericDalitz()); modelist.registerModel(new EvtBcVHad); modelist.registerModel(new Evtbs2llGammaMNT); modelist.registerModel(new Evtbs2llGammaISRFSR); modelist.registerModel(new EvtbTosllMS); modelist.registerModel(new EvtbTosllMSExt); modelist.registerModel(new EvtLb2plnuLQCD); modelist.registerModel(new EvtLb2plnuLCSR); modelist.registerModel(new EvtLb2Baryonlnu); modelist.registerModel(new EvtBToDiBaryonlnupQCD); modelist.registerModel(new EvtFlatSqDalitz); modelist.registerModel(new EvtPhspFlatLifetime); } diff --git a/src/EvtGenModels/EvtRareLbToLll.cpp b/src/EvtGenModels/EvtRareLbToLll.cpp index 2be3e01..c810719 100644 --- a/src/EvtGenModels/EvtRareLbToLll.cpp +++ b/src/EvtGenModels/EvtRareLbToLll.cpp @@ -1,498 +1,504 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 2003 Caltech, UCSB // // Module: EvtRareLbToLll // -// Description: +// Description: // Implements the rare Lb --> Lambda^(*) ell ell models described in // http://arxiv.org/pdf/1108.6129.pdf // // Modification history: // // T. Blake November 2013 Created module // //------------------------------------------------------------------------ -// +// + - #include #include "EvtGenModels/EvtRareLbToLll.hh" #include "EvtGenModels/EvtRareLbToLllFF.hh" #include "EvtGenModels/EvtRareLbToLllFFGutsche.hh" #include "EvtGenModels/EvtRareLbToLllFFlQCD.hh" #include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtVector4C.hh" #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtDiracParticle.hh" #include "EvtGenBase/EvtRaritaSchwinger.hh" #include "EvtGenBase/EvtSpinDensity.hh" #include "EvtGenBase/EvtPDL.hh" -EvtRareLbToLll::EvtRareLbToLll() : m_maxProbability( 0 ), ffmodel_( 0 ), wcmodel_( 0 ) {} +EvtRareLbToLll::EvtRareLbToLll() : m_maxProbability( 0 ), ffmodel_( 0 ), wcmodel_( 0 ) {} EvtRareLbToLll::~EvtRareLbToLll() { if ( wcmodel_ ) delete wcmodel_; if ( ffmodel_ ) delete ffmodel_; } // The module name specification std::string EvtRareLbToLll::getName( ) { return "RareLbToLll" ; } // The implementation of the clone() method EvtDecayBase* EvtRareLbToLll::clone(){ return new EvtRareLbToLll; } void EvtRareLbToLll::init(){ checkNArg(1); - + // check that there are 3 daughteres checkNDaug(3); // Parent should be spin 1/2 Lambda_b0 const EvtSpinType::spintype spin = EvtPDL::getSpinType(getDaug(0)); - + if ( !( spin == EvtSpinType::DIRAC || spin == EvtSpinType::RARITASCHWINGER ) ) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << " EvtRareLbToLll expects DIRAC or RARITASWINGER daughter " << std::endl; } - - // We expect that the second and third daughters - // are the ell+ and ell- + + // We expect that the second and third daughters + // are the ell+ and ell- checkSpinDaughter(1,EvtSpinType::DIRAC); checkSpinDaughter(2,EvtSpinType::DIRAC); std::string model = getArgStr(0); if ( model == "Gutsche" ) { ffmodel_ = new EvtRareLbToLllFFGutsche(); } else if ( model == "LQCD" ) { ffmodel_ = new EvtRareLbToLllFFlQCD(); } else if ( model == "MR" ) { ffmodel_ = new EvtRareLbToLllFF(); } else { - EvtGenReport(EVTGEN_ERROR ,"EvtGen") << " Unknown form-factor model, valid options are MR, LQCD, Gutsche." << std::endl; + EvtGenReport(EVTGEN_INFO,"EvtGen") << " Unknown form-factor model, valid options are MR, LQCD, Gutsche." << std::endl; ::abort(); } wcmodel_ = new EvtRareLbToLllWC(); ffmodel_->init(); - + return; } void EvtRareLbToLll::initProbMax(){ EvtGenReport(EVTGEN_INFO,"EvtGen") << " EvtRareLbToLll is finding maximum probability ... " << std::endl; m_maxProbability=0; if(m_maxProbability==0){ EvtDiracParticle *parent = new EvtDiracParticle; parent->noLifeTime(); parent->init(getParentId(),EvtVector4R(EvtPDL::getMass(getParentId()),0,0,0)); parent->setDiagonalSpinDensity(); EvtAmp amp; EvtId daughters[3] = {getDaug(0),getDaug(1),getDaug(2)}; amp.init(getParentId(),3,daughters); parent->makeDaughters(3,daughters); EvtParticle *lambda = parent->getDaug(0); EvtParticle *lep1 = parent->getDaug(1); EvtParticle *lep2 = parent->getDaug(2); lambda -> noLifeTime(); lep1 -> noLifeTime(); lep2 -> noLifeTime(); EvtSpinDensity rho; rho.setDiag(parent->getSpinStates()); double M0 = EvtPDL::getMass(getParentId()); double mL = EvtPDL::getMass(getDaug(0)); double m1 = EvtPDL::getMass(getDaug(1)); double m2 = EvtPDL::getMass(getDaug(2)); double q2,pstar,elambda,theta; double q2min = (m1+m2)*(m1+m2); double q2max = (M0-mL)*(M0-mL); EvtVector4R p4lambda,p4lep1,p4lep2,boost; EvtGenReport(EVTGEN_INFO,"EvtGen") << " EvtRareLbToLll is probing whole phase space ..." << std::endl; int i,j; double prob=0; for(i=0;i<=100;i++){ q2 = q2min+i*(q2max-q2min)/100.; elambda = (M0*M0+mL*mL-q2)/2/M0; if(i==0) pstar = 0; else pstar = sqrt(q2-(m1+m2)*(m1+m2))*sqrt(q2-(m1-m2)*(m1-m2))/2/sqrt(q2); boost.set(M0-elambda,0,0,+sqrt(elambda*elambda-mL*mL)); if ( i != 100 ) { p4lambda.set(elambda,0,0,-sqrt(elambda*elambda-mL*mL)); } else { - p4lambda.set(mL,0,0,0); + p4lambda.set(mL,0,0,0); } for(j=0;j<=45;j++){ theta = j*EvtConst::pi/45; p4lep1.set(sqrt(pstar*pstar+m1*m1),0,+pstar*sin(theta),+pstar*cos(theta)); p4lep2.set(sqrt(pstar*pstar+m2*m2),0,-pstar*sin(theta),-pstar*cos(theta)); //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " pstar: " << pstar << std::endl; if ( i != 100 ) // At maximal q2 we are already in correct frame as Lambda and W/Zvirtual are at rest { p4lep1 = boostTo(p4lep1,boost); p4lep2 = boostTo(p4lep2,boost); } lambda -> init(getDaug(0),p4lambda); lep1 -> init(getDaug(1),p4lep1 ); lep2 -> init(getDaug(2),p4lep2 ); calcAmp(amp,parent); prob = rho.normalizedProb(amp.getSpinDensity()); //std::cout << "q2: " << q2 << " \t theta: " << theta << " \t prob: " << prob << std::endl; //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " q2-q2min: " << q2-(m1+m2)*(m1+m2) << std::endl; if(prob>m_maxProbability){ - EvtGenReport(EVTGEN_INFO,"EvtGen") << " - probability " << prob - << " found at q2 = " << q2 << " (" << 100*(q2-q2min)/(q2max-q2min) - << " %) and theta = " << theta*180/EvtConst::pi << std::endl; + EvtGenReport(EVTGEN_INFO,"EvtGen") << " - probability " << prob + << " found at q2 = " << q2 << " (" << 100*(q2-q2min)/(q2max-q2min) + << " %) and theta = " << theta*180/EvtConst::pi << std::endl; m_maxProbability=prob; } } //::abort(); } //m_poleSize = 0.04*q2min; m_maxProbability *= 1.2; delete parent; } setProbMax(m_maxProbability); EvtGenReport(EVTGEN_INFO,"EvtGen") << " EvtRareLbToLll set up maximum probability to " << m_maxProbability << std::endl; } void EvtRareLbToLll::decay( EvtParticle *parent ){ parent->initializePhaseSpace( getNDaug(),getDaugs() ); calcAmp( _amp2,parent ); } -bool EvtRareLbToLll::isParticle( EvtParticle *parent ) const +bool EvtRareLbToLll::isParticle( EvtParticle *parent ) const { static EvtIdSet partlist("Lambda_b0"); - + return partlist.contains( parent->getId() ); } void EvtRareLbToLll::calcAmp( EvtAmp& amp, EvtParticle *parent ){ //parent->setDiagonalSpinDensity(); EvtParticle* lambda = parent->getDaug(0); static EvtIdSet leptons("e-","mu-","tau-"); - + const bool isparticle = isParticle( parent ); - + EvtParticle* lp = 0; EvtParticle* lm = 0; - + if ( leptons.contains(parent->getDaug(1)->getId()) ) { lp = parent->getDaug(1); lm = parent->getDaug(2); } - else + else { lp = parent->getDaug(2); lm = parent->getDaug(1); - } - + } + EvtVector4R P; - P.set(parent->mass(), 0.0,0.0,0.0); - + P.set(parent->mass(), 0.0,0.0,0.0); + EvtVector4R q = lp->getP4() + lm->getP4(); double qsq = q.mass2(); - + // Leptonic currents - EvtVector4C LV[2][2]; // \bar{\ell} \gamma^{\mu} \ell + EvtVector4C LV[2][2]; // \bar{\ell} \gamma^{\mu} \ell EvtVector4C LA[2][2]; // \bar{\ell} \gamma^{\mu} \gamma^{5} \ell for ( int i =0 ; i < 2; ++i ){ for ( int j = 0; j < 2; ++j ){ if ( isparticle ) { LV[i][j] = EvtLeptonVCurrent( lp->spParent(i), lm->spParent(j) ); LA[i][j] = EvtLeptonACurrent( lp->spParent(i), lm->spParent(j) ); } - else + else { LV[i][j] = EvtLeptonVCurrent( lp->spParent(1-i), lm->spParent(1-j) ); LA[i][j] = EvtLeptonACurrent( lp->spParent(1-i), lm->spParent(1-j) ); } } } - + EvtRareLbToLllFF::FormFactors FF; //F, G, FT and GT ffmodel_->getFF( parent, lambda, FF ); - + EvtComplex C7eff = wcmodel_-> GetC7Eff(qsq); EvtComplex C9eff = wcmodel_-> GetC9Eff(qsq); EvtComplex C10eff = wcmodel_->GetC10Eff(qsq); - + EvtComplex AC[4]; EvtComplex BC[4]; EvtComplex DC[4]; EvtComplex EC[4]; // check to see if particle is same or opposite parity to Lb const int parity = ffmodel_->isNatural( lambda ) ? 1 : -1 ; // Lambda spin type const EvtSpinType::spintype spin = EvtPDL::getSpinType(lambda->getId()); - + static const double mb = 5.209; - + // Eq. 48 + 49 for ( unsigned int i = 0; i < 4; ++i ) { if ( parity > 0 ) { - AC[i] = -2.*mb*C7eff*FF.FT_[i]/qsq + parity*C9eff*FF.F_[i]; - BC[i] = -2.*mb*C7eff*FF.GT_[i]/qsq - parity*C9eff*FF.G_[i]; + AC[i] = -2.*mb*C7eff*FF.FT_[i]/qsq + C9eff*FF.F_[i]; + BC[i] = -2.*mb*C7eff*FF.GT_[i]/qsq - C9eff*FF.G_[i]; DC[i] = C10eff*FF.F_[i]; EC[i] = -C10eff*FF.G_[i]; } - else + else { AC[i] = -2.*mb*C7eff*FF.GT_[i]/qsq - C9eff*FF.G_[i]; BC[i] = -2.*mb*C7eff*FF.FT_[i]/qsq + C9eff*FF.F_[i]; DC[i] = -C10eff*FF.G_[i]; EC[i] = C10eff*FF.F_[i]; - } + } } - + // handle particle -> antiparticle in Hadronic currents const double cv = ( isparticle > 0 ) ? 1.0 : -1.0*parity; - const double ca = ( isparticle > 0 ) ? 1.0 : 1.0*parity; - const double cs = ( isparticle > 0 ) ? 1.0 : 1.0*parity; + const double ca = ( isparticle > 0 ) ? 1.0 : +1.0*parity; + const double cs = ( isparticle > 0 ) ? 1.0 : +1.0*parity; const double cp = ( isparticle > 0 ) ? 1.0 : -1.0*parity; - - + + if (EvtSpinType::DIRAC == spin ){ EvtVector4C H1[2][2]; // vector current EvtVector4C H2[2][2]; // axial-vector - + EvtVector4C T[6]; // Hadronic currents for ( int i =0 ; i < 2; ++i ){ - for ( int j = 0; j < 2; ++j ){ - HadronicAmp( parent, lambda, T, i, j ); - - H1[i][j] = ( cv*AC[0]*T[0] + ca*BC[0]*T[1] + cs*AC[1]*T[2] + - cp*BC[1]*T[3] + cs*AC[2]*T[4] + cp*BC[2]*T[5] ); - - H2[i][j] = ( cv*DC[0]*T[0] + ca*EC[0]*T[1] + cs*DC[1]*T[2] + - cp*EC[1]*T[3] + cs*DC[2]*T[4] + cp*EC[2]*T[5] ); - } + for ( int j = 0; j < 2; ++j ){ + + HadronicAmp( parent, lambda, T, i, j ); + + + H1[i][j] = ( cv*AC[0]*T[0] + ca*BC[0]*T[1] + cs*AC[1]*T[2] + + cp*BC[1]*T[3] + cs*AC[2]*T[4] + cp*BC[2]*T[5] ); + + H2[i][j] = ( cv*DC[0]*T[0] + ca*EC[0]*T[1] + cs*DC[1]*T[2] + + cp*EC[1]*T[3] + cs*DC[2]*T[4] + cp*EC[2]*T[5] ); + } } // Spin sums int spins[4]; for ( int i =0; i < 2; ++i ) { for ( int ip = 0; ip < 2; ++ip ) { for ( int j = 0; j < 2; ++j ) { for ( int jp = 0; jp < 2; ++jp ) { spins[0] = i; spins[1] = ip; spins[2] = j; spins[3] = jp; - - amp.vertex( spins, H1[i][ip]*LV[j][jp] + H2[i][ip]*LA[j][jp] ); + + EvtComplex M = H1[i][ip]*LV[j][jp] + H2[i][ip]*LA[j][jp]; + + amp.vertex( spins, M ); } - } + } } } } else if ( EvtSpinType::RARITASCHWINGER == spin ) { EvtVector4C T[8]; - - EvtVector4C H1[4][2]; // vector current - EvtVector4C H2[4][2]; // axial-vector - + + EvtVector4C H1[2][4]; // vector current // swaped + EvtVector4C H2[2][4]; // axial-vector + // Build hadronic amplitude - for ( int i =0 ; i < 4; ++i ){ - for ( int j = 0; j < 2; ++j ){ - HadronicAmpRS( parent, lambda, T, i, j ); - - H1[i][j] = ( cv*AC[0]*T[0] + ca*BC[0]*T[1] + + for ( int i =0 ; i < 2; ++i ){ + for ( int j = 0; j < 4; ++j ){ + + H1[i][j] = ( cv*AC[0]*T[0] + ca*BC[0]*T[1] + cs*AC[1]*T[2] + cp*BC[1]*T[3] + - cs*AC[2]*T[4] + cp*BC[2]*T[5] + + cs*AC[2]*T[4] + cp*BC[2]*T[5] + cs*AC[3]*T[6] + cp*BC[3]*T[7] ); - H2[i][j] = ( cv*DC[0]*T[0] + ca*EC[0]*T[1] + + H2[i][j] = ( cv*DC[0]*T[0] + ca*EC[0]*T[1] + cs*DC[1]*T[2] + cp*EC[1]*T[3] + - cs*DC[2]*T[4] + cp*EC[2]*T[5] + + cs*DC[2]*T[4] + cp*EC[2]*T[5] + cs*DC[3]*T[6] + cp*EC[3]*T[7] ); } } - + // Spin sums int spins[4]; - - for ( int i =0; i < 4; ++i ) { - for ( int ip = 0; ip < 2; ++ip ) { + + for ( int i =0; i < 2; ++i ) { + for ( int ip = 0; ip < 4; ++ip ) { for ( int j = 0; j < 2; ++j ) { for ( int jp = 0; jp < 2; ++jp ) { spins[0] = i; spins[1] = ip; spins[2] = j; spins[3] = jp; - - amp.vertex( spins, H1[i][ip]*LV[j][jp] + H2[i][ip]*LA[j][jp] ); + + EvtComplex M = H1[i][ip]*LV[j][jp] + H2[i][ip]*LA[j][jp]; + + amp.vertex( spins, M ); } - } + } } } } - else + else { EvtGenReport(EVTGEN_ERROR,"EvtGen" ) << " EvtRareLbToLll expects DIRAC or RARITASWINGER daughter " << std::endl; } - + return; } // spin 1/2 daughters -void EvtRareLbToLll::HadronicAmp( EvtParticle* parent, - EvtParticle* lambda, +void EvtRareLbToLll::HadronicAmp( EvtParticle* parent, + EvtParticle* lambda, EvtVector4C* T, - const int i, + const int i, const int j ) { - - const EvtDiracSpinor Sfinal = lambda->spParent(i); - const EvtDiracSpinor Sinit = parent->sp(j); - + + const EvtDiracSpinor Sfinal = lambda->spParent(j); + const EvtDiracSpinor Sinit = parent->sp(i); + const EvtVector4R L = lambda->getP4(); - + EvtVector4R P; P.set(parent->mass(), 0.0,0.0,0.0); - + const double Pm = parent->mass(); const double Lm = lambda->mass(); - + // \bar{u} \gamma^{\mu} u T[0] = EvtLeptonVCurrent( Sfinal, Sinit ); // \bar{u} \gamma^{\mu}\gamma^{5} u T[1] = EvtLeptonACurrent( Sfinal, Sinit ); - + // \bar{u} v^{\mu} u T[2] = EvtLeptonSCurrent( Sfinal, Sinit )*( P/Pm ); - + // \bar{u} v^{\mu} \gamma^{5} u T[3] = EvtLeptonPCurrent( Sfinal, Sinit )*( P/Pm ); - + // \bar{u} v^{\prime\mu} u T[4] = EvtLeptonSCurrent( Sfinal, Sinit )*( L/Lm ); - - // \bar{u} v^{\prime\mu} \gamma^{5} + + // \bar{u} v^{\prime\mu} \gamma^{5} T[5] = EvtLeptonPCurrent( Sfinal, Sinit )*( L/Lm ); - - // v = p_{\Lambda_b}/m_{\Lambda_b} - // v^{\prime} = p_{\Lambda}/m_{\Lambda} - + + // Where: + // v = p_{\Lambda_b}/m_{\Lambda_b} + // v^{\prime} = p_{\Lambda}/m_{\Lambda} + return; } // spin 3/2 daughters -void EvtRareLbToLll::HadronicAmpRS( EvtParticle* parent, - EvtParticle* lambda, +void EvtRareLbToLll::HadronicAmpRS( EvtParticle* parent, + EvtParticle* lambda, EvtVector4C* T, - const int i, + const int i, const int j ) { - const EvtRaritaSchwinger Sfinal = lambda->spRSParent(i); - const EvtDiracSpinor Sinit = parent->sp(j); - + const EvtRaritaSchwinger Sfinal = lambda->spRSParent(j); + const EvtDiracSpinor Sinit = parent->sp(i); + EvtVector4R P; P.set(parent->mass(), 0.0,0.0,0.0); - + const EvtVector4R L = lambda->getP4(); - + EvtTensor4C ID; ID.setdiag(1.0,1.0,1.0,1.0); EvtDiracSpinor Sprime; - + for(int i=0 ; i<4 ; i++ ){ Sprime.set_spinor(i,Sfinal.getVector(i)*P); } const double Pmsq = P.mass2(); const double Pm = parent->mass(); const double PmLm = Pm*lambda->mass(); - - + + EvtVector4C V1,V2; - + for(int i=0;i<4;i++){ V1.set(i,EvtLeptonSCurrent(Sfinal.getSpinor(i),Sinit)); V2.set(i,EvtLeptonPCurrent(Sfinal.getSpinor(i),Sinit)); } // \bar{u}_{alpha} v^{\alpha} \gamma^{\mu} u - T[0] = EvtLeptonVCurrent(Sprime, Sinit)*(1/Pm); + T[0] = EvtLeptonVCurrent(Sprime, Sinit)*(1/Pm); - // \bar{u}_{alpha} v^{\alpha} \gamma^{\mu} \gamma^{5} u + // \bar{u}_{alpha} v^{\alpha} \gamma^{\mu} \gamma^{5} u T[1] = EvtLeptonACurrent(Sprime, Sinit)*(1/Pm); // \bar{u}_{\alpha} v^{\alpha} v^{\mu} u T[2] = EvtLeptonSCurrent(Sprime, Sinit)*(P/Pmsq); // \bar{u}_{\alpha} v^{\alpha} v^{\mu} \gamma^{5} u T[3] = EvtLeptonPCurrent(Sprime, Sinit)*(P/Pmsq); - + // \bar{u}_{\alpha} v^{\alpha} v^{\prime \mu} u T[4] = EvtLeptonSCurrent(Sprime, Sinit)*(L/PmLm); - + // \bar{u}_{\alpha} v^{\alpha} v^{\prime \mu} \gamma^{5} u T[5] = EvtLeptonPCurrent(Sprime, Sinit)*(L/PmLm); - + // \bar{u}_{\alpha} g^{\alpha\mu} u T[6] = ID.cont2(V1); // \bar{u}_{\alpha} g^{\alpha\mu} \gamma^{5} u T[7] = ID.cont2(V2); - + // Where: - // v = p_{\Lambda_b}/m_{\Lambda_b} - // v^{\prime} = p_{\Lambda}/m_{\Lambda} - + // v = p_{\Lambda_b}/m_{\Lambda_b} + // v^{\prime} = p_{\Lambda}/m_{\Lambda} + return; } diff --git a/src/EvtGenModels/EvtSLDiBaryonAmp.cpp b/src/EvtGenModels/EvtSLDiBaryonAmp.cpp index dc08ed0..bfc4a5f 100644 --- a/src/EvtGenModels/EvtSLDiBaryonAmp.cpp +++ b/src/EvtGenModels/EvtSLDiBaryonAmp.cpp @@ -1,325 +1,403 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 1998 Caltech, UCSB // // Module: EvtSLDiBaryonAmp.cc // // Description: Routine to implement semileptonic decays to dibaryonic final // state. Details of amplitude calculation to be found in arXiv:1107.0801. // // Modification history: // // Mark Smith July 18, 2017 Module created // Heavily adapted from the EvtSLBaryonAmp module // Ryan Newcombe May 2018 Added capability for Rarita-Schwinger // daughters. Indexing convention follows EvtSLBaryonAmp // John B Oct 2018 Optimise amplitude calculation // //-------------------------------------------------------------------------- #include "EvtGenBase/EvtPatches.hh" -#include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtGammaMatrix.hh" -#include "EvtGenBase/EvtId.hh" +#include "EvtGenBase/EvtIdSet.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtRaritaSchwinger.hh" #include "EvtGenBase/EvtReport.hh" -#include "EvtGenBase/EvtVector4C.hh" -#include "EvtGenBase/EvtVector4R.hh" +#include "EvtGenBase/EvtTensor4C.hh" -#include "EvtGenModels/EvtBToDiBaryonlnupQCDFF.hh" #include "EvtGenModels/EvtSLDiBaryonAmp.hh" using std::endl; -EvtSLDiBaryonAmp::EvtSLDiBaryonAmp(EvtBToDiBaryonlnupQCDFF* formFactors) : - ffModel(formFactors) +EvtSLDiBaryonAmp::EvtSLDiBaryonAmp(const EvtBToDiBaryonlnupQCDFF& formFactors) : + ffModel_(formFactors) { } void EvtSLDiBaryonAmp::CalcAmp(EvtParticle *parent, EvtAmp& amp) const { static EvtId EM = EvtPDL::getId("e-"); static EvtId MUM = EvtPDL::getId("mu-"); static EvtId TAUM = EvtPDL::getId("tau-"); static EvtId EP = EvtPDL::getId("e+"); static EvtId MUP = EvtPDL::getId("mu+"); static EvtId TAUP = EvtPDL::getId("tau+"); - // B charge (x3) to find parity factor and baryon daughter ordering + // The amplitude assumes B- -> p+ p- l- nubar ordering + // i.e. the B- decay is the "particle" mode + + // B charge (x3) to check for antiparticle mode and baryon daughter ordering EvtId BId = parent->getId(); int qB3 = EvtPDL::chg3(BId); - // The baryon, charged lepton and neutrino daughters. - // Make sure the first baryon has a charge opposite to the B + bool particleMode(true); + // Check if we have B+ instead (antiparticle mode) + if (qB3 > 0) {particleMode = false;} + + // The baryon, charged lepton and neutrino daughters + + // Make sure the first baryon has a charge opposite to the B, since the + // amplitude expressions assume this order EvtParticle* baryon1 = parent->getDaug(0); EvtParticle* baryon2 = parent->getDaug(1); + // Check if we need to reverse the baryon ordering if (EvtPDL::chg3(baryon1->getId()) == qB3) { - // Reverse baryon ordering baryon1 = parent->getDaug(1); baryon2 = parent->getDaug(0); } EvtParticle* lepton = parent->getDaug(2); EvtParticle* neutrino = parent->getDaug(3); - // 4-momenta + // 4-momenta in B rest frame EvtVector4R p0(parent->mass(), 0.0, 0.0, 0.0); EvtVector4R p1 = baryon1->getP4(); EvtVector4R p2 = baryon2->getP4(); EvtVector4R pSum = p1 + p2; EvtVector4R p = p0 - pSum; EvtVector4R pDiff = p2 - p1; // Particle id's: retrieve 1st baryon again in case order has changed EvtId Id1 = baryon1->getId(); EvtId Id2 = baryon2->getId(); EvtId l_num = lepton->getId(); - // Parity = 1 for B- mode; amplitude is written assuming we have B- - double parity(1.0); - // Parity = -1 factor for the vector terms for the B+ decay - if (qB3 > 0) {parity = -1.0;} - EvtSpinType::spintype type1 = EvtPDL::getSpinType(Id1); EvtSpinType::spintype type2 = EvtPDL::getSpinType(Id2); + + // Parity of B+- = -1. Check if the parity of the dibaryon state is the same. + // If so, set the sameParity integer to 1. Otherwise set it to -1, + // i.e. the dibaryon system has opposite parity to the B meson + int J1 = EvtSpinType::getSpin2(type1); + int J2 = EvtSpinType::getSpin2(type2); + int sameParity = this->checkDibaryonParity(Id1, Id2, J1, J2); // Number of chiral components of the baryon spinors int N1 = EvtSpinType::getSpinStates(type1); int N2 = EvtSpinType::getSpinStates(type2); // Invariant mass of the two baryon particle system double m_dibaryon = sqrt(pSum.mass2()); // Complex number i EvtComplex I(0, 1); // Lepton currents, same for all baryon options EvtVector4C l1, l2; if (l_num == EM || l_num == MUM || l_num == TAUM) { // B- l1 = EvtLeptonVACurrent(lepton->spParent(0), neutrino->spParentNeutrino()); l2 = EvtLeptonVACurrent(lepton->spParent(1), neutrino->spParentNeutrino()); } else if (l_num == EP || l_num == MUP || l_num == TAUP) { // B+ l1 = EvtLeptonVACurrent(neutrino->spParentNeutrino(), lepton->spParent(0)); l2 = EvtLeptonVACurrent(neutrino->spParentNeutrino(), lepton->spParent(1)); } else { EvtGenReport(EVTGEN_ERROR,"EvtSLDiBaryonAmp") << "Wrong lepton number"< sigmaVect; - for (int k = 0; k < 4; k++) { - - EvtGammaMatrix sigmaSum; - for (int index = 0; index < 4; index++) { - sigmaSum += EvtGammaMatrix::sigmaLower(k, index)*p.get(index); - } - - sigmaVect.push_back(sigmaSum); - - } - // Vector and axial-vector terms; these get reset within the loop over k - EvtVector4C term1; - EvtVector4C term2; + // Parity multiplication factors for the antiparticle mode hadronic currents + double sign1 = (particleMode == true) ? 1.0: 1.0*sameParity; + double sign2 = (particleMode == true) ? 1.0: 1.0*sameParity; + double sign3 = (particleMode == true) ? 1.0: -1.0*sameParity; + double sign4 = (particleMode == true) ? 1.0: -1.0*sameParity; + double sign5 = (particleMode == true) ? 1.0: -1.0*sameParity; + double sign6 = (particleMode == true) ? 1.0: 1.0*sameParity; - // Handle case of two Dirac-type daughters, e.g. ppbar, pN(1440) + // Define form factor coeff variables + double f1(0.0), f2(0.0), f3(0.0), f4(0.0), f5(0.0); + double g1(0.0), g2(0.0), g3(0.0), g4(0.0), g5(0.0); + + // Handle case of two Dirac-type daughters, e.g. p pbar, p N(1440) if (type1 == EvtSpinType::DIRAC && type2 == EvtSpinType::DIRAC) { // Form factor parameters EvtBToDiBaryonlnupQCDFF::FormFactors FF; - ffModel->getDiracFF(parent, m_dibaryon, FF); + ffModel_.getDiracFF(parent, m_dibaryon, FF); + + if (sameParity == 1) { + f1 = FF.F1; f2 = FF.F2; f3 = FF.F3; f4 = FF.F4; f5 = FF.F5; + g1 = FF.G1; g2 = FF.G2; g3 = FF.G3; g4 = FF.G4; g5 = FF.G5; + } else { + // Swap coeffs: f_i <--> g_i + f1 = FF.G1; f2 = FF.G2; f3 = FF.G3; f4 = FF.G4; f5 = FF.G5; + g1 = FF.F1; g2 = FF.F2; g3 = FF.F3; g4 = FF.F4; g5 = FF.F5; + } + + EvtVector4R gMtmTerms = g3*p + g4*pSum + g5*pDiff; + EvtVector4R fMtmTerms = f3*p + f4*pSum + f5*pDiff; // First baryon for (int i = 0; i < N1; i++) { // Get the baryon spinor in the B rest frame. Also just use u and not i*u, // since the imaginary constant factor is not needed for the probability EvtDiracSpinor u = baryon1->spParent(i); // Second baryon for(int j = 0; j < N2; j++) { EvtDiracSpinor v = baryon2->spParent(j); - EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v; - - // Calculate and set the 4 complex amplitude components - for (int k = 0; k < 4; k++) { - EvtGammaMatrix sigmaSum = sigmaVect[k]; - - // Need two terms owing to the requirements of using the product operators - EvtDiracSpinor vgTermA = (FF.G1*EvtGammaMatrix::g(k) + I*FF.G2*sigmaSum)*g5v; - EvtDiracSpinor vgTermB = (FF.G3*p.get(k) + FF.G4*pSum.get(k) + FF.G5*pDiff.get(k))*g5v; + // Hadronic currents + std::vector hadCurrents = this->getHadronicCurrents(u, v, p, gMtmTerms, + fMtmTerms); - // S current = ubar adjoint*vTerms - term1.set(k, EvtLeptonSCurrent(u, vgTermA+vgTermB)); + // First amplitude terms: 3rd current already has the form factor coeffs applied (gMtmTerms) + EvtVector4C amp1 = g1*sign1*hadCurrents[0] + g2*sign2*hadCurrents[1] + sign3*hadCurrents[2]; - EvtDiracSpinor vfTermA = (FF.F1*EvtGammaMatrix::g(k) + I*FF.F2*sigmaSum)*v; - EvtDiracSpinor vfTermB = (FF.F3*p.get(k) + FF.F4*pSum.get(k) + FF.F5*pDiff.get(k))*v; - - // S current = ubar adjoint*vTerms - term2.set(k, EvtLeptonSCurrent(u, vfTermA+vfTermB)); + // Second amplitude terms: 6th current already has the form factor coeffs applied (fMtmTerms) + EvtVector4C amp2 = f1*sign4*hadCurrents[3] + f2*sign5*hadCurrents[4] + sign6*hadCurrents[5]; + EvtVector4C hadAmp; + if (sameParity == 1) { + hadAmp = amp1 - amp2; + } else { + hadAmp = amp2 - amp1; } - // Set the V-A decay amplitude element - EvtVector4C term = parity*term1 - term2; - amp.vertex(i,j,0,l1*term); - amp.vertex(i,j,1,l2*term); + amp.vertex(i, j, 0, l1*hadAmp); + amp.vertex(i, j, 1, l2*hadAmp); } // j } // i } else if ((type1 == EvtSpinType::DIRAC && type2 == EvtSpinType::RARITASCHWINGER) || (type1 == EvtSpinType::RARITASCHWINGER && type2 == EvtSpinType::DIRAC)) { // Handle the case of one Dirac-type daughter (not including the leptons), e.g. one proton, and one // Rarita-Schwinger-type (spin 3/2) daughter e.g. B -> p N(1520) l nu // Form factor parameters EvtBToDiBaryonlnupQCDFF::FormFactors FF; - ffModel->getRaritaFF(parent, m_dibaryon, FF); + ffModel_.getRaritaFF(parent, m_dibaryon, FF); + + if (sameParity == 1) { + f1 = FF.F1; f2 = FF.F2; f3 = FF.F3; f4 = FF.F4; f5 = FF.F5; + g1 = FF.G1; g2 = FF.G2; g3 = FF.G3; g4 = FF.G4; g5 = FF.G5; + } else { + // Swap coeffs: f_i <--> g_i + f1 = FF.G1; f2 = FF.G2; f3 = FF.G3; f4 = FF.G4; f5 = FF.G5; + g1 = FF.F1; g2 = FF.F2; g3 = FF.F3; g4 = FF.F4; g5 = FF.F5; + } + + EvtVector4R gMtmTerms = g3*p + g4*pSum + g5*pDiff; + EvtVector4R fMtmTerms = f3*p + f4*pSum + f5*pDiff; if (type1 == EvtSpinType::DIRAC) { // First baryon is Dirac for (int i = 0; i < N1; i++) { // Get the baryon spinor in the B rest frame. Also just use u and not i*u, // since the imaginary constant factor is not needed for the probability EvtDiracSpinor u = baryon1->spParent(i); // Second baryon is RS-type for (int j = 0; j < N2; j++) { EvtRaritaSchwinger vRS = baryon2->spRSParent(j); - - // Store products of g5 with the spinor components as well as - // EvtDiracSpinors to limit constructor calls inside k loop - std::vector g5vVect, vVect; - for (int index = 0; index < 4; index++) { - - EvtDiracSpinor v = vRS.getSpinor(index); - EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v; - - g5vVect.push_back(g5v); - vVect.push_back(v); - + + EvtDiracSpinor v; + for (int k = 0; k < 4; k++) { + v.set_spinor(k, vRS.getVector(k)*p0); } - // Calculate and set the 4 complex amplitude components - for (int k = 0; k < 4; k++) { - - EvtGammaMatrix sigmaSum = sigmaVect[k]; - EvtDiracSpinor g5v = g5vVect[k]; - EvtDiracSpinor v = vVect[k]; + // Hadronic currents + std::vector hadCurrents = this->getHadronicCurrents(u, v, p, gMtmTerms, + fMtmTerms); - // Need two terms owing to the requirements of using the product operators - EvtDiracSpinor vgTermA = (FF.G1*EvtGammaMatrix::g(k) + I*FF.G2*sigmaSum)*g5v; - EvtDiracSpinor vgTermB = (FF.G3*p.get(k) + FF.G4*pSum.get(k) + FF.G5*pDiff.get(k))*g5v; + // First amplitude terms: 3rd current already has the form factor coeffs applied (gMtmTerms) + EvtVector4C amp1 = g1*sign1*hadCurrents[0] + g2*sign2*hadCurrents[1] + sign3*hadCurrents[2]; - // S current = ubar adjoint*vTerms - term1.set(k, EvtLeptonSCurrent(u, vgTermA+vgTermB)); - - EvtDiracSpinor vfTermA = (FF.F1*EvtGammaMatrix::g(k) + I*FF.F2*sigmaSum)*v; - EvtDiracSpinor vfTermB = (FF.F3*p.get(k) + FF.F4*pSum.get(k) + FF.F5*pDiff.get(k))*v; + // Second amplitude terms: 6th current already has the form factor coeffs applied (fMtmTerms) + EvtVector4C amp2 = f1*sign4*hadCurrents[3] + f2*sign5*hadCurrents[4] + sign6*hadCurrents[5]; - // S current = ubar adjoint*vTerms - term2.set(k, EvtLeptonSCurrent(u, vfTermA+vfTermB)); - + EvtVector4C hadAmp; + if (sameParity == 1) { + hadAmp = amp1 - amp2; + } else { + hadAmp = amp2 - amp1; } - // Set the V-A decay amplitude element - EvtVector4C term = parity*term1 - term2; - amp.vertex(i,j,0,l1*term); - amp.vertex(i,j,1,l2*term); - + amp.vertex(i, j, 0, l1*hadAmp); + amp.vertex(i, j, 1, l2*hadAmp); + } // j } // i } else if (type2 == EvtSpinType::DIRAC) { // Same as before, but where the first daughter is RS-type, e.g. B -> N(1520) p l nu // First baryon is RS for (int i = 0; i < N1; i++) { // Get the baryon spinor in the B rest frame EvtRaritaSchwinger uRS = baryon1->spRSParent(i); - - // Store EvtDiracSpinors to reduce constructor calls inside k loop - std::vector uVect; - for (int index = 0; index < 4; index++) { - - // Just use u and not i*u, since the imaginary constant factor is not needed - // for the probability - EvtDiracSpinor u = uRS.getSpinor(index); - uVect.push_back(u); - + + EvtDiracSpinor u; + for (int k = 0; k < 4; k++) { + u.set_spinor(k, uRS.getVector(k)*p0); } - + // Second baryon is Dirac for (int j = 0; j < N2; j++) { EvtDiracSpinor v = baryon2->spParent(j); - EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v; - - // Calculate and set the 4 complex amplitude components - for (int k = 0; k < 4; k++) { - - EvtGammaMatrix sigmaSum = sigmaVect[k]; - EvtDiracSpinor u = uVect[k]; - - // Need two terms owing to the requirements of using the product operators - EvtDiracSpinor vgTermA = (FF.G1*EvtGammaMatrix::g(k) + I*FF.G2*sigmaSum)*g5v; - EvtDiracSpinor vgTermB = (FF.G3*p.get(k) + FF.G4*pSum.get(k) + FF.G5*pDiff.get(k))*g5v; - - // S current = ubar adjoint*vTerms - term1.set(k, EvtLeptonSCurrent(u, vgTermA+vgTermB)); - - EvtDiracSpinor vfTermA = (FF.F1*EvtGammaMatrix::g(k) + I*FF.F2*sigmaSum)*v; - EvtDiracSpinor vfTermB = (FF.F3*p.get(k) + FF.F4*pSum.get(k) + FF.F5*pDiff.get(k))*v; - // S current = ubar adjoint*vTerms - term2.set(k, EvtLeptonSCurrent(u, vfTermA+vfTermB)); + // Hadronic currents + std::vector hadCurrents = this->getHadronicCurrents(u, v, p, gMtmTerms, + fMtmTerms); + + // First amplitude terms: 3rd current already has the form factor coeffs applied (gMtmTerms) + EvtVector4C amp1 = g1*sign1*hadCurrents[0] + g2*sign2*hadCurrents[1] + sign3*hadCurrents[2]; + // Second amplitude terms: 6th current already has the form factor coeffs applied (fMtmTerms) + EvtVector4C amp2 = f1*sign4*hadCurrents[3] + f2*sign5*hadCurrents[4] + sign6*hadCurrents[5]; + + EvtVector4C hadAmp; + if (sameParity == 1) { + hadAmp = amp1 - amp2; + } else { + hadAmp = amp2 - amp1; } - // Set the V-A decay amplitude element - EvtVector4C term = parity*term1 - term2; - amp.vertex(i,j,0,l1*term); - amp.vertex(i,j,1,l2*term); - + amp.vertex(i, j, 0, l1*hadAmp); + amp.vertex(i, j, 1, l2*hadAmp); + } // j } // i } // RS daughter check } // Have Dirac and RS baryons } + +std::vector EvtSLDiBaryonAmp::getHadronicCurrents(const EvtDiracSpinor& u, const EvtDiracSpinor& v, + const EvtVector4R& p, const EvtVector4R& gMtmTerms, + const EvtVector4R& fMtmTerms) const { + + // Store the currents used in Eq 6 (in order of appearance) + std::vector currents; + currents.reserve(6); + + EvtDiracSpinor g5v = EvtGammaMatrix::g5()*v; + + // ubar*gamma*gamma5*v + EvtVector4C current1 = EvtLeptonACurrent(u, v); + currents.push_back(current1); + + // ubar*sigma*p*gamma5*v -> [ubar*sigma*(gamma5*v)]*p + EvtTensor4C TC1 = EvtLeptonTCurrent(u, g5v); + // Contract tensor with 4-momentum + EvtVector4C current2 = TC1.cont2(p); + currents.push_back(current2); + + // ubar*p*gamma5*v; "p" = p, pSum and pDiff + EvtComplex PC1 = EvtLeptonPCurrent(u, v); + EvtVector4C current3 = PC1*gMtmTerms; + currents.push_back(current3); + + // ubar*gamma*v + EvtVector4C current4 = EvtLeptonVCurrent(u, v); + currents.push_back(current4); + + // ubar*sigma*p*v -> [ubar*sigma*v]*p + EvtTensor4C TC2 = EvtLeptonTCurrent(u, v); + // Contract tensor with 4-momentum + EvtVector4C current5 = TC2.cont2(p); + currents.push_back(current5); + + // ubar*p*v; "p" = p, pSum and pDiff + EvtComplex S1 = EvtLeptonSCurrent(u, v); + EvtVector4C current6 = S1*fMtmTerms; + currents.push_back(current6); + + return currents; + +} + +int EvtSLDiBaryonAmp::checkDibaryonParity(const EvtId& id1, const EvtId& id2, + const int J1, const int J2) const { + + // Get intrisic parities of the two baryons, then multiply by (-1)^|J1 - J2|. + // Note here that the J1 and J2 function arguments = 2*spin + int par1 = this->getBaryonParity(id1); + int par2 = this->getBaryonParity(id2); + + // mult should be either 0 or 1 for allowed Dirac/RS baryon pairs + int mult = static_cast(pow(-1.0, 0.5*fabs(J1 - J2))); + + int dbParity = par1*par2*mult; + + // Initialise result to 1, i.e. dibaryon parity = B parity = -1 + int result(1); + + // Dibaryon parity is opposite to the negative B parity + if (dbParity > 0) {result = -1;} + + return result; + +} + +int EvtSLDiBaryonAmp::getBaryonParity(const EvtId& id) const { + + // Initialise parity to +1 + int parity(1); + + // List of baryons with parity = +1 + static EvtIdSet posParity("p+", "Delta+", "Lambda_c+", "anti-Lambda_c(2593)-", + "anti-Lambda_c(2625)-", "N(1440)+", "anti-N(1520)-", + "anti-N(1535)-", "anti-N(1650)-", "anti-N(1700)-", + "N(1710)+", "N(1720)+"); + + // If the baryon id is not in the list, set the parity to -1 + if (!posParity.contains(id)) {parity = -1;} + + return parity; + +} diff --git a/src/EvtGenModels/EvtWHad.cpp b/src/EvtGenModels/EvtWHad.cpp index dbf8c33..3b1865c 100644 --- a/src/EvtGenModels/EvtWHad.cpp +++ b/src/EvtGenModels/EvtWHad.cpp @@ -1,190 +1,299 @@ //-------------------------------------------------------------------------- // // Environment: -// This software is part of the EvtGen package developed jointly -// for the BaBar and CLEO collaborations. If you use all or part +// This software is part of the EvtGen package. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT -// Copyright (C) 1998 Caltech, UCSB // // Module: EvtWHad.cc // // Description: Routine to calculate W -> (n pi) + (m K) current // according to [Kuhn, Was, Acta.Phys.Polon B39 (2008) 147] // // Modification history: -// AVL 20 Jan, 2013 Module created +// A V Luchinsky 20 Jan, 2013 Module created // //------------------------------------------------------------------------ -// #include "EvtGenModels/EvtWHad.hh" -using namespace std; +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtTensor4C.hh" + +EvtWHad::EvtWHad() : + mRho_(), + gamma0_(), + cK_(0) +{ + // cK coefficients from Eur. Phys. J. C39, 41 (2005), arXiv:hep-ph/0409080 [hep-ph] + + // rho(770) + mRho_.push_back(0.77511); + gamma0_.push_back(0.1491); + cK_.push_back(1.195); + + // rho(1450) + mRho_.push_back(1.465); + gamma0_.push_back(0.400); + cK_.push_back(-0.112); + + // rho(1700) + mRho_.push_back(1.72); + gamma0_.push_back(0.250); // rho(1700) + cK_.push_back(-0.083); + + // rho(2150), PRD 76 092005 + mRho_.push_back(2.150); + gamma0_.push_back(0.310); + cK_.push_back(0.0); +} + +EvtComplex EvtWHad::BWKK(double s, int i) const { + + double m2 = pow(mRho_[i], 2); + EvtComplex qs = pcm(s); + EvtComplex qm = pcm(m2); + if (abs(qm) < 1e-10) {return 0;} + + EvtComplex rat = qs/qm; + EvtComplex rat3 = rat*rat*rat; + if (abs(s) < 1e-10) {return 0;} + + EvtComplex gamma = m2*rat3*gamma0_[i]/s; + EvtComplex I(0.0, 1.0); + + EvtComplex denBW = m2 - s - I*sqrt(s)*gamma; + if (abs(denBW) < 1e-10) {return 0;} + + return cK_[i]*m2/denBW; +} + +EvtVector4C EvtWHad::WCurrent_KSK(const EvtVector4R& pKS, const EvtVector4R& pKplus) const { + + double s = (pKS + pKplus).mass2(); + EvtComplex f = BWKK(s, 0) + BWKK(s, 1) + BWKK(s, 2); + + return f*(pKS - pKplus); + +} + +EvtComplex EvtWHad::pcm(double s) const { + + double mpi2 = pow(0.140, 2); + if (abs(s) < 1e-10) return 0; + + double pcm2 = 1.0 - 4.0*mpi2/s; + EvtComplex result; + + if (pcm2 >= 0.0) { + result = EvtComplex(sqrt(pcm2), 0.0); + } else { + result = EvtComplex(0.0, sqrt(-pcm2)); + } + + return result; +} // =================== W+ -> pi_ current ======================================== -EvtVector4C EvtWHad::WCurrent(EvtVector4R q1) { - return q1; + +EvtVector4C EvtWHad::WCurrent(const EvtVector4R& q1) const { + return q1; } //====================== W+ -> pi+ pi0 current ========================================= -EvtVector4C EvtWHad::WCurrent(EvtVector4R q1, EvtVector4R q2) { - return BWr(q1+q2)*(q1-q2); + +EvtVector4C EvtWHad::WCurrent(const EvtVector4R& q1, const EvtVector4R& q2) const { + return BWr(q1 + q2)*(q1 - q2); } //========================= W+ -> pi+ pi+ pi- current ============================================== -EvtVector4C EvtWHad::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3) { - EvtVector4R Q=q1+q2+q3; - double Q2=Q.mass2(); - return BWa(Q)*( (q1-q3) - (Q*(Q*(q1-q3))/Q2)*BWr(q2+q3) + - (q2-q3) - (Q*(Q*(q2-q3))/Q2)*BWr(q1+q3) ); + +EvtVector4C EvtWHad::WCurrent(const EvtVector4R& q1, const EvtVector4R& q2, const EvtVector4R& q3) const { + + EvtVector4R Q = q1 + q2 + q3; + EvtVector4R q13 = q1 - q3, q23 = q2 - q3; + double Q2 = Q.mass2(); + + return BWa(Q)*(q13 - (Q*(Q*q13)/Q2) * BWr(q2 + q3) + q23 - (Q*(Q*q23)/Q2) * BWr(q1 + q3)); + } // ================= W+ -> pi+ pi+ pi- pi- pi+ current with symmetrization ================================ -EvtVector4C EvtWHad::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3, EvtVector4R q4, EvtVector4R q5) { - EvtVector4C V = JB(q1, q2, q3, q4, q5) + JB(q5, q2, q3, q4, q1) + JB(q1, q5, q3, q4, q2) + - JB(q1,q2,q4,q3,q5)+JB(q5,q2,q4,q3,q1)+JB(q1,q5,q4,q3,q2); - return V; -} +EvtVector4C EvtWHad::WCurrent(const EvtVector4R& q1, const EvtVector4R& q2, const EvtVector4R& q3, + const EvtVector4R& q4, const EvtVector4R& q5) const { -// =========================W+ -> K+ K- pi+ current ================================================== -EvtVector4C EvtWHad::WCurrent_KKP(EvtVector4R pKplus, EvtVector4R pKminus, EvtVector4R pPiPlus) { - const double mK_892 = 0.892, gammaK_892 = 0.051; - const double mA1 = 1.239, gammaA1 = 0.600; - - EvtVector4R q = pKplus + pKminus + pPiPlus; - double q2=q.mass2(); - EvtVector4R pK = pKminus + pPiPlus; - double pK2 = pK.mass2(); - - EvtComplex I(0., 1.), den1, den2; - - den1=1/(q2-mA1*mA1+I*mA1*gammaA1); - den2=1/(pK2-mK_892*mK_892+I*mK_892*gammaK_892); - - EvtTensor4C ten = EvtTensor4C::g()-(1./q2)*EvtGenFunctions::directProd(q,q); - EvtVector4C vec=den1*den2*(pKminus-pPiPlus); - vec = ten.cont2(vec); - return vec; + EvtVector4C term1 = JB(q1, q2, q3, q4, q5); + EvtVector4C term2 = JB(q5, q2, q3, q4, q1); + EvtVector4C term3 = JB(q1, q5, q3, q4, q2); + EvtVector4C term4 = JB(q1, q2, q4, q3, q5); + EvtVector4C term5 = JB(q5, q2, q4, q3, q1); + EvtVector4C term6 = JB(q1, q5, q4, q3, q2); + + EvtVector4C V = term1 + term2 + term3 + term4 + term5 + term6; + return V; } +// =========================W+ -> K+ K- pi+ current ================================================== -// hadronic current W+ -> K+ pi+ pi- -EvtVector4C EvtWHad::WCurrent_KPP(EvtVector4R pKplus, EvtVector4R pPiPlus, EvtVector4R pPiMinus) { - static double const cK1p=0.210709, cK1r=-0.0152997, cK2p=0.0945309, cK2r=0.504315; - const double mK1_1270 = 1.270, gammaK1_1270 = 0.090, gK1270_Krho = 2.71, gK1270_KsPi = 0.792; - const double mK1_1400 = 1.400, gammaK1_1400 = 0.174, gK11400_Krho = 0.254, gK11400_KsPi = 2.509; - const double mK_892 = 0.892, gammaK_892 = 0.051, gK892_Kpi = 3.26; - const double mRho = 0.770, gammaRho = 0.150, gRho_PiPi = 6.02; - - EvtVector4R q = pKplus + pPiPlus + pPiMinus; - double q2 = q.mass2(), pp2; - - EvtVector4C curr(0,0,0,0),curr1; - EvtComplex I(0., 1.), den1, den2; - - // W+ -> K1+(1270) -> K+ rho0 -> K+ pi+ pi- - den1 = gK1270_Krho/(q2 - mK1_1270*mK1_1270 + I*mK1_1270*gammaK1_1270); - pp2 = (pPiPlus+pPiMinus).mass2(); - den2 = gRho_PiPi/(pp2 - mRho*mRho + I*mRho*gammaRho); - curr1 = (pPiPlus - pPiMinus)*den1*den2; - curr = curr + cK1r*curr1; - - // W+ -> K1+(1270) -> K*(892)0 pi+ -> K+ pi- pi- - den1 = gK1270_KsPi/(q2 - mK1_1270*mK1_1270 + I*mK1_1270*gammaK1_1270); - pp2 = (pKplus+pPiMinus).mass2(); - den2 = gK892_Kpi/(pp2 - mK_892*mK_892 + I*mK_892*gammaK_892); - curr1 = (pKplus - pPiMinus)*den1*den2; - curr = curr + cK1p*curr1; - - // W+ -> K1+(1400) -> K+ rho0 -> K+ pi+ pi- - den1 = gK11400_Krho/(q2 - mK1_1400*mK1_1400 + I*mK1_1400*gammaK1_1400); - pp2 = (pPiMinus+pPiPlus).mass2(); - den2 = gRho_PiPi/(pp2 - mRho*mRho + I*mRho*gammaRho); - curr1 = (pPiPlus - pPiMinus)*den1*den2; - curr = curr + cK2r*curr1; - - // W+ -> K1+(1400) -> K*(892)0 pi+ -> K+ pi- pi+ - den1 = gK11400_KsPi/(q2 - mK1_1400*mK1_1400 + I*mK1_1400*gammaK1_1400); - pp2 = (pKplus+pPiMinus).mass2(); - den2 = gK892_Kpi/(pp2 - mK_892*mK_892 + I*mK_892*gammaK_892); - curr1 = (pKplus - pPiPlus)*den1*den2; - curr = curr + cK2p*curr1; - - EvtTensor4C ten = EvtTensor4C::g()-(1./q2)*EvtGenFunctions::directProd(q,q); - curr = ten.cont2(curr); - - return curr; +EvtVector4C EvtWHad::WCurrent_KKP(const EvtVector4R& pKplus, const EvtVector4R& pKminus, + const EvtVector4R& pPiPlus) const { + + const double mK_892(0.892), gammaK_892(0.051); + const double mA1(1.239), gammaA1(0.600); + + EvtVector4R q = pKplus + pKminus + pPiPlus; + double q2 = q.mass2(); + EvtVector4R pK = pKminus + pPiPlus; + double pK2 = pK.mass2(); + + EvtComplex I(0.0, 1.0), den1, den2; + + den1 = 1.0/(q2 - mA1*mA1 + I*mA1*gammaA1); + den2 = 1.0/(pK2 - mK_892*mK_892 + I*mK_892*gammaK_892); + + EvtTensor4C ten = EvtTensor4C::g() - (1.0/q2)*EvtGenFunctions::directProd(q, q); + + EvtVector4C vec = den1*den2*(pKminus - pPiPlus); + vec = ten.cont2(vec); + + return vec; } +// hadronic current W+ -> K+ pi+ pi- + +EvtVector4C EvtWHad::WCurrent_KPP(const EvtVector4R& pKplus, const EvtVector4R& pPiPlus, + const EvtVector4R& pPiMinus) const { + + const double cK1p = 0.210709, cK1r = -0.0152997, cK2p = 0.0945309, cK2r = 0.504315; + const double mK1_1270 = 1.270, gammaK1_1270 = 0.090, gK1270_Krho = 2.71, gK1270_KsPi = 0.792; + const double mK1_1400 = 1.400, gammaK1_1400 = 0.174, gK11400_Krho = 0.254, gK11400_KsPi = 2.509; + const double mK_892 = 0.892, gammaK_892 = 0.051, gK892_Kpi = 3.26; + const double mRho = 0.770, gammaRho = 0.150, gRho_PiPi = 6.02; + + EvtVector4R q = pKplus + pPiPlus + pPiMinus; + double q2 = q.mass2(), pp2(0.0); + + EvtVector4C curr(0, 0, 0, 0), curr1; + EvtComplex I(0.0, 1.0), den1, den2; + + // W+ -> K1+(1270) -> K+ rho0 -> K+ pi+ pi- + den1 = gK1270_Krho / (q2 - mK1_1270*mK1_1270 + I*mK1_1270*gammaK1_1270); + pp2 = (pPiPlus + pPiMinus).mass2(); + den2 = gRho_PiPi / (pp2 - mRho*mRho + I*mRho*gammaRho); + curr1 = (pPiPlus - pPiMinus) * den1*den2; + curr = curr + cK1r*curr1; + + // W+ -> K1+(1270) -> K*(892)0 pi+ -> K+ pi- pi- + den1 = gK1270_KsPi / (q2 - mK1_1270*mK1_1270 + I*mK1_1270*gammaK1_1270); + pp2 = (pKplus + pPiMinus).mass2(); + den2 = gK892_Kpi / (pp2 - mK_892*mK_892 + I*mK_892*gammaK_892); + curr1 = (pKplus - pPiMinus) * den1*den2; + curr = curr + cK1p*curr1; + + // W+ -> K1+(1400) -> K+ rho0 -> K+ pi+ pi- + den1 = gK11400_Krho / (q2 - mK1_1400*mK1_1400 + I*mK1_1400*gammaK1_1400); + pp2 = (pPiMinus + pPiPlus).mass2(); + den2 = gRho_PiPi / (pp2 - mRho*mRho + I*mRho*gammaRho); + curr1 = (pPiPlus - pPiMinus) * den1*den2; + curr = curr + cK2r*curr1; + + // W+ -> K1+(1400) -> K*(892)0 pi+ -> K+ pi- pi+ + den1 = gK11400_KsPi / (q2 - mK1_1400*mK1_1400 + I*mK1_1400*gammaK1_1400); + pp2 = (pKplus + pPiMinus).mass2(); + den2 = gK892_Kpi / (pp2 - mK_892*mK_892 + I*mK_892*gammaK_892); + curr1 = (pKplus - pPiPlus) * den1*den2; + curr = curr + cK2p*curr1; + + EvtTensor4C ten = EvtTensor4C::g() - (1.0/q2)*EvtGenFunctions::directProd(q, q); + curr = ten.cont2(curr); + + return curr; +} // a1 -> pi+ pi+ pi- BW -EvtComplex EvtWHad::BWa(EvtVector4R q) { - double const _mA1=1.26, _GA1=0.4; - EvtComplex I(0,1); - double Q2 = q.mass2(); - double GA1=_GA1*pi3G(Q2)/pi3G(_mA1*_mA1); - EvtComplex denBA1(_mA1*_mA1 - Q2,-1.*_mA1*GA1); - return _mA1*_mA1 / denBA1; + +EvtComplex EvtWHad::BWa(const EvtVector4R& q) const { + + double _mA1(1.26), _GA1(0.4); + double _mA1Sq = _mA1*_mA1; + EvtComplex I(0.0, 1.0); + double Q2 = q.mass2(); + double GA1 = _GA1*pi3G(Q2)/pi3G(_mA1Sq); + + EvtComplex denBA1(_mA1Sq - Q2, -1.0*_mA1*GA1); + + return _mA1Sq/denBA1; } +EvtComplex EvtWHad::BWf(const EvtVector4R& q) const { -EvtComplex EvtWHad::BWf(EvtVector4R q) { - double const mf=0.8, Gf=0.6; - EvtComplex I(0,1); - double Q2 = q.mass2(); - return mf*mf/(mf*mf-Q2-I*mf*Gf); + double mf(0.8), Gf(0.6); + double mfSq = mf*mf; + EvtComplex I(0.0, 1.0); + double Q2 = q.mass2(); + return mfSq/(mfSq - Q2 - I*mf*Gf); } -EvtComplex EvtWHad::BWr(EvtVector4R q) { - double _mRho = 0.775, _gammaRho=0.149, _mRhopr=1.364, _gammaRhopr=0.400, _beta=-0.108; - double m1=EvtPDL::getMeanMass(EvtPDL::getId("pi+")), m2=EvtPDL::getMeanMass(EvtPDL::getId("pi+")); - double mQ2=q.mass2(); - - // momenta in the rho->pipi decay - double dRho= _mRho*_mRho - m1*m1 - m2*m2; - double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2); - - double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2; - double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2); - - double dQ= mQ2 - m1*m1 - m2*m2; - double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2); - - - double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3); - EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho); - EvtComplex BRho= _mRho*_mRho / BRhoDem; - - double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3); - EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr); - EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem; - - return (BRho + _beta*BRhopr)/(1+_beta); +EvtComplex EvtWHad::BWr(const EvtVector4R& q) const { + + double _mRho(0.775), _gammaRho(0.149), _mRhopr(1.364), _gammaRhopr(0.400), _beta(-0.108); + double m1 = EvtPDL::getMeanMass(EvtPDL::getId("pi+")), m2 = EvtPDL::getMeanMass(EvtPDL::getId("pi+")); + double mQ2 = q.mass2(); + + double m1Sq = m1*m1, m2Sq = m2*m2, m12Sq = m1Sq*m2Sq; + + // momenta in the rho->pipi decay + double dRho = _mRho*_mRho - m1Sq - m2Sq; + double pPiRho = (1.0/_mRho) * sqrt((dRho*dRho)/4.0 - m12Sq); + + double dRhopr = _mRhopr*_mRhopr - m1Sq - m2Sq; + double pPiRhopr = (1.0/_mRhopr) * sqrt((dRhopr*dRhopr)/4.0 - m12Sq); + + double dQ = mQ2 - m1Sq - m2Sq; + double pPiQ = (1.0/sqrt(mQ2)) * sqrt((dQ*dQ)/4.0 - m12Sq); + + double gammaRho = _gammaRho*_mRho/sqrt(mQ2) * pow((pPiQ/pPiRho), 3); + EvtComplex BRhoDem(_mRho*_mRho - mQ2, -1.0*_mRho*gammaRho); + EvtComplex BRho = _mRho*_mRho/BRhoDem; + + double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2) * pow((pPiQ/pPiRhopr), 3); + EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2, -1.0*_mRho*gammaRhopr); + EvtComplex BRhopr = _mRhopr*_mRhopr/BRhoprDem; + + return (BRho + _beta*BRhopr)/(1.0 + _beta); } -double EvtWHad::pi3G(double m2) { - double mPi = EvtPDL::getMeanMass(EvtPDL::getId("pi+")); - double _mRho = 0.775; - if ( m2 > (_mRho+mPi) ) { - return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2)); - } - else { - double t1=m2-9.0*mPi*mPi; - return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1); - }; +double EvtWHad::pi3G(double m2) const { + + double mPi = EvtPDL::getMeanMass(EvtPDL::getId("pi+")); + double mRho(0.775); + double val(0.0); + + if (m2 > (mRho + mPi)) { + val = m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2)); + } else { + double t1 = m2 - 9.0 * mPi*mPi; + val = 4.1*pow(t1, 3.0)*(1.0 - 3.3*t1 + 5.8*t1*t1); + } + + return val; + } -EvtVector4C EvtWHad::JB( EvtVector4R p1, EvtVector4R p2, EvtVector4R p3, EvtVector4R p4, EvtVector4R p5) { - EvtVector4R Qtot = p1+p2+p3+p4+p5, Qa=p1+p2+p3; - EvtTensor4C T= (1/Qtot.mass2())*EvtGenFunctions::directProd(Qtot,Qtot) - EvtTensor4C::g(); - EvtVector4R V13 = Qa*( p2*(p1-p3) )/Qa.mass2() - (p1-p3); - EvtVector4R V23 = Qa*( p1*(p2-p3) )/Qa.mass2() - (p2-p3); - return BWa(Qtot)*BWa(Qa)*BWf(p4+p5)*( - T.cont1(V13)*BWr(p1+p3) + T.cont1(V23)*BWr(p2+p3) - ); +EvtVector4C EvtWHad::JB(const EvtVector4R& p1, const EvtVector4R& p2, const EvtVector4R& p3, + const EvtVector4R& p4, const EvtVector4R& p5) const { + + EvtVector4R Qtot = p1 + p2 + p3 + p4 + p5, Qa = p1 + p2 + p3; + EvtTensor4C T = (1.0/Qtot.mass2()) * EvtGenFunctions::directProd(Qtot, Qtot) - EvtTensor4C::g(); + + EvtVector4R p13 = p1 - p3, p23 = p2 - p3; + EvtVector4R V13 = Qa*(p2*p13)/Qa.mass2() - p13; + EvtVector4R V23 = Qa*(p1*p23)/Qa.mass2() - p23; + + return BWa(Qtot)*BWa(Qa)*BWf(p4 + p5)*(T.cont1(V13)*BWr(p1 + p3) + T.cont1(V23)*BWr(p2 + p3)); + }