diff --git a/History.md b/History.md
--- a/History.md
+++ b/History.md
@@ -11,6 +11,9 @@
 ===
 ## R02-02-0X
 
+16 Dec 2022 John Back
+* D89: Added probabilities for B_c -> V pi+ and V pi+ pi0 EvtBcVHad modes.
+
 16 Nov 2022 Tom Latham
 * T123: Provide documention of how to contribute bug reports, feature requests, new/modified code
 
diff --git a/src/EvtGenModels/EvtBcVHad.cpp b/src/EvtGenModels/EvtBcVHad.cpp
--- a/src/EvtGenModels/EvtBcVHad.cpp
+++ b/src/EvtGenModels/EvtBcVHad.cpp
@@ -126,7 +126,33 @@
 
 void EvtBcVHad::initProbMax()
 {
-    if ( out_code == 3 ) {
+    if ( out_code == 1 ) {
+        if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
+             getNDaug() == 2 )
+            setProbMax( 500. );
+        else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
+                  whichfit == 2 && getNDaug() == 2 )
+            setProbMax( 300. );
+        else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
+                  whichfit == 1 && getNDaug() == 2 )
+            setProbMax( 17. );
+        else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
+                  whichfit == 2 && getNDaug() == 2 )
+            setProbMax( 40. );
+    } else if ( out_code == 2 ) {
+        if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
+             getNDaug() == 3 )
+            setProbMax( 10950. );
+        else if ( idVector == EvtPDL::getId( "J/psi" ).getId() &&
+                  whichfit == 2 && getNDaug() == 3 )
+            setProbMax( 4200. );
+        else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
+                  whichfit == 1 && getNDaug() == 3 )
+            setProbMax( 500. );
+        else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
+                  whichfit == 2 && getNDaug() == 3 )
+            setProbMax( 700. );
+    } else if ( out_code == 3 ) {
         if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
              getNDaug() == 4 )
             setProbMax( 42000. );
@@ -139,7 +165,6 @@
         else if ( idVector == EvtPDL::getId( "psi(2S)" ).getId() &&
                   whichfit == 2 && getNDaug() == 4 )
             setProbMax( 2600. );
-
     } else if ( out_code == 5 ) {
         if ( idVector == EvtPDL::getId( "J/psi" ).getId() && whichfit == 1 &&
              getNDaug() == 6 )
diff --git a/test/jsonFiles/BC_VHAD=VLL__Bc_Jpsipipi0_mumu.json b/test/jsonFiles/BC_VHAD=VLL__Bc_Jpsipipi0_mumu.json
new file mode 100644
--- /dev/null
+++ b/test/jsonFiles/BC_VHAD=VLL__Bc_Jpsipipi0_mumu.json
@@ -0,0 +1,38 @@
+{
+    "parent" : "B_c+",
+    "daughters" : ["J/psi", "pi+", "pi0"],
+    "grand_daughters" : [["mu+", "mu-"], [], []],
+    "models" : ["BC_VHAD", "VLL", "", ""],
+    "parameters" : [["1"], [], [], []],
+    "extras" : ["noPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "Prob(BC_VHAD)",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "prob_daug1", "title" : "Prob(VLL)",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "mass", "title" : "m(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.0, "xmax" : 3.2},
+        {"variable" : "mass", "title" : "m(Jpsi pi+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0, "xmax" : 6.5},
+        {"variable" : "mass", "title" : "m(Jpsi pi0)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 3.0, "xmax" : 6.5},
+        {"variable" : "mass", "title" : "m(pi+ pi0)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 3.0},
+        {"variable" : "pSumSq", "title" : "qSq(pi+ pi0)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 4.0},
+        {"variable" : "p", "title" : "p(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.5, "xmax" : 2.5},
+        {"variable" : "pz", "title" : "pz(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -3.0, "xmax" : 3.0},
+        {"variable" : "cosHel", "title" : "cosHel(Jpsi,pi+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel", "title" : "cosHel(pi+,pi0)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel", "title" : "cosHel(Jpsi,pi0)", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel_daug1", "title" : "cosHel(mu+,Jpsi)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel_daug1", "title" : "cosHel(mu-,Jpsi)", "d1" : 2, "d2" : 1, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosAcoplanarityAngle", "title" : "cosAcopl", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "acoplanarityAngle", "title" : "acopl", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0, "xmax" : 180.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "E", "title" : "E(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.5, "xmax" : 4.2},
+	{"variable" : "phi", "title" : "phi(Jpsi)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+	{"variable" : "pLab_daug1", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
+	{"variable" : "pLab_daug1", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
+	{"variable" : "cosTheta_daug1", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "phi_daug1", "title" : "phi(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+	{"variable" : "cosTheta_daug1", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "phi_daug1", "title" : "phi(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
+    ],
+    "outfile" : "BC_VHAD=VLL__Bc_Jpsipipi0_mumu.root",
+    "reference" : "RefBC_VHAD=VLL__Bc_JpsiKpipi0_mumu.root"
+}
diff --git a/test/jsonFiles/BC_VHAD=VLL__Bc_psi2Spi_mumu.json b/test/jsonFiles/BC_VHAD=VLL__Bc_psi2Spi_mumu.json
new file mode 100644
--- /dev/null
+++ b/test/jsonFiles/BC_VHAD=VLL__Bc_psi2Spi_mumu.json
@@ -0,0 +1,29 @@
+{
+    "parent" : "B_c+",
+    "daughters" : ["psi(2S)", "pi+"],
+    "grand_daughters" : [["mu+", "mu-"], []],
+    "models" : ["BC_VHAD", "VLL", ""],
+    "parameters" : [["2"], [], []],
+    "extras" : ["noPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "Prob(BC_VHAD)",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "prob_daug1", "title" : "Prob(VLL)",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "mass", "title" : "m(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 3.68, "xmax" : 3.695},
+        {"variable" : "mass", "title" : "m(psi2S pi+)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 6.0, "xmax" : 6.5},
+        {"variable" : "p", "title" : "p(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 2.0, "xmax" : 2.1},
+        {"variable" : "pz", "title" : "pz(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -2.5, "xmax" : 2.5},
+        {"variable" : "cosHel", "title" : "cosHel(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "E", "title" : "E(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 4.0, "xmax" : 4.5},
+	{"variable" : "phi", "title" : "phi(psi2S)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+	{"variable" : "pLab_daug1", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
+	{"variable" : "pLab_daug1", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.5, "xmax" : 3.5},
+	{"variable" : "cosTheta_daug1", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "phi_daug1", "title" : "phi(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+	{"variable" : "cosTheta_daug1", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "phi_daug1", "title" : "phi(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
+    ],
+    "outfile" : "BC_VHAD=VLL__Bc_psi2Spi_mumu.root",
+    "reference" : "RefBC_VHAD=VLL__Bc_psi2Spi_mumu.root"
+}