From 59b975368c57973706e1f13d663d67ccc4d82ccd Mon Sep 17 00:00:00 2001 From: Dmitrii Pereima <dmitrii.pereima@cern.ch> Date: Tue, 3 Jan 2023 15:42:26 +0100 Subject: [PATCH 01/14] Add BC_VPPHAD model --- Gen/EvtGen/src/EvtWPPHad.cpp | 145 +++++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) create mode 100644 Gen/EvtGen/src/EvtWPPHad.cpp diff --git a/Gen/EvtGen/src/EvtWPPHad.cpp b/Gen/EvtGen/src/EvtWPPHad.cpp new file mode 100644 index 000000000..d209e13a4 --- /dev/null +++ b/Gen/EvtGen/src/EvtWPPHad.cpp @@ -0,0 +1,145 @@ +//-------------------------------------------------------------------------- +// +// 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: EvtWPPHad.cc +// +// Description: Routine to calculate W -> p ~p + (n pi) + (m K) current +// +// Modification history: +// A V Luchinsky Dec 2022, Module created +//------------------------------------------------------------------------ + +#include "EvtGenModels/EvtWPPHad.hh" +#include "EvtGenBase/EvtTensor4C.hh" + +#include "EvtGenBase/EvtPDL.hh" +#include "EvtGenBase/EvtTensor4C.hh" + +EvtWPPHad::EvtWPPHad() : + 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); +} + +double EvtWPPHad::pi3G(double Q2) const { + double mpi = EvtPDL::getMeanMass(EvtPDL::getId("pi+")); + double mpi2 = pow( mpi,2.); + double _mRho(0.775); + if (Q2 < pow(_mRho + mpi, 2.)) { + double arg = Q2 - 9. * mpi2; + return 4.1 * pow(arg, 3.) * (1. - 3.3 * arg + 5.8 * pow(arg, 2.)); + } else + return Q2 * (1.623 + 10.38 / Q2 - 9.32 / pow(Q2, 2.) + 0.65 / pow(Q2, 3.)); + +} + + +// a1 -> pi+ pi+ pi- BW +EvtComplex EvtWPPHad::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); + if(abs(denBA1)<1e-10) return 0; + return _mA1Sq/denBA1; +} + + +#define sp(p1, p2) ((p1)*(p2)) + + +EvtVector4C EvtWPPHad::WCurrent_ppPi(EvtVector4R p1, EvtDiracSpinor sp1, EvtVector4R p2, EvtDiracSpinor sp2, EvtVector4R k) { + EvtVector4C current; + double a = 1.25; + EvtVector4R q = p1 + p2 + k; + double q2 = q.mass2(); + double mp = p1.mass(), mp2 = mp*mp, mpi = k.mass(), mpi2 = mpi*mpi; + double mn = 939.656e-3, mn2 = mn*mn; + EvtComplex II(0, 1); + + double kp1=(k*p1); + double kp2=(k*p2); + double p1p2=(p1*p2); + + + double f1 = 1; + double f2 = 3.7/(2+mp); + double g1 = 1.25; + double g3 = 2*mp*g1/(p1p2+mp2); + + EvtComplex curS = EvtLeptonSCurrent(sp1, sp2); + EvtComplex curP = EvtLeptonPCurrent(sp1, sp2); + EvtVector4C curV = EvtLeptonVCurrent(sp1, sp2); + EvtVector4C curA = EvtLeptonACurrent(sp1, sp2); + EvtTensor4C curT = EvtLeptonTCurrent(sp1, sp2); + + double D1 = 1/(2*kp2 + mp2 - mn2 + mpi2); // (k+p2)^2-mn^2 + if(abs(D1)>1e1) { + return current; + } + +// 1 + current += curA*(-(f2*II)); +// 2 + current += (D1*g1*II)*curT.cont2(k); +// 3 + current += -(-0.5*(D1*f1))*dual(curT).cont2(k); +// 4 + current += -(D1*f2*II*mp)*dual(curT).cont2(k); +// 5 + current += k*(-(D1*g1))*curS; +// 6 + current += k*(-(D1*f1))*curP; +// 7 + current += k*(2*D1*f2*II*mp)*curP; +// 8 + current += k*(curV*k)*(D1*g3); +// 9 + current += k*(curA*k)*(D1*f2*II); +// 10 + current += p1*(curV*k)*(D1*g3); +// 11 + current += p1*(curA*k)*(-(D1*f2*II)); +// 12 + current += p2*(curV*k)*(D1*g3); +// 13 + current += p2*(curA*k)*(D1*f2*II); + + EvtTensor4C JT = (1/q2)*EvtGenFunctions::directProd(q, q) - EvtTensor4C::g(); + current = JT.cont2(current); + + // std::cout << abs2(current*current.conj()) << std::endl; + return BWa(q)*current; +} + -- GitLab From 3c4e193d7889b31bc6258753f1220f36687eb1e7 Mon Sep 17 00:00:00 2001 From: Dmitrii Pereima <dmitrii.pereima@cern.ch> Date: Tue, 3 Jan 2023 15:44:09 +0100 Subject: [PATCH 02/14] Add BC_VPPHAD model --- Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh | 72 ++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh new file mode 100644 index 000000000..44a97f590 --- /dev/null +++ b/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh @@ -0,0 +1,72 @@ +//-------------------------------------------------------------------------- +// +// 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: EvtBcVPPHad.hh +// +// Description: Module to implement Bc -> psi + p ~p + (n pi) + (m K) decays +// +// Modification history: +// +// A V Luchinsky Dec 2022, Created +//------------------------------------------------------------------------ + +#ifndef EvtBcVPPHad_HH +#define EvtBcVPPHad_HH + +#include "EvtGenBase/EvtVector4C.hh" +#include "EvtGenBase/EvtDecayAmp.hh" + +#include <string> + +class EvtBCVFF2; +class EvtParticle; +class EvtWPPHad; + +class EvtBcVPPHad : public EvtDecayAmp { +public: + + EvtBcVPPHad(); + virtual ~EvtBcVPPHad(); + + std::string getName(); + EvtDecayBase* clone(); + void initProbMax(); + void init(); + void decay(EvtParticle *p); + +protected: + + // Hadronic current function + EvtVector4C hardCurrPP(EvtParticle *root_particle, int i1, int i2) 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- + // 8 - K_S0 K+ + int out_code; + + EvtBCVFF2 *ffmodel; + EvtWPPHad *wcurr; + +}; + +#endif -- GitLab From 82af63da5757afac897404ab6d267b91f678b44a Mon Sep 17 00:00:00 2001 From: Dmitrii Pereima <dmitrii.pereima@cern.ch> Date: Tue, 3 Jan 2023 15:46:09 +0100 Subject: [PATCH 03/14] Add BC_VPPHAD model --- Gen/EvtGen/EvtGenModels/EvtWPPHad.hh | 50 ++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 Gen/EvtGen/EvtGenModels/EvtWPPHad.hh diff --git a/Gen/EvtGen/EvtGenModels/EvtWPPHad.hh b/Gen/EvtGen/EvtGenModels/EvtWPPHad.hh new file mode 100644 index 000000000..b00eb5886 --- /dev/null +++ b/Gen/EvtGen/EvtGenModels/EvtWPPHad.hh @@ -0,0 +1,50 @@ +//-------------------------------------------------------------------------- +// +// 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: EvtWPPHad.hh +// +// Description: Routine to calculate W -> p ~p + (n pi) + (m K) current +// +// Modification history: +// A V Luchinsky 20 Jan, 2013 Module created +// +//------------------------------------------------------------------------ + +#ifndef EvtWPPHad_HH +#define EvtWPPHad_HH + +#include "EvtGenBase/EvtComplex.hh" +#include "EvtGenBase/EvtVector4C.hh" +#include "EvtGenBase/EvtVector4R.hh" +#include "EvtGenBase/EvtTensor4C.hh" +#include "EvtGenBase/EvtDiracSpinor.hh" + +#include <vector> + +class EvtWPPHad { + +public: + + EvtWPPHad(); + + // a1 -> p+ p- pi+ + EvtVector4C WCurrent_ppPi(EvtVector4R p1, EvtDiracSpinor sp1, EvtVector4R p2, EvtDiracSpinor sp2, EvtVector4R k); + +protected: + + double pi3G(double Q2) const; + EvtComplex BWa(const EvtVector4R& q) const; + + +private: + + std::vector<double> mRho_, gamma0_, cK_; + +}; + +#endif -- GitLab From 3a1f00c4ccc46fe72d6a70fb8f004953b035f836 Mon Sep 17 00:00:00 2001 From: Dmitrii Pereima <dmitrii.pereima@cern.ch> Date: Tue, 3 Jan 2023 15:48:37 +0100 Subject: [PATCH 04/14] Add BC_VPPHAD model --- Gen/EvtGen/src/EvtBcVPPHad.cpp | 185 +++++++++++++++++++++++++++++++++ 1 file changed, 185 insertions(+) create mode 100644 Gen/EvtGen/src/EvtBcVPPHad.cpp diff --git a/Gen/EvtGen/src/EvtBcVPPHad.cpp b/Gen/EvtGen/src/EvtBcVPPHad.cpp new file mode 100644 index 000000000..fe8bd5f80 --- /dev/null +++ b/Gen/EvtGen/src/EvtBcVPPHad.cpp @@ -0,0 +1,185 @@ +//-------------------------------------------------------------------------- +// +// 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: EvtBcVPPHad.cc +// +// Description: Module to implement Bc -> psi + p ~p (n pi) + (m K) decays +// +// Modification history: +// +// A V Luchinsky Dec 2022 Module created +// +//------------------------------------------------------------------------ + + +// The following decays are supported: +// out_code=1: B_c+ -> V p ~p pi+ + + +#include "EvtGenModels/EvtBcVPPHad.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 "EvtGenBase/EvtDiracSpinor.hh" +#include "EvtGenModels/EvtBCVFF2.hh" +#include "EvtGenModels/EvtWPPHad.hh" + +#include <iostream> + +EvtBcVPPHad::EvtBcVPPHad() : + whichfit(0), + idVector(0), + out_code(0), + ffmodel(0), + wcurr(0) +{ +} + +EvtBcVPPHad::~EvtBcVPPHad() +{ + + if (ffmodel) {delete ffmodel;} + if (wcurr) {delete wcurr;} + +} + +std::string EvtBcVPPHad::getName() { + return "BC_VPPHAD"; +} + +EvtDecayBase* EvtBcVPPHad::clone() { + return new EvtBcVPPHad; +} + +//====================================================== + +void EvtBcVPPHad::init() { + + checkNArg(1); + checkSpinParent(EvtSpinType::SCALAR); + checkSpinDaughter(0, EvtSpinType::VECTOR); + checkSpinDaughter(1, EvtSpinType::DIRAC); + checkSpinDaughter(2, EvtSpinType::DIRAC); + for (int i = 3; i <= (getNDaug() - 1); i++) { + checkSpinDaughter(i, EvtSpinType::SCALAR); + } + + + idVector = getDaug(0).getId(); + whichfit = int(getArg(0) + 0.1); + ffmodel = new EvtBCVFF2(idVector, whichfit); + + wcurr = new EvtWPPHad(); + + // determine the code of final hadronic state + EvtIdSet thePis("pi+", "pi-", "pi0"); + EvtIdSet theK("K+", "K-", "K_S0"); + EvtIdSet thePs("p+", "anti-p-"); + + if(getNDaug() == 4 && thePs.contains(getDaug(1)) && thePs.contains(getDaug(2)) && thePis.contains(getDaug(3))) { + out_code = 1; // p+ p- pi+ + } + else + { + std::cout << "EvtBcPPHad::init --- unkonown decay" << std::endl; + }; + std::cout << "============ EvtBcVPPHad: out_code = " << out_code << " whichfit = "<< whichfit << std::endl; +} + +//====================================================== + +void EvtBcVPPHad::initProbMax() { + if(out_code == 1) { // p ~p pi + if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1) setProbMax(0.2e4); + else if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2) setProbMax(1000); + else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1) setProbMax(2000); + else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) setProbMax(7); + } + else { + std::cout << "probax: Have not yet implemented this final state in BC_VPPHAD model, out_code = " << out_code << std::endl; + ::abort(); + } +} + +//====================================================== +EvtVector4C EvtBcVPPHad::hardCurrPP(EvtParticle *root_particle, int i1, int i2) const { + EvtVector4C hardCur; + EvtVector4R p1 = root_particle->getDaug(1)->getP4(); + EvtDiracSpinor sp1 = root_particle->getDaug(1)->spParent(i1); + + EvtVector4R p2 = root_particle->getDaug(2) -> getP4(); + EvtDiracSpinor sp2 = root_particle->getDaug(2)->spParent(i2); + + if(out_code==1) { + EvtVector4R k = root_particle->getDaug(3) -> getP4(); + hardCur = wcurr ->WCurrent_ppPi(p1, sp1, p2, sp2, k); + } else { + std::cout << "hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; + for (int id = 0; id < (getNDaug() - 1); id++) { + // report(ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + } + ::abort(); + }; + return hardCur; +}; + + +//====================================================== +void EvtBcVPPHad::decay(EvtParticle *root_particle) { + + root_particle->initializePhaseSpace(getNDaug(), getDaugs()); + + + EvtParticle* Jpsi = root_particle->getDaug(0); + + EvtVector4R + 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 hardCur, Heps, eps; + EvtComplex amp; + for(int i1=0; i1<2; ++i1) { + for(int i2=0; i2<2; ++i2) { + hardCur = hardCurrPP(root_particle, i1, i2); + Heps = H.cont2(hardCur); + for (int i = 0; i < 4; i++) { + eps = Jpsi->epsParent(i).conj(); // psi-meson polarization vector + amp = eps*Heps; + vertex(i, i1, i2, amp); + } + }; + }; +}; + -- GitLab From aa3f276b12b6ae5a953c039c6d052f3c075bf9a1 Mon Sep 17 00:00:00 2001 From: Dmitrii Pereima <dmitrii.pereima@cern.ch> Date: Tue, 3 Jan 2023 15:51:00 +0100 Subject: [PATCH 05/14] Update Gen/EvtGen/src/EvtModelReg.cpp --- Gen/EvtGen/src/EvtModelReg.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Gen/EvtGen/src/EvtModelReg.cpp b/Gen/EvtGen/src/EvtModelReg.cpp index 8f07bb4b9..d44ba3d71 100755 --- a/Gen/EvtGen/src/EvtModelReg.cpp +++ b/Gen/EvtGen/src/EvtModelReg.cpp @@ -170,6 +170,7 @@ #include "EvtGenModels/EvtBToDiBaryonlnupQCD.hh" #include "EvtGenModels/EvtBcVHad.hh" +#include "EvtGenModels/EvtBcVPPHad.hh" #include "EvtGenModels/Evtbs2llGammaMNT.hh" #include "EvtGenModels/Evtbs2llGammaISRFSR.hh" @@ -354,6 +355,7 @@ EvtModelReg::EvtModelReg(const std::list<EvtDecayBase*>* extraModels) modelist.registerModel(new EvtGenericDalitz()); modelist.registerModel(new EvtBcVHad); + modelist.registerModel(new EvtBcVPPHad); modelist.registerModel(new Evtbs2llGammaMNT); modelist.registerModel(new Evtbs2llGammaISRFSR); -- GitLab From b8f981768fe40b1c01b8a600f69947a61937c864 Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Sat, 11 Feb 2023 22:35:37 +0100 Subject: [PATCH 06/14] Bug fix --- Gen/EvtGen/src/EvtBcVHad.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Gen/EvtGen/src/EvtBcVHad.cpp b/Gen/EvtGen/src/EvtBcVHad.cpp index 5ead598ee..b3d4b7c6a 100644 --- a/Gen/EvtGen/src/EvtBcVHad.cpp +++ b/Gen/EvtGen/src/EvtBcVHad.cpp @@ -230,7 +230,7 @@ EvtVector4C EvtBcVHad::hardCurr(EvtParticle *root_particle) const { } 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(), + hardCur = wcurr->WCurrent_5pi(root_particle->getDaug(1)->getP4(), root_particle->getDaug(2)->getP4(), root_particle->getDaug(3)->getP4(), root_particle->getDaug(4)->getP4(), -- GitLab From b4cd56de7d26896f487d591e739819781bb6017b Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Sat, 11 Feb 2023 22:43:11 +0100 Subject: [PATCH 07/14] Suggestions from Jhon --- Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh | 4 +- Gen/EvtGen/EvtGenModels/EvtWPPHad.hh | 50 --------- Gen/EvtGen/src/EvtBcVPPHad.cpp | 16 +-- Gen/EvtGen/src/EvtWHad.cpp | 72 +++++++++++- Gen/EvtGen/src/EvtWPPHad.cpp | 145 ------------------------- 5 files changed, 78 insertions(+), 209 deletions(-) delete mode 100644 Gen/EvtGen/EvtGenModels/EvtWPPHad.hh delete mode 100644 Gen/EvtGen/src/EvtWPPHad.cpp diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh index 44a97f590..b928c3928 100644 --- a/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtBcVPPHad.hh @@ -25,7 +25,7 @@ class EvtBCVFF2; class EvtParticle; -class EvtWPPHad; +class EvtWHad; class EvtBcVPPHad : public EvtDecayAmp { public: @@ -65,7 +65,7 @@ private: int out_code; EvtBCVFF2 *ffmodel; - EvtWPPHad *wcurr; + EvtWHad *wcurr; }; diff --git a/Gen/EvtGen/EvtGenModels/EvtWPPHad.hh b/Gen/EvtGen/EvtGenModels/EvtWPPHad.hh deleted file mode 100644 index b00eb5886..000000000 --- a/Gen/EvtGen/EvtGenModels/EvtWPPHad.hh +++ /dev/null @@ -1,50 +0,0 @@ -//-------------------------------------------------------------------------- -// -// 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: EvtWPPHad.hh -// -// Description: Routine to calculate W -> p ~p + (n pi) + (m K) current -// -// Modification history: -// A V Luchinsky 20 Jan, 2013 Module created -// -//------------------------------------------------------------------------ - -#ifndef EvtWPPHad_HH -#define EvtWPPHad_HH - -#include "EvtGenBase/EvtComplex.hh" -#include "EvtGenBase/EvtVector4C.hh" -#include "EvtGenBase/EvtVector4R.hh" -#include "EvtGenBase/EvtTensor4C.hh" -#include "EvtGenBase/EvtDiracSpinor.hh" - -#include <vector> - -class EvtWPPHad { - -public: - - EvtWPPHad(); - - // a1 -> p+ p- pi+ - EvtVector4C WCurrent_ppPi(EvtVector4R p1, EvtDiracSpinor sp1, EvtVector4R p2, EvtDiracSpinor sp2, EvtVector4R k); - -protected: - - double pi3G(double Q2) const; - EvtComplex BWa(const EvtVector4R& q) const; - - -private: - - std::vector<double> mRho_, gamma0_, cK_; - -}; - -#endif diff --git a/Gen/EvtGen/src/EvtBcVPPHad.cpp b/Gen/EvtGen/src/EvtBcVPPHad.cpp index fe8bd5f80..3099ce568 100644 --- a/Gen/EvtGen/src/EvtBcVPPHad.cpp +++ b/Gen/EvtGen/src/EvtBcVPPHad.cpp @@ -32,7 +32,7 @@ #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenModels/EvtBCVFF2.hh" -#include "EvtGenModels/EvtWPPHad.hh" +#include "EvtGenModels/EvtWHad.hh" #include <iostream> @@ -79,7 +79,7 @@ void EvtBcVPPHad::init() { whichfit = int(getArg(0) + 0.1); ffmodel = new EvtBCVFF2(idVector, whichfit); - wcurr = new EvtWPPHad(); + wcurr = new EvtWHad(); // determine the code of final hadronic state EvtIdSet thePis("pi+", "pi-", "pi0"); @@ -91,9 +91,9 @@ void EvtBcVPPHad::init() { } else { - std::cout << "EvtBcPPHad::init --- unkonown decay" << std::endl; + report(ERROR, "EvtBcPPHad::init --- unkonown decay"); }; - std::cout << "============ EvtBcVPPHad: out_code = " << out_code << " whichfit = "<< whichfit << std::endl; + report(INFO, "========== EvtBcVPPHad: out_code = ")<<out_code<<", whichfit = "<<whichfit<<std::endl; } //====================================================== @@ -106,7 +106,7 @@ void EvtBcVPPHad::initProbMax() { else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) setProbMax(7); } else { - std::cout << "probax: Have not yet implemented this final state in BC_VPPHAD model, out_code = " << out_code << std::endl; + report(ERROR,"probax: Have not yet implemented this final state in BC_VPPHAD model, out_code = ")<< out_code << std::endl; ::abort(); } } @@ -124,11 +124,7 @@ EvtVector4C EvtBcVPPHad::hardCurrPP(EvtParticle *root_particle, int i1, int i2) EvtVector4R k = root_particle->getDaug(3) -> getP4(); hardCur = wcurr ->WCurrent_ppPi(p1, sp1, p2, sp2, k); } else { - std::cout << "hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; - for (int id = 0; id < (getNDaug() - 1); id++) { - // report(ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; - } - ::abort(); + report(ERROR, "hardCurr: Have not yet implemented this final state in BC_VPPHAD model"); }; return hardCur; }; diff --git a/Gen/EvtGen/src/EvtWHad.cpp b/Gen/EvtGen/src/EvtWHad.cpp index ce32a77e0..b995dc2bd 100644 --- a/Gen/EvtGen/src/EvtWHad.cpp +++ b/Gen/EvtGen/src/EvtWHad.cpp @@ -258,7 +258,7 @@ EvtVector4C EvtWHad::WCurrent_7pi_nosymm(const EvtVector4R& p1, const EvtVector4 const EvtVector4R& p3, const EvtVector4R& p4, const EvtVector4R& p5, const EvtVector4R& p6, const EvtVector4R& p7) const { EvtVector4R qTot = p1+p2+p3+p4+p5+p6+p7; - EvtVector4R qa = p1+p2+p3+p5+p6; + //EvtVector4R qa = p1+p2+p3+p5+p6; EvtVector4C eps1 = EvtWHad::WCurrent_5pi(p1, p2, p3, p5, p6); // pi+ pi+ pi+ pi- p- EvtVector4R pf0 = p4 + p7; return eps1*BWa(qTot)*BWf(pf0); @@ -402,7 +402,7 @@ EvtVector4C EvtWHad::WCurrent_K4pi_nosymm(const EvtVector4R& p1, const EvtVector const EvtVector4R& p3, const EvtVector4R& p4, const EvtVector4R& p5) const { const EvtComplex I(0, 1); - EvtVector4R pKstar = p1 + p4, pa1 = p2+p3+p5, pTot = pKstar + pa1; + EvtVector4R pKstar = p1 + p4, pa1 = p2+p3+p5;//, pTot = pKstar + pa1; EvtComplex denKstar = pKstar*pKstar - mK_[0]*mK_[0] + I*mK_[0]*gammaK_[0]; //K(892) if(abs(denKstar)<1e-10) { denKstar = 1e10; @@ -414,3 +414,71 @@ EvtVector4C EvtWHad::WCurrent_K4pi_nosymm(const EvtVector4R& p1, const EvtVector return eps; } +#define sp(p1, p2) ((p1)*(p2)) + + +EvtVector4C EvtWHad::WCurrent_ppPi(const EvtVector4R& p1, const EvtDiracSpinor& sp1, const EvtVector4R& p2, const EvtDiracSpinor& sp2, const EvtVector4R& k) const { + EvtVector4C current; + //double a = 1.25; + EvtVector4R q = p1 + p2 + k; + double q2 = q.mass2(); + double mp = p1.mass(), mp2 = mp*mp, mpi = k.mass(), mpi2 = mpi*mpi; + const double mn = EvtPDL::getMeanMass( EvtPDL::getId( "n0" ) ) ; + double mn2 = mn*mn; + EvtComplex II(0, 1); + + //double kp1=(k*p1); + double kp2=(k*p2); + double p1p2=(p1*p2); + + + double f1 = 1; + double f2 = 3.7/(2+mp); + double g1 = 1.25; + double g3 = 2*mp*g1/(p1p2+mp2); + + EvtComplex curS = EvtLeptonSCurrent(sp1, sp2); + EvtComplex curP = EvtLeptonPCurrent(sp1, sp2); + EvtVector4C curV = EvtLeptonVCurrent(sp1, sp2); + EvtVector4C curA = EvtLeptonACurrent(sp1, sp2); + EvtTensor4C curT = EvtLeptonTCurrent(sp1, sp2); + + double D1 = 1/(2*kp2 + mp2 - mn2 + mpi2); // (k+p2)^2-mn^2 + if(abs(D1)>1e1) { + return current; + } + +// 1 + current += curA*(-(f2*II)); +// 2 + current += (D1*g1*II)*curT.cont2(k); +// 3 + current += -(-0.5*(D1*f1))*dual(curT).cont2(k); +// 4 + current += -(D1*f2*II*mp)*dual(curT).cont2(k); +// 5 + current += k*(-(D1*g1))*curS; +// 6 + current += k*(-(D1*f1))*curP; +// 7 + current += k*(2*D1*f2*II*mp)*curP; +// 8 + current += k*(curV*k)*(D1*g3); +// 9 + current += k*(curA*k)*(D1*f2*II); +// 10 + current += p1*(curV*k)*(D1*g3); +// 11 + current += p1*(curA*k)*(-(D1*f2*II)); +// 12 + current += p2*(curV*k)*(D1*g3); +// 13 + current += p2*(curA*k)*(D1*f2*II); + + EvtTensor4C JT = (1/q2)*EvtGenFunctions::directProd(q, q) - EvtTensor4C::g(); + current = JT.cont2(current); + + // std::cout << abs2(current*current.conj()) << std::endl; + return BWa(q)*current; +} + diff --git a/Gen/EvtGen/src/EvtWPPHad.cpp b/Gen/EvtGen/src/EvtWPPHad.cpp deleted file mode 100644 index d209e13a4..000000000 --- a/Gen/EvtGen/src/EvtWPPHad.cpp +++ /dev/null @@ -1,145 +0,0 @@ -//-------------------------------------------------------------------------- -// -// 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: EvtWPPHad.cc -// -// Description: Routine to calculate W -> p ~p + (n pi) + (m K) current -// -// Modification history: -// A V Luchinsky Dec 2022, Module created -//------------------------------------------------------------------------ - -#include "EvtGenModels/EvtWPPHad.hh" -#include "EvtGenBase/EvtTensor4C.hh" - -#include "EvtGenBase/EvtPDL.hh" -#include "EvtGenBase/EvtTensor4C.hh" - -EvtWPPHad::EvtWPPHad() : - 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); -} - -double EvtWPPHad::pi3G(double Q2) const { - double mpi = EvtPDL::getMeanMass(EvtPDL::getId("pi+")); - double mpi2 = pow( mpi,2.); - double _mRho(0.775); - if (Q2 < pow(_mRho + mpi, 2.)) { - double arg = Q2 - 9. * mpi2; - return 4.1 * pow(arg, 3.) * (1. - 3.3 * arg + 5.8 * pow(arg, 2.)); - } else - return Q2 * (1.623 + 10.38 / Q2 - 9.32 / pow(Q2, 2.) + 0.65 / pow(Q2, 3.)); - -} - - -// a1 -> pi+ pi+ pi- BW -EvtComplex EvtWPPHad::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); - if(abs(denBA1)<1e-10) return 0; - return _mA1Sq/denBA1; -} - - -#define sp(p1, p2) ((p1)*(p2)) - - -EvtVector4C EvtWPPHad::WCurrent_ppPi(EvtVector4R p1, EvtDiracSpinor sp1, EvtVector4R p2, EvtDiracSpinor sp2, EvtVector4R k) { - EvtVector4C current; - double a = 1.25; - EvtVector4R q = p1 + p2 + k; - double q2 = q.mass2(); - double mp = p1.mass(), mp2 = mp*mp, mpi = k.mass(), mpi2 = mpi*mpi; - double mn = 939.656e-3, mn2 = mn*mn; - EvtComplex II(0, 1); - - double kp1=(k*p1); - double kp2=(k*p2); - double p1p2=(p1*p2); - - - double f1 = 1; - double f2 = 3.7/(2+mp); - double g1 = 1.25; - double g3 = 2*mp*g1/(p1p2+mp2); - - EvtComplex curS = EvtLeptonSCurrent(sp1, sp2); - EvtComplex curP = EvtLeptonPCurrent(sp1, sp2); - EvtVector4C curV = EvtLeptonVCurrent(sp1, sp2); - EvtVector4C curA = EvtLeptonACurrent(sp1, sp2); - EvtTensor4C curT = EvtLeptonTCurrent(sp1, sp2); - - double D1 = 1/(2*kp2 + mp2 - mn2 + mpi2); // (k+p2)^2-mn^2 - if(abs(D1)>1e1) { - return current; - } - -// 1 - current += curA*(-(f2*II)); -// 2 - current += (D1*g1*II)*curT.cont2(k); -// 3 - current += -(-0.5*(D1*f1))*dual(curT).cont2(k); -// 4 - current += -(D1*f2*II*mp)*dual(curT).cont2(k); -// 5 - current += k*(-(D1*g1))*curS; -// 6 - current += k*(-(D1*f1))*curP; -// 7 - current += k*(2*D1*f2*II*mp)*curP; -// 8 - current += k*(curV*k)*(D1*g3); -// 9 - current += k*(curA*k)*(D1*f2*II); -// 10 - current += p1*(curV*k)*(D1*g3); -// 11 - current += p1*(curA*k)*(-(D1*f2*II)); -// 12 - current += p2*(curV*k)*(D1*g3); -// 13 - current += p2*(curA*k)*(D1*f2*II); - - EvtTensor4C JT = (1/q2)*EvtGenFunctions::directProd(q, q) - EvtTensor4C::g(); - current = JT.cont2(current); - - // std::cout << abs2(current*current.conj()) << std::endl; - return BWa(q)*current; -} - -- GitLab From 770e079d1a4f6447c4e228e1b537d08ab26e9fc8 Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Sun, 12 Feb 2023 14:25:50 +0100 Subject: [PATCH 08/14] Update EvtWHad --- Gen/EvtGen/EvtGenModels/EvtWHad.hh | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Gen/EvtGen/EvtGenModels/EvtWHad.hh b/Gen/EvtGen/EvtGenModels/EvtWHad.hh index 21fd7cf6e..2b2198a34 100644 --- a/Gen/EvtGen/EvtGenModels/EvtWHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtWHad.hh @@ -23,6 +23,7 @@ #include "EvtGenBase/EvtVector4C.hh" #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenBase/EvtTensor4C.hh" +#include "EvtGenBase/EvtDiracSpinor.hh" #include <vector> @@ -64,6 +65,9 @@ public: // 1=K+ 2 = pi+ 3 = pi+ 4 = pi- 5 = pi- with symmetrizatiom EvtVector4C WCurrent_K4pi(const EvtVector4R& p1, const EvtVector4R& p2, const EvtVector4R& p3, const EvtVector4R& p4, const EvtVector4R& p5) const; + + // a1 -> p+ p- pi+ + EvtVector4C WCurrent_ppPi(const EvtVector4R& p1, const EvtDiracSpinor& sp1, const EvtVector4R& p2, const EvtDiracSpinor& sp2, const EvtVector4R& k) const; -- GitLab From 73b979d6831346c02dd90bef584e1cd87faf374b Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Tue, 14 Feb 2023 17:22:11 +0100 Subject: [PATCH 09/14] Code polishing fro John --- Gen/EvtGen/src/EvtBcVPPHad.cpp | 124 ++++++++++++++++----------------- Gen/EvtGen/src/EvtWHad.cpp | 109 ++++++++++++----------------- 2 files changed, 107 insertions(+), 126 deletions(-) diff --git a/Gen/EvtGen/src/EvtBcVPPHad.cpp b/Gen/EvtGen/src/EvtBcVPPHad.cpp index 3099ce568..bff73cd7a 100644 --- a/Gen/EvtGen/src/EvtBcVPPHad.cpp +++ b/Gen/EvtGen/src/EvtBcVPPHad.cpp @@ -47,10 +47,8 @@ EvtBcVPPHad::EvtBcVPPHad() : EvtBcVPPHad::~EvtBcVPPHad() { - if (ffmodel) {delete ffmodel;} if (wcurr) {delete wcurr;} - } std::string EvtBcVPPHad::getName() { @@ -74,7 +72,6 @@ void EvtBcVPPHad::init() { checkSpinDaughter(i, EvtSpinType::SCALAR); } - idVector = getDaug(0).getId(); whichfit = int(getArg(0) + 0.1); ffmodel = new EvtBCVFF2(idVector, whichfit); @@ -82,100 +79,101 @@ void EvtBcVPPHad::init() { wcurr = new EvtWHad(); // determine the code of final hadronic state - EvtIdSet thePis("pi+", "pi-", "pi0"); - EvtIdSet theK("K+", "K-", "K_S0"); - EvtIdSet thePs("p+", "anti-p-"); + const EvtIdSet thePis("pi+", "pi-", "pi0"); + const EvtIdSet theK("K+", "K-", "K_S0"); + const EvtIdSet thePs("p+", "anti-p-"); - if(getNDaug() == 4 && thePs.contains(getDaug(1)) && thePs.contains(getDaug(2)) && thePis.contains(getDaug(3))) { - out_code = 1; // p+ p- pi+ + if ( getNDaug() == 4 && thePs.contains(getDaug(1)) && + thePs.contains(getDaug(2)) && thePis.contains(getDaug(3)) ) { + out_code = 1; // p+ p- pi+ + } else { + report(ERROR, "EvtBcPPHad::init has unknown decay mode"); } - else - { - report(ERROR, "EvtBcPPHad::init --- unkonown decay"); - }; - report(INFO, "========== EvtBcVPPHad: out_code = ")<<out_code<<", whichfit = "<<whichfit<<std::endl; + report(INFO, "EvtBcVPPHad")<<": out_code = "<<out_code<<", whichfit = "<<whichfit<<std::endl; } //====================================================== void EvtBcVPPHad::initProbMax() { - if(out_code == 1) { // p ~p pi - if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1) setProbMax(0.2e4); - else if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2) setProbMax(1000); - else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1) setProbMax(2000); - else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) setProbMax(7); - } - else { - report(ERROR,"probax: Have not yet implemented this final state in BC_VPPHAD model, out_code = ")<< out_code << std::endl; - ::abort(); - } + + if (out_code == 1) { + // p ~p pi + if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1) { + setProbMax(2000.0); + } else if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2) { + setProbMax(1000.0); + } else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1) { + setProbMax(2000.0); + } else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) { + setProbMax(7.0); + } + } else { + report(ERROR,"probax: Have not yet implemented this final state in BC_VPPHAD model, out_code = ") + << out_code << std::endl; + ::abort(); + } } //====================================================== EvtVector4C EvtBcVPPHad::hardCurrPP(EvtParticle *root_particle, int i1, int i2) const { + EvtVector4C hardCur; - EvtVector4R p1 = root_particle->getDaug(1)->getP4(); - EvtDiracSpinor sp1 = root_particle->getDaug(1)->spParent(i1); + const EvtVector4R p1 = root_particle->getDaug(1)->getP4(); + const EvtDiracSpinor sp1 = root_particle->getDaug(1)->spParent(i1); - EvtVector4R p2 = root_particle->getDaug(2) -> getP4(); - EvtDiracSpinor sp2 = root_particle->getDaug(2)->spParent(i2); + const EvtVector4R p2 = root_particle->getDaug(2)->getP4(); + const EvtDiracSpinor sp2 = root_particle->getDaug(2)->spParent(i2); - if(out_code==1) { - EvtVector4R k = root_particle->getDaug(3) -> getP4(); - hardCur = wcurr ->WCurrent_ppPi(p1, sp1, p2, sp2, k); - } else { - report(ERROR, "hardCurr: Have not yet implemented this final state in BC_VPPHAD model"); - }; + if (out_code == 1) { + const EvtVector4R k = root_particle->getDaug(3)->getP4(); + hardCur = wcurr->WCurrent_ppPi(p1, sp1, p2, sp2, k); + } return hardCur; -}; - +} //====================================================== void EvtBcVPPHad::decay(EvtParticle *root_particle) { - root_particle->initializePhaseSpace(getNDaug(), getDaugs()); + root_particle->initializePhaseSpace(getNDaug(), getDaugs()); + EvtParticle* Jpsi = root_particle->getDaug(0); - EvtParticle* Jpsi = root_particle->getDaug(0); - - EvtVector4R - 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(); + const EvtVector4R + p4b(root_particle->mass(), 0., 0., 0.), // Bc momentum + p4meson = Jpsi->getP4(), // J/psi momenta + Q = p4b - p4meson, + p4Sum = p4meson + p4b; + const 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; + const double m_meson = Jpsi->mass(); + const double m_b = root_particle->mass(); + const 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; + const 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 hardCur, Heps, eps; - EvtComplex amp; - for(int i1=0; i1<2; ++i1) { - for(int i2=0; i2<2; ++i2) { - hardCur = hardCurrPP(root_particle, i1, i2); - Heps = H.cont2(hardCur); - for (int i = 0; i < 4; i++) { - eps = Jpsi->epsParent(i).conj(); // psi-meson polarization vector - amp = eps*Heps; - vertex(i, i1, i2, amp); - } - }; - }; -}; + for (int i1 = 0; i1 < 2; ++i1) { + for (int i2 = 0; i2 < 2; ++i2) { + const EvtVector4C hardCur = hardCurrPP(root_particle, i1, i2); + const EvtVector4C Heps = H.cont2(hardCur); + for (int i = 0; i < 4; i++) { + // psi-meson polarization vector + const EvtVector4C eps = Jpsi->epsParent(i).conj(); + const EvtComplex amp = eps*Heps; + vertex(i, i1, i2, amp); + } + } + } +} diff --git a/Gen/EvtGen/src/EvtWHad.cpp b/Gen/EvtGen/src/EvtWHad.cpp index b995dc2bd..ffe545ec8 100644 --- a/Gen/EvtGen/src/EvtWHad.cpp +++ b/Gen/EvtGen/src/EvtWHad.cpp @@ -262,7 +262,7 @@ EvtVector4C EvtWHad::WCurrent_7pi_nosymm(const EvtVector4R& p1, const EvtVector4 EvtVector4C eps1 = EvtWHad::WCurrent_5pi(p1, p2, p3, p5, p6); // pi+ pi+ pi+ pi- p- EvtVector4R pf0 = p4 + p7; return eps1*BWa(qTot)*BWf(pf0); -}; +} // hadronic current W+ -> K+ pi+ pi- @@ -414,71 +414,54 @@ EvtVector4C EvtWHad::WCurrent_K4pi_nosymm(const EvtVector4R& p1, const EvtVector return eps; } -#define sp(p1, p2) ((p1)*(p2)) +EvtVector4C EvtWHad::WCurrent_ppPi(const EvtVector4R& p1, const EvtDiracSpinor& sp1, const EvtVector4R& p2, + const EvtDiracSpinor& sp2, const EvtVector4R& k) const { + const EvtVector4R q = p1 + p2 + k; + const double q2 = q.mass2(); + const double mp = p1.mass(), mp2 = mp*mp, mpi = k.mass(), mpi2 = mpi*mpi; + const double mn = EvtPDL::getMeanMass( EvtPDL::getId( "n0" ) ); + const double mn2 = mn*mn; + const EvtComplex II(0, 1); + + const double kp2 = k*p2; + const double p1p2 = p1*p2; + + const double f1 = 1.0; + const double f2 = 3.7/(2.0 + mp); + const double g1 = 1.25; + const double g3 = 2.0*mp*g1/(p1p2 + mp2); + + const EvtComplex curS = EvtLeptonSCurrent(sp1, sp2); + const EvtComplex curP = EvtLeptonPCurrent(sp1, sp2); + const EvtVector4C curV = EvtLeptonVCurrent(sp1, sp2); + const EvtVector4C curA = EvtLeptonACurrent(sp1, sp2); + const EvtTensor4C curT = EvtLeptonTCurrent(sp1, sp2); + + const double D1 = 1/(2*kp2 + mp2 - mn2 + mpi2); //(k+p2)^2-mn^2 + // Amplitude: ~U(p1).GA5.Vpp(alpha).prop(-p2-k).V(p2) + permutations + // U() and V() are proton and antiproton spinors, prop() is the proton's propagator, + // and Vpp(alpha) is the W->pp vertex. + // Expand terms to use basic spinor currents (Scalar, Vector, Axial, etc) -EvtVector4C EvtWHad::WCurrent_ppPi(const EvtVector4R& p1, const EvtDiracSpinor& sp1, const EvtVector4R& p2, const EvtDiracSpinor& sp2, const EvtVector4R& k) const { EvtVector4C current; - //double a = 1.25; - EvtVector4R q = p1 + p2 + k; - double q2 = q.mass2(); - double mp = p1.mass(), mp2 = mp*mp, mpi = k.mass(), mpi2 = mpi*mpi; - const double mn = EvtPDL::getMeanMass( EvtPDL::getId( "n0" ) ) ; - double mn2 = mn*mn; - EvtComplex II(0, 1); - - //double kp1=(k*p1); - double kp2=(k*p2); - double p1p2=(p1*p2); - - - double f1 = 1; - double f2 = 3.7/(2+mp); - double g1 = 1.25; - double g3 = 2*mp*g1/(p1p2+mp2); - - EvtComplex curS = EvtLeptonSCurrent(sp1, sp2); - EvtComplex curP = EvtLeptonPCurrent(sp1, sp2); - EvtVector4C curV = EvtLeptonVCurrent(sp1, sp2); - EvtVector4C curA = EvtLeptonACurrent(sp1, sp2); - EvtTensor4C curT = EvtLeptonTCurrent(sp1, sp2); - - double D1 = 1/(2*kp2 + mp2 - mn2 + mpi2); // (k+p2)^2-mn^2 - if(abs(D1)>1e1) { - return current; - } - -// 1 - current += curA*(-(f2*II)); -// 2 - current += (D1*g1*II)*curT.cont2(k); -// 3 - current += -(-0.5*(D1*f1))*dual(curT).cont2(k); -// 4 - current += -(D1*f2*II*mp)*dual(curT).cont2(k); -// 5 - current += k*(-(D1*g1))*curS; -// 6 - current += k*(-(D1*f1))*curP; -// 7 - current += k*(2*D1*f2*II*mp)*curP; -// 8 - current += k*(curV*k)*(D1*g3); -// 9 - current += k*(curA*k)*(D1*f2*II); -// 10 - current += p1*(curV*k)*(D1*g3); -// 11 - current += p1*(curA*k)*(-(D1*f2*II)); -// 12 - current += p2*(curV*k)*(D1*g3); -// 13 - current += p2*(curA*k)*(D1*f2*II); - - EvtTensor4C JT = (1/q2)*EvtGenFunctions::directProd(q, q) - EvtTensor4C::g(); - current = JT.cont2(current); - - // std::cout << abs2(current*current.conj()) << std::endl; + current += curA*(-(f2*II)); + current += (D1*g1*II)*curT.cont2(k); + current += -(-0.5*(D1*f1))*dual(curT).cont2(k); + current += -(D1*f2*II*mp)*dual(curT).cont2(k); + current += k*(-(D1*g1))*curS; + current += k*(-(D1*f1))*curP; + current += k*(2*D1*f2*II*mp)*curP; + current += k*(curV*k)*(D1*g3); + current += k*(curA*k)*(D1*f2*II); + current += p1*(curV*k)*(D1*g3); + current += p1*(curA*k)*(-(D1*f2*II)); + current += p2*(curV*k)*(D1*g3); + current += p2*(curA*k)*(D1*f2*II); + + const EvtTensor4C JT = (1.0/q2)*EvtGenFunctions::directProd(q, q) - EvtTensor4C::g(); + current = JT.cont2(current); + return BWa(q)*current; } -- GitLab From b7874563eaa2ed2582a8004f83b4712499be9ee6 Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Sun, 7 May 2023 09:29:54 +0200 Subject: [PATCH 10/14] charge order fix --- Gen/EvtGen/EvtGenModels/EvtBcVHad.hh | 12 +- Gen/EvtGen/src/EvtBcVHad.cpp | 269 +++++++++++++++++---------- 2 files changed, 184 insertions(+), 97 deletions(-) diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh index c21fe0a24..9f2780933 100644 --- a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh @@ -45,6 +45,7 @@ protected: // Hadronic current function EvtVector4C hardCurr(EvtParticle *root_particle) const; + void parseDecay(void); private: // whichfit --- code of the Bc -> VW formfactor set: @@ -66,9 +67,18 @@ private: // 8 - K_S0 K+ int out_code; + + EvtBCVFF2 *ffmodel; EvtWHad *wcurr; + int iPiPlus[4] = {-1, -1, -1, -1}, + iPiMinus[4] = {-1, -1, -1, -1}, + iPiZero[4] = {-1,-1,-1,-1}, + iKPlus[4] = {-1,-1,-1,-1}, + iKMinus[4] = {-1,-1,-1,-1}; + + }; -#endif +#endif \ No newline at end of file diff --git a/Gen/EvtGen/src/EvtBcVHad.cpp b/Gen/EvtGen/src/EvtBcVHad.cpp index b3d4b7c6a..1ebb101de 100644 --- a/Gen/EvtGen/src/EvtBcVHad.cpp +++ b/Gen/EvtGen/src/EvtBcVHad.cpp @@ -21,15 +21,15 @@ // The following decays are supported: // out_code=1: B_c+ -> V pi+ -// out_code=2: B_c+ -> V pi+ pi0 -// out_code=3: B_c+ -> V pi+ pi+ pi- -// out_code=5: B_c+ -> V pi+ pi+ pi- pi- pi+ -// out_code=6: B_c+ -> V K+ K- pi+ -// out_code=7: B_c+ -> V K+ pi+ pi- +// out_code=2: B_c+ -> V pi+ pi0 (auto sym) +// out_code=3: B_c+ -> V + 2pi+ pi- (auto sym) +// out_code=5: B_c+ -> V + 3pi+ + 2 pi- (auto sym) +// out_code=6: B_c+ -> V K+ K- pi+ (auto sym) +// out_code=7: B_c+ -> V K+ pi+ pi- (auto sym) // out_code=8: B_c+ -> V K_S0 K+ -// out_code=9: B_c+ -> V K+ K- pi+ pi+ pi- -// out_code=10: B_c+ -> V pi+ pi+ pi+ pi+ pi- pi- pi- (7 charged pi) -// out_code=11: B_c+ -> V K+ pi+ pi+ pi- pi- +// out_code=9: B_c+ -> V K+ K- pi+ pi+ pi- (auto sym) +// out_code=10: B_c+ -> V + 4 pi+ + 3 pi- (auto sym) +// out_code=11: B_c+ -> V K+ pi+ pi+ pi- pi- (auto sym) #include "EvtGenModels/EvtBcVHad.hh" @@ -46,7 +46,7 @@ #include "EvtGenModels/EvtWHad.hh" #include <iostream> - + EvtBcVHad::EvtBcVHad() : whichfit(0), idVector(0), @@ -72,70 +72,147 @@ EvtDecayBase* EvtBcVHad::clone() { return new EvtBcVHad; } -//====================================================== - -void EvtBcVHad::init() { - - checkNArg(1); - checkSpinParent(EvtSpinType::SCALAR); - checkSpinDaughter(0, EvtSpinType::VECTOR); - for (int i = 1; i <= (getNDaug() - 1); i++) { - checkSpinDaughter(i, EvtSpinType::SCALAR); - } - - - 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 +void EvtBcVHad::parseDecay() +{ EvtIdSet thePis("pi+", "pi-", "pi0"); EvtIdSet theK("K+", "K-", "K_S0"); - if (getNDaug() == 2 && thePis.contains(getDaug(1))) { + EvtIdSet *PiPlusID, *PiMinusID, *PiZeroID, *KPlusID, *KMinusID, + *BcPlusID = new EvtIdSet("B_c+"), *BcMinusID = new EvtIdSet("B_c-"); + report(INFO,"getParentId()=")<<getParentId()<<std::endl; + + if(BcPlusID->contains(getParentId())) + { + PiPlusID = new EvtIdSet("pi+"); + PiMinusID = new EvtIdSet("pi-"); + PiZeroID = new EvtIdSet("pi0"); + KPlusID = new EvtIdSet("K+"); + KMinusID = new EvtIdSet("K-"); + } + else if(BcMinusID->contains(getParentId())) + { + PiPlusID = new EvtIdSet("pi-"); + PiMinusID = new EvtIdSet("pi+"); + PiZeroID = new EvtIdSet("pi0"); + KPlusID = new EvtIdSet("K-"); + KMinusID = new EvtIdSet("K+"); + }; + + + int PiPlusFound = 0, PiMinusFound = 0, PiZeroFound = 0, KPlusFound = 0, KMinusFound = 0; + for(int iDaughter = 0; iDaughter<getNDaug(); ++iDaughter) + { + report(INFO, "iDaughter=")<<iDaughter<<" id="<<getDaug(iDaughter).getName()<<std::endl; + if(PiPlusID->contains(getDaug(iDaughter))) + { + iPiPlus[PiPlusFound] = iDaughter; + PiPlusFound++; + } + else if(PiMinusID->contains(getDaug(iDaughter))) + { + iPiMinus[PiMinusFound] = iDaughter; + PiMinusFound++; + } + else if(PiZeroID->contains(getDaug(iDaughter))) + { + iPiZero[PiZeroFound] = iDaughter; + PiZeroFound++; + } + else if(KPlusID->contains(getDaug(iDaughter))) + { + iKPlus[KPlusFound] = iDaughter; + KPlusFound++; + } + else if(KMinusID->contains(getDaug(iDaughter))) + { + iKMinus[KMinusFound] = iDaughter; + KMinusFound++; + }; + }; + + + if (getNDaug() == 2 && PiPlusFound == 1) + { out_code = 1; // pi+ - } else if (getNDaug() == 3 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2))) { + } + else if (getNDaug() == 3 && PiPlusFound == 1 && PiZeroFound==1) + { out_code = 2; // pi+ pi0 - } else if (getNDaug() == 4 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) && - thePis.contains(getDaug(3))) { - out_code = 3;// pi+ pi+ pi- - } else if (getNDaug() == 5 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) && - thePis.contains(getDaug(3)) && thePis.contains(getDaug(4))) { - out_code = 4; //4pi - } else if (getNDaug() == 6 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) && - thePis.contains(getDaug(3)) && thePis.contains(getDaug(4)) && thePis.contains(getDaug(5))) { + } + else if (getNDaug() == 4 && PiPlusFound == 2 && PiMinusFound == 1) + { + out_code = 3;// pi+ pi+ pi- + } + else if (getNDaug() == 5 && PiPlusFound == 2 && PiMinusFound==1 && PiZeroFound==1) + { + out_code = 4; // pi+ pi+ pi- pi0 + } + else if (getNDaug() == 6 && PiPlusFound==3 && PiMinusFound == 2) + { out_code = 5; //5pi - } else if (getNDaug() == 4 && theK.contains(getDaug(1)) && theK.contains(getDaug(2)) && - thePis.contains(getDaug(3))) { + } else if (getNDaug() == 4 && KPlusFound == 1 && KMinusFound==1 && PiPlusFound == 1) + { out_code = 6; //KKpi - } else if (getNDaug() == 4 && theK.contains(getDaug(1)) && thePis.contains(getDaug(2)) && - thePis.contains(getDaug(3))) { - out_code = 7; // Kpi pi + } else if (getNDaug() == 4 && KPlusFound == 1 && PiPlusFound == 1 && PiMinusFound==1) + { + out_code = 7; // K+ pi+ pi- } else if (getNDaug() == 3 && theK.contains(getDaug(1)) && theK.contains(getDaug(2))) { out_code = 8; // KK - } else if (getNDaug() == 6 && - theK.contains(getDaug(1)) && theK.contains(getDaug(2)) - && thePis.contains(getDaug(3)) && thePis.contains(getDaug(4)) && thePis.contains(getDaug(5)) - ) + } else if (getNDaug() == 6 && KPlusFound == 1 && KMinusFound == 1 && PiPlusFound == 2 && PiMinusFound == 1) + { + out_code = 9; // K+ K- pi+ pi+ pi- + } else if (getNDaug() == 8 && PiPlusFound == 4 && PiMinusFound == 3) { - out_code = 9; // KK+3pi - } else if (getNDaug() == 8 && thePis.contains(getDaug(1)) && thePis.contains(getDaug(2)) && - thePis.contains(getDaug(3)) && thePis.contains(getDaug(4)) && thePis.contains(getDaug(5)) - && thePis.contains(getDaug(6)) && thePis.contains(getDaug(7)) ) { out_code = 10; //7pi - } else if (getNDaug() == 6 && theK.contains(getDaug(1)) && thePis.contains(getDaug(2)) && - thePis.contains(getDaug(3)) && thePis.contains(getDaug(4)) && thePis.contains(getDaug(5)) - ) + } else if (getNDaug() == 6 && KPlusFound == 1 && PiPlusFound == 2 && PiMinusFound == 2) { out_code = 11; // K+ pi+ pi+ pi- pi- } else { + report(ERROR,"EvtBcHad")<<"Init: unknown decay"<<std::endl; report(ERROR,"EvtBcHad")<<"Init: unknown decay"<<std::endl; }; report(INFO,"EvtBcHad")<<"out_code = "<<out_code<<", whichfit = "<<whichfit<<std::endl; + for(int i=0; i<4; ++i) + { + report(INFO, "i= " )<< i; + report(INFO, "iPiPlus= " )<< iPiPlus[i]; + report(INFO, "iPiMinus= ")<< iPiZero[i]; + report(INFO, "iPiZero= " )<< iPiMinus[i]; + report(INFO, "iKPlus= " )<< iKPlus[i]; + report(INFO, "iKMinus= " )<< iKMinus[i] << std::endl; + }; + report(INFO, "PiPlusFound= " ) << PiPlusFound; + report(INFO, "PiMinusFound= " ) << PiMinusFound; + report(INFO, "PiZeroFound= " ) << PiZeroFound; + report(INFO, "KPlusFound= " ) << KPlusFound; + report(INFO, "KMinusFound= " ) << KMinusFound << std::endl; + +} + + + +//====================================================== +void EvtBcVHad::init() { + + checkNArg(1); + checkSpinParent(EvtSpinType::SCALAR); + checkSpinDaughter(0, EvtSpinType::VECTOR); + for (int i = 1; i <= (getNDaug() - 1); i++) { + checkSpinDaughter(i, EvtSpinType::SCALAR); + } + + + idVector = getDaug(0).getId(); + whichfit = int(getArg(0) + 0.1); + ffmodel = new EvtBCVFF2(idVector, whichfit); + + wcurr = new EvtWHad(); + parseDecay(); + + // determine the code of final hadronic state + } //====================================================== @@ -212,43 +289,46 @@ EvtVector4C EvtBcVHad::hardCurr(EvtParticle *root_particle) const { if (out_code == 1) { // pi+ - hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4()); + hardCur = wcurr->WCurrent(root_particle->getDaug(iPiPlus[0])->getP4()); } else if (out_code == 2) { // pi+ pi0 - hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(), - root_particle->getDaug(2)->getP4()); + hardCur = wcurr->WCurrent( + root_particle->getDaug(iPiPlus[0])->getP4(), + root_particle->getDaug(iPiZero[0])->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()); - + hardCur = wcurr->WCurrent( + root_particle->getDaug(iPiPlus[0])->getP4(), + root_particle->getDaug(iPiPlus[1])->getP4(), + root_particle->getDaug(iPiMinus[0])->getP4() + ); } else if (out_code == 5) { - // Bc -> psi pi+ pi+ pi- pi- pi+ from Kuhn, Was, hep-ph/0602162 - hardCur = wcurr->WCurrent_5pi(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()); + hardCur = wcurr->WCurrent_5pi( + root_particle->getDaug(iPiPlus[0])->getP4(), + root_particle->getDaug(iPiPlus[1])->getP4(), + root_particle->getDaug(iPiPlus[2])->getP4(), + root_particle->getDaug(iPiMinus[0])->getP4(), + root_particle->getDaug(iPiMinus[1])->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()); + hardCur = wcurr->WCurrent_KKP( + root_particle->getDaug(iKPlus[0])->getP4(), + root_particle->getDaug(iKMinus[0])->getP4(), + root_particle->getDaug(iPiPlus[0])->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()); + hardCur = wcurr->WCurrent_KPP( + root_particle->getDaug(iKPlus[0])->getP4(), + root_particle->getDaug(iPiPlus[0])->getP4(), + root_particle->getDaug(iPiMinus[0])->getP4()); } else if (out_code == 8) { @@ -260,40 +340,38 @@ EvtVector4C EvtBcVHad::hardCurr(EvtParticle *root_particle) const { else if (out_code == 9) { // K+ K- pi+ pi+ pi- hardCur = wcurr -> WCurrent_KKPPP( - root_particle->getDaug(1)->getP4(), // K+ - root_particle->getDaug(2)->getP4(), // K- - root_particle->getDaug(3)->getP4(), // pi+ - root_particle->getDaug(4)->getP4(), // pi+ - root_particle->getDaug(5)->getP4() // pi- + root_particle->getDaug(iKPlus[0])->getP4(), // K+ + root_particle->getDaug(iKMinus[0])->getP4(), // K- + root_particle->getDaug(iPiPlus[0])->getP4(), // pi+ + root_particle->getDaug(iPiPlus[1])->getP4(), // pi+ + root_particle->getDaug(iPiMinus[0])->getP4() // pi- ); } else if (out_code == 10) { // 1=pi+ 2=pi+ 3=pi+ 4=pi+ 5=pi- 6=pi- 7=pi- with symmetrization of the identical particles hardCur = wcurr -> WCurrent_7pi( - root_particle->getDaug(1)->getP4(), // pi+ - root_particle->getDaug(2)->getP4(), // pi+ - root_particle->getDaug(3)->getP4(), // pi+ - root_particle->getDaug(4)->getP4(), // pi+ - root_particle->getDaug(5)->getP4(), // pi- - root_particle->getDaug(6)->getP4(), // pi- - root_particle->getDaug(7)->getP4() // pi- + root_particle->getDaug(iPiPlus[0])->getP4(), // pi+ + root_particle->getDaug(iPiPlus[1])->getP4(), // pi+ + root_particle->getDaug(iPiPlus[2])->getP4(), // pi+ + root_particle->getDaug(iPiPlus[3])->getP4(), // pi+ + root_particle->getDaug(iPiMinus[0])->getP4(), // pi- + root_particle->getDaug(iPiMinus[1])->getP4(), // pi- + root_particle->getDaug(iPiMinus[2])->getP4() // pi- ); } else if (out_code == 11) { // 1=K+ 2 = pi+ 3 = pi+ 4 = pi- 5 = pi- with symmetrizatiom hardCur = wcurr -> WCurrent_K4pi( - root_particle->getDaug(1)->getP4(), // pi+ - root_particle->getDaug(2)->getP4(), // pi+ - root_particle->getDaug(3)->getP4(), // pi+ - root_particle->getDaug(4)->getP4(), // pi- - root_particle->getDaug(5)->getP4() // pi- + root_particle->getDaug(iKPlus[0])->getP4(), // pi+ + root_particle->getDaug(iPiPlus[0])->getP4(), // pi+ + root_particle->getDaug(iPiPlus[1])->getP4(), // pi+ + root_particle->getDaug(iPiMinus[0])->getP4(), // pi- + root_particle->getDaug(iPiMinus[1])->getP4() // pi- ); } else { report(ERROR,"EvtBcHad")<<"hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; - // report(ERROR, "EvtGen") << "Have not yet implemented this final state in BC_VHAD model" << std::endl; - // report(ERROR, "EvtGen") << "Ndaug=" << getNDaug() << std::endl; for (int id = 0; id < (getNDaug() - 1); id++) { - // report(ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + //report(ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; } ::abort(); @@ -347,4 +425,3 @@ void EvtBcVHad::decay(EvtParticle *root_particle) { } } - -- GitLab From 3caf15aa7262365d1d410c00c1a51769e2cea70f Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Wed, 24 May 2023 21:01:17 +0200 Subject: [PATCH 11/14] Code polishing from John --- Gen/EvtGen/EvtGenModels/EvtBcVHad.hh | 20 ++-- Gen/EvtGen/src/EvtBcVHad.cpp | 153 ++++++++++++--------------- 2 files changed, 79 insertions(+), 94 deletions(-) diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh index 9f2780933..32179442e 100644 --- a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh @@ -23,6 +23,7 @@ #include "EvtGenBase/EvtVector4C.hh" #include "EvtGenBase/EvtDecayAmp.hh" +#include <array> #include <string> class EvtBCVFF2; @@ -59,26 +60,23 @@ private: // out_code: code of the hadronic final state // 1 - pi+ // 2 - pi+ pi0 - // 3 - pi+ pi+ pi- + // 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; - - + int out_code; EvtBCVFF2 *ffmodel; EvtWHad *wcurr; - int iPiPlus[4] = {-1, -1, -1, -1}, - iPiMinus[4] = {-1, -1, -1, -1}, - iPiZero[4] = {-1,-1,-1,-1}, - iKPlus[4] = {-1,-1,-1,-1}, - iKMinus[4] = {-1,-1,-1,-1}; - + std::array<int, 4> iPiPlus{-1, -1, -1, -1}; + std::array<int, 4> iPiMinus{-1, -1, -1, -1}; + std::array<int, 4> iPiZero{-1, -1, -1, -1}; + std::array<int, 4> iKPlus{-1, -1, -1, -1}; + std::array<int, 4> iKMinus{-1, -1, -1, -1}; }; -#endif \ No newline at end of file +#endif diff --git a/Gen/EvtGen/src/EvtBcVHad.cpp b/Gen/EvtGen/src/EvtBcVHad.cpp index 1ebb101de..64927eeab 100644 --- a/Gen/EvtGen/src/EvtBcVHad.cpp +++ b/Gen/EvtGen/src/EvtBcVHad.cpp @@ -74,125 +74,116 @@ EvtDecayBase* EvtBcVHad::clone() { void EvtBcVHad::parseDecay() { - EvtIdSet thePis("pi+", "pi-", "pi0"); - EvtIdSet theK("K+", "K-", "K_S0"); - EvtIdSet *PiPlusID, *PiMinusID, *PiZeroID, *KPlusID, *KMinusID, - *BcPlusID = new EvtIdSet("B_c+"), *BcMinusID = new EvtIdSet("B_c-"); - report(INFO,"getParentId()=")<<getParentId()<<std::endl; - - if(BcPlusID->contains(getParentId())) - { - PiPlusID = new EvtIdSet("pi+"); - PiMinusID = new EvtIdSet("pi-"); - PiZeroID = new EvtIdSet("pi0"); - KPlusID = new EvtIdSet("K+"); - KMinusID = new EvtIdSet("K-"); - } - else if(BcMinusID->contains(getParentId())) - { - PiPlusID = new EvtIdSet("pi-"); - PiMinusID = new EvtIdSet("pi+"); - PiZeroID = new EvtIdSet("pi0"); - KPlusID = new EvtIdSet("K-"); - KMinusID = new EvtIdSet("K+"); - }; - + const EvtIdSet BcPlusID("B_c+"), BcMinusID("B_c-"); + const EvtIdSet theK("K+", "K-", "K_S0"); + const int cMode = BcMinusID.contains(getParentId()); + const EvtIdSet PiPlusID(cMode == 0 ? "pi+" : "pi-"); + const EvtIdSet PiMinusID(cMode == 0 ? "pi-" : "pi+"); + const EvtIdSet PiZeroID("pi0"); + const EvtIdSet KPlusID(cMode == 0 ? "K+" : "K-"); + const EvtIdSet KMinusID(cMode == 0 ? "K-" : "K+"); + EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "parentId = " << getParentId() << std::endl; + int PiPlusFound = 0, PiMinusFound = 0, PiZeroFound = 0, KPlusFound = 0, KMinusFound = 0; - for(int iDaughter = 0; iDaughter<getNDaug(); ++iDaughter) + for (int iDaughter = 0; iDaughter < getNDaug(); ++iDaughter) { - report(INFO, "iDaughter=")<<iDaughter<<" id="<<getDaug(iDaughter).getName()<<std::endl; - if(PiPlusID->contains(getDaug(iDaughter))) + const EvtId daugId = getDaug(iDaughter); + EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "iDaughter = " << iDaughter + << " id = " <<getDaug(iDaughter).getName() <<std::endl; + if (PiPlusID.contains(daugId) && PiPlusFound < 4) { iPiPlus[PiPlusFound] = iDaughter; PiPlusFound++; } - else if(PiMinusID->contains(getDaug(iDaughter))) + else if (PiMinusID.contains(daugId) && PiMinusFound < 4) { iPiMinus[PiMinusFound] = iDaughter; PiMinusFound++; } - else if(PiZeroID->contains(getDaug(iDaughter))) + else if (PiZeroID.contains(daugId) && PiZeroFound < 4) { iPiZero[PiZeroFound] = iDaughter; PiZeroFound++; } - else if(KPlusID->contains(getDaug(iDaughter))) + else if (KPlusID.contains(daugId) && KPlusFound < 4) { iKPlus[KPlusFound] = iDaughter; KPlusFound++; } - else if(KMinusID->contains(getDaug(iDaughter))) + else if (KMinusID.contains(daugId) && KMinusFound < 4) { iKMinus[KMinusFound] = iDaughter; KMinusFound++; - }; - }; - + } + } - if (getNDaug() == 2 && PiPlusFound == 1) + if (getNDaug() == 2 && PiPlusFound == 1) { out_code = 1; // pi+ } - else if (getNDaug() == 3 && PiPlusFound == 1 && PiZeroFound==1) + else if (getNDaug() == 3 && PiPlusFound == 1 && PiZeroFound == 1) { out_code = 2; // pi+ pi0 } else if (getNDaug() == 4 && PiPlusFound == 2 && PiMinusFound == 1) { - out_code = 3;// pi+ pi+ pi- + out_code = 3; // pi+ pi+ pi- } - else if (getNDaug() == 5 && PiPlusFound == 2 && PiMinusFound==1 && PiZeroFound==1) + else if (getNDaug() == 5 && PiPlusFound == 2 && PiMinusFound == 1 && PiZeroFound == 1) { out_code = 4; // pi+ pi+ pi- pi0 } - else if (getNDaug() == 6 && PiPlusFound==3 && PiMinusFound == 2) + else if (getNDaug() == 6 && PiPlusFound == 3 && PiMinusFound == 2) { - out_code = 5; //5pi - } else if (getNDaug() == 4 && KPlusFound == 1 && KMinusFound==1 && PiPlusFound == 1) + out_code = 5; // 5pi + } + else if (getNDaug() == 4 && KPlusFound == 1 && KMinusFound == 1 && PiPlusFound == 1) { - out_code = 6; //KKpi - } else if (getNDaug() == 4 && KPlusFound == 1 && PiPlusFound == 1 && PiMinusFound==1) + out_code = 6; // KKpi + } + else if (getNDaug() == 4 && KPlusFound == 1 && PiPlusFound == 1 && PiMinusFound == 1) { out_code = 7; // K+ pi+ pi- - } else if (getNDaug() == 3 && theK.contains(getDaug(1)) && theK.contains(getDaug(2))) { + } + else if (getNDaug() == 3 && theK.contains(getDaug(1)) && theK.contains(getDaug(2))) + { out_code = 8; // KK - } else if (getNDaug() == 6 && KPlusFound == 1 && KMinusFound == 1 && PiPlusFound == 2 && PiMinusFound == 1) + } + else if (getNDaug() == 6 && KPlusFound == 1 && KMinusFound == 1 && PiPlusFound == 2 && PiMinusFound == 1) { out_code = 9; // K+ K- pi+ pi+ pi- - } else if (getNDaug() == 8 && PiPlusFound == 4 && PiMinusFound == 3) + } + else if (getNDaug() == 8 && PiPlusFound == 4 && PiMinusFound == 3) { - out_code = 10; //7pi - } else if (getNDaug() == 6 && KPlusFound == 1 && PiPlusFound == 2 && PiMinusFound == 2) + out_code = 10; // 7pi + } + else if (getNDaug() == 6 && KPlusFound == 1 && PiPlusFound == 2 && PiMinusFound == 2) { out_code = 11; // K+ pi+ pi+ pi- pi- - } + } else { - report(ERROR,"EvtBcHad")<<"Init: unknown decay"<<std::endl; - report(ERROR,"EvtBcHad")<<"Init: unknown decay"<<std::endl; - }; + EvtGenReport(EVTGEN_ERROR, "EvtBcVHad")<<"Init: unknown decay"<<std::endl; + } - report(INFO,"EvtBcHad")<<"out_code = "<<out_code<<", whichfit = "<<whichfit<<std::endl; - for(int i=0; i<4; ++i) + EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "out_code = " <<out_code << ", whichfit = " << whichfit << std::endl; + for (int i = 0; i < 4; ++i) { - report(INFO, "i= " )<< i; - report(INFO, "iPiPlus= " )<< iPiPlus[i]; - report(INFO, "iPiMinus= ")<< iPiZero[i]; - report(INFO, "iPiZero= " )<< iPiMinus[i]; - report(INFO, "iKPlus= " )<< iKPlus[i]; - report(INFO, "iKMinus= " )<< iKMinus[i] << std::endl; - }; - report(INFO, "PiPlusFound= " ) << PiPlusFound; - report(INFO, "PiMinusFound= " ) << PiMinusFound; - report(INFO, "PiZeroFound= " ) << PiZeroFound; - report(INFO, "KPlusFound= " ) << KPlusFound; - report(INFO, "KMinusFound= " ) << KMinusFound << std::endl; - + EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << " i = "<< i + << ", iPiPlus = " << iPiPlus[i] + << ", iPiMinus = " << iPiMinus[i] + << ", iPiZero = " << iPiZero[i] + << ", iKPlus = " << iKPlus[i] + << ", iKMinus = " << iKMinus[i] << std::endl; + } + EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "PiPlusFound = " << PiPlusFound + << ", PiMinusFound = " << PiMinusFound + << ", PiZeroFound = " << PiZeroFound + << ", KPlusFound = " << KPlusFound + << ", KMinusFound = " << KMinusFound << std::endl; } - - //====================================================== void EvtBcVHad::init() { @@ -203,16 +194,14 @@ void EvtBcVHad::init() { checkSpinDaughter(i, EvtSpinType::SCALAR); } - 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 parseDecay(); - // determine the code of final hadronic state - } //====================================================== @@ -264,18 +253,16 @@ void EvtBcVHad::initProbMax() { else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1) setProbMax(1.5e5); else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) setProbMax(1e5); } else if(out_code == 11) { - if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1) setProbMax(2.5e8); //setProbMax(1.4e7); + if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1) setProbMax(2.5e8); else if (idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2) setProbMax(1.4e7); else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1) setProbMax(2e6); else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) setProbMax(8e4); } else { - report(ERROR,"EvtBcHad")<<"probax: Have not yet implemented this final state in BC_VHAD model, out_code = " << out_code << std::endl; - // report(ERROR, "EvtGen") << "Have not yet implemented this final state in BC_VHAD model" << std::endl; - // report(ERROR, "EvtGen") << "Ndaug=" << getNDaug() << std::endl; - for (int id = 0; id < (getNDaug() - 1); id++) { - // report(ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; - } - ::abort(); + EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "probax: Have not yet implemented this final state in BC_VHAD model, out_code = " << out_code << std::endl; + for (int id = 0; id < (getNDaug() - 1); id++) { + EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + } + ::abort(); } } @@ -369,10 +356,10 @@ EvtVector4C EvtBcVHad::hardCurr(EvtParticle *root_particle) const { root_particle->getDaug(iPiMinus[1])->getP4() // pi- ); } else { - report(ERROR,"EvtBcHad")<<"hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; - for (int id = 0; id < (getNDaug() - 1); id++) { - //report(ERROR, "EvtGen") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; - } + EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; + for (int id = 0; id < (getNDaug() - 1); id++) { + EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + } ::abort(); } -- GitLab From 9734d91e2b6d7ef6051e011601688c5d6a4b4a37 Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Thu, 22 Jun 2023 21:04:42 +0200 Subject: [PATCH 12/14] Bug fixes --- Gen/EvtGen/EvtGenModels/EvtBcVHad.hh | 12 ++++++------ Gen/EvtGen/src/EvtBcVHad.cpp | 20 ++++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh index 32179442e..4d75d638c 100644 --- a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh @@ -53,7 +53,7 @@ private: // 1 - SR // 2 - PM int whichfit; - +S // idVector --- final vector particle code int idVector; @@ -71,11 +71,11 @@ private: EvtBCVFF2 *ffmodel; EvtWHad *wcurr; - std::array<int, 4> iPiPlus{-1, -1, -1, -1}; - std::array<int, 4> iPiMinus{-1, -1, -1, -1}; - std::array<int, 4> iPiZero{-1, -1, -1, -1}; - std::array<int, 4> iKPlus{-1, -1, -1, -1}; - std::array<int, 4> iKMinus{-1, -1, -1, -1}; + std::array<int, 4> iPiPlus ={-1, -1, -1, -1}; + std::array<int, 4> iPiMinus={-1, -1, -1, -1}; + std::array<int, 4> iPiZero ={-1, -1, -1, -1}; + std::array<int, 4> iKPlus ={-1, -1, -1, -1}; + std::array<int, 4> iKMinus ={-1, -1, -1, -1}; }; diff --git a/Gen/EvtGen/src/EvtBcVHad.cpp b/Gen/EvtGen/src/EvtBcVHad.cpp index 64927eeab..37f3430f1 100644 --- a/Gen/EvtGen/src/EvtBcVHad.cpp +++ b/Gen/EvtGen/src/EvtBcVHad.cpp @@ -83,13 +83,13 @@ void EvtBcVHad::parseDecay() const EvtIdSet PiZeroID("pi0"); const EvtIdSet KPlusID(cMode == 0 ? "K+" : "K-"); const EvtIdSet KMinusID(cMode == 0 ? "K-" : "K+"); - EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "parentId = " << getParentId() << std::endl; + report(INFO, "EvtBcVHad") << "parentId = " << getParentId() << std::endl; int PiPlusFound = 0, PiMinusFound = 0, PiZeroFound = 0, KPlusFound = 0, KMinusFound = 0; for (int iDaughter = 0; iDaughter < getNDaug(); ++iDaughter) { const EvtId daugId = getDaug(iDaughter); - EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "iDaughter = " << iDaughter + report(INFO, "EvtBcVHad") << "iDaughter = " << iDaughter << " id = " <<getDaug(iDaughter).getName() <<std::endl; if (PiPlusID.contains(daugId) && PiPlusFound < 4) { @@ -164,20 +164,20 @@ void EvtBcVHad::parseDecay() } else { - EvtGenReport(EVTGEN_ERROR, "EvtBcVHad")<<"Init: unknown decay"<<std::endl; + report(ERROR, "EvtBcVHad")<<"Init: unknown decay"<<std::endl; } - EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "out_code = " <<out_code << ", whichfit = " << whichfit << std::endl; + report(INFO, "EvtBcVHad") << "out_code = " <<out_code << ", whichfit = " << whichfit << std::endl; for (int i = 0; i < 4; ++i) { - EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << " i = "<< i + report(INFO, "EvtBcVHad") << " i = "<< i << ", iPiPlus = " << iPiPlus[i] << ", iPiMinus = " << iPiMinus[i] << ", iPiZero = " << iPiZero[i] << ", iKPlus = " << iKPlus[i] << ", iKMinus = " << iKMinus[i] << std::endl; } - EvtGenReport(EVTGEN_INFO, "EvtBcVHad") << "PiPlusFound = " << PiPlusFound + report(INFO, "EvtBcVHad") << "PiPlusFound = " << PiPlusFound << ", PiMinusFound = " << PiMinusFound << ", PiZeroFound = " << PiZeroFound << ", KPlusFound = " << KPlusFound @@ -258,9 +258,9 @@ void EvtBcVHad::initProbMax() { else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1) setProbMax(2e6); else if (idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2) setProbMax(8e4); } else { - EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "probax: Have not yet implemented this final state in BC_VHAD model, out_code = " << out_code << std::endl; + report(ERROR, "EvtBcVHad") << "probax: Have not yet implemented this final state in BC_VHAD model, out_code = " << out_code << std::endl; for (int id = 0; id < (getNDaug() - 1); id++) { - EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + report(ERROR, "EvtBcVHad") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; } ::abort(); } @@ -356,9 +356,9 @@ EvtVector4C EvtBcVHad::hardCurr(EvtParticle *root_particle) const { root_particle->getDaug(iPiMinus[1])->getP4() // pi- ); } else { - EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; + report(ERROR, "EvtBcVHad") << "hardCurr: Have not yet implemented this final state in BC_VHAD model" << std::endl; for (int id = 0; id < (getNDaug() - 1); id++) { - EvtGenReport(EVTGEN_ERROR, "EvtBcVHad") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; + report(ERROR, "EvtBcVHad") << "Daug " << id << " " << EvtPDL::name(getDaug(id)).c_str() << std::endl; } ::abort(); -- GitLab From 8a4d9f7587bb98c1c8b503e1f999ab19946b211e Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Thu, 22 Jun 2023 21:06:38 +0200 Subject: [PATCH 13/14] Bug fixes --- Gen/EvtGen/EvtGenModels/EvtBcVHad.hh | 1 - 1 file changed, 1 deletion(-) diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh index 4d75d638c..3a8403b38 100644 --- a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh @@ -53,7 +53,6 @@ private: // 1 - SR // 2 - PM int whichfit; -S // idVector --- final vector particle code int idVector; -- GitLab From b57dd856995685b49937f09ba17205576b76fa6f Mon Sep 17 00:00:00 2001 From: dpereima <dpereima dmitrii.pereima@cern.ch> Date: Fri, 23 Jun 2023 09:41:25 +0200 Subject: [PATCH 14/14] Bug fixes --- Gen/EvtGen/EvtGenModels/EvtBcVHad.hh | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh index 3a8403b38..a5bd55ab6 100644 --- a/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh +++ b/Gen/EvtGen/EvtGenModels/EvtBcVHad.hh @@ -70,11 +70,11 @@ private: EvtBCVFF2 *ffmodel; EvtWHad *wcurr; - std::array<int, 4> iPiPlus ={-1, -1, -1, -1}; - std::array<int, 4> iPiMinus={-1, -1, -1, -1}; - std::array<int, 4> iPiZero ={-1, -1, -1, -1}; - std::array<int, 4> iKPlus ={-1, -1, -1, -1}; - std::array<int, 4> iKMinus ={-1, -1, -1, -1}; + std::array<int, 4> iPiPlus ={{-1, -1, -1, -1}}; + std::array<int, 4> iPiMinus={{-1, -1, -1, -1}}; + std::array<int, 4> iPiZero ={{-1, -1, -1, -1}}; + std::array<int, 4> iKPlus ={{-1, -1, -1, -1}}; + std::array<int, 4> iKMinus ={{-1, -1, -1, -1}}; }; -- GitLab