diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000..ce947e0
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,82 @@
+stages:
+  - build
+  - test
+
+variables:
+  LCG_VERSION: "98"
+
+.production_image:
+  variables:
+    LCG_OS: x86_64-centos7
+  image: gitlab-registry.cern.ch/ci-tools/ci-worker:cc7
+  tags:
+    - cvmfs
+
+.lcg_setup:
+  before_script:
+    - set +e && source /cvmfs/sft.cern.ch/lcg/views/setupViews.sh LCG_$LCG_VERSION $LCG_OS-$LCG_COMPILER; set -e
+
+.build_template:
+  stage: build
+  extends:
+    - .lcg_setup
+  script:
+    - mkdir install
+    - mkdir build && cd build
+    - cmake -DCMAKE_BUILD_TYPE:STRING=$BUILD_TYPE -DCMAKE_INSTALL_PREFIX:PATH=$CI_PROJECT_DIR/install -DEVTGEN_HEPMC3:BOOL=ON -DEVTGEN_PYTHIA:BOOL=ON -DEVTGEN_PHOTOS:BOOL=ON -DEVTGEN_TAUOLA:BOOL=ON -DEVTGEN_BUILD_TESTS:BOOL=$BUILD_TESTS $CI_PROJECT_DIR
+    - cmake --build .
+    - cmake --build . --target install
+
+build_clang10_opt:
+  variables:
+    LCG_COMPILER: "clang10-opt"
+    BUILD_TYPE: "Release"
+    BUILD_TESTS: "OFF"
+  extends:
+    - .production_image
+    - .build_template
+
+build_gcc10_opt:
+  variables:
+    LCG_COMPILER: "gcc10-opt"
+    BUILD_TYPE: "Release"
+    BUILD_TESTS: "OFF"
+  extends:
+    - .production_image
+    - .build_template
+  allow_failure: true
+
+build_gcc9_opt:
+  variables:
+    LCG_COMPILER: "gcc9-opt"
+    BUILD_TYPE: "Release"
+    BUILD_TESTS: "OFF"
+  extends:
+    - .production_image
+    - .build_template
+
+build_gcc9_dbg:
+  variables:
+    LCG_COMPILER: "gcc9-dbg"
+    BUILD_TYPE: "Debug"
+    BUILD_TESTS: "ON"
+  extends:
+    - .production_image
+    - .build_template
+  artifacts:
+    paths:
+      - install
+    expire_in: 1 day
+
+test_placeholder:
+  stage: test
+  variables:
+    LCG_COMPILER: "gcc9-dbg"
+  extends:
+    - .production_image
+    - .lcg_setup
+  dependencies:
+    - build_gcc9_dbg
+  script:
+    - cd install/share/EvtGen/test
+    - ./testDecayModel jsonFiles/BppmunuQCD.json
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 18ca882..05a24cc 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,55 +1,60 @@
 
 ########################################################################
 # Copyright 1998-2020 CERN for the benefit of the EvtGen authors       #
 #                                                                      #
 # This file is part of EvtGen.                                         #
 #                                                                      #
 # EvtGen is free software: you can redistribute it and/or modify       #
 # it under the terms of the GNU General Public License as published by #
 # the Free Software Foundation, either version 3 of the License, or    #
 # (at your option) any later version.                                  #
 #                                                                      #
 # EvtGen is distributed in the hope that it will be useful,            #
 # but WITHOUT ANY WARRANTY; without even the implied warranty of       #
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        #
 # GNU General Public License for more details.                         #
 #                                                                      #
 # You should have received a copy of the GNU General Public License    #
 # along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     #
 ########################################################################
 
 
 # the test executables require ROOT for histogramming etc.
 find_package(ROOT QUIET COMPONENTS Hist)
+find_package(jsoncpp)
 if (ROOT_FOUND)
     add_library(root_interface INTERFACE IMPORTED)
     target_include_directories(root_interface INTERFACE ${ROOT_INCLUDE_DIRS})
     target_link_libraries(root_interface INTERFACE ${ROOT_LIBRARIES})
 
     # build each of the executables
-    foreach( test_exe evt_dalitz evtgenlhc_test1 example1 exampleWriteHepMC )
+    foreach( test_exe evt_dalitz evtgenlhc_test1 example1 exampleWriteHepMC testDecayModel )
         add_executable(${test_exe} ${test_exe}.cc)
         target_link_libraries(${test_exe} PRIVATE EvtGen)
         if (EVTGEN_HEPMC3)
             target_compile_definitions(${test_exe} PRIVATE EVTGEN_HEPMC3)
         endif()
         if( EVTGEN_PYTHIA OR EVTGEN_PHOTOS OR EVTGEN_TAUOLA )
             target_compile_definitions(${test_exe} PRIVATE EVTGEN_EXTERNAL)
             target_link_libraries(${test_exe} PRIVATE EvtGenExternal)
         endif()
         target_link_libraries(${test_exe} PRIVATE root_interface)
     endforeach()
 
+    # testDecayModel needs jsoncpp as well
+    target_link_libraries(testDecayModel PRIVATE jsoncpp_lib)
+
     # install the executables
-    install(TARGETS evt_dalitz evtgenlhc_test1 example1 exampleWriteHepMC
+    install(TARGETS evt_dalitz evtgenlhc_test1 example1 exampleWriteHepMC testDecayModel
         DESTINATION ${CMAKE_INSTALL_DATADIR}/EvtGen/test
         )
 
     # install the decay files, macros, scripts, etc.
     install(DIRECTORY exampleFiles DESTINATION ${CMAKE_INSTALL_DATADIR}/EvtGen/test)
+    install(DIRECTORY jsonFiles DESTINATION ${CMAKE_INSTALL_DATADIR}/EvtGen/test)
     install(FILES do_tests
         DESTINATION ${CMAKE_INSTALL_DATADIR}/EvtGen/test
         PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE)
 else()
     message(WARNING "Could not find ROOT, cannot build test executables")
 endif()
diff --git a/test/jsonFiles/BdeltapmunuQCD.json b/test/jsonFiles/BdeltapmunuQCD.json
new file mode 100644
index 0000000..68313b5
--- /dev/null
+++ b/test/jsonFiles/BdeltapmunuQCD.json
@@ -0,0 +1,17 @@
+{
+    "parent" : "B-",
+    "daughters" : ["Delta+", "anti-p-", "mu-", "anti-nu_mu"],
+    "model" : "BToDiBaryonlnupQCD",
+    "parameters" : [67.7, -280.0, -187.3, -840.1, -10.1, -157.0, 7.0e6],
+    "events" : 10000,
+    "histograms" : [
+	{"variable" : "mass", "title" : "mPP", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 1.5, "xmax" : 3.5},
+	{"variable" : "mass", "title" : "mMuNu", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 3.5},
+	{"variable" : "psumsq", "title" : "qSqMuNu", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 12.0},
+	{"variable" : "coshel", "title" : "cosThetaB", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "coshel", "title" : "cosThetaL", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}
+    ],
+    "outfile" : "BDeltapmunuQCD.root",
+    "reference" : "RefBDeltapmunuQCD.root"
+}
diff --git a/test/jsonFiles/BppmunuQCD.json b/test/jsonFiles/BppmunuQCD.json
new file mode 100644
index 0000000..2cea582
--- /dev/null
+++ b/test/jsonFiles/BppmunuQCD.json
@@ -0,0 +1,17 @@
+{
+    "parent" : "B-",
+    "daughters" : ["p+", "anti-p-", "mu-", "anti-nu_mu"],
+    "model" : "BToDiBaryonlnupQCD",
+    "parameters" : [67.7, -280.0, -187.3, -840.1, -10.1, -157.0, 3.0e6],
+    "events" : 10000,
+    "histograms" : [
+	{"variable" : "mass", "title" : "mPP", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 1.5, "xmax" : 3.5},
+	{"variable" : "mass", "title" : "mMuNu", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 3.5},
+	{"variable" : "psumsq", "title" : "qSqMuNu", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 12.0},
+	{"variable" : "coshel", "title" : "cosThetaB", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "coshel", "title" : "cosThetaL", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}
+    ],
+    "outfile" : "BppmunuQCD.root",
+    "reference" : "RefBppmunuQCD.root"
+}
diff --git a/test/jsonFiles/EtaEEPiPi.json b/test/jsonFiles/EtaEEPiPi.json
new file mode 100644
index 0000000..434e83d
--- /dev/null
+++ b/test/jsonFiles/EtaEEPiPi.json
@@ -0,0 +1,18 @@
+{
+    "parent" : "eta'",
+    "daughters" : ["e+", "e-", "pi+", "pi-"],
+    "model" : "ETA_LLPIPI",
+    "parameters" : [],
+    "events" : 10000,
+    "histograms" : [
+	{"variable" : "mass", "title" : "mPiPi", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+	{"variable" : "mass", "title" : "mLL"  , "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+	{"variable" : "mass", "title" : "mPiL1", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.85},
+	{"variable" : "mass", "title" : "mPiL2", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.85},
+	{"variable" : "coshel", "title" : "cosL", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "coshel", "title" : "cosPi", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}
+    ],
+    "outfile" : "EtaEEPiPi.root",
+    "reference" : "RefEtaEEPiPi.root"
+}
diff --git a/test/jsonFiles/EtaMuMuPiPi.json b/test/jsonFiles/EtaMuMuPiPi.json
new file mode 100644
index 0000000..42eeffe
--- /dev/null
+++ b/test/jsonFiles/EtaMuMuPiPi.json
@@ -0,0 +1,18 @@
+{
+    "parent" : "eta'",
+    "daughters" : ["mu+", "mu-", "pi+", "pi-"],
+    "model" : "ETA_LLPIPI",
+    "parameters" : [],
+    "events" : 10000,
+    "histograms" : [
+	{"variable" : "mass", "title" : "mPiPi", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+	{"variable" : "mass", "title" : "mLL"  , "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+	{"variable" : "mass", "title" : "mPiL1", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.85},
+	{"variable" : "mass", "title" : "mPiL2", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.85},
+	{"variable" : "coshel", "title" : "cosL", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "coshel", "title" : "cosPi", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}
+    ],
+    "outfile" : "EtaMuMuPiPi.root",
+    "reference" : "RefEtaMuMuPiPi.root"
+}
diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
new file mode 100644
index 0000000..817e8c0
--- /dev/null
+++ b/test/testDecayModel.cc
@@ -0,0 +1,357 @@
+
+/***********************************************************************
+* Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
+*                                                                      *
+* This file is part of EvtGen.                                         *
+*                                                                      *
+* EvtGen is free software: you can redistribute it and/or modify       *
+* it under the terms of the GNU General Public License as published by *
+* the Free Software Foundation, either version 3 of the License, or    *
+* (at your option) any later version.                                  *
+*                                                                      *
+* EvtGen is distributed in the hope that it will be useful,            *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
+* GNU General Public License for more details.                         *
+*                                                                      *
+* You should have received a copy of the GNU General Public License    *
+* along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
+***********************************************************************/
+
+#include "testDecayModel.hh"
+
+#include "EvtGen/EvtGen.hh"
+
+#include "EvtGenBase/EvtAbsRadCorr.hh"
+#include "EvtGenBase/EvtDecayBase.hh"
+#include "EvtGenBase/EvtId.hh"
+#include "EvtGenBase/EvtMTRandomEngine.hh"
+#include "EvtGenBase/EvtParticle.hh"
+#include "EvtGenBase/EvtParticleFactory.hh"
+#include "EvtGenBase/EvtPDL.hh"
+#include "EvtGenBase/EvtVector4R.hh"
+
+#ifdef EVTGEN_EXTERNAL
+#include "EvtGenExternal/EvtExternalGenList.hh"
+#endif
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <list>
+#include <vector>
+
+testDecayModel::testDecayModel( const std::string& jsonFile ) :
+    m_jsonFile( jsonFile ),
+    m_infoVect(),
+    m_histVect(),
+    m_nHistos( 0 )
+{
+}
+
+void testDecayModel::run() {
+
+    Json::Value config;
+    std::ifstream inputStr( m_jsonFile.c_str() );
+    inputStr >> config;
+    inputStr.close();
+
+    const std::string parentName = config["parent"].asString();
+
+    Json::Value jDaughters = config["daughters"];
+    std::vector<std::string> daugNames;
+    for ( auto daugN : jDaughters ) {
+	daugNames.push_back( daugN.asString() );
+    }
+
+    const std::string modelName = config["model"].asString();
+
+    Json::Value jParameters = config["parameters"];
+    std::vector<double> modelPars;
+    for ( auto par : jParameters ) {
+	modelPars.push_back( par.asDouble() );
+    }
+
+    // Create decay file based on json file input
+    const std::string decFile = createDecFile( parentName, daugNames,
+					       modelName, modelPars );
+
+    int nEvents = config["events"].asInt();
+    std::cout << "Number of events = " << nEvents << std::endl;
+
+    // ROOT output file
+    const std::string outFileName = config["outfile"].asString();
+    TFile* outFile = TFile::Open( outFileName.c_str(), "recreate" );
+
+    // Define histograms
+    defineHistos( config, outFile );
+
+    // Generate events and fill histograms
+    generateEvents( decFile, parentName, nEvents );
+
+    // Normalize histograms
+    for ( auto hist : m_histVect ) {
+	double area = hist->Integral();
+	if ( area > 0.0 ) { hist->Scale( 1.0/area ); }
+    }
+    
+    // Compare with reference histograms
+    const std::string refFileName = config["reference"].asString();
+    compareHistos( refFileName );
+    
+    // Write output
+    outFile->cd();
+    for ( auto hist : m_histVect ) {
+	hist->Write();
+    }
+    outFile->Close();
+
+}
+
+std::string testDecayModel::createDecFile( const std::string& parent,
+					   const std::vector<std::string>& daugNames,
+					   const std::string& modelName,
+					   const std::vector<double>& parameters ) const
+{
+
+    // Create (or overwrite) the decay file
+    std::string decName( modelName ); decName += ".dec";
+    std::ofstream decFile;
+    decFile.open( decName.c_str() );
+
+    decFile << "Decay " << parent << std::endl;
+    decFile << "1.0";
+    for ( auto daug : daugNames ) {
+	decFile << " " << daug;
+    }
+    decFile << " " << modelName;
+    for ( auto par : parameters ) {
+	decFile << " " << par;
+    }
+    decFile << ";" << std::endl;
+    decFile << "Enddecay" << std::endl;
+    decFile << "End" << std::endl;
+
+    decFile.close();
+    
+    return decName;
+    
+}
+
+void testDecayModel::defineHistos( Json::Value& config, TFile* outFile ) {
+
+    // Histogram information
+    Json::Value jHistos = config["histograms"];
+    size_t nHistos = jHistos.size();
+    
+    m_infoVect.reserve( nHistos );
+    m_histVect.reserve( nHistos );
+    
+    for ( auto hInfo : jHistos ) {
+
+	const std::string varName = hInfo["variable"].asString();
+	// Integer values that define what particles need to be used
+	// for invariant mass combinations or helicity angles etc
+	int d1 = hInfo["d1"].asInt();
+	int d2 = hInfo["d2"].asInt();
+	
+	testInfo info( varName, d1, d2 );
+	m_infoVect.push_back( info );
+	
+	const std::string varTitle = hInfo["title"].asString();
+	int nBins = hInfo["nbins"].asInt();
+	double xmin = hInfo["xmin"].asDouble();
+	double xmax = hInfo["xmax"].asDouble();
+
+	TString histName( varName.c_str() );
+	if ( d1 != 0 ) { histName += "_"; histName += d1; }
+	if ( d2 != 0 ) { histName += "_"; histName += d2; }
+	TH1D* hist = new TH1D( histName.Data(), varTitle.c_str(), nBins, xmin, xmax );
+	hist->SetDirectory( outFile );
+	m_histVect.push_back( hist );
+
+    }	
+
+    m_nHistos = m_histVect.size();
+
+}
+
+void testDecayModel::generateEvents( const std::string& decFile, const std::string& parentName,
+				     int nEvents ) {
+    
+    // Define the random number generator
+    auto randomEngine = std::make_unique<EvtMTRandomEngine>();
+
+    EvtAbsRadCorr* radCorrEngine = nullptr;
+    std::list<EvtDecayBase*> extraModels;
+
+#ifdef EVTGEN_EXTERNAL
+    bool convertPythiaCodes( false );
+    bool useEvtGenRandom( true );
+    EvtExternalGenList genList( convertPythiaCodes, "", "gamma", useEvtGenRandom );
+    radCorrEngine = genList.getPhotosModel();
+    extraModels = genList.getListOfModels();
+#endif
+    
+    EvtGen theGen( "../DECAY.DEC", "../evt.pdl", randomEngine.get(),
+		   radCorrEngine, &extraModels );
+
+    theGen.readUDecay( decFile.c_str() );
+
+    // Generate the decays
+    EvtId parId = EvtPDL::getId( parentName.c_str() );
+    for ( int i = 0; i < nEvents; i++ ) {
+
+	if ( i%1000 == 0 ) { std::cout << "Event " << nEvents - i << std::endl; }
+	
+	// Initial 4-momentum and particle
+	EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 );
+	EvtParticle* parent = EvtParticleFactory::particleFactory( parId, pInit );
+
+	theGen.generateDecay( parent );
+
+	// Store information
+	for ( size_t i = 0; i < m_nHistos; i++ ) {
+
+	    testInfo info = m_infoVect[i];
+	    double value = getValue( parent, info.getName(), info.getd1(), info.getd2() );
+	    
+	    TH1D* hist = m_histVect[i];
+	    if ( hist ) {
+		hist->Fill( value );
+	    }
+	    
+	}
+
+	//parent->printTree();
+	
+	// Cleanup
+	parent->deleteTree();
+	
+    }
+    
+}
+
+double testDecayModel::getValue( EvtParticle* parent, const std::string& varName,
+				 const int d1, const int d2 ) const {
+
+    double value( 0.0 );
+    if ( !parent ) { return value; }
+
+    // Daughters
+    const EvtParticle* par1 = d1 > 0 ? parent->getDaug( d1 - 1 ) : nullptr;
+    const EvtParticle* par2 = d2 > 0 ? parent->getDaug( d2 - 1 ) : nullptr;
+
+    // 4-momenta (in parent rest frame)
+    const EvtVector4R p1 = par1 != nullptr ? par1->getP4() : EvtVector4R();
+    const EvtVector4R p2 = par2 != nullptr ? par2->getP4() : EvtVector4R();
+
+    if ( !varName.compare( "mass" ) ) {
+
+	// Invariant mass
+	if ( d2 != 0 ) {
+	    // Invariant 4-mass combination of particles d1 and d2
+	    value = ( p1 + p2 ).mass();
+	} else {
+	    // Invariant mass of particle d1 only
+	    value = p1.mass();
+	}
+	
+    } else if ( !varName.compare( "psumsq" ) ) {
+
+	// Invariant momentum sum squared
+	value = ( p1 + p2 ).mass2();
+
+    } else if ( !varName.compare( "pdiffsq" ) ) {
+
+	// Invariant momentum difference squared
+	value = ( p1 - p2 ).mass2();
+
+    } else if ( !varName.compare( "mtm" ) ) {
+
+	// Momentum of particle d1
+	value = p1.d3mag();
+	
+    } else if ( !varName.compare( "coshel" ) ) {
+
+	// Cosine of helicity angle
+
+	// Resonance center-of-mass system (d1 and d2)
+	const EvtVector4R p12 = p1 + p2;
+	// Boost vector
+	const EvtVector4R boost( p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ), -p12.get( 3 ) );
+	// Momentum of particle d1 in resonance frame, p1Res
+	const EvtVector4R p1Res = boostTo( p1, boost );
+	// Cosine of angle between p1Res and momentum of resonance in parent frame
+	double p1ResMag = p1Res.d3mag();
+	double p12Mag = p12.d3mag();
+	if ( p1ResMag > 0.0 && p12Mag > 0.0 ) {
+	    value = p1Res.dot( p12 ) / ( p1ResMag * p12Mag );
+	}
+
+    } else if ( !varName.compare( "prob" ) ) {
+
+	// Decay probability
+	double* dProb = parent->decayProb();
+        if ( dProb ) {
+	    value = *dProb;
+        }
+	
+    }
+    
+    return value;
+    
+}
+
+void testDecayModel::compareHistos( const std::string& refFileName ) const {
+
+    // Compare histograms with the same name, calculating the chi-squared
+    TFile* refFile = TFile::Open( refFileName.c_str(), "read" );
+
+    if ( !refFile ) {
+	std::cout << "Could not open reference file " << refFileName << std::endl;
+	return;
+    }
+    
+    for ( auto hist : m_histVect ) {
+
+	const std::string histName = hist->GetName();
+	// Get equivalent reference histogram
+	const TH1* refHist = dynamic_cast<TH1*>( refFile->Get ( histName.c_str() ) );
+
+	if ( refHist ) {
+
+	    double chiSq( 0.0 );
+	    int nDof( 0 );
+	    int iGood( 0 );
+	    double pValue = refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" );
+	    std::cout << "Histogram " << histName << " chiSq/nDof = "
+		      << chiSq << "/" << nDof << ", pValue = " << pValue << std::endl;
+
+	} else {
+
+	    std::cout << "Could not find reference histogram " << histName << std::endl;
+
+	}
+	
+    }
+
+    refFile->Close();
+    
+}
+
+int main( int argc, char* argv[] ) {
+    
+    if ( argc != 2 ) {
+	std::cout << "Expecting one argument: json input file" << std::endl;
+	return -1;
+    }
+
+    const std::string jsonFile = argv[1];
+
+    testDecayModel test( jsonFile );
+    test.run();
+
+    return 0;
+
+}
diff --git a/test/testDecayModel.hh b/test/testDecayModel.hh
new file mode 100644
index 0000000..10b2196
--- /dev/null
+++ b/test/testDecayModel.hh
@@ -0,0 +1,87 @@
+
+/***********************************************************************
+* Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
+*                                                                      *
+* This file is part of EvtGen.                                         *
+*                                                                      *
+* EvtGen is free software: you can redistribute it and/or modify       *
+* it under the terms of the GNU General Public License as published by *
+* the Free Software Foundation, either version 3 of the License, or    *
+* (at your option) any later version.                                  *
+*                                                                      *
+* EvtGen is distributed in the hope that it will be useful,            *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
+* GNU General Public License for more details.                         *
+*                                                                      *
+* You should have received a copy of the GNU General Public License    *
+* along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
+***********************************************************************/
+
+#ifndef TEST_DECAY_MODEL_HH
+#define TEST_DECAY_MODEL_HH
+
+#include "json/json.h"
+
+#include "TFile.h"
+#include "TH1D.h"
+
+#include <string>
+#include <vector>
+
+class EvtParticle;
+
+class testInfo {
+
+  public:
+    testInfo( const std::string& name, int d1, int d2 ) :
+	m_name( name ),
+	m_d1( d1 ),
+	m_d2( d2 )
+    {;}
+    
+    const std::string getName() const { return m_name; }
+    int getd1() const { return m_d1; }
+    int getd2() const { return m_d2; }
+    
+  private:
+    std::string m_name;
+    int m_d1;
+    int m_d2;
+
+};
+
+class testDecayModel {
+
+  public:
+    
+    testDecayModel( const std::string& jsonFile );
+
+    void run();    
+
+  private:
+
+    std::string createDecFile( const std::string& parent,
+			       const std::vector<std::string>& daugNames,
+			       const std::string& model,
+			       const std::vector<double>& parameters ) const;
+
+    void defineHistos( Json::Value& config, TFile* theFile );
+
+    void generateEvents( const std::string& decFile, const std::string& parentName,
+			 int nEvents );
+
+    double getValue( EvtParticle* rootPart, const std::string& varName,
+		     const int d1, const int d2 ) const;
+    
+    void compareHistos( const std::string& refFileName ) const;
+    
+    std::string m_jsonFile;
+    std::vector<testInfo> m_infoVect;
+    std::vector<TH1D*> m_histVect;
+    size_t m_nHistos;
+
+    
+};
+
+#endif