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