From 5c7d701300f1a7c2d860cceaad1070aa9ee8ab8c Mon Sep 17 00:00:00 2001
From: Ulascan Sarica <ulascan90@gmail.com>
Date: Sun, 18 Dec 2016 18:11:44 +0100
Subject: [PATCH 1/4] WW bugfix and JJQCD missing elements

---
 MELA/interface/TMCFMUtils.hh |  2 ++
 MELA/src/TMCFMUtils.cc       | 67 ++++++++++++++++++++++++++++++++++--
 MELA/test/loadMELA.C         |  5 +--
 MELA/test/testME_v2.c        | 22 ++++++------
 4 files changed, 78 insertions(+), 18 deletions(-)

diff --git a/MELA/interface/TMCFMUtils.hh b/MELA/interface/TMCFMUtils.hh
index 936657cb..fbc5a228 100755
--- a/MELA/interface/TMCFMUtils.hh
+++ b/MELA/interface/TMCFMUtils.hh
@@ -33,6 +33,8 @@ namespace TMCFMUtils{
 
   void AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rSel, int sSel, int order[2]);
   std::vector<intQuad_t> Hash_QQVVQQAny();
+  std::vector<intQuad_t> Hash_QQVVQQ();
+  std::vector<intQuad_t> Hash_QQVVQQStrong();
 }
 
 #endif
diff --git a/MELA/src/TMCFMUtils.cc b/MELA/src/TMCFMUtils.cc
index 74254abc..1521a5cc 100755
--- a/MELA/src/TMCFMUtils.cc
+++ b/MELA/src/TMCFMUtils.cc
@@ -86,6 +86,14 @@ void TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rS
   if (!outFound){ for (unsigned int ip=0; ip<2; ip++) order[ip]=-1; }
 }
 std::vector<TMCFMUtils::intQuad_t> TMCFMUtils::Hash_QQVVQQAny(){
+  std::vector<TMCFMUtils::intQuad_t> pcfg;
+  std::vector<TMCFMUtils::intQuad_t> hash_qqvvqq = TMCFMUtils::Hash_QQVVQQ();
+  std::vector<TMCFMUtils::intQuad_t> hash_qqvvqqstrong = TMCFMUtils::Hash_QQVVQQStrong();
+  for (unsigned int c=0; c<hash_qqvvqq.size(); c++) pcfg.push_back(hash_qqvvqq.at(c));
+  for (unsigned int c=0; c<hash_qqvvqqstrong.size(); c++) pcfg.push_back(hash_qqvvqqstrong.at(c));
+  return pcfg;
+}
+std::vector<TMCFMUtils::intQuad_t> TMCFMUtils::Hash_QQVVQQ(){
   // FIXME: NEED TO INCLUDE LEPTONS AS WELL FOR J=4-8 AND 9-12
   /*
   Based on the following cases in MCFM:
@@ -163,9 +171,62 @@ std::vector<TMCFMUtils::intQuad_t> TMCFMUtils::Hash_QQVVQQAny(){
   // Uncommenting the lines below prints out the hash when the library is loaded.
   /*
   for (unsigned int ic=0; ic<pcfg.size(); ic++) std::cout
-    << "TMCFMUtils::Hash_QQVVQQAny: Hash configuration " << ic << " requests ids=( "
-    << pcfg.at(ic)[0] << ", " << pcfg.at(ic)[1] << ", " << pcfg.at(ic)[2] << ", " << pcfg.at(ic)[3]
-    << ")" << std::endl;
+  << "TMCFMUtils::Hash_QQVVQQAny: Hash configuration " << ic << " requests ids=( "
+  << pcfg.at(ic)[0] << ", " << pcfg.at(ic)[1] << ", " << pcfg.at(ic)[2] << ", " << pcfg.at(ic)[3]
+  << ")" << std::endl;
+  */
+  return pcfg;
+}
+std::vector<TMCFMUtils::intQuad_t> TMCFMUtils::Hash_QQVVQQStrong(){
+  /*
+  Based on the following cases in MCFM:
+  call qq4lggampf(1,2,3,4,5,6,7,8,3,4,za,zb,msqgg)
+  msq(1,-1)=msq(1,-1)+stat*aveqq*msqgg(1)
+  call qq4lggampf(2,1,3,4,5,6,7,8,3,4,za,zb,msqgg)
+  msq(-1,1)=msq(-1,1)+stat*aveqq*msqgg(1)
+  call qq4lggampf(7,8,3,4,5,6,1,2,3,4,za,zb,msqgg)
+  msq(0,0)=msq(0,0)+avegg*(3d0*msqgg(1)+2d0*msqgg(2))
+  call qq4lggampf(7,2,3,4,5,6,1,8,3,4,za,zb,msqgg)
+  msq(0,-1)=msq(0,-1)+aveqg*msqgg(1)
+  call qq4lggampf(7,1,3,4,5,6,2,8,3,4,za,zb,msqgg)
+  msq(-1,0)=msq(-1,0)+aveqg*msqgg(1)
+  call qq4lggampf(2,8,3,4,5,6,1,7,3,4,za,zb,msqgg)
+  msq(0,1)=msq(0,1)+aveqg*msqgg(1)
+  call qq4lggampf(1,8,3,4,5,6,2,7,3,4,za,zb,msqgg)
+  msq(1,0)=msq(1,0)+aveqg*msqgg(1)
+  */
+  std::vector<TMCFMUtils::intQuad_t> base_cfg;
+  // Start with qqb_gg
+  for (int iq=1; iq<=5; iq++) base_cfg.push_back(TMCFMUtils::intQuad_t(iq, -iq, 21, 21));
+
+  std::vector<TMCFMUtils::intQuad_t> jcfg;
+  jcfg.push_back(TMCFMUtils::intQuad_t(0, 1, 2, 3)); // qqb->gg
+  jcfg.push_back(TMCFMUtils::intQuad_t(1, 0, 2, 3)); // qbq->gg
+  jcfg.push_back(TMCFMUtils::intQuad_t(2, 3, 0, 1)); // gg->qbq
+  jcfg.push_back(TMCFMUtils::intQuad_t(2, 1, 0, 3)); // gqb->qbg
+  jcfg.push_back(TMCFMUtils::intQuad_t(1, 2, 0, 3)); // qbg->qbg
+  jcfg.push_back(TMCFMUtils::intQuad_t(3, 0, 2, 1)); // gq->gq
+  jcfg.push_back(TMCFMUtils::intQuad_t(0, 3, 2, 1)); // qg->gq
+
+  std::vector<TMCFMUtils::intQuad_t> pcfg;
+  for (unsigned int j=0; j<jcfg.size(); j++){
+    for (unsigned int p=0; p<base_cfg.size(); p++){
+      TMCFMUtils::intQuad_t cfg;
+      for (unsigned int ipos=0; ipos<4; ipos++){
+        int idpos = jcfg.at(j)[ipos];
+        int idAssigned = base_cfg.at(p)[ipos];
+        if (((idpos<2 && ipos>=2) || (idpos>=2 && ipos<2)) && !PDGHelpers::isAGluon(idAssigned)) idAssigned = -idAssigned;
+        cfg[jcfg.at(j)[ipos]] = idAssigned;
+      }
+      pcfg.push_back(cfg);
+    }
+  }
+  // Uncommenting the lines below prints out the hash when the library is loaded.
+  /*
+  for (unsigned int ic=0; ic<pcfg.size(); ic++) std::cout
+  << "TMCFMUtils::Hash_QQVVQQAny: Hash configuration " << ic << " requests ids=( "
+  << pcfg.at(ic)[0] << ", " << pcfg.at(ic)[1] << ", " << pcfg.at(ic)[2] << ", " << pcfg.at(ic)[3]
+  << ")" << std::endl;
   */
   return pcfg;
 }
diff --git a/MELA/test/loadMELA.C b/MELA/test/loadMELA.C
index 7d356fc7..a12f36d1 100755
--- a/MELA/test/loadMELA.C
+++ b/MELA/test/loadMELA.C
@@ -1,8 +1,5 @@
-//
-// MELA root package loader - see testKD.C for instructions
-//
 {
-  TString LIBMCFM="libmcfm_702.so";
+  TString LIBMCFM="libmcfm_702b2.so";
   gSystem->AddIncludePath("-I$ROOFITSYS/include/");
   gSystem->AddIncludePath("-I$CMSSW_BASE/interface/");
   gSystem->AddIncludePath("-I$CMSSW_BASE/src/");
diff --git a/MELA/test/testME_v2.c b/MELA/test/testME_v2.c
index 5ed9b19a..d428a5e4 100755
--- a/MELA/test/testME_v2.c
+++ b/MELA/test/testME_v2.c
@@ -1294,16 +1294,16 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
 
   int idMother[2]={ 0 };
   // VBF ZZ(+)WW
-  if (motherflavor==1){ idMother[0]=2; idMother[1]=1; }
-  else if (motherflavor==2){ idMother[0]=-2; idMother[1]=-1; }
+  if (motherflavor==1){ idMother[0]=2; idMother[1]=1; } // Passed: (1,1,0). (1,2,0) SM ok, fa3 not (fa3 ratio less than SM ratio, no swap seen).
+  else if (motherflavor==2){ idMother[0]=-2; idMother[1]=-1; } // Passed: (2,1,0). (2,2,0) SM ok, fa3 not (fa3 ratio less than SM ratio, no swap seen).
   // VBF ZZ-only(+)WH
   else if (motherflavor==3){ idMother[0]=2; idMother[1]=-1; } // Passed: (3,1,0);(3,2,1) -> (3,0,0);(3,0,1) -> Check!
-  else if (motherflavor==4){ idMother[0]=-2; idMother[1]=1; }
+  else if (motherflavor==4){ idMother[0]=-2; idMother[1]=1; } // Passed: (4,1,0);(4,2,1) -> (4,0,0);(4,0,1) -> Check!
   // VBF ZZ(+)ZH or WW(+)ZH
-  else if (motherflavor==5){ idMother[0]=2; idMother[1]=-2; } // Passed: (5,1,0)
-  else if (motherflavor==6){ idMother[0]=-2; idMother[1]=2; }
-  else if (motherflavor==7){ idMother[0]=1; idMother[1]=-1; }
-  else if (motherflavor==8){ idMother[0]=-1; idMother[1]=1; }
+  else if (motherflavor==5){ idMother[0]=2; idMother[1]=-2; } // Passed: (5,1,0). (5,2,0) SM ok, fa3 not (fa3 ratio less than SM ratio, no swap seen, same as 1,2,0). (5,1,1) does not give correct ratios, qqb and qbq initial states have opposite uub, ddb less-more trend.
+  else if (motherflavor==6){ idMother[0]=-2; idMother[1]=2; } // Passed: (6,1,0)
+  else if (motherflavor==7){ idMother[0]=1; idMother[1]=-1; } // Passed: (7,1,0)
+  else if (motherflavor==8){ idMother[0]=-1; idMother[1]=1; } // Passed: (8,1,0)
 
   bool doEval=true;
   TVar::VerbosityLevel verbosity = TVar::DEBUG;
@@ -1467,9 +1467,6 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
     double invckm_ud=1, invckm_us=0, invckm_cd=0, invckm_cs=1, invckm_ts=0, invckm_tb=1, invckm_ub=0, invckm_cb=0, invckm_td=0;
     TUtil::SetCKMElements(&invckm_ud, &invckm_us, &invckm_cd, &invckm_cs, &invckm_ts, &invckm_tb, &invckm_ub, &invckm_cb, &invckm_td);
 
-    if (isZZWW==2) spinzerohiggs_anomcoupl_.AnomalCouplDK=0;
-    else spinzerohiggs_anomcoupl_.AnomalCouplDK=1;
-
     /***** JHUGen *****/
 
     mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, TVar::ZZINDEPENDENT);
@@ -1538,7 +1535,7 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
     p_prod_0minus_VAJHU *= propagator;
 
     p_prod_0mplus_dec_0mplus_VAJHU=p_prod_0mplus_VAJHU*p_dec_0mplus_VAJHU; for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) pArray_prod_0mplus_dec_0mplus_VAJHU[ii][jj] = pArray_prod_0mplus_VAJHU[ii][jj]*p_dec_0mplus_VAJHU; }
-    if (prod==TVar::Had_WH_S && isZZWW==2){ // MCFM setting above turns off ZZ anomalous couplings if WH with WW couplings tested, so we should do the same here. 
+    if (isZZWW==2){ // MCFM setting above turns off ZZ anomalous couplings if WH with WW couplings tested, so we should do the same here. 
       p_prod_0minus_dec_0minus_VAJHU=p_prod_0minus_VAJHU*p_dec_0mplus_VAJHU; for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) pArray_prod_0minus_dec_0minus_VAJHU[ii][jj] = pArray_prod_0minus_VAJHU[ii][jj]*p_dec_0mplus_VAJHU; }
     }
     else{
@@ -1547,6 +1544,9 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
 
     /***** MCFM *****/
 
+    if (isZZWW==2) spinzerohiggs_anomcoupl_.AnomalCouplDK=0;
+    else spinzerohiggs_anomcoupl_.AnomalCouplDK=1;
+
     mela.setProcess(TVar::HSMHiggs, TVar::MCFM, prod);
 
     if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }

From 87ee1eda68e6863280232a037f0295fe3fef3bf5 Mon Sep 17 00:00:00 2001
From: Ulascan Sarica <ulascan90@gmail.com>
Date: Tue, 20 Dec 2016 22:02:04 +0100
Subject: [PATCH 2/4] testME_v2.c

---
 MELA/test/testME_v2.c | 196 ++++++++++++++++++++++++++++++------------
 1 file changed, 139 insertions(+), 57 deletions(-)

diff --git a/MELA/test/testME_v2.c b/MELA/test/testME_v2.c
index d428a5e4..f431bdd8 100755
--- a/MELA/test/testME_v2.c
+++ b/MELA/test/testME_v2.c
@@ -1283,27 +1283,60 @@ void testME_ProdDec_MCFM_Ordering(int iSel, int jSel, int rSel, int sSel){
     << order[1] << endl;
 }
 
-void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0 /*1==ZZ, 2==WW*/, int vbfvhchannel=0 /*0==VBF, 1==VH*/){
+void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0 /*1==ZZ, 2==WW*/, int vbfvhchannel=0 /*0==VBF, 1==VH*/, int hasInterf=0){
   ofstream tout("testME_ProdDec_MCFM_JHUGen_Comparison_Ping.out");
   streambuf* coutbuf = cout.rdbuf();
   //cout.rdbuf(tout.rdbuf());
 
+  struct mcfmme{
+    float proddecme;
+    double mearray[nmsq][nmsq];
+    mcfmme(){
+      proddecme=0;
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) mearray[ii][jj]=0; }
+    }
+    void multiplyarray(const float val){
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) mearray[ii][jj]*=val; }
+    }
+    void printarray(){
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << mearray[ii][jj]; cout << endl; }
+    }
+  };
+  struct jhume{
+    float prodme;
+    float decme;
+    float proddecme;
+    double mearray[nmsq][nmsq];
+    jhume(){
+      prodme=0;
+      decme=0;
+      proddecme=0;
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) mearray[ii][jj]=0; }
+    }
+    void multiplyarray(const float val){
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) mearray[ii][jj]*=val; }
+    }
+    void printarray(){
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << mearray[ii][jj]; cout << endl; }
+    }
+  };
+
   int erg_tev=13;
   float mPOLE=125.0;
   float wPOLE=4.07e-3;
 
   int idMother[2]={ 0 };
   // VBF ZZ(+)WW
-  if (motherflavor==1){ idMother[0]=2; idMother[1]=1; } // Passed: (1,1,0). (1,2,0) SM ok, fa3 not (fa3 ratio less than SM ratio, no swap seen).
-  else if (motherflavor==2){ idMother[0]=-2; idMother[1]=-1; } // Passed: (2,1,0). (2,2,0) SM ok, fa3 not (fa3 ratio less than SM ratio, no swap seen).
+  if (motherflavor==1){ idMother[0]=2; idMother[1]=1; } // Passed: (1,1,0);(1,2,0) -> Check!
+  else if (motherflavor==2){ idMother[0]=-2; idMother[1]=-1; } // Passed: (2,1,0);(2,2,0) -> Check!
   // VBF ZZ-only(+)WH
   else if (motherflavor==3){ idMother[0]=2; idMother[1]=-1; } // Passed: (3,1,0);(3,2,1) -> (3,0,0);(3,0,1) -> Check!
   else if (motherflavor==4){ idMother[0]=-2; idMother[1]=1; } // Passed: (4,1,0);(4,2,1) -> (4,0,0);(4,0,1) -> Check!
   // VBF ZZ(+)ZH or WW(+)ZH
-  else if (motherflavor==5){ idMother[0]=2; idMother[1]=-2; } // Passed: (5,1,0). (5,2,0) SM ok, fa3 not (fa3 ratio less than SM ratio, no swap seen, same as 1,2,0). (5,1,1) does not give correct ratios, qqb and qbq initial states have opposite uub, ddb less-more trend.
-  else if (motherflavor==6){ idMother[0]=-2; idMother[1]=2; } // Passed: (6,1,0)
-  else if (motherflavor==7){ idMother[0]=1; idMother[1]=-1; } // Passed: (7,1,0)
-  else if (motherflavor==8){ idMother[0]=-1; idMother[1]=1; } // Passed: (8,1,0)
+  else if (motherflavor==5){ idMother[0]=2; idMother[1]=-2; } // Passed: (5,1,0);(5,2,0);(5,1,1) -> Check
+  else if (motherflavor==6){ idMother[0]=-2; idMother[1]=2; } // Passed: (6,1,0);(6,2,0);(6,1,1) -> Check
+  else if (motherflavor==7){ idMother[0]=1; idMother[1]=-1; } // Passed: (7,1,0);(7,2,0);(7,1,1) -> Check
+  else if (motherflavor==8){ idMother[0]=-1; idMother[1]=1; } // Passed: (8,1,0);(8,2,0);(8,1,1) -> Check
 
   bool doEval=true;
   TVar::VerbosityLevel verbosity = TVar::DEBUG;
@@ -1319,10 +1352,13 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
   if (doEval){
     Mela mela(erg_tev, mPOLE, verbosity);
 
-    float p_prod_0mplus_dec_0mplus_VAJHU; double pArray_prod_0mplus_dec_0mplus_VAJHU[nmsq][nmsq]={ 0 };
-    float p_prod_0mplus_dec_0mplus_VAMCFM; double pArray_prod_0mplus_dec_0mplus_VAMCFM[nmsq][nmsq]={ 0 };
-    float p_prod_0minus_dec_0minus_VAJHU; double pArray_prod_0minus_dec_0minus_VAJHU[nmsq][nmsq]={ 0 };
-    float p_prod_0minus_dec_0minus_VAMCFM; double pArray_prod_0minus_dec_0minus_VAMCFM[nmsq][nmsq]={ 0 };
+    jhume p_prod_0mplus_dec_0mplus_VAJHU;
+    jhume p_prod_0minus_dec_0minus_VAJHU;
+    jhume p_prod_fa3_dec_fa3_VAJHU;
+
+    mcfmme p_prod_0mplus_dec_0mplus_VAMCFM;
+    mcfmme p_prod_0minus_dec_0minus_VAMCFM;
+    mcfmme p_prod_fa3_dec_fa3_VAMCFM;
 
     float mzz = 0;
     float mjj = 0;
@@ -1351,8 +1387,14 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
     int GenLep1Id=0, GenLep2Id=0, GenLep3Id=0, GenLep4Id=0;
     GenLep1Id=13;
     GenLep2Id=-13;
-    GenLep3Id=11;
-    GenLep4Id=-11;
+    if (hasInterf==0){
+      GenLep3Id=11;
+      GenLep4Id=-11;
+    }
+    else{
+      GenLep3Id=GenLep1Id;
+      GenLep4Id=GenLep2Id;
+    }
     mela.setCandidateDecayMode(TVar::CandidateDecay_ZZ);
     int idOrdered[4] ={ GenLep1Id, GenLep2Id, GenLep3Id, GenLep4Id };
 
@@ -1468,32 +1510,35 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
     TUtil::SetCKMElements(&invckm_ud, &invckm_us, &invckm_cd, &invckm_cs, &invckm_ts, &invckm_tb, &invckm_ub, &invckm_cb, &invckm_td);
 
     /***** JHUGen *****/
-
     mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, TVar::ZZINDEPENDENT);
-    float p_dec_0mplus_VAJHU;
+
     mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1;
-    mela.computeP(p_dec_0mplus_VAJHU, false);
+    mela.computeP(p_prod_0mplus_dec_0mplus_VAJHU.decme, false);
 
-    float p_dec_0minus_VAJHU;
     mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1;
-    mela.computeP(p_dec_0minus_VAJHU, false);
+    mela.computeP(p_prod_0minus_dec_0minus_VAJHU.decme, false);
+
+    mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1;
+    mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1;
+    mela.computeP(p_prod_fa3_dec_fa3_VAJHU.decme, false);
 
-    float p_prod_0mplus_VAJHU=0; double pArray_prod_0mplus_VAJHU[nmsq][nmsq]={ 0 };
-    float p_prod_0minus_VAJHU=0; double pArray_prod_0minus_VAJHU[nmsq][nmsq]={ 0 };
     if (vbfvhchannel==0){
       mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, TVar::JJVBF);
 
-      mela.setVerbosity(TVar::DEBUG_VERBOSE);
       if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }
       else{ mela.selfDHwwcoupl[0][gHIGGS_VV_1][0]=1; mela.differentiate_HWW_HZZ=true; }
-      mela.computeProdP(p_prod_0mplus_VAJHU, false);
-      mela.getIORecord()->getUnweightedMEArray(pArray_prod_0mplus_VAJHU);
-      mela.setVerbosity(TVar::DEBUG);
+      mela.computeProdP(p_prod_0mplus_dec_0mplus_VAJHU.prodme, false);
+      mela.getIORecord()->getUnweightedMEArray(p_prod_0mplus_dec_0mplus_VAJHU.mearray);
 
       if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }
       else{ mela.selfDHwwcoupl[0][gHIGGS_VV_4][0]=1; mela.differentiate_HWW_HZZ=true; }
-      mela.computeProdP(p_prod_0minus_VAJHU, false);
-      mela.getIORecord()->getUnweightedMEArray(pArray_prod_0minus_VAJHU);
+      mela.computeProdP(p_prod_0minus_dec_0minus_VAJHU.prodme, false);
+      mela.getIORecord()->getUnweightedMEArray(p_prod_0minus_dec_0minus_VAJHU.mearray);
+
+      if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }
+      else{ mela.selfDHwwcoupl[0][gHIGGS_VV_1][0]=1; mela.selfDHwwcoupl[0][gHIGGS_VV_4][0]=1; mela.differentiate_HWW_HZZ=true; }
+      mela.computeProdP(p_prod_fa3_dec_fa3_VAJHU.prodme, false);
+      mela.getIORecord()->getUnweightedMEArray(p_prod_fa3_dec_fa3_VAJHU.mearray);
     }
     else{
       if (prod==TVar::Had_ZH_S){
@@ -1501,26 +1546,36 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
 
         if (isZZWW!=2) mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1;
         else mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=0;
-        mela.computeProdP(p_prod_0mplus_VAJHU, false);
-        mela.getIORecord()->getUnweightedMEArray(pArray_prod_0mplus_VAJHU);
+        mela.computeProdP(p_prod_0mplus_dec_0mplus_VAJHU.prodme, false);
+        mela.getIORecord()->getUnweightedMEArray(p_prod_0mplus_dec_0mplus_VAJHU.mearray);
 
         if (isZZWW!=2) mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1;
         else mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=0;
-        mela.computeProdP(p_prod_0minus_VAJHU, false);
-        mela.getIORecord()->getUnweightedMEArray(pArray_prod_0minus_VAJHU);
+        mela.computeProdP(p_prod_0minus_dec_0minus_VAJHU.prodme, false);
+        mela.getIORecord()->getUnweightedMEArray(p_prod_0minus_dec_0minus_VAJHU.mearray);
+
+        if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; }
+        else{ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=0; mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=0; }
+        mela.computeProdP(p_prod_fa3_dec_fa3_VAJHU.prodme, false);
+        mela.getIORecord()->getUnweightedMEArray(p_prod_fa3_dec_fa3_VAJHU.mearray);
       }
       else if (prod==TVar::Had_WH_S){
         mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, TVar::Had_WH);
 
         if (isZZWW!=1) mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1;
         else mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=0;
-        mela.computeProdP(p_prod_0mplus_VAJHU, false);
-        mela.getIORecord()->getUnweightedMEArray(pArray_prod_0mplus_VAJHU);
+        mela.computeProdP(p_prod_0mplus_dec_0mplus_VAJHU.prodme, false);
+        mela.getIORecord()->getUnweightedMEArray(p_prod_0mplus_dec_0mplus_VAJHU.mearray);
 
         if (isZZWW!=1) mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1;
         else mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=0;
-        mela.computeProdP(p_prod_0minus_VAJHU, false);
-        mela.getIORecord()->getUnweightedMEArray(pArray_prod_0minus_VAJHU);
+        mela.computeProdP(p_prod_0minus_dec_0minus_VAJHU.prodme, false);
+        mela.getIORecord()->getUnweightedMEArray(p_prod_0minus_dec_0minus_VAJHU.mearray);
+
+        if (isZZWW!=1){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; }
+        else{ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=0; mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=0; }
+        mela.computeProdP(p_prod_fa3_dec_fa3_VAJHU.prodme, false);
+        mela.getIORecord()->getUnweightedMEArray(p_prod_fa3_dec_fa3_VAJHU.mearray);
       }
     }
 
@@ -1531,15 +1586,18 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
       mela.getIORecord()->getHiggsMassWidth(mh, gah, 0);
       propagator /= 1./(pow(pow(mzz, 2)-pow(mh, 2), 2) + pow(mh*gah, 2));
     }
-    p_prod_0mplus_VAJHU *= propagator;
-    p_prod_0minus_VAJHU *= propagator;
+    p_prod_0mplus_dec_0mplus_VAJHU.prodme *= propagator; p_prod_0mplus_dec_0mplus_VAJHU.multiplyarray(propagator);
+    p_prod_0minus_dec_0minus_VAJHU.prodme *= propagator; p_prod_0minus_dec_0minus_VAJHU.multiplyarray(propagator);
+    p_prod_fa3_dec_fa3_VAJHU.prodme *= propagator; p_prod_fa3_dec_fa3_VAJHU.multiplyarray(propagator);
 
-    p_prod_0mplus_dec_0mplus_VAJHU=p_prod_0mplus_VAJHU*p_dec_0mplus_VAJHU; for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) pArray_prod_0mplus_dec_0mplus_VAJHU[ii][jj] = pArray_prod_0mplus_VAJHU[ii][jj]*p_dec_0mplus_VAJHU; }
+    p_prod_0mplus_dec_0mplus_VAJHU.proddecme = p_prod_0mplus_dec_0mplus_VAJHU.prodme*p_prod_0mplus_dec_0mplus_VAJHU.decme; p_prod_0mplus_dec_0mplus_VAJHU.multiplyarray(p_prod_0mplus_dec_0mplus_VAJHU.decme);
     if (isZZWW==2){ // MCFM setting above turns off ZZ anomalous couplings if WH with WW couplings tested, so we should do the same here. 
-      p_prod_0minus_dec_0minus_VAJHU=p_prod_0minus_VAJHU*p_dec_0mplus_VAJHU; for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) pArray_prod_0minus_dec_0minus_VAJHU[ii][jj] = pArray_prod_0minus_VAJHU[ii][jj]*p_dec_0mplus_VAJHU; }
+      p_prod_0minus_dec_0minus_VAJHU.proddecme = p_prod_0minus_dec_0minus_VAJHU.prodme*p_prod_0mplus_dec_0mplus_VAJHU.decme; p_prod_0minus_dec_0minus_VAJHU.multiplyarray(p_prod_0mplus_dec_0mplus_VAJHU.decme);
+      p_prod_fa3_dec_fa3_VAJHU.proddecme = p_prod_fa3_dec_fa3_VAJHU.prodme*p_prod_0mplus_dec_0mplus_VAJHU.decme; p_prod_fa3_dec_fa3_VAJHU.multiplyarray(p_prod_0mplus_dec_0mplus_VAJHU.decme);
     }
     else{
-      p_prod_0minus_dec_0minus_VAJHU=p_prod_0minus_VAJHU*p_dec_0minus_VAJHU; for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) pArray_prod_0minus_dec_0minus_VAJHU[ii][jj] = pArray_prod_0minus_VAJHU[ii][jj]*p_dec_0minus_VAJHU; }
+      p_prod_0minus_dec_0minus_VAJHU.proddecme = p_prod_0minus_dec_0minus_VAJHU.prodme*p_prod_0minus_dec_0minus_VAJHU.decme; p_prod_0minus_dec_0minus_VAJHU.multiplyarray(p_prod_0minus_dec_0minus_VAJHU.decme);
+      p_prod_fa3_dec_fa3_VAJHU.proddecme = p_prod_fa3_dec_fa3_VAJHU.prodme*p_prod_fa3_dec_fa3_VAJHU.decme; p_prod_fa3_dec_fa3_VAJHU.multiplyarray(p_prod_fa3_dec_fa3_VAJHU.decme);
     }
 
     /***** MCFM *****/
@@ -1551,13 +1609,18 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
 
     if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }
     else{ mela.selfDHwwcoupl[0][gHIGGS_VV_1][0]=1; mela.differentiate_HWW_HZZ=true; }
-    mela.computeProdDecP(p_prod_0mplus_dec_0mplus_VAMCFM, false);
-    mela.getIORecord()->getUnweightedMEArray(pArray_prod_0mplus_dec_0mplus_VAMCFM);
+    mela.computeProdDecP(p_prod_0mplus_dec_0mplus_VAMCFM.proddecme, false);
+    mela.getIORecord()->getUnweightedMEArray(p_prod_0mplus_dec_0mplus_VAMCFM.mearray);
 
     if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }
     else{ mela.selfDHwwcoupl[0][gHIGGS_VV_4][0]=1; mela.differentiate_HWW_HZZ=true; }
-    mela.computeProdDecP(p_prod_0minus_dec_0minus_VAMCFM, false);
-    mela.getIORecord()->getUnweightedMEArray(pArray_prod_0minus_dec_0minus_VAMCFM);
+    mela.computeProdDecP(p_prod_0minus_dec_0minus_VAMCFM.proddecme, false);
+    mela.getIORecord()->getUnweightedMEArray(p_prod_0minus_dec_0minus_VAMCFM.mearray);
+
+    if (isZZWW!=2){ mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; if (isZZWW==1) mela.differentiate_HWW_HZZ=true; }
+    else{ mela.selfDHwwcoupl[0][gHIGGS_VV_1][0]=1; mela.selfDHwwcoupl[0][gHIGGS_VV_4][0]=1; mela.differentiate_HWW_HZZ=true; }
+    mela.computeProdDecP(p_prod_fa3_dec_fa3_VAMCFM.proddecme, false);
+    mela.getIORecord()->getUnweightedMEArray(p_prod_fa3_dec_fa3_VAMCFM.mearray);
 
     cout << "Production variables:\n";
     cout << "\tmJJ = " << mjj << endl;
@@ -1565,40 +1628,59 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
     cout << "\tJHUGen (mass, width): (" << mh << ", " << gah << ")" << endl;
     cout << "\tJHUGen propagator: " << propagator << endl;
     cout << "0mplus" << endl;
-    cout << "\tJHUGen decay-alone: " << p_dec_0mplus_VAJHU << endl;
-    cout << "\tJHUGen prod.-alone: " << p_prod_0mplus_VAJHU << endl;
-    cout << "\tJHUGen ME: " << p_prod_0mplus_dec_0mplus_VAJHU << endl;
-    cout << "\tMCFM ME: " << p_prod_0mplus_dec_0mplus_VAMCFM << endl;
+    cout << "\tJHUGen decay-alone: " << p_prod_0mplus_dec_0mplus_VAJHU.decme << endl;
+    cout << "\tJHUGen prod.-alone: " << p_prod_0mplus_dec_0mplus_VAJHU.prodme << endl;
+    cout << "\tJHUGen ME: " << p_prod_0mplus_dec_0mplus_VAJHU.proddecme << endl;
+    cout << "\tMCFM ME: " << p_prod_0mplus_dec_0mplus_VAMCFM.proddecme << endl;
     cout << "0minus" << endl;
-    cout << "\tJHUGen decay-alone: " << p_dec_0minus_VAJHU << endl;
-    cout << "\tJHUGen prod.-alone: " << p_prod_0minus_VAJHU << endl;
-    cout << "\tJHUGen ME: " << p_prod_0minus_dec_0minus_VAJHU << endl;
-    cout << "\tMCFM ME: " << p_prod_0minus_dec_0minus_VAMCFM << endl;
+    cout << "\tJHUGen decay-alone: " << p_prod_0minus_dec_0minus_VAJHU.decme << endl;
+    cout << "\tJHUGen prod.-alone: " << p_prod_0minus_dec_0minus_VAJHU.prodme << endl;
+    cout << "\tJHUGen ME: " << p_prod_0minus_dec_0minus_VAJHU.proddecme << endl;
+    cout << "\tMCFM ME: " << p_prod_0minus_dec_0minus_VAMCFM.proddecme << endl;
+    cout << "fa3" << endl;
+    cout << "\tJHUGen decay-alone: " << p_prod_fa3_dec_fa3_VAJHU.decme << endl;
+    cout << "\tJHUGen prod.-alone: " << p_prod_fa3_dec_fa3_VAJHU.prodme << endl;
+    cout << "\tJHUGen ME: " << p_prod_fa3_dec_fa3_VAJHU.proddecme << endl;
+    cout << "\tMCFM ME: " << p_prod_fa3_dec_fa3_VAMCFM.proddecme << endl;
 
     cout << "Arrays:\n0mplus" << endl;
     cout << "\tJHUGen" << endl;
-    for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << pArray_prod_0mplus_dec_0mplus_VAJHU[ii][jj]; cout << endl; }
+    p_prod_0mplus_dec_0mplus_VAJHU.printarray();
     cout << "\tMCFM" << endl;
-    for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << pArray_prod_0mplus_dec_0mplus_VAMCFM[ii][jj]; cout << endl; }
+    p_prod_0mplus_dec_0mplus_VAMCFM.printarray();
     cout << "\tJHUGen/MCFM Ratio" << endl;
     for (int ii=0; ii<nmsq; ii++){
       for (int jj=0; jj<nmsq; jj++){
         cout << '\t';
-        if (pArray_prod_0mplus_dec_0mplus_VAMCFM[ii][jj]!=0.) cout << pArray_prod_0mplus_dec_0mplus_VAJHU[ii][jj]/pArray_prod_0mplus_dec_0mplus_VAMCFM[ii][jj];
+        if (p_prod_0mplus_dec_0mplus_VAMCFM.mearray[ii][jj]!=0.) cout << p_prod_0mplus_dec_0mplus_VAJHU.mearray[ii][jj]/p_prod_0mplus_dec_0mplus_VAMCFM.mearray[ii][jj];
         else cout << 0;
       }
       cout << endl;
     }
     cout << "0minus" << endl;
     cout << "\tJHUGen" << endl;
-    for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << pArray_prod_0minus_dec_0minus_VAJHU[ii][jj]; cout << endl; }
+    p_prod_0minus_dec_0minus_VAJHU.printarray();
+    cout << "\tMCFM" << endl;
+    p_prod_0minus_dec_0minus_VAMCFM.printarray();
+    cout << "\tJHUGen/MCFM Ratio" << endl;
+    for (int ii=0; ii<nmsq; ii++){
+      for (int jj=0; jj<nmsq; jj++){
+        cout << '\t';
+        if (p_prod_0minus_dec_0minus_VAMCFM.mearray[ii][jj]!=0.) cout << p_prod_0minus_dec_0minus_VAJHU.mearray[ii][jj]/p_prod_0minus_dec_0minus_VAMCFM.mearray[ii][jj];
+        else cout << 0;
+      }
+      cout << endl;
+    }
+    cout << "fa3" << endl;
+    cout << "\tJHUGen" << endl;
+    p_prod_fa3_dec_fa3_VAJHU.printarray();
     cout << "\tMCFM" << endl;
-    for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << pArray_prod_0minus_dec_0minus_VAMCFM[ii][jj]; cout << endl; }
+    p_prod_fa3_dec_fa3_VAMCFM.printarray();
     cout << "\tJHUGen/MCFM Ratio" << endl;
     for (int ii=0; ii<nmsq; ii++){
       for (int jj=0; jj<nmsq; jj++){
         cout << '\t';
-        if (pArray_prod_0minus_dec_0minus_VAMCFM[ii][jj]!=0.) cout << pArray_prod_0minus_dec_0minus_VAJHU[ii][jj]/pArray_prod_0minus_dec_0minus_VAMCFM[ii][jj];
+        if (p_prod_fa3_dec_fa3_VAMCFM.mearray[ii][jj]!=0.) cout << p_prod_fa3_dec_fa3_VAJHU.mearray[ii][jj]/p_prod_fa3_dec_fa3_VAMCFM.mearray[ii][jj];
         else cout << 0;
       }
       cout << endl;

From 959dbf5e1460d16de94ada4a7f8b9b8e4d451705 Mon Sep 17 00:00:00 2001
From: Ulascan Sarica <ulascan90@gmail.com>
Date: Tue, 20 Dec 2016 22:02:34 +0100
Subject: [PATCH 3/4] Fix lepton interference from MELA

---
 MELA/src/TUtil.cc     | 14 ++++++++++++++
 MELA/test/testME_v2.c |  2 +-
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/MELA/src/TUtil.cc b/MELA/src/TUtil.cc
index 49abc1dc..a8e74157 100755
--- a/MELA/src/TUtil.cc
+++ b/MELA/src/TUtil.cc
@@ -1586,6 +1586,13 @@ bool TUtil::MCFM_chooser(
     breit_.mass3 =masses_mcfm_.zmass;
     breit_.width3=masses_mcfm_.zwidth;
 
+    vsymfact_.vsymfact=1.0;
+    interference_.interference=false;
+    if (hasZZ4fInterf && (leptonInterf==TVar::DefaultLeptonInterf || leptonInterf==TVar::InterfOn)){
+      vsymfact_.vsymfact=0.5;
+      interference_.interference=true;
+    }
+
     if (verbosity>=TVar::DEBUG) cout << "TUtil::MCFM_chooser: Setup is (production, process)=(" << TVar::ProductionName(production) << ", " << TVar::ProcessName(process) << ")" << endl;
 
   }
@@ -1691,6 +1698,13 @@ bool TUtil::MCFM_chooser(
     breit_.mass3 =masses_mcfm_.zmass;
     breit_.width3=masses_mcfm_.zwidth;
 
+    vsymfact_.vsymfact=1.0;
+    interference_.interference=false;
+    if (hasZZ4fInterf && (leptonInterf==TVar::DefaultLeptonInterf || leptonInterf==TVar::InterfOn)){
+      vsymfact_.vsymfact=0.5;
+      interference_.interference=true;
+    }
+
     if (verbosity>=TVar::DEBUG) cout << "TUtil::MCFM_chooser: Setup is (production, process)=(" << TVar::ProductionName(production) << ", " << TVar::ProcessName(process) << ")" << endl;
 
   }
diff --git a/MELA/test/testME_v2.c b/MELA/test/testME_v2.c
index f431bdd8..bd2b3291 100755
--- a/MELA/test/testME_v2.c
+++ b/MELA/test/testME_v2.c
@@ -1283,7 +1283,7 @@ void testME_ProdDec_MCFM_Ordering(int iSel, int jSel, int rSel, int sSel){
     << order[1] << endl;
 }
 
-void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0 /*1==ZZ, 2==WW*/, int vbfvhchannel=0 /*0==VBF, 1==VH*/, int hasInterf=0){
+void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0 /*1==ZZ, 2==WW*/, int vbfvhchannel=0 /*0==VBF, 1==VH*/, int hasInterf=0 /*0==2l2l, 1==4l*/){
   ofstream tout("testME_ProdDec_MCFM_JHUGen_Comparison_Ping.out");
   streambuf* coutbuf = cout.rdbuf();
   //cout.rdbuf(tout.rdbuf());

From cbd6a855a3875cb43cc0c287dbfcbe38c96342ab Mon Sep 17 00:00:00 2001
From: Ulascan Sarica <ulascan90@gmail.com>
Date: Wed, 21 Dec 2016 22:27:37 +0100
Subject: [PATCH 4/4] Upgrade to MCFM 702b3

---
 MELA/BuildFile.xml     |   4 +-
 MELA/src/TMCFMUtils.cc |  12 ++--
 MELA/src/TUtil.cc      |  16 ++++-
 MELA/test/loadMELA.C   |   2 +-
 MELA/test/testME_v2.c  | 151 ++++++++++++++++++++++++++++++++++++++++-
 5 files changed, 171 insertions(+), 14 deletions(-)

diff --git a/MELA/BuildFile.xml b/MELA/BuildFile.xml
index b3b5f3d3..54a0f00c 100755
--- a/MELA/BuildFile.xml
+++ b/MELA/BuildFile.xml
@@ -3,7 +3,7 @@
 <flags FFLAGS="-O3 -ffree-form -ffree-line-length-none -Dcompiler=2 -Dmela -fno-automatic -fno-f2c  -fPIC -g"/>
 
 <!--Some magic to retrieve the mcfm library over http, and link it hardcoding its path in the resulting library, so that LD_LIBRARY_PATH is not required. -->
-<Flags LDFLAGS="-Wl,-rpath=$(CMSSW_BASE)/src/ZZMatrixElement/MELA/data/$(SCRAM_ARCH) $(shell $(CMSSW_BASE)/src/ZZMatrixElement/MELA/data/retrieve.csh $(SCRAM_ARCH) mcfm_702b2) -L$(CMSSW_BASE)/src/ZZMatrixElement/MELA/data/$(SCRAM_ARCH) -lmcfm_702b2 -ljhugenmela" />
+<Flags LDFLAGS="-Wl,-rpath=$(CMSSW_BASE)/src/ZZMatrixElement/MELA/data/$(SCRAM_ARCH) $(shell $(CMSSW_BASE)/src/ZZMatrixElement/MELA/data/retrieve.csh $(SCRAM_ARCH) mcfm_702b3) -L$(CMSSW_BASE)/src/ZZMatrixElement/MELA/data/$(SCRAM_ARCH) -lmcfm_702b3 -ljhugenmela" />
 <Flags CPPDEFINES="_def_melatools_">
 
 </flags>
@@ -18,7 +18,7 @@
 
 <export>
   <lib name="1"/>
-  <use name="mcfm_702b2"/>   <!--For this to work, one needs to specify mcfm as a scram tool.-->
+  <use name="mcfm_702b3"/>   <!--For this to work, one needs to specify mcfm as a scram tool.-->
 </export>
 
 <export>
diff --git a/MELA/src/TMCFMUtils.cc b/MELA/src/TMCFMUtils.cc
index 1521a5cc..ad6959a0 100755
--- a/MELA/src/TMCFMUtils.cc
+++ b/MELA/src/TMCFMUtils.cc
@@ -30,7 +30,7 @@ void TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rS
         order[0]=0;
         order[1]=1;
         outFound=true;
-        cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", unswapped r, s = " << rSel << " " << sSel << endl;
+        //cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", unswapped r, s = " << rSel << " " << sSel << endl;
       }
       else if (
         (PDGHelpers::isAnUnknownJet(rSel) || rSel==hash.at(ih)[3])
@@ -40,7 +40,7 @@ void TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rS
         order[0]=1;
         order[1]=0;
         outFound=true;
-        cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", swapped r, s = " << rSel << " " << sSel << endl;
+        //cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", swapped r, s = " << rSel << " " << sSel << endl;
       }
     }
     // Final particles l/nu
@@ -61,7 +61,7 @@ void TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rS
         order[0]=0;
         order[1]=1;
         outFound=true;
-        cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", unswapped r, s = " << rSel << " " << sSel << endl;
+        //cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", unswapped r, s = " << rSel << " " << sSel << endl;
       }
       else if (
         (
@@ -77,7 +77,7 @@ void TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rS
         order[0]=1;
         order[1]=0;
         outFound=true;
-        cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", swapped r, s = " << rSel << " " << sSel << endl;
+        //cout << "Hash requested outgoing "<< hash.at(ih)[2] << " " << hash.at(ih)[3] << ", swapped r, s = " << rSel << " " << sSel << endl;
       }
       if (PDGHelpers::getCoupledVertex(rSel, sSel)!=PDGHelpers::getCoupledVertex(hash.at(ih)[2], hash.at(ih)[3])) outFound=false;
     }
@@ -204,8 +204,8 @@ std::vector<TMCFMUtils::intQuad_t> TMCFMUtils::Hash_QQVVQQStrong(){
   jcfg.push_back(TMCFMUtils::intQuad_t(1, 0, 2, 3)); // qbq->gg
   jcfg.push_back(TMCFMUtils::intQuad_t(2, 3, 0, 1)); // gg->qbq
   jcfg.push_back(TMCFMUtils::intQuad_t(2, 1, 0, 3)); // gqb->qbg
-  jcfg.push_back(TMCFMUtils::intQuad_t(1, 2, 0, 3)); // qbg->qbg
-  jcfg.push_back(TMCFMUtils::intQuad_t(3, 0, 2, 1)); // gq->gq
+  jcfg.push_back(TMCFMUtils::intQuad_t(2, 0, 1, 3)); // qbg->qbg
+  jcfg.push_back(TMCFMUtils::intQuad_t(1, 3, 2, 0)); // gq->gq
   jcfg.push_back(TMCFMUtils::intQuad_t(0, 3, 2, 1)); // qg->gq
 
   std::vector<TMCFMUtils::intQuad_t> pcfg;
diff --git a/MELA/src/TUtil.cc b/MELA/src/TUtil.cc
index a8e74157..27d42221 100755
--- a/MELA/src/TUtil.cc
+++ b/MELA/src/TUtil.cc
@@ -2542,7 +2542,13 @@ bool TUtil::MCFM_SetupParticleCouplings(
       // All of these require ZZ ordering since they use either ZZ ME or VV ME
       if (process==TVar::HSMHiggs || process==TVar::bkgZZ_SMHiggs || process==TVar::bkgZZ){
         ordering = pZOrder;
-        if (!hasZ1 || !hasZ2) result = false;
+        if (
+          !(
+          (hasZ1 && hasZ2)
+          ||
+          (process==TVar::HSMHiggs && hasW1 && hasW2)
+          )
+          ) result = false;
       }
       else if (process==TVar::HSMHiggs_WWZZ || process==TVar::bkgWWZZ_SMHiggs || process==TVar::bkgWWZZ){
         ordering = pZOrder;
@@ -2572,7 +2578,13 @@ bool TUtil::MCFM_SetupParticleCouplings(
       // All of these require ZZ ordering since they use either ZZ ME or VV ME
       else if (process==TVar::HSMHiggs || process==TVar::bkgZZ_SMHiggs || process==TVar::bkgZZ){
         ordering = pZOrder;
-        if (!hasZ1 || !hasZ2) result = false;
+        if (
+          !(
+          (hasZ1 && hasZ2)
+          ||
+          (process==TVar::HSMHiggs && hasW1 && hasW2)
+          )
+          ) result = false;
       }
       else if (process==TVar::HSMHiggs_WWZZ || process==TVar::bkgWWZZ_SMHiggs || process==TVar::bkgWWZZ){
         ordering = pZOrder;
diff --git a/MELA/test/loadMELA.C b/MELA/test/loadMELA.C
index a12f36d1..dc8bebbd 100755
--- a/MELA/test/loadMELA.C
+++ b/MELA/test/loadMELA.C
@@ -1,5 +1,5 @@
 {
-  TString LIBMCFM="libmcfm_702b2.so";
+  TString LIBMCFM="libmcfm_702b3.so";
   gSystem->AddIncludePath("-I$ROOFITSYS/include/");
   gSystem->AddIncludePath("-I$CMSSW_BASE/interface/");
   gSystem->AddIncludePath("-I$CMSSW_BASE/src/");
diff --git a/MELA/test/testME_v2.c b/MELA/test/testME_v2.c
index bd2b3291..5c2502b8 100755
--- a/MELA/test/testME_v2.c
+++ b/MELA/test/testME_v2.c
@@ -1283,8 +1283,8 @@ void testME_ProdDec_MCFM_Ordering(int iSel, int jSel, int rSel, int sSel){
     << order[1] << endl;
 }
 
-void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0 /*1==ZZ, 2==WW*/, int vbfvhchannel=0 /*0==VBF, 1==VH*/, int hasInterf=0 /*0==2l2l, 1==4l*/){
-  ofstream tout("testME_ProdDec_MCFM_JHUGen_Comparison_Ping.out");
+void testME_ProdDec_MCFM_JHUGen_WBFZZ_Comparison_Ping(int motherflavor=0, int isZZWW=0 /*1==ZZ, 2==WW*/, int vbfvhchannel=0 /*0==VBF, 1==VH*/, int hasInterf=0 /*0==2l2l, 1==4l*/){
+  ofstream tout("testME_ProdDec_MCFM_JHUGen_WBFZZ_Comparison_Ping.out");
   streambuf* coutbuf = cout.rdbuf();
   //cout.rdbuf(tout.rdbuf());
 
@@ -1326,6 +1326,10 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
   float wPOLE=4.07e-3;
 
   int idMother[2]={ 0 };
+  // FIX ME:
+  // Noticed that MCFM passes (3,2,0) == WW fusion for udb. The passing MEs are (-2,2)/(2,-2)/(-2,4)/(4,-2) -> (2,-1)
+  // Noticed that MCFM passes (4,2,0) == WW fusion for dub. The passing MEs are (-1,1)/(1,-1)/(-1,3)/(3,-1) -> (1,-2)
+  // Noticed that MCFM passes (5,1,0) == ZZ fusion for uub. The passing MEs are (-2,1..3.4.5)/(2,-1..-3.-4.-5)/(-4,3)/(3,-4) -> (2,-2). Initial (-2,4)==(-2,2) etc; (-2,1)==(-2,3) etc.
   // VBF ZZ(+)WW
   if (motherflavor==1){ idMother[0]=2; idMother[1]=1; } // Passed: (1,1,0);(1,2,0) -> Check!
   else if (motherflavor==2){ idMother[0]=-2; idMother[1]=-1; } // Passed: (2,1,0);(2,2,0) -> Check!
@@ -1339,7 +1343,7 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
   else if (motherflavor==8){ idMother[0]=-1; idMother[1]=1; } // Passed: (8,1,0);(8,2,0);(8,1,1) -> Check
 
   bool doEval=true;
-  TVar::VerbosityLevel verbosity = TVar::DEBUG;
+  TVar::VerbosityLevel verbosity = TVar::ERROR;
   TVar::Production prod;
   if (vbfvhchannel==0) prod=TVar::JJVBF_S;
   else if (vbfvhchannel==1){
@@ -1696,6 +1700,147 @@ void testME_ProdDec_MCFM_JHUGen_Comparison_Ping(int motherflavor=0, int isZZWW=0
 }
 
 
+void testME_ProdDec_MCFM_JHUGen_JJQCDZZ_Comparison_Ping(int motherflavor=0, int hasInterf=0 /*0==2l2l, 1==4l*/){
+  ofstream tout("testME_ProdDec_MCFM_JHUGen_JJQCDZZ_Comparison_Ping.out");
+  streambuf* coutbuf = cout.rdbuf();
+  //cout.rdbuf(tout.rdbuf());
+
+  struct mcfmme{
+    float proddecme;
+    double mearray[nmsq][nmsq];
+    mcfmme(){
+      proddecme=0;
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) mearray[ii][jj]=0; }
+    }
+    void multiplyarray(const float val){
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) mearray[ii][jj]*=val; }
+    }
+    void printarray(){
+      for (int ii=0; ii<nmsq; ii++){ for (int jj=0; jj<nmsq; jj++) cout << '\t' << mearray[ii][jj]; cout << endl; }
+    }
+  };
+
+  int erg_tev=13;
+  float mPOLE=125.0;
+  float wPOLE=4.07e-3;
+
+  int idMother[2]={ 0 };
+  // Ordering logic that involves gluons needs to be revised: case 19 only gives qqb/qbq/gq/qg, not qbg or gqb
+  if (motherflavor==1){ idMother[0]=2; idMother[1]=-2; }
+  else if (motherflavor==2){ idMother[0]=-2; idMother[1]=2; }
+  else if (motherflavor==3){ idMother[0]=1; idMother[1]=-1; }
+  else if (motherflavor==4){ idMother[0]=-1; idMother[1]=1; }
+
+  else if (motherflavor==5){ idMother[0]=1; idMother[1]=21; }
+  else if (motherflavor==6){ idMother[0]=21; idMother[1]=1; }
+  else if (motherflavor==7){ idMother[0]=-1; idMother[1]=21; }
+  else if (motherflavor==8){ idMother[0]=21; idMother[1]=-1; }
+  else if (motherflavor==9){ idMother[0]=2; idMother[1]=21; }
+  else if (motherflavor==10){ idMother[0]=21; idMother[1]=2; }
+  else if (motherflavor==11){ idMother[0]=-2; idMother[1]=21; }
+  else if (motherflavor==12){ idMother[0]=21; idMother[1]=-2; }
+  else if (motherflavor==13){ idMother[0]=21; idMother[1]=21; }
+
+  else if (motherflavor==14){ idMother[0]=2; idMother[1]=0; }
+  else if (motherflavor==15){ idMother[0]=-2; idMother[1]=0; }
+  else if (motherflavor==16){ idMother[0]=1; idMother[1]=0; }
+  else if (motherflavor==17){ idMother[0]=-1; idMother[1]=0; }
+  else if (motherflavor==19){ idMother[0]=21; idMother[1]=0; }
+
+  else if (motherflavor==20){ idMother[0]=0; idMother[1]=0; }
+
+  TVar::VerbosityLevel verbosity = TVar::ERROR;
+  TVar::Production prod = TVar::JJQCD;
+
+  Mela mela(erg_tev, mPOLE, verbosity);
+
+  mcfmme p_prod_JJQCD_VAMCFM;
+
+  float mzz = 0;
+  float mjj = 0;
+
+  float pingMom[8][4]={
+    { 0, 0, -865.37881546721542, -865.37881546721542 },
+    { 0, 0, 624.03396598421773, -624.03396598421773 },
+    { 7.6145299215002638, -17.259247740062808, 9.4660586470659975, 21.106135714241464 },
+    { 90.901719112641416, -69.683681833050798, 32.066319224729980, 118.94194752090492 },
+    { 78.476352131782917, -35.264818847819797, -8.8615639484695272, 86.490881645951262 },
+    { 191.68369742375290, -197.85205601463366, 100.99437243828194, 293.40746273989180 },
+    { -131.59521398083137, 330.56000090294270, 437.01695094737875, 563.53440884737279 },
+    { -237.08108460884614, -10.500196467375645, -329.33728782598945, 405.93194498307093 }
+  };
+
+  int GenLep1Id=0, GenLep2Id=0, GenLep3Id=0, GenLep4Id=0;
+  GenLep1Id=13;
+  GenLep2Id=-13;
+  if (hasInterf==0){
+    GenLep3Id=11;
+    GenLep4Id=-11;
+  }
+  else{
+    GenLep3Id=GenLep1Id;
+    GenLep4Id=GenLep2Id;
+  }
+  mela.setCandidateDecayMode(TVar::CandidateDecay_ZZ);
+  int idOrdered[4] ={ GenLep1Id, GenLep2Id, GenLep3Id, GenLep4Id };
+
+  SimpleParticleCollection_t mothers;
+  for (unsigned int ip=0; ip<2; ip++){
+    mothers.push_back(
+      SimpleParticle_t(
+      0,
+      TLorentzVector(pingMom[ip][0], pingMom[ip][1], pingMom[ip][2], pingMom[ip][3])
+      )
+      );
+  };
+  SimpleParticleCollection_t daughters;
+  for (unsigned int ip=2; ip<6; ip++){
+    daughters.push_back(
+      SimpleParticle_t(
+      idOrdered[ip-2],
+      TLorentzVector(pingMom[ip][0], pingMom[ip][1], pingMom[ip][2], pingMom[ip][3])
+      )
+      );
+  };
+  SimpleParticleCollection_t associated;
+  for (unsigned int ip=6; ip<8; ip++){
+    associated.push_back(
+      SimpleParticle_t(
+      idMother[ip-6], // Is this wrong? No, not really. We want to check all initial particles.
+      TLorentzVector(pingMom[ip][0], pingMom[ip][1], pingMom[ip][2], pingMom[ip][3])
+      )
+      );
+  };
+  mjj = (associated.at(0).second+associated.at(1).second).M();
+  mzz = (daughters.at(0).second+daughters.at(1).second+daughters.at(2).second+daughters.at(3).second).M();
+
+  mela.setInputEvent(&daughters, &associated, &mothers, true);
+
+  /***** MCFM *****/
+
+  mela.setProcess(TVar::bkgZZ, TVar::MCFM, prod);
+
+  mela.computeProdDecP(p_prod_JJQCD_VAMCFM.proddecme, false);
+  mela.getIORecord()->getUnweightedMEArray(p_prod_JJQCD_VAMCFM.mearray);
+
+  cout << "Production variables:\n";
+  cout << "\tmZZ = " << mzz << endl;
+  cout << "\tmJJ = " << mjj << endl;
+  cout << "JJQCD ZZ" << endl;
+  cout << "\tMCFM ME: " << p_prod_JJQCD_VAMCFM.proddecme << endl;
+  cout << "\tArray:" << endl;
+  p_prod_JJQCD_VAMCFM.printarray();
+
+  TUtil::PrintCandidateSummary(mela.getCurrentCandidate());
+
+  mela.resetInputEvent();
+
+
+  cout.rdbuf(coutbuf);
+  tout.close();
+}
+
+
 void testME_ProdDec_MCFM_Ping(int flavor=2){
   ofstream tout("testME_ProdDec_MCFM_Ping.out");
   streambuf* coutbuf = cout.rdbuf();