Merge branch 'RHadron_Cleanup_Master' into 'master'

R-hadron cleanup for master

See merge request atlas/athena!22034

Generators/EvgenJobTransforms/python/EvgenConfig.py

-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+#  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 
-## Generators providing input events via the LHEF format (MC\@NLO produces LHEF
+## Generators providing input events via the LHEF format (MC@NLO produces LHEF
 ## when used in connection with Herwig++, which is fine since this variable is
 ## just used to determine the input file dummy-naming strategy for C++ generators)
 lhefGenerators = ["Lhef", # generic name: prefer to use the names below
                   "aMcAtNlo", "McAtNlo", "Powheg", "PowHel", "MadGraph", "CompHep", "CalcHep",
-                  "Whizard", "MCFM", "JHU", "MEtop", "Charybdis", "Charybdis2",
+                  "Whizard", "MCFM", "JHU", "MEtop", "Charybdis", "Charybdis2","BCVEGPY", "Dire4Pythia8",
                   "BlackMax", "QBH", "gg2ww", "gg2zz", "gg2vv", "HvyN", "VBFNLO", "FPMC", "ProtosLHEF"]
 
 ## A more general list of generators which provide partonic input, including non-LHEF ones

Generators/EvgenJobTransforms/python/evgenTrfArgs.py

-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+#  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 
 import PyJobTransforms.trfArgClasses as trfArgClasses
 
@@ -83,7 +83,7 @@ def addStdEvgenArgs(parser):
 
     parser.add_argument("--outputTXTFile", group="Evgen",
                         help="optional output TXT file for LHEF events, default is None",
-                        type=trfArgClasses.argFactory(trfArgClasses.argFile))
+                        type=trfArgClasses.argFactory(trfArgClasses.argLHEFile, io='output', runarg=True))
     parser.add_argument("--generatorRunMode", group="Evgen",
                         help="Run mode to pass on to the event generator",
                         default=trfArgClasses.argString("run", runarg=True),
@@ -99,4 +99,20 @@ def addStdEvgenArgs(parser):
                         default=trfArgClasses.argInt(0, runarg=True),
                         type=trfArgClasses.argFactory(trfArgClasses.argInt, runarg=True))
 
+    parser.add_argument('--lheOnly', '--lheOnly', group='Evgen',
+                        default=trfArgClasses.argInt(0, runarg=True),
+                        type=trfArgClasses.argFactory(trfArgClasses.argInt, runarg=True),
+                        help='Set to 1 if Evgen file should be removed')
+
+    parser.add_argument("--localPath", group="Evgen",
+                        help="Name of local path to JOs",
+                        type=trfArgClasses.argFactory(trfArgClasses.argString, runarg=True))
 
+    parser.add_argument('--saveList', '--saveList', group='Evgen',
+                        help="a comma-separated list of files to be saved in output or a single digit if the whole output directory is to be saved",
+                        type=trfArgClasses.argFactory(trfArgClasses.argSubstepList, runarg=True))
+
+    parser.add_argument('--cleanOut', '--cleanOut', group='Evgen',
+                        default=trfArgClasses.argInt(0, runarg=True),
+                        type=trfArgClasses.argFactory(trfArgClasses.argInt, runarg=True),
+                        help='Set to any digit !=0 if the output files should be cleaned up --saveList may help to save selected files')

Generators/EvgenJobTransforms/python/jo_proxy.py

+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+# Helper functions for setting up a job options proxy
+import os,shutil
+
+def get_immediate_subdirectories(a_dir):
+    return [name for name in os.listdir(a_dir)
+            if os.path.isdir(os.path.join(a_dir, name))]
+
+def mk_jo_proxy(targetbasepath, pkgname, proxypath, addtosearch=True):
+    "Make a JO proxy dir such that the MCxxJobOptions/dddd dirs contents are found via include(MCxxJobOptions/yyyy)"
+    if proxypath:
+        if os.path.exists(proxypath):
+            shutil.rmtree(proxypath)
+        os.mkdir(proxypath)
+    os.environ['LOCAL_INSTALL_DIR'] = (os.environ['JOBOPTSEARCHPATH']).split(":")[0]
+    os.environ['LOCAL_DATA_DIR'] = (os.environ['DATAPATH']).split(":")[0]
+
+    dirlist =  get_immediate_subdirectories(targetbasepath)
+    subdirlist=dirlist;
+    for dd in (dirlist):
+      if (('.svn' not in dd) and ('cmt' not in dd) and ('_joproxy' not in dd)):
+        deepdir = os.path.join(targetbasepath, dd)
+        subdirlist1 = get_immediate_subdirectories(deepdir)
+        subdirlist = subdirlist+["%s" % dd+"/%s" % item for item in subdirlist1]
+    for d in (subdirlist):
+        # TODO: we could _maybe_ add the appropriate share/DSIDxxxx/ dir to the path based on the jobConfig arg... too much magic?
+        if (('.svn' not in d) and ('cmt' not in d) and ('_joproxy' not in d) and ('share/' not in d)):
+          dpath = os.path.join(proxypath, d)
+          if proxypath:
+                   os.mkdir(dpath)
+
+          if ('nonStandard' in dpath):
+                dpath_ex = os.path.join(dpath,pkgname)
+                if proxypath:
+                   os.mkdir(dpath_ex)
+
+                os.symlink(os.path.join(targetbasepath, d), os.path.join(dpath_ex, "nonStandard"))
+          else :
+                os.symlink(os.path.join(targetbasepath, d), os.path.join(dpath, pkgname))
+          if addtosearch:
+            os.environ["JOBOPTSEARCHPATH"] = dpath+":"+os.environ["JOBOPTSEARCHPATH"]
+            os.environ["DATAPATH"] =os.path.join(targetbasepath, d)+":"+os.environ["DATAPATH"]
+
+    os.environ["JOBOPTSEARCHPATH"] = os.environ['LOCAL_INSTALL_DIR']+":"+os.environ["JOBOPTSEARCHPATH"]
+    os.environ["DATAPATH"] = os.environ['LOCAL_DATA_DIR']+":"+os.environ["DATAPATH"]
\ No newline at end of file

Generators/EvgenJobTransforms/scripts/AfterBurner_tf.py

-#! /usr/bin/env python
+#  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 
-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+#! /usr/bin/env python
 
 """
 Run event simulation and produce an EVNT file.
@@ -14,6 +14,7 @@ from PyJobTransforms.trfArgs import addAthenaArguments
 from PyJobTransforms.trfDecorators import stdTrfExceptionHandler, sigUsrStackTrace
 from EvgenJobTransforms.evgenTrfArgs import addStdEvgenArgs
 import PyJobTransforms.trfArgClasses as trfArgClasses
+from EvgenJobTransforms.jo_proxy import mk_jo_proxy
 
 
 ## Prodsys1 hack...
@@ -46,36 +47,6 @@ class EvgenExecutor(athenaExecutor):
                 zf.extractall()
                 zf.close()
 
-        def get_immediate_subdirectories(a_dir):
-            return [name for name in os.listdir(a_dir)
-                    if os.path.isdir(os.path.join(a_dir, name))]                
-
-        def mk_jo_proxy(targetbasepath, pkgname, proxypath, addtosearch=True):
-            "Make a JO proxy dir such that the MCxxJobOptions/dddd dirs contents are found via include(MCxxJobOptions/yyyy)"
-            if proxypath:
-                if os.path.exists(proxypath):
-                    shutil.rmtree(proxypath)
-                os.mkdir(proxypath)
-            os.environ['LOCAL_INSTALL_DIR'] = (os.environ['JOBOPTSEARCHPATH']).split(":")[0]
-            comdir = os.path.join(targetbasepath, "common")
-            subdirlist = get_immediate_subdirectories(comdir)
-            subdirlist1 = ['common/%s' % item for item in subdirlist]
-            dirlist = ['common','share','gencontrol','susycontrol']
-            for d in (dirlist+subdirlist1):
-                # TODO: we could _maybe_ add the appropriate share/DSIDxxxx/ dir to the path based on the jobConfig arg... too much magic?
-                if (d != 'common/.svn'):
-                  dpath = os.path.join(proxypath, d)
-
-                  if proxypath:
-                    os.mkdir(dpath)
-                    os.symlink(os.path.join(targetbasepath, d), os.path.join(dpath, pkgname))
-                  if addtosearch:
-                    os.environ["JOBOPTSEARCHPATH"] = dpath+":"+os.environ["JOBOPTSEARCHPATH"]
-                    os.environ["DATAPATH"] =os.path.join(targetbasepath, d)+":"+os.environ["DATAPATH"]
-
-            os.environ["JOBOPTSEARCHPATH"] = os.environ['LOCAL_INSTALL_DIR']+":"+os.environ["JOBOPTSEARCHPATH"]
-        
-
         ## Handle locating of evgen job options / fragments, either from a tarball or CVMFS
         if "evgenJobOpts" in self._trf.argdict: ## Use a specified JO tarball
             tarball = self._trf.argdict["evgenJobOpts"].value

Generators/EvgenJobTransforms/scripts/Generate_tf.py

 #! /usr/bin/env python
 
-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-
 """
+#  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 Run event simulation and produce an EVNT file.
 """
 
@@ -14,11 +13,11 @@ from PyJobTransforms.trfArgs import addAthenaArguments
 from PyJobTransforms.trfDecorators import stdTrfExceptionHandler, sigUsrStackTrace
 from EvgenJobTransforms.evgenTrfArgs import addStdEvgenArgs
 import PyJobTransforms.trfArgClasses as trfArgClasses
-
+from EvgenJobTransforms.jo_proxy import mk_jo_proxy
 
 ## Prodsys1 hack...
 # TODO: Remove!
-ListOfDefaultPositionalKeys=['--AMIConfig', '--AMITag', '--argJSON', '--asetup', '--athena', '--athenaMPMergeTargetSize', '--athenaopts', '--attempt', '--checkEventCount', '--command', '--dumpJSON', '--dumpPickle', '--ecmEnergy', '--env', '--eventAcceptanceEfficiency', '--evgenJobOpts', '--execOnly', '--fileValidation', '--firstEvent', '--ignoreErrors', '--ignoreFiles', '--ignorePatterns', '--imf', '--inputEVNT_PreFile', '--inputFileValidation', '--inputGenConfFile', '--inputGeneratorFile', '--jobConfig', '--jobid', '--maxEvents', '--orphanKiller', '--outputEVNTFile', '--outputEVNT_PreFile', '--outputFileValidation', '--outputNTUP_TRUTHFile', '--outputTXTFile', '--parallelFileValidation', '--postExec', '--postInclude', '--preExec', '--preInclude', '--printEvts', '--randomSeed', '--reportName', '--reportType', '--rivetAnas', '--runNumber', '--showGraph', '--showPath', '--showSteps', '--skipEvents', '--skipFileValidation', '--skipInputFileValidation', '--skipOutputFileValidation', '--steering', '--taskid', '--tcmalloc', '--valgrind', '--valgrindbasicopts', '--valgrindextraopts']
+ListOfDefaultPositionalKeys=['--AMIConfig', '--AMITag', '--argJSON', '--asetup', '--athena', '--athenaMPMergeTargetSize', '--athenaopts', '--attempt', '--checkEventCount', '--command', '--dumpJSON', '--dumpPickle', '--ecmEnergy', '--env', '--eventAcceptanceEfficiency', '--evgenJobOpts', '--execOnly', '--fileValidation', '--firstEvent', '--ignoreErrors', '--ignoreFiles', '--ignorePatterns', '--imf', '--inputEVNT_PreFile', '--inputFileValidation', '--inputGenConfFile', '--inputGeneratorFile', '--jobConfig', '--jobid', '--maxEvents', '--orphanKiller', '--outputEVNTFile', '--outputEVNT_PreFile', '--outputFileValidation', '--outputNTUP_TRUTHFile', '--outputTXTFile', '--parallelFileValidation', '--postExec', '--postInclude', '--preExec', '--preInclude', '--printEvts', '--randomSeed', '--reportName', '--reportType', '--rivetAnas', '--runNumber', '--showGraph', '--showPath', '--showSteps', '--skipEvents', '--skipFileValidation', '--skipInputFileValidation', '--skipOutputFileValidation', '--steering', '--taskid', '--tcmalloc', '--valgrind', '--valgrindbasicopts', '--valgrindextraopts', '--lheOnly', '--localPath', '--cleanOut', '--saveList']
 
 class EvgenExecutor(athenaExecutor):
     "Specialised trf executor class for event generation jobs"
@@ -45,49 +44,6 @@ class EvgenExecutor(athenaExecutor):
                 zf.extractall()
                 zf.close()
 
-        def get_immediate_subdirectories(a_dir):
-            return [name for name in os.listdir(a_dir)
-                    if os.path.isdir(os.path.join(a_dir, name))]                
-
-        def mk_jo_proxy(targetbasepath, pkgname, proxypath, addtosearch=True):
-            "Make a JO proxy dir such that the MCxxJobOptions/dddd dirs contents are found via include(MCxxJobOptions/yyyy)"
-            if proxypath:
-                if os.path.exists(proxypath):
-                    shutil.rmtree(proxypath)
-                os.mkdir(proxypath)
-            os.environ['LOCAL_INSTALL_DIR'] = (os.environ['JOBOPTSEARCHPATH']).split(":")[0]
-            os.environ['LOCAL_DATA_DIR'] = (os.environ['DATAPATH']).split(":")[0]
-
-            dirlist =  get_immediate_subdirectories(targetbasepath)
-            subdirlist=dirlist;
-            for dd in (dirlist):
-              if (('.svn' not in dd) and ('cmt' not in dd) and ('_joproxy' not in dd)):
-                deepdir = os.path.join(targetbasepath, dd)
-                subdirlist1 = get_immediate_subdirectories(deepdir)
-                subdirlist = subdirlist+["%s" % dd+"/%s" % item for item in subdirlist1]
-            for d in (subdirlist):
-                # TODO: we could _maybe_ add the appropriate share/DSIDxxxx/ dir to the path based on the jobConfig arg... too much magic?
-                if (('.svn' not in d) and ('cmt' not in d) and ('_joproxy' not in d) and ('share/' not in d)):
-                  dpath = os.path.join(proxypath, d)
-                  if proxypath:
-                           os.mkdir(dpath)
-
-                  if ('nonStandard' in dpath):
-                        dpath_ex = os.path.join(dpath,pkgname)
-                        if proxypath:
-                           os.mkdir(dpath_ex)
-                            
-                        os.symlink(os.path.join(targetbasepath, d), os.path.join(dpath_ex, "nonStandard"))
-                  else :      
-                        os.symlink(os.path.join(targetbasepath, d), os.path.join(dpath, pkgname))
-                  if addtosearch:
-                    os.environ["JOBOPTSEARCHPATH"] = dpath+":"+os.environ["JOBOPTSEARCHPATH"]
-                    os.environ["DATAPATH"] =os.path.join(targetbasepath, d)+":"+os.environ["DATAPATH"]
-
-            os.environ["JOBOPTSEARCHPATH"] = os.environ['LOCAL_INSTALL_DIR']+":"+os.environ["JOBOPTSEARCHPATH"]
-            os.environ["DATAPATH"] = os.environ['LOCAL_DATA_DIR']+":"+os.environ["DATAPATH"]
-        
-
         ## Handle locating of evgen job options / fragments, either from a tarball or CVMFS
         # read the JO used
         joparam = self._trf.argdict["jobConfig"].value
@@ -132,9 +88,16 @@ class EvgenExecutor(athenaExecutor):
             else:    
                 msg.info('ATLAS_SW_BASE path: %s' %sw_base)
 
+            local_path = None
+            if ("localPath" in self._trf.argdict ):
+                local_path = self._trf.argdict["localPath"].value
+                print("local path",local_path)
             cvmfs_path = os.path.join(sw_base, "atlas.cern.ch")
             
-            if os.path.exists(cvmfs_path):
+            if ((local_path is not None) and (os.path.exists(local_path))) :
+              mk_jo_proxy(local_path, "MC15JobOptions","_joproxy15")
+              print("JO fragments taken from local path i.e. ",local_path)             
+            elif os.path.exists(cvmfs_path):
                 # TODO: Make the package name configurable
                 if "MC14" in str(joparam):
                     cvmfs_mc14 = os.path.join(cvmfs_path, "repo/sw/Generators/MC14JobOptions/latest/")
@@ -154,10 +117,29 @@ class EvgenExecutor(athenaExecutor):
 
         ## Expand tarball input event and generator conf files, if provided
         if "inputGeneratorFile" in self._trf.argdict:
-            expand_if_archive(self._trf.argdict["inputGeneratorFile"].value)
+#            expand_if_archive(self._trf.argdict["inputGeneratorFile"].value)
+            myinputfiles = self._trf.argdict["inputGeneratorFile"].value
+            genInputFiles = myinputfiles.split(',')
+            for file in genInputFiles:
+               expand_if_archive(file)
         if "inputGenConfFile" in self._trf.argdict:
             expand_if_archive(self._trf.argdict["inputGenConfFile"].value)
 
+def move_files(main_dir,tmp_dir,whitelist):
+    files = os.listdir(tmp_dir)
+    print("list of files ",files)
+    print("white list ",whitelist)    
+    files.sort()
+    for f in files:
+       for i in whitelist:
+            if i in f:
+                src = tmp_dir+"/"+f
+                dest = main_dir+"/"+f
+#                print("src ",src)
+#                print("dest ",dest)
+                os.rename(src,dest)
+                break
+
 
 def getTransform():
     exeSet = set()
@@ -175,11 +157,67 @@ def getTransform():
 @sigUsrStackTrace
 def main():
     msg.info("This is %s" % sys.argv[0])
+    main_dir = os.getcwd()
     trf = getTransform()
     trf.parseCmdLineArgs(sys.argv[1:])
+#    print("main ",main_dir)
+    if (("cleanOut" in trf.argdict) and (trf.argdict["cleanOut"].value != 0)):
+      name_tmpdir = "tmprun"
+      tmp_dir = os.path.join(main_dir, name_tmpdir)
+#    print("tmpdir ",tmp_dir)
+      if os.path.isdir(tmp_dir):
+         shutil.rmtree(tmp_dir, ignore_errors=True)
+      os.mkdir("tmprun")
+      os.chdir("tmprun")
+      tmp_dir = os.getcwd()
+      whitelist_in = ['MC','group','TXT']
+      move_files(tmp_dir,main_dir,whitelist_in)
+      
+#    trf = getTransform()
+#    trf.parseCmdLineArgs(sys.argv[1:])
     trf.execute()
     trf.generateReport()
     msg.info("%s stopped at %s, trf exit code %d" % (sys.argv[0], time.asctime(), trf.exitCode))
+
+
+#    if "outputTXTFile" in trf.argdict:
+#      whitelist_out.append('TXT')
+
+# read files/dirs that should be saved and if present in cwd - remove
+
+    if (("cleanOut" in trf.argdict) and (trf.argdict["cleanOut"].value!=0)):
+#       print("mamy cleanOut !!!",trf.argdict["cleanOut"].value) 
+       whitelist_out = ['log.generate','.root']
+       if "outputTXTFile" in trf.argdict:
+          whitelist_out.append('TXT')
+       if "saveList" in trf.argdict:
+         saveList_dic= trf.argdict["saveList"].value
+         saveList_str= str(saveList_dic)
+         saveList_str=saveList_str[10:-3]
+         saveList= saveList_str.split(",")
+         for item in saveList:
+           test_ex = os.path.join(main_dir,str(item))
+           if os.path.isdir(test_ex):
+             shutil.rmtree(test_ex, ignore_errors=True)
+#           print("directory removed ", test_ex)
+           elif os.path.isfile(test_ex):
+             os.remove(test_ex)
+#           print("file removed ", test_ex)        
+#       new_list = list(list(int(a) for a in b) for b in T1 if a.isdigit())
+         if not saveList[0].isdigit():
+             whitelist_out=whitelist_out+saveList 
+
+       move_files(main_dir,tmp_dir,whitelist_out)
+       os.chdir(main_dir)
+       if "saveList" not in trf.argdict:
+         shutil.rmtree(tmp_dir, ignore_errors=True)
+#    elif firstSaveList!=1:
+       elif not saveList[0].isdigit():
+         shutil.rmtree(tmp_dir, ignore_errors=True)  
+
+    if (("lheOnly" in trf.argdict ) and (trf.argdict["lheOnly"].value == 1)):
+          outputName = ''.join(trf.argdict["outputEVNTFile"].value)
+          os.remove(outputName)
     sys.exit(trf.exitCode)
 
 

Generators/Pythia8_i/Pythia8_i/Pythia8_i.h

@@ -79,6 +79,8 @@ public:
   double pythiaVersion()const;
   
   static std::string    pythia_stream;
+
+  static std::string xmlpath();
     
 protected:
 
@@ -89,9 +91,7 @@ protected:
   HepMC::Pythia8ToHepMC m_pythiaToHepMC;
 
 private:
-  
-  static std::string xmlpath();
-  
+
   static std::string findValue(const std::string &command, const std::string &key);
   
   int m_internal_event_number;

Simulation/G4Extensions/G4ExternalDecay/CMakeLists.txt

@@ -7,23 +7,41 @@ atlas_subdir( G4ExternalDecay )
 
 # Declare the package's dependencies:
 atlas_depends_on_subdirs( PRIVATE
-                          Generators/Pythia_i )
+                          Generators/Pythia8_i )
 
 # External dependencies:
 find_package( CLHEP )
 find_package( Geant4 )
 find_package( HepMC )
 find_package( Lhapdf )
-find_package( Pythia6 )
+find_package( Pythia8 )
 find_package( XercesC )
+find_package( HEPUtils )
+find_package( MCUtils )
+
+set( extra_includes )
+set( extra_libs )
+if( HEPUTILS_FOUND )
+   list( APPEND extra_includes ${HEPUTILS_INCLUDE_DIRS} )
+   list( APPEND extra_libs ${HEPUTILS_LIBRARIES} )
+endif()
+if( MCUTILS_FOUND )
+   list( APPEND extra_includes ${MCUTILS_INCLUDE_DIRS} )
+   list( APPEND extra_libs ${MCUTILS_LIBRARIES} )
+endif()
+if( HEPMC_FOUND )
+   list( APPEND extra_includes ${HEPMC_INCLUDE_DIRS} )
+   list( APPEND extra_libs ${HEPMC_LIBRARIES} )
+endif()
+
 
 # Component(s) in the package:
 atlas_add_library( G4ExternalDecay
                    src/*.cxx
                    PUBLIC_HEADERS G4ExternalDecay
-                   INCLUDE_DIRS ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA6_INCLUDE_DIRS} ${LHAPDF_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
+                   INCLUDE_DIRS ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${LHAPDF_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS} ${extra_includes}
                    PRIVATE_INCLUDE_DIRS ${HEPMC_INCLUDE_DIRS}
                    DEFINITIONS ${CLHEP_DEFINITIONS}
-                   LINK_LIBRARIES ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA6_LIBRARIES} ${LHAPDF_LIBRARIES} ${GEANT4_LIBRARIES} Pythia_iLib
+                   LINK_LIBRARIES ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} ${LHAPDF_LIBRARIES} ${GEANT4_LIBRARIES} Pythia_iLib Pythia8_iLib ${extra_libs}
                    PRIVATE_LINK_LIBRARIES ${HEPMC_LIBRARIES} )
 

Simulation/G4Extensions/G4ExternalDecay/G4ExternalDecay/Pythia8ForDecays.h

+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Abused from Geant4 version of Pythai6.hh from extended examples
+
+#ifndef Pythia8ForDecays_H
+#define Pythia8ForDecays_H
+
+// For std::vector
+#include <vector>
+// G4 track function parameter
+#include "G4Track.hh"
+// For std::pair
+#include <utility>
+// For all the various Pythia8 classes used here
+#include "Pythia8_i/Pythia8_i.h"
+// For unique_ptr and once flag
+#include <memory>
+#include <mutex>
+
+class G4DynamicParticle;
+class G4ParticleDefinition;
+
+class Pythia8ForDecays
+{
+  public:
+
+   virtual ~Pythia8ForDecays() = default;
+
+   static Pythia8ForDecays *Instance();
+
+   /// Function that decays the RHadron; returns products in G4 format
+   void Py1ent(const G4Track&, std::vector<G4DynamicParticle*> &);
+
+  private:
+
+   /// Private constructor for singleton pattern
+   Pythia8ForDecays();
+
+   /// Helper for getting G4ParticleDefinition from PDG ID
+   G4ParticleDefinition* GetParticleDefinition(const int) const;
+
+   /// Fill a Pythia8 event with the information from a G4Track
+   void fillParticle(const G4Track&, Pythia8::Event& event) const;
+
+   /// Get the quarks from a gluino R-hadron.  From Pythia8 code.
+   std::pair<int,int> fromIdWithGluino( int idRHad, Pythia8::Rndm* rndmPtr) const;
+   std::pair<int,int> fromIdWithSquark( int idRHad) const;
+   bool isGluinoRHadron(int pdgId) const;
+
+   /// My own class; singleton pattern; thread safe for future-proofing
+   static std::unique_ptr<Pythia8ForDecays> s_instance;
+   static std::once_flag m_onceFlag;
+
+   /// The instance of Pythia8 that will do the work
+   std::unique_ptr<Pythia8::Pythia> m_pythia;
+};
+
+#endif
+

Simulation/G4Extensions/G4ExternalDecay/G4ExternalDecay/PythiaForDecays.h

-/*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// Abused from Geant4 version of Pythai6.hh from extended examples
-
-#ifndef PythiaForDecays_H
-#define PythiaForDecays_H
-
-#include <vector>
-#include "G4DynamicParticle.hh"
-
-class PythiaForDecays
-{
-  public:
-
-   PythiaForDecays();
-   virtual ~PythiaForDecays();
-
-   static PythiaForDecays *Instance();
-
-   int  Pycomp(int);
-   void Py1ent(int, int, double, double, double);
-
-   void ImportParticles( std::vector<G4DynamicParticle*> &);
-
-   void DecayRhadrons(int pdgid);
-
-  private:
-
-   G4ParticleDefinition* GetParticleDefinition(const int) const;
-
-   static  PythiaForDecays* s_instance;
-};
-
-#endif
-

Simulation/G4Extensions/G4ExternalDecay/src/Pythia8ForDecays.cxx

+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Abused from Pythia6.cc in Geant4 extended examples
+
+// My own class definition
+#include "G4ExternalDecay/Pythia8ForDecays.h"
+
+// Helper functions for Pythia8 and Pythia8 classes
+#include "Pythia8_i/Pythia8_i.h"
+
+// Pythia8 RHadrons code that we call into
+#include "Pythia8/RHadrons.h"
+
+// HepMC for translation into format Pythia likes
+#include "HepMC/IO_HEPEVT.h"
+#include "HepMC/GenEvent.h"
+#include "HepMC/GenParticle.h"
+
+// G4 classes for translation into G4 format
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+#include "G4ThreeVector.hh"
+#include "TLorentzVector.h"
+#include "MCUtils/HepMCUtils.h"
+
+// STL includes
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <fstream>
+
+std::unique_ptr<Pythia8ForDecays> Pythia8ForDecays::s_instance(nullptr);
+std::once_flag Pythia8ForDecays::m_onceFlag;
+
+Pythia8ForDecays* Pythia8ForDecays::Instance() {
+  std::call_once(m_onceFlag,
+        [] {
+            s_instance.reset(new Pythia8ForDecays);
+    });
+  return s_instance.get();
+}
+
+Pythia8ForDecays::Pythia8ForDecays()
+  {
+  // Pythia instance where RHadrons can decay
+  std::string docstring = Pythia8_i::xmlpath();
+  m_pythia = std::make_unique<Pythia8::Pythia>(docstring);
+  m_pythia->readString("SLHA:file = SLHA_INPUT.DAT");
+  m_pythia->readString("ProcessLevel:all = off");
+  m_pythia->readString("Init:showChangedSettings = off");
+  m_pythia->readString("RHadrons:allow = on");
+  m_pythia->readString("RHadrons:allowDecay = on");
+  m_pythia->readString("RHadrons:probGluinoball = 0.1");
+  m_pythia->readString("PartonLevel:FSR = off");
+  m_pythia->readString("Init:showAllParticleData = on");
+
+  // Process the file of commands left for us by the python layer
+  std::string line;
+  std::ifstream command_stream ("PYTHIA8_COMMANDS.TXT");
+  while(getline(command_stream,line)){
+    // Leaving it to the top-level to get this file right
+    m_pythia->readString(line);
+  }
+  command_stream.close();
+
+  m_pythia->init();
+
+}
+
+G4ParticleDefinition* Pythia8ForDecays::GetParticleDefinition(const int pdgEncoding) const
+{
+  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
+  G4ParticleDefinition* particleDefinition(nullptr);
+  if (pdgEncoding != 0) particleDefinition = particleTable->FindParticle(pdgEncoding);
+  return particleDefinition;
+}
+
+
+std::pair<int,int> Pythia8ForDecays::fromIdWithSquark( int idRHad) const
+{
+
+  // Find squark flavour content.
+  int idRSb            = m_pythia->settings.mode("RHadrons:idSbottom");
+  int idRSt            = m_pythia->settings.mode("RHadrons:idStop");
+  int idLight = (abs(idRHad) - 1000000) / 10;
+  int idSq    = (idLight < 100) ? idLight/10 : idLight/100;
+  int id1     = (idSq == 6) ? idRSt : idRSb;
+  if (idRHad < 0) id1 = -id1;
+
+  // Find light (di)quark flavour content.
+  int id2     =  (idLight < 100) ? idLight%10 : idLight%100;
+  if (id2 > 10) id2 = 100 * id2 + abs(idRHad)%10;
+  if ((id2 < 10 && idRHad > 0) || (id2 > 10 && idRHad < 0)) id2 = -id2;
+
+  // Done.
+  return std::make_pair( id1, id2);
+
+}
+
+// TODO: Would be nice for this to be a public function in Pythia8::RHadrons.hh
+std::pair<int,int> Pythia8ForDecays::fromIdWithGluino( int idRHad, Pythia8::Rndm* rndmPtr) const
+{
+  // Find light flavour content of R-hadron.
+  int idLight = (abs(idRHad) - 1000000) / 10;
+  int id1, id2, idTmp, idA, idB, idC;
+  double diquarkSpin1RH = 0.5;
+
+  // Gluinoballs: split g into d dbar or u ubar.
+  if (idLight < 100) {
+    id1 = (rndmPtr->flat() < 0.5) ? 1 : 2;
+    id2 = -id1;
+
+  // Gluino-meson: split into q + qbar.
+  } else if (idLight < 1000) {
+    id1 = (idLight / 10) % 10;
+    id2 = -(idLight % 10);
+    // Flip signs when first quark of down-type.
+    if (id1%2 == 1) {
+      idTmp = id1;
+      id1   = -id2;
+      id2   = -idTmp;
+    }
+
+  // Gluino-baryon: split to q + qq (diquark).
+  // Pick diquark at random, except if c or b involved.
+  } else {
+    idA = (idLight / 100) % 10;
+    idB = (idLight / 10) % 10;
+    idC = idLight % 10;
+    double rndmQ = 3. * rndmPtr->flat();
+    if (idA > 3) rndmQ = 0.5;
+    if (rndmQ < 1.) {
+      id1 = idA;
+      id2 = 1000 * idB + 100 * idC + 3;
+      if (idB != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
+    } else if (rndmQ < 2.) {
+      id1 = idB;
+      id2 = 1000 * idA + 100 * idC + 3;
+      if (idA != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
+    } else {
+      id1 = idC;
+      id2 = 1000 * idA + 100 * idB +3;
+      if (idA != idB && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
+    }
+  }
+  // Flip signs for anti-R-hadron.
+  if (idRHad < 0) {
+    idTmp = id1;
+    id1   = -id2;
+    id2   = -idTmp;
+  }
+
+  // Done.
+  return std::make_pair( id1, id2);
+
+}
+
+/// Add a G4Track to a Pythia8 event to make it a single-particle gun. The particle must be a colour singlet.
+/// Input: particle, Pythia8 event
+void Pythia8ForDecays::fillParticle(const G4Track& aTrack, Pythia8::Event& event) const
+{
+  // Reset event record to allow for new event.
+  event.reset();
+
+  // Select particle mass; where relevant according to Breit-Wigner.
+  double mm = aTrack.GetDynamicParticle()->GetMass();
+
+  // Store the particle in the event record.
+  event.append( aTrack.GetDefinition()->GetPDGEncoding(), 1, 0, 0, aTrack.GetMomentum().x()/CLHEP::GeV, aTrack.GetMomentum().y()/CLHEP::GeV,
+                aTrack.GetMomentum().z()/CLHEP::GeV, aTrack.GetDynamicParticle()->GetTotalEnergy()/CLHEP::GeV, mm/CLHEP::GeV);
+  // Note: this function returns an int, but we don't need or use its output
+}
+
+bool Pythia8ForDecays::isGluinoRHadron(int pdgId) const{
+  // Checking what kind of RHadron this is based on the digits in its PDGID
+  const unsigned short digitValue_q1 = MCUtils::PID::_digit(MCUtils::PID::Location::nq1,pdgId);
+  const unsigned short digitValue_l = MCUtils::PID::_digit(MCUtils::PID::Location::nl,pdgId);
+
+  // Gluino R-Hadrons have the form 109xxxx or 1009xxx
+  if (digitValue_l == 9 || (digitValue_l==0 && digitValue_q1 == 9) ){
+    // This is a gluino R-Hadron
+    return true;
+  }
+
+  // Special case : R-gluinoball
+  if (pdgId==1000993) return true;
+
+  // This is not a gluino R-Hadron (probably a squark R-Hadron)
+  return false;
+
+}
+
+void Pythia8ForDecays::Py1ent(const G4Track& aTrack, std::vector<G4DynamicParticle*> & particles)
+{
+
+  // Get members from Pythia8 instance where RHadrons can decay
+  Pythia8::Event& event      = m_pythia->event;
+  Pythia8::ParticleData& pdt = m_pythia->particleData;
+
+  // Use pythiaDecay information to fill event with the input particle
+  fillParticle(aTrack, event);
+
+  // Copy and paste of RHadron decay code
+  int    iRNow  = 1;
+  int    idRHad = event[iRNow].id();
+  double mRHad  = event[iRNow].m();
+  int    iR0    = 0;
+  int    iR2    = 0;
+
+  bool isTriplet = !isGluinoRHadron(idRHad);
+
+  // Find flavour content of squark or gluino R-hadron.
+  std::pair<int,int> idPair = (isTriplet) ? fromIdWithSquark( idRHad) : fromIdWithGluino( idRHad, &(m_pythia->rndm));
+  int id1 = idPair.first;
+  int id2 = idPair.second;
+
+  // Sharing of momentum: the squark/gluino should be restored
+  // to original mass, but error if negative-mass spectators.
+  int idRSb            = m_pythia->settings.mode("RHadrons:idSbottom");
+  int idRSt            = m_pythia->settings.mode("RHadrons:idStop");
+  int idRGo            = m_pythia->settings.mode("RHadrons:idGluino");
+  int idLight = (abs(idRHad) - 1000000) / 10;
+  int idSq    = (idLight < 100) ? idLight/10 : idLight/100;
+  int idRSq     = (idSq == 6) ? idRSt : idRSb;
+
+  // Handling R-Hadrons with anti-squarks
+  idRSq = idRSq * std::copysign(1, idRHad);
+
+  int idRBef = isTriplet ? idRSq : idRGo;
+
+  // Mass of the underlying sparticle
+  double mRBef = pdt.mSel(idRBef);
+
+  // Fraction of the RHadron mass given by the sparticle
+  double fracR = mRBef / mRHad;
+  int counter=0;
+  while (fracR>=1.){
+    if (counter==10){
+      G4cout << "Needed more than 10 attempts with constituent " << idRBef << " mass (" << mRBef << " above R-Hadron " << idRHad << " mass " << mRHad << G4endl;
+    } else if (counter>100){
+      G4cout << "Pythia8ForDecays::Py1ent ERROR   Failed >100 times. Constituent " << idRBef << " mass (" << mRBef << " above R-Hadron " << idRHad << " mass " << mRHad << G4endl;
+      return;
+    }
+    mRBef = pdt.mSel(idRBef);
+    fracR = mRBef / mRHad;
+    counter++;
+  }
+
+  // Squark case
+  if(isTriplet){
+    int col    = (event[idRBef].col() != 0) ? event.nextColTag() : 0;
+    int tmpSparticleColor = id1>0 ? col : 0;
+    int tmpSparticleAnticolor = id1>0 ? 0 : col;
+
+    // Store the constituents of a squark R-hadron.
+
+    // Sparticle
+    // (id, status, mother1, mother2, daughter1, daughter2, col, acol, px, py, pz, e, m=0., scaleIn=0., polIn=9.)
+    iR0 = event.append( id1, 106, iRNow, 0, 0, 0, tmpSparticleColor, tmpSparticleAnticolor, fracR * event[iRNow].p(), fracR * mRHad, 0.);
+    // Spectator quark
+    iR2 = event.append( id2, 106, iRNow, 0, 0, 0, tmpSparticleAnticolor, tmpSparticleColor, (1. - fracR) * event[iRNow].p(), (1. - fracR) * mRHad, 0.);
+  }
+  // Gluino case
+  else{
+    double mOffsetCloudRH = 0.2; // could be read from internal data?
+    double m1Eff  = pdt.constituentMass(id1) + mOffsetCloudRH;
+    double m2Eff  = pdt.constituentMass(id2) + mOffsetCloudRH;
+    double frac1 = (1. - fracR) * m1Eff / ( m1Eff + m2Eff);
+    double frac2 = (1. - fracR) * m2Eff / ( m1Eff + m2Eff);
+
+    // Two new colours needed in the breakups.
+    int col1 = event.nextColTag();
+    int col2 = event.nextColTag();
+
+    // Store the constituents of a gluino R-hadron.
+    iR0 = event.append( idRBef, 106, iRNow, 0, 0, 0, col2, col1, fracR * event[iRNow].p(), fracR * mRHad, 0.);
+    event.append( id1, 106, iRNow, 0, 0, 0, col1, 0, frac1 * event[iRNow].p(), frac1 * mRHad, 0.);
+    iR2 = event.append( id2, 106, iRNow, 0, 0, 0, 0, col2, frac2 * event[iRNow].p(), frac2 * mRHad, 0.);
+  }
+
+  // Mark R-hadron as decayed and update history.
+  event[iRNow].statusNeg();
+  event[iRNow].daughters( iR0, iR2);
+
+  // Generate events. Quit if failure.
+  if (!m_pythia->next()) {
+    m_pythia->forceRHadronDecays();
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Add the particles from the Pythia event into the GEANT particle vector
+  ///////////////////////////////////////////////////////////////////////////
+  particles.clear();
+
+  for(int i=0; i<m_pythia->event.size(); i++){
+    if ( m_pythia->event[i].status()<0 ) continue; // stable only
+    G4ThreeVector momentum( m_pythia->event[i].px() , m_pythia->event[i].py() , m_pythia->event[i].pz() );
+    momentum*=1000.0;//GeV to MeV
+    const G4ParticleDefinition * particleDefinition = GetParticleDefinition( m_pythia->event[i].id() );
+
+    if (!particleDefinition){
+      G4cout << "I don't know a definition for pdgid "<<m_pythia->event[i].id()<<"! Skipping it..." << G4endl;
+      continue;
+    }
+
+    G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, momentum);
+    particles.push_back( dynamicParticle );
+  }
+
+}
+

Simulation/G4Extensions/G4ExternalDecay/src/PythiaForDecays.cxx

-/*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// Abused from Pythia6.cc in Geant4 extended examples
-
-#include "G4ExternalDecay/PythiaForDecays.h"
-#include "Pythia_i/Pythia.h"
-#include "Pythia_i/Pydat3.h"
-#include "Pythia_i/Pydat2.h"
-#include "Pythia_i/Pydat1.h"
-#include "HepMC/IO_HEPEVT.h"
-#include "HepMC/GenEvent.h"
-#include "HepMC/GenParticle.h"
-
-#include "G4ParticleDefinition.hh"
-#include "G4ParticleTable.hh"
-#include "G4ThreeVector.hh"
-
-#include <iostream>
-#include <cstdlib>
-#include <cstring>
-#include <cmath>
-
-struct PPYJETS{
-  int   n;      
-  int npad;
-  int   k[5][4000];
-  double p[5][4000];
-  double v[5][4000];
-};
-
-extern "C" {
-  int  pycomp_(int*);
-  void py1ent_(int&, int&, double&, double&, double&);
-  void pyhepc_(int*);
-  PPYJETS pyjets_;
-  void pylist_(int*);
-  void pyexec_();
-}
-void pygldc_() {}//ACH for gluino decay
-void pysqdc_(int*) {}//ACH for stop/sbottom decay
-
-
-PythiaForDecays* PythiaForDecays::s_instance = 0;
-
-PythiaForDecays* PythiaForDecays::Instance() 
-{
-   if ( ! s_instance ) s_instance = new PythiaForDecays();
-   return s_instance;
-}
-
-PythiaForDecays::PythiaForDecays()  
-{
-   // Protect against multiple objects.   All access should be via the
-   // Instance member function. 
-   if ( s_instance ) {
-      std::cerr << "There's already an instance of PythiaForDecays" << std::endl;
-      exit (1);
-   }   
-}
-
-PythiaForDecays::~PythiaForDecays()
-{;} // It's ok little class, I never owned anything either...
-
-int PythiaForDecays::Pycomp(int kf) 
-{
-/// Interface with fortran routine pycomp
-   return pycomp_(&kf);
-}
-
-void PythiaForDecays::Py1ent(int ip, int kf, double pe, double theta, double phi)
-{
-/// Add one entry to the event record, i.e. either a parton or a
-/// particle. 
-///
-///  IP:   normally line number for the parton/particle. There are two
-///        exceptions:
-/// 
-///        If IP = 0: line number 1 is used and PYEXEC is called. 
-///        If IP < 0: line -IP is used, with status code K(-IP,2)=2
-///                   rather than 1; thus a parton system may be built
-///                   up by filling all but the last parton of the
-///                   system with IP < 0.   
-///  KF:   parton/particle flavour code (PDG code)
-///  PE:   parton/particle energy. If PE is smaller than the mass,
-///        the parton/particle is taken to be at rest.  
-///  THETA:
-///  PHI:  polar and azimuthal angle for the momentum vector of the
-///        parton/particle. 
-  
-  //std::cout<<"ACH start Py1ent 1"<<std::endl;
-  py1ent_(ip, kf, pe, theta, phi);
-  //int pyl=3; pylist_(&pyl);
-}
-
-void PythiaForDecays::DecayRhadrons(int pdgid){
-  //std::cout<<"ACH start DecayRhadrons for pdgid "<<pdgid<<std::endl;
-  int pyl=3;//level to list to output text
-
-  //only decay daughter particles from the decay that are "near" the hadron
-  Pydat1 pydat1;
-  pydat1.mstj(21)=2;
-  pydat1.mstj(22)=3;
-  pydat1.parj(72)=0.001;
-
-  //Deal with standard hadrons first...
-  if (abs(pdgid)<1000000){
-    //std::cout<<"ACH decay standard hadron for pdgid "<<pdgid<<std::endl;
-    pyjets_.k[0][0]=5;//make it decay right away!
-    //pylist_(&pyl);
-    pyexec_();
-    //pylist_(&pyl);
-    return;
-  }
-
-  //Real Rhadron stuff...
-  pdgid=abs(pdgid);
-  pyjets_.k[0][0]=6;
-  Pydat3 pydat3;
-  Pydat2 pydat2;
-  //pylist_(&pyl);
-  if ( pdgid/10000==109 || pdgid/1000==1009 || pdgid/100==10009 ){
-    //std::cout<<"ACH decay gluino rhadron .."<<std::endl;
-    int kfgl = 1000021;//gluino
-    pydat3.mdcy(pycomp_(&kfgl),1)=1;// MDCY(PYCOMP(KFGL),1)=1 //unstable gluino
-    float oldl=pydat2.pmas(pycomp_(&kfgl),4);//remember lifetime
-    pydat2.pmas(pycomp_(&kfgl),4)=0;//decay it promptly!
-    pygldc_();//gluino DECAY
-    pydat2.pmas(pycomp_(&kfgl),4)=oldl;//set lifetime back to what it was
-    //pylist_(&pyl);
-  }
-  else if ( pdgid/1000==1006 || pdgid/100==10006 ){
-    //std::cout<<"ACH decay stop rhadron .."<<std::endl;
-    Pydat3 pydat3;
-    int kfgl = 1000006;//stop
-    pydat3.mdcy(pycomp_(&kfgl),1)=1;// MDCY(PYCOMP(KFGL),1)=1 //unstable stop
-    float oldl=pydat2.pmas(pycomp_(&kfgl),4);//remember lifetime
-    pydat2.pmas(pycomp_(&kfgl),4)=0;//decay it promptly!
-    pysqdc_(&kfgl);//stop DECAY!
-    pydat2.pmas(pycomp_(&kfgl),4)=oldl;//set lifetime back to what it was
-    //pylist_(&pyl);
-  }
-  else if ( pdgid/1000==1005 || pdgid/100==10005 ){
-    //std::cout<<"ACH decay sbottom rhadron .."<<std::endl;
-    Pydat3 pydat3;
-    int kfgl = 1000005;//sbottom
-    pydat3.mdcy(pycomp_(&kfgl),1)=1;// MDCY(PYCOMP(KFGL),1)=1 //unstable sbottom
-    float oldl=pydat2.pmas(pycomp_(&kfgl),4);//remember lifetime
-    pydat2.pmas(pycomp_(&kfgl),4)=0;//decay it promptly!    
-    pysqdc_(&kfgl);//sbottom DECAY!
-    pydat2.pmas(pycomp_(&kfgl),4)=oldl;//set lifetime back to what it was
-    //pylist_(&pyl);
-  }
-  else{
-    std::cout<<"ACH ERROR !!! Can't decay unknown particle with abspdgid="<<pdgid<<std::endl;
-    pylist_(&pyl);
-    //exit(-19);
-  }
-
-}
-
-//______________________________________________________________________________
-void PythiaForDecays::ImportParticles( std::vector<G4DynamicParticle*> & particles)
-{
-  particles.clear();
-
-  int mconv=1;
-  pyhepc_(&mconv);
-
-  // Fill up the HEP Event
-  HepMC::GenEvent* evt = new HepMC::GenEvent(1,1);
-  HepMC::IO_HEPEVT hepio;
-  hepio.set_print_inconsistency_errors(0);
-  hepio.fill_next_event(evt);
-
-  for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p){
-    if (! (*p) ) continue;
-
-    /*
-    double r=0;
-    if ((*p)->production_vertex()) {
-      double x=(*p)->production_vertex()->position().x();
-      double y=(*p)->production_vertex()->position().y();
-      r=sqrt(x*x+y*y);
-    }
-    std::cout<<"ImportParticles: "<<(*p)->pdg_id()<<", vertex r "<<r<<", status "<<(*p)->status()<<std::endl;
-    */
-
-    if ( (*p)->status()!=1 ) continue; // stable only
-    //if ( std::abs((*p)->pdg_id())==12 || std::abs((*p)->pdg_id())==14 || std::abs((*p)->pdg_id())==16 ) continue; // skip neutrinos
-    G4ThreeVector momentum( (*p)->momentum().x() , (*p)->momentum().y() , (*p)->momentum().z() );
-    momentum*=1000.0;//GeV to MeV
-    const G4ParticleDefinition * particleDefinition = GetParticleDefinition( (*p)->pdg_id() );
-    if (!particleDefinition){
-      std::cout<<"WARNING: I don't know a definition for pdgid "<<(*p)->pdg_id()<<"! Skipping it..."<<std::endl;
-      continue;
-    }
-    G4DynamicParticle* dynamicParticle = new G4DynamicParticle(particleDefinition, momentum);
-    particles.push_back( dynamicParticle );
-  }
-}
-
-G4ParticleDefinition* PythiaForDecays::GetParticleDefinition(const int pdgEncoding) const
-{
-  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
-  G4ParticleDefinition* particleDefinition = 0;
-  if (pdgEncoding != 0) particleDefinition = particleTable->FindParticle(pdgEncoding);
-  return particleDefinition;
-}
-

Simulation/G4Extensions/RHadrons/CMakeLists.txt

@@ -5,33 +5,32 @@
 # Declare the package name:
 atlas_subdir( RHadrons )
 
-# Declare the package's dependencies:
-atlas_depends_on_subdirs( PUBLIC
-                          GaudiKernel
-                          PRIVATE
-                          Control/AthenaBaseComps
-			  Control/CxxUtils
-                          Simulation/G4Atlas/G4AtlasInterfaces
-                          Simulation/G4Atlas/G4AtlasTools
-                          Simulation/G4Extensions/G4ExternalDecay
-                          Simulation/G4Sim/SimHelpers )
+if( NOT GENERATIONBASE )
+    # Declare the package's dependencies:
+    atlas_depends_on_subdirs( PUBLIC
+                              GaudiKernel
+                              PRIVATE
+                              Control/AthenaBaseComps
+                              Control/CxxUtils
+                              Simulation/G4Atlas/G4AtlasInterfaces
+                              Simulation/G4Atlas/G4AtlasTools
+                              Simulation/G4Extensions/G4ExternalDecay
+                              Simulation/G4Sim/SimHelpers )
 
-# External dependencies:
-find_package( CLHEP )
-find_package( Geant4 )
-find_package( XercesC )
+    # External dependencies:
+    find_package( CLHEP )
+    find_package( Geant4 )
+    find_package( XercesC )
+    find_package( Pythia8 )
 
-# Component(s) in the package:
-atlas_add_component( RHadrons
-                     src/*.cxx
-                     src/components/*.cxx
-                     INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
-                     LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib G4ExternalDecay SimHelpers )
+    # Component(s) in the package:
+    atlas_add_component( RHadrons
+                         src/*.cxx
+                         src/components/*.cxx
+                         INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
+                         LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib G4ExternalDecay SimHelpers Pythia8_iLib )
+endif()
 
 # Install files from the package:
 atlas_install_python_modules( python/*.py )
 atlas_install_joboptions( share/*.py )
-atlas_install_generic( share/*.txt
-                       DESTINATION share
-                       EXECUTABLE )
-

Simulation/G4Extensions/RHadrons/python/RHadronsConfig.py

@@ -11,5 +11,14 @@ def getSG_StepNtupleTool(name="G4UA::SG_StepNtupleTool", **kwargs):
         #from AthenaCommon.AppMgr import theApp
         #theApp.exit(1)
         return False
+    # Get the PDG IDs for RHadrons
+    from RHadronMasses import offset_options
+    kwargs.setdefault('RHadronPDGIDList',offset_options.keys())
+    if name in simFlags.UserActionConfig.get_Value().keys():
+        for prop,value in simFlags.UserActionConfig.get_Value()[name].iteritems():
+            kwargs.setdefault(prop,value)
     from RHadrons.RHadronsConf import G4UA__SG_StepNtupleTool
     return G4UA__SG_StepNtupleTool(name, **kwargs)
+
+def getRHadronsPhysicsTool(name='RHadronsPhysicsTool', **kwargs):
+    return CfgMgr.RHadronsPhysicsTool(name,**kwargs)

Simulation/G4Extensions/RHadrons/python/RHadronsConfigDb.py

@@ -3,4 +3,4 @@
 from AthenaCommon.CfgGetter import addTool
 
 addTool("RHadrons.RHadronsConfig.getSG_StepNtupleTool", "G4UA::SG_StepNtupleTool")
-
+addTool("RHadrons.RHadronsConfig.getRHadronsPhysicsTool", "RHadronsPhysicsTool")

Simulation/G4Extensions/RHadrons/share/PhysicsConfiguration_generic.txt

-Resonant = 0.
-XsecMultiplier = 1.
-ResonanceEnergy = 200.
-Gamma = 0.1
-Amplitude = 100.
-ReggeSuppression = 0.
-HadronLifeTime = 100.
-ReggeModel = 0.
-Mixing = 0.

Simulation/G4Extensions/RHadrons/share/PhysicsConfiguration_regge.txt

-Resonant = 0.
-XsecMultiplier = 1.
-ResonanceEnergy = 200.
-Gamma = 0.1
-Amplitude = 100.
-ReggeSuppression = 0.
-HadronLifeTime = 100.
-ReggeModel = 1.
-Mixing = 1.

Simulation/G4Extensions/RHadrons/share/PhysicsConfiguration_regge_nomix.txt

-Resonant = 0.
-XsecMultiplier = 1.
-ResonanceEnergy = 200.
-Gamma = 0.1
-Amplitude = 100.
-ReggeSuppression = 0.
-HadronLifeTime = 100.
-ReggeModel = 1.
-Mixing = 0.

Simulation/G4Extensions/RHadrons/share/ProcessList_generic_stop.txt

-~T0 # neutron # ~T0 # neutron 
-~T0 # neutron # ~T0 # neutron # pi0
-~T0 # neutron # ~T0 # proton # pi-
-~T0 # neutron # ~T_dd10 # pi0
-~T0 # neutron # ~T_dd10 # pi0 # pi0
-~T0 # neutron # ~T_dd10 # pi- # pi+
-~T0 # neutron # ~T+ # neutron # pi-
-~T0 # neutron # ~T_ud0+ # pi-
-~T0 # neutron # ~T_ud0+ # pi- # pi0
-~T0 # proton # ~T0 # neutron # pi+
-~T0 # proton # ~T0 # proton 
-~T0 # proton # ~T0 # proton # pi0
-~T0 # proton # ~T+ # neutron 
-~T0 # proton # ~T+ # neutron # pi0
-~T0 # proton # ~T+ # proton # pi-
-~T0 # proton # ~T_ud0+ # pi0
-~T0 # proton # ~T_ud0+ # pi0 # pi0
-~T0 # proton # ~T_ud0+ # pi- # pi+
-~T0 # proton # ~T_uu1++ # pi-
-~T0 # proton # ~T_uu1++ # pi- # pi0
-~Tbar0 # neutron # ~Tbar0 # neutron
-~Tbar0 # neutron # ~Tbar0 # neutron # pi0
-~Tbar0 # neutron # ~Tbar0 # proton # pi-
-~Tbar0 # neutron # ~Tbar- # neutron # pi+
-~Tbar0 # neutron # ~Tbar- # proton
-~Tbar0 # neutron # ~Tbar- # proton # pi0
-~Tbar0 # proton # ~Tbar0 # neutron # pi+
-~Tbar0 # proton # ~Tbar0 # proton
-~Tbar0 # proton # ~Tbar0 # proton # pi0
-~Tbar0 # proton # ~Tbar- # proton # pi+
-~Tbar_dd10 # neutron # ~Tbar0 # pi0
-~Tbar_dd10 # neutron # ~Tbar_dd10 # neutron 
-~Tbar_dd10 # neutron # ~Tbar_dd10 # neutron # pi0 
-~Tbar_dd10 # neutron # ~Tbar- # pi0
-~Tbar_dd10 # neutron # ~Tbar_ud0- # proton
-~Tbar_dd10 # proton # ~Tbar0 # pi+
-~Tbar_dd10 # proton # ~Tbar_dd10 # proton 
-~Tbar_dd10 # proton # ~Tbar_dd10 # proton # pi0 
-~Tbar- # neutron # ~Tbar0 # neutron # pi-
-~Tbar- # neutron # ~Tbar- # neutron
-~Tbar- # neutron # ~Tbar- # neutron # pi0
-~Tbar- # neutron # ~Tbar- # proton # pi- 
-~Tbar- # proton # ~Tbar0 # neutron 
-~Tbar- # proton # ~Tbar0 # neutron # pi0
-~Tbar- # proton # ~Tbar0 # proton # pi-
-~Tbar- # proton # ~Tbar- # neutron # pi+
-~Tbar- # proton # ~Tbar- # proton
-~Tbar- # proton # ~Tbar- # proton # pi0
-~Tbar_ud0- # neutron # ~Tbar0 # pi-
-~Tbar_ud0- # neutron # ~Tbar_dd10 # neutron # pi-
-~Tbar_ud0- # neutron # ~Tbar- # pi0
-~Tbar_ud0- # neutron # ~Tbar_ud0- # neutron 
-~Tbar_ud0- # neutron # ~Tbar_ud0- # neutron # pi0 
-~Tbar_ud0- # neutron # ~Tbar_ud0- # proton # pi-
-~Tbar_ud0- # neutron # ~Tbar_uu1-- # neutron # pi+
-~Tbar_ud0- # neutron # ~Tbar_uu1-- # proton
-~Tbar_ud0- # neutron # ~Tbar_uu1-- # proton # pi0
-~Tbar_ud0- # proton # ~Tbar0 # pi0
-~Tbar_ud0- # proton # ~Tbar_dd10 # neutron
-~Tbar_ud0- # proton # ~Tbar_dd10 # neutron # pi0
-~Tbar_ud0- # proton # ~Tbar_dd10 # proton # pi-
-~Tbar_ud0- # proton # ~Tbar- # pi+
-~Tbar_ud0- # proton # ~Tbar_ud0- # neutron # pi+
-~Tbar_ud0- # proton # ~Tbar_ud0- # proton 
-~Tbar_ud0- # proton # ~Tbar_ud0- # proton # pi0 
-~Tbar_ud0- # proton # ~Tbar_uu1-- # proton # pi+
-~Tbar_uu1-- # neutron # ~Tbar- # pi-
-~Tbar_uu1-- # neutron # ~Tbar_ud0- # neutron # pi-
-~Tbar_uu1-- # neutron # ~Tbar_uu1-- # neutron 
-~Tbar_uu1-- # neutron # ~Tbar_uu1-- # neutron # pi0 
-~Tbar_uu1-- # neutron # ~Tbar_uu1-- # proton # pi-
-~Tbar_uu1-- # proton # ~Tbar0 # pi-
-~Tbar_uu1-- # proton # ~Tbar_dd10 # neutron # pi-
-~Tbar_uu1-- # proton # ~Tbar- # pi0
-~Tbar_uu1-- # proton # ~Tbar_ud0- # neutron
-~Tbar_uu1-- # proton # ~Tbar_ud0- # neutron # pi0
-~Tbar_uu1-- # proton # ~Tbar_ud0- # proton # pi-
-~Tbar_uu1-- # proton # ~Tbar_uu1-- # neutron # pi+
-~Tbar_uu1-- # proton # ~Tbar_uu1-- # proton 
-~Tbar_uu1-- # proton # ~Tbar_uu1-- # proton # pi0 
-~T_dd10 # neutron # ~T_dd10 # neutron
-~T_dd10 # neutron # ~T_dd10 # neutron # pi0
-~T_dd10 # neutron # ~T_dd10 # proton # pi-
-~T_dd10 # neutron # ~T_ud0+ # neutron # pi-
-~T_dd10 # proton # ~T_dd10 # neutron # pi+
-~T_dd10 # proton # ~T_dd10 # proton
-~T_dd10 # proton # ~T_dd10 # proton # pi0
-~T_dd10 # proton # ~T_ud0+ # neutron
-~T_dd10 # proton # ~T_ud0+ # neutron # pi0
-~T_dd10 # proton # ~T_ud0+ # proton # pi- 
-~T_dd10 # proton # ~T_uu1++ # neutron # pi-
-~T+ # neutron # ~T0 # neutron # pi-
-~T+ # neutron # ~T0 # proton 
-~T+ # neutron # ~T0 # proton # pi0
-~T+ # neutron # ~T_dd10 # pi+
-~T+ # neutron # ~T_dd10 # pi+ # pi0
-~T+ # neutron # ~T+ # neutron 
-~T+ # neutron # ~T+ # neutron # pi0
-~T+ # neutron # ~T+ # proton # pi-
-~T+ # neutron # ~T_ud0+ # pi0
-~T+ # neutron # ~T_ud0+ # pi0 # pi0
-~T+ # neutron # ~T_ud0+ # pi- # pi+
-~T+ # proton # ~T0 # proton # pi+
-~T+ # proton # ~T+ # neutron # pi+
-~T+ # proton # ~T+ # proton 
-~T+ # proton # ~T+ # proton # pi0
-~T+ # proton # ~T_ud0+ # pi+
-~T+ # proton # ~T_ud0+ # pi+ # pi0
-~T+ # proton # ~T_uu1++ # pi0
-~T+ # proton # ~T_uu1++ # pi0 # pi0
-~T+ # proton # ~T_uu1++ # pi+ # pi-
-~T_ud0+ # neutron # ~T_dd10 # neutron # pi+
-~T_ud0+ # neutron # ~T_dd10 # proton
-~T_ud0+ # neutron # ~T_dd10 # proton # pi0
-~T_ud0+ # neutron # ~T_ud0+ # neutron
-~T_ud0+ # neutron # ~T_ud0+ # neutron # pi0
-~T_ud0+ # neutron # ~T_ud0+ # proton # pi-
-~T_ud0+ # proton # ~T_dd10 # proton # pi+
-~T_ud0+ # proton # ~T_ud0+ # neutron # pi+
-~T_ud0+ # proton # ~T_ud0+ # proton
-~T_ud0+ # proton # ~T_ud0+ # proton # pi0
-~T_ud0+ # proton # ~T_uu1++ # neutron
-~T_ud0+ # proton # ~T_uu1++ # neutron # pi0
-~T_ud0+ # proton # ~T_uu1++ # proton # pi-
-~T_uu1++ # neutron # ~T_ud0+ # neutron # pi+
-~T_uu1++ # neutron # ~T_ud0+ # proton
-~T_uu1++ # neutron # ~T_ud0+ # proton # pi0
-~T_uu1++ # neutron # ~T_uu1++ # neutron
-~T_uu1++ # neutron # ~T_uu1++ # neutron # pi0
-~T_uu1++ # neutron # ~T_uu1++ # proton # pi-
-~T_uu1++ # proton # ~T_ud0+ # proton # pi+
-~T_uu1++ # proton # ~T_uu1++ # neutron # pi+
-~T_uu1++ # proton # ~T_uu1++ # proton
-~T_uu1++ # proton # ~T_uu1++ # proton # pi0
-~T_s+ # neutron # ~T0 # neutron # pi-
-~T_s+ # neutron # ~T0 # proton 
-~T_s+ # neutron # ~T0 # proton # pi0
-~T_s+ # neutron # ~T_dd10 # pi+
-~T_s+ # neutron # ~T_dd10 # pi+ # pi0
-~T_s+ # neutron # ~T+ # neutron 
-~T_s+ # neutron # ~T+ # neutron # pi0
-~T_s+ # neutron # ~T+ # proton # pi-
-~T_s+ # neutron # ~T_ud0+ # pi0
-~T_s+ # neutron # ~T_ud0+ # pi0 # pi0
-~T_s+ # neutron # ~T_ud0+ # pi- # pi+
-~T_s+ # proton # ~T0 # proton # pi+
-~T_s+ # proton # ~T+ # neutron # pi+
-~T_s+ # proton # ~T+ # proton 
-~T_s+ # proton # ~T+ # proton # pi0
-~T_s+ # proton # ~T_ud0+ # pi+
-~T_s+ # proton # ~T_ud0+ # pi+ # pi0
-~T_s+ # proton # ~T_uu1++ # pi0
-~T_s+ # proton # ~T_uu1++ # pi0 # pi0
-~T_s+ # proton # ~T_uu1++ # pi+ # pi-
-~T_sd00 # neutron # ~T_dd10 # neutron
-~T_sd00 # neutron # ~T_dd10 # neutron # pi0
-~T_sd00 # neutron # ~T_dd10 # proton # pi-
-~T_sd00 # neutron # ~T_ud0+ # neutron # pi-
-~T_sd00 # proton # ~T_dd10 # neutron # pi+
-~T_sd00 # proton # ~T_dd10 # proton
-~T_sd00 # proton # ~T_dd10 # proton # pi0
-~T_sd00 # proton # ~T_ud0+ # neutron
-~T_sd00 # proton # ~T_ud0+ # neutron # pi0
-~T_sd00 # proton # ~T_ud0+ # proton # pi- 
-~T_sd00 # proton # ~T_uu1++ # neutron # pi-
-~T_sd10 # neutron # ~T_dd10 # neutron
-~T_sd10 # neutron # ~T_dd10 # neutron # pi0
-~T_sd10 # neutron # ~T_dd10 # proton # pi-
-~T_sd10 # neutron # ~T_ud0+ # neutron # pi-
-~T_sd10 # proton # ~T_dd10 # neutron # pi+
-~T_sd10 # proton # ~T_dd10 # proton
-~T_sd10 # proton # ~T_dd10 # proton # pi0
-~T_sd10 # proton # ~T_ud0+ # neutron
-~T_sd10 # proton # ~T_ud0+ # neutron # pi0
-~T_sd10 # proton # ~T_ud0+ # proton # pi- 
-~T_sd10 # proton # ~T_uu1++ # neutron # pi-
-~T_su0+ # neutron # ~T_dd10 # neutron # pi+
-~T_su0+ # neutron # ~T_dd10 # proton
-~T_su0+ # neutron # ~T_dd10 # proton # pi0
-~T_su0+ # neutron # ~T_ud0+ # neutron
-~T_su0+ # neutron # ~T_ud0+ # neutron # pi0
-~T_su0+ # neutron # ~T_ud0+ # proton # pi-
-~T_su0+ # proton # ~T_dd10 # proton # pi+
-~T_su0+ # proton # ~T_ud0+ # neutron # pi+
-~T_su0+ # proton # ~T_ud0+ # proton
-~T_su0+ # proton # ~T_ud0+ # proton # pi0
-~T_su0+ # proton # ~T_uu1++ # neutron
-~T_su0+ # proton # ~T_uu1++ # neutron # pi0
-~T_su0+ # proton # ~T_uu1++ # proton # pi-
-~T_su1+ # neutron # ~T_dd10 # neutron # pi+
-~T_su1+ # neutron # ~T_dd10 # proton
-~T_su1+ # neutron # ~T_dd10 # proton # pi0
-~T_su1+ # neutron # ~T_ud0+ # neutron
-~T_su1+ # neutron # ~T_ud0+ # neutron # pi0
-~T_su1+ # neutron # ~T_ud0+ # proton # pi-
-~T_su1+ # proton # ~T_dd10 # proton # pi+
-~T_su1+ # proton # ~T_ud0+ # neutron # pi+
-~T_su1+ # proton # ~T_ud0+ # proton
-~T_su1+ # proton # ~T_ud0+ # proton # pi0
-~T_su1+ # proton # ~T_uu1++ # neutron
-~T_su1+ # proton # ~T_uu1++ # neutron # pi0
-~T_su1+ # proton # ~T_uu1++ # proton # pi-
-~T_ss10 # neutron # ~T_dd10 # neutron
-~T_ss10 # neutron # ~T_dd10 # neutron # pi0
-~T_ss10 # neutron # ~T_dd10 # proton # pi-
-~T_ss10 # neutron # ~T_ud0+ # neutron # pi-
-~T_ss10 # proton # ~T_dd10 # neutron # pi+
-~T_ss10 # proton # ~T_dd10 # proton
-~T_ss10 # proton # ~T_dd10 # proton # pi0
-~T_ss10 # proton # ~T_ud0+ # neutron
-~T_ss10 # proton # ~T_ud0+ # neutron # pi0
-~T_ss10 # proton # ~T_ud0+ # proton # pi- 
-~T_ss10 # proton # ~T_uu1++ # neutron # pi-
-~Tbar_s- # neutron # ~Tbar0 # neutron # pi-
-~Tbar_s- # neutron # ~Tbar- # neutron
-~Tbar_s- # neutron # ~Tbar- # neutron # pi0
-~Tbar_s- # neutron # ~Tbar- # proton # pi- 
-~Tbar_s- # proton # ~Tbar0 # neutron 
-~Tbar_s- # proton # ~Tbar0 # neutron # pi0
-~Tbar_s- # proton # ~Tbar0 # proton # pi-
-~Tbar_s- # proton # ~Tbar- # neutron # pi+
-~Tbar_s- # proton # ~Tbar- # proton
-~Tbar_s- # proton # ~Tbar- # proton # pi0
-~Tbar_sd00 # neutron # ~Tbar0 # pi0
-~Tbar_sd00 # neutron # ~Tbar_dd10 # neutron 
-~Tbar_sd00 # neutron # ~Tbar_dd10 # neutron # pi0 
-~Tbar_sd00 # neutron # ~Tbar- # pi0
-~Tbar_sd00 # neutron # ~Tbar_ud0- # proton
-~Tbar_sd00 # proton # ~Tbar0 # pi+
-~Tbar_sd00 # proton # ~Tbar_dd10 # proton 
-~Tbar_sd00 # proton # ~Tbar_dd10 # proton # pi0 
-~Tbar_sd10 # neutron # ~Tbar0 # pi0
-~Tbar_sd10 # neutron # ~Tbar_dd10 # neutron 
-~Tbar_sd10 # neutron # ~Tbar_dd10 # neutron # pi0 
-~Tbar_sd10 # neutron # ~Tbar- # pi0
-~Tbar_sd10 # neutron # ~Tbar_ud0- # proton
-~Tbar_sd10 # proton # ~Tbar0 # pi+
-~Tbar_sd10 # proton # ~Tbar_dd10 # proton 
-~Tbar_sd10 # proton # ~Tbar_dd10 # proton # pi0 
-~Tbar_su0- # neutron # ~Tbar0 # pi-
-~Tbar_su0- # neutron # ~Tbar_dd10 # neutron # pi-
-~Tbar_su0- # neutron # ~Tbar- # pi0
-~Tbar_su0- # neutron # ~Tbar_ud0- # neutron 
-~Tbar_su0- # neutron # ~Tbar_ud0- # neutron # pi0 
-~Tbar_su0- # neutron # ~Tbar_ud0- # proton # pi-
-~Tbar_su0- # neutron # ~Tbar_uu1-- # neutron # pi+
-~Tbar_su0- # neutron # ~Tbar_uu1-- # proton
-~Tbar_su0- # neutron # ~Tbar_uu1-- # proton # pi0
-~Tbar_su0- # proton # ~Tbar0 # pi0
-~Tbar_su0- # proton # ~Tbar_dd10 # neutron
-~Tbar_su0- # proton # ~Tbar_dd10 # neutron # pi0
-~Tbar_su0- # proton # ~Tbar_dd10 # proton # pi-
-~Tbar_su0- # proton # ~Tbar- # pi+
-~Tbar_su0- # proton # ~Tbar_ud0- # neutron # pi+
-~Tbar_su0- # proton # ~Tbar_ud0- # proton 
-~Tbar_su0- # proton # ~Tbar_ud0- # proton # pi0 
-~Tbar_su0- # proton # ~Tbar_uu1-- # proton # pi+
-~Tbar_su1- # neutron # ~Tbar0 # pi-
-~Tbar_su1- # neutron # ~Tbar_dd10 # neutron # pi-
-~Tbar_su1- # neutron # ~Tbar- # pi0
-~Tbar_su1- # neutron # ~Tbar_ud0- # neutron 
-~Tbar_su1- # neutron # ~Tbar_ud0- # neutron # pi0 
-~Tbar_su1- # neutron # ~Tbar_ud0- # proton # pi-
-~Tbar_su1- # neutron # ~Tbar_uu1-- # neutron # pi+
-~Tbar_su1- # neutron # ~Tbar_uu1-- # proton
-~Tbar_su1- # neutron # ~Tbar_uu1-- # proton # pi0
-~Tbar_su1- # proton # ~Tbar0 # pi0
-~Tbar_su1- # proton # ~Tbar_dd10 # neutron
-~Tbar_su1- # proton # ~Tbar_dd10 # neutron # pi0
-~Tbar_su1- # proton # ~Tbar_dd10 # proton # pi-
-~Tbar_su1- # proton # ~Tbar- # pi+
-~Tbar_su1- # proton # ~Tbar_ud0- # neutron # pi+
-~Tbar_su1- # proton # ~Tbar_ud0- # proton 
-~Tbar_su1- # proton # ~Tbar_ud0- # proton # pi0 
-~Tbar_su1- # proton # ~Tbar_uu1-- # proton # pi+
-~Tbar_ss10 # neutron # ~Tbar0 # pi0
-~Tbar_ss10 # neutron # ~Tbar_dd10 # neutron 
-~Tbar_ss10 # neutron # ~Tbar_dd10 # neutron # pi0 
-~Tbar_ss10 # neutron # ~Tbar- # pi0
-~Tbar_ss10 # neutron # ~Tbar_ud0- # proton
-~Tbar_ss10 # proton # ~Tbar0 # pi+
-~Tbar_ss10 # proton # ~Tbar_dd10 # proton 
-~Tbar_ss10 # proton # ~Tbar_dd10 # proton # pi0