diff --git a/Generators/TruthRivetTools/CMakeLists.txt b/Generators/TruthRivetTools/CMakeLists.txt
index 96953559755c0262402a270e07409d9985e0578d..18b97de981f37ee93a09acb2fabe58e81f801193 100644
--- a/Generators/TruthRivetTools/CMakeLists.txt
+++ b/Generators/TruthRivetTools/CMakeLists.txt
@@ -1,3 +1,4 @@
+# $Id: CMakeLists.txt 769074 2016-08-22 10:04:26Z krasznaa $
 ################################################################################
 # Package: TruthRivetTools
 ################################################################################
@@ -8,9 +9,9 @@ atlas_subdir( TruthRivetTools )
 # Declare the dependencies of the package:
 atlas_depends_on_subdirs(
    PUBLIC
-   Control/AthToolSupport/AsgTools 
+	Control/AthToolSupport/AsgTools 
    PRIVATE
-   GaudiKernel )
+	GaudiKernel )
 
 # External dependencies:
 find_package( HepMC )
@@ -24,11 +25,12 @@ atlas_disable_as_needed()
 
 # Component(s) in the package:
 atlas_add_library( TruthRivetToolsLib
-   TruthRivetTools/*.h Root/*.cxx
+	Root/*.cxx TruthRivetTools/HiggsTemplateCrossSections.cc TruthRivetTools/*.h
+   SHARED
    PUBLIC_HEADERS TruthRivetTools
    INCLUDE_DIRS ${HEPMC_INCLUDE_DIRS} ${RIVET_INCLUDE_DIRS} ${YODA_INCLUDE_DIRS}
    ${FASTJET_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}
-   LINK_LIBRARIES ${HEPMC_LIBRARIES} ${RIVET_LIBRARIES} ${YODA_LIBRARIES}
+	LINK_LIBRARIES ${HEPMC_LIBRARIES} ${RIVET_LIBRARIES} ${YODA_LIBRARIES}
    ${FASTJET_LIBRARIES} ${ROOT_LIBRARIES} AsgTools )
 
 atlas_add_component( TruthRivetTools
diff --git a/Generators/TruthRivetTools/Root/HiggsTemplateCrossSections.cxx b/Generators/TruthRivetTools/Root/HiggsTemplateCrossSections.cxx
deleted file mode 100644
index 9c2562c6573d3ccc1ca5714114123bbbb6a083a8..0000000000000000000000000000000000000000
--- a/Generators/TruthRivetTools/Root/HiggsTemplateCrossSections.cxx
+++ /dev/null
@@ -1,670 +0,0 @@
-/*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// System include(s):
-#include <cmath>
-#include <cstdlib>
-
-// HepMC include(s):
-#include <HepMC/GenEvent.h>
-#include <HepMC/GenParticle.h>
-
-// Rivet include(s):
-#include <Rivet/Analysis.hh>
-#include <Rivet/Particle.hh>
-#include <Rivet/Projections/FastJets.hh>
-
-// Local include(s):
-#include "TruthRivetTools/HiggsTemplateCrossSections.h"
-
-namespace {
-
-   /// Check whether the input particle has a child with a given PDGID
-   bool hasChild( const HepMC::GenParticle& ptcl, int pdgID ) {
-
-      for( auto child : Rivet::Particle( ptcl ).children() ) {
-         if( child.pdgId() == pdgID ) {
-            return true;
-         }
-      }
-      return false;
-   }
-
-   /// Check whether the input particle has a parent with a given PDGID
-   bool hasParent( const HepMC::GenParticle& ptcl, int pdgID ) {
-
-      for( auto parent : Rivet::particles( ptcl.production_vertex(),
-                                           HepMC::parents ) ) {
-         if( parent->pdg_id() == pdgID ) {
-            return true;
-         }
-      }
-      return false;
-   }
-
-   /// Follow a "propagating" particle and return its last instance
-   const Rivet::Particle& getLastInstance( const Rivet::Particle& ptcl ) {
-
-      if( ptcl.genParticle()->end_vertex() ) {
-         if( ! hasChild( *( ptcl.genParticle() ), ptcl.pdgId() ) ) {
-            return ptcl;
-         } else {
-            return getLastInstance( ptcl.children()[ 0 ] );
-         }
-      }
-      return ptcl;
-   }
-
-   /// Check whether particle p originate from any of the ptcls
-   bool originateFrom( const Rivet::Particle& p,
-                       const Rivet::Particles& ptcls ) {
-
-      const HepMC::GenVertex* prodVtx = p.genParticle()->production_vertex();
-      if( ! prodVtx ) {
-         return false;
-      }
-
-      // for each ancestor, check if it matches any of the input particles
-      for( auto ancestor : Rivet::particles( prodVtx, HepMC::ancestors ) ) {
-         for( auto part : ptcls ) {
-            if( ancestor == part.genParticle() ) {
-               return true;
-            }
-         }
-      }
-
-      // if we get here, no ancetor matched any input particle
-      return false; 
-   }
-
-   /// Check whether particle p originates from p2
-   bool originateFrom( const Rivet::Particle& p, const Rivet::Particle& p2 ) {
-
-      const Rivet::Particles ptcls = { p2 };
-      return originateFrom( p, ptcls );
-   }
-
-   /// Return true is particle decays to quarks
-   bool quarkDecay( const Rivet::Particle& p ) {
-
-      for( auto child : p.children() ) {
-         if( Rivet::PID::isQuark( child.pdgId() ) ) {
-            return true;
-         }
-      }
-      return false;
-   }
-
-   /// Return bin index of x given the provided bin edges
-   /// 0=first bin, -1=underflow bin.
-   int getBin( double x, const std::vector< double >& bins ) {
-      if (bins.size()==0||x<bins[0]) return -1; // should not happen!
-      for (size_t i=1;i<bins.size();++i) {
-         if (x<bins[i]) {
-            return i-1;
-         }
-      }
-      return bins.size()-1;
-   }
-
-   /// VBF topolog selection
-   ///
-   /// 0 = fail loose selction: m_jj > 400 GeV and Dy_jj > 2.8
-   /// 1 pass loose, but fail additional cut pT(Hjj)<25. 2 pass tight selection
-   ///
-   int vbfTopology( const Rivet::Jets& jets, const Rivet::Particle& higgs ) {
-
-      if (jets.size()<2) {
-         return 0;
-      }
-      const Rivet::FourMomentum &j1=jets[0].momentum(), &j2=jets[1].momentum();
-      bool VBFtopo = ( (j1+j2).mass() > 400.0 && 
-                       std::abs(j1.rapidity()-j2.rapidity()) > 2.8 );
-      return VBFtopo ? (j1+j2+higgs.momentum()).pt()<25 ? 2 : 1 : 0;
-   }
-
-   /// Whether the Higgs is produced in association with a vector boson (VH)
-   bool isVH( HTXS::HiggsProdMode p ) {
-      return p==HTXS::WH || p==HTXS::QQ2ZH || p==HTXS::GG2ZH;
-   }
-
-   /// Stage-0 HTXS categorization
-   HTXS::Stage0::Category getStage0Category( const HTXS::HiggsProdMode prodMode,
-                                             const Rivet::Particle& higgs,
-                                             const Rivet::Particle& V ) {
-
-      using namespace HTXS::Stage0;
-      int ctrlHiggs = std::abs(higgs.rapidity())<2.5;
-      // Special cases first, qq→Hqq
-      if ( (prodMode==HTXS::WH||prodMode==HTXS::QQ2ZH) && quarkDecay(V) ) {
-         return ctrlHiggs ? VH2HQQ : VH2HQQ_FWDH;
-      } else if ( prodMode==HTXS::GG2ZH && quarkDecay(V) ) {
-         return Category(HTXS::GGF*10 + ctrlHiggs);
-      }
-      // General case after
-      return  Category(prodMode*10 + ctrlHiggs);
-   }
-
-   /// Stage-1 categorization
-   HTXS::Stage1::Category getStage1Category( const HTXS::HiggsProdMode prodMode,
-                                             const Rivet::Particle &higgs,
-                                             const Rivet::Jets &jets,
-                                             const Rivet::Particle &V ) {
-      using namespace HTXS::Stage1;
-      int Njets=jets.size();
-      int ctrlHiggs = std::abs(higgs.rapidity())<2.5;
-      int fwdHiggs = !ctrlHiggs;
-      double pTj1 = jets.size() ? jets[0].momentum().pt() : 0;
-      int vbfTopo = vbfTopology(jets,higgs);
-
-      // 1. GGF Stage 1 categories
-      //    Following YR4 write-up: XXXXX
-      if (prodMode==HTXS::GGF || (prodMode==HTXS::GG2ZH && quarkDecay(V)) ) {
-         if (fwdHiggs) {
-            return GG2H_FWDH;
-         }
-         if (Njets==0) {
-            return GG2H_0J;
-         } else if (Njets==1) {
-            return Category(GG2H_1J_PTH_0_60+getBin(higgs.pt(),{0,60,120,200}));
-         } else if (Njets>=2) {
-            // events with pT_H>200 get priority over VBF cuts
-            if(higgs.pt()<=200){
-               if      (vbfTopo==2) return GG2H_VBFTOPO_JET3VETO;
-               else if (vbfTopo==1) return GG2H_VBFTOPO_JET3;
-            }
-            // Njets >= 2jets without VBF topology
-            return Category(GG2H_GE2J_PTH_0_60+getBin(higgs.pt(),
-                                                      {0,60,120,200}));
-         }
-      }
-      // 2. Electroweak qq->Hqq Stage 1 categories
-      else if (prodMode==HTXS::VBF || ( isVH(prodMode) && quarkDecay(V)) ) {
-         if (std::abs(higgs.rapidity())>2.5) {
-            return QQ2HQQ_FWDH;
-         }
-      	if (pTj1>200) {
-            return QQ2HQQ_PTJET1_GT200;
-         }
-         if (vbfTopo==2) {
-            return QQ2HQQ_VBFTOPO_JET3VETO;
-         }
-         if (vbfTopo==1) {
-            return QQ2HQQ_VBFTOPO_JET3;
-         }
-         const double mjj = ( jets.size()>1 ?
-                              (jets[0].mom()+jets[1].mom()).mass():0 );
-         if ( 60 < mjj && mjj < 120 ) {
-            return QQ2HQQ_VH2JET;
-         }
-      	return QQ2HQQ_REST;
-      }
-      // 3. WH->Hlv categories
-      else if (prodMode==HTXS::WH) {
-         if (fwdHiggs) {
-            return QQ2HLNU_FWDH;
-         } else if (V.pt()<150) {
-            return QQ2HLNU_PTV_0_150;
-         } else if (V.pt()>250) {
-            return QQ2HLNU_PTV_GT250;
-         }
-      	// 150 < pTV/GeV < 250
-      	return ( jets.size()==0 ? QQ2HLNU_PTV_150_250_0J :
-                  QQ2HLNU_PTV_150_250_GE1J );
-      }
-      // 4. qq->ZH->llH categories
-      else if (prodMode==HTXS::QQ2ZH) {
-         if (fwdHiggs) {
-            return QQ2HLL_FWDH;
-      	} else if (V.pt()<150) {
-            return QQ2HLL_PTV_0_150;
-      	} else if (V.pt()>250) {
-            return QQ2HLL_PTV_GT250;
-         }
-      	// 150 < pTV/GeV < 250
-      	return jets.size()==0 ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
-      }
-      // 5. gg->ZH->llH categories
-      else if (prodMode==HTXS::GG2ZH ) {
-         if (fwdHiggs) {
-            return GG2HLL_FWDH;
-         }
-         if (V.pt()<150) {
-            return GG2HLL_PTV_0_150;
-         } else if (jets.size()==0) {
-            return GG2HLL_PTV_GT150_0J;
-         }
-         return GG2HLL_PTV_GT150_GE1J;
-      }
-      // 6.ttH,bbH,tH categories
-      else if (prodMode==HTXS::TTH) {
-         return Category(TTH_FWDH+ctrlHiggs);
-      } else if (prodMode==HTXS::BBH) {
-         return Category(BBH_FWDH+ctrlHiggs);
-      } else if (prodMode==HTXS::TH ) {
-         return Category(TH_FWDH+ctrlHiggs);
-      }
-      return UNKNOWN;
-    }
-
-} // private namespace
-
-namespace Rivet {
-
-   HiggsTemplateCrossSections::HiggsTemplateCrossSections()
-      : Analysis( "HiggsTemplateCrossSections" ),
-        m_HiggsProdMode( HTXS::UNKNOWN ) {
-
-   }
-
-   HiggsClassification
-   HiggsTemplateCrossSections::error( HiggsClassification& cat,
-                                      HTXS::ErrorCode err, 
-                                      std::string msg,
-                                      int NmaxWarnings ) {
-
-      // Set the error, and keep statistics
-      cat.errorCode = err;
-      ++m_errorCount[err];
-
-      // Print warning message to the screen/log
-      static int Nwarnings = 0;
-      if( msg!="" && ++Nwarnings < NmaxWarnings ) {
-         MSG_WARNING( msg );
-      }
-
-      return cat;
-   }
-
-   HiggsClassification
-   HiggsTemplateCrossSections::classifyEvent( const Event& event,
-                                              const HTXS::HiggsProdMode prodMode ) {
-
-      if( m_HiggsProdMode == HTXS::UNKNOWN ) {
-         m_HiggsProdMode = prodMode; 
-      }
-
-      // the classification object
-      HiggsClassification cat;
-      cat.prodMode = prodMode;
-      cat.errorCode = HTXS::UNDEFINED;
-      cat.stage0_cat = HTXS::Stage0::UNKNOWN;
-      cat.stage1_cat_pTjet25GeV = HTXS::Stage1::UNKNOWN;
-      cat.stage1_cat_pTjet30GeV = HTXS::Stage1::UNKNOWN;
-
-      if( prodMode == HTXS::UNKNOWN ) {
-         return error( cat, HTXS::PRODMODE_DEFINED,
-                       "Unkown Higgs production mechanism. Cannot classify "
-                       "event. Classification for all events will most likely "
-                       "fail.");
-      }
-
-      /*****
-       * Step 1. 
-       *  Idenfify the Higgs boson and the hard scatter vertex
-       *  There should be only one of each.
-       */
-
-      GenVertex *HSvtx = event.genEvent()->signal_process_vertex();
-      int Nhiggs=0;
-      for ( const GenParticle *ptcl : Rivet::particles( event.genEvent() ) ) {
-
-         // a) Reject all non-Higgs particles
-         if ( !PID::isHiggs(ptcl->pdg_id()) ) {
-            continue;
-         }
-         // b) select only the final Higgs boson copy, prior to decay
-         if ( ptcl->end_vertex() && !hasChild(*ptcl,PID::HIGGS) ) {
-            cat.higgs = Particle(ptcl); ++Nhiggs;
-         }
-         // c) if HepMC::signal_proces_vertex is missing
-         //    set hard-scatter vertex based on first Higgs boson
-         if ( HSvtx==nullptr && ptcl->production_vertex() &&
-              !hasParent(*ptcl,PID::HIGGS) ) {
-            HSvtx = ptcl->production_vertex();
-         }
-      }
-
-      // Make sure things are in order so far
-      if( Nhiggs != 1 ) {
-         return error( cat, HTXS::HIGGS_IDENTIFICATION,
-                       "Current event has " + std::to_string(Nhiggs) +
-                       " Higgs bosons. There must be only one." );
-      }
-      if( cat.higgs.children().size() < 2 ) {
-         return error( cat, HTXS::HIGGS_DECAY_IDENTIFICATION,
-                       "Could not identify Higgs boson decay products." );
-      }
-
-      if( HSvtx == nullptr ) {
-         return error( cat, HTXS::HS_VTX_IDENTIFICATION,
-                       "Cannot find hard-scatter vertex of current event." );
-      }
-
-      /*****
-       * Step 2. 
-       *   Identify associated vector bosons
-       */
- 
-      // Find associated vector bosons
-      bool is_uncatdV = false;
-      Particles uncatV_decays;
-      FourMomentum uncatV_p4(0,0,0,0);
-      FourVector uncatV_v4(0,0,0,0);
-      int nWs=0, nZs=0, nTop=0;
-      if ( isVH(prodMode) ) {
-         for (auto ptcl:particles(HSvtx,HepMC::children)) {
-            if (PID::isW(ptcl->pdg_id())) { ++nWs; cat.V=Particle(ptcl); }
-            if (PID::isZ(ptcl->pdg_id())) { ++nZs; cat.V=Particle(ptcl); }
-         }
-         if(nWs+nZs>0) {
-            cat.V = getLastInstance(cat.V);
-         } else {
-            for (auto ptcl:particles(HSvtx,HepMC::children)) {
-               if (!PID::isHiggs(ptcl->pdg_id())) {
-                  uncatV_decays += Particle(ptcl);
-                  uncatV_p4 += Particle(ptcl).momentum();
-                  uncatV_v4 += Particle(ptcl).origin();
-               }
-            }
-            is_uncatdV = true; cat.V = Particle(24,uncatV_p4,uncatV_v4);
-         }
-      }
-
-      if ( !is_uncatdV ){
-
-         if( isVH(prodMode) && !cat.V.genParticle()->end_vertex() ) {
-            return error( cat, HTXS::VH_DECAY_IDENTIFICATION,
-                          "Vector boson does not decay!" );
-         }
-
-         if( isVH(prodMode) && cat.V.children().size()<2 ) {
-            return error( cat, HTXS::VH_DECAY_IDENTIFICATION,
-                          "Vector boson does not decay!" );
-         }
-	
-         if ( ( prodMode==HTXS::WH && (nZs>0||nWs!=1) ) ||
-              ( (prodMode==HTXS::QQ2ZH||prodMode==HTXS::GG2ZH) &&
-                (nZs!=1||nWs>0) ) ) {
-            return error( cat, HTXS::VH_IDENTIFICATION,
-                          "Found " + std::to_string(nWs) + " W-bosons and " +
-                          std::to_string(nZs) +
-                          " Z-bosons. Inconsitent with VH expectation." );
-         }
-      }
-
-      // Find and store the W-bosons from ttH->WbWbH
-      Particles Ws;
-      if ( prodMode==HTXS::TTH || prodMode==HTXS::TH ){
-         // loop over particles produced in hard-scatter vertex
-      	for ( auto ptcl : particles(HSvtx,HepMC::children) ) {
-            if ( !PID::isTop(ptcl->pdg_id()) ) {
-               continue;
-            }
-            ++nTop;
-            Particle top = getLastInstance(Particle(ptcl));
-            if ( top.genParticle()->end_vertex() ) {
-               for (auto child:top.children()) {
-                  if ( PID::isW(child.pdgId()) ) {
-                     Ws += child;
-                  }
-               }
-            }
-         }
-      }
-
-      // Make sure result make sense
-      if ( (prodMode==HTXS::TTH && Ws.size()<2) ||
-           (prodMode==HTXS::TH && Ws.size()<1 ) ) {
-         return error( cat, HTXS::TOP_W_IDENTIFICATION,
-                       "Failed to identify W-boson(s) from t-decay!" );
-      }
-
-      /*****
-       * Step 3.
-       *   Build jets
-       *   Make sure all stable particles are present
-       */
-
-      // Create a list of the vector bosons that decay leptonically
-      // Either the vector boson produced in association with the Higgs boson,
-      // or the ones produced from decays of top quarks produced with the Higgs
-      Particles leptonicVs;
-      if ( !is_uncatdV ){
-         if ( isVH(prodMode) && !quarkDecay(cat.V) ) {
-            leptonicVs += cat.V;
-         }
-      } else {
-         leptonicVs = uncatV_decays;
-      }
-      for ( auto W:Ws ) {
-         if ( W.genParticle()->end_vertex() && !quarkDecay(W) ) {
-            leptonicVs += W;
-         }
-      }
-
-      // Obtain all stable, final-state particles
-      const ParticleVector FS =
-         applyProjection<FinalState>(event, "FS").particles();
-      Particles hadrons;
-      FourMomentum sum(0,0,0,0), vSum(0,0,0,0), hSum(0,0,0,0);
-      for ( const Particle &p : FS ) {
-         // Add up the four momenta of all stable particles as a cross check
-         sum += p.momentum();
-         // ignore particles from the Higgs boson
-         if ( originateFrom(p,cat.higgs) ) { hSum += p.momentum(); continue; }
-         // Cross-check the V decay products for VH
-         if ( isVH(prodMode) && !is_uncatdV && originateFrom(p,Ws) ) {
-            vSum += p.momentum();
-         }
-         // ignore final state particles from leptonic V decays
-         if ( leptonicVs.size() && originateFrom(p,leptonicVs) ) {
-            continue;
-         }
-         // All particles reaching here are considered hadrons and will be used
-         // to build jets
-         hadrons += p;
-      }
-      
-      cat.p4decay_higgs = hSum;
-      cat.p4decay_V = vSum;
-
-      FinalState fps_temp;
-      FastJets jets(fps_temp, FastJets::ANTIKT, 0.4 );
-      jets.calc(hadrons);
-
-      cat.jets25 = jets.jetsByPt( Cuts::pT > 25.0 );
-      cat.jets30 = jets.jetsByPt( Cuts::pT > 30.0 );
- 
-      // check that four mometum sum of all stable particles satisfies
-      // momentum consevation
-      if ( sum.pt()>0.1 ) {
-         return error( cat, HTXS::MOMENTUM_CONSERVATION,
-                       "Four vector sum does not amount to pT=0, m=E=sqrt(s), "
-                       "but pT=" + std::to_string( sum.pt() ) +
-                       " GeV and m = " + std::to_string( sum.mass() ) +
-                       " GeV" );
-      }
-      
-      // check if V-boson was not included in the event record but decay
-      // particles were EFT contact interaction: return UNKNOWN for category but
-      // set all event/particle kinematics
-      if(is_uncatdV) {
-         return error( cat, HTXS::VH_IDENTIFICATION,
-                       "Failed to identify associated V-boson!" );
-      }
-
-      /*****
-       * Step 4.
-       *   Classify and save output
-       */
-      
-      // Apply the categorization categorization
-      cat.stage0_cat = getStage0Category(prodMode,cat.higgs,cat.V);
-      cat.stage1_cat_pTjet25GeV =
-         getStage1Category(prodMode,cat.higgs,cat.jets25,cat.V);
-      cat.stage1_cat_pTjet30GeV =
-         getStage1Category(prodMode,cat.higgs,cat.jets30,cat.V);      
-      cat.errorCode = HTXS::SUCCESS; ++m_errorCount[HTXS::SUCCESS];
-
-      return cat;
-   }
-
-   void
-   HiggsTemplateCrossSections::setHiggsProdMode( HTXS::HiggsProdMode prodMode ) {
-
-      m_HiggsProdMode = prodMode;
-      return;
-   }
-
-   void HiggsTemplateCrossSections::init() {
-
-      MSG_INFO( "==============================================================" );
-      MSG_INFO( "========  HiggsTemplateCrossSections Initialization  =========" );
-      MSG_INFO( "==============================================================" );
-
-      // check that the production mode has been set
-      // if running in standalone Rivet the production mode is set through an env
-      // variable
-      if( m_HiggsProdMode == HTXS::UNKNOWN ) {
-         const char *pm_env = getenv("HIGGSPRODMODE");
-         std::string pm( pm_env == nullptr ? "" : pm_env );
-
-         if      ( pm == "GGF"   ) m_HiggsProdMode = HTXS::GGF;
-         else if ( pm == "VBF"   ) m_HiggsProdMode = HTXS::VBF;
-         else if ( pm == "WH"    ) m_HiggsProdMode = HTXS::WH;
-         else if ( pm == "ZH"    ) m_HiggsProdMode = HTXS::QQ2ZH;
-         else if ( pm == "QQ2ZH" ) m_HiggsProdMode = HTXS::QQ2ZH;
-         else if ( pm == "GG2ZH" ) m_HiggsProdMode = HTXS::GG2ZH;
-         else if ( pm == "TTH"   ) m_HiggsProdMode = HTXS::TTH;
-         else if ( pm == "BBH"   ) m_HiggsProdMode = HTXS::BBH;
-         else if ( pm == "TH"    ) m_HiggsProdMode = HTXS::TH;
-         else {
-            MSG_WARNING( "No HIGGSPRODMODE shell variable found. "
-                         "Needed when running Rivet stand-alone." );
-         }
-      }
-
-      // Projections for final state particles
-      const FinalState FS; 
-      addProjection( FS, "FS" );
-
-      // initialize the histograms with for each of the stages
-      initializeHistos();
-      m_sumw = 0.0;
-      MSG_INFO( "==============================================================" );
-      MSG_INFO( "========             Higgs prod mode " << m_HiggsProdMode
-                << "              =========" );
-      MSG_INFO( "========          Sucessful Initialization           =========" );
-      MSG_INFO( "==============================================================" );
-
-      return;
-   }
-
-   void HiggsTemplateCrossSections::analyze( const Event& event ) {
-
-      // get the classification
-      HiggsClassification cat = classifyEvent( event, m_HiggsProdMode );
-
-      // Fill histograms: categorization --> linerize the categories
-      const double weight = event.weight();
-      m_sumw += weight;
-
-      int F=cat.stage0_cat%10, P=cat.stage1_cat_pTjet30GeV/100;
-      hist_stage0->fill( cat.stage0_cat/10*2+F, weight );
-      
-      // Stage 1 enum offsets for each production mode:
-      //   GGF=12, VBF=6, WH= 5, QQ2ZH=5, GG2ZH=4, TTH=2, BBH=2, TH=2
-      static vector<int> offset({0,1,13,19,24,29,33,35,37,39});
-      int off = offset[P];
-      hist_stage1_pTjet25->fill(cat.stage1_cat_pTjet25GeV%100 + off, weight);
-      hist_stage1_pTjet30->fill(cat.stage1_cat_pTjet30GeV%100 + off, weight);
-
-      // Fill histograms: variables used in the categorization
-      hist_pT_Higgs->fill(cat.higgs.pT(),weight);
-      hist_y_Higgs->fill(cat.higgs.rapidity(),weight);
-      hist_pT_V->fill(cat.V.pT(),weight);
-
-      hist_Njets25->fill(cat.jets25.size(),weight);
-      hist_Njets30->fill(cat.jets30.size(),weight);
-
-      // Jet variables. Use jet collection with pT threshold at 30 GeV
-      if (cat.jets30.size()) hist_pT_jet1->fill(cat.jets30[0].pt(),weight);
-      if (cat.jets30.size()>=2) {
-         const FourMomentum &j1 = cat.jets30[0].momentum(),
-            &j2 = cat.jets30[1].momentum();
-         hist_deltay_jj->fill(std::abs(j1.rapidity()-j2.rapidity()),weight);
-         hist_dijet_mass->fill((j1+j2).mass(),weight);
-         hist_pT_Hjj->fill((j1+j2+cat.higgs.momentum()).pt(),weight);
-      }
-
-      return;
-   }
-
-   void HiggsTemplateCrossSections::finalize() {
-
-      printClassificationSummary();
-      double sf = m_sumw>0?1.0/m_sumw:1.0;
-      for( auto hist:{hist_stage0,hist_stage1_pTjet25,hist_stage1_pTjet30,
-               hist_Njets25,hist_Njets30, hist_pT_Higgs,hist_y_Higgs,
-               hist_pT_V,hist_pT_jet1,hist_deltay_jj,hist_dijet_mass,
-               hist_pT_Hjj}) {
-         scale(hist, sf);
-      }
-
-      return;
-   }
-
-   void HiggsTemplateCrossSections::printClassificationSummary() {
-
-      MSG_INFO (" ====================================================== ");
-      MSG_INFO ("      Higgs Template X-Sec Categorization Tool          ");
-      MSG_INFO ("                Status Code Summary                     ");
-      MSG_INFO (" ====================================================== ");
-      bool allSuccess = (numEvents()==m_errorCount[HTXS::SUCCESS]);
-      if ( allSuccess ) {
-         MSG_INFO ("     >>>> All "<< m_errorCount[HTXS::SUCCESS]
-                   <<" events successfully categorized!");
-      } else {
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::SUCCESS]
-                   <<" events successfully categorized");
-         MSG_INFO ("     >>>> --> the following errors occured:");
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::PRODMODE_DEFINED]
-                   <<" had an undefined Higgs production mode.");
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::MOMENTUM_CONSERVATION]
-                   <<" failed momentum conservation.");
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::HIGGS_IDENTIFICATION]
-                   <<" failed to identify a valid Higgs boson.");
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::HS_VTX_IDENTIFICATION]
-                   <<" failed to identify the hard scatter vertex.");
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::VH_IDENTIFICATION]
-                   <<" VH: to identify a valid V-boson.");
-         MSG_INFO ("     >>>> "<< m_errorCount[HTXS::TOP_W_IDENTIFICATION]
-                   <<" failed to identify valid Ws from top decay.");
-      }
-      MSG_INFO (" ====================================================== ");
-      MSG_INFO (" ====================================================== ");
-
-      return;
-   }
-
-   void HiggsTemplateCrossSections::initializeHistos(){
-
-      hist_stage0          = bookHisto1D("HTXS_stage0",20,0,20);
-      hist_stage1_pTjet25  = bookHisto1D("HTXS_stage1_pTjet25",40,0,40);
-      hist_stage1_pTjet30  = bookHisto1D("HTXS_stage1_pTjet30",40,0,40);
-      hist_pT_Higgs    = bookHisto1D("pT_Higgs",80,0,400);
-      hist_y_Higgs     = bookHisto1D("y_Higgs",80,-4,4);
-      hist_pT_V        = bookHisto1D("pT_V",80,0,400);
-      hist_pT_jet1     = bookHisto1D("pT_jet1",80,0,400);
-      hist_deltay_jj   = bookHisto1D("deltay_jj",50,0,10);
-      hist_dijet_mass  = bookHisto1D("m_jj",50,0,2000);
-      hist_pT_Hjj      = bookHisto1D("pT_Hjj",50,0,250);
-      hist_Njets25     = bookHisto1D("Njets25",10,0,10);
-      hist_Njets30     = bookHisto1D("Njets30",10,0,10);
-
-      return;
-   }
-
-} // namespace Rivet    
diff --git a/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx b/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx
index 7ea0b66b92a5d428af1d5e6b0c41fea4e7a46f67..753b0b1aad41f238db48a43d33fe4ef556bce556 100644
--- a/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx
+++ b/Generators/TruthRivetTools/Root/HiggsTruthCategoryTool.cxx
@@ -2,110 +2,61 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
-// Local include(s):
+// TruthRivetTools includes
 #include "TruthRivetTools/HiggsTruthCategoryTool.h"
 
-namespace {
-
-   template< class T >
-   TLorentzVector makeTLV( const T& p4 ) {
-
-      return TLorentzVector( p4.px(), p4.py(), p4.pz(), p4.E() );
-   }
-
-   template< class T >
-   std::vector< TLorentzVector > makeTLVs( const T& p4vec ) {
-
-      std::vector< TLorentzVector > result;
-      for( const auto& p4 : p4vec ) {
-         result.push_back( makeTLV( p4 ) );
-      }
-      return result;
-   }
-
-   HTXS::HiggsClassification
-   rivet2Root( const Rivet::HiggsClassification& cat ) {
-
-      HTXS::HiggsClassification ret;
-
-      ret.prodMode = cat.prodMode;
-      ret.errorCode = cat.errorCode;
-      ret.higgs = makeTLV( cat.higgs );
-      ret.V = makeTLV( cat.V );
-      ret.p4decay_higgs = makeTLV( cat.p4decay_higgs );
-      ret.p4decay_V = makeTLV( cat.p4decay_V );
-      ret.jets25 = makeTLVs( cat.jets25 );
-      ret.jets30 = makeTLVs( cat.jets30 );
-      ret.stage0_cat = cat.stage0_cat;
-      ret.stage1_cat_pTjet25GeV = cat.stage1_cat_pTjet25GeV;
-      ret.stage1_cat_pTjet30GeV = cat.stage1_cat_pTjet30GeV;
-
-      return ret;
-    }
-
-}
-
 HiggsTruthCategoryTool::HiggsTruthCategoryTool( const std::string& name) 
-  : asg::AsgTool( name ),
-    m_isInitialized( false ),
-    m_anaHandler(),
-    m_higgsTCS() {
-
-   // cannot be set to true until the issue with the beam protons in the truth
-   // event record is fixed..
-   // see JIRA ticket: https://its.cern.ch/jira/browse/ATLASRECTS-3072?filter=-2
-   declareProperty( "OutputHistos", m_outHistos = false );
-   declareProperty( "ProductionMode", m_prodMode = HTXS::UNKNOWN );
+: asg::AsgTool( name ),
+  rivetAnaHandler(nullptr),
+  higgsTemplateCrossSections(nullptr)
+{
+  // cannot be set to true until the issue with the beam protons in the truth event record is fixed..
+  // see JIRA ticket: https://its.cern.ch/jira/browse/ATLASRECTS-3072?filter=-2
+  m_outHistos = false; 
+  m_isInitialized = false;
 }
 
 
 StatusCode HiggsTruthCategoryTool::initialize() {
-
-   // Greet the user:
-   ATH_MSG_INFO( "Initializing ..." );
-
-   // Set up the helper objects:
-   const HTXS::HiggsProdMode prodMode =
-      static_cast< HTXS::HiggsProdMode >( m_prodMode );
-   m_higgsTCS.setHiggsProdMode( prodMode );
-   m_anaHandler.addAnalysis( &m_higgsTCS );
-
-   // Return gracefully:
-   return StatusCode::SUCCESS;
+  ATH_MSG_INFO ("Initializing " << name() << "...");
+  // Rivet analysis :: Higgs truth event classifier class
+  higgsTemplateCrossSections = new Rivet::HiggsTemplateCrossSections();
+  // crreate an instance of the Rivet analysis handler
+  rivetAnaHandler = new Rivet::AnalysisHandler();
+  // Add the Higgs truth classifier class to the handler
+  rivetAnaHandler->addAnalysis(&(*higgsTemplateCrossSections));
+  return StatusCode::SUCCESS;
 }
 
 StatusCode HiggsTruthCategoryTool :: finalize () {
+  ATH_MSG_INFO (" ====================================================== ");
+  ATH_MSG_INFO (" ---- Finalizing" << name() << "...");
+  ATH_MSG_INFO (" ====================================================== ");
+  if ( !m_outHistos )
+    higgsTemplateCrossSections->printClassificationSummary( );
+  else{
+    // TODO:: update the tool properly deal with output files/paths
+    rivetAnaHandler->finalize();
+    rivetAnaHandler->writeData("HiggsTruthCategoryTool.yoda");
+  }
+  ATH_MSG_INFO (" ====================================================== ");
+  return StatusCode::SUCCESS;  
+}
 
-   // Finalise the helper object(s):
-   if( m_outHistos ) {
-      m_anaHandler.finalize();
-      m_anaHandler.writeData( "HiggsTruthCategoryTool.yoda" );
-   } else {
-      m_higgsTCS.finalize();
-   }
-
-   // Return gracefully:
-   return StatusCode::SUCCESS;  
+HTXS::HiggsClassification HiggsTruthCategoryTool :: getHiggsTruthCategoryObject (const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode){
+  if ( !m_isInitialized ) {
+    higgsTemplateCrossSections->setHiggsProdMode(prodMode); 
+    rivetAnaHandler->init(HepMCEvent); 
+    m_isInitialized=true;
+  }
+  // fill histos if flag is specified
+  if ( m_outHistos ) rivetAnaHandler->analyze(HepMCEvent);
+  // get the category output object containing the template cross section category,
+  // and Higgs, V-boson, jets 4-vectors
+  const Rivet::HiggsClassification htxs_cat_rivet = higgsTemplateCrossSections->classifyEvent(HepMCEvent,prodMode);  
+  return HTXS::Rivet2Root(htxs_cat_rivet);
 }
 
-HTXS::HiggsClassification
-HiggsTruthCategoryTool::getHiggsTruthCategoryObject( const HepMC::GenEvent& ev ) {
 
-  if( ! m_isInitialized ) {
-     m_anaHandler.init( ev ); 
-     m_isInitialized = true;
-  }
 
-  // fill histos if flag is specified
-  if( m_outHistos ) {
-     m_anaHandler.analyze( ev );
-  }
 
-  // get the category output object containing the template cross section
-  // category, and Higgs, V-boson, jets 4-vectors
-  const HTXS::HiggsProdMode prodMode =
-     static_cast< HTXS::HiggsProdMode >( m_prodMode );
-  const Rivet::HiggsClassification htxs_cat_rivet =
-     m_higgsTCS.classifyEvent( ev, prodMode );
-  return rivet2Root( htxs_cat_rivet );
-}
diff --git a/Generators/TruthRivetTools/TruthRivetTools/Enums.h b/Generators/TruthRivetTools/TruthRivetTools/Enums.h
deleted file mode 100644
index 9845f0d6e877c9d0fa65409c386477a49d670dac..0000000000000000000000000000000000000000
--- a/Generators/TruthRivetTools/TruthRivetTools/Enums.h
+++ /dev/null
@@ -1,187 +0,0 @@
-// Dear emacs, this is -*- c++ -*-
-
-/*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef TRUTHRIVETTOOLS_ENUMS_H
-#define TRUTHRIVETTOOLS_ENUMS_H
-
-/// Higgs Template Cross Section namespace
-namespace HTXS {
-
-   /// Error code: whether the classification was successful or failed
-   enum ErrorCode {
-
-      /// Undefined/unknown error
-      UNDEFINED=-99,
-      /// Successful classification
-      SUCCESS = 0,
-      /// Production mode not defined
-      PRODMODE_DEFINED = 1,
-      /// Failed momentum conservation
-      MOMENTUM_CONSERVATION = 2,
-      /// Failed to identify Higgs boson
-      HIGGS_IDENTIFICATION = 3,
-      /// Failed to identify Higgs boson decay products
-      HIGGS_DECAY_IDENTIFICATION = 4,
-      /// Failed to identify hard scatter vertex
-      HS_VTX_IDENTIFICATION = 5,
-      /// Failed to identify associated vector boson
-      VH_IDENTIFICATION = 6,
-      /// Failed to identify associated vector boson decay products
-      VH_DECAY_IDENTIFICATION = 7,
-      /// Failed to identify top decay
-      TOP_W_IDENTIFICATION = 8
-
-   }; // enum ErrorCode
-
-   /// Higgs production modes, corresponding to input sample
-   enum HiggsProdMode {
-
-      UNKNOWN = 0,
-      GGF = 1,
-      VBF = 2,
-      WH = 3,
-      QQ2ZH = 4,
-      GG2ZH = 5,
-      TTH = 6,
-      BBH = 7,
-      TH = 8 
-
-   }; // enum HiggsProdMode
-  
-   /// Additional identifier flag for TH production modes
-   enum tH_type {
-
-      noTH=0,
-      THQB=1,
-      TWH=2
-
-   }; // enum tH_type
-
-   /// Namespace for Stage-0 categorization
-   namespace Stage0 {
-
-      /// Stage-0 categorization
-      ///
-      /// Two-digit number of format PF, with P for process,
-      /// and F being 0 for |yH|>2.5 and 1 for |yH|<2.5.
-      ///
-      enum Category {
-
-         UNKNOWN  = 0,
-         GG2H_FWDH = 10,
-         GG2H = 11,
-         VBF_FWDH = 20,
-         VBF = 21,
-         VH2HQQ_FWDH = 22,
-         VH2HQQ = 23,
-         QQ2HLNU_FWDH = 30,
-         QQ2HLNU = 31,
-         QQ2HLL_FWDH = 40,
-         QQ2HLL = 41,
-         GG2HLL_FWDH = 50,
-         GG2HLL = 51,
-         TTH_FWDH = 60,
-         TTH = 61,
-         BBH_FWDH = 70,
-         BBH = 71,
-         TH_FWDH = 80,
-         TH = 81
-
-      }; // enum Category
-
-   } // namespace Stage0
- 
-  /// Namespace for Stage-1 categorization
-  namespace Stage1 {
-
-     /// Stage-1 categorization
-     ///
-     /// Three digit integer of format PF.
-     /// Where P is a digit representing the process
-     /// F is a unique integer ( F < 99 ) corresponding to each Stage1
-     /// phase-space region (bin)
-     ///
-     enum Category {
-
-        UNKNOWN  = 0,
-
-        /// @name Gluon fusion
-        /// @{
-        GG2H_FWDH = 100,
-        GG2H_VBFTOPO_JET3VETO = 101,
-        GG2H_VBFTOPO_JET3 = 102,
-        GG2H_0J = 103,
-        GG2H_1J_PTH_0_60 = 104,
-        GG2H_1J_PTH_60_120 = 105,
-        GG2H_1J_PTH_120_200 = 106,
-        GG2H_1J_PTH_GT200 = 107,
-        GG2H_GE2J_PTH_0_60 = 108,
-        GG2H_GE2J_PTH_60_120 = 109,
-        GG2H_GE2J_PTH_120_200 = 110,
-        GG2H_GE2J_PTH_GT200 = 111,
-        /// @}
-
-        /// @name VBF
-        /// @{
-        QQ2HQQ_FWDH = 200,
-        QQ2HQQ_VBFTOPO_JET3VETO = 201,
-        QQ2HQQ_VBFTOPO_JET3 = 202,
-        QQ2HQQ_VH2JET = 203,
-        QQ2HQQ_REST = 204,
-        QQ2HQQ_PTJET1_GT200 = 205,
-        /// @}
-
-        /// @name qq -> WH
-        /// @{
-        QQ2HLNU_FWDH = 300,
-        QQ2HLNU_PTV_0_150 = 301,
-        QQ2HLNU_PTV_150_250_0J = 302,
-        QQ2HLNU_PTV_150_250_GE1J = 303,
-        QQ2HLNU_PTV_GT250 = 304,
-        /// @}
-
-        /// @name qq -> ZH
-        /// @{
-        QQ2HLL_FWDH = 400,
-        QQ2HLL_PTV_0_150 = 401,
-        QQ2HLL_PTV_150_250_0J = 402,
-        QQ2HLL_PTV_150_250_GE1J = 403,
-        QQ2HLL_PTV_GT250 = 404,
-        /// @}
-
-        /// @name gg -> ZH
-        /// @{
-        GG2HLL_FWDH = 500,
-        GG2HLL_PTV_0_150 = 501,
-        GG2HLL_PTV_GT150_0J = 502,
-        GG2HLL_PTV_GT150_GE1J = 503,
-        /// @}
-
-        /// @name ttH
-        /// @{
-        TTH_FWDH = 600,
-        TTH = 601,
-        /// @}
-
-        /// @name bbH
-        /// @{
-        BBH_FWDH = 700,
-        BBH = 701,
-        /// @}
-
-        /// @name tH
-        /// @{
-        TH_FWDH = 800,
-        TH = 801
-        /// @}
-
-     }; // enum Category
-
-  } // namespace Stage1
-
-} // namespace HTXS
-
-#endif // TRUTHRIVETTOOLS_ENUMS_H
diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsClassification.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsClassification.h
deleted file mode 100644
index 4f0d2288b3ee59b5bcad70f6684915f022d0cc75..0000000000000000000000000000000000000000
--- a/Generators/TruthRivetTools/TruthRivetTools/HiggsClassification.h
+++ /dev/null
@@ -1,61 +0,0 @@
-// Dear emacs, this is -*- c++ -*-
-
-/*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef TRUTHRIVETTOOLS_HIGGSCLASSIFICATION_H
-#define TRUTHRIVETTOOLS_HIGGSCLASSIFICATION_H
-
-// Rivet include(s):
-#include <Rivet/Particle.hh>
-#include <Rivet/Jet.hh>
-
-// Local include(s):
-#include "TruthRivetTools/Enums.h"
-
-/// Rivet namespace
-namespace Rivet {
-
-   /// Structure holding information about the current event
-   ///
-   /// Four-momenta and event classification according to the
-   /// Higgs Template Cross Section
-   ///
-   struct HiggsClassification {
-
-      /// Higgs production mode
-      HTXS::HiggsProdMode prodMode;
-
-      /// The Higgs boson
-      Rivet::Particle higgs;
-      /// Vector boson produced in asscoiation with the Higgs
-      Rivet::Particle V;
-      /// The four momentum sum of all stable decay products orignating from
-      /// the Higgs boson
-      Rivet::FourMomentum p4decay_higgs;
-      /// The four momentum sum of all stable decay products orignating from
-      /// the vector boson in associated production
-      Rivet::FourMomentum p4decay_V;
-      /// Jets built ignoring Higgs decay products and leptons from V decays,
-      /// pT thresholds at 25 GeV and 30 GeV
-      Rivet::Jets jets25, jets30;
-
-      /// Stage-0 HTXS event classifcation
-      /// See: https://cds.cern.ch/record/2138079
-      HTXS::Stage0::Category stage0_cat;
-      /// Stage-1 HTXS event classifcation
-      /// See: https://cds.cern.ch/record/2138079
-      HTXS::Stage1::Category stage1_cat_pTjet25GeV;
-      /// Stage-1 HTXS event classifcation
-      /// See: https://cds.cern.ch/record/2138079
-      HTXS::Stage1::Category stage1_cat_pTjet30GeV;
-
-      /// Error code: classification was succesful or some error occured
-      HTXS::ErrorCode errorCode;
-
-   }; // class HiggsClassification
-
-} // namespace Rivet
-
-#endif // TRUTHRIVETTOOLS_HIGGSCLASSIFICATION_H
diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h
index f8db1e989f3e068c8bf7688321b5a3f806d4b34b..46d3ea0b60d6a816dac435d7c5a7508463330af8 100644
--- a/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h
+++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTemplateCrossSections.h
@@ -1,5 +1,3 @@
-// Dear emacs, this is -*- c++ -*-
-
 /*
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
@@ -7,95 +5,235 @@
 #ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H
 #define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H 1
 
-// Rivet include(s):
-#include "Rivet/Analysis.hh"
-#include "Rivet/Event.hh"
-
-// Local include(s):
-#include "TruthRivetTools/HiggsClassification.h"
-#include "TruthRivetTools/Enums.h"
-
+/// Higgs Template Cross Section namespace
+namespace HTXS {
+
+  /// Error code: whether the classification was successful or failed
+  enum ErrorCode {
+    UNDEFINED=-99,
+    SUCCESS = 0,               ///< successful classification
+    PRODMODE_DEFINED = 1,      ///< production mode not defined
+    MOMENTUM_CONSERVATION = 2, ///< failed momentum conservation
+    HIGGS_IDENTIFICATION = 3,  ///< failed to identify Higgs boson
+    HIGGS_DECAY_IDENTIFICATION = 4,  ///< failed to identify Higgs boson decay products
+    HS_VTX_IDENTIFICATION = 5, ///< failed to identify hard scatter vertex
+    VH_IDENTIFICATION = 6,     ///< failed to identify associated vector boson
+    VH_DECAY_IDENTIFICATION = 7,     ///< failed to identify associated vector boson decay products
+    TOP_W_IDENTIFICATION = 8   ///< failed to identify top decay
+  };
+
+  /// Higgs production modes, corresponding to input sample
+  enum HiggsProdMode {
+    UNKNOWN = 0,
+    GGF = 1, VBF = 2, WH = 3, QQ2ZH = 4, GG2ZH = 5,
+    TTH = 6, BBH = 7, TH = 8 
+  };
+  
+  /// Additional identifier flag for TH production modes
+  enum tH_type { noTH=0, THQB=1, TWH=2 };
+  
+  ///   Two digit number of format PF
+  ///   P is digit for the physics process
+  ///   and F is 0 for |yH|>2.5 and 11 for |yH|<2.5 ("in fiducial")
+
+  /// Namespace for Stage0 categorization
+  namespace Stage0 {
+    /// @enum Stage-0 ategorization: Two-digit number of format PF, with P for process and F being 0 for |yH|>2.5 and 1 for |yH|<2.5
+    enum Category {
+      UNKNOWN  = 0, GG2H_FWDH = 10, GG2H = 11, VBF_FWDH = 20, VBF = 21, VH2HQQ_FWDH = 22, VH2HQQ = 23,
+      QQ2HLNU_FWDH = 30, QQ2HLNU = 31, QQ2HLL_FWDH = 40, QQ2HLL = 41, GG2HLL_FWDH = 50, GG2HLL = 51,
+      TTH_FWDH = 60, TTH = 61, BBH_FWDH = 70, BBH = 71, TH_FWDH = 80, TH = 81 };
+  }
+ 
+  /// Categorization Stage 1:
+  /// Three digit integer of format PF
+  /// Where P is a digit representing the process
+  /// F is a unique integer ( F < 99 ) corresponding to each Stage1 phase-space region (bin)
+  namespace Stage1 {
+    enum Category {
+      UNKNOWN  = 0,
+      // Gluon fusion
+      GG2H_FWDH = 100,
+      GG2H_VBFTOPO_JET3VETO = 101, GG2H_VBFTOPO_JET3 = 102,
+      GG2H_0J   = 103,
+      GG2H_1J_PTH_0_60 = 104,      GG2H_1J_PTH_60_120 = 105,
+      GG2H_1J_PTH_120_200 = 106,   GG2H_1J_PTH_GT200 = 107,
+      GG2H_GE2J_PTH_0_60 = 108,      GG2H_GE2J_PTH_60_120 = 109,
+      GG2H_GE2J_PTH_120_200 = 110,   GG2H_GE2J_PTH_GT200 = 111,
+      // "VBF"
+      QQ2HQQ_FWDH = 200,
+      QQ2HQQ_VBFTOPO_JET3VETO = 201, QQ2HQQ_VBFTOPO_JET3 = 202,
+      QQ2HQQ_VH2JET = 203, QQ2HQQ_REST = 204, QQ2HQQ_PTJET1_GT200 = 205,
+      // qq -> WH
+      QQ2HLNU_FWDH = 300,
+      QQ2HLNU_PTV_0_150 = 301,
+      QQ2HLNU_PTV_150_250_0J = 302,
+      QQ2HLNU_PTV_150_250_GE1J = 303,
+      QQ2HLNU_PTV_GT250 = 304,
+      // qq -> ZH
+      QQ2HLL_FWDH = 400,
+      QQ2HLL_PTV_0_150 = 401,
+      QQ2HLL_PTV_150_250_0J = 402,
+      QQ2HLL_PTV_150_250_GE1J = 403,
+      QQ2HLL_PTV_GT250 = 404,
+      // gg -> ZH
+      GG2HLL_FWDH = 500,
+      GG2HLL_PTV_0_150 = 501,
+      GG2HLL_PTV_GT150_0J = 502,
+      GG2HLL_PTV_GT150_GE1J = 503,
+      // ttH
+      TTH_FWDH = 600, TTH = 601,
+      // bbH
+      BBH_FWDH = 700, BBH = 701,
+      // tH
+      TH_FWDH = 800, TH = 801
+    };
+  } // namespace Stage1
+
+  
+#ifdef ROOT_TLorentzVector
+
+    typedef TLorentzVector TLV;
+    typedef std::vector<TLV> TLVs;
+    
+    template <class vec4>
+      TLV MakeTLV(vec4 const p) { return TLV(p.px(),p.py(),p.pz(),p.E()); }
+    
+    template <class Vvec4>
+      inline TLVs MakeTLVs(Vvec4 const &rivet_jets){ 
+      TLVs jets; for ( auto jet:rivet_jets ) jets.push_back(MakeTLV(jet)); 
+      return jets; 
+    }
+    
+    // Structure holding information about the current event:
+    // Four-momenta and event classification according to the
+    // Higgs Template Cross Section
+    struct HiggsClassification {
+      // Higgs production mode
+      HTXS::HiggsProdMode prodMode;
+      // The Higgs boson
+      TLV higgs;
+      // The Higgs boson decay products
+      TLV p4decay_higgs;
+      // Associated vector bosons
+      TLV V;
+      // The V-boson decay products
+      TLV p4decay_V;
+      // Jets are built ignoring Higgs decay products and leptons from V decays
+      // jets with pT > 25 GeV and 30 GeV
+      TLVs jets25, jets30;
+      // Event categorization according to YR4 wrtietup
+      // https://cds.cern.ch/record/2138079
+      HTXS::Stage0::Category stage0_cat;
+      HTXS::Stage1::Category stage1_cat_pTjet25GeV;
+      HTXS::Stage1::Category stage1_cat_pTjet30GeV;
+      // Error code :: classification was succesful or some error occured
+      HTXS::ErrorCode errorCode;
+    };
+    
+    template <class category>
+      inline HTXS::HiggsClassification Rivet2Root(category const &htxs_cat_rivet){
+      HTXS::HiggsClassification cat;
+      cat.prodMode = htxs_cat_rivet.prodMode;
+      cat.errorCode = htxs_cat_rivet.errorCode;
+      cat.higgs = MakeTLV(htxs_cat_rivet.higgs);
+      cat.V = MakeTLV(htxs_cat_rivet.V);
+      cat.p4decay_higgs = MakeTLV(htxs_cat_rivet.p4decay_higgs);
+      cat.p4decay_V = MakeTLV(htxs_cat_rivet.p4decay_V);
+      cat.jets25 = MakeTLVs(htxs_cat_rivet.jets25);
+      cat.jets30 = MakeTLVs(htxs_cat_rivet.jets30);
+      cat.stage0_cat = htxs_cat_rivet.stage0_cat;
+      cat.stage1_cat_pTjet25GeV = htxs_cat_rivet.stage1_cat_pTjet25GeV;
+      cat.stage1_cat_pTjet30GeV = htxs_cat_rivet.stage1_cat_pTjet30GeV;
+      return cat;    
+    }
+    
+
+    
+    inline int HTXSstage1_to_HTXSstage1FineIndex(HTXS::Stage1::Category stage1, 
+						 HiggsProdMode prodMode, tH_type tH) {
+
+      if(stage1==HTXS::Stage1::Category::UNKNOWN) return 0;
+      int P = (int)(stage1 / 100);
+      int F = (int)(stage1 % 100);
+      // 1.a spit tH categories
+      if (prodMode==HiggsProdMode::TH) {
+	// check that tH splitting is valid for Stage-1 FineIndex
+	// else return unknown category
+	if(tH==tH_type::noTH) return 0;
+	// check if forward tH
+	int fwdH = F==0?0:1;
+	return (49 + 2*(tH-1) +fwdH);
+      }
+      // 1.b QQ2HQQ --> split into VBF, WH, ZH -> HQQ
+      // offset vector 1: input is the Higgs prodMode 
+      // first two indicies are dummies, given that this is only called for prodMode=2,3,4 
+      std::vector<int> pMode_offset = {0,0,13,19,25};
+      if (P==2) return (F + pMode_offset[prodMode]);
+      // 1.c remaining categories
+      // offset vector 2: input is the Stage-1 category P 
+      // third index is dummy, given that this is called for category P=0,1,3,4,5,6,7
+      std::vector<int> catP_offset = {0,1,0,31,36,41,45,47};
+      return (F + catP_offset[P]);
+    }
+
+    inline int HTXSstage1_to_HTXSstage1FineIndex(const HiggsClassification &stxs, 
+						 tH_type tH=noTH, bool jets_pT25 = false) {
+      HTXS::Stage1::Category stage1 = 
+	jets_pT25==false?stxs.stage1_cat_pTjet30GeV:
+	stxs.stage1_cat_pTjet25GeV;
+      return HTXSstage1_to_HTXSstage1FineIndex(stage1,stxs.prodMode,tH);
+    }
+    
+    inline int HTXSstage1_to_index(HTXS::Stage1::Category stage1) {
+      // the Stage-1 categories
+      int P = (int)(stage1 / 100);
+      int F = (int)(stage1 % 100);    
+      std::vector<int> offset{0,1,13,19,24,29,33,35,37,39};
+      // convert to linear values
+      return ( F + offset[P] );
+    }
+
+     
+#endif 
+
+} // namespace HTXS
+
+
+#ifdef RIVET_Particle_HH
+//#ifdef HIGGSTRUTHCLASSIFIER_HIGGSTRUTHCLASSIFIER_CC
+//#include "Rivet/Particle.hh"
 namespace Rivet {
 
-   /// Rivet routine for classifying MC events according to the Higgs template
-   /// cross section categories
-   ///
-   /// @author Jim Lacey (DESY) <james.lacey@cern.ch,jlacey@desy.de>
-   /// @author Dag Gillberg (Carleton University) <dag.gillberg@cern.ch>
-   ///
-   class HiggsTemplateCrossSections : public Analysis {
-
-   public:
-      /// Constructor
-      HiggsTemplateCrossSections();
-
-      /// @name Utility methods
-      /// Methods to identify the Higgs boson and
-      /// associated vector boson and to build jets
-      /// @{
-
-      /// Returns the classification object with the error code set.
-      /// Prints an warning message, and keeps track of number of errors
-      HiggsClassification
-      error( HiggsClassification& cat, HTXS::ErrorCode err, 
-             std::string msg = "", int NmaxWarnings = 20 );
-
-      /// @}
-
-      /// Main classificaion method.
-      HiggsClassification classifyEvent( const Event& event,
-                                         const HTXS::HiggsProdMode prodMode );
-
-      /// @name Default Rivet analysis methods and steering methods
-      /// @{
-
-      /// Sets the Higgs production mode
-      void setHiggsProdMode( HTXS::HiggsProdMode prodMode );
-
-      /// Default Rivet Analysis::init method
-      ///
-      /// Booking of histograms, initializing Rivet projection
-      /// Extracts Higgs production mode from shell variable if not set
-      /// manually using setHiggsProdMode
-      void init() override;
-
-      /// Perform the per-event analysis
-      void analyze( const Event& event ) override;
-
-      /// Finalize the analysis
-      void finalize() override;
-
-      /// @}
-
-   private:
-      /// @name Functions used during initialisation and finalisation
-      /// @{
-
-      void printClassificationSummary();
-      void initializeHistos();
-
-      /// @}
-
-      double m_sumw;
-      HTXS::HiggsProdMode m_HiggsProdMode;
-      std::map<HTXS::ErrorCode,size_t> m_errorCount;
-      Histo1DPtr hist_stage0;
-      Histo1DPtr hist_stage1_pTjet25, hist_stage1_pTjet30;
-      Histo1DPtr hist_pT_Higgs, hist_y_Higgs;
-      Histo1DPtr hist_pT_V, hist_pT_jet1;
-      Histo1DPtr hist_deltay_jj, hist_dijet_mass, hist_pT_Hjj;
-      Histo1DPtr hist_Njets25, hist_Njets30;
-
-   }; // class HiggsTemplateCrossSections
-
-   // the PLUGIN only needs to be decleared when running standalone Rivet
-   // and causes compilation / linking issues if included in Athena / RootCore
-   //check for Rivet environment variable RIVET_ANALYSIS_PATH
-#ifdef RIVET_ANALYSIS_PATH
-   // The hook for the plugin system
-   DECLARE_RIVET_PLUGIN( HiggsTemplateCrossSections );
+  /// @struct HiggsClassification
+  /// @brief Structure holding information about the current event:
+  ///        Four-momenta and event classification according to the
+  ///        Higgs Template Cross Section
+  struct HiggsClassification {
+    /// Higgs production mode
+    HTXS::HiggsProdMode prodMode;
+    /// The Higgs boson
+    Rivet::Particle higgs;
+    /// Vector boson produced in asscoiation with the Higgs
+    Rivet::Particle V;
+    /// The four momentum sum of all stable decay products orignating from the Higgs boson
+    Rivet::FourMomentum p4decay_higgs;
+    /// The four momentum sum of all stable decay products orignating from the vector boson in associated production
+    Rivet::FourMomentum p4decay_V;
+    /// Jets built ignoring Higgs decay products and leptons from V decays, pT thresholds at 25 GeV and 30 GeV
+    Rivet::Jets jets25, jets30;
+    /// Stage-0 HTXS event classifcation, see: https://cds.cern.ch/record/2138079
+    HTXS::Stage0::Category stage0_cat;
+    /// Stage-1 HTXS event classifcation, see: https://cds.cern.ch/record/2138079
+    HTXS::Stage1::Category stage1_cat_pTjet25GeV;
+    /// Stage-1 HTXS event classifcation, see: https://cds.cern.ch/record/2138079
+    HTXS::Stage1::Category stage1_cat_pTjet30GeV;
+    /// Error code: Whether classification was succesful or some error occured
+    HTXS::ErrorCode errorCode;
+  };
+} // namespace Rivet
 #endif
 
-} // namespace Rivet
 
-#endif // TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_H
+
+#endif
diff --git a/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h b/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h
index 90d4be1be86dc80117814a7588b7127aa8b61ad3..df538af59c8fe82f9fd3c7b00131724dd065676d 100644
--- a/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h
+++ b/Generators/TruthRivetTools/TruthRivetTools/HiggsTruthCategoryTool.h
@@ -1,85 +1,41 @@
-// Dear emacs, this is -*- c++ -*-
-
 /*
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
+/*
+ * Dual-use tool interface for Rivet routine for classifying MC events according to the Higgs template cross section categories
+ * Authors: Jim Lacey (Carleton University)
+ * <james.lacey@cern.ch,jlacey@physics.carleton.ca>
+ */
+
 #ifndef TRUTHRIVETTOOLS_HIGGSTRUTHCATEGORYTOOL_H
 #define TRUTHRIVETTOOLS_HIGGSTRUTHCATEGORYTOOL_H 1
 
-// ASG include(s):
-#include "AsgTools/AsgTool.h"
-
-// Rivet include(s):
+#include "TLorentzVector.h"
 #include "Rivet/AnalysisHandler.hh"
+#include "TruthRivetTools/HiggsTemplateCrossSections.cc"
 
-// Local include(s):
+#include "AsgTools/AsgTool.h"
 #include "TruthRivetTools/IHiggsTruthCategoryTool.h"
-#include "TruthRivetTools/HiggsTemplateCrossSections.h"
 
-/**
- * Dual-use tool for a Rivet routine for classifying MC events according to
- * the Higgs template cross section categories
- *
- * @author Jim Lacey (Carleton University) <james.lacey@cern.ch,jlacey@physics.carleton.ca>
- */
-class HiggsTruthCategoryTool : public asg::AsgTool, 
-                               public virtual IHiggsTruthCategoryTool {
-
-public:
-   /// Create an "Athena constructor" for the tool
-   ASG_TOOL_CLASS( HiggsTruthCategoryTool, IHiggsTruthCategoryTool )
 
-   /// AsgTool constructor
+class HiggsTruthCategoryTool 
+: public asg::AsgTool, 
+  public virtual IHiggsTruthCategoryTool 
+{ 
+ public: 
+   ASG_TOOL_CLASS( HiggsTruthCategoryTool , IHiggsTruthCategoryTool )
    HiggsTruthCategoryTool( const std::string& name );
-
-   /// @name Functions inherited from AsgTool
-   /// @{
-
-   /// Function initialising the tool
-   StatusCode initialize() override;
-   /// Function finalising the tool
-   StatusCode finalize() override;
-
-   /// @}
-
-   /// @name Function(s) inherited from IHiggsTruthCategoryTool
-   /// @{
-
-   HTXS::HiggsClassification
-   getHiggsTruthCategoryObject( const HepMC::GenEvent& HepMCEvent ) override;
-
-   /// @}
-
-private:
-   /// @name Tool properties
-   /// @{
-
-   /// Flag for writing out data during finalisation
-   bool m_outHistos;
-   /// Higgs production mode to use
-   int m_prodMode;
-
-   /// @}
-
-   /// @name Status flag(s)
-   /// @{
-
-   /// Flag showing whether the internal analysis handler was initialised
+   ~HiggsTruthCategoryTool() { };
+ public:
+   Rivet::AnalysisHandler *rivetAnaHandler; //!
+   Rivet::HiggsTemplateCrossSections *higgsTemplateCrossSections; //!
+   virtual StatusCode  initialize() override;
+   StatusCode finalize () override;
+   HTXS::HiggsClassification getHiggsTruthCategoryObject(const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode) override;
+ private:
    bool m_isInitialized;
+   bool m_outHistos;
+};
 
-   /// @}
-
-   /// @name Helper tool(s)
-   /// @{
-
-   /// The analysis handler object
-   Rivet::AnalysisHandler m_anaHandler;
-   /// The analysis to run
-   Rivet::HiggsTemplateCrossSections m_higgsTCS;
-
-   /// @}
-
-}; // class HiggsTruthCategoryTool
-
-#endif // !HIGGSTRUTHCLASSIFIER_HIGGSTRUTHCATEGORYTOOL_H
+#endif //> !HIGGSTRUTHCLASSIFIER_HIGGSTRUTHCATEGORYTOOL_H
diff --git a/Generators/TruthRivetTools/TruthRivetTools/IHiggsTruthCategoryTool.h b/Generators/TruthRivetTools/TruthRivetTools/IHiggsTruthCategoryTool.h
index 4a9ecad2fad286c3d5069f171caedd3c91a50cdd..77e6478002d89765379aca4fc52540353d549a81 100644
--- a/Generators/TruthRivetTools/TruthRivetTools/IHiggsTruthCategoryTool.h
+++ b/Generators/TruthRivetTools/TruthRivetTools/IHiggsTruthCategoryTool.h
@@ -1,84 +1,29 @@
-// Dear emacs, this is -*- c++ -*-
-
 /*
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
+/*
+ * Dual-use tool interface for Rivet routine for classifying MC events according to the Higgs template cross section categories
+ * Authors: Jim Lacey (Carleton University)
+ * <james.lacey@cern.ch,jlacey@physics.carleton.ca>
+ */
+
 #ifndef TRUTHRIVETTOOLS_IHIGGSTRUTHCATEGORYTOOL_H
 #define TRUTHRIVETTOOLS_IHIGGSTRUTHCATEGORYTOOL_H 1
 
-// ROOT include(s):
-#include <TLorentzVector.h>
-
-// ASG include(s):
 #include "AsgTools/IAsgTool.h"
+#include "TLorentzVector.h"
+#include "HiggsTemplateCrossSections.h"
+#include "HepMC/GenEvent.h"
 
-// Local include(s):
-#include "TruthRivetTools/Enums.h"
-
-// Forward declaration(s):
-namespace HepMC {
-   class GenEvent;
-}
-
-/// Higgs Template Cross Section namespace
-namespace HTXS {
-
-   /// Structure holding information about the current event
-   ///
-   /// Four-momenta and event classification according to the
-   /// Higgs Template Cross Section.
-   ///
-   struct HiggsClassification {
-
-      /// Higgs production mode
-      HTXS::HiggsProdMode prodMode;
-
-      /// The Higgs boson
-      TLorentzVector higgs;
-      /// The Higgs boson decay products
-      TLorentzVector p4decay_higgs;
-      /// Associated vector bosons
-      TLorentzVector V;
-      /// The V-boson decay products
-      TLorentzVector p4decay_V;
-      /// Jets are built ignoring Higgs decay products and leptons from V decays
-      /// jets with pT > 25 GeV and 30 GeV
-      std::vector< TLorentzVector > jets25, jets30;
-
-      /// Stage-0 HTXS event classifcation
-      /// See: https://cds.cern.ch/record/2138079
-      HTXS::Stage0::Category stage0_cat;
-      /// Stage-1 HTXS event classifcation
-      /// See: https://cds.cern.ch/record/2138079
-      HTXS::Stage1::Category stage1_cat_pTjet25GeV;
-      /// Stage-1 HTXS event classifcation
-      /// See: https://cds.cern.ch/record/2138079
-      HTXS::Stage1::Category stage1_cat_pTjet30GeV;
-
-      /// Error code: classification was succesful or some error occured
-      HTXS::ErrorCode errorCode;
-
-   }; // struct HiggsClassification
-
-} // namespace HTXS
-
-/**
- * Dual-use tool interface for Rivet routine for classifying MC events according
- * to the Higgs template cross section categories.
- *
- * @authors Jim Lacey (Carleton University) <james.lacey@cern.ch,jlacey@physics.carleton.ca>
- */
 class IHiggsTruthCategoryTool : public virtual asg::IAsgTool {
-
-public:
-   /// Declare the interface to Athena
-   ASG_TOOL_INTERFACE( IHiggsTruthCategoryTool )
-
-   /// Get the classification for the specifed generated event
-   virtual HTXS::HiggsClassification
-   getHiggsTruthCategoryObject( const HepMC::GenEvent& HepMCEvent ) = 0;
-
-}; // class IHiggsTruthCategoryTool
-
-#endif // !HIGGSTRUTHCLASSIFIER_IHIGGSTRUTHCATEGORYTOOL_H
+ public:
+  ASG_TOOL_INTERFACE( IHiggsTruthCategoryTool ) //declares the interface to athena
+    virtual ~IHiggsTruthCategoryTool() {};
+ public:
+  virtual StatusCode initialize() = 0;
+  virtual StatusCode finalize () = 0;  
+  virtual HTXS::HiggsClassification getHiggsTruthCategoryObject(const HepMC::GenEvent& HepMCEvent, const HTXS::HiggsProdMode prodMode)=0;
+};
+
+#endif //> !HIGGSTRUTHCLASSIFIER_IHIGGSTRUTHCATEGORYTOOL_H
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D1ExtraContent.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D1ExtraContent.py
index 59aaca5da87ab5ab479fddc2563a2c16b72dba11..50f11cc7c7f0cd2abb2eddd2f49d2dc43946951c 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D1ExtraContent.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D1ExtraContent.py
@@ -26,4 +26,5 @@ HIGG2D1ExtraContainersTruth=[
     "TruthParticles",
     "TruthVertices",
     "AntiKt4TruthJets",
+    "AntiKt4TruthWZJets",
     "MuonTruthParticles"]
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D2ExtraContent.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D2ExtraContent.py
index e81363502a0c7150b398d24354778b1b13c48f21..759af597b9b361689f60e0f86b45f4f9919c1a60 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D2ExtraContent.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/python/HIGG2D2ExtraContent.py
@@ -26,4 +26,5 @@ HIGG2D2ExtraContainersTruth=[
     "TruthParticles",
     "TruthVertices",
     "AntiKt4TruthJets",
+    "AntiKt4TruthWZJets",
     "MuonTruthParticles"]
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG1D1.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG1D1.py
index 252c45af5813bd14c44ebf2b2a6e1084589a89bc..b9164cf7d2f2640c69657ca6ee0b9c888d5698da 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG1D1.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG1D1.py
@@ -255,6 +255,7 @@ HIGG1D1Stream.AddItem("xAOD::EventShapeAuxInfo_v1#*")
 HIGG1D1SlimmingHelper.SmartCollections = ["Electrons",
                                           "Photons",
                                           "Muons",
+                                          "TauJets",
                                           "MET_Reference_AntiKt4EMTopo",
                                           "AntiKt4EMTopoJets",
                                           "AntiKt4EMPFlowJets",
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D1.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D1.py
index 3b84830366b1d4f18e4969fc10e65c2f7a9c84d5..fce98cc954afa1f99f518090a3b750f070342924 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D1.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D1.py
@@ -171,7 +171,7 @@ else:
 from AthenaCommon.GlobalFlags import globalflags
 print "HIGG2D1.py globalflags.DataSource()", globalflags.DataSource()
 
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     ToolSvc += HIGG2D1TruthThinningTool
     thinningTools.append(HIGG2D1TruthThinningTool)
 print "HIGG2D1.py thinningTools", thinningTools
@@ -198,7 +198,7 @@ if jobproperties.Beam.energy()==4000000.0:
     electronMuonTriggerRequirement=["EF_e12Tvh_medium1_mu8", "EF_e24vhi_loose1_mu8"]
 triggerRequirement=singleElectronTriggerRequirement+diElectronTriggerRequirement+singleMuonTriggerRequirement+diMuonTriggerRequirement+electronMuonTriggerRequirement
 # 8 TeV MC does not have trigger information
-SkipTriggerRequirement=((globalflags.DataSource()=='geant4') and (jobproperties.Beam.energy()==4000000.0))
+SkipTriggerRequirement=(DerivationFrameworkIsMonteCarlo and (jobproperties.Beam.energy()==4000000.0))
 print "HIGG2D1.py SkipTriggerRequirement", SkipTriggerRequirement
 if SkipTriggerRequirement:
     triggerRequirement=[]
@@ -243,16 +243,34 @@ if Do4LVertexing:
     ToolSvc += SkimmingToolHIGG2D1
     augmentationTools.append(SkimmingToolHIGG2D1)
 
+#=======================================
+# CREATE PRIVATE SEQUENCE
+#=======================================
+higg2d1Seq = CfgMgr.AthSequencer("HIGG2D1Sequence")
+
 #====================================================================
 # CREATE THE DERIVATION KERNEL ALGORITHM AND PASS THE ABOVE TOOLS  
 #====================================================================
 
 # The name of the kernel (LooseSkimKernel in this case) must be unique to this derivation
 from DerivationFrameworkCore.DerivationFrameworkCoreConf import DerivationFramework__DerivationKernel
-DerivationFrameworkJob += CfgMgr.DerivationFramework__DerivationKernel("HIGG2D1Kernel",
-                                                                       SkimmingTools = [SkimmingToolHIGG2D1],
-                                                                       ThinningTools = thinningTools, 
-                                                                       AugmentationTools = augmentationTools)
+higg2d1Seq += CfgMgr.DerivationFramework__DerivationKernel("HIGG2D1Kernel",
+                                                           SkimmingTools = [SkimmingToolHIGG2D1],
+                                                           ThinningTools = thinningTools, 
+                                                           AugmentationTools = augmentationTools)
+
+#====================================================================
+# Standard jets
+#====================================================================
+
+if not "HIGG2D1Jets" in OutputJets:
+    OutputJets["HIGG2D1Jets"] = []
+    reducedJetList = []
+    if jetFlags.useTruth:
+        reducedJetList += ["AntiKt4TruthJets", "AntiKt4TruthWZJets"]
+    replaceAODReducedJets(reducedJetList, higg2d1Seq, "HIGG2D1Jets")
+
+DerivationFrameworkJob += higg2d1Seq
 
 #====================================================================
 # Add the containers to the output stream - slimming done here
@@ -275,7 +293,7 @@ HIGG2D1SlimmingHelper.SmartCollections = ["Electrons",
 
 HIGG2D1SlimmingHelper.ExtraVariables = HIGG2D1ExtraContent
 HIGG2D1SlimmingHelper.AllVariables = HIGG2D1ExtraContainers
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     HIGG2D1SlimmingHelper.ExtraVariables += HIGG2D1ExtraContentTruth
     HIGG2D1SlimmingHelper.AllVariables += HIGG2D1ExtraContainersTruth
 
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D2.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D2.py
index 2704c8907ba803594d1ac0bf0fdeba967d8d11ed..53c648944f191f87686876396e5bdba73dffcf77 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D2.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D2.py
@@ -161,7 +161,7 @@ else:
 from AthenaCommon.GlobalFlags import globalflags
 print "HIGG2D2.py globalflags.DataSource()", globalflags.DataSource()
 
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     ToolSvc += HIGG2D2TruthThinningTool
     thinningTools.append(HIGG2D2TruthThinningTool)
 print "HIGG2D2.py thinningTools", thinningTools
@@ -188,7 +188,7 @@ if jobproperties.Beam.energy()==4000000.0:
     electronMuonTriggerRequirement=["EF_e12Tvh_medium1_mu8", "EF_e24vhi_loose1_mu8"]
 triggerRequirement=singleElectronTriggerRequirement+diElectronTriggerRequirement+singleMuonTriggerRequirement+diMuonTriggerRequirement+electronMuonTriggerRequirement
 # 8 TeV MC does not have trigger information
-SkipTriggerRequirement=((globalflags.DataSource()=='geant4') and (jobproperties.Beam.energy()==4000000.0))
+SkipTriggerRequirement=(DerivationFrameworkIsMonteCarlo and (jobproperties.Beam.energy()==4000000.0))
 print "HIGG2D2.py SkipTriggerRequirement", SkipTriggerRequirement
 if SkipTriggerRequirement:
     triggerRequirement=[]
@@ -228,16 +228,34 @@ if Do4LVertexing:
     ToolSvc += SkimmingToolHIGG2D2
     augmentationTools.append(SkimmingToolHIGG2D2)
 
+#=======================================
+# CREATE PRIVATE SEQUENCE
+#=======================================
+higg2d2Seq = CfgMgr.AthSequencer("HIGG2D2Sequence")
+
 #====================================================================
 # CREATE THE DERIVATION KERNEL ALGORITHM AND PASS THE ABOVE TOOLS  
 #====================================================================
 
 # The name of the kernel (LooseSkimKernel in this case) must be unique to this derivation
 from DerivationFrameworkCore.DerivationFrameworkCoreConf import DerivationFramework__DerivationKernel
-DerivationFrameworkJob += CfgMgr.DerivationFramework__DerivationKernel("HIGG2D2Kernel",
-                                                                       SkimmingTools = [SkimmingToolHIGG2D2],
-                                                                       ThinningTools = thinningTools,
-                                                                       AugmentationTools = augmentationTools)
+higg2d2Seq += CfgMgr.DerivationFramework__DerivationKernel("HIGG2D2Kernel",
+                                                           SkimmingTools = [SkimmingToolHIGG2D2],
+                                                           ThinningTools = thinningTools,
+                                                           AugmentationTools = augmentationTools)
+
+#====================================================================
+# Standard jets
+#====================================================================
+
+if not "HIGG2D2Jets" in OutputJets:
+    OutputJets["HIGG2D2Jets"] = []
+    reducedJetList = []
+    if jetFlags.useTruth:
+        reducedJetList += ["AntiKt4TruthJets", "AntiKt4TruthWZJets"]
+    replaceAODReducedJets(reducedJetList, higg2d2Seq, "HIGG2D2Jets")
+
+DerivationFrameworkJob += higg2d2Seq
 
 #====================================================================
 # Add the containers to the output stream - slimming done here
@@ -262,7 +280,7 @@ HIGG2D2SlimmingHelper.SmartCollections = ["Electrons",
 
 HIGG2D2SlimmingHelper.ExtraVariables = HIGG2D2ExtraContent
 HIGG2D2SlimmingHelper.AllVariables = HIGG2D2ExtraContainers
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     HIGG2D2SlimmingHelper.ExtraVariables += HIGG2D2ExtraContentTruth
     HIGG2D2SlimmingHelper.AllVariables += HIGG2D2ExtraContainersTruth
 
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D4.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D4.py
index 707553e53e16e8b148d5d76bdf093668e08f5689..369d00ffa33e37f8d280d1cb081aecc4678edbbd 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D4.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D4.py
@@ -162,7 +162,7 @@ else:
 from AthenaCommon.GlobalFlags import globalflags
 print "HIGG2D4.py globalflags.DataSource()", globalflags.DataSource()
 
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     ToolSvc += HIGG2D4TruthThinningTool
     thinningTools.append(HIGG2D4TruthThinningTool)
 print "HIGG2D4.py thinningTools", thinningTools
@@ -189,7 +189,7 @@ if jobproperties.Beam.energy()==4000000.0:
     electronMuonTriggerRequirement=["EF_e12Tvh_medium1_mu8", "EF_e24vhi_loose1_mu8"]
 triggerRequirement=singleElectronTriggerRequirement+diElectronTriggerRequirement+singleMuonTriggerRequirement+diMuonTriggerRequirement+electronMuonTriggerRequirement
 # 8 TeV MC does not have trigger information
-SkipTriggerRequirement=((globalflags.DataSource()=='geant4') and (jobproperties.Beam.energy()==4000000.0))
+SkipTriggerRequirement=(DerivationFrameworkIsMonteCarlo and (jobproperties.Beam.energy()==4000000.0))
 print "HIGG2D4.py SkipTriggerRequirement", SkipTriggerRequirement
 if SkipTriggerRequirement:
     triggerRequirement=[]
@@ -268,6 +268,8 @@ if not "HIGG2D4Jets" in OutputJets:
     OutputJets["HIGG2D4Jets"] = []
 
     reducedJetList = ["AntiKt2PV0TrackJets", "AntiKt4PV0TrackJets", "AntiKt10LCTopoJets"]
+    if jetFlags.useTruth:
+        reducedJetList += ["AntiKt4TruthJets", "AntiKt4TruthWZJets"]
     replaceAODReducedJets(reducedJetList, higg2d4Seq, "HIGG2D4Jets")
 
 #====================================================================
@@ -275,11 +277,9 @@ if not "HIGG2D4Jets" in OutputJets:
 #====================================================================
 
     if jetFlags.useTruth:
-        OutputJets["HIGG2D4Jets"].append("AntiKt4TruthJets")
-        OutputJets["HIGG2D4Jets"].append("AntiKt4TruthWZJets")
-        addTrimmedJets("AntiKt", 1.0, "TruthWZ", rclus=0.2, ptfrac=0.05, includePreTools=False, algseq=higg2d4Seq,outputGroup="HIGG2D4Jets")
+        addTrimmedJets("AntiKt", 1.0, "TruthWZ", rclus=0.2, ptfrac=0.05, mods="groomed", includePreTools=False, algseq=higg2d4Seq,outputGroup="HIGG2D4Jets")
 
-    addTrimmedJets("AntiKt", 1.0, "LCTopo", rclus=0.2, ptfrac=0.05, includePreTools=False, algseq=higg2d4Seq,outputGroup="HIGG2D4Jets")
+    addTrimmedJets("AntiKt", 1.0, "LCTopo", rclus=0.2, ptfrac=0.05, mods="lctopo_groomed", includePreTools=False, algseq=higg2d4Seq,outputGroup="HIGG2D4Jets")
 
 #====================================================================
 # Create variable-R trackjets and dress AntiKt10LCTopo with ghost VR-trkjet 
@@ -381,7 +381,7 @@ HIGG2D4SlimmingHelper.SmartCollections = ["Electrons",
 
 HIGG2D4SlimmingHelper.ExtraVariables = HIGG2D4ExtraContent
 HIGG2D4SlimmingHelper.AllVariables = HIGG2D4ExtraContainers
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     HIGG2D4SlimmingHelper.ExtraVariables += HIGG2D4ExtraContentTruth
     HIGG2D4SlimmingHelper.AllVariables += HIGG2D4ExtraContainersTruth
 HIGG2D4SlimmingHelper.ExtraVariables += JetTagConfig.GetExtraPromptVariablesForDxAOD()
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D5.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D5.py
index 859f89117d3963e0e9213bea2c95ed441db92ba9..9a579ff709972eb8b44dc0c8ab8e2f9a672b3dac 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D5.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG2D5.py
@@ -430,7 +430,7 @@ else:
 from AthenaCommon.GlobalFlags import globalflags
 print "HIGG2D5.py globalflags.DataSource()", globalflags.DataSource()
 
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     ToolSvc += HIGG2D5TruthThinningTool
     thinningTools.append(HIGG2D5TruthThinningTool)
 print "HIGG2D5.py thinningTools", thinningTools
@@ -441,10 +441,7 @@ print "HIGG2D5.py thinningTools", thinningTools
 
 from AthenaCommon.BeamFlags import jobproperties
 print "HIGG2D5.py jobproperties.Beam.energy()", jobproperties.Beam.energy()
-# SkipTriggerRequirement=((globalflags.DataSource()=='geant4') and (jobproperties.Beam.energy()==4000000.0))
-#8 TeV MC does not have trigger information
-#SkipTriggerRequirement=True # Temporally disabled (2015-05-28)
-SkipTriggerRequirement= (globalflags.DataSource()=='geant4') # or (jobproperties.Beam.energy()==4000000.0))
+SkipTriggerRequirement= DerivationFrameworkIsMonteCarlo
 # no need to apply trigger requirements on MC 
 print "HIGG2D5.py SkipTriggerRequirement", SkipTriggerRequirement
 TriggerJPSI= []
@@ -576,7 +573,7 @@ HIGG2D5SlimmingHelper.SmartCollections = ["Electrons",
 
 HIGG2D5SlimmingHelper.ExtraVariables = HIGG2D5ExtraContent
 HIGG2D5SlimmingHelper.AllVariables = HIGG2D5ExtraContainers
-if globalflags.DataSource()=='geant4':
+if DerivationFrameworkIsMonteCarlo:
     HIGG2D5SlimmingHelper.ExtraVariables += HIGG2D5ExtraContentTruth
     HIGG2D5SlimmingHelper.AllVariables += HIGG2D5ExtraContainersTruth
 
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG3D1.py b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG3D1.py
index e1ae659dd6d75535fad7395a4bfab50318761453..3978a963a2658e854a9ee224112633cc8b1fe8bc 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG3D1.py
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/share/HIGG3D1.py
@@ -6,6 +6,7 @@
 from DerivationFrameworkCore.DerivationFrameworkMaster import *
 from DerivationFrameworkInDet.InDetCommon import *
 from DerivationFrameworkJetEtMiss.JetCommon import *
+from DerivationFrameworkJetEtMiss.ExtendedJetCommon import *
 from DerivationFrameworkJetEtMiss.METCommon import *
 from DerivationFrameworkEGamma.EGammaCommon import *
 from DerivationFrameworkMuons.MuonsCommon import *
@@ -168,26 +169,15 @@ higg3d1Seq = CfgMgr.AthSequencer("HIGG3d1Sequence")
 #          ghostArea = 0 , ptmin = 2000, ptminFilter = 7000,
 #          variableRMinRadius = 0.02, variableRMassScale = 30000, calibOpt = "none")
 
-#====================================================================
-# Add AntiKt2PV0TrackJets and AntiKt4PV0TrackJets
-#====================================================================
-
-#if not "HIGG3D1Jets" in OutputJets:
-    #OutputJets["HIGG3D1Jets"] = []
-    #AntiKt4PV0TrackJets
-    #addStandardJets("AntiKt", 0.4, "PV0Track", 2000, mods="track_ungroomed", algseq=higg3d1Seq, outputGroup="HIGG3D1Jets")
-    #AntiKt2PV0TrackJets
-    #addStandardJets("AntiKt", 0.2, "PV0Track", 2000, mods="track_ungroomed", algseq=higg3d1Seq, outputGroup="HIGG3D1Jets")
-
 #====================================================================
 # RESTORE JET COLLECTIONS REMOVED BETWEEN r20 AND r21
 #====================================================================
 
-from DerivationFrameworkJetEtMiss.ExtendedJetCommon import replaceAODReducedJets
 reducedJetList = ["AntiKt2PV0TrackJets","AntiKt4PV0TrackJets"]
+if jetFlags.useTruth:
+   reducedJetList += ["AntiKt4TruthJets", "AntiKt4TruthWZJets"]
 replaceAODReducedJets(reducedJetList, higg3d1Seq,"HIGG3D1Jets")
 
-
 #===================================================================
 # Run b-tagging
 #===================================================================
diff --git a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx
index 602e8bde9354d56b541e63b9ce4474e93672609b..2605c2b2928a9996ad74328e57a19a3a1f4fc791 100644
--- a/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx
+++ b/PhysicsAnalysis/DerivationFramework/DerivationFrameworkHiggs/src/TruthCategoriesDecorator.cxx
@@ -190,7 +190,7 @@ namespace DerivationFramework {
     }
 
     // classify event according to simplified template cross section
-    HTXS::HiggsClassification htxs =  m_higgsTruthCatTool->getHiggsTruthCategoryObject(hepmc_evts[0]);
+    HTXS::HiggsClassification htxs =  m_higgsTruthCatTool->getHiggsTruthCategoryObject(hepmc_evts[0],prodMode);
     
     // Decorate the enums
     eventInfo->auxdecor<int>("HTXS_prodMode")   = (int)htxs.prodMode;
@@ -200,9 +200,8 @@ namespace DerivationFramework {
     eventInfo->auxdecor<int>("HTXS_Stage1_Category_pTjet25") = (int)htxs.stage1_cat_pTjet25GeV;
     eventInfo->auxdecor<int>("HTXS_Stage1_Category_pTjet30") = (int)htxs.stage1_cat_pTjet30GeV;
 
-    // JRC - NEXT TWO LINES TEMPORARILY COMMENTED IN GIT TRANSITION
-    //eventInfo->auxdecor<int>("HTXS_Stage1_FineIndex_pTjet30") = HTXSstage1_to_HTXSstage1FineIndex(htxs,th_type);
-    //eventInfo->auxdecor<int>("HTXS_Stage1_FineIndex_pTjet25") = HTXSstage1_to_HTXSstage1FineIndex(htxs,th_type,true);
+    eventInfo->auxdecor<int>("HTXS_Stage1_FineIndex_pTjet30") = HTXSstage1_to_HTXSstage1FineIndex(htxs,th_type);
+    eventInfo->auxdecor<int>("HTXS_Stage1_FineIndex_pTjet25") = HTXSstage1_to_HTXSstage1FineIndex(htxs,th_type,true);
 
     eventInfo->auxdecor<int>("HTXS_Njets_pTjet25")  = (int)htxs.jets25.size();
     eventInfo->auxdecor<int>("HTXS_Njets_pTjet30")  = (int)htxs.jets30.size();