From 3dc8cf84015ec89d02df5b7b7dd4ad09062a12d4 Mon Sep 17 00:00:00 2001
From: James William Monk <james.william.monk@cern.ch>
Date: Wed, 27 Nov 2013 15:51:06 +0100
Subject: [PATCH] tagging for Pythia 8.180 compatibility
(Pythia8_i-00-09-10-01)
---
Generators/Pythia8_i/Pythia8_i/Pythia8_i.h | 129 +++++
.../Pythia8_i/Pythia8_i/UserHooksFactory.h | 69 +++
.../Pythia8_i/Pythia8_i/UserProcessFactory.h | 70 +++
.../Pythia8_i/UserResonanceFactory.h | 72 +++
Generators/Pythia8_i/cmt/requirements | 35 ++
.../Pythia8_i/share/Powheg.ZMu.MC11.events | 208 ++++++++
.../Pythia8_i/share/enhanceMPI_example.py | 12 +
.../Pythia8_i/share/jobOptions.pythia8.py | 12 +
Generators/Pythia8_i/share/lhef_example.py | 32 ++
.../Pythia8_i/share/qcdVetoedShowerExample.py | 36 ++
.../Pythia8_i/share/suppressMPI_example.py | 12 +
.../Pythia8_i/share/userProcess_example.py | 36 ++
.../Pythia8_i/share/vetoedISR_example.py | 35 ++
Generators/Pythia8_i/src/Pythia8_i.cxx | 410 ++++++++++++++
.../Pythia8_i/src/UserHooks/EnhanceMPI.cxx | 104 ++++
.../Pythia8_i/src/UserHooks/GravFlat.cxx | 117 ++++
.../src/UserHooks/ISRVetoedShower.cxx | 173 ++++++
.../src/UserHooks/PTRelVetoedShower.cxx | 219 ++++++++
.../src/UserHooks/PoWHEGVetoedShower.cxx | 180 +++++++
.../src/UserHooks/QCDVetoedShower.cxx | 138 +++++
.../Pythia8_i/src/UserHooks/SuppressMPI.cxx | 85 +++
.../Pythia8_i/src/UserHooks/UserHooksUtils.h | 116 ++++
.../Pythia8_i/src/UserHooks/VetoedShower.cxx | 151 ++++++
.../src/UserHooks/WZVetoedShower.cxx | 128 +++++
.../Pythia8_i/src/UserHooks/WprimeFlat.cxx | 60 +++
.../Pythia8_i/src/UserHooks/WprimeWZFlat.cxx | 68 +++
Generators/Pythia8_i/src/UserHooks/main31.cxx | 505 ++++++++++++++++++
Generators/Pythia8_i/src/UserHooksFactory.cxx | 27 +
.../Pythia8_i/src/UserProcessFactory.cxx | 28 +
.../src/UserProcesses/Sigma2qqbar2emu.cxx | 127 +++++
.../UserProcesses/Sigma2qqbar2lStarlBar.cxx | 187 +++++++
.../Sigma2qqbar2lStarlStarBar.cxx | 186 +++++++
.../Pythia8_i/src/UserResonanceFactory.cxx | 33 ++
.../src/UserResonances/ResonanceExcitedCI.cxx | 114 ++++
.../src/components/Pythia8_i_entries.cxx | 11 +
.../src/components/Pythia8_i_load.cxx | 7 +
36 files changed, 3932 insertions(+)
create mode 100644 Generators/Pythia8_i/Pythia8_i/Pythia8_i.h
create mode 100644 Generators/Pythia8_i/Pythia8_i/UserHooksFactory.h
create mode 100644 Generators/Pythia8_i/Pythia8_i/UserProcessFactory.h
create mode 100644 Generators/Pythia8_i/Pythia8_i/UserResonanceFactory.h
create mode 100644 Generators/Pythia8_i/cmt/requirements
create mode 100644 Generators/Pythia8_i/share/Powheg.ZMu.MC11.events
create mode 100644 Generators/Pythia8_i/share/enhanceMPI_example.py
create mode 100644 Generators/Pythia8_i/share/jobOptions.pythia8.py
create mode 100644 Generators/Pythia8_i/share/lhef_example.py
create mode 100644 Generators/Pythia8_i/share/qcdVetoedShowerExample.py
create mode 100644 Generators/Pythia8_i/share/suppressMPI_example.py
create mode 100644 Generators/Pythia8_i/share/userProcess_example.py
create mode 100644 Generators/Pythia8_i/share/vetoedISR_example.py
create mode 100644 Generators/Pythia8_i/src/Pythia8_i.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/EnhanceMPI.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/GravFlat.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/ISRVetoedShower.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/PTRelVetoedShower.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/PoWHEGVetoedShower.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/QCDVetoedShower.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/SuppressMPI.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h
create mode 100644 Generators/Pythia8_i/src/UserHooks/VetoedShower.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/WZVetoedShower.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/WprimeFlat.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/WprimeWZFlat.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooks/main31.cxx
create mode 100644 Generators/Pythia8_i/src/UserHooksFactory.cxx
create mode 100644 Generators/Pythia8_i/src/UserProcessFactory.cxx
create mode 100644 Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2emu.cxx
create mode 100644 Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlBar.cxx
create mode 100644 Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlStarBar.cxx
create mode 100644 Generators/Pythia8_i/src/UserResonanceFactory.cxx
create mode 100644 Generators/Pythia8_i/src/UserResonances/ResonanceExcitedCI.cxx
create mode 100644 Generators/Pythia8_i/src/components/Pythia8_i_entries.cxx
create mode 100644 Generators/Pythia8_i/src/components/Pythia8_i_load.cxx
diff --git a/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h b/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h
new file mode 100644
index 00000000000..6fadee22fc9
--- /dev/null
+++ b/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h
@@ -0,0 +1,129 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GENERATOR_PYTHIA8_H
+#define GENERATOR_PYTHIA8_H
+
+#include "GeneratorModules/GenModule.h"
+
+#include "Pythia8/Pythia.h"
+#include "Pythia8/Pythia8ToHepMC.h"
+
+// calls to fortran routines
+#include "CLHEP/Random/RandFlat.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+
+/**
+ * Author: James Monk (jmonk@hep.ucl.ac.uk)
+*/
+
+using std::string;
+
+class IAtRndmGenSvc;
+
+namespace Pythia8{
+ class Sigma2Process;
+}
+
+
+class customRndm : public Pythia8::RndmEngine {
+public:
+
+ // Constructor.
+ customRndm():
+ m_engine(0), m_RndmGenSvc(0) {}
+
+ // Routine for generating a random number.
+ inline double flat(){
+ m_engine=m_RndmGenSvc->GetEngine(m_stream);
+ return CLHEP::RandFlat::shoot(m_engine);
+ };
+
+ // Initialisation Routine
+ inline void init(IAtRndmGenSvc &RndmGenSvc, std::string RndmStream){m_RndmGenSvc=&RndmGenSvc; m_stream=RndmStream;};
+ std::string m_stream;
+
+private:
+ CLHEP::HepRandomEngine* m_engine;
+ IAtRndmGenSvc* m_RndmGenSvc;
+};
+
+
+//namespace Generator{
+class Pythia8_i: public GenModule{
+
+public:
+ Pythia8_i(const string &name, ISvcLocator *pSvcLocator);
+
+ ~Pythia8_i();
+
+ virtual StatusCode genInitialize();
+ virtual StatusCode callGenerator();
+ virtual StatusCode fillEvt(GenEvent *evt);
+ virtual StatusCode genFinalize();
+
+ double pythiaVersion()const;
+
+ static std::string pythia_stream;
+
+protected:
+
+ // make these protected so that Pythia8B can access them
+ Pythia8::Pythia m_pythia;
+ HepMC::Pythia8ToHepMC m_pythiaToHepMC;
+
+private:
+
+ static std::string xmlpath();
+
+// StoreGateSvc* m_sGateService;
+// string m_mcEventKey;
+
+ int m_internal_event_number;
+
+ double m_version;
+
+ std::vector<std::string> m_commands;
+
+ enum PDGID {PROTON=2212, ANTIPROTON=-2212, ELECTRON=11, POSITRON=-11, INVALID=0};
+
+ double m_collisionEnergy;
+ bool m_useRndmGenSvc;
+ customRndm *m_atlasRndmEngine;
+
+
+ std::string m_beam1;
+ std::string m_beam2;
+
+ std::string m_lheFile;
+
+ unsigned int m_maxFailures;
+ unsigned int m_failureCount;
+
+ std::map<std::string, PDGID> m_particleIDs;
+
+ std::vector<long int> m_seeds;
+
+ std::string m_userProcess;
+
+ // ptr to possible user process
+ Pythia8::Sigma2Process *m_procPtr;
+
+ std::string m_userHook;
+
+ Pythia8::UserHooks *m_userHookPtr;
+
+ std::string m_userResonances;
+
+ vector<Pythia8::ResonanceWidths*> m_userResonancePtrs;
+
+ bool m_useLHAPDF;
+
+ std::string m_particleDataFile;
+ std::string m_outputParticleDataFile;
+
+};
+
+#endif
diff --git a/Generators/Pythia8_i/Pythia8_i/UserHooksFactory.h b/Generators/Pythia8_i/Pythia8_i/UserHooksFactory.h
new file mode 100644
index 00000000000..a9cebaaf0ef
--- /dev/null
+++ b/Generators/Pythia8_i/Pythia8_i/UserHooksFactory.h
@@ -0,0 +1,69 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GENERATOR_PYTHIA8_USER_HOOKS_FACTORY_H
+#define GENERATOR_PYTHIA8_USER_HOOKS_FACTORY_H
+
+#include "Pythia8/UserHooks.h"
+
+#include <string>
+
+namespace Pythia8_UserHooks{
+
+ using Pythia8::UserHooks;
+ using std::string;
+
+
+ class UserHooksFactory{
+
+ public:
+
+ static UserHooks* create(const string &hookName);
+
+ private:
+
+ UserHooksFactory(){};
+
+ class ICreator{
+ public:
+ virtual UserHooks *create() const = 0;
+ virtual ~ICreator(){};
+ };
+
+ public:
+
+ template <class T>
+ class Creator: public ICreator{
+
+ public:
+
+ Creator(const string &name){
+ m_name = name;
+ UserHooksFactory::s_creators()[name] = this;
+ }
+
+ ~Creator(){
+ if(s_creators()[m_name] == this){
+ s_creators().erase(m_name);
+ }
+ }
+
+ UserHooks *create()const{
+ return new T;
+ }
+
+ private:
+
+ string m_name;
+ };
+
+ private:
+
+ static map<string, const ICreator*> &s_creators();
+
+ };
+}
+
+
+#endif
diff --git a/Generators/Pythia8_i/Pythia8_i/UserProcessFactory.h b/Generators/Pythia8_i/Pythia8_i/UserProcessFactory.h
new file mode 100644
index 00000000000..fe3bb8d825c
--- /dev/null
+++ b/Generators/Pythia8_i/Pythia8_i/UserProcessFactory.h
@@ -0,0 +1,70 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GENERATOR_PYTHIA8_USER_PROCESS_FACTORY_H
+#define GENERATOR_PYTHIA8_USER_PROCESS_FACTORY_H
+
+#include "Pythia8/SigmaProcess.h"
+
+#include <boost/smart_ptr.hpp>
+
+#include <map>
+#include <string>
+
+namespace Pythia8_UserProcess{
+
+ using Pythia8::Sigma2Process;
+
+ using std::string;
+ using std::map;
+
+ class UserProcessFactory{
+
+ public:
+
+ static Sigma2Process* create(const string &procName);
+
+ private:
+
+ UserProcessFactory(){};
+
+ class ICreator{
+ public:
+ virtual Sigma2Process *create() const = 0;
+ virtual ~ICreator(){};
+ };
+
+ public:
+
+ template <class T>
+ class Creator: public ICreator{
+
+ public:
+ Creator(const string &name){
+ m_name = name;
+ UserProcessFactory::s_creators()[name] = this;
+ }
+
+ ~Creator(){
+ if(s_creators()[m_name] == this){
+ s_creators().erase(m_name);
+ }
+ }
+
+ Sigma2Process *create()const{
+ return new T;
+ }
+
+ private:
+
+ string m_name;
+
+ };
+
+ private:
+ static map<string, const ICreator*> &s_creators();
+
+ };
+}
+#endif
diff --git a/Generators/Pythia8_i/Pythia8_i/UserResonanceFactory.h b/Generators/Pythia8_i/Pythia8_i/UserResonanceFactory.h
new file mode 100644
index 00000000000..32736c4bb00
--- /dev/null
+++ b/Generators/Pythia8_i/Pythia8_i/UserResonanceFactory.h
@@ -0,0 +1,72 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GENERATOR_PYTHIA8_USER_RESONANCE_FACTORY_H
+#define GENERATOR_PYTHIA8_USER_RESONANCE_FACTORY_H
+
+#include "Pythia8/ResonanceWidths.h"
+
+#include <string>
+
+namespace Pythia8_UserResonance{
+
+ using Pythia8::ResonanceWidths;
+ using std::string;
+ using std::map;
+
+ class UserResonanceFactory{
+
+ public:
+
+ /**
+ * Call this with the name of the ResonanceWidth and PDG ID to which it will be applied
+ * e.g. create("MyResonance", 23) will return a MyResonance instance that will be applied to the Z
+ *
+ */
+ static ResonanceWidths* create(const string &name, int pdgid);
+
+ private:
+
+ UserResonanceFactory(){};
+
+ class ICreator{
+ public:
+ virtual ResonanceWidths *create(int idResIn)const = 0;
+ virtual ~ICreator(){};
+ };
+
+ public:
+
+ template <class T>
+ class Creator: public ICreator{
+
+ public:
+ Creator(const string &name){
+ m_name = name;
+ UserResonanceFactory::s_creators()[name] = this;
+ }
+
+ ~Creator(){
+ if(s_creators()[m_name] == this){
+ s_creators().erase(m_name);
+ }
+ }
+
+ ResonanceWidths *create(int idResIn)const{
+ return new T(idResIn);
+ }
+
+ private:
+
+ string m_name;
+ };
+
+ private:
+
+ static map<string, const ICreator*> &s_creators();
+
+ };
+}
+
+#endif
diff --git a/Generators/Pythia8_i/cmt/requirements b/Generators/Pythia8_i/cmt/requirements
new file mode 100644
index 00000000000..b9b83a15d70
--- /dev/null
+++ b/Generators/Pythia8_i/cmt/requirements
@@ -0,0 +1,35 @@
+package Pythia8_i
+
+author James Monk <jmonk@hep.ucl.ac.uk>
+
+use AtlasPolicy AtlasPolicy-*
+use AtlasBoost AtlasBoost-* External
+use Lhapdf Lhapdf-* External
+use Pythia8 Pythia8-* External
+use GeneratorModules GeneratorModules-* Generators
+use AthenaKernel AthenaKernel-* Control
+use AtlasCLHEP AtlasCLHEP-* External
+
+private
+
+use GaudiInterface GaudiInterface-* External
+###use AtlasHepMC AtlasHepMC-* External
+use GeneratorObjects GeneratorObjects-* Generators
+
+macro_append Pythia8_i_pp_cppflags ' -DPY8VERSION=\"$(Pythia8_version)\" '
+macro_append fflags "" Linux ""
+
+# The dependencies of the various generator packages need to be sorted out
+# so that they work in full asNeeded/noUndefined mode. Until that time, this
+# package explicitly sets the allowUndefined tag
+#apply_tag allowUndefined
+apply_tag notAsNeeded
+end_private
+
+macro_append pythia8_ifiles 'Pythia8_i.cxx UserProcessFactory.cxx UserHooksFactory.cxx UserResonanceFactory.cxx'
+macro_append UserProcessFiles 'UserProcesses/Sigma2qqbar2emu.cxx UserProcesses/Sigma2qqbar2lStarlBar.cxx UserProcesses/Sigma2qqbar2lStarlStarBar.cxx'
+macro_append UserHookFiles 'UserHooks/WZVetoedShower.cxx UserHooks/QCDVetoedShower.cxx UserHooks/PoWHEGVetoedShower.cxx UserHooks/GravFlat.cxx UserHooks/SuppressMPI.cxx UserHooks/EnhanceMPI.cxx UserHooks/ISRVetoedShower.cxx UserHooks/PTRelVetoedShower.cxx UserHooks/WprimeFlat.cxx UserHooks/WprimeWZFlat.cxx UserHooks/main31.cxx'
+macro_append UserResonanceFiles 'UserResonances/ResonanceExcitedCI.cxx'
+
+apply_pattern named_dual_use_library library="Pythia8_i" files="$(pythia8_ifiles) $(UserProcessFiles) $(UserHookFiles) $(UserResonanceFiles)"
+apply_pattern declare_joboptions files="*.py"
diff --git a/Generators/Pythia8_i/share/Powheg.ZMu.MC11.events b/Generators/Pythia8_i/share/Powheg.ZMu.MC11.events
new file mode 100644
index 00000000000..43e4e334f74
--- /dev/null
+++ b/Generators/Pythia8_i/share/Powheg.ZMu.MC11.events
@@ -0,0 +1,208 @@
+<LesHouchesEvents version="1.0">
+<!--
+file generated with POWHEG-BOX version 1.0
+Input file powheg.input contained:
+!Single vector boson production parameters
+idvecbos 23 ! PDG code for vector boson to be produced ( W+:24 W-:-24 )
+vdecaymode 2 ! PDG code for charged decay product of the vector boson (11:e-; -11:e+; ...)
+
+numevts 10 ! number of events to be generated
+ih1 1 ! hadron 1 (1 for protons, -1 for antiprotons)
+ih2 1 ! hadron 2 (1 for protons, -1 for antiprotons)
+ndns1 131 ! pdf set for hadron 1 (mlm numbering)
+ndns2 131 ! pdf set for hadron 2 (mlm numbering)
+ebeam1 3500d0 ! energy of beam 1
+ebeam2 3500d0 ! energy of beam 2
+
+! To be set only if using LHA pdfs
+lhans1 10800 ! pdf set for hadron 1 (LHA numbering)
+lhans2 10800 ! pdf set for hadron 2 (LHA numbering)
+! To be set only if using different pdf sets for the two incoming hadrons
+! QCDLambda5 0.25 ! for not equal pdf sets
+
+! Parameters to allow or not the use of stored data
+use-old-grid 1 ! if 1 use old grid if file pwggrids.dat is present (<> 1 regenerate)
+use-old-ubound 1 ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; <>
+
+ncall1 120000 ! number of calls for initializing the integration grid
+itmx1 5 ! number of iterations for initializing the integration grid
+ncall2 250000 ! number of calls for computing the integral and finding upper bound
+itmx2 5 ! number of iterations for computing the integral and finding upper bound
+foldcsi 1 ! number of folds on csi integration
+foldy 1 ! number of folds on y integration
+foldphi 1 ! number of folds on phi integration
+nubound 20000 ! number of bbarra calls to setup norm of upper bounding function
+icsimax 1 ! <= 100, number of csi subdivision when computing the upper bounds
+iymax 1 ! <= 100, number of y subdivision when computing the upper bounds
+xupbound 2d0 ! increase upper bound for radiation generation
+
+! OPTIONAL PARAMETERS
+
+!SM parameters
+
+Zmass 91.1876 ! Z mass in GeV
+Zwidth 2.4952 ! Z width in GeV
+sthw2 0.23113 ! sin**2 theta w
+!alphaem 0.0072973 ! em coupling 1/137
+alphaem 0.00781653 ! em coupling 1/127
+Wmass 80.399 ! W mass in GeV
+Wwidth 2.085 ! W width in GeV
+
+CKM_Vud 0.97428
+CKM_Vus 0.2253
+CKM_Vub 0.00347
+CKM_Vcd 0.2252
+CKM_Vcs 0.97345
+CKM_Vcb 0.0410
+CKM_Vtd 0.00862
+CKM_Vts 0.0403
+CKM_Vtb 0.999152
+
+masswindow_low 15. ! M Z > Zmass - masswindow low * Zwidth
+masswindow_high 5000. ! M Z < Zmass + masswindow high * Zwidth
+
+runningscale 1 ! choice for ren and fac scales in Bbar integration
+ !0: fixed scale M Z
+ !1: running scale inv mass Z
+ !2: running scale transverse mass Z
+renscfact 1d0 ! (default 1d0) ren scale factor: muren = muref * renscfact
+facscfact 1d0 ! (default 1d0) fac scale factor: mufact = muref * facscfact
+# withdamp 1 ! (default 0, do not use) use Born-zero damping factor
+pdfreweight 1 ! (default 0) write extra pdf infos on LHEF
+
+! RANDOM NUMBER SEEDS
+
+iseed 1 ! initialize random number sequence
+rand1 -1 ! initialize random number sequence
+rand2 -1 ! initialize random number sequence
+
+
+#ptsupp 0d0 ! (default 0d0) mass param for Born suppression factor (generation cut) If < 0 su
+#bornonly 0 ! (default 0) if 1 do Born only
+#smartsig 1 ! (default 1) remember equal amplitudes (0 do not remember)
+#withsubtr 0 ! (default 1) subtract real counterterms (0 do not subtract)
+#ptsqmin 0.8 ! (default 0.8 GeV) minimum pt for generation of radiation
+#charmthr 1.5 ! (default 1.5 GeV) charm treshold for gluon splitting
+#bottomthr 5.0 ! (default 5.0 GeV) bottom treshold for gluon splitting
+#testplots 1 ! (default 0, do not) do NLO and PWHG distributions
+#hfact 100d0 ! (default no dumping factor) dump factor for high-pt radiation: > 0 dumpfac=h**2
+#testsuda 1 ! (default 0, do not test) test Sudakov form factor
+#radregion 1 ! (default all regions) only generate radiation in the selected singular region
+#charmthrpdf 1.5 ! (default 1.5 GeV) pdf charm treshold
+#bottomthrpdf 5.0 ! (default 5.0 GeV) pdf bottom treshold
+
+#iupperisr 1 ! (default 1) choice of ISR upper bounding functional form
+#iupperfsr 2 ! (default 2) choice of FSR upper bounding functional form
+
+#manyseeds 1 ! (default 0) allow for the generation of different statistically independent samples (
+End of powheg.input content
+ Random number generator initialized with: 1 33627373 0
+-->
+<init>
+ 2212 2212 3.50000E+03 3.50000E+03 -1 -1 -1 -1 3 1
+ 9.67831E+02 2.67947E-01 1.00000E+00 10013
+</init>
+<event>
+ 6 10013 1.00000E+00 2.82379E+00 -1.00000E+00 3.04932E-01
+ 1 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 3.607458198E+02 3.607458198E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ -1 -1 0 0 0 511 0.000000000E+00 0.000000000E+00 -6.454449689E+00 6.454449689E+00 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 -2.323316383E-01 2.814217670E+00 3.520085741E+02 3.635691599E+02 9.090942750E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -3.948222587E+00 3.481745721E+00 -4.569218708E+00 6.971365587E+00 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 3.715890949E+00 -6.675280510E-01 3.565777928E+02 3.565977943E+02 1.056583669E-01 0.00000E+00 9.000E+00
+ 21 1 1 2 501 511 2.323316383E-01 -2.814217670E+00 2.282796027E+00 3.631109580E+00 6.615340996E-08 0.00000E+00 9.000E+00
+#pdf 1 -1 0.10217611E+00 0.16507161E-02 0.90909428E+02 0.35811953E+00 0.11922570E+01 1.
+</event>
+<event>
+ 6 10013 1.00000E+00 1.56447E+00 -1.00000E+00 4.21196E-01
+ 3 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 9.713038824E+00 9.713038824E+00 0.000000000E+00 0.00000E+00 9.000E+00
+ -3 -1 0 0 0 511 0.000000000E+00 0.000000000E+00 -2.407279087E+02 2.407279087E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 8.179690066E-01 1.333599350E+00 -2.032460040E+02 2.226280464E+02 9.083975624E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 2.384785824E+01 -2.242676161E+01 -2.461580450E+01 4.095886929E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 -2.302988924E+01 2.376036096E+01 -1.786301995E+02 1.816691771E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ 21 1 1 2 501 511 -8.179690066E-01 -1.333599350E+00 -2.776886588E+01 2.781290119E+01 4.128186458E-07 0.00000E+00 9.000E+00
+#pdf 3 -3 0.27684526E-02 0.60830129E-01 0.90839756E+02 0.91714285E+00 0.12179079E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 1.12054E+02 -1.00000E+00 1.21623E-01
+ 2 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 8.018857281E+02 8.018857281E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ 21 -1 0 0 511 501 0.000000000E+00 0.000000000E+00 -2.822868203E+01 2.822868203E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 9.287337735E+01 -6.269504394E+01 2.224186047E+02 2.676022307E+02 9.790191727E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -8.657041539E+00 -4.183373359E+01 4.445776523E+01 6.165638399E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 1.015304189E+02 -2.086131035E+01 1.779608394E+02 2.059458467E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ 2 1 1 2 511 0 -9.287337735E+01 6.269504394E+01 5.512384414E+02 5.625121795E+02 6.125084199E-06 0.00000E+00 9.000E+00
+#pdf 2 -2 0.46058511E-01 0.42469429E-02 0.97901917E+02 0.63591894E+00 0.79446277E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 1.00369E+01 -1.00000E+00 1.95514E-01
+ 2 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 2.010003947E+02 2.010003947E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ -2 -1 0 0 0 511 0.000000000E+00 0.000000000E+00 -1.529721111E+01 1.529721111E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 3.141951762E+00 -9.532469809E+00 1.839949162E+02 2.061163476E+02 9.235301690E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 1.232415980E+01 -3.569065842E+01 2.632963771E+01 4.603226041E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 -9.182208041E+00 2.615818861E+01 1.576652784E+02 1.600840872E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ 21 1 1 2 501 511 -3.141951762E+00 9.532469809E+00 1.708267440E+00 1.018125823E+01 2.122164147E-07 0.00000E+00 9.000E+00
+#pdf 2 -2 0.55403944E-01 0.31417051E-02 0.92353017E+02 0.62989306E+00 0.90292466E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 1.04202E+01 -1.00000E+00 1.93656E-01
+ 3 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 1.332282132E+02 1.332282132E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ -3 -1 0 0 0 511 0.000000000E+00 0.000000000E+00 -2.326270033E+01 2.326270033E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 -1.021549498E+01 2.055254461E+00 1.126918126E+02 1.457199749E+02 9.179480396E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -2.308841974E+00 1.220137428E+01 -1.380625449E+01 1.856954798E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 -7.906653006E+00 -1.014611982E+01 1.264980671E+02 1.271504269E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ 21 1 1 2 501 511 1.021549498E+01 -2.055254461E+00 -2.726299782E+00 1.077093864E+01 1.933845970E-07 0.00000E+00 9.000E+00
+#pdf 3 -3 0.36680396E-01 0.46881997E-02 0.91794804E+02 0.20111262E+00 0.72101262E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 1.59194E+01 -1.00000E+00 1.74937E-01
+ 2 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 6.771182630E+01 6.771182630E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 21 -1 0 0 511 501 0.000000000E+00 0.000000000E+00 -4.157412921E+01 4.157412921E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 2.738398956E+00 1.568213068E+01 3.627724513E+01 9.041167445E+01 8.126994716E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -2.886571268E+00 -1.222412057E+01 -2.212432259E+01 2.544127023E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 5.624970224E+00 2.790625125E+01 5.840156771E+01 6.497040422E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ 2 1 1 2 511 0 -2.738398956E+00 -1.568213068E+01 -1.013954804E+01 1.887428107E+01 2.469673610E-07 0.00000E+00 9.000E+00
+#pdf 2 -2 0.17760880E-01 0.75892595E-02 0.81269947E+02 0.68683024E+00 0.59404329E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 2.66750E+01 -1.00000E+00 1.56633E-01
+ -3 -1 0 0 0 511 0.000000000E+00 0.000000000E+00 4.528832990E+01 4.528832990E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 3 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 -1.351504985E+02 1.351504985E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 2.433492024E+01 -1.092564348E+01 -1.111498417E+02 1.463107893E+02 9.132908482E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -2.101434561E-01 -2.802960983E+01 -1.129147169E+02 1.163419427E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 2.454506369E+01 1.710396635E+01 1.764875200E+00 2.996884657E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ 21 1 1 2 501 511 -2.433492024E+01 1.092564348E+01 2.128767307E+01 3.412803910E+01 6.425450858E-07 0.00000E+00 9.000E+00
+#pdf -3 3 0.48215418E-02 0.35304998E-01 0.91329085E+02 0.71114296E+00 0.20783819E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 1.47555E+01 -1.00000E+00 1.78011E-01
+ -1 -1 0 0 0 501 0.000000000E+00 0.000000000E+00 7.623749254E+01 7.623749254E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 21 -1 0 0 501 511 0.000000000E+00 0.000000000E+00 -4.136031431E+01 4.136031431E+01 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 -1.208561607E+01 8.465327255E+00 2.106795380E+01 9.738843872E+01 9.393042871E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -3.299092274E+01 -5.810500577E+00 -2.702423006E+01 4.304048181E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 2.090530668E+01 1.427582783E+01 4.809218386E+01 5.434795691E+01 1.056583668E-01 0.00000E+00 9.000E+00
+ -1 1 1 2 0 511 1.208561607E+01 -8.465327255E+00 1.380922443E+01 2.020936814E+01 3.537064561E-07 0.00000E+00 9.000E+00
+#pdf -1 1 0.16717332E-01 0.10770840E-01 0.93930429E+02 0.41769435E+00 0.66171399E+00
+</event>
+<event>
+ 6 10013 1.00000E+00 2.16022E+00 -1.00000E+00 3.47947E-01
+ 1 -1 0 0 511 0 0.000000000E+00 0.000000000E+00 7.250822326E+02 7.250822326E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ -1 -1 0 0 0 501 0.000000000E+00 0.000000000E+00 -4.300923230E+00 4.300923230E+00 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 -2.140464289E+00 -2.915204767E-01 4.678954473E+02 4.764880672E+02 9.005588318E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 2.359202157E+01 3.683460625E+01 1.920555977E+02 1.969739461E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 -2.573248586E+01 -3.712612672E+01 2.758398496E+02 2.795141211E+02 1.056583668E-01 0.00000E+00 9.000E+00
+ 21 1 1 2 511 501 2.140464289E+00 2.915204767E-01 2.528858621E+02 2.528950886E+02 1.066240300E-06 0.00000E+00 9.000E+00
+#pdf 1 -1 0.13487313E+00 0.12271641E-02 0.90055883E+02 0.31305398E+00 0.13403693E+01
+</event>
+<event>
+ 6 10013 1.00000E+00 3.43348E+00 -1.00000E+00 2.80022E-01
+ 2 -1 0 0 501 0 0.000000000E+00 0.000000000E+00 8.958829797E+02 8.958829797E+02 0.000000000E+00 0.00000E+00 9.000E+00
+ 21 -1 0 0 511 501 0.000000000E+00 0.000000000E+00 -7.087092274E+00 7.087092274E+00 0.000000000E+00 0.00000E+00 9.000E+00
+ 23 2 1 2 0 0 1.159842391E+00 3.231643577E+00 8.929604436E+02 8.975726390E+02 9.081023980E+01 0.00000E+00 9.000E+00
+ 13 1 3 3 0 0 -3.567196758E+01 1.692535072E+01 6.710009277E+02 6.721615973E+02 1.056583670E-01 0.00000E+00 9.000E+00
+ -13 1 3 3 0 0 3.683180997E+01 -1.369370714E+01 2.219595159E+02 2.254110417E+02 1.056583669E-01 0.00000E+00 9.000E+00
+ 2 1 1 2 511 0 -1.159842391E+00 -3.231643577E+00 -4.164556187E+00 5.397432984E+00 4.029991083E-07 0.00000E+00 9.000E+00
+#pdf 2 -2 0.25560780E+00 0.65841460E-03 0.90810240E+02 0.41014918E+00 0.17034656E+01
+</event>
+</LesHouchesEvents>
+ #Random number generator exit values: 1 34256472 0
+
+
diff --git a/Generators/Pythia8_i/share/enhanceMPI_example.py b/Generators/Pythia8_i/share/enhanceMPI_example.py
new file mode 100644
index 00000000000..01a6e458edb
--- /dev/null
+++ b/Generators/Pythia8_i/share/enhanceMPI_example.py
@@ -0,0 +1,12 @@
+# example JO for using the EnhanceMPI plugin
+# For testing only - not validated for production yet!!
+
+evgenConfig.description = "Example QCD jets with MPI enhanced with UserHook"
+evgenConfig.keywords = ["QCD", "MPI"]
+include("MC12JobOptions/Pythia8_AU2_CTEQ6L1_Common.py")
+
+topAlg.Pythia8.Commands += ["HardQCD:all=on",
+ "PhaseSpace:pTHatMin=60"]
+topAlg.Pythia8.UserHook = "EnhanceMPI"
+
+
diff --git a/Generators/Pythia8_i/share/jobOptions.pythia8.py b/Generators/Pythia8_i/share/jobOptions.pythia8.py
new file mode 100644
index 00000000000..66c34e32e04
--- /dev/null
+++ b/Generators/Pythia8_i/share/jobOptions.pythia8.py
@@ -0,0 +1,12 @@
+include("GeneratorUtils/StdEvgenSetup.py")
+
+svcMgr.MessageSvc.OutputLevel = INFO
+svcMgr.AtRndmGenSvc.Seeds = ["PYTHIA8 4789899 989240512",
+ "PYTHIA8_INIT 820021 2347532"]
+
+from Pythia8_i.Pythia8_iConf import Pythia8_i
+topAlg += Pythia8_i()
+topAlg.Pythia8_i.CollisionEnergy = 7000
+topAlg.Pythia8_i.Commands += ['HardQCD:all = on']
+
+include("GeneratorUtils/postJO.DumpMC.py")
diff --git a/Generators/Pythia8_i/share/lhef_example.py b/Generators/Pythia8_i/share/lhef_example.py
new file mode 100644
index 00000000000..20a33472ca9
--- /dev/null
+++ b/Generators/Pythia8_i/share/lhef_example.py
@@ -0,0 +1,32 @@
+#import AthenaCommon.AtlasUnixGeneratorJob
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+
+#--------------------------------------------------------------
+# Event related parameters
+#--------------------------------------------------------------
+# Number of events to be processed (default is 10)
+theApp.EvtMax = 10
+#--------------------------------------------------------------
+# Algorithms Private Options
+#--------------------------------------------------------------
+from AthenaServices.AthenaServicesConf import AtRndmGenSvc
+ServiceMgr += AtRndmGenSvc()
+ServiceMgr.AtRndmGenSvc.Seeds = ["PYTHIA8 4789899 989240512",
+ "PYTHIA8_INIT 820021 2347532"]
+
+svcMgr.MessageSvc.OutputLevel = INFO
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+from Pythia8_i.Pythia8_iConf import Pythia8_i
+
+topSequence+=Pythia8_i()
+topSequence.Pythia8_i.CollisionEnergy = 7000
+
+topSequence.Pythia8_i.LHEFile = "Powheg.ZMu.MC11.events"
+
+
+from TruthExamples.TruthExamplesConf import DumpMC
+topSequence += DumpMC()
+
+
diff --git a/Generators/Pythia8_i/share/qcdVetoedShowerExample.py b/Generators/Pythia8_i/share/qcdVetoedShowerExample.py
new file mode 100644
index 00000000000..8f0c1e39163
--- /dev/null
+++ b/Generators/Pythia8_i/share/qcdVetoedShowerExample.py
@@ -0,0 +1,36 @@
+#import AthenaCommon.AtlasUnixGeneratorJob
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+
+#--------------------------------------------------------------
+# Event related parameters
+#--------------------------------------------------------------
+# Number of events to be processed (default is 10)
+theApp.EvtMax = 10
+#--------------------------------------------------------------
+# Algorithms Private Options
+#--------------------------------------------------------------
+from AthenaServices.AthenaServicesConf import AtRndmGenSvc
+ServiceMgr += AtRndmGenSvc()
+ServiceMgr.AtRndmGenSvc.Seeds = ["PYTHIA8 4789899 989240512",
+ "PYTHIA8_INIT 820021 2347532"]
+
+svcMgr.MessageSvc.OutputLevel = INFO
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+from Pythia8_i.Pythia8_iConf import Pythia8_i
+
+topSequence+=Pythia8_i()
+topSequence.Pythia8_i.CollisionEnergy = 7000
+
+#topSequence.Pythia8_i.Commands += ['WeakSingleBoson:ffbar2gmZ = on']
+topSequence.Pythia8_i.LHEFile = "Powheg.ZMu.MC11.events"
+
+topSequence.Pythia8_i.Commands += ['SpaceShower:pTmaxMatch = 2']
+topSequence.Pythia8_i.Commands += ['TimeShower:pTmaxMatch = 2']
+topSequence.Pythia8_i.UserHook = "QCDVetoedShower"
+
+from TruthExamples.TruthExamplesConf import DumpMC
+topSequence += DumpMC()
+
+
diff --git a/Generators/Pythia8_i/share/suppressMPI_example.py b/Generators/Pythia8_i/share/suppressMPI_example.py
new file mode 100644
index 00000000000..5c64b3ab561
--- /dev/null
+++ b/Generators/Pythia8_i/share/suppressMPI_example.py
@@ -0,0 +1,12 @@
+# example JO for using the SuppressMPI plugin
+# For testing only - not validated for production yet!!
+
+evgenConfig.description = "Example QCD jets with MPI suppressed with UserHook"
+evgenConfig.keywords = ["QCD", "MPI"]
+include("MC12JobOptions/Pythia8_AU2_CTEQ6L1_Common.py")
+
+topAlg.Pythia8.Commands += ["HardQCD:all=on",
+ "PhaseSpace:pTHatMin=60"]
+topAlg.Pythia8.UserHook = "SuppressMPI"
+
+
diff --git a/Generators/Pythia8_i/share/userProcess_example.py b/Generators/Pythia8_i/share/userProcess_example.py
new file mode 100644
index 00000000000..b5cca22532f
--- /dev/null
+++ b/Generators/Pythia8_i/share/userProcess_example.py
@@ -0,0 +1,36 @@
+#import AthenaCommon.AtlasUnixGeneratorJob
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+
+#--------------------------------------------------------------
+# Event related parameters
+#--------------------------------------------------------------
+# Number of events to be processed (default is 10)
+theApp.EvtMax = 3
+#--------------------------------------------------------------
+# Algorithms Private Options
+#--------------------------------------------------------------
+from AthenaServices.AthenaServicesConf import AtRndmGenSvc
+ServiceMgr += AtRndmGenSvc()
+ServiceMgr.AtRndmGenSvc.Seeds = ["PYTHIA8 4789899 989240512",
+ "PYTHIA8_INIT 820021 2347532"]
+
+svcMgr.MessageSvc.OutputLevel = INFO
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+from Pythia8_i.Pythia8_iConf import Pythia8_i
+
+
+topSequence+=Pythia8_i()
+topSequence.Pythia8_i.useRndmGenSvc = True
+topSequence.Pythia8_i.UserProcess = "qqbar2emu"
+
+#topSequence.Pythia8_i.Commands += ['Random:setSeed = on']
+#topSequence.Pythia8_i.Commands += ['Random:seed = 123456789']
+
+#topSequence.Pythia8_i.CollisionEnergy = 10000
+
+from TruthExamples.TruthExamplesConf import DumpMC
+topSequence += DumpMC()
+
+
diff --git a/Generators/Pythia8_i/share/vetoedISR_example.py b/Generators/Pythia8_i/share/vetoedISR_example.py
new file mode 100644
index 00000000000..f4eb2a344e1
--- /dev/null
+++ b/Generators/Pythia8_i/share/vetoedISR_example.py
@@ -0,0 +1,35 @@
+#import AthenaCommon.AtlasUnixGeneratorJob
+from AthenaCommon.AppMgr import ServiceMgr as svcMgr
+
+#--------------------------------------------------------------
+# Event related parameters
+#--------------------------------------------------------------
+# Number of events to be processed (default is 10)
+theApp.EvtMax = 10
+#--------------------------------------------------------------
+# Algorithms Private Options
+#--------------------------------------------------------------
+from AthenaServices.AthenaServicesConf import AtRndmGenSvc
+ServiceMgr += AtRndmGenSvc()
+ServiceMgr.AtRndmGenSvc.Seeds = ["PYTHIA8 4789899 989240512",
+ "PYTHIA8_INIT 820021 2347532"]
+
+svcMgr.MessageSvc.OutputLevel = INFO
+
+from AthenaCommon.AlgSequence import AlgSequence
+topSequence = AlgSequence()
+from Pythia8_i.Pythia8_iConf import Pythia8_i
+
+topSequence+=Pythia8_i()
+topSequence.Pythia8_i.CollisionEnergy = 7000
+
+topSequence.Pythia8_i.LHEFile = "Powheg.ZMu.MC11.events"
+
+topSequence.Pythia8_i.Commands += ['SpaceShower:pTmaxMatch = 2',
+ 'TimeShower:pTmaxMatch = 1']
+topSequence.Pythia8_i.UserHook = "ISRVetoedShower"
+
+from TruthExamples.TruthExamplesConf import DumpMC
+topSequence += DumpMC()
+
+
diff --git a/Generators/Pythia8_i/src/Pythia8_i.cxx b/Generators/Pythia8_i/src/Pythia8_i.cxx
new file mode 100644
index 00000000000..6a2bb8cf3c8
--- /dev/null
+++ b/Generators/Pythia8_i/src/Pythia8_i.cxx
@@ -0,0 +1,410 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/Pythia8_i.h"
+#include "Pythia8_i/UserProcessFactory.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include "Pythia8_i/UserResonanceFactory.h"
+
+#include "GeneratorObjects/McEventCollection.h"
+#include <boost/algorithm/string.hpp>
+#include <boost/lexical_cast.hpp>
+
+// calls to fortran routines
+#include "CLHEP/Random/RandFlat.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+// Name of AtRndmGenSvc stream
+std::string Pythia8_i::pythia_stream = "PYTHIA8_INIT";
+
+/**
+ * author: James Monk (jmonk@cern.ch)
+ */
+////////////////////////////////////////////////////////////////////////////////
+Pythia8_i::Pythia8_i(const string &name, ISvcLocator *pSvcLocator)
+: GenModule(name, pSvcLocator),
+m_pythia(xmlpath()),
+m_internal_event_number(0),
+m_version(-1.),
+m_atlasRndmEngine(0),
+m_failureCount(0),
+m_procPtr(0),
+m_userHookPtr(0)
+{
+ declareProperty("Commands", m_commands);
+ declareProperty("CollisionEnergy", m_collisionEnergy = 14000.0);
+ declareProperty("useRndmGenSvc", m_useRndmGenSvc = true);
+ declareProperty("Beam1", m_beam1 = "PROTON");
+ declareProperty("Beam2", m_beam2 = "PROTON");
+ declareProperty("LHEFile", m_lheFile = "");
+ declareProperty("MaxFailures", m_maxFailures = 10);//the max number of consecutive failures
+ declareProperty("UserProcess", m_userProcess="");
+ declareProperty("UserHook", m_userHook="");
+ declareProperty("UserResonances", m_userResonances="");
+ declareProperty("UseLHAPDF", m_useLHAPDF=true);
+ declareProperty("ParticleData", m_particleDataFile="");
+ declareProperty("OutputParticleData",m_outputParticleDataFile="ParticleData.local.xml");
+
+ m_particleIDs["PROTON"] = PROTON;
+ m_particleIDs["ANTIPROTON"] = ANTIPROTON;
+ m_particleIDs["ELECTRON"] = ELECTRON;
+ m_particleIDs["POSITRON"] = POSITRON;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+Pythia8_i::~Pythia8_i() {
+
+ delete m_atlasRndmEngine;
+
+ if(m_procPtr != 0) delete m_procPtr;
+ if(m_userHookPtr != 0) delete m_userHookPtr;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+StatusCode Pythia8_i::genInitialize() {
+
+ ATH_MSG_DEBUG("Pythia8_i from genInitialize()");
+
+ m_version = m_pythia.settings.parm("Pythia:versionNumber");
+
+ std::string pythiaVersion = boost::lexical_cast<std::string>(m_version + 0.00000000001);
+ pythiaVersion.erase(5);
+ std::string libVersion = "8." + std::string(PY8VERSION);
+
+ if(pythiaVersion != libVersion){
+ ATH_MSG_ERROR("Version of Pythia in xmldoc (" + pythiaVersion + ") does not matched linked library version (" + libVersion + ")");
+ return StatusCode::FAILURE;
+ }
+
+ Pythia8_i::pythia_stream = "PYTHIA8_INIT";
+
+ // Default LEP tune is LEP 1 (previously Tune:ee = 3).
+ // Not setting it should be identical, so we don't anymore!
+
+ // LEP1 tune for hadronisation. Default from version 8.125 anyway
+ // m_pythia.readString("Tune:ee = 3");
+
+ // We do explicitly set tune 4C, since it is the starting point for many other tunes
+ // Tune 4C for pp collisions (default from 8.150 anyway)
+ m_pythia.readString("Tune:pp = 5");
+
+ // but also use LHAPDF and CTEQ6L1 unless otherwise instructed in JO
+ if(m_useLHAPDF){
+ m_pythia.readString("PDF:useLHAPDF = on");
+ m_pythia.readString("PDF:LHAPDFset = cteq6ll.LHpdf");
+ }else{
+ m_pythia.readString("PDF:useLHAPDF = off");
+ }
+
+ // Now apply the settings from the JO
+ for(vector<string>::const_iterator cmd = m_commands.begin();
+ cmd != m_commands.end(); ++cmd){
+ bool read = m_pythia.readString(*cmd);
+
+ if(!read){
+ ATH_MSG_ERROR("Pythia could not understand the command '"<< *cmd<<"'");
+ return StatusCode::FAILURE;
+ }
+ }
+
+ PDGID beam1 = m_particleIDs[m_beam1];
+ PDGID beam2 = m_particleIDs[m_beam2];
+
+ ATH_MSG_INFO("Beam1 = "<<beam1);
+ ATH_MSG_INFO("Beam2 = "<<beam2);
+
+ if(beam1 == 0 || beam2 == 0){
+ ATH_MSG_ERROR("Invalid beam particle!");
+ return StatusCode::FAILURE;
+ }
+
+ if(m_useRndmGenSvc){
+
+ ATH_MSG_INFO(" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
+ ATH_MSG_INFO(" THE ATHENA SERVICE AtRndmGenSvc IS USED.");
+ ATH_MSG_INFO(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
+
+ if(m_atlasRndmEngine) delete m_atlasRndmEngine;
+ m_atlasRndmEngine = new customRndm();
+
+ m_atlasRndmEngine->init(atRndmGenSvc(),Pythia8_i::pythia_stream);
+ m_pythia.setRndmEnginePtr(m_atlasRndmEngine);
+
+ // Save the PYTHIA_INIT stream seeds....
+ CLHEP::HepRandomEngine* engine = atRndmGenSvc().GetEngine(Pythia8_i::pythia_stream);
+
+ const long* sip = engine->getSeeds();
+ long int si1 = sip[0];
+ long int si2 = sip[1];
+
+ atRndmGenSvc().CreateStream(si1, si2, Pythia8_i::pythia_stream);
+ Pythia8_i::pythia_stream = "PYTHIA8";
+ m_atlasRndmEngine->m_stream=Pythia8_i::pythia_stream;
+
+ }else{
+ ATH_MSG_INFO(" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
+ ATH_MSG_INFO(" THE STANDARD PYTHIA8 RANDOM NUMBER SERVICE IS USED.");
+ ATH_MSG_INFO(" THE ATHENA SERVICE AtRndmGenSvc IS ***NOT*** USED.");
+ ATH_MSG_INFO(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
+ }
+
+ m_procPtr = 0;
+
+ if(m_userProcess != ""){
+ m_procPtr = Pythia8_UserProcess::UserProcessFactory::create(m_userProcess);
+ }
+
+ bool canInit = true;
+
+ m_userHookPtr = 0;
+
+ if(m_userHook != ""){
+ m_userHookPtr = Pythia8_UserHooks::UserHooksFactory::create(m_userHook);
+ if(!m_pythia.setUserHooksPtr(m_userHookPtr)){
+ ATH_MSG_ERROR("Unable to set requested user hook: " + m_userHook + " !!");
+ ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
+ canInit = false;
+ }
+ }
+
+ if(m_userResonances != ""){
+
+ vector<string> resonanceArgs;
+
+ boost::split(resonanceArgs, m_userResonances, boost::is_any_of(":"));
+ if(resonanceArgs.size() != 2){
+ ATH_MSG_ERROR("Cannot Understand UserResonance job option!");
+ ATH_MSG_ERROR("You should specify it as a 'name:id1,id2,id3...'");
+ ATH_MSG_ERROR("Where name is the name of your UserResonance, and id1,id2,id3 are a comma separated list of PDG IDs to which it is applied");
+ canInit = false;
+ }
+ vector<string> resonanceIds;
+ boost::split(resonanceIds, resonanceArgs.back(), boost::is_any_of(","));
+ if(resonanceIds.size()==0){
+ ATH_MSG_ERROR("You did not specifiy any PDG ids to which your user resonance width should be applied!");
+ ATH_MSG_ERROR("You should specify a list as 'name:id1,id2,id3...'");
+ ATH_MSG_ERROR("Where name is the name of your UserResonance, and id1,id2,id3 are a comma separated list of PDG IDs to which it is applied");
+ canInit=false;
+ }
+
+ for(vector<string>::const_iterator sId = resonanceIds.begin();
+ sId != resonanceIds.end(); ++sId){
+ int idResIn = boost::lexical_cast<int>(*sId);
+ m_userResonancePtrs.push_back(Pythia8_UserResonance::UserResonanceFactory::create(resonanceArgs.front(), idResIn));
+ }
+
+ for(vector<Pythia8::ResonanceWidths*>::const_iterator resonance = m_userResonancePtrs.begin();
+ resonance != m_userResonancePtrs.end(); ++resonance){
+ m_pythia.setResonancePtr(*resonance);
+ }
+
+ }
+
+ if(m_particleDataFile != "") {
+ if(!m_pythia.particleData.reInit(m_particleDataFile, true)){
+ ATH_MSG_ERROR("Unable to read requested particle data table: " + m_particleDataFile + " !!");
+ ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
+ canInit = false;
+ }
+ }
+
+ if(m_lheFile != ""){
+ if(!m_pythia.init(m_lheFile)){
+ ATH_MSG_ERROR("Unable to read requested LHE file: " + m_lheFile + " !!");
+ ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
+ canInit = false;
+ }
+ }else{
+ if(m_procPtr != 0){
+ if(!m_pythia.setSigmaPtr(m_procPtr)){
+ ATH_MSG_ERROR("Unable to set requested user process: " + m_userProcess + " !!");
+ ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
+ canInit = false;
+ }
+ }
+ if(!m_pythia.init(beam1, beam2, m_collisionEnergy)){
+ ATH_MSG_ERROR("Unable to initialise Pythia beams !!");
+ ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
+ canInit = false;
+ }
+ }
+
+ StatusCode returnCode = SUCCESS;
+
+ if(!canInit){
+ returnCode = StatusCode::FAILURE;
+ ATH_MSG_ERROR(" *** Unable to initialise Pythia !! ***");
+ }
+
+ m_pythia.settings.listChanged();
+ m_pythia.particleData.listChanged();
+ m_pythia.particleData.listXML(m_outputParticleDataFile);
+
+ //counter for event failures;
+ m_failureCount = 0;
+
+ m_internal_event_number = 0;
+
+ m_pythiaToHepMC.set_store_pdf(true);
+
+ return returnCode;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+StatusCode Pythia8_i::callGenerator(){
+
+ ATH_MSG_DEBUG(">>> Pythia8_i from callGenerator");
+
+ if(m_useRndmGenSvc){
+ // Save the random number seeds in the event
+ CLHEP::HepRandomEngine* engine = atRndmGenSvc().GetEngine(Pythia8_i::pythia_stream);
+ const long* s = engine->getSeeds();
+ m_seeds.clear();
+ m_seeds.push_back(s[0]);
+ m_seeds.push_back(s[1]);
+ }
+
+ bool status = m_pythia.next();
+
+ StatusCode returnCode = StatusCode::SUCCESS;
+
+ if(!status){
+ ++m_failureCount;
+ if(m_failureCount >= m_maxFailures){
+ ATH_MSG_ERROR("Exceeded the max number of consecutive event failures.");
+ returnCode = StatusCode::FAILURE;
+ }else{
+ ATH_MSG_INFO("Event generation failed - re-trying.");
+ returnCode = this->callGenerator();
+ }
+ }
+
+ m_failureCount = 0;
+
+ // some CKKWL merged events have zero weight (or unfilled event).
+ // start again with such events
+ if(fabs(m_pythia.info.mergingWeight()*m_pythia.info.weight()) < 1.e-18 ||
+ m_pythia.event.size() < 2)
+ returnCode = this->callGenerator();
+
+ ++m_internal_event_number;
+
+ return returnCode;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+StatusCode Pythia8_i::fillEvt(GenEvent *evt){
+
+ ATH_MSG_DEBUG(">>> Pythia8_i from fillEvt");
+
+ evt->set_event_number(m_internal_event_number);
+
+ if(m_pythia.event.size() < 2){
+ ATH_MSG_ERROR("Something wrong with this event - it contains fewer than 2 particles!");
+ ATH_MSG_ERROR("internal event number is "<<m_internal_event_number);
+ return StatusCode::FAILURE;
+ }
+
+ m_pythiaToHepMC.fill_next_event(m_pythia, evt, m_internal_event_number);
+
+ // in debug mode you can check whether the pdf information is stored
+ if(evt->pdf_info()){
+ ATH_MSG_DEBUG("PDFinfo id1:" << evt->pdf_info()->id1());
+ ATH_MSG_DEBUG("PDFinfo id2:" << evt->pdf_info()->id2());
+ ATH_MSG_DEBUG("PDFinfo x1:" << evt->pdf_info()->x1());
+ ATH_MSG_DEBUG("PDFinfo x2:" << evt->pdf_info()->x2());
+ ATH_MSG_DEBUG("PDFinfo scalePDF:" << evt->pdf_info()->scalePDF());
+ ATH_MSG_DEBUG("PDFinfo pdf1:" << evt->pdf_info()->pdf1());
+ ATH_MSG_DEBUG("PDFinfo pdf2:" << evt->pdf_info()->pdf2());
+ }
+ else
+ ATH_MSG_DEBUG("No PDF information available in HepMC::GenEvent!");
+
+ // set the randomseeds
+ if(m_useRndmGenSvc)evt->set_random_states(m_seeds);
+
+ double phaseSpaceWeight = m_pythia.info.weight();
+ double mergingWeight = m_pythia.info.mergingWeight();
+ double eventWeight = phaseSpaceWeight*mergingWeight;
+
+ ATH_MSG_DEBUG("Event weights: phase space weight, merging weight, total weight = "<<phaseSpaceWeight<<", "<<mergingWeight<<", "<<eventWeight);
+
+ // set the event weight
+ evt->weights().clear();
+ evt->weights().push_back(eventWeight);
+
+ // Units correction
+ #ifdef HEPMC_HAS_UNITS
+ //std::cout << (evt->momentum_unit() == HepMC::Units::GEV ? "GeV" : "MeV") << std::endl;
+ /// @todo We shouldn't be having to rescale these events if they're already in MeV :S Where's the screw-up: HepMC or Py8?
+ /// Hopefully this is permanently fixed in version 8.170 onwards
+ if(pythiaVersion() < 8.170 ){
+ GeVToMeV(evt);
+ }
+ #else
+ ATH_MSG_DEBUG("Not scaling HepMC energies for Py8 > 8.153!");
+ /// @todo We *should* be having to rescale these events! Again, is HepMC or Py8 wrong? Or us?
+ //GeVToMeV(evt);
+ #endif
+
+#ifndef HEPMC_HAS_UNITS
+ ATH_MSG_DEBUG();
+#endif
+
+ //HepMC::GenEvent *evtCopy = new HepMC::GenEvent(*evt);
+
+ return SUCCESS;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+StatusCode Pythia8_i::genFinalize(){
+
+ ATH_MSG_INFO(">>> Pythia8_i from genFinalize");
+
+ m_pythia.stat();
+
+ Pythia8::Info info = m_pythia.info;
+ double xs = info.sigmaGen();// in mb
+ xs *= 1000. * 1000.;//convert to nb
+
+ std::cout << "MetaData: cross-section (nb)= " << xs <<std::endl;
+ std::cout << "MetaData: generator= Pythia 8." << PY8VERSION <<std::endl;
+
+ return SUCCESS;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Pythia8_i::pythiaVersion()const{
+ return m_version;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+string Pythia8_i::xmlpath(){
+
+ char *cmtpath = getenv("CMTPATH");
+ char *cmtconfig = getenv("CMTCONFIG");
+
+ string foundpath = "";
+
+ if(cmtpath != 0 && cmtconfig != 0){
+
+ vector<string> cmtpaths;
+ boost::split(cmtpaths, cmtpath, boost::is_any_of(string(":")));
+
+ string installPath = "/InstallArea/" + string(cmtconfig) + "/share/Pythia8/xmldoc";
+
+ for(vector<string>::const_iterator path = cmtpaths.begin();
+ path != cmtpaths.end() && foundpath == ""; ++path){
+ string testPath = *path + installPath;
+ std::ifstream testFile(testPath.c_str());
+ if(testFile.good()) foundpath = testPath;
+ testFile.close();
+ }
+
+ }
+
+ return foundpath;
+}
diff --git a/Generators/Pythia8_i/src/UserHooks/EnhanceMPI.cxx b/Generators/Pythia8_i/src/UserHooks/EnhanceMPI.cxx
new file mode 100644
index 00000000000..6c771551649
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/EnhanceMPI.cxx
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include "boost/lexical_cast.hpp"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class EnhanceMPI;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::EnhanceMPI> EnhanceMPICreator("EnhanceMPI");
+
+namespace Pythia8{
+
+ class EnhanceMPI : public UserHooks{
+
+ public:
+
+ EnhanceMPI(): m_pTCut(10.), m_nMPIVeto(3){
+
+ std::cout<<"**********************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Enhancing MPI emissions with UserHook! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"**********************************************************"<<std::endl;
+
+ m_passedEvent = 0;
+ }
+
+ ~EnhanceMPI(){}
+
+ bool doVetoMPIStep(int nMPI, const Event &event){
+
+ // MPI 1 is the hard process. We do not veto that!
+ if(nMPI < 2){
+ m_passedEvent = false;
+ return false;
+ }
+
+ if(m_passedEvent) return false;
+
+ // start at the end of the event record and work back
+ // This is prior to showering, so there should be at most 2 new MPI emissions
+ // event[0] is documentation, so stop before that.
+ size_t nEmissions=0;
+ for(int ii=event.size()-1; ii > 0 && nEmissions != 2; --ii){
+ if(event[ii].status() != 33) continue;
+ if(event[ii].pT() > m_pTCut){
+ m_passedEvent = true;
+ return false;
+ }
+
+ ++nEmissions;
+ }
+
+ if(nMPI == m_nMPIVeto && !m_passedEvent){
+ // std::cout<<"Vetoing event on too-little MPI"<<std::endl;
+ return true;
+ }
+
+ return false;
+ }
+
+ bool doVetoPartonLevel(const Event&){
+ if(m_passedEvent) return false;
+
+ return true;
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep three times
+ /// First is the hard process
+ /// second is first MPI emission
+ /// *or*
+ /// the second part of a double diffractive event
+ /// *or*
+ /// the second hard process if there is on.
+ /// Therefore check up to 3
+ int numberVetoMPIStep(){return m_nMPIVeto;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return false;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return false;}
+ /// Check the event after the generation of the parton level but before hadronisation
+ bool canVetoPartonLevel(){return true;}
+
+ private:
+
+ double m_pTCut;
+
+ bool m_passedEvent;
+
+ int m_nMPIVeto;
+
+ };
+
+
+
+}
diff --git a/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx b/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx
new file mode 100644
index 00000000000..f83217d1a39
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx
@@ -0,0 +1,117 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserHooksFactory.h"
+#include "Pythia8/PhaseSpace.h"
+
+namespace Pythia8{
+ class GravFlat;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::GravFlat> gravFlatCreator("GravFlat");
+
+namespace Pythia8 {
+
+ class GravFlat : public UserHooks {
+
+ public:
+
+ // Constructor.
+ GravFlat(){}
+
+ // Destructor.
+ ~GravFlat(){}
+
+ // Allow process cross section to be modified...
+ virtual bool canModifySigma() {return true;}
+
+ // ...which gives access to the event at the trial level, before selection.
+ virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr,
+ const PhaseSpace* phaseSpacePtr,
+ bool /* inEvent */) {
+ // All events should be 2 -> 1, but kill them if not.
+ if (sigmaProcessPtr->nFinal() != 1) return 0.;
+
+ // Weight cross section with BW propagator, i.e. to remove it.
+ // (inEvent = false for initialization).
+ // No inEvent criteria, want weight both for cross section
+ // and MC generation.
+ int idRes = sigmaProcessPtr->resonanceA();
+ double mRes = particleDataPtr->m0(idRes);
+ double wRes = particleDataPtr->mWidth(idRes);
+ double m2Res = mRes*mRes;
+ double GamMRat = wRes/mRes;
+ double sHat = phaseSpacePtr->sHat();
+
+ double weightBW = pow2(sHat - m2Res) + pow2(sHat * GamMRat);
+
+ double weightPL = 1;
+ double rH = sqrt(sHat);
+
+ if (rH >= 50 && rH < 2000){
+
+ // Completely Flat M=2 TeV, Above 0.05 TeV, Below 2 TeV
+ // 1e5
+ double p0 = 12.4176;
+ double p1 = -0.364611;
+ double p2 = 0.00293827;
+ double p3 = 4.3106e-07;
+ double p4 = -2.61662e-09;
+ double p5 = 1.33037e-12;
+ double p6 = -2.07848e-16;
+
+ weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6)));
+
+ }
+ else if (rH >= 2000 && rH < 4000){
+
+ // Completely Flat M=2 TeV, Above 2 TeV, Below 4 TeV
+ // 1e5
+ double p0 = 41237.9;
+ double p1 = -95.9745;
+ double p2 = 0.0979667;
+ double p3 = -5.14644e-05;
+ double p4 = 1.45857e-08;
+ double p5 = -2.13169e-12;
+ double p6 = 1.26473e-16;
+
+ weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6)));
+
+ }
+ else if (rH >= 4000 && rH < 5500){
+
+ // Completely Flat M=2 TeV, Above 4 TeV, Below 5.5 TeV
+ // 1e5
+ double p0 = 1.11676e+06;
+ double p1 = -925.647;
+ double p2 = 0.311595;
+ double p3 = -5.38465e-05;
+ double p4 = 4.92286e-09;
+ double p5 = -2.1452e-13;
+ double p6 = 2.97112e-18;
+
+ weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6)));
+
+ }
+ else if (rH >= 5500 && rH <= 6500){
+
+ // Completely Flat M=2 TeV, Above 5.5 TeV, Below 6.5 TeV
+ // 1e5
+ double p0 = 9.70964e+13;
+ double p1 = -4.17142e-03;
+ double p2 = -3.06415e-03;
+
+ weightPL = 1/((p0*exp(p1*rH))+(p2*rH));
+
+ }
+
+ double weight = weightBW * weightPL;
+
+ return weight;
+
+ }
+
+ };
+
+} // end namespace Pythia8
diff --git a/Generators/Pythia8_i/src/UserHooks/ISRVetoedShower.cxx b/Generators/Pythia8_i/src/UserHooks/ISRVetoedShower.cxx
new file mode 100644
index 00000000000..8aa845e09e2
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/ISRVetoedShower.cxx
@@ -0,0 +1,173 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class ISRVetoedShower;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::ISRVetoedShower> ISRVetoedShowerCreator("ISRVetoedShower");
+
+namespace Pythia8{
+
+ /**
+ * This UserHook should be used when showering QCD jets generated with PoWHEG-box.
+ * The following settings must be activated in Pythia:
+ *
+ * SpaceShower:pTMaxMatch = 2
+ * TimeShower:pTMaxMatch = 1
+ *
+ * These set the limit for emission from initial state radiation to the kinematic (beam) energy.
+ * The veto hook then vetos any emission that would be above the appropriate scale.
+ *
+ * This Hook determines a new (lower) veto scale according to the following prescription:
+ * Compute the pT of every PoWHEG leg relative to the beam, and the pT of every outgoing PoWHEG leg relative to all
+ * the other legs in the CMS frame. Use the lowest of these values as the new veto scale
+ *
+ * This hook is ONLY TESTED FOR QCD JET PRODUCTION!!
+ */
+ class ISRVetoedShower : public UserHooks{
+
+ public:
+
+ ISRVetoedShower(): m_powhegScale(0.){
+
+ std::cout<<"*****************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Using new scale definition for QCD production! *"<<std::endl;
+ std::cout<<"* New showering for ISR ONLY - FSR as normal! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"*****************************************************"<<std::endl;
+
+ }
+
+ ~ISRVetoedShower(){}
+
+ /**
+ * doVetoMPIStep is called immediately after the MPI generation
+ * In this case it never actually vetoes the MPI, but since it is called
+ * before the ISR veto check this is a convenient place to find the PoWHEG
+ * scale from the LHEF event
+ */
+ bool doVetoMPIStep(int nMPI, const Event &evt){
+
+ // Only do anything on the first call.
+ if(nMPI > 1) return false;
+
+ m_powhegScale = infoPtr->QRen();;
+
+ // momentum components to boost to CMS frame
+ double pxCMS = 0.;
+ double pyCMS = 0.;
+ double pzCMS = 0.;
+ double eCMS = 0.;
+
+ vector<Particle> powhegLegs;
+
+ // Find the entries corresponding to outgoing legs from PoWHEG
+ // start the loop at 1, since entry 0 represents the event as a whole
+ for(int ii=1; ii != evt.size(); ++ii){
+
+ // status -21 is the incoming hard partons
+ if(evt[ii].status() == -21){
+ pxCMS += evt[ii].px();
+ pyCMS += evt[ii].py();
+ pzCMS += evt[ii].pz();
+ eCMS += evt[ii].e();
+ }
+
+ // here Event::isFinal refers to whether the particle can still decay/radiate, or whether it is an internal leg
+ if(evt[ii].isFinal()){
+ powhegLegs.push_back(Particle(evt[ii]));
+ }
+ }
+
+ // compare the pT of each leg to the powheg scale.
+ // Set the scale to the lowest (or leave the scale unchanged if it is already lower)
+ for(vector<Particle>::const_iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ if(leg->pT() < m_powhegScale )m_powhegScale = leg->pT();
+ }
+
+ double norm = -1./eCMS;
+ pxCMS *= norm;
+ pyCMS *= norm;
+ pzCMS *= norm;
+
+ // boost all outgoing legs to the CMS frame
+ for(vector<Particle>::iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ leg->bst(pxCMS, pyCMS, pzCMS);
+ }
+
+ for(vector<Particle>::const_iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+
+ if(leg->pT() < m_powhegScale )m_powhegScale = leg->pT();
+
+ // calculate the pT relative to each other leg
+ // if any such pT is lower than the current scale, reset the scale to that value
+ for(vector<Particle>::const_iterator otherLeg = powhegLegs.begin();
+ otherLeg != powhegLegs.end(); ++otherLeg){
+ if(otherLeg == leg) continue;
+
+ double pTLeg = Pythia8_UserHooks::pTLeg(*leg, *otherLeg);
+
+ if(pTLeg < m_powhegScale) m_powhegScale = pTLeg;
+ }
+ }
+
+// std::cout<<m_powhegScale<<" "<<infoPtr->QRen()<<std::endl;
+
+
+ return false;
+ }
+
+ /**
+ * This is called after the generation of each new ISR emission
+ * Can use it to test if the last generated emission is above the
+ * veto scale
+ *
+ */
+ bool doVetoISREmission(int, const Event &evt, int iSys){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastISREmission(evt);
+
+ // Veto emission above the veto scale
+ if(evt[emission].pT() > m_powhegScale) return true;
+
+ return false;
+ }
+ /**
+ * This is similar to the ISR veto
+ *
+ */
+ bool doVetoFSREmission(int, const Event&, int, bool){
+
+ return false;
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep once
+ int numberVetoMPIStep(){return 1;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return true;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return false;}
+
+ private:
+
+ double m_powhegScale;
+
+ };
+}
+
diff --git a/Generators/Pythia8_i/src/UserHooks/PTRelVetoedShower.cxx b/Generators/Pythia8_i/src/UserHooks/PTRelVetoedShower.cxx
new file mode 100644
index 00000000000..7fa73342dc4
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/PTRelVetoedShower.cxx
@@ -0,0 +1,219 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include "boost/lexical_cast.hpp"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class PTRelVetoedShower;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::PTRelVetoedShower> pTRelVetoedShowerCreator("PTRelVetoedShower");
+
+namespace Pythia8{
+
+ /**
+ * This user hook is derived from the main31 example distributed with Pythia 8 and
+ * is for use with QCD events generated with PoWHEG. The following settings must be
+ * used in Pythia:
+ *
+ * SpaceShower:pTMaxMatch = 2
+ * TimeShower:pTMaxMatch = 2
+ *
+ * These set the limit for emission from both initial and final state radiation to the kinematic (beam) energy.
+ * The veto hook then vetos any emission that would be above the appropriate scale.
+ *
+ * One key difference between this veto and the one implemented in main31 is the determination of the
+ * boost to the CMS frame in which the FSR veto is tested. In main31, that boost is determined
+ * on each shower emission, with the result that it evolves due to recoil of the incoming partons
+ * against ISR. In this veto, the boost is determined once in each event and does not evolve with recoil
+ */
+
+ class PTRelVetoedShower : public UserHooks{
+
+ public:
+
+ PTRelVetoedShower(): m_nPoWHEGFinal(2), m_powhegScale(0.){
+
+ std::cout<<"*************************************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Using PTRel vetoed shower for PoWHEG QCD production! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"*************************************************************************"<<std::endl;
+
+ m_pxCMS = 0;
+ m_pyCMS = 0;
+ m_pzCMS = 0;
+
+ }
+
+ ~PTRelVetoedShower(){}
+
+ /**
+ * doVetoMPIStep is called immediately after the MPI generation
+ * In this case it never actually vetoes the MPI, but since it is called
+ * before the ISR veto check this is a convenient place to find the PoWHEG
+ * scale from the LHEF event
+ */
+ bool doVetoMPIStep(int nMPI, const Event &evt){
+
+ // Only do anything on the first call.
+ if(nMPI > 1) return false;
+
+ m_powhegScale = infoPtr->QRen();
+
+ size_t nPwgOutgoing = 0;
+
+ for(int ii = evt.size()-1; ii > 0; --ii){
+ if(! (evt[ii].isFinal())) break;
+ ++nPwgOutgoing;
+ }
+
+ if(nPwgOutgoing == m_nPoWHEGFinal) return false;
+
+ if(nPwgOutgoing != m_nPoWHEGFinal + 1){
+ throw std::runtime_error("Wrong number of final state PoWHEG legs: " + boost::lexical_cast<string>(nPwgOutgoing));
+ }
+
+ // momentum components to boost to CMS frame
+ m_pxCMS = 0.;
+ m_pyCMS = 0.;
+ m_pzCMS = 0.;
+ double eCMS = 0.;
+
+ // The outgoing powheg legs
+ vector<Particle> powhegLegs;
+
+ // Find the entries corresponding to outgoing legs from PoWHEG
+ // start the loop at 1, since entry 0 represents the event as a whole
+ for(int ii=1; ii != evt.size(); ++ii){
+
+ // status -21 is the incoming hard partons
+ if(evt[ii].status() == -21){
+ m_pxCMS += evt[ii].px();
+ m_pyCMS += evt[ii].py();
+ m_pzCMS += evt[ii].pz();
+ eCMS += evt[ii].e();
+ }
+
+ // here Event::isFinal refers to whether the particle can still decay/radiate, or whether it is an internal leg
+ if(evt[ii].isFinal()){
+ powhegLegs.push_back(Particle(evt[ii]));
+ }
+ }
+
+ // Start the search for the powheg scale above the collision energy,
+ // which is guaranteed to be above the veto scale.
+
+ m_powhegScale = 1.1*infoPtr->s();
+
+ // compare the pT of each leg to the powheg scale.
+ // Set the scale to the lowest
+ for(vector<Particle>::const_iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ double pTTmp = leg->pT();
+ if(pTTmp < m_powhegScale )m_powhegScale = pTTmp;
+ }
+
+ // normalise the boost vector to the CMS frame...
+ double norm = -1./eCMS;
+ m_pxCMS *= norm;
+ m_pyCMS *= norm;
+ m_pzCMS *= norm;
+
+ // ...and boost all outgoing legs to that frame
+ for(vector<Particle>::iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ leg->bst(m_pxCMS, m_pyCMS, m_pzCMS);
+ }
+
+ for(vector<Particle>::const_iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ // calculate the pT relative to each other leg
+ // if any such pT is lower than the current scale, reset the scale to that value
+ for(vector<Particle>::const_iterator otherLeg = powhegLegs.begin();
+ otherLeg != powhegLegs.end(); ++otherLeg){
+ if(otherLeg == leg) continue;
+
+ double pTLeg = Pythia8_UserHooks::pTLeg(*leg, *otherLeg);
+
+ if(pTLeg < m_powhegScale) m_powhegScale = pTLeg;
+ }
+ }
+
+ // barf if the scale was not found properly!
+ if(m_powhegScale > infoPtr->s()){
+ throw std::runtime_error("Veto scale could not be determined!");
+ }
+
+ return false;
+ }
+
+ /**
+ * This is called after the generation of each new ISR emission
+ * Can use it to test if the last generated emission is above the
+ * veto scale
+ *
+ */
+ bool doVetoISREmission(int, const Event &evt, int iSys){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastISREmission(evt);
+
+ if(evt[emission].pT() > m_powhegScale) return true;
+
+ return false;
+ }
+
+ /**
+ * This is similar to the ISR veto
+ *
+ */
+ bool doVetoFSREmission(int, const Event &evt, int iSys, bool){
+ if(iSys != 0) return false;
+
+ size_t emissionIndex = Pythia8_UserHooks::findLastFSREmission(evt);
+ size_t radiatorIndex = Pythia8_UserHooks::findLastFSRRadiator(evt);
+
+ Particle emission(evt[emissionIndex]);
+ Particle radiator(evt[radiatorIndex]);
+
+ emission.bst(m_pxCMS, m_pyCMS, m_pzCMS);
+ radiator.bst(m_pxCMS, m_pyCMS, m_pzCMS);
+
+ double pTRel = Pythia8_UserHooks::pTLeg(emission, radiator);
+ if(pTRel > m_powhegScale) return true;
+
+ return false;
+ }
+
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep once
+ int numberVetoMPIStep(){return 1;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return true;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return true;}
+
+ private:
+
+ const size_t m_nPoWHEGFinal;
+
+ double m_powhegScale;
+
+ // defines the boost vector to the CMS frame
+ double m_pxCMS;
+ double m_pyCMS;
+ double m_pzCMS;
+
+ };
+}
+
diff --git a/Generators/Pythia8_i/src/UserHooks/PoWHEGVetoedShower.cxx b/Generators/Pythia8_i/src/UserHooks/PoWHEGVetoedShower.cxx
new file mode 100644
index 00000000000..cce7c4451af
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/PoWHEGVetoedShower.cxx
@@ -0,0 +1,180 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class PoWHEGVetoedShower;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::PoWHEGVetoedShower> powhegVetoedShowerCreator("PoWHEGVetoedShower");
+
+namespace Pythia8{
+
+ /**
+ * This UserHook should be used when showering QCD jets generated with PoWHEG-box.
+ * The following settings must be activated in Pythia:
+ *
+ * SpaceShower:pTMaxMatch = 2
+ * TimeShower:pTMaxMatch = 2
+ *
+ * These set the limit for emission from both initial and final state radiation to the kinematic (beam) energy.
+ * The veto hook then vetos any emission that would be above the appropriate scale.
+ *
+ * This Hook determines a new (lower) veto scale according to the following prescription:
+ * Compute the pT of every PoWHEG leg relative to the beam, and the pT of every outgoing PoWHEG leg relative to all
+ * the other legs in the CMS frame. Use the lowest of these values as the new veto scale
+ *
+ * This hook is ONLY SUITABLE FOR QCD JET PRODUCTION!!
+ */
+ class PoWHEGVetoedShower : public UserHooks{
+
+ public:
+
+ PoWHEGVetoedShower(): m_powhegScale(0.){
+
+ std::cout<<"*******************************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Using new PoWHEG author scale definition for QCD production! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"*******************************************************************"<<std::endl;
+
+ }
+
+ ~PoWHEGVetoedShower(){}
+
+ /**
+ * doVetoMPIStep is called immediately after the MPI generation
+ * In this case it never actually vetoes the MPI, but since it is called
+ * before the ISR veto check this is a convenient place to find the PoWHEG
+ * scale from the LHEF event
+ */
+ bool doVetoMPIStep(int nMPI, const Event &evt){
+
+ // Only do anything on the first call.
+ if(nMPI > 1) return false;
+
+ m_powhegScale = infoPtr->QRen();;
+
+ // momentum components to boost to CMS frame
+ double pxCMS = 0.;
+ double pyCMS = 0.;
+ double pzCMS = 0.;
+ double eCMS = 0.;
+
+ vector<Particle> powhegLegs;
+
+ // Find the entries corresponding to outgoing legs from PoWHEG
+ // start the loop at 1, since entry 0 represents the event as a whole
+ for(int ii=1; ii != evt.size(); ++ii){
+
+ // status -21 is the incoming hard partons
+ if(evt[ii].status() == -21){
+ pxCMS += evt[ii].px();
+ pyCMS += evt[ii].py();
+ pzCMS += evt[ii].pz();
+ eCMS += evt[ii].e();
+ }
+
+ // here Event::isFinal refers to whether the particle can still decay/radiate, or whether it is an internal leg
+ if(evt[ii].isFinal()){
+ powhegLegs.push_back(Particle(evt[ii]));
+ }
+ }
+
+ // compare the pT of each leg to the powheg scale.
+ // Set the scale to the lowest (or leave the scale unchanged if it is already lower)
+ for(vector<Particle>::const_iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ if(leg->pT() < m_powhegScale )m_powhegScale = leg->pT();
+ }
+
+ double norm = -1./eCMS;
+ pxCMS *= norm;
+ pyCMS *= norm;
+ pzCMS *= norm;
+
+ // boost all outgoing legs to the CMS frame
+ for(vector<Particle>::iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+ leg->bst(pxCMS, pyCMS, pzCMS);
+ }
+
+ for(vector<Particle>::const_iterator leg=powhegLegs.begin();
+ leg != powhegLegs.end(); ++leg){
+
+ if(leg->pT() < m_powhegScale )m_powhegScale = leg->pT();
+
+ // calculate the pT relative to each other leg
+ // if any such pT is lower than the current scale, reset the scale to that value
+ for(vector<Particle>::const_iterator otherLeg = powhegLegs.begin();
+ otherLeg != powhegLegs.end(); ++otherLeg){
+ if(otherLeg == leg) continue;
+
+ double pTLeg = Pythia8_UserHooks::pTLeg(*leg, *otherLeg);
+
+ if(pTLeg < m_powhegScale) m_powhegScale = pTLeg;
+ }
+ }
+
+// std::cout<<m_powhegScale<<" "<<infoPtr->QRen()<<std::endl;
+
+
+ return false;
+ }
+
+ /**
+ * This is called after the generation of each new ISR emission
+ * Can use it to test if the last generated emission is above the
+ * veto scale
+ *
+ */
+ bool doVetoISREmission(int, const Event &evt, int iSys){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastISREmission(evt);
+
+ // Veto emission above the veto scale
+ if(evt[emission].pT() > m_powhegScale) return true;
+
+ return false;
+ }
+ /**
+ * This is similar to the ISR veto
+ *
+ */
+ bool doVetoFSREmission(int, const Event &evt, int iSys, bool){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastFSREmission(evt);
+
+ // Veto if above the POWHEG scale
+ if (evt[emission].pT() > m_powhegScale) return true;
+
+ return false;
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep once
+ int numberVetoMPIStep(){return 1;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return true;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return true;}
+
+ private:
+
+ double m_powhegScale;
+
+ };
+}
+
diff --git a/Generators/Pythia8_i/src/UserHooks/QCDVetoedShower.cxx b/Generators/Pythia8_i/src/UserHooks/QCDVetoedShower.cxx
new file mode 100644
index 00000000000..fa2c01aac40
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/QCDVetoedShower.cxx
@@ -0,0 +1,138 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class QCDVetoedShower;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::QCDVetoedShower> qcdVetoedShowerCreator("QCDVetoedShower");
+
+namespace Pythia8{
+
+ /**
+ * This UserHook should be used when showering QCD jets generated with PoWHEG-box.
+ * The following settings must be activated in Pythia:
+ *
+ * SpaceShower:pTMaxMatch = 2
+ * TimeShower:pTMaxMatch = 2
+ *
+ * These set the limit for emission from both initial and final state radiation to the kinematic (beam) energy.
+ * The veto hook then vetos any emission that would be above the appropriate scale.
+ *
+ * In QCD, emissions could come from any PoWHEG leg, therefore each candidate Pythia emission has its pT determined
+ * relative to each PoWHEG leg, including the beams. If any of these pTs is larger than the PoWHEG veto scale then
+ * that proposed emission is rejected.
+ *
+ * this differs from colourless resonance production, which requires a different scheme and therefore this hook...
+ *
+ * ...SHOULD NOT BE USED FOR W/Z OR OTHER COLOURLESS RESONANCE PRODUCTION!!
+ */
+ class QCDVetoedShower : public UserHooks{
+
+ public:
+
+ QCDVetoedShower(): m_powhegScale(0.){
+
+ std::cout<<"*******************************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Using vetoed shower for PoWHEG QCD production! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"*******************************************************************"<<std::endl;
+
+ }
+
+ ~QCDVetoedShower(){}
+
+ /**
+ * doVetoMPIStep is called immediately after the MPI generation
+ * In this case it never actually vetoes the MPI, but since it is called
+ * before the ISR veto check this is a convenient place to find the PoWHEG
+ * scale from the LHEF event
+ */
+ bool doVetoMPIStep(int nMPI, const Event &evt){
+
+ // Only do anything on the first call.
+ if(nMPI > 1) return false;
+
+ m_powhegScale = infoPtr->QRen();
+
+ m_powhegLegs.clear();
+ // start the loop at 1, since entry 0 represents the event as a whole
+ for(int ii=1; ii != evt.size(); ++ii){
+ // here Event::isFinal refers to whether the particle can still decay/radiate, or whether it is an internal leg
+ if(evt[ii].isFinal()){
+ m_powhegLegs.push_back(ii);
+ }
+ }
+
+ return false;
+ }
+
+ /**
+ * This is called after the generation of each new ISR emission
+ * Can use it to test if the last generated emission is above the
+ * veto scale
+ *
+ */
+ bool doVetoISREmission(int, const Event &evt, int iSys){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastISREmission(evt);
+
+ return checkEmission(emission, evt);
+
+ }
+
+ /**
+ * This is similar to the ISR veto
+ *
+ */
+ bool doVetoFSREmission(int, const Event &evt, int iSys, bool){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastFSREmission(evt);
+
+ return checkEmission(emission, evt);
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep once
+ int numberVetoMPIStep(){return 1;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return true;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return true;}
+
+ private:
+
+ bool checkEmission(size_t emission, const Event &evt){
+ // Veto if above the POWHEG scale
+ if (evt[emission].pT() > m_powhegScale) return true;
+
+ // determine pT relative to each outgoing PoWHEG leg and veto if above veto scale
+
+ for(vector<size_t>::const_iterator legIndex = m_powhegLegs.begin();
+ legIndex != m_powhegLegs.end(); ++legIndex){
+
+ if(Pythia8_UserHooks::pTLeg(emission, *legIndex, evt) > m_powhegScale)return true;
+ }
+ return false;
+ }
+
+ double m_powhegScale;
+ // Indices of the PoWHEG legs in the event record
+ vector<size_t> m_powhegLegs;
+ };
+}
+
diff --git a/Generators/Pythia8_i/src/UserHooks/SuppressMPI.cxx b/Generators/Pythia8_i/src/UserHooks/SuppressMPI.cxx
new file mode 100644
index 00000000000..e43691982b3
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/SuppressMPI.cxx
@@ -0,0 +1,85 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include "boost/lexical_cast.hpp"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class SuppressMPI;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::SuppressMPI> SuppressMPICreator("SuppressMPI");
+
+namespace Pythia8{
+
+ class SuppressMPI : public UserHooks{
+
+ public:
+
+ SuppressMPI(): m_pTCut(10.){
+
+ std::cout<<"************************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Suppressing MPI emissions with UserHook! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"************************************************************"<<std::endl;
+
+ m_nMPIVeto = 0;
+ }
+
+ ~SuppressMPI(){}
+
+ bool doVetoMPIStep(int nMPI, const Event &event){
+ // MPI 1 is the hard process. We do not veto that!
+ if(nMPI < 2){
+ return false;
+ }
+
+ // start at the end of the event record and work back
+ // This is prior to showering, so there should be at most 2 new MPI emissions
+ // event[0] is documentation, so stop before that.
+ size_t nEmissions=0;
+ for(int ii=event.size()-1; ii > 0 && nEmissions != 2; --ii){
+ if(event[ii].status() != 33) continue;
+ if(event[ii].pT() > m_pTCut){
+ return true;
+ }
+
+ ++nEmissions;
+ }
+
+ return false;
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep three times
+ /// First is the hard process
+ /// second is first MPI emission
+ /// *or*
+ /// the second part of a double diffractive event
+ /// *or*
+ /// the second hard process if there is on.
+ /// Therefore check up to 3
+ int numberVetoMPIStep(){return m_nMPIVeto;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return false;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return false;}
+
+ private:
+
+ double m_pTCut;
+
+ int m_nMPIVeto;
+
+
+ };
+
+
+
+}
diff --git a/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h b/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h
new file mode 100644
index 00000000000..0d73d312521
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/UserHooksUtils.h
@@ -0,0 +1,116 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PYTHIA8_USERHOOKS_USERHOOKSUTILS_H
+#define PYTHIA8_USERHOOKS_USERHOOKSUTILS_H
+#include "Pythia8/Event.h"
+#include <stdexcept>
+
+/**
+ * Some common functions for determining pTs and navigating event records for the PoWHEG + Pythia user hooks
+ */
+
+namespace Pythia8_UserHooks{
+
+ /**
+ * \return the dot product of \param leg and \param comparison
+ */
+ inline double pTProj(const Pythia8::Particle &leg, const Pythia8::Particle &comparison){
+ return leg.px()*comparison.px() + leg.py()*comparison.py() + leg.pz()*comparison.pz();
+ }
+
+ /**
+ * \return the dot product of the legs in the \param evt with indices
+ * \param legIndex and \param comparedIndex
+ */
+ inline double pTProj(size_t legIndex, size_t comparedIndex, const Pythia8::Event &evt){
+ return pTProj(evt[legIndex], evt[comparedIndex]);
+ }
+
+ /**
+ * \return the pT squared of \param leg relative to \param comparison
+ */
+ inline double pT2Leg(const Pythia8::Particle &leg, const Pythia8::Particle &comparison){
+ double proj = pTProj(leg, comparison);
+ return leg.pAbs2() - (proj*proj) / comparison.pAbs2();
+ }
+
+ /**
+ * \return the pT squared of the leg with index \param legIndex compared to the leg with index
+ * \param comparedIndex in Event \param evt
+ */
+ inline double pT2Leg(size_t legIndex, size_t comparedIndex, const Pythia8::Event &evt){
+ return pT2Leg(evt[legIndex], evt[comparedIndex]);
+ }
+
+ /**
+ * \return the pT of \param leg relative to \param comparison
+ */
+ inline double pTLeg(const Pythia8::Particle &leg, const Pythia8::Particle &comparison){
+ return sqrt(pT2Leg(leg, comparison));
+ }
+
+ /**
+ * \return the pT of the leg with index \param legIndex compared to the leg with index
+ * \param comparedIndex in Event \param evt
+ */
+ inline double pTLeg(size_t legIndex, size_t comparedIndex, const Pythia8::Event &evt){
+ return sqrt(pT2Leg(legIndex, comparedIndex, evt));
+ }
+
+ /**
+ * Return the index of the most recent emission in a \param evt with a given \param status
+ */
+ inline size_t findLastEmission(const Pythia8::Event &evt, int status){
+ size_t emission = evt.size() - 1;
+
+ while(emission != 0){
+ if (evt[emission].isFinal() && evt[emission].status() == status) return emission;
+
+ --emission;
+ }
+
+ return 0;
+ }
+
+ /**
+ * \return the index of the most recent ISR emission in an \param evt
+ */
+ inline size_t findLastISREmission(const Pythia8::Event &evt){
+ size_t emission = findLastEmission(evt, 43);
+ if(emission == 0) throw std::runtime_error("findLastISREmission:: Could not find ISR emission");
+ return emission;
+ }
+
+ /**
+ * \return the index of the most recent FSR emission in an \param evt
+ */
+ inline size_t findLastFSREmission(const Pythia8::Event &evt){
+ size_t emission = findLastEmission(evt, 51);
+ if(emission == 0) throw std::runtime_error("findLastFSREmission:: Could not find FSR emission");
+ return emission;
+ }
+
+ /**
+ * \return the index of the most recent particle to radiate an ISR emission in an \param evt
+ */
+
+ inline size_t findLastISRRadiator(const Pythia8::Event &evt){
+ size_t radiator = findLastEmission(evt, -41);
+ if(radiator == 0) throw std::runtime_error("findLastISRRadiator:: Could not find ISR radiator");
+ return radiator;
+ }
+
+ /**
+ * \return the index of the most recent particle to radiate a FSR emission in an \param evt
+ */
+
+ inline size_t findLastFSRRadiator(const Pythia8::Event &evt){
+ size_t emitted = findLastFSREmission(evt);
+ return evt[emitted].mother1();
+ }
+
+}
+
+#endif
diff --git a/Generators/Pythia8_i/src/UserHooks/VetoedShower.cxx b/Generators/Pythia8_i/src/UserHooks/VetoedShower.cxx
new file mode 100644
index 00000000000..9282b58dbc5
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/VetoedShower.cxx
@@ -0,0 +1,151 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserHooksFactory.h"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class VetoedShower;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::VetoedShower> vetoedShowerCreator("VetoedShower");
+
+namespace Pythia8{
+
+
+ class VetoedShower : public UserHooks{
+
+ public:
+
+ VetoedShower(): m_pTVeto(0.), m_doneISR(false){
+ std::cout<<std::endl;
+ std::cout<<"************************************************************************************************"<<std::endl;
+ std::cout<<"Pythia8::VetoedShower: *** WARNING: THIS USER HOOK IS NOT VALIDATED FOR QCD JET PRODUCTION!! ***"<<std::endl;
+ std::cout<<"************************************************************************************************"<<std::endl;
+ std::cout<<std::endl;
+ }
+
+ ~VetoedShower(){}
+
+ /**
+ * doVetoMPIStep is called immediately after the MPI generation
+ * In this case it nver actually vetoes the MPI, but since it is called
+ * before the ISR veto check this is a convenient place to find the veto
+ * scale from the LHEF event
+ */
+ bool doVetoMPIStep(int, const Event&){
+
+ m_pTVeto = infoPtr->QFac();
+ m_doneISR = false;
+
+// std::cout<<" ******** veto scale = "<<m_pTVeto<<" *********"<<std::endl;
+
+ return false;
+ }
+
+ /**
+ * This is called after the generation of each new ISR emission
+ * Can use it to test if the last generated emission is above the
+ * veto scale
+ *
+ * Since the ISR is ordered in pT (with the same definition relative to the
+ * beam as Powheg), once we find an emission below the veto scale we know that
+ * all subsequent emissions will also be below the scale. Therefore m_doneISR
+ * gets set to true at that point and nothing more needs doing
+ *
+ */
+ bool doVetoISREmission(int, const Event &e, int){
+ // Note ignoring new iSys param for now (addition required to fix compilation against >= 8.157
+
+ if(m_doneISR) return false;
+
+ int i = e.size() - 1;
+
+ while(i != -1){
+
+ if(e[i].isFinal() && e[i].status() == 43) break;
+
+ --i;
+
+ if(i == 0){
+ e.list();
+ throw std::runtime_error("Pythia8::VetoedShower::doVetoISREmission: Could not find ISR emission");
+ }
+ }
+
+ //check the leg is connected to the hard system
+ int iMother = e[i].mother1();
+ if(e[iMother].status() != -41)
+ throw std::runtime_error("Pythia8::VetoedShower::doVetoISREmission: Unexpected status code in ISR");
+
+ iMother = e[iMother].daughter2();
+ if(iMother != partonSystemsPtr->getInA(0) &&
+ iMother != partonSystemsPtr->getInB(0)) return false;
+
+ // Veto emission above the veto scale
+ if(e[i].pT() > m_pTVeto) return true;
+
+ m_doneISR = true;
+ return false;
+ }
+ /**
+ * This is similar to the ISR veto, with the exception that since the ordering is in
+ * a slightly different pT definition to Powheg the veto must always be checked, rather
+ * than stopping once an emission is allowed
+ *
+ */
+ bool doVetoFSREmission(int, const Event &e, int, bool){
+ // Note ignoring new iSys and isResonance params for now (addition required to fix compilation against >= 8.157
+
+ int i = e.size() - 1;
+
+ while(i != -1){
+ if (e[i].isFinal() && e[i].status() == 51) break;
+
+ --i;
+
+ if(i == 0){
+ e.list();
+ throw std::runtime_error("Pythia8::VetoedShower::doVetoFSREmission: Could not find FSR emission");
+ }
+ }
+
+ // Make sure radiation from the hard system
+ int iMother = e[i].mother1();
+ int sysSize = partonSystemsPtr->sizeOut(0);
+ bool hardSys = false;
+
+ for (int j = 0; j < sysSize; ++j){
+ int iOut = partonSystemsPtr->getOut(0, j);
+ if (iOut == iMother) {
+ hardSys = true;
+ break;
+ }
+ }
+
+ if (!hardSys) return false;
+
+ // Veto if above the POWHEG scale
+ if (e[i].pT() > m_pTVeto) return true;
+
+ return false;
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep once
+ int numberVetoMPIStep(){return 1;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return true;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return true;}
+
+ private:
+
+ double m_pTVeto;
+ bool m_doneISR;
+ };
+}
+
diff --git a/Generators/Pythia8_i/src/UserHooks/WZVetoedShower.cxx b/Generators/Pythia8_i/src/UserHooks/WZVetoedShower.cxx
new file mode 100644
index 00000000000..efe72b1d32c
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/WZVetoedShower.cxx
@@ -0,0 +1,128 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class WZVetoedShower;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::WZVetoedShower> wzVetoedShowerCreator("WZVetoedShower");
+
+namespace Pythia8{
+
+ /**
+ * This UserHook should be used when showering W, Z or other colourless resonances generated with PoWHEG.
+ * The following setting must be activated in Pythia:
+ *
+ * SpaceShower:pTMaxMatch = 2
+ * TimeShower:pTMaxMatch = 2
+ *
+ * These set the limit for emission from both initial and final state radiation to the kinematic (beam) energy.
+ * The veto hook then vetos any emission that would be above the appropriate scale.
+ *
+ * In the cases of colourless resonance production there should be no FSR emission from PoWHEG, therefore the veto
+ * determines the pT of every Pythia emission relative to the beam, and vetos if it is above the PoWHEG scale.
+ *
+ * This is ok because
+ *
+ * 1) Although you might think the pT should also be determined relative to any PoWHEG ISR emission, there is no way
+ * PoWHEG could have generated an emission off its own ISR.
+ * 2) There can be FSR off what was a PoWHEG ISR leg, but again, PoWHEG could never have generated that emission, so you
+ * only need to compare to the incoming (beam) legs.
+ *
+ * This is in contrast to QCD production, which requires a different veto scheme. As such, this hook...
+ *
+ * ...SHOULD NOT BE USED FOR QCD JET PRODUCTION!!
+ *
+ */
+
+ class WZVetoedShower : public UserHooks{
+
+ public:
+
+ WZVetoedShower(): m_powhegScale(0.){
+
+ std::cout<<"*******************************************************************"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"* Using vetoed shower for PoWHEG colourless resonance production! *"<<std::endl;
+ std::cout<<"* *"<<std::endl;
+ std::cout<<"*******************************************************************"<<std::endl;
+
+ }
+
+ ~WZVetoedShower(){}
+
+ /**
+ * doVetoMPIStep is called immediately after the MPI generation
+ * In this case it never actually vetoes the MPI, but since it is called
+ * before the ISR veto check this is a convenient place to find the PoWHEG
+ * scale from the LHEF event
+ */
+ bool doVetoMPIStep(int nMPI, const Event &){
+
+ // Only do anything on the first call.
+ if(nMPI > 1) return false;
+
+ m_powhegScale = infoPtr->QRen();
+// m_powhegScale = infoPtr->QFac();
+
+// std::cout<<" ******** veto scale = "<<m_powhegScale<<", "<<infoPtr->QRen() <<" *********"<<std::endl;
+
+ return false;
+ }
+
+ /**
+ * This is called after the generation of each new ISR emission
+ * Can use it to test if the last generated emission is above the
+ * veto scale
+ *
+ */
+ bool doVetoISREmission(int, const Event &evt, int iSys){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastISREmission(evt);
+
+ // Veto emission above the veto scale
+ if(evt[emission].pT() > m_powhegScale) return true;
+
+ return false;
+ }
+ /**
+ * This is similar to the ISR veto
+ *
+ */
+ bool doVetoFSREmission(int, const Event &evt, int iSys, bool){
+
+ // only veto emissions from the hard system
+ if(iSys != 0) return false;
+
+ size_t emission = Pythia8_UserHooks::findLastFSREmission(evt);
+
+ // Veto if above the POWHEG scale
+ if (evt[emission].pT() > m_powhegScale) return true;
+
+ return false;
+ }
+
+ /// Switch on calling of doVetoMPIStep
+ bool canVetoMPIStep(){return true;}
+ /// Call doVetoMIStep once
+ int numberVetoMPIStep(){return 1;}
+ /// Switch on veto of ISR
+ bool canVetoISREmission(){return true;}
+ /// Switch off veto of FSR
+ bool canVetoFSREmission(){return true;}
+
+ private:
+
+ double m_powhegScale;
+ };
+}
+
diff --git a/Generators/Pythia8_i/src/UserHooks/WprimeFlat.cxx b/Generators/Pythia8_i/src/UserHooks/WprimeFlat.cxx
new file mode 100644
index 00000000000..a1b4b934756
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/WprimeFlat.cxx
@@ -0,0 +1,60 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserHooksFactory.h"
+#include "Pythia8/PhaseSpace.h"
+
+namespace Pythia8{
+ class WprimeFlat;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::WprimeFlat> WprimeFlatCreator("WprimeFlat");
+
+namespace Pythia8 {
+
+ class WprimeFlat : public UserHooks {
+
+ public:
+
+ // Constructor.
+ WprimeFlat(){}
+
+ // Destructor.
+ ~WprimeFlat(){}
+
+ // Allow process cross section to be modified...
+ virtual bool canModifySigma() {return true;}
+
+ // ...which gives access to the event at the trial level, before selection.
+ virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr,
+ const PhaseSpace* phaseSpacePtr,
+ bool /* inEvent */) {
+ // All events should be 2 -> 1, but kill them if not.
+ if (sigmaProcessPtr->nFinal() != 1) return 0.;
+
+ // Weight cross section with BW propagator, i.e. to remove it.
+ // (inEvent = false for initialization).
+ // No inEvent criteria, want weight both for cross section
+ // and MC generation.
+ int idRes = sigmaProcessPtr->resonanceA();
+ double mRes = particleDataPtr->m0(idRes);
+ double wRes = particleDataPtr->mWidth(idRes);
+ double m2Res = mRes*mRes;
+ double GamMRat = wRes/mRes;
+ double sHat = phaseSpacePtr->sHat();
+ double weightBW = pow2(sHat - m2Res) + pow2(sHat * GamMRat);
+ double m = sqrt(sHat)/8000.0;
+ if(m < 0.0375)
+ weightBW *= 121.88e-12*exp(13.0*m);
+ else
+ weightBW *= 1.0e-12*exp(18.5*m-1.4*log(m));
+
+ return weightBW;
+ }
+
+ };
+
+} // end namespace Pythia8
+
+
diff --git a/Generators/Pythia8_i/src/UserHooks/WprimeWZFlat.cxx b/Generators/Pythia8_i/src/UserHooks/WprimeWZFlat.cxx
new file mode 100644
index 00000000000..c27743d2839
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/WprimeWZFlat.cxx
@@ -0,0 +1,68 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserHooksFactory.h"
+#include "Pythia8/PhaseSpace.h"
+
+namespace Pythia8{
+ class WprimeWZFlat;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::WprimeWZFlat> WprimeWZFlatCreator("WprimeWZFlat");
+
+namespace Pythia8 {
+
+ class WprimeWZFlat : public UserHooks {
+
+ public:
+
+ // Constructor.
+ WprimeWZFlat(){}
+
+ // Destructor.
+ ~WprimeWZFlat(){}
+
+ // Allow process cross section to be modified...
+ virtual bool canModifySigma() {return true;}
+
+ // ...which gives access to the event at the trial level, before selection.
+ virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr,
+ const PhaseSpace* phaseSpacePtr,
+ bool /* inEvent */) {
+ // All events should be 2 -> 1, but kill them if not.
+ if (sigmaProcessPtr->nFinal() != 1) return 0.;
+
+ // Weight cross section with BW propagator, i.e. to remove it.
+ // (inEvent = false for initialization).
+ // No inEvent criteria, want weight both for cross section
+ // and MC generation.
+ int idRes = sigmaProcessPtr->resonanceA();
+ double mRes = particleDataPtr->m0(idRes);
+ double wRes = particleDataPtr->mWidth(idRes);
+ double m2Res = mRes*mRes;
+ double gamMRat = wRes/mRes;
+ double sHat = phaseSpacePtr->sHat();
+ double weightBW = m2Res*m2Res + sHat*sHat*(1 + gamMRat*gamMRat) - 2.*sHat*m2Res;
+
+ return weightBW * breitWignerDenom(sqrt(sHat)/8000.0);
+ }
+
+ private:
+
+ double breitWignerDenom(double mFrac){
+
+ if(mFrac < 0.025) return breitWignerDenom(0.025);
+ if(mFrac > 0.6) return breitWignerDenom(0.6);
+
+ if(mFrac < 0.0425) return 1e-12/(-1.293+1.098e+2*mFrac-2.800e+3*mFrac*mFrac+2.345e+4*mFrac*mFrac*mFrac);
+ if(mFrac < 0.073) return 1.248e-12*(exp(1.158+18.34*mFrac));
+
+ return 5.733e-10*pow(mFrac,-3.798-0.6555*log(mFrac))/pow(1.427-mFrac,30.017);
+ }
+
+ };
+
+} // end namespace Pythia8
+
+
diff --git a/Generators/Pythia8_i/src/UserHooks/main31.cxx b/Generators/Pythia8_i/src/UserHooks/main31.cxx
new file mode 100644
index 00000000000..3898f605ae0
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooks/main31.cxx
@@ -0,0 +1,505 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "UserHooksUtils.h"
+#include "Pythia8_i/UserHooksFactory.h"
+#include "boost/lexical_cast.hpp"
+#include <stdexcept>
+#include <iostream>
+
+namespace Pythia8{
+ class main31;
+}
+
+Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::main31> main31Creator("Main31");
+
+namespace Pythia8{
+
+ // Use userhooks to veto PYTHIA emissions above the POWHEG scale.
+
+ class main31 : public UserHooks {
+
+ public:
+
+ // Constructor and destructor.
+ main31() : nFinal(2),
+ vetoMode(1), vetoCount(3),
+ pThardMode(2), pTemtMode(0),
+ emittedMode(0), pTdefMode(1),
+ MPIvetoMode(0),
+ pThard(0), pTMPI(0),
+ accepted(0),
+ nAcceptSeq(0),
+ nISRveto(0), nFSRveto(0) {};
+ ~main31() {}
+
+ //--------------------------------------------------------------------------
+
+ // Routines to calculate the pT (according to pTdefMode) in a splitting:
+ // ISR: i (radiator after) -> j (emitted after) k (radiator before)
+ // FSR: i (radiator before) -> j (emitted after) k (radiator after)
+ // For the Pythia pT definition, a recoiler (after) must be specified.
+
+ // Compute the Pythia pT separation. Based on pTLund function in History.cc
+ double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch,
+ int RecAfterBranch, bool FSR) {
+
+ // Convenient shorthands for later
+ Vec4 radVec = e[RadAfterBranch].p();
+ Vec4 emtVec = e[EmtAfterBranch].p();
+ Vec4 recVec = e[RecAfterBranch].p();
+ int radID = e[RadAfterBranch].id();
+
+ // Calculate virtuality of splitting
+ double sign = (FSR) ? 1. : -1.;
+ Vec4 Q(radVec + sign * emtVec);
+ double Qsq = sign * Q.m2Calc();
+
+ // Mass term of radiator
+ double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
+ pow2(particleDataPtr->m0(radID)) : 0.;
+
+ // z values for FSR and ISR
+ double z, pTnow;
+ if (FSR) {
+ // Construct 2 -> 3 variables
+ Vec4 sum = radVec + recVec + emtVec;
+ double m2Dip = sum.m2Calc();
+ double x1 = 2. * (sum * radVec) / m2Dip;
+ double x3 = 2. * (sum * emtVec) / m2Dip;
+ z = x1 / (x1 + x3);
+ pTnow = z * (1. - z);
+
+ } else {
+ // Construct dipoles before/after splitting
+ Vec4 qBR(radVec - emtVec + recVec);
+ Vec4 qAR(radVec + recVec);
+ z = qBR.m2Calc() / qAR.m2Calc();
+ pTnow = (1. - z);
+ }
+
+ // Virtuality with correct sign
+ pTnow *= (Qsq - sign * m2Rad);
+
+ // Can get negative pT for massive splittings
+ if (pTnow < 0.) {
+ cout << "Warning: pTpythia was negative" << endl;
+ return -1.;
+ }
+
+#ifdef DBGOUTPUT
+ cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
+ << EmtAfterBranch << ", rec = " << RecAfterBranch
+ << ", pTnow = " << sqrt(pTnow) << endl;
+#endif
+
+ // Return pT
+ return sqrt(pTnow);
+ }
+
+ // Compute the POWHEG pT separation between i and j
+ double pTpowheg(const Event &e, int i, int j, bool FSR) {
+
+ // pT value for FSR and ISR
+ double pTnow = 0.;
+ if (FSR) {
+ // POWHEG d_ij (in CM frame). Note that the incoming beams have not
+ // been updated in the parton systems pointer yet (i.e. prior to any
+ // potential recoil).
+ int iInA = partonSystemsPtr->getInA(0);
+ int iInB = partonSystemsPtr->getInB(0);
+ double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
+ ( e[iInA].e() + e[iInB].e() );
+
+ Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
+ iVecBst.bst(0., 0., betaZ);
+ jVecBst.bst(0., 0., betaZ);
+ pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
+ iVecBst.e() * jVecBst.e() /
+ pow2(iVecBst.e() + jVecBst.e()) );
+
+ } else {
+ // POWHEG pT_ISR is just kinematic pT
+ pTnow = e[j].pT();
+ }
+
+ // Check result
+ if (pTnow < 0.) {
+ cout << "Warning: pTpowheg was negative" << endl;
+ return -1.;
+ }
+
+#ifdef DBGOUTPUT
+ cout << "pTpowheg: i = " << i << ", j = " << j
+ << ", pTnow = " << pTnow << endl;
+#endif
+
+ return pTnow;
+ }
+
+ // Calculate pT for a splitting based on pTdefMode.
+ // If j is -1, all final-state partons are tried.
+ // If i, k, r and xSR are -1, then all incoming and outgoing
+ // partons are tried.
+ // xSR set to 0 means ISR, while xSR set to 1 means FSR
+ double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
+
+ // Loop over ISR and FSR if necessary
+ double pTemt = -1., pTnow;
+ int xSR1 = (xSRin == -1) ? 0 : xSRin;
+ int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
+ for (int xSR = xSR1; xSR < xSR2; xSR++) {
+ // FSR flag
+ bool FSR = (xSR == 0) ? false : true;
+
+ // If all necessary arguments have been given, then directly calculate.
+ // POWHEG ISR and FSR, need i and j.
+ if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
+ pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
+
+ // Pythia ISR, need i, j and r.
+ } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
+ pTemt = pTpythia(e, i, j, r, FSR);
+
+ // Pythia FSR, need k, j and r.
+ } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
+ pTemt = pTpythia(e, k, j, r, FSR);
+
+ // Otherwise need to try all possible combintations.
+ } else {
+ // Start by finding incoming legs to the hard system after
+ // branching (radiator after branching, i for ISR).
+ // Use partonSystemsPtr to find incoming just prior to the
+ // branching and track mothers.
+ int iInA = partonSystemsPtr->getInA(0);
+ int iInB = partonSystemsPtr->getInB(0);
+ while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
+ while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
+
+ // If we do not have j, then try all final-state partons
+ int jNow = (j > 0) ? j : 0;
+ int jMax = (j > 0) ? j + 1 : e.size();
+ for (; jNow < jMax; jNow++) {
+
+ // Final-state and coloured jNow only
+ if (!e[jNow].isFinal() || e[jNow].colType() == 0) continue;
+
+ // POWHEG
+ if (pTdefMode == 0 || pTdefMode == 1) {
+
+ // ISR - only done once as just kinematical pT
+ if (!FSR) {
+ pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
+
+ // FSR - try all outgoing partons from system before branching
+ // as i. Note that for the hard system, there is no
+ // "before branching" information.
+ } else {
+
+ int outSize = partonSystemsPtr->sizeOut(0);
+ for (int iMem = 0; iMem < outSize; iMem++) {
+ int iNow = partonSystemsPtr->getOut(0, iMem);
+
+ // Coloured only, i != jNow and no carbon copies
+ if (iNow == jNow || e[iNow].colType() == 0) continue;
+ if (jNow == e[iNow].daughter1()
+ && jNow == e[iNow].daughter2()) continue;
+
+ pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
+ ? false : FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0)
+ ? pTnow : min(pTemt, pTnow);
+ } // for (iMem)
+
+ } // if (!FSR)
+
+ // Pythia
+ } else if (pTdefMode == 2) {
+
+ // ISR - other incoming as recoiler
+ if (!FSR) {
+ pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
+ pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
+
+ // FSR - try all final-state coloured partons as radiator
+ // after emission (k).
+ } else {
+ for (int kNow = 0; kNow < e.size(); kNow++) {
+ if (kNow == jNow || !e[kNow].isFinal() ||
+ e[kNow].colType() == 0) continue;
+
+ // For this kNow, need to have a recoiler.
+ // Try two incoming.
+ pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0)
+ ? pTnow : min(pTemt, pTnow);
+ pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0)
+ ? pTnow : min(pTemt, pTnow);
+
+ // Try all other outgoing.
+ for (int rNow = 0; rNow < e.size(); rNow++) {
+ if (rNow == kNow || rNow == jNow ||
+ !e[rNow].isFinal() || e[rNow].colType() == 0) continue;
+ pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
+ if (pTnow > 0.) pTemt = (pTemt < 0)
+ ? pTnow : min(pTemt, pTnow);
+ } // for (rNow)
+
+ } // for (kNow)
+ } // if (!FSR)
+ } // if (pTdefMode)
+ } // for (j)
+ }
+ } // for (xSR)
+
+#ifdef DBGOUTPUT
+ cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k
+ << ", r = " << r << ", xSR = " << xSRin
+ << ", pTemt = " << pTemt << endl;
+#endif
+
+ return pTemt;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // Extraction of pThard based on the incoming event.
+ // Assume that all the final-state particles are in a continuous block
+ // at the end of the event and the final entry is the POWHEG emission.
+ // If there is no POWHEG emission, then pThard is set to Qfac.
+
+ bool canVetoMPIStep() { return true; }
+ int numberVetoMPIStep() { return 1; }
+ bool doVetoMPIStep(int nMPI, const Event &e) {
+ // Extra check on nMPI
+ if (nMPI > 1) return false;
+
+ // Find if there is a POWHEG emission. Go backwards through the
+ // event record until there is a non-final particle. Also sum pT and
+ // find pT_1 for possible MPI vetoing
+ int count = 0;
+ double pT1 = 0., pTsum = 0.;
+ for (int i = e.size() - 1; i > 0; i--) {
+ if (e[i].isFinal()) {
+ count++;
+ pT1 = e[i].pT();
+ pTsum += e[i].pT();
+ } else break;
+ }
+ // Extra check that we have the correct final state
+ if (count != nFinal && count != nFinal + 1) {
+ cout << "Error: wrong number of final state particles in event" << endl;
+ exit(1);
+ }
+ // Flag if POWHEG radiation present and index
+ bool isEmt = (count == nFinal) ? false : true;
+ int iEmt = (isEmt) ? e.size() - 1 : -1;
+
+ // If there is no radiation or if pThardMode is 0 then set pThard to Qfac.
+ if (!isEmt || pThardMode == 0) {
+ pThard = infoPtr->QFac();
+
+ // If pThardMode is 1 then the pT of the POWHEG emission is checked against
+ // all other incoming and outgoing partons, with the minimal value taken
+ } else if (pThardMode == 1) {
+ pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
+
+ // If pThardMode is 2, then the pT of all final-state partons is checked
+ // against all other incoming and outgoing partons, with the minimal value
+ // taken
+ } else if (pThardMode == 2) {
+ pThard = pTcalc(e, -1, -1, -1, -1, -1);
+
+ }
+
+ // Find MPI veto pT if necessary
+ if (MPIvetoMode == 1) {
+ pTMPI = (isEmt) ? pTsum / 2. : pT1;
+ }
+
+#ifdef DBGOUTPUT
+ cout << "doVetoMPIStep: Qfac = " << infoPtr->QFac()
+ << ", pThard = " << pThard << endl << endl;
+#endif
+
+// std::cout<<"vetoScale = "<<pThard<<std::endl;
+
+ // Initialise other variables
+ accepted = false;
+ nAcceptSeq = nISRveto = nFSRveto = 0;
+
+ // Do not veto the event
+ return false;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // ISR veto
+
+ bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
+ bool doVetoISREmission(int, const Event &e, int iSys) {
+ // Must be radiation from the hard system
+ if (iSys != 0) return false;
+
+ // If we already have accepted 'vetoCount' emissions in a row, do nothing
+ if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
+
+ // Pythia radiator after, emitted and recoiler after.
+ int iRadAft = -1, iEmt = -1, iRecAft = -1;
+ for (int i = e.size() - 1; i > 0; i--) {
+ if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
+ else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
+ else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
+ if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
+ }
+ if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
+ e.list();
+ cout << "Error: couldn't find Pythia ISR emission" << endl;
+ exit(1);
+ }
+
+ // pTemtMode == 0: pT of emitted w.r.t. radiator
+ // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
+ // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
+ int xSR = (pTemtMode == 0) ? 0 : -1;
+ int i = (pTemtMode == 0) ? iRadAft : -1;
+ int j = (pTemtMode != 2) ? iEmt : -1;
+ int k = -1;
+ int r = (pTemtMode == 0) ? iRecAft : -1;
+ double pTemt = pTcalc(e, i, j, k, r, xSR);
+
+#ifdef DBGOUTPUT
+ cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl;
+#endif
+
+ // Veto if pTemt > pThard
+ if (pTemt > pThard) {
+ nAcceptSeq = 0;
+ nISRveto++;
+ return true;
+ }
+
+ // Else mark that an emission has been accepted and continue
+ nAcceptSeq++;
+ accepted = true;
+ return false;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // FSR veto
+
+ bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
+ bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
+ // Must be radiation from the hard system
+ if (iSys != 0) return false;
+
+ // If we already have accepted 'vetoCount' emissions in a row, do nothing
+ if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
+
+ // Pythia radiator (before and after), emitted and recoiler (after)
+ int iRecAft = e.size() - 1;
+ int iEmt = e.size() - 2;
+ int iRadAft = e.size() - 3;
+ int iRadBef = e[iEmt].mother1();
+ if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
+ e[iEmt].status() != 51 || e[iRadAft].status() != 51) {
+ e.list();
+ cout << "Error: couldn't find Pythia FSR emission" << endl;
+ exit(1);
+ }
+
+ // Behaviour based on pTemtMode:
+ // 0 - pT of emitted w.r.t. radiator before
+ // 1 - min(pT of emitted w.r.t. all incoming/outgoing)
+ // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
+ int xSR = (pTemtMode == 0) ? 1 : -1;
+ int i = (pTemtMode == 0) ? iRadBef : -1;
+ int k = (pTemtMode == 0) ? iRadAft : -1;
+ int r = (pTemtMode == 0) ? iRecAft : -1;
+
+ // When pTemtMode is 0 or 1, iEmt has been selected
+ double pTemt = 0.;
+ if (pTemtMode == 0 || pTemtMode == 1) {
+ // Which parton is emitted, based on emittedMode:
+ // 0 - Pythia definition of emitted
+ // 1 - Pythia definition of radiated after emission
+ // 2 - Random selection of emitted or radiated after emission
+ // 3 - Try both emitted and radiated after emission
+ int j = iRadAft;
+ if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
+
+ for (int jLoop = 0; jLoop < 2; jLoop++) {
+ if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
+ else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
+
+ // For emittedMode == 3, have tried iRadAft, now try iEmt
+ if (emittedMode != 3) break;
+ if (k != -1) swap(j, k); else j = iEmt;
+ }
+
+ // If pTemtMode is 2, then try all final-state partons as emitted
+ } else if (pTemtMode == 2) {
+ pTemt = pTcalc(e, i, -1, k, r, xSR);
+
+ }
+
+#ifdef DBGOUTPUT
+ cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl;
+#endif
+
+ // Veto if pTemt > pThard
+ if (pTemt > pThard) {
+ nAcceptSeq = 0;
+ nFSRveto++;
+ return true;
+ }
+
+ // Else mark that an emission has been accepted and continue
+ nAcceptSeq++;
+ accepted = true;
+ return false;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // MPI veto
+
+ bool canVetoMPIEmission() { return (MPIvetoMode == 0) ? false : true; }
+ bool doVetoMPIEmission(int, const Event &e) {
+ if (MPIvetoMode == 1) {
+ if (e[e.size() - 1].pT() > pTMPI) {
+#ifdef DBGOUTPUT
+ cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
+ << ", pTMPI = " << pTMPI << endl << endl;
+#endif
+ return true;
+ }
+ }
+ return false;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // Functions to return information
+
+ int getNISRveto() { return nISRveto; }
+ int getNFSRveto() { return nFSRveto; }
+
+ private:
+ int nFinal, vetoMode, vetoCount, pThardMode, pTemtMode,
+ emittedMode, pTdefMode, MPIvetoMode;
+ double pThard, pTMPI;
+ bool accepted;
+ // The number of accepted emissions (in a row)
+ int nAcceptSeq;
+ // Statistics on vetos
+ unsigned long int nISRveto, nFSRveto;
+
+ };
+}
diff --git a/Generators/Pythia8_i/src/UserHooksFactory.cxx b/Generators/Pythia8_i/src/UserHooksFactory.cxx
new file mode 100644
index 00000000000..1bddd3c8c27
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserHooksFactory.cxx
@@ -0,0 +1,27 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserHooksFactory.h"
+
+#include <stdexcept>
+
+namespace Pythia8_UserHooks{
+
+ UserHooks *UserHooksFactory::create(const string &name){
+ map<string, const ICreator*>::const_iterator it = s_creators().find(name);
+ if(it == s_creators().end()){
+ //eek!
+ throw std::runtime_error("Pythia 8 UserHooksFactory: cannot create user hook " + name);
+ }
+ return it->second->create();
+ }
+
+ ///static function to instantiate map of string name Vs. creator object on first use
+ map<string, const UserHooksFactory::ICreator*> &UserHooksFactory::s_creators(){
+ static map<string, const UserHooksFactory::ICreator*> creators;
+ return creators;
+ }
+
+}
+
diff --git a/Generators/Pythia8_i/src/UserProcessFactory.cxx b/Generators/Pythia8_i/src/UserProcessFactory.cxx
new file mode 100644
index 00000000000..5dbacc6888b
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserProcessFactory.cxx
@@ -0,0 +1,28 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserProcessFactory.h"
+
+#include <stdexcept>
+
+namespace Pythia8_UserProcess{
+
+ Sigma2Process *UserProcessFactory::create(const string &name){
+ map<string, const ICreator*>::const_iterator it = s_creators().find(name);
+ if(it == s_creators().end()){
+ //eek!
+ throw std::runtime_error("Pythia 8 UserProcessFactory: cannot create user process " + name);
+ }
+
+ return it->second->create();
+ }
+
+ ///static function to instantiate map of string name Vs. creator object on first use
+ map<string, const UserProcessFactory::ICreator*> &UserProcessFactory::s_creators(){
+ static map<string, const UserProcessFactory::ICreator*> creators;
+ return creators;
+ }
+
+}
+
diff --git a/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2emu.cxx b/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2emu.cxx
new file mode 100644
index 00000000000..f9fb4b72e4d
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2emu.cxx
@@ -0,0 +1,127 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+ * Lepton flavour violating differential cross sections
+ */
+
+#include "Pythia8_i/UserProcessFactory.h"
+
+// calls to fortran routines
+#include "CLHEP/Random/RandFlat.h"
+
+#include <iostream>
+#include <vector>
+
+namespace Pythia8{
+ class Sigma2qqbar2emu;
+}
+
+// This one line registers this process with the factory
+// An instance of this process can be created with
+// UserProcessFactory::create("qqbar2emu")
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2emu> qqbar2emuCreator("qqbar2emu");
+
+namespace Pythia8 {
+
+ /**
+ *
+ * Sigma1ql2LFV class.
+ * Cross section for q qbar -> e mu (t:squark).
+ */
+
+ class Sigma2qqbar2emu : public Sigma2Process{
+
+ public:
+ //Constructor.
+ Sigma2qqbar2emu()
+ : lambdaCoup(0.), sigmaTemp(0.), m_stop(0.),
+ lambdaSqddbar(0.), lambdaSqdsbar(0.), lambdaSqsdbar(0.), lambdaSqssbar(0.),
+ Nc(0), myRand(0)
+ {}
+
+
+ /// Initialize process.
+ void initProc() {
+
+ // Store parameters.
+ lambdaCoup = 0.2;
+ Nc = 3;
+ m_stop = particleDataPtr->m0(1000006); //95.5;
+ lambdaSqddbar = 0.065;
+ lambdaSqdsbar = 3.8e-5;
+ lambdaSqsdbar = 3.8e-5;
+ lambdaSqssbar = 0.065;
+ myRand = new Rndm(1234567);
+
+ return;
+ }
+
+ //*********
+
+ /// Evaluate sigmaHat(sHat), part independent of incoming flavour.
+ void sigmaKin() {
+
+ // The 2 out front is for e- mu+ plus e+ mu -
+ sigmaTemp = 2./((float)Nc * 64. * M_PI);
+
+ return;
+ }
+
+ //*********
+
+ /// Evaluate sigmaHat(sHat), part dependent of incoming flavour.
+ double sigmaHat() {
+
+ double lambdaSq = 0.;
+ if((id1 == 1 && id2 == -1) || (id1 == -1 && id2 == 1)) lambdaSq = lambdaSqddbar;
+ else if((id1 == 1 && id2 == -3) || (id1 == -3 && id2 == 1)) lambdaSq = lambdaSqdsbar;
+ else if((id1 == 3 && id2 == -1) || (id1 == -1 && id2 == 3)) lambdaSq = lambdaSqsdbar;
+ else if((id1 == 3 && id2 == -3) || (id1 == -3 && id2 == 3)) lambdaSq = lambdaSqssbar;
+ double sigma = pow2(lambdaSq)*sigmaTemp*tH2/(sH2*pow2(tH-pow2(m_stop)));
+ return sigma;
+ }
+
+ //*********
+
+ // Select identity, colour and anticolour.
+
+ void setIdColAcol() {
+
+ // Flavours.
+ double x = myRand->flat();
+ if(x < 0.5) setId( id1, id2, 11, -13);
+ else setId( id1, id2, -11, 13);
+
+ // tH defined between f and LQ: must swap tHat <-> uHat if qbar q in.
+ swapTU = (id1 < 0);
+
+ //cab I don't understand what's done here.
+ // Colour flow topologies.
+ if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
+ else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
+
+ return;
+ }
+
+ // Info on the subprocess.
+ string name() const {return "q qbar -> e mu (t:squark)";}
+ //cab Need Codes for q qbar -> e+- mu-+
+ int code() const {return 10001;}
+ string inFlux() const {return "qq";}
+ //cab It looks like here is where the final state is specified
+ //cab electrons are +-11, muons are +-13
+ int id3Mass() const {return 11;}
+ int id4Mass() const {return 13;}
+
+ private:
+
+ // Parameters set at initialization or for current kinematics.
+ double lambdaCoup, sigmaTemp, m_stop;
+ double lambdaSqddbar, lambdaSqdsbar, lambdaSqsdbar, lambdaSqssbar;
+ int Nc;
+ Rndm* myRand;
+
+ };
+}
diff --git a/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlBar.cxx b/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlBar.cxx
new file mode 100644
index 00000000000..b3d362e651c
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlBar.cxx
@@ -0,0 +1,187 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Single excited lepton production - similar to original pythia8 code,
+// but polarization of final state is reversed [to match Calchep distributions
+
+#include "Pythia8_i/UserProcessFactory.h"
+
+#include <iostream>
+#include <vector>
+
+namespace Pythia8 {
+ class Sigma2qqbar2lStarlBar;
+ class Sigma2qqbar2nueStarnueBar;
+ class Sigma2qqbar2numuStarnumuBar;
+ class Sigma2qqbar2nutauStarnutauBar;
+ class Sigma2qqbar2eStareBar;
+ class Sigma2qqbar2muStarmuBar;
+ class Sigma2qqbar2tauStartauBar;
+}
+
+// This one line registers this process with the factory
+// An instance of this process can be created with
+// UserProcessFactory::create("qqbar2lStarlBar")
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2nueStarnueBar> qqbar2nueStarnueBarCreator( "qqbar2nueStarnueBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2numuStarnumuBar> qqbar2numuStarnumuBarCreator( "qqbar2numuStarnumuBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2nutauStarnutauBar> qqbar2nutauStarnutauBarCreator("qqbar2nutauStarnutauBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2eStareBar> qqbar2eStareBarCreator( "qqbar2eStareBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2muStarmuBar> qqbar2muStarmuBarCreator( "qqbar2muStarmuBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2tauStartauBar> qqbar2tauStartauBarCreator( "qqbar2tauStartauBar");
+
+namespace Pythia8 {
+
+ //**************************************************************************
+
+ // Sigma2qqbar2lStarlBar class.
+ // Cross section for q qbar -> lStar lBar
+
+ class Sigma2qqbar2lStarlBar: public Sigma2Process {
+
+ public:
+ //Constructor.
+ Sigma2qqbar2lStarlBar(int idlIn) : idl(idlIn),
+ idRes(0), codeSave(0),
+ Lambda(0.), preFac(0.), openFracPos(0.), openFracNeg(0.), sigma(0.), m2ResTimes4(0.){};
+
+ //*********
+
+ // Initialize process.
+
+ void initProc() {
+
+ // Set up process properties from the chosen lepton flavour.
+ idRes = 4000000 + idl;
+ codeSave = 4060 + idl;
+ if (idl == 11) nameSave = "q qbar -> e^*+- e^-+";
+ else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_ebar";
+ else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^-+";
+ else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mubar";
+ else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^-+";
+ else nameSave = "q qbar -> nu_tau^* nu_taubar";
+
+ // Secondary open width fractions.
+ openFracPos = particleDataPtr->resOpenFrac( idRes);
+ openFracNeg = particleDataPtr->resOpenFrac(-idRes);
+
+ // Locally stored properties and couplings.
+ Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
+ preFac = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 3.; // orig pythia
+ //preFac = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 6.;
+
+ }
+
+ //*********
+ // Evaluate sigmaHat(sHat).
+ virtual double sigmaHat() {return sigma;}
+
+
+ // Evaluate sigmaHat(sHat), part independent of incoming flavour.
+
+ void sigmaKin() {
+
+ // later this gets multiplied further with sH*(1-s3/sH)
+
+ sigma = preFac * (-uH) * (sH + tH) / sH2; // orig pythia
+ //sigma = preFac * (2. + s3/sH) *(1 - s3/sH) ;
+ //Baur sigma = preFac * ( 1 + 1./3.*(sH - s3)/(sH+s3) ) * (1 - s3/sH) * (1 + s3/sH);
+ }
+
+ //*********
+
+ // Select identity, colour and anticolour.
+
+ void setIdColAcol() {
+
+ // Flavours: either lepton or antilepton may be excited.
+ if (rndmPtr->flat() * (openFracPos + openFracNeg) < openFracPos) {
+ setId( id1, id2, idRes, -idl);
+ if (id1 < 0) swapTU = true;
+ } else {
+ setId( id1, id2, -idRes, idl);
+ if (id1 > 0) swapTU = true;
+ }
+
+ // Colour flow trivial.
+ if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
+ else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
+
+ }
+
+ //**************************************************************************
+
+ double weightDecay( Event& process, int iResBeg,
+ int iResEnd) {
+
+ // l* should sit in entry 5. and 6 Sequential Z/W decay assumed isotropic.
+ if (iResBeg != 5 && iResEnd != 5 && iResBeg != 6 && iResEnd != 6 ) return 1.;
+
+ // Phase space factors.
+ double mr1 = pow2(process[7].m() / process[5].m());
+ double mr2 = pow2(process[8].m() / process[5].m());
+
+ // Reconstruct decay angle in l* CoM frame.
+ int idAbs3 = process[7].idAbs();
+ Vec4 pLStarCom = (idAbs3 < 20) ? process[8].p() : process[7].p();
+ //Vec4 pLStarCom = (idAbs3 < 20) ? process[7].p() : process[8].p(); // orig pythia
+ pLStarCom.bstback(process[5].p());
+ double cosThe = costheta(pLStarCom, process[5].p());
+ double wt = 1.;
+
+ // Decay, l* -> l + gamma/Z^0/W^+-).
+ int idBoson = (idAbs3 < 20) ? process[8].idAbs() : process[7].idAbs();
+ if (idBoson == 22) {
+ wt = 0.5 * (1. + cosThe);
+ } else if (idBoson == 23 || idBoson == 24) {
+ double mrB = (idAbs3 < 20) ? mr2 : mr1;
+ double kTrm = 0.5 * (mrB * (1. - cosThe));
+ wt = (1. + cosThe + kTrm) / (2 + mrB);
+ }
+
+ // Done.
+ return wt;
+ }
+
+ // Info on the subprocess.
+ virtual string name() const {return nameSave;}
+ virtual int code() const {return codeSave;}
+ virtual string inFlux() const {return "qqbarSame";}
+ virtual int id3Mass() const {return abs(idRes);}
+
+ //virtual bool convertM2() const {return false;} //true is only suitable for resonance production pp->A->ll
+ //virtual bool convert2mb() const {return true;}
+
+ // Parameters set at initialization or for current kinematics.
+ int idl, idRes, codeSave;
+ string nameSave;
+ double Lambda, preFac, openFracPos, openFracNeg, sigma;
+ double m2ResTimes4;
+ };
+ //**************************************************************************
+ class Sigma2qqbar2nueStarnueBar: public Sigma2qqbar2lStarlBar {
+ public:
+ Sigma2qqbar2nueStarnueBar() : Sigma2qqbar2lStarlBar(12) {}
+ };
+ class Sigma2qqbar2numuStarnumuBar: public Sigma2qqbar2lStarlBar {
+ public:
+ Sigma2qqbar2numuStarnumuBar() : Sigma2qqbar2lStarlBar(14) {}
+ };
+ class Sigma2qqbar2nutauStarnutauBar: public Sigma2qqbar2lStarlBar {
+ public:
+ Sigma2qqbar2nutauStarnutauBar() : Sigma2qqbar2lStarlBar(16) {}
+ };
+ class Sigma2qqbar2eStareBar: public Sigma2qqbar2lStarlBar {
+ public:
+ Sigma2qqbar2eStareBar() : Sigma2qqbar2lStarlBar(11) {}
+ };
+ class Sigma2qqbar2muStarmuBar: public Sigma2qqbar2lStarlBar {
+ public:
+ Sigma2qqbar2muStarmuBar() : Sigma2qqbar2lStarlBar(13) {}
+ };
+ class Sigma2qqbar2tauStartauBar: public Sigma2qqbar2lStarlBar {
+ public:
+ Sigma2qqbar2tauStartauBar() : Sigma2qqbar2lStarlBar(15) {}
+ };
+
+}
diff --git a/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlStarBar.cxx b/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlStarBar.cxx
new file mode 100644
index 00000000000..02c0d32dfa6
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserProcesses/Sigma2qqbar2lStarlStarBar.cxx
@@ -0,0 +1,186 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// Double excited lepton production
+
+#include "Pythia8_i/UserProcessFactory.h"
+
+#include <iostream>
+#include <vector>
+
+namespace Pythia8 {
+ class Sigma2qqbar2lStarlStarBar;
+ class Sigma2qqbar2nueStarnueStarBar;
+ class Sigma2qqbar2numuStarnumuStarBar;
+ class Sigma2qqbar2nutauStarnutauStarBar;
+ class Sigma2qqbar2eStareStarBar;
+ class Sigma2qqbar2muStarmuStarBar;
+ class Sigma2qqbar2tauStartauStarBar;
+}
+
+// This one line registers this process with the factory
+// An instance of this process can be created with
+// UserProcessFactory::create("qqbar2lStarlStarBar")
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2nueStarnueStarBar> qqbar2nueStarlnuetarBarCreator( "qqbar2nueStarnueStarBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2numuStarnumuStarBar> qqbar2numuStarlnumutarBarCreator( "qqbar2numuStarnumuStarBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2nutauStarnutauStarBar> qqbar2nutauStarlnutautarBarCreator("qqbar2nutauStarnutauStarBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2eStareStarBar> qqbar2eStarletarBarCreator( "qqbar2eStareStarBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2muStarmuStarBar> qqbar2muStarlmutarBarCreator( "qqbar2muStarmuStarBar");
+Pythia8_UserProcess::UserProcessFactory::Creator<Pythia8::Sigma2qqbar2tauStartauStarBar> qqbar2tauStarltautarBarCreator( "qqbar2tauStartauStarBar");
+
+namespace Pythia8 {
+
+ //**************************************************************************
+
+ // Sigma2qqbar2lStarlStarBar class.
+ // Cross section for q qbar -> lStar lStarBar
+
+ // A derived class for q qbar -> lStar lStarBar
+
+ class Sigma2qqbar2lStarlStarBar: public Sigma2Process {
+
+ public:
+ //Constructor.
+ Sigma2qqbar2lStarlStarBar(int idlIn) : idl(idlIn),
+ idRes(0), codeSave(0),
+ Lambda(0.), preFac(0.), openFracPos(0.), openFracNeg(0.), sigma(0.),
+ m2ResTimes4(0.){};
+
+ // Initialize process.
+
+ void initProc() {
+
+ // Set up process properties from the chosen lepton flavour.
+ idRes = 4000000 + idl;
+ codeSave = 4040 + idl;
+ if (idl == 11) nameSave = "q qbar -> e^*+- e^*-+";
+ else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_e^*bar";
+ else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^*-+";
+ else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mu^*bar";
+ else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^*-+";
+ else nameSave = "q qbar -> nu_tau^* nu_tau^*bar";
+
+ // Secondary open width fractions.
+ openFracPos = particleDataPtr->resOpenFrac( idRes);
+ openFracNeg = particleDataPtr->resOpenFrac(-idRes);
+
+ // Locally stored properties and couplings.
+ Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
+ preFac = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 9.;
+
+ double mRes = particleDataPtr->m0(idRes);
+ m2ResTimes4 = mRes*mRes*4.;
+ }
+
+ //*********
+
+ // Evaluate sigmaHat(sHat).
+ virtual double sigmaHat() {return sigma;}
+
+ // Evaluate sigmaHat(sHat), part independent of incoming flavour.
+
+ void sigmaKin() {
+
+ // use assumption that sigma is multiplied later with sH*(1-s3/sH)
+
+ //sigma = preFac * sqrtpos( 1.- m2ResTimes4/sH );
+ sigma = preFac * sqrtpos( 1.- 4.*s3/sH );
+
+ }
+
+ //*********
+
+ // Select identity, colour and anticolour.
+
+ void setIdColAcol() {
+
+ // Flavours: both lepton and antilepton are excited.
+ setId( id1, id2, idRes, -idRes);
+ //if (id1 < 0) swapTU = true; // OI how is this used?
+
+ // Colour flow trivial.
+ if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
+ else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
+
+ }
+
+ //**************************************************************************
+
+ double weightDecay( Event& process, int iResBeg,
+ int iResEnd) {
+
+ // l* should sit in entry 5. and 6 Sequential Z/W decay assumed isotropic.
+ if (iResBeg != 5 && iResEnd != 5 && iResBeg != 6 && iResEnd != 6 ) return 1.;
+
+ // Phase space factors.
+ double mr1 = pow2(process[7].m() / process[5].m());
+ double mr2 = pow2(process[8].m() / process[5].m());
+
+ // Reconstruct decay angle in l* CoM frame.
+ int idAbs3 = process[7].idAbs();
+ Vec4 pLStarCom = (idAbs3 < 20) ? process[8].p() : process[7].p();
+ pLStarCom.bstback(process[5].p());
+ double cosThe = costheta(pLStarCom, process[5].p());
+ double wt = 1.;
+
+ // Decay, l* -> l + gamma/Z^0/W^+-).
+ int idBoson = (idAbs3 < 20) ? process[8].idAbs() : process[7].idAbs();
+ if (idBoson == 22) {
+ wt = 0.5 * (1. + cosThe);
+ } else if (idBoson == 23 || idBoson == 24) {
+ double mrB = (idAbs3 < 20) ? mr2 : mr1;
+ double kTrm = 0.5 * (mrB * (1. - cosThe));
+ wt = (1. + cosThe + kTrm) / (2 + mrB);
+ }
+
+ // Done.
+ return wt;
+ }
+
+ // Info on the subprocess.
+ virtual string name() const {return nameSave;}
+ virtual int code() const {return codeSave;}
+ virtual string inFlux() const {return "qqbarSame";}
+ virtual int id3Mass() const {return abs(idRes);}
+ virtual int id4Mass() const {return abs(idRes);}
+
+ virtual bool convertM2() const {return false;} //true is only suitable for resonance production pp->A->ll
+ virtual bool convert2mb() const {return true;}
+
+ // Parameters set at initialization or for current kinematics.
+ int idl, idRes, codeSave;
+ string nameSave;
+ double Lambda, preFac, openFracPos, openFracNeg, sigma;
+ double m2ResTimes4;
+ };
+
+
+ //**************************************************************************
+ class Sigma2qqbar2nueStarnueStarBar: public Sigma2qqbar2lStarlStarBar {
+ public:
+ Sigma2qqbar2nueStarnueStarBar() : Sigma2qqbar2lStarlStarBar(12) {}
+ };
+ class Sigma2qqbar2numuStarnumuStarBar: public Sigma2qqbar2lStarlStarBar {
+ public:
+ Sigma2qqbar2numuStarnumuStarBar() : Sigma2qqbar2lStarlStarBar(14) {}
+ };
+ class Sigma2qqbar2nutauStarnutauStarBar: public Sigma2qqbar2lStarlStarBar {
+ public:
+ Sigma2qqbar2nutauStarnutauStarBar() : Sigma2qqbar2lStarlStarBar(16) {}
+ };
+ class Sigma2qqbar2eStareStarBar: public Sigma2qqbar2lStarlStarBar {
+ public:
+ Sigma2qqbar2eStareStarBar() : Sigma2qqbar2lStarlStarBar(11) {}
+ };
+ class Sigma2qqbar2muStarmuStarBar: public Sigma2qqbar2lStarlStarBar {
+ public:
+ Sigma2qqbar2muStarmuStarBar() : Sigma2qqbar2lStarlStarBar(13) {}
+ };
+
+ class Sigma2qqbar2tauStartauStarBar: public Sigma2qqbar2lStarlStarBar {
+ public:
+ Sigma2qqbar2tauStartauStarBar() : Sigma2qqbar2lStarlStarBar(15) {}
+ };
+
+}
diff --git a/Generators/Pythia8_i/src/UserResonanceFactory.cxx b/Generators/Pythia8_i/src/UserResonanceFactory.cxx
new file mode 100644
index 00000000000..04d513c862d
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserResonanceFactory.cxx
@@ -0,0 +1,33 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "Pythia8_i/UserResonanceFactory.h"
+
+#include <boost/lexical_cast.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <vector>
+#include <stdexcept>
+
+namespace Pythia8_UserResonance{
+
+ using std::vector;
+
+ ResonanceWidths *UserResonanceFactory::create(const string &name, int pdgid){
+
+ map<string, const ICreator*>::const_iterator it = s_creators().find(name);
+ if(it == s_creators().end()){
+ //eek!
+ throw std::runtime_error("Pythia 8 UserResonanceFactory: cannot create user resonance " + name);
+ }
+ return it->second->create(pdgid);
+ }
+
+ ///static function to instantiate map of string name Vs. creator object on first use
+ map<string, const UserResonanceFactory::ICreator*> &UserResonanceFactory::s_creators(){
+ static map<string, const UserResonanceFactory::ICreator*> creators;
+ return creators;
+ }
+
+}
\ No newline at end of file
diff --git a/Generators/Pythia8_i/src/UserResonances/ResonanceExcitedCI.cxx b/Generators/Pythia8_i/src/UserResonances/ResonanceExcitedCI.cxx
new file mode 100644
index 00000000000..9d9a45f40f3
--- /dev/null
+++ b/Generators/Pythia8_i/src/UserResonances/ResonanceExcitedCI.cxx
@@ -0,0 +1,114 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// The ResonanceExcited class handles excited-fermion resonances.
+// This is a copy of existing pythia ResonanceExcited class
+// with extention to include decays going via contact interactions (4-fermion vertexes)
+//
+// author Olya Igonkina, modified by James Monk for Pythia8_i
+
+#include "Pythia8_i/UserResonanceFactory.h"
+
+namespace Pythia8{
+ class ResonanceExcitedCI;
+}
+
+Pythia8_UserResonance::UserResonanceFactory::Creator<Pythia8::ResonanceExcitedCI> resonanceExcitedCreator("ExcitedCI");
+
+namespace Pythia8{
+
+ class ResonanceExcitedCI : public ResonanceWidths {
+
+ public:
+
+ // Constructor.
+ ResonanceExcitedCI(int idResIn):
+ m_lambda(0.), m_coupF(0.), m_coupFprime(0.), m_coupFcol(0.),
+ m_sin2tW(0.), m_cos2tW(0.){
+ initBasic(idResIn); std::cout << " ResonanceExcitedCI constructor\n";
+ }
+
+ private:
+
+ void initConstants() {
+
+ // Locally stored properties and couplings.
+ m_lambda = settingsPtr->parm("ExcitedFermion:Lambda");
+ m_coupF = settingsPtr->parm("ExcitedFermion:coupF");
+ m_coupFprime = settingsPtr->parm("ExcitedFermion:coupFprime");
+ m_coupFcol = settingsPtr->parm("ExcitedFermion:coupFcol");
+ m_sin2tW = couplingsPtr->sin2thetaW();
+ m_cos2tW = 1. - m_sin2tW;
+
+ return;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // Calculate various common prefactors for the current mass.
+
+ void calcPreFac(bool) {
+ // Common coupling factors.
+ alpEM = couplingsPtr->alphaEM(mHat * mHat);
+ alpS = couplingsPtr->alphaS(mHat * mHat);
+ preFac = pow3(mHat) / pow2(m_lambda);
+
+ return;
+ }
+
+ //--------------------------------------------------------------------------
+
+ // Calculate width for currently considered channel.
+
+ void calcWidth(bool) {
+
+ // Check that above threshold.
+ if (ps == 0.) return;
+
+ // f^* -> f g.
+ if (id1Abs == 21) widNow = preFac * alpS * pow2(m_coupFcol) / 3.;
+
+ // f^* -> f gamma.
+ else if (id1Abs == 22) {
+ double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
+ double chgY = (id2Abs < 9) ? 1. / 6. : -0.5;
+ double chg = chgI3 * m_coupF + chgY * m_coupFprime;
+ widNow = preFac * alpEM * pow2(chg) / 4.;
+ }
+ // f^* -> f Z^0.
+ else if (id1Abs == 23) {
+ double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
+ double chgY = (id2Abs < 9) ? 1. / 6. : -0.5;
+ double chg = chgI3 * m_cos2tW * m_coupF - chgY * m_sin2tW * m_coupFprime;
+ widNow = preFac * (alpEM * pow2(chg) / (8. * m_sin2tW * m_cos2tW))
+ * ps*ps * (2. + mr1);
+ }
+
+ // f^* -> f' W^+-.
+ else if (id1Abs == 24) widNow = preFac * (alpEM * pow2(m_coupF)
+ / (16. * m_sin2tW)) * ps*ps * (2. + mr1);
+ // idxAbs are sorted according to their abs value ; use idRes
+
+ // f^* -> f f' \bar{f'}
+ else {
+ //std::cout << " hei! in ResonanceExcitedCI processing F* -> f f'f'bar decays "<< id1Abs <<" "<<id2Abs<<" "<<id3Abs<< "\n";
+
+ if( id1Abs < 17 && id2Abs < 17 && id3Abs>0 && id3Abs < 17 ){ // f* -> f f'\bar(f'} decay
+ widNow = preFac * pow2(mHat) / pow2(m_lambda) / 96. / M_PI ;
+ if( id3Abs < 10 ) widNow *= 3.; // quarks in final state. id1 is highest pdgId of f, f'
+ if( id1Abs == id2Abs && id1Abs == id3Abs ){
+ if( abs(idRes)-4000000 < 10 ) widNow *= 4./3.;
+ else widNow *= 2.;
+ }
+ }
+ }
+
+ return;
+ }
+
+ // Locally stored properties and couplings.
+ double m_lambda, m_coupF, m_coupFprime, m_coupFcol, m_sin2tW, m_cos2tW;
+
+ };
+}
diff --git a/Generators/Pythia8_i/src/components/Pythia8_i_entries.cxx b/Generators/Pythia8_i/src/components/Pythia8_i_entries.cxx
new file mode 100644
index 00000000000..562b20b314b
--- /dev/null
+++ b/Generators/Pythia8_i/src/components/Pythia8_i_entries.cxx
@@ -0,0 +1,11 @@
+#include "Pythia8_i/Pythia8_i.h"
+#include "GaudiKernel/DeclareFactoryEntries.h"
+
+/**
+ * author: James Monk (jmonk@hep.ucl.ac.uk)
+ */
+
+DECLARE_ALGORITHM_FACTORY( Pythia8_i )
+DECLARE_FACTORY_ENTRIES( Pythia8_i ){
+ DECLARE_ALGORITHM( Pythia8_i )
+}
diff --git a/Generators/Pythia8_i/src/components/Pythia8_i_load.cxx b/Generators/Pythia8_i/src/components/Pythia8_i_load.cxx
new file mode 100644
index 00000000000..4a1c8d8acff
--- /dev/null
+++ b/Generators/Pythia8_i/src/components/Pythia8_i_load.cxx
@@ -0,0 +1,7 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+/**
+ * author: James Monk (jmonk@hep.ucl.ac.uk)
+ */
+
+LOAD_FACTORY_ENTRIES( Pythia8_i )
--
GitLab