From 076f790d18359cae10da98f06404f812df06be5a Mon Sep 17 00:00:00 2001
From: Sebastien Binet <sebastien.binet@cern.ch>
Date: Fri, 26 Aug 2011 08:20:06 +0200
Subject: [PATCH] migration to HepMC v2.06: removed the IO_xxx classes
(McParticleUtils-00-07-00)
---
.../McParticleUtils/DecayParser.h | 134 ++++
.../McParticleUtils/McParticleUtils/McUtils.h | 32 +
.../McParticleUtils/McVtxFilter.h | 295 +++++++
.../McParticleUtils/cmt/requirements | 44 +
.../McParticleUtils/python/DecayParser.py | 52 ++
.../McParticleUtils/python/__init__.py | 3 +
.../McParticleUtils/src/DecayParser.cxx | 419 ++++++++++
.../McParticleUtils/src/McUtils.cxx | 104 +++
.../McParticleUtils/src/McVtxFilter.cxx | 508 ++++++++++++
.../McParticleUtils/test/McParticleUtils.sh | 13 +
.../test/McVtxFilterTest_CppUnit.cxx | 752 ++++++++++++++++++
11 files changed, 2356 insertions(+)
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/DecayParser.h
create mode 100644 PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McUtils.h
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McVtxFilter.h
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/cmt/requirements
create mode 100644 PhysicsAnalysis/TruthParticleID/McParticleUtils/python/DecayParser.py
create mode 100644 PhysicsAnalysis/TruthParticleID/McParticleUtils/python/__init__.py
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/src/DecayParser.cxx
create mode 100644 PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McUtils.cxx
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McVtxFilter.cxx
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McParticleUtils.sh
create mode 100755 PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McVtxFilterTest_CppUnit.cxx
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/DecayParser.h b/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/DecayParser.h
new file mode 100755
index 000000000000..c61858937508
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/DecayParser.h
@@ -0,0 +1,134 @@
+/////////////////// -*- C++ -*- ////////////////////////////////////
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// DecayParser.h
+// Header file for class Object
+// Author: S.Binet<binet@cern.ch>
+///////////////////////////////////////////////////////////////////
+#ifndef MCPARTICLEUTILS_DECAYPARSER_H
+#define MCPARTICLEUTILS_DECAYPARSER_H
+
+/** DecayParser is the class which parses and models decay patterns.
+ * You can describe <i>via</i> a std::string the decay pattern you are
+ * looking for.
+ * <i>Eg:</i> you just have to say "Z0[e+;e-]" to create a pattern which
+ * can then be used to filter for vertices where such a decay occurs.
+ */
+
+// STL includes
+#include <string>
+#include <vector>
+
+// EventKernel includes
+#include "EventKernel/PdtPdg.h"
+
+// fwd declare
+struct _object;
+typedef _object PyObject;
+
+namespace McUtils {
+ typedef std::vector<std::string> Strings;
+}
+
+class DecayParser
+{
+ ///////////////////////////////////////////////////////////////////
+ // Public methods:
+ ///////////////////////////////////////////////////////////////////
+ public:
+
+ /// Constructor with parameters:
+ DecayParser( const std::string& cmd );
+
+ /// Destructor:
+ virtual ~DecayParser();
+
+ ///////////////////////////////////////////////////////////////////
+ // Const methods:
+ ///////////////////////////////////////////////////////////////////
+ const std::vector<McUtils::Strings>& getParents() const;
+ const std::vector<McUtils::Strings>& getChildren() const;
+
+ void dump() const;
+
+ PDG::pidType pdgId( const std::string& pdgIdString ) const;
+
+ ///////////////////////////////////////////////////////////////////
+ // Non-const methods:
+ ///////////////////////////////////////////////////////////////////
+ void parse( const std::string& cmd );
+
+
+ ///////////////////////////////////////////////////////////////////
+ // Protected methods:
+ ///////////////////////////////////////////////////////////////////
+ protected:
+
+ /// Default constructor:
+ DecayParser(); //> not implemented
+
+ /// Copy constructor:
+ DecayParser( const DecayParser& rhs ); //> not implemented
+
+ /// Assignment operator:
+ DecayParser &operator=(const DecayParser &obj); //> not implemented
+
+ ///////////////////////////////////////////////////////////////////
+ // Protected data:
+ ///////////////////////////////////////////////////////////////////
+ protected:
+
+ /** Print the content of a vector of McUtils::Strings to std::cout
+ */
+ void printMcUtilsStrings( const std::vector<McUtils::Strings>& list ) const;
+
+ ///////////////////////////////////////////////////////////////////
+ // Protected data:
+ ///////////////////////////////////////////////////////////////////
+ protected:
+
+ /** python function to parse the input string modeling the decay pattern
+ * to look for.
+ */
+ PyObject *m_parseFct;
+
+ /** List of parents : each slot of the vector is a list of candidates
+ * So one could have something like : [ ["22","23"] ] to model a decay
+ * of a gamma or a Z0
+ */
+ std::vector<McUtils::Strings> m_parents;
+
+ /** List of children : each slot of the vector is a list of candidates
+ * So one could have something like : [ ["11","13"], ["-11","-13" ] ]
+ * to model a decay into a pair of electrons or muons
+ */
+ std::vector<McUtils::Strings> m_children;
+
+};
+
+/// I/O operators
+//////////////////////
+//MsgStream& operator<<( MsgStream & msg, const DecayParser &obj );
+
+///////////////////////////////////////////////////////////////////
+/// Inline methods:
+///////////////////////////////////////////////////////////////////
+
+inline
+const std::vector<McUtils::Strings>&
+DecayParser::getParents() const
+{
+ return m_parents;
+}
+
+inline
+const std::vector<McUtils::Strings>&
+DecayParser::getChildren() const
+{
+ return m_children;
+}
+
+#endif //> MCPARTICLEUTILS_DECAYPARSER_H
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McUtils.h b/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McUtils.h
new file mode 100644
index 000000000000..0b31016b43ef
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McUtils.h
@@ -0,0 +1,32 @@
+///////////////////////// -*- C++ -*- /////////////////////////////
+
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// McUtils.h
+// Header file for various Mc utils
+// Author: S.Binet<binet@cern.ch>
+///////////////////////////////////////////////////////////////////
+#ifndef MCPARTICLEUTILS_MCUTILS_H
+#define MCPARTICLEUTILS_MCUTILS_H
+
+// STL includes
+
+// HepMC includes
+
+// Forward declaration
+namespace HepPDT { class ParticleDataTable; }
+
+namespace McUtils {
+
+ /** @brief Helper function to get the charge of a particle given its PDG Id.
+ * Also handles (some) particles which are not in the PDGTABLE.MeV file
+ * This function internally uses HepPDT.
+ */
+ double chargeFromPdgId( const int pdgId,
+ const HepPDT::ParticleDataTable* pdt );
+
+} //> end namespace McUtils
+
+#endif //> MCPARTICLEUTILS_MCUTILS_H
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McVtxFilter.h b/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McVtxFilter.h
new file mode 100755
index 000000000000..3c1687e879a9
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/McParticleUtils/McVtxFilter.h
@@ -0,0 +1,295 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// McVtxFilter.h
+// Header file for class McVtxFilter
+// Author: S.Binet<binet@cern.ch>
+///////////////////////////////////////////////////////////////////
+#ifndef MCPARTICLEUTILS_MCVTXFILTER_H
+#define MCPARTICLEUTILS_MCVTXFILTER_H
+
+/** McVtxFilter allows to select HepMC::GenParticle based on the decay pattern.
+ * Hopefully one will be able to also select the particles by pt (in a coming
+ * patch).
+ * It uses internally the DecayPattern class to build the filtering decision
+ * The basic <i>modus operandi</i> is to loop on HepMC::GenVertex (of a given
+ * HepMC::GenEvent) to select the ones which are of interest (and fulfill
+ * the decay pattern).
+ */
+
+// STL includes
+#include <sstream>
+#include <ostream>
+
+// DataModel
+#include "DataModel/DataVector.h"
+
+// Framework includes
+#include "AthenaKernel/getMessageSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/GaudiException.h"
+
+// HepMC includes
+#include "HepMC/GenEvent.h"
+#include "HepMC/GenParticle.h"
+#include "HepMC/GenVertex.h"
+
+// EventKernel includes
+#include "EventKernel/PdtPdg.h"
+
+// AnalysisUtils includes
+#include "AnalysisUtils/IFilterCuts.h"
+#include "AnalysisUtils/ParticleCandidateList.h"
+
+// McParticleUtils includes
+#include "McParticleUtils/DecayParser.h"
+
+// Forward declaration
+
+class McVtxFilter : virtual public IFilterCuts
+{
+
+ ///////////////////////////////////////////////////////////////////
+ // Public methods:
+ ///////////////////////////////////////////////////////////////////
+ public:
+
+ /** Default constructor:
+ */
+ McVtxFilter();
+
+ /** Copy constructor:
+ */
+ McVtxFilter( const McVtxFilter& rhs );
+
+ /** Constructor with parameters:
+ */
+ McVtxFilter( const std::string& decayPattern,
+ const bool matchSign = false,
+ const bool matchBranches = false );
+
+ /** Constructor with parameters:
+ */
+ McVtxFilter( const bool matchSign, const bool matchBranches );
+
+ /** Destructor:
+ */
+ virtual ~McVtxFilter();
+
+ /** Assignment operator:
+ */
+ McVtxFilter &operator=(const McVtxFilter &rhs);
+
+ ///////////////////////////////////////////////////////////////////
+ // Const methods:
+ ///////////////////////////////////////////////////////////////////
+
+ /** Main filter method. This method takes the decision whether or not
+ * the filter has to accept the vertex.
+ * It checks if the vertex has the required number of ingoing particles
+ * the required number of outgoing particles. Then for each branch (parent
+ * child) it checks if the particle can match an asked for particle.
+ */
+ virtual bool isAccepted( const HepMC::GenVertex * vtx ) const;
+
+ /// Dump to std::ostream (default = std::cout) the decay pattern the
+ /// filter one is looking for
+ void dump( std::ostream& out = std::cout ) const;
+
+ /** returns true if the vertex filter requires a parent AND a child list
+ * <i>ie:</i> returns false if the filter looks only for a lonely particle
+ * (eg: a photon or an electron, see TruthParticleAlgs for a concrete example)
+ */
+ bool isFullVtx() const;
+
+ /** return true if the filter discriminates the sign of the particles
+ */
+ bool matchSign() const;
+
+ /** return true if the filter strictly matchs the number of output particles
+ * eg: a vertex t->Wbgg won't be selected if matchBranches == true and
+ * decay pattern = "#id 6[24,5]"
+ */
+ bool matchBranches() const;
+
+ /** return the decay pattern used by this filter
+ */
+ const std::string& decayPattern() const;
+
+ /** return the list of particle candidates for the input particles (parent)
+ */
+ const DataVector<const ParticleCandidateList>& parentList() const;
+
+ /** return the list of particle candidates for the output particles (child)
+ */
+ const DataVector<const ParticleCandidateList>& childList() const;
+
+ ///////////////////////////////////////////////////////////////////
+ // Non-const methods:
+ ///////////////////////////////////////////////////////////////////
+
+ /** Set filter cut options
+ */
+ virtual void setFilter( const IFilterCuts * filter );
+
+ /**Set if one wants to discriminate the sign of the particles
+ */
+ void setMatchSign( const bool matchSign );
+
+ /**Set if one wants to strictly match the number of output particles
+ */
+ void setMatchBranches( const bool matchVtxBranches );
+
+ /**Set the decay pattern to look for : it correctly setups the parent
+ * and children according to the LaTeX-like string which is given as
+ * argument
+ */
+ void setDecayPattern( const std::string& decayPattern );
+
+ ///////////////////////////////////////////////////////////////////
+ // Protected methods:
+ ///////////////////////////////////////////////////////////////////
+ protected:
+
+ /**Check if the parent branch fulfills the requirements
+ */
+ bool checkParentBranch( const HepMC::GenVertex * vtx ) const;
+
+ /**Check if the child branch fulfills the requirements
+ */
+ bool checkChildBranch ( const HepMC::GenVertex * vtx ) const;
+
+ /// Check the branches for the special case of a 1->2 body decay
+ /// (there is some room for optimisation in that particular case)
+ bool checkTwoBodyDecay( const HepMC::GenVertex * vtx ) const;
+
+ /// Populate vtx by HepPDT::ParticleData
+ //void addParent( const HepPDT::ParticleData& pdt );
+ //void addChild ( const HepPDT::ParticleData& pdt );
+
+ /**Populate in-going particles' vtx by particle identification
+ */
+ //void addParent( const PDG::pidType& parentPID );
+ /**Populate out-going particles' vtx by particle identification
+ */
+ //void addChild ( const PDG::pidType& childPID );
+
+ /// Populate vtx by particle name
+ //void addParent( const std::string& parentName );
+ //void addChild ( const std::string& childName );
+
+ /** Populate in-going particles' vtx list by group name
+ */
+ //void addParentList( const std::string& parentClass );
+
+ /** Populate out-going particles' vtx list by group name
+ */
+ //void addChildList ( const std::string& childClass );
+
+ /// Populate parents' vtx list by copy
+ /// McVtxFilter takes ownership of the given list
+ //void addParentList( const ParticleCandidateList* parentList );
+ /// Populate children's vtx list by copy :
+ /// McVtxFilter takes ownership of the given list
+ //void addChildList ( const ParticleCandidateList* childList );
+
+ ///////////////////////////////////////////////////////////////////
+ // Protected data:
+ ///////////////////////////////////////////////////////////////////
+ protected:
+
+ /// MsgStream instance (a std::cout like with print-out levels)
+ mutable MsgStream m_msg;
+
+ /// Discrimination between particles and conjugates
+ bool m_matchSign;
+
+ /// Tell if one wants a strict vertex branches matching.
+ /// This is to cope with radiations in the output branches
+ /// If true then a decay t->Wbggg will not be kept if one asked for
+ /// a top into W+b decay...
+ /// Default : false
+ bool m_matchBranches;
+
+ /// The decay pattern one looks for.
+ /// Ex: 23 -> -11 + 11 to model a Z0 -> e+ e- decay
+ std::string m_decayPattern;
+
+ /// The list of the matching particles Ids : each item in DataVector
+ /// stands for an entering branch to the vertex
+ DataVector<const ParticleCandidateList> m_parentList;
+
+ /// The list of the matching particles Ids : each item in DataVector
+ /// stands for an out-going branch to the vertex
+ DataVector<const ParticleCandidateList> m_childList;
+
+};
+
+/// I/O operators
+//////////////////////
+MsgStream& operator<<( MsgStream & msg, const McVtxFilter &obj );
+std::ostream& operator<<( std::ostream& out, const McVtxFilter &obj );
+
+///////////////////////////////////////////////////////////////////
+// Inline methods:
+///////////////////////////////////////////////////////////////////
+
+inline McVtxFilter::~McVtxFilter(){}
+
+///////////////////////////////////////////////////////////////////
+// Const methods:
+///////////////////////////////////////////////////////////////////
+
+inline bool McVtxFilter::isFullVtx() const
+{
+ if ( m_parentList.empty() || m_childList.empty() ) {
+ return false;
+ } else {
+ return true;
+ }
+}
+
+inline bool McVtxFilter::matchSign() const
+{
+ return m_matchSign;
+}
+
+inline bool McVtxFilter::matchBranches() const
+{
+ return m_matchBranches;
+}
+
+inline const std::string& McVtxFilter::decayPattern() const
+{
+ return m_decayPattern;
+}
+
+inline
+const DataVector<const ParticleCandidateList>& McVtxFilter::parentList() const
+{
+ return m_parentList;
+}
+
+inline
+const DataVector<const ParticleCandidateList>& McVtxFilter::childList() const
+{
+ return m_childList;
+}
+
+///////////////////////////////////////////////////////////////////
+// Non-const methods:
+///////////////////////////////////////////////////////////////////
+
+inline void McVtxFilter::setMatchSign( const bool matchSign )
+{
+ m_matchSign = matchSign;
+}
+
+inline void McVtxFilter::setMatchBranches( const bool matchBranches )
+{
+ m_matchBranches = matchBranches;
+}
+
+#endif //> MCPARTICLEUTILS_MCVTXFILTER_H
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/cmt/requirements b/PhysicsAnalysis/TruthParticleID/McParticleUtils/cmt/requirements
new file mode 100755
index 000000000000..c7e77212516a
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/cmt/requirements
@@ -0,0 +1,44 @@
+package McParticleUtils
+
+author Sebastien Binet <binet@cern.ch>
+
+use AtlasPolicy AtlasPolicy-*
+
+use GaudiInterface GaudiInterface-* External
+use AtlasHepMC AtlasHepMC-* External
+
+use EventKernel EventKernel-* Event
+
+use DataModel DataModel-* Control
+use AthenaKernel AthenaKernel-* Control
+use DataModel DataModel-* Control
+
+use AnalysisUtils AnalysisUtils-* PhysicsAnalysis/AnalysisCommon
+
+## private uses
+private
+use AtlasBoost AtlasBoost-* External
+use AtlasPython AtlasPython-* External
+use AtlasCLHEP AtlasCLHEP-* External
+use HepPDT v* LCG_Interfaces
+end_private
+
+branches src McParticleUtils doc test run share
+
+library McParticleUtils *.cxx
+
+apply_pattern installed_library
+
+#########################################################
+# Unit tests: CppUnit tests for McParticleUtils classes #
+#########################################################
+private
+use TestPolicy TestPolicy-* -no_auto_imports
+use TestTools TestTools-* AtlasTest -no_auto_imports
+use StoreGate StoreGate-* Control
+use AtlasCppUnit AtlasCppUnit-* External -no_auto_imports
+apply_pattern CppUnit name=cppUnit_mcVtxFilterTest files="-s=${McParticleUtils_root}/test McVtxFilterTest_CppUnit.cxx"
+macro_append HepPDT_linkopts " -lHepPID"
+
+apply_pattern declare_python_modules files="*.py"
+end_private
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/python/DecayParser.py b/PhysicsAnalysis/TruthParticleID/McParticleUtils/python/DecayParser.py
new file mode 100644
index 000000000000..e5bc7e2ee8a0
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/python/DecayParser.py
@@ -0,0 +1,52 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+# @file: McParticleUtils/python/DecayParser.py
+# @purpose: implement the parsing/tokenization of decay-patterns' strings
+# @author: Sebastien Binet <binet@cern.ch>
+
+import re
+
+_slot_separator = "+"
+_candidate_sep = "|"
+_decay_arrow = "->"
+_wild_card = "*"
+
+def py_parse (cmd):
+ # remove all spaces to ease further parsing
+ cmd = re.sub(pattern=" ",
+ repl="",
+ string=cmd)
+## print "=> [%s]"%cmd
+
+ if cmd == _decay_arrow:
+ # special case of a single arrow without any parent nor child:
+ # this decay pattern will select every single vertex
+ return (0, None, None)
+ buf = cmd.split (_decay_arrow)
+ assert (len(buf)==2)
+
+## print "==> buf:",buf
+ parents = process_block (buf[0])
+ children= process_block (buf[1])
+
+ return (0, parents, children)
+
+def process_block (cmd):
+ if cmd == "":
+ return None
+ candidates = cmd.split (_slot_separator)
+ return [list(set(c.split(_candidate_sep))) for c in candidates]
+
+
+if __name__ == "__main__":
+ print ":"*80
+ print "::: tests..."
+ for cmd in ("-1|-2|-3|-4|-5|-6|21 + 1|2|3|4|5|6|21 -> ",
+ "-> 91|92|94"
+ ):
+ print "::: testing [%s]..." % cmd
+ print "::: result:",py_parse (cmd)
+
+ print "::: bye."
+ print ":"*80
+
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/python/__init__.py b/PhysicsAnalysis/TruthParticleID/McParticleUtils/python/__init__.py
new file mode 100644
index 000000000000..8056e12aa60d
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/python/__init__.py
@@ -0,0 +1,3 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+# hook for McParticleUtils python package
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/DecayParser.cxx b/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/DecayParser.cxx
new file mode 100755
index 000000000000..4115e4d6e484
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/DecayParser.cxx
@@ -0,0 +1,419 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// DecayParser.cxx
+// Implementation file for class DecayParser
+// Author: S.Binet<binet@cern.ch>
+///////////////////////////////////////////////////////////////////
+
+// Python includes
+#include "Python.h"
+
+// STL includes
+#include <iostream>
+#include <list>
+#include <stdexcept>
+#include <sstream>
+
+// boost includes
+// #include <boost/algorithm/string/split.hpp>
+// #include <boost/algorithm/string/erase.hpp>
+// #include <boost/algorithm/string/regex.hpp>
+// #include <boost/algorithm/string/trim.hpp>
+// #include <boost/algorithm/string/compare.hpp>
+
+// McParticleUtils includes
+#include "McParticleUtils/DecayParser.h"
+
+// namespace balg = boost::algorithm;
+
+namespace {
+ PyObject *fetch_py_parse_fct();
+ bool py_to_cpp (PyObject* candidates,
+ std::vector<McUtils::Strings>& parsed);
+}
+
+///////////////////////////////////////////////////////////////////
+// Static data
+///////////////////////////////////////////////////////////////////
+
+// const std::string DecayParser::m_slotSep = "+" ;
+// const boost::regex DecayParser::m_candidateSep ( "\\|" );
+// const boost::regex DecayParser::m_decayArrow ( "->" );
+// const boost::regex DecayParser::m_wildCardParticle ( "\\*" );
+
+///////////////////////////////////////////////////////////////////
+/// Public methods:
+///////////////////////////////////////////////////////////////////
+
+/// Constructors
+////////////////
+// DecayParser::DecayParser() :
+// m_parseFct (0),
+// m_parents ( ),
+// m_children ( )
+// {}
+
+DecayParser::DecayParser( const std::string& cmd ) :
+ m_parseFct (0),
+ m_parents ( ),
+ m_children ( )
+{
+ m_parseFct = ::fetch_py_parse_fct();
+ parse(cmd);
+}
+
+// DecayParser::DecayParser( const DecayParser& rhs ) :
+// m_parseFct (0),
+// m_parents ( rhs.m_parents ),
+// m_children ( rhs.m_children )
+// {}
+
+/// Destructor
+///////////////
+DecayParser::~DecayParser()
+{
+ Py_XDECREF (m_parseFct);
+}
+
+///////////////////////////////////////////////////////////////////
+/// Const methods:
+///////////////////////////////////////////////////////////////////
+void DecayParser::dump() const
+{
+ std::cout << "--- Parents ---" << std::endl;
+ printMcUtilsStrings( m_parents );
+
+ std::cout << "--- Children ---" << std::endl;
+ printMcUtilsStrings( m_children );
+}
+
+PDG::pidType DecayParser::pdgId( const std::string& pdgIdString ) const
+{
+ PDG::pidType pdgID = PDG::null;
+
+ /// We have been setup to directly use PdgIds
+ /// so here we just convert a PDG-ID string into a bare number
+ /// Note: if we were in 1990's we could have used atoi(char*)
+ /// but, hey, let us use ANSI-C++ : (need to go through a temporary integer
+ /// because explicit cast of PDG::pidType into integer does not work
+ /// with int std::stringstream::operator>>() ==> Why ? //FIXME
+ int iPDG = 0;
+ std::stringstream( pdgIdString ) >> iPDG;
+ pdgID = static_cast<PDG::pidType>(iPDG);
+
+ return pdgID;
+}
+///////////////////////////////////////////////////////////////////
+/// Non-const methods:
+///////////////////////////////////////////////////////////////////
+void DecayParser::parse( const std::string& inputCmd )
+{
+ if ( inputCmd.empty() ) {
+ return;
+ }
+
+ // Reset the parents and children lists
+ m_parents.clear();
+ m_children.clear();
+
+// /// Test if there is an arrow to enforce our syntax
+// if ( inputCmd.find (m_decayArrow.str()) == std::string::npos ) {
+// // Malformed command
+// std::string error = "";
+// error += "Malformed command : <"+inputCmd+"> !!\n";
+// error += "You should enter an arrow... ex: 23 -> -13 +13 or 22-> 11 + -11";
+// throw std::invalid_argument(error);
+// }
+
+// std::string cmd = inputCmd;
+
+// /// remove spaces to ease further parsing
+// rmSpaces(cmd);
+
+// if ( cmd == m_decayArrow.str() ) {
+// }
+
+ // real parsing takes place now.
+ PyObject *res = PyObject_CallFunction (m_parseFct,
+ (char*)"s",
+ inputCmd.c_str());
+ if (!res) {
+ Py_XDECREF (res);
+ std::string error = "problem while parsing command [" + inputCmd +"]";
+ throw std::runtime_error (error);
+ }
+
+ if (!PyTuple_Check (res)) {
+ Py_DECREF (res);
+ std::string error = "expected a python tuple";
+ throw std::runtime_error (error);
+ }
+
+ if (PyTuple_GET_SIZE (res) != 3) {
+ Py_DECREF (res);
+ std::string error = "expected a python tuple of size 3";
+ throw std::runtime_error (error);
+ }
+
+ PyObject *sc = PyTuple_GET_ITEM (res, 0);
+ Py_XINCREF (sc);
+ if (!sc || !PyInt_Check (sc)) {
+ Py_XDECREF (sc);
+ Py_DECREF (res);
+ std::string error = "corrupted return code";
+ throw std::runtime_error (error);
+ }
+
+ Py_ssize_t status = PyInt_AsSsize_t (sc);
+ if (status != 0) {
+ Py_DECREF (sc);
+ Py_DECREF (res);
+ std::string error = "failed to parse command ["+inputCmd+"]";
+ throw std::runtime_error (error);
+ }
+ Py_DECREF (sc);
+
+ PyObject *parents = PyTuple_GET_ITEM (res, 1);
+ Py_XINCREF (parents);
+ if (!parents) {
+ Py_DECREF (res);
+ std::string error = "corrupted parents' list";
+ throw std::runtime_error (error);
+ }
+
+ PyObject *children= PyTuple_GET_ITEM (res, 2);
+ Py_XINCREF (children);
+ if (!children) {
+ Py_DECREF (parents);
+ Py_DECREF (res);
+ std::string error = "corrupted children' list";
+ throw std::runtime_error (error);
+ }
+ Py_DECREF (res);
+
+ if (parents==Py_None && children==Py_None) {
+ // special case of a single arrow without any parent nor child :
+ // this decay pattern will select every single vertex
+ Py_DECREF (parents);
+ Py_DECREF (children);
+ return;
+ }
+
+ if (!py_to_cpp (parents, m_parents)) {
+ Py_DECREF (parents);
+ Py_DECREF (children);
+ std::string error = "could not translate parents' list";
+ throw std::runtime_error (error);
+ }
+
+ if (!py_to_cpp (children, m_children)) {
+ Py_DECREF (parents);
+ Py_DECREF (children);
+ std::string error = "could not translate children' list";
+ throw std::runtime_error (error);
+ }
+
+ return;
+}
+
+
+///////////////////////////////////////////////////////////////////
+/// Protected methods:
+///////////////////////////////////////////////////////////////////
+
+// void
+// DecayParser::parseCandidates( const std::string& inString,
+// McUtils::Strings& outList ) const
+// {
+// McUtils::Strings tmpList;
+// balg::split_regex (tmpList, inString, m_candidateSep);
+// outList.clear();
+// outList.reserve (tmpList.size());
+// const McUtils::Strings::const_iterator tmpListEnd = tmpList.end();
+// for ( McUtils::Strings::const_iterator i = tmpList.begin();
+// i != tmpListEnd;
+// ++i ) {
+// if ( !i->empty() ) {
+// outList.push_back( *i );
+// }
+// }
+// return;
+// }
+
+void
+DecayParser::printMcUtilsStrings( const std::vector<McUtils::Strings>& list ) const
+{
+ unsigned int iSlot = 0;
+ for( std::vector<McUtils::Strings>::const_iterator itr = list.begin();
+ itr != list.end();
+ ++itr,++iSlot ) {
+ std::stringstream iSlotStr;
+ iSlotStr << iSlot;
+ const McUtils::Strings::const_iterator candEnd = itr->end();
+ std::cout << "slot #" << iSlotStr.str() << ": candidates= [ ";
+ for( McUtils::Strings::const_iterator candidate = itr->begin();
+ candidate != candEnd;
+ ++candidate ) {
+ std::cout << *candidate;
+ if ( candidate+1 != candEnd ) {
+ std::cout << " | ";
+ }
+ }
+ std::cout << " ]" << std::endl;
+ }
+ return;
+}
+
+// McUtils::Strings
+// DecayParser::tokenizeString( const std::string& stringList ) const
+// {
+// std::cerr << "==DecayParser== tokenizeString (" << stringList << ")...\n";
+// // std::cerr << " ==> slotSep: [" << m_slotSep << "]\n"
+// // << " ==> candSep: [" << m_candidateSep << "]\n"
+// // << " ==> dcyArrw: [" << m_decayArrow << "]\n"
+// // << " ==> wildCrd: [" << m_wildCardParticle << "]\n";
+
+// std::string nameList = stringList;
+// /// Removes all the spaces in the strings : " " -> ""
+// rmSpaces(nameList);
+
+// McUtils::Strings out;
+// balg::split (out, nameList, balg::is_any_of(m_slotSep));
+// return out;
+// }
+
+// void DecayParser::fillToken( const McUtils::Strings& token,
+// std::vector<McUtils::Strings>& tokens ) const
+// {
+// std::cerr << "==DecayParser== fillToken (" << token.size() << ")...\n";
+// for ( McUtils::Strings::const_iterator itr = token.begin();
+// itr != token.end();
+// ++itr ) {
+// //std::cout << "Parsing : <" << *itr << ">..." << std::endl;
+// /// this is the list of matching candidates for each child
+// McUtils::Strings candidates;
+// /// parse the child string to unmangle the possible choices for each
+// /// of the children ( resolves that e-|mu- can be either a muon or
+// /// an electron )
+// parseCandidates( *itr, candidates );
+// if ( !candidates.empty() ) {
+// tokens.push_back( candidates );
+// }
+// }
+// return;
+// }
+
+// void DecayParser::processBlock( const std::string& tokenString,
+// std::vector<McUtils::Strings>& tokens ) const
+// {
+// std::cerr << "==DecayParser== processBlock ("
+// << tokenString <<", tokens=" << tokens.size()
+// << ")...\n";
+
+// // cut the tokenString into slots
+// McUtils::Strings tokenNames = tokenizeString( tokenString );
+
+// // cut each slot into a list of candidates
+// fillToken( tokenNames, tokens );
+// }
+
+///////////////////////////////////////////////////////////////////
+// Operators:
+///////////////////////////////////////////////////////////////////
+DecayParser & DecayParser::operator=(const DecayParser& rhs )
+{
+ if ( this != &rhs ) {
+ m_parents = rhs.m_parents;
+ m_children = rhs.m_children;
+ }
+ return *this;
+}
+
+namespace {
+
+PyObject*
+fetch_py_parse_fct()
+{
+ // need to ensure the python interpreter has been initialized...
+ if (!Py_IsInitialized()) {
+ Py_Initialize();
+ }
+
+ const std::string n = "McParticleUtils.DecayParser";
+ PyObject *module = PyImport_ImportModule (const_cast<char*>(n.c_str()));
+ if (!module || !PyModule_Check (module)) {
+ Py_XDECREF (module);
+ std::string error = "could not import module ["+n+"]";
+ throw std::runtime_error (error);
+ }
+
+ const std::string fct_name = "py_parse";
+ PyObject *fct = PyDict_GetItemString (PyModule_GetDict (module),
+ const_cast<char*>(fct_name.c_str()));
+ // borrowed ref.
+ Py_XINCREF (fct);
+ // don't need the module anymore
+ Py_DECREF (module);
+
+ if (!fct || !PyFunction_Check (fct)) {
+ std::string error = "could not get '"+fct_name+"' from module ["+n+"] or not a function";
+ throw std::runtime_error (error);
+ }
+
+ return fct;
+}
+
+bool
+py_to_cpp (PyObject* candidates,
+ std::vector<McUtils::Strings>& parsed)
+{
+ bool all_good = true;
+ if (candidates==Py_None) {
+ // nothing to do
+ return true;
+ }
+
+ if (!PySequence_Check (candidates)) {
+ return false;
+ }
+ Py_ssize_t isz = PySequence_Size (candidates);
+ if (isz==-1) {
+ return false;
+ }
+ parsed.resize (isz);
+
+ for (Py_ssize_t i = 0; i!=isz; ++i) {
+ PyObject *cand = PySequence_GetItem(candidates, i);
+ if (!cand) {
+ return false;
+ }
+ if (!PySequence_Check (cand)) {
+ Py_DECREF (cand);
+ return false;
+ }
+ Py_ssize_t jsz = PySequence_Size (cand);
+ if (jsz==-1) {
+ Py_DECREF (cand);
+ return false;
+ }
+
+ parsed[i].resize(jsz);
+
+ for (Py_ssize_t j = 0; j!=jsz; ++j) {
+ PyObject *pdgid = PySequence_GetItem(cand, j);
+ if (!pdgid) {
+ Py_DECREF (cand);
+ return false;
+ }
+
+ parsed[i][j] = std::string(PyString_AS_STRING(pdgid));
+ }
+
+ Py_DECREF (cand);
+ }
+ return all_good;
+}
+} //> anon-namespace
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McUtils.cxx b/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McUtils.cxx
new file mode 100644
index 000000000000..7670d5a96fa2
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McUtils.cxx
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// McUtils.cxx
+// Implementation file for various Mc helpers and methods
+// Author: S.Binet<binet@cern.ch>
+///////////////////////////////////////////////////////////////////
+
+// STL includes
+#include <cmath>
+#include <stdexcept>
+#include <cstring>
+#include <string>
+
+// boost includes
+#include "boost/array.hpp"
+
+// FrameWork includes
+#include "GaudiKernel/IPartPropSvc.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// CLHEP/HepMC includes
+#include "HepPDT/ParticleData.hh"
+#include "HepPDT/ParticleDataTable.hh"
+
+// McParticleUtils includes
+#include "McParticleUtils/McUtils.h"
+
+namespace {
+
+/** 3*charge for basic pdgId codes -- used to parse unknown id's
+ Fix from Frank for the charge of the MC Truth Particle */
+static const boost::array<int, 100> qcharge =
+ {{+0, -1, +2, -1, +2, -1, +2, -1, +2, +0, // 0- 9
+ +0, -3, +0, -3, +0, -3, +0, -3, +0, +0, // 10-19
+ +0, +0, +0, +3, +0, +0, +0, +0, +0, +0, // 20-29
+ +0, +0, +0, +3, +0, +0, +0, +3, +0, +0, // 30-39
+ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0,
+ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0,
+ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0,
+ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0,
+ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0,
+ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0}};
+
+}
+
+
+double McUtils::chargeFromPdgId( const int pdgId,
+ const HepPDT::ParticleDataTable* pdt )
+{
+ using std::abs;
+ if ( !pdt ) {
+ throw std::invalid_argument( "Null pointer to HepPDT::ParticleDataTable !" );
+ }
+
+ /** FIXME : should use HepPDT::ParticleDataTable to retrieve
+ the charge of the MC particle
+ Problem : how do we handle particles which are not in the
+ PDGTABLE.MeV file ? */
+ const HepPDT::ParticleData* ap = pdt->particle( abs( pdgId ) );
+ if ( ap ) {
+ return ( pdgId < 0 )
+ ? -(ap->charge())
+ : ap->charge();
+ } else if ( 0 == pdgId ) {
+ // we are trying to compute a charge for a non existing particle...
+ //throw std::invalid_argument( "Non existing PDG-id" );
+ return -999.;
+ } else {
+ /** Set charge using PDG convention:
+ id = nnnnijkl
+ i == 0, j == 0: see qcharge[100]
+ i == 0: meson, j kbar quarks l = 2*spin+1
+ i != 0: baryon, i j k quarks l = 2*spin+1
+ Default is 0; */
+ const int idmod = abs(pdgId) % 10000;
+ const int q1 = (idmod/10) % 10;
+ const int q2 = (idmod/100) % 10;
+ const int q3 = (idmod/1000) % 10;
+ double q = 0;
+
+ if (abs(pdgId) >= 1000000000) {
+ // Seems to be a nucleus: 100pppnnn0
+ q = (abs(pdgId) / 10000) % 1000;
+ }
+ else if( idmod < 100 ) {
+ q = qcharge[idmod]/3.;
+ }
+ else if ( idmod < 1000 ) {
+ q = (qcharge[q1]-qcharge[q2])/3.;
+ if ( qcharge[q2] == 2 ) {
+ q *= -1.;
+ }
+ }
+ else if( idmod < 10000 ) {
+ q = (qcharge[q3]+qcharge[q2]+qcharge[q1])/3.;
+ }
+ if (q == 0) q = 0; // Change -0 -> 0.
+ return (pdgId < 0) ? -q : q;
+ }
+}
+
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McVtxFilter.cxx b/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McVtxFilter.cxx
new file mode 100755
index 000000000000..84cae9ef52e9
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/src/McVtxFilter.cxx
@@ -0,0 +1,508 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// McVtxFilter.cxx
+// Implementation file for class McVtxFilter
+// Author: S.Binet<binet@cern.ch>
+///////////////////////////////////////////////////////////////////
+
+
+// STL includes
+
+// Framework includes
+#include "AthenaKernel/getMessageSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/GaudiException.h"
+
+// HepMC includes
+#include "HepMC/GenParticle.h"
+#include "HepMC/GenVertex.h"
+
+// AnalysisUtils includes
+#include "AnalysisUtils/AnalysisPermutation.h"
+
+// McParticleUtils includes
+#include "McParticleUtils/McVtxFilter.h"
+
+///////////////////////////////////////////////////////////////////
+/// Public methods:
+///////////////////////////////////////////////////////////////////
+
+/// Constructors
+////////////////
+McVtxFilter::McVtxFilter() :
+ IFilterCuts(),
+ m_msg ( Athena::getMessageSvc(), "McVtxFilter" ),
+ m_matchSign( false ),
+ m_matchBranches( false ),
+ m_decayPattern( "" ),
+ m_parentList(),
+ m_childList()
+{}
+
+McVtxFilter::McVtxFilter( const std::string& decayPattern,
+ const bool matchSign,
+ const bool matchBranches ) :
+ IFilterCuts(),
+ m_msg ( Athena::getMessageSvc(), "McVtxFilter" ),
+ m_matchSign( matchSign ),
+ m_matchBranches( matchBranches ),
+ m_decayPattern( decayPattern ),
+ m_parentList(),
+ m_childList()
+{
+ setDecayPattern(decayPattern);
+}
+
+McVtxFilter::McVtxFilter( const bool matchSign,
+ const bool matchBranches ) :
+ IFilterCuts(),
+ m_msg ( Athena::getMessageSvc(), "McVtxFilter" ),
+ m_matchSign( matchSign ),
+ m_matchBranches( matchBranches ),
+ m_decayPattern( "" ),
+ m_parentList(),
+ m_childList()
+{}
+
+McVtxFilter::McVtxFilter( const McVtxFilter& rhs ) :
+ IFilterCuts ( rhs ),
+ m_msg ( rhs.m_msg ),
+ m_matchSign ( rhs.m_matchSign ),
+ m_matchBranches( rhs.m_matchBranches ),
+ m_decayPattern ( rhs.m_decayPattern ),
+ m_parentList(),
+ m_childList()
+{
+ // deep copy of the parent branch
+ m_parentList.reserve( rhs.m_parentList.size() );
+ for ( DataVector<const ParticleCandidateList>::const_iterator itr =
+ rhs.m_parentList.begin();
+ itr != rhs.m_parentList.end();
+ ++itr ) {
+ m_parentList.push_back( new ParticleCandidateList(**itr) );
+ }
+
+ // deep copy of the child branch
+ m_childList.reserve( rhs.m_childList.size() );
+ for ( DataVector<const ParticleCandidateList>::const_iterator itr =
+ rhs.m_childList.begin();
+ itr != rhs.m_childList.end();
+ ++itr ) {
+ m_childList.push_back( new ParticleCandidateList(**itr) );
+ }
+}
+
+McVtxFilter & McVtxFilter::operator=(const McVtxFilter &rhs)
+{
+ if ( this != &rhs ) {
+ IFilterCuts::operator=(rhs);
+ m_matchSign = rhs.m_matchSign;
+ m_matchBranches = rhs.m_matchBranches;
+ m_decayPattern = rhs.m_decayPattern;
+
+ // deep copy of the parent branch
+ m_parentList.reserve( rhs.m_parentList.size() );
+ for ( DataVector<const ParticleCandidateList>::const_iterator itr =
+ rhs.m_parentList.begin();
+ itr != rhs.m_parentList.end();
+ ++itr ) {
+ m_parentList.push_back( new ParticleCandidateList(**itr) );
+ }
+
+ // deep copy of the child branch
+ m_childList.reserve( rhs.m_childList.size() );
+ for ( DataVector<const ParticleCandidateList>::const_iterator itr =
+ rhs.m_childList.begin();
+ itr != rhs.m_childList.end();
+ ++itr ) {
+ m_childList.push_back( new ParticleCandidateList(**itr) );
+ }
+ }
+ return *this;
+}
+
+/// Destructor
+///////////////
+
+///////////////////////////////////////////////////////////////////
+/// Const methods:
+///////////////////////////////////////////////////////////////////
+bool McVtxFilter::isAccepted( const HepMC::GenVertex * vtx ) const
+{
+ m_msg << MSG::VERBOSE
+ << "In McVtxFilter::isAccepted(...)" << endreq;
+
+ ////////////////////////////////////////////////////////////////
+ /// First handle special case where the filter has only 1 child
+ /// and no parent : one is looking for a stable particle with
+ /// no end_vertex
+ if ( m_childList.size() == static_cast<unsigned int>( 1 ) &&
+ m_parentList.size() == static_cast<unsigned int>( 0 ) &&
+ vtx->particles_out_size() == static_cast<unsigned int>( 1 ) ) {
+ const HepMC::GenParticle * part = *(vtx->particles_out_const_begin());
+ const ParticleCandidateList * item = *( m_childList.begin() );
+ if ( item->hasInList( static_cast<PDG::pidType>(part->pdg_id()),
+ m_matchSign ) ) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ /// Skip vertices which don't match the number of branches
+ /// we are looking for : check if number of parents
+ /// and of children is OK
+ if ( !m_matchBranches ) {
+
+ //
+ // Test only if there is enough *output* branches to match for
+ //
+ if ( static_cast<unsigned int>(vtx->particles_in_size()) < m_parentList.size() ||
+ static_cast<unsigned int>(vtx->particles_out_size()) < m_childList.size() ) {
+ return false;
+ }
+
+ } else {
+ //
+ // Strict match of output branches
+ //
+ if ( static_cast<unsigned int>(vtx->particles_in_size()) != m_parentList.size() ||
+ static_cast<unsigned int>(vtx->particles_out_size()) != m_childList.size() ) {
+ return false;
+ }
+ }
+
+ //////////////////////////////////////////////////////////////
+ /// Then handle the case we are looking for a 1->2 body decay
+ ///
+ if ( m_matchBranches &&
+ m_parentList.size() == static_cast<unsigned int>(1) &&
+ m_childList.size() == static_cast<unsigned int>(2) &&
+ vtx->particles_out_size() >= 2 ) {
+ return checkTwoBodyDecay( vtx );
+ } //> two-body decay
+
+
+ ////////////////////////////////
+ /// Handle other generic cases
+ ///
+ m_msg << MSG::VERBOSE << "trying checkParentBranch(...)" << endreq;
+ if ( checkParentBranch( vtx ) == false ) return false;
+
+ m_msg << MSG::VERBOSE << "trying checkChildBranch(...)" << endreq;
+ if ( checkChildBranch ( vtx ) == false ) return false;
+
+ m_msg << MSG::VERBOSE << "McVtxFilter::isAccepted(...) => DONE" << endreq;
+ return true;
+}
+
+void McVtxFilter::dump( std::ostream& out ) const
+{
+ out << ">>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<" << std::endl;
+ out << ">>> Parents : " << std::endl;
+ for ( DataVector<const ParticleCandidateList>::const_iterator itr = m_parentList.begin();
+ itr != m_parentList.end();
+ ++itr ) {
+ out << "\nList for : " << *itr
+ << " (" << (*itr)->size() << ")" << std::endl;
+ (*itr)->dropList();
+ }
+
+ out << ">>> Children : " << std::endl;
+ for ( DataVector<const ParticleCandidateList>::const_iterator itr = m_childList.begin();
+ itr != m_childList.end();
+ ++itr ) {
+ out << "\nList for : " << *itr
+ << " (" << (*itr)->size() << ")" << std::endl;
+ (*itr)->dropList();
+ }
+ out << ">>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<" << std::endl;
+
+ return;
+}
+
+
+///////////////////////////////////////////////////////////////////
+/// Non-const methods:
+///////////////////////////////////////////////////////////////////
+
+/// Set filter cut options
+void McVtxFilter::setFilter( const IFilterCuts * filter )
+{
+ if ( filter ) {
+ const McVtxFilter * vtxFilter =
+ dynamic_cast<const McVtxFilter*>(filter);
+
+ if ( vtxFilter ) {
+ operator=(*vtxFilter);
+ } else {
+ m_msg << MSG::ERROR
+ << "Can't dynamic_cast " << typeid(filter).name()
+ << " to a McVtxFilter"
+ << endreq;
+ }
+ } //> filter is a valid pointer
+}
+
+/// Set the filter from a decay pattern
+void McVtxFilter::setDecayPattern( const std::string& decayPattern )
+{
+ // clear parent and child candidates
+ m_parentList.clear();
+ m_childList.clear();
+
+ m_decayPattern = decayPattern;
+
+ DecayParser parser( m_decayPattern );
+
+ //std::cout << "Populate parent list" << std::endl;
+ std::vector<std::vector<std::string> > parents = parser.getParents();
+ for(std::vector<std::vector<std::string> >::const_iterator itr = parents.begin();
+ itr != parents.end();
+ ++itr ) {
+ ParticleCandidateList * list = new ParticleCandidateList();
+ const std::vector<std::string>::const_iterator candEnd = itr->end();
+ for( std::vector<std::string>::const_iterator candidate = itr->begin();
+ candidate != candEnd;
+ ++candidate ) {
+ PDG::pidType pdgID = parser.pdgId( *candidate );
+ list->push_back( pdgID );
+ }
+ if ( ! list->empty() ) {
+ m_parentList.push_back( list );
+ } else {
+ delete list;
+ }
+ }
+ //std::cout << "Populate children list" << std::endl;
+ std::vector<std::vector<std::string> > children = parser.getChildren();
+ for(std::vector<std::vector<std::string> >::const_iterator itr = children.begin();
+ itr != children.end();
+ ++itr ) {
+ ParticleCandidateList * list = new ParticleCandidateList();
+ const std::vector<std::string>::const_iterator candEnd = itr->end();
+ for( std::vector<std::string>::const_iterator candidate = itr->begin();
+ candidate != candEnd;
+ ++candidate ) {
+ PDG::pidType pdgID = parser.pdgId( *candidate );
+ list->push_back( pdgID );
+ }
+ if ( ! list->empty() ) {
+ m_childList.push_back( list );
+ } else {
+ delete list;
+ }
+ }
+
+ return;
+}
+
+
+///////////////////////////////////////////////////////////////////
+/// Protected methods:
+///////////////////////////////////////////////////////////////////
+
+bool McVtxFilter::checkParentBranch( const HepMC::GenVertex * vtx ) const
+{
+ m_msg << MSG::VERBOSE << "In checkParentBranch..." << endreq;
+
+ /// Check we aren't in the "any particle" case
+ if ( m_parentList.empty() ) {
+ return true;
+ }
+
+ if ( m_msg.level() <= MSG::VERBOSE ) {
+ vtx->print();
+ }
+
+ /// Check if number of parents is OK
+ if ( static_cast<unsigned int>(vtx->particles_in_size()) < m_parentList.size() ){
+ return false;
+ }
+
+ if ( m_msg.level() <= MSG::VERBOSE ) {
+ m_msg << MSG::VERBOSE
+ << "Number of list of parents : "
+ << m_parentList.size()
+ << endreq;
+ m_parentList.front()->dropList();
+ }
+
+ std::vector<int> parentIds;
+ for ( HepMC::GenVertex::particles_in_const_iterator itrPart = vtx->particles_in_const_begin();
+ itrPart != vtx->particles_in_const_end();
+ ++itrPart ) {
+ parentIds.push_back( (*itrPart)->pdg_id() );
+ }
+
+ AnalysisUtils::Permutation<std::vector<int> > permute( &parentIds,
+ m_parentList.size() );
+ std::vector<int> parents;
+
+ bool accepted = false;
+ while ( permute.get(parents) && !accepted ) {
+ accepted = true;
+ const unsigned int iMax = std::min( m_parentList.size(), parentIds.size() );
+ for ( unsigned int i = 0; i != iMax; ++i ) {
+ const bool hasInList =
+ m_parentList[i]->hasInList( static_cast<PDG::pidType>(parents[i]),
+ m_matchSign );
+ if ( !hasInList ) {
+ // this permutation is not suiting, going to the next one
+ // (if any)
+ accepted = false;
+ break;
+ }
+ }
+
+ // everything went fine for this permutation !
+ if ( accepted ) {
+ break;
+ }
+ }
+
+ m_msg << MSG::VERBOSE << ">>> CheckParentBranch is DONE : "
+ << ( accepted ? "accept" : "reject" )
+ << " vtx= " << vtx->barcode()
+ << endreq;
+ return accepted;
+}
+
+bool McVtxFilter::checkChildBranch( const HepMC::GenVertex * vtx ) const
+{
+ m_msg << MSG::VERBOSE << "In checkChildBranch..." << endreq;
+
+ if ( m_msg.level() <= MSG::VERBOSE ) {
+ vtx->print();
+ }
+
+ /// Check we aren't in the "any particle" case
+ if ( m_childList.empty() ) {
+ return true;
+ }
+
+ /// Check there is enough outgoing particles in the current vertex
+ if ( static_cast<unsigned int>(vtx->particles_out_size()) < m_childList.size() ) return false;
+
+ m_msg << MSG::VERBOSE << "Number of list of children : "
+ << m_childList.size() << endreq;
+
+ std::vector<int> childIds;
+ for ( HepMC::GenVertex::particles_out_const_iterator itrPart = vtx->particles_out_const_begin();
+ itrPart != vtx->particles_out_const_end();
+ ++itrPart ) {
+ childIds.push_back( (*itrPart)->pdg_id() );
+ }
+
+ AnalysisUtils::Permutation<std::vector<int> > permute( &childIds,
+ m_childList.size() );
+ std::vector<int> children;
+
+ bool accepted = false;
+ while ( permute.get(children) && !accepted ) {
+ accepted = true;
+ const unsigned int iMax = std::min( m_childList.size(), childIds.size() );
+ for ( unsigned int i = 0; i != iMax; ++i ) {
+ const bool hasInList =
+ m_childList[i]->hasInList( static_cast<PDG::pidType>(children[i]),
+ m_matchSign );
+ if ( !hasInList ) {
+ // this permutation is not suiting, going to the next one
+ // (if any)
+ accepted = false;
+ break;
+ }
+ }
+ }
+
+ m_msg << MSG::VERBOSE << ">>> CheckChildBranch is DONE : "
+ << ( accepted ? "accept" : "reject" )
+ << " vtx= " << vtx->barcode()
+ << endreq;
+ return accepted;
+}
+
+bool McVtxFilter::checkTwoBodyDecay( const HepMC::GenVertex * vtx ) const
+{
+ m_msg << MSG::VERBOSE << "In checkTwoBodyDecay..." << endreq;
+
+ /// First check the parent branch matching decision
+ /// if it doesn't fulfil our requirements, it is not worth
+ /// continuing the process
+ const bool parentTree = checkParentBranch( vtx );
+ if ( parentTree == false ) {
+ return false;
+ }
+
+ /// Now, handle the child branch
+ m_msg << MSG::VERBOSE << ">>> Check child branch" << endreq;
+
+ /// Cache the 2 particle candidate lists
+ const ParticleCandidateList * children1 = m_childList[0];
+ const ParticleCandidateList * children2 = m_childList[1];
+
+ /// Cache the id of the outgoing particles of the vertex being analysed
+ HepMC::GenVertex::particles_out_const_iterator
+ itrPart = vtx->particles_out_const_begin();
+ const PDG::pidType pdgId1 = static_cast<PDG::pidType>((*itrPart)->pdg_id());
+ ++itrPart;
+ const PDG::pidType pdgId2 = static_cast<PDG::pidType>((*itrPart)->pdg_id());
+
+ /// Loop over candidates for the 1st child
+ for( ParticleCandidateList::const_iterator itr1 = children1->begin();
+ itr1 != children1->end();
+ ++itr1 ) {
+ /// Loop over candidate for the 2nd child
+ for( ParticleCandidateList::const_iterator itr2 = children2->begin();
+ itr2 != children2->end();
+ ++itr2 ) {
+ m_msg << MSG::VERBOSE << "Checking the pair : "
+ << (*itr1) << "/" << (*itr2)
+ << endreq;
+
+ /// If the strict match sign has been required, we check if
+ /// the PDG ids are matching
+ if ( m_matchSign &&
+ ( ( (*itr1) == pdgId1 && (*itr2) == pdgId2 ) ||
+ ( (*itr1) == pdgId2 && (*itr2) == pdgId1 ) ) ) {
+ m_msg << MSG::VERBOSE << "Strict sign matching found !" << endreq;
+ return true;
+ /// if we are in a loose sign request, we check only that the
+ /// absolute values of the pdg IDs are matching
+ } else if ( ( std::abs(*itr1) == std::abs(pdgId1) &&
+ std::abs(*itr2) == std::abs(pdgId2) ) ||
+ ( std::abs(*itr1) == std::abs(pdgId2) &&
+ std::abs(*itr2) == std::abs(pdgId1) ) ) {
+ m_msg << MSG::VERBOSE << "Loose sign matching found !" << endreq;
+ return true;
+ }
+ m_msg << MSG::VERBOSE << "Checking next pair" << endreq;
+ }//> loop over 2nd child's candidates
+ }//> loop over 1st child's candidates
+
+ /// If we are here, then the vertex candidate didn't fulfil the
+ /// requirements we have setup
+ m_msg << MSG::VERBOSE << ">>> CheckTwoBodyDecay is DONE." << endreq;
+ return false;
+}
+
+///////////////////////////////////////////////////////////////////
+// Operators
+///////////////////////////////////////////////////////////////////
+
+MsgStream& operator<<( MsgStream & msg, const McVtxFilter &obj )
+{
+ std::stringstream out;
+ obj.dump( out );
+ msg << out.str() << endreq;
+ return msg;
+}
+
+std::ostream& operator<<( std::ostream& out, const McVtxFilter &obj )
+{
+ obj.dump( out );
+ return out;
+}
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McParticleUtils.sh b/PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McParticleUtils.sh
new file mode 100755
index 000000000000..408b03eebb4f
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McParticleUtils.sh
@@ -0,0 +1,13 @@
+#!/bin/sh
+# This just selects the directory to run cmt commands
+# But with project builds and as this package (McParticleUtils) can be moved
+# from some project to another, we wildcard the AtlasXyzRelease directory
+cd ${NIGHTLYAREA}/Atlas*Release/cmt
+#
+cmt broadcast -select=McParticleUtils make CppUnit
+stat=$?
+if [ "$stat" != "0" ]; then
+ echo " -------------------------------- "
+ echo " FAILURE : test McParticleUtils "
+ echo " -------------------------------- "
+fi
diff --git a/PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McVtxFilterTest_CppUnit.cxx b/PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McVtxFilterTest_CppUnit.cxx
new file mode 100755
index 000000000000..3a7853fb446d
--- /dev/null
+++ b/PhysicsAnalysis/TruthParticleID/McParticleUtils/test/McVtxFilterTest_CppUnit.cxx
@@ -0,0 +1,752 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/// Class to test the McVtxFilter class
+/// Author : S.Binet<binet@cern.ch>
+
+// CppUnit includes
+#include<cppunit/extensions/HelperMacros.h>
+#include<cppunit/Exception.h>
+
+// STL includes
+#include <iostream>
+#include <stdexcept>
+#include <cmath>
+
+// Framework includes
+#include "StoreGate/StoreGateSvc.h"
+
+// CLHEP includes
+#include "CLHEP/Vector/LorentzVector.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+// HepMC includes
+#include "HepMC/GenEvent.h"
+#include "HepMC/GenVertex.h"
+#include "HepMC/GenParticle.h"
+
+// DataModel includes
+#include "DataModel/DataVector.h"
+
+// McParticleUtils includes
+#include "McParticleUtils/McVtxFilter.h"
+
+class McVtxFilterTest : public CppUnit::TestFixture
+{
+ /// Definition of the unit test suite "FilterTest"
+ CPPUNIT_TEST_SUITE( McVtxFilterTest );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testCopyMcVtxFilter );
+ CPPUNIT_TEST( testSettersAndGetters );
+
+ CPPUNIT_TEST( testFullVtx );
+
+ CPPUNIT_TEST( testZeeVertices );
+ CPPUNIT_TEST( testTopWbggVertices );
+ CPPUNIT_TEST( testWlnuVertices );
+ CPPUNIT_TEST( testQuarkVertices );
+
+ CPPUNIT_TEST( testStableParticle );
+
+ /// end the definition test suite
+ CPPUNIT_TEST_SUITE_END();
+
+private:
+ HepMC::GenEvent * m_evt;
+ int m_bcZee;
+ int m_bcZgee;
+ int m_bcTopWbgg;
+ typedef HepLorentzVector HLV_t;
+
+public:
+
+ /// Set up the data members
+ void setUp()
+ {
+ const int signalProcessId = 1000082;
+ const int evtNbr = 1;
+ m_evt = new HepMC::GenEvent( signalProcessId, evtNbr );
+ m_evt->set_event_scale( -1 );
+ m_evt->set_alphaQCD( -1 );
+ m_evt->set_alphaQED( -1 );
+ std::vector<double> weights(3);
+ weights[0] = 1;
+ weights[1] = 1;
+ weights[2] = 1;
+ std::vector<long> rdmStates(2);
+ rdmStates[0] = 85909879;
+ rdmStates[1] = 9707499;
+ m_evt->weights() = weights;
+ m_evt->set_random_states( rdmStates );
+
+ // Add 2 vertices
+ HepMC::GenVertex * v1 = new HepMC::GenVertex();
+ m_evt->add_vertex( v1 );
+ v1->add_particle_in( new HepMC::GenParticle( HLV_t(0,0,
+ 7000*GeV,
+ 7000*GeV),
+ 2212, 3 ) );
+ HepMC::GenVertex* v2 = new HepMC::GenVertex;
+ m_evt->add_vertex( v2 );
+ v2->add_particle_in( new HepMC::GenParticle( HLV_t(0,0,
+ -7000*GeV,
+ +7000*GeV),
+ 2212, 3 ) );
+
+ //
+ // create the outgoing particles of v1 and v2
+ HepMC::GenParticle* p3 =
+ new HepMC::GenParticle( HLV_t(.750*GeV,
+ -1.569*GeV,
+ 32.191*GeV,
+ 32.238*GeV), 1, 3 );
+ v1->add_particle_out( p3 );
+ HepMC::GenParticle* p4 =
+ new HepMC::GenParticle( HLV_t( -3.047*GeV,
+ -19.*GeV,
+ -54.629*GeV,
+ 57.920*GeV), -2, 3 );
+ v2->add_particle_out( p4 );
+
+ //
+ // create v3
+ HepMC::GenVertex* v3 = new HepMC::GenVertex();
+ m_evt->add_vertex( v3 );
+ v3->add_particle_in( p3 );
+ v3->add_particle_in( p4 );
+ v3->add_particle_out(
+ new HepMC::GenParticle( HLV_t(-3.813,0.113,-1.833,4.233 ),
+ 22, 1 )
+ );
+ HepMC::GenParticle* p5 =
+ new HepMC::GenParticle( HLV_t(1.517,-20.68,-20.605,85.925),
+ -24,3);
+ v3->add_particle_out( p5 );
+
+ //
+ // tell the event which vertex is the signal process vertex
+ m_evt->set_signal_process_vertex( v3 );
+
+ //
+ // add a Z0->4gammas,e+,e-
+ //
+ HepMC::GenVertex * vZgee = new HepMC::GenVertex;
+ m_evt->add_vertex( vZgee );
+ // Z0
+ vZgee->add_particle_in( new HepMC::GenParticle( HLV_t( -4.49e+04,
+ +8.36e+03,
+ -2.70e+05,
+ +2.89e+05 ),
+ 23, 2 ) );
+ // Gammas
+ vZgee->add_particle_out( new HepMC::GenParticle( HLV_t( -1.28e+03,
+ +1.03e+03,
+ -5.47e+03,
+ +5.71e+03 ),
+ 22, 1 ) );
+ vZgee->add_particle_out( new HepMC::GenParticle( HLV_t( +3.89e+02,
+ -3.16e+02,
+ -6.69e+03,
+ +6.70e+03 ),
+ 22, 1 ) );
+ vZgee->add_particle_out( new HepMC::GenParticle( HLV_t( +7.34e+00,
+ -2.71e+01,
+ -4.12e+01,
+ +4.98e+01 ),
+ 22, 1 ) );
+ vZgee->add_particle_out( new HepMC::GenParticle( HLV_t( -1.36e+02,
+ +9.38e+01,
+ -5.62e+02,
+ +5.86e+02 ),
+ 22, 1 ) );
+ // Electrons
+ vZgee->add_particle_out( new HepMC::GenParticle( HLV_t( +8.01e+03,
+ -2.96e+04,
+ -4.50e+04,
+ +5.44e+04 ),
+ -11, 1 ) );
+ vZgee->add_particle_out( new HepMC::GenParticle( HLV_t( -5.19e+04,
+ +3.72e+04,
+ -2.13e+05,
+ +2.22e+05 ),
+ 11, 1 ) );
+ // store its barcode for later use
+ m_bcZgee = vZgee->barcode();
+
+ //
+ // Add a Z->e+e-
+ HepMC::GenVertex * vZee = new HepMC::GenVertex;
+ m_evt->add_vertex( vZee );
+ vZee->add_particle_in( new HepMC::GenParticle( HLV_t( +7.29e+03,
+ +2.34e+04,
+ +2.81e+05,
+ +2.96e+05 ),
+ 23, 2 ) );
+
+ vZee->add_particle_out( new HepMC::GenParticle( HLV_t( +2.74e+04,
+ -1.83e+04,
+ +4.70e+04,
+ +5.74e+04 ),
+ 11, 1 ) );
+ vZee->add_particle_out( new HepMC::GenParticle( HLV_t( -2.01e+04,
+ +4.17e+04,
+ +2.34e+05,
+ +2.38e+05 ),
+ -11, 1 ) );
+ m_bcZee = vZee->barcode();
+
+ //
+ // Add a t->W+bgg
+ HepMC::GenVertex * vtWbgg = new HepMC::GenVertex;
+ m_evt->add_vertex( vtWbgg );
+ // top
+ vtWbgg->add_particle_in(new HepMC::GenParticle(HLV_t(-2.35e+05,
+ +7.34e+04,
+ +3.60e+04,
+ +3.04e+05),
+ 6, 3 ) );
+
+ // Wbgg
+ vtWbgg->add_particle_out(new HepMC::GenParticle(HLV_t(-1.09e+05,
+ +6.99e+04,
+ -3.86e+04,
+ +1.57e+05),
+ 24, 2 ) );
+ vtWbgg->add_particle_out(new HepMC::GenParticle(HLV_t(-9.23e+04,
+ +2.54e+03,
+ +5.32e+04,
+ +1.07e+05),
+ 5, 2 ) );
+ vtWbgg->add_particle_out(new HepMC::GenParticle(HLV_t(-4.76e+03,
+ +6.72e+02,
+ +2.90e+03,
+ +5.62e+03),
+ 21, 2 ) );
+ vtWbgg->add_particle_out(new HepMC::GenParticle(HLV_t(-2.93e+04,
+ +2.13e+02,
+ +1.85e+04,
+ +3.46e+04),
+ 21, 2 ) );
+ m_bcTopWbgg = vtWbgg->barcode();
+
+ }
+
+ /// destroy any on-the-heap-created data member
+ void tearDown()
+ {
+ if ( m_evt ) {
+ delete m_evt;
+ m_evt = 0;
+ }
+ }
+
+ /// Test the McVtxFilter constructors
+ void testConstructor()
+ {
+ const unsigned int defaultCand = 0;
+ McVtxFilter filter;
+ CPPUNIT_ASSERT( false == filter.matchSign() );
+ CPPUNIT_ASSERT( false == filter.matchBranches() );
+ CPPUNIT_ASSERT( "" == filter.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == filter.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == filter.childList().size() );
+
+ McVtxFilter f2( false, false );
+ CPPUNIT_ASSERT( false == f2.matchSign() );
+ CPPUNIT_ASSERT( false == f2.matchBranches() );
+ CPPUNIT_ASSERT( "" == f2.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == f2.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == f2.childList().size() );
+
+ McVtxFilter f3( true, false );
+ CPPUNIT_ASSERT( true == f3.matchSign() );
+ CPPUNIT_ASSERT( false == f3.matchBranches() );
+ CPPUNIT_ASSERT( "" == f3.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == f3.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == f3.childList().size() );
+
+ McVtxFilter f4( false, true );
+ CPPUNIT_ASSERT( false == f4.matchSign() );
+ CPPUNIT_ASSERT( true == f4.matchBranches() );
+ CPPUNIT_ASSERT( "" == f4.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == f4.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == f4.childList().size() );
+
+ McVtxFilter f5( f4 );
+ CPPUNIT_ASSERT( false == f5.matchSign() );
+ CPPUNIT_ASSERT( true == f5.matchBranches() );
+ CPPUNIT_ASSERT( "" == f5.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == f5.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == f5.childList().size() );
+
+ McVtxFilter f6 = f4;
+ CPPUNIT_ASSERT( false == f6.matchSign() );
+ CPPUNIT_ASSERT( true == f6.matchBranches() );
+ CPPUNIT_ASSERT( "" == f6.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == f6.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == f6.childList().size() );
+
+ {
+ McVtxFilter f( "6 -> 24 + -5" );
+ CPPUNIT_ASSERT( false == f.matchSign() );
+ CPPUNIT_ASSERT( false == f.matchBranches() );
+ CPPUNIT_ASSERT( "6 -> 24 + -5" == f.decayPattern() );
+ CPPUNIT_ASSERT( 1 == f.parentList().size() );
+ CPPUNIT_ASSERT( 2 == f.childList().size() );
+ CPPUNIT_ASSERT( f.isFullVtx() );
+ }
+ {
+ McVtxFilter f( "6 -> 24 + -5", true, false );
+ CPPUNIT_ASSERT( true == f.matchSign() );
+ CPPUNIT_ASSERT( false == f.matchBranches() );
+ CPPUNIT_ASSERT( "6 -> 24 + -5" == f.decayPattern() );
+ CPPUNIT_ASSERT( 1 == f.parentList().size() );
+ CPPUNIT_ASSERT( 2 == f.childList().size() );
+ CPPUNIT_ASSERT( f.isFullVtx() );
+ }
+ {
+ McVtxFilter f( "6 -> 24 + -5", false, true );
+ CPPUNIT_ASSERT( false == f.matchSign() );
+ CPPUNIT_ASSERT( true == f.matchBranches() );
+ CPPUNIT_ASSERT( "6 -> 24 + -5" == f.decayPattern() );
+ CPPUNIT_ASSERT( 1 == f.parentList().size() );
+ CPPUNIT_ASSERT( 2 == f.childList().size() );
+ CPPUNIT_ASSERT( f.isFullVtx() );
+ }
+ {
+ McVtxFilter f( "6 -> ", false, true );
+ CPPUNIT_ASSERT( false == f.matchSign() );
+ CPPUNIT_ASSERT( true == f.matchBranches() );
+ CPPUNIT_ASSERT( "6 -> " == f.decayPattern() );
+ CPPUNIT_ASSERT( 1 == f.parentList().size() );
+ CPPUNIT_ASSERT( 0 == f.childList().size() );
+ CPPUNIT_ASSERT( !f.isFullVtx() );
+ }
+ }
+
+ /// Test that McVtxFilter cuts are well copied
+ void testCopyMcVtxFilter()
+ {
+ const unsigned int defaultCand = 0;
+ McVtxFilter filter( false, true );
+ CPPUNIT_ASSERT( false == filter.matchSign() );
+ CPPUNIT_ASSERT( true == filter.matchBranches() );
+ CPPUNIT_ASSERT( "" == filter.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == filter.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == filter.childList().size() );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+
+ filter.setDecayPattern( "6 -> 24 + -5" );
+ CPPUNIT_ASSERT( "6 -> 24 + -5" == filter.decayPattern() );
+ CPPUNIT_ASSERT( 1 == filter.parentList().size() );
+ CPPUNIT_ASSERT( 2 == filter.childList().size() );
+ CPPUNIT_ASSERT( true == filter.isFullVtx() );
+
+ McVtxFilter copy;
+ CPPUNIT_ASSERT( false == copy.matchSign() );
+ CPPUNIT_ASSERT( false == copy.matchBranches() );
+ CPPUNIT_ASSERT( "" == copy.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == copy.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == copy.childList().size() );
+ CPPUNIT_ASSERT( false == copy.isFullVtx() );
+
+ copy.setFilter( &filter );
+ CPPUNIT_ASSERT( false == copy.matchSign() );
+ CPPUNIT_ASSERT( true == copy.matchBranches() );
+ CPPUNIT_ASSERT( "6 -> 24 + -5" == copy.decayPattern() );
+ CPPUNIT_ASSERT( 1 == copy.parentList().size() );
+ CPPUNIT_ASSERT( 2 == copy.childList().size() );
+ CPPUNIT_ASSERT( true == copy.isFullVtx() );
+ }
+
+ /// Test setters and getters
+ void testSettersAndGetters()
+ {
+ std::cout << std::endl;
+ //
+ // Print event
+ //
+ m_evt->print();
+
+ const unsigned int defaultCand = 0;
+ McVtxFilter filter( false, true );
+ CPPUNIT_ASSERT( false == filter.matchSign() );
+ CPPUNIT_ASSERT( true == filter.matchBranches() );
+ CPPUNIT_ASSERT( "" == filter.decayPattern() );
+ CPPUNIT_ASSERT( defaultCand == filter.parentList().size() );
+ CPPUNIT_ASSERT( defaultCand == filter.childList().size() );
+
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( true == filter.matchSign() );
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( false == filter.matchSign() );
+
+ filter.setMatchBranches( false );
+ CPPUNIT_ASSERT( false == filter.matchBranches() );
+ filter.setMatchBranches( true );
+ CPPUNIT_ASSERT( true == filter.matchBranches() );
+
+ filter.setDecayPattern( "6 -> 24 + -5" );
+ CPPUNIT_ASSERT( "6 -> 24 + -5" == filter.decayPattern() );
+ CPPUNIT_ASSERT( 1 == filter.parentList().size() );
+ CPPUNIT_ASSERT( 2 == filter.childList().size() );
+
+ filter.setDecayPattern( "23 -> -11|-13 + 11|13|15" );
+ CPPUNIT_ASSERT( "23 -> -11|-13 + 11|13|15" == filter.decayPattern() );
+ CPPUNIT_ASSERT( 1 == filter.parentList().size() );
+ CPPUNIT_ASSERT( 2 == filter.childList().size() );
+
+ bool matchSign = false;
+ CPPUNIT_ASSERT( (filter.parentList()[0])->hasInList(PDG::Z0, matchSign) );
+ matchSign = true;
+ CPPUNIT_ASSERT( (filter.parentList()[0])->hasInList(PDG::Z0, matchSign) );
+
+ matchSign = false;
+ CPPUNIT_ASSERT((filter.childList()[0])->hasInList(PDG::e_plus, matchSign));
+ CPPUNIT_ASSERT((filter.childList()[0])->hasInList(PDG::mu_plus,matchSign));
+ matchSign = true;
+ CPPUNIT_ASSERT((filter.childList()[0])->hasInList(PDG::e_plus, matchSign));
+ CPPUNIT_ASSERT((filter.childList()[0])->hasInList(PDG::mu_plus,matchSign));
+
+ matchSign = true;
+ CPPUNIT_ASSERT(!(filter.childList()[0])->hasInList( PDG::e_minus,
+ matchSign));
+ CPPUNIT_ASSERT(!(filter.childList()[0])->hasInList( PDG::mu_minus,
+ matchSign));
+
+ CPPUNIT_ASSERT( 1 == (filter.parentList()[0])->size() );
+ std::vector<PDG::pidType> parent;
+ for ( std::list<PDG::pidType>::const_iterator itr =
+ (filter.parentList()[0])->begin();
+ itr != (filter.parentList()[0])->end();
+ ++itr ) {
+ parent.push_back( *itr );
+ }
+
+ CPPUNIT_ASSERT( PDG::Z0 == parent[0] );
+
+ CPPUNIT_ASSERT( 2 == (filter.childList()[0])->size() );
+ std::vector<PDG::pidType> child1;
+ for ( std::list<PDG::pidType>::const_iterator itr =
+ (filter.childList()[0])->begin();
+ itr != (filter.childList()[0])->end();
+ ++itr ) {
+ child1.push_back( *itr );
+ }
+ CPPUNIT_ASSERT( 2 == child1.size() );
+ CPPUNIT_ASSERT( PDG::e_plus == child1[0] );
+ CPPUNIT_ASSERT( PDG::mu_plus == child1[1] );
+
+ CPPUNIT_ASSERT( 3 == (filter.childList()[1])->size() );
+ std::vector<PDG::pidType> child2;
+ for ( std::list<PDG::pidType>::const_iterator itr =
+ (filter.childList()[1])->begin();
+ itr != (filter.childList()[1])->end();
+ ++itr ) {
+ child2.push_back( *itr );
+ }
+ CPPUNIT_ASSERT( 3 == child2.size() );
+
+ CPPUNIT_ASSERT( PDG::e_minus == child2[0] );
+ CPPUNIT_ASSERT( PDG::mu_minus == child2[1] );
+ CPPUNIT_ASSERT( PDG::tau_minus == child2[2] );
+
+ }
+
+ /// Test full vtx
+ void testFullVtx()
+ {
+ HepMC::GenVertex * vtx = new HepMC::GenVertex;
+ m_evt->add_vertex( vtx );
+ vtx->add_particle_in( new HepMC::GenParticle( HLV_t( -2.45e+04,
+ +1.88e+04,
+ -8.65e+05,
+ +8.65e+05 ),
+ 22, 3 ) );
+ vtx->add_particle_out( new HepMC::GenParticle( HLV_t( -2.45e+04,
+ +1.88e+04,
+ -8.65e+05,
+ +8.65e+05 ),
+ 22, 1 ) );
+
+ McVtxFilter filter;
+ filter.setDecayPattern( "22->" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "->22" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "24->-24" );
+ CPPUNIT_ASSERT( true == filter.isFullVtx() );
+ CPPUNIT_ASSERT( !filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "-> -24 + 24" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( !filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "22 -> -11 + 11" );
+ CPPUNIT_ASSERT( true == filter.isFullVtx() );
+ CPPUNIT_ASSERT( !filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "22 -> -11" );
+ CPPUNIT_ASSERT( true == filter.isFullVtx() );
+ CPPUNIT_ASSERT( !filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "22 -> 22" );
+ CPPUNIT_ASSERT( true == filter.isFullVtx() );
+ CPPUNIT_ASSERT( filter.isAccepted( vtx ) );
+
+ }
+
+ /// Test some vertices
+ void testZeeVertices()
+ {
+ McVtxFilter filter;
+ filter.setDecayPattern( "23 -> -11 + 11" );
+
+ const HepMC::GenVertex * vtx = m_evt->barcode_to_vertex(m_bcZee);
+ CPPUNIT_ASSERT( 0 != vtx );
+
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setMatchBranches( true );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ // Retrieve a vtx with FSR showering
+ vtx = m_evt->barcode_to_vertex(m_bcZgee);
+ CPPUNIT_ASSERT( 0 != vtx );
+
+ filter.setMatchSign( false );
+ filter.setMatchBranches( true );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtx) );
+
+ filter.setMatchSign( false );
+ filter.setMatchBranches( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setMatchSign( true );
+ filter.setMatchBranches( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setMatchSign( true );
+ filter.setMatchBranches( true );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtx) );
+ }
+
+ /// Test some vertices
+ void testTopWbggVertices()
+ {
+ const HepMC::GenVertex * vtx = m_evt->barcode_to_vertex(m_bcTopWbgg);
+ CPPUNIT_ASSERT( 0 != vtx );
+
+ {
+ McVtxFilter filter;
+ filter.setDecayPattern( "6 -> 24 + 5" );
+
+ filter.setMatchSign( true );
+ filter.setMatchBranches( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setMatchSign( true );
+ filter.setMatchBranches( true );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtx) );
+
+ filter.setMatchSign( false );
+ filter.setMatchBranches( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setMatchSign( false );
+ filter.setMatchBranches( true );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtx) );
+ }
+ }
+
+ /// Test some vertices
+ void testWlnuVertices()
+ {
+
+ HepMC::GenVertex * vtx = new HepMC::GenVertex;
+ m_evt->add_vertex( vtx );
+ vtx->add_particle_in( new HepMC::GenParticle( HLV_t( -6.76e+04,
+ +4.85e+03,
+ -1.46e+03,
+ +9.51e+04 ),
+ -24, 2 ) );
+ vtx->add_particle_out( new HepMC::GenParticle( HLV_t( -7.14e+04,
+ -6.17e+03,
+ +1.67e+04,
+ +7.36e+04 ),
+ 13, 1 ) );
+
+ vtx->add_particle_out( new HepMC::GenParticle( HLV_t( +3.75e+03,
+ +1.10e+04,
+ -1.81e+04,
+ +2.15e+04 ),
+ -14, 1 ) );
+ //const int bcWlnu = vtx->barcode();
+
+ McVtxFilter filter;
+ filter.setDecayPattern( "24 -> -13 + -14" );
+
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtx) );
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ filter.setDecayPattern( "-24 -> 13 + -14" );
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtx) );
+
+ }
+
+ /*
+GenVertex: -5 ID: 0 (X,cT):0
+ I: 2 5 21 +1.23e+04,-4.79e+03,+6.52e+04,+6.65e+04 3 -5
+ 6 21 +1.14e+02,+1.35e+04,-6.42e+04,+6.56e+04 3 -5
+ O: 2 7 5 +4.56e+04,+1.53e+04,-5.08e+04,+7.01e+04 3 -8
+ 8 -5 -3.32e+04,-6.56e+03,+5.18e+04,+6.20e+04 3 -9
+GenVertex: -39 ID: 0 (X,cT):0
+ I: 1 158 -5 -3.21e+04,-6.19e+03,+5.05e+04,+6.38e+04 2 -39
+ O: 2 159 -5 -1.29e+04,+1.12e+03,+3.50e+04,+3.77e+04 2 -92
+ 160 21 -1.92e+04,-7.30e+03,+1.55e+04,+2.61e+04 2 -40
+ */
+
+ /// Test some vertices
+ void testQuarkVertices()
+ {
+ // create a b->g+b vertex
+ HepMC::GenVertex * vtxgb = new HepMC::GenVertex;
+ m_evt->add_vertex( vtxgb );
+ vtxgb->add_particle_in( new HepMC::GenParticle( HLV_t( -3.21e+04,
+ -6.19e+03,
+ +5.05e+04,
+ +6.38e+04 ),
+ -5, 2 ) );
+ vtxgb->add_particle_out( new HepMC::GenParticle( HLV_t( -1.29e+04,
+ +1.12e+03,
+ +3.50e+04,
+ +3.77e+04 ),
+ -5, 2 ) );
+
+ vtxgb->add_particle_out( new HepMC::GenParticle( HLV_t( -1.92e+04,
+ -7.30e+03,
+ +1.55e+04,
+ +2.61e+04 ),
+ 21, 2 ) );
+
+ // create a gg->b+bbar vertex
+ HepMC::GenVertex * vtxbb = new HepMC::GenVertex;
+ m_evt->add_vertex( vtxbb );
+ vtxbb->add_particle_in( new HepMC::GenParticle( HLV_t( +1.23e+04,
+ -4.79e+03,
+ +6.52e+04,
+ +6.65e+04 ),
+ 21, 3 ) );
+ vtxbb->add_particle_in( new HepMC::GenParticle( HLV_t( +1.14e+02,
+ +1.35e+04,
+ -6.42e+04,
+ +6.56e+04 ),
+ 21, 3 ) );
+ vtxbb->add_particle_out( new HepMC::GenParticle( HLV_t( +4.56e+04,
+ +1.53e+04,
+ -5.08e+04,
+ +7.01e+04 ),
+ 5, 3 ) );
+
+ vtxbb->add_particle_out( new HepMC::GenParticle( HLV_t( -3.32e+04,
+ -6.56e+03,
+ +5.18e+04,
+ +6.20e+04 ),
+ -5, 3 ) );
+
+ McVtxFilter filter;
+ filter.setDecayPattern( " -> -5 + 5" );
+
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtxbb) );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtxgb) );
+
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( filter.isAccepted(vtxbb) );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtxgb) );
+
+ filter.setDecayPattern( "21+21 -> -5 + 5" );
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtxbb) );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtxgb) );
+
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( filter.isAccepted(vtxbb) );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtxgb) );
+
+ filter.setMatchBranches(true);
+ filter.setMatchSign( false );
+ CPPUNIT_ASSERT( filter.isAccepted(vtxbb) );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtxgb) );
+
+ filter.setMatchSign( true );
+ CPPUNIT_ASSERT( filter.isAccepted(vtxbb) );
+ CPPUNIT_ASSERT( !filter.isAccepted(vtxgb) );
+
+ }
+
+ /// Test Stable particles (no end_vertex)
+ void testStableParticle()
+ {
+ HepMC::GenVertex * vtx = new HepMC::GenVertex;
+ m_evt->add_vertex( vtx );
+ vtx->add_particle_in( new HepMC::GenParticle( HLV_t( -2.45e+04,
+ +1.88e+04,
+ -8.65e+05,
+ +8.65e+05 ),
+ 23, 3 ) );
+ HepMC::GenParticle * photon = 0;
+ photon = new HepMC::GenParticle( HLV_t( -2.45e+04,
+ +1.88e+04,
+ -8.65e+05,
+ +8.65e+05 ),
+ 22, 1 );
+ vtx->add_particle_out( photon );
+
+ // no end_vertex so : particle is stable
+ CPPUNIT_ASSERT( photon->end_vertex() == 0 );
+
+ McVtxFilter filter;
+ filter.setDecayPattern( "23->" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "22->" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( !filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "->23" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( !filter.isAccepted( vtx ) );
+
+ filter.setDecayPattern( "->22" );
+ CPPUNIT_ASSERT( false == filter.isFullVtx() );
+ CPPUNIT_ASSERT( filter.isAccepted( vtx ) );
+
+ }
+
+};
+
+/// Registration of the test suite "McVtxFilterTest" so it will be
+/// recognized by the main test program which drives the different tests
+CPPUNIT_TEST_SUITE_REGISTRATION( McVtxFilterTest );
+
+/// CppUnit test-driver common for all the cppunit test classes.
+/// In ATLAS sw release it is located in TestPolicy package
+#include <TestPolicy/CppUnit_SGtestdriver.cxx>
--
GitLab