Page MenuHomeHEPForge
Files F7879389

No OneTemporary
Actions

View Options

diff --git a/Contrib/FxFx/FxFxHandler.h b/Contrib/FxFx/FxFxHandler.h
--- a/Contrib/FxFx/FxFxHandler.h
+++ b/Contrib/FxFx/FxFxHandler.h
@@ -1,640 +1,640 @@
// -*- C++ -*-
#ifndef HERWIG_FxFxHandler_H
#define HERWIG_FxFxHandler_H
//
// This is the declaration of the FxFxHandler class.
//
#include "Herwig/Shower/QTilde/QTildeShowerHandler.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/Config/Pointers.h"
-#include "Herwig/Shower/QTilde/Couplings/ShowerAlphaQCD.h"
+#include "Herwig/Shower/Core/Couplings/ShowerAlpha.h"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "ThePEG/Utilities/CompSelector.h"
#include "ThePEG/Utilities/XSecStat.h"
namespace Herwig {
class FxFxHandler;
}
//declaration of thepeg ptr
namespace ThePEG {
ThePEG_DECLARE_POINTERS(Herwig::FxFxHandler,FxFxHandlerPtr);
}
namespace Herwig {
using namespace ThePEG;
typedef vector< string > split_vector_type;
/**
* Here is the documentation of the FxFxHandler class.
*
* @see \ref FxFxHandlerInterfaces "The interfaces"
* defined for FxFxHandler.
*/
class FxFxHandler: public QTildeShowerHandler {
/**
* FxFxHandler should have access to our private parts.
*/
friend class FxFxEventHandler;
friend class FxFxReader;
public:
/**
* The default constructor.
*/
FxFxHandler();
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Finalize the object
*/
virtual void dofinish();
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
public:
/**
* Hook to allow vetoing of event after showering hard sub-process
* as in e.g. MLM merging.
*/
virtual bool showerHardProcessVeto() const;
/**
* information for FxFx merging
*/
mutable int npLO_;
mutable int npNLO_;
/**
* information for tree-level merging
*/
mutable vector<double> ptclust_;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/*
* whether a heavy quark has been found in the merging
*/
mutable bool hvqfound = false;
/*
* Run MLM jet-parton matching on the 'extra' jets.
*/
bool lightJetPartonVeto();
/*
* Function that calculates deltaR between a parton and a jet
*/
double partonJetDeltaR(ThePEG::tPPtr partonptr, LorentzMomentum jetmom) const;
/*
* Function that calculates deltaR between two jets
*/
double partonJetDeltaR(LorentzMomentum jetmom1, LorentzMomentum jetmom2) const;
/**
* Find jets using the FastJet package on particlesToCluster_.
*/
void getFastJets(double rjet, Energy ejcut, double etajcut) const;
/**
* Find jets using the FastJet package on particlesToCluster_.
*/
void getFastJetsToMatch(double rjet, Energy ejcut, double etajcut) const;
/**
* Deletes particles from partonsToMatch_ and particlesToCluster_
* vectors so that these contain only the partons to match to the
* jets and the particles used to build jets respectively. By and
* large the candidates for deletion are: vector bosons and their
* decay products, Higgs bosons, photons as well as _primary_, i.e.
* present in the lowest multiplicity process, heavy quarks and
* any related decay products.
*/
void caldel_m() const;
/**
* Deletes particles from partonsToMatch_ and particlesToCluster_
* vectors so that these contain only the partons to match to the
* jets and the particles used to build jets respectively. The candidates
* are chosen according to the information passed from madgraph.
*/
void caldel_mg() const;
/**
* c++ translation of subroutine of same name from alpsho.f.
* Label all particles with status between ISTLO and ISTHI
* (until a particle with status ISTOP is found) as final-state,
* call calsim_m and then put labels back to normal. This
* version keeps only all IST=1 particles rejected by caldel as
* daughters of vetoed heavy-quark mothers: jets complementary
* to those reconstructed by caldel.
*/
void caldel_hvq() const;
/**
* get the MG5_aMC information required for FxFx merging
*/
void getnpFxFx() const;
/**
* get the MG5_aMC information required for FxFx merging
*/
void getECOM() const;
/**
* get the MG5_aMC information required for tree-level merging
*/
void getptclust() const;
/**
* Erases all occurences of a substring from a string
*/
void erase_substr(std::string& subject, const std::string& search) const;
/**
* Get the particles from lastXCombPtr filling the pair
* preshowerISPs_ and particle pointer vector preshowerFSPs_.
*/
void getPreshowerParticles() const;
/**
* Get the particles from eventHandler()->currentEvent()->...
* filling the particle pairs showeredISHs_, showeredISPs_,
* showeredRems_ and the particle pointer vector showeredFSPs_.
*/
void getShoweredParticles() const;
/**
* Allows printing of debug output and sanity checks like
* total momentum consrvation to be carried out.
* debugLevel = -1, 0, ...5
* = no debugging, minimal debugging, ... verbose.
*/
void doSanityChecks(int debugLevel) const;
/**
* Given a pointer to a particle this finds all its final state
* descendents.
*/
void getDescendents(PPtr theParticle) const;
/**
* Accumulates all descendents of tops down to the b and W
* but not including them.
*/
void getTopRadiation(PPtr theParticle) const;
/**
* Sorts a given vector of particles by descending pT or ETJET
*/
ParticleVector pTsort(ParticleVector unsortedVec);
pair< vector<Energy>, vector<Lorentz5Momentum> > ETsort(vector<Energy> unsortedetjet, vector<Lorentz5Momentum> unsortedVec);
/*
* A function that prints a vector of Lorentz5Momenta in a fancy way
*/
void printMomVec(vector<Lorentz5Momentum> momVec);
/*
* A probability function for varying etclus_ about the mean value
*/
Energy etclusran_(double petc) const;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FxFxHandler> initFxFxHandler;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FxFxHandler & operator=(const FxFxHandler &);
private:
/**
* Initial-state incoming partons prior to showering
* (i.e. from lastXCombPtr).
*/
mutable PPair preshowerISPs_;
/**
* Final-state outgoing partICLEs prior to showering
* (i.e. from lastXCombPtr).
*/
mutable ParticleVector preshowerFSPs_;
/**
* Final-state outgoing partICLEs prior to showering _to_be_removed_
* from preShowerFSPs_ prior to the light-parton-light-jet matching
* step. This same list is the starting point for determining
* partonsToMatch_ for the case of merging in heavy quark production.
*/
mutable ParticleVector preshowerFSPsToDelete_;
/**
* Initial-state incoming hadrons after shower of hard process
* (eventHandler()->currentEvent()->incoming()).
*/
mutable PPair showeredISHs_;
/**
* Initial-state incoming partons after shower of hard process
* (look for partonic children of showeredISHs_).
*/
mutable PPair showeredISPs_;
/**
* Final-state outgoing partICLEs after shower of hard process
* (eventHandler()->currentEvent()->getFinalState()).
*/
mutable tPVector showeredFSPs_;
/**
* Final-state outgoing partICLEs after shower of hard process
* _to_be_removed_ from showeredFSPs_ prior to the
* light-parton-light-jet matching step. This same list is the
* starting point for determining particlesToCluster_ for the
* case of merging in heavy quark production.
*/
mutable ParticleVector showeredFSPsToDelete_;
/**
* ONLY the final-state partons from preshowerFSPs_ that are
* supposed to enter the jet-parton matching.
*/
mutable ParticleVector partonsToMatch_;
/*
* The shower progenitors
*/
mutable PPtr theProgenitor;
mutable PPtr theLastProgenitor;
/**
* ONLY the final-state particles from showeredFSPs_ (and maybe
* also showeredRems_) that are supposed to go for jet clustering.
*/
mutable tPVector particlesToCluster_;
/**
* Final-state remnants after shower of hard process
* (look for remnants initially in showeredFSPs_).
*/
mutable PPair showeredRems_;
/**
* the COM of the incoming hadrons
*/
mutable double ECOM_;
/**
* Pointer to the object calculating the strong coupling
*/
ShowerAlphaPtr alphaS_;
/**
* Information extracted from the XComb object
*/
//@{
/**
* The fixed factorization scale used in the MEs.
*/
Energy pdfScale_;
/**
* Centre of mass energy
*/
Energy2 sHat_;
/**
* Constant alphaS used to generate LH events - if not already
* using CKKW scale (ickkw = 1 in AlpGen for example).
*/
double alphaSME_;
//@}
/*
* Number of rapidity segments of the calorimeter.
*/
unsigned int ncy_;
/*
* Number of phi segments of the calorimeter.
*/
unsigned int ncphi_;
/*
* Heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t).
*/
int ihvy_;
/*
* Number of photons in the AlpGen process.
*/
int nph_;
/*
* Number of higgses in the AlpGen process.
*/
int nh_;
/*
* Jet ET cut to apply in jet clustering (in merging).
*/
mutable Energy etclus_;
/*
* Mean Jet ET cut to apply in jet clustering (in merging).
*/
Energy etclusmean_;
/*
* maximum deviation from mean Jet ET cut to apply in jet clustering (in merging).
*/
Energy epsetclus_;
/*
* Cone size used in jet clustering (in merging).
*/
double rclus_;
/*
* Max |eta| for jets in clustering (in merging).
*/
double etaclmax_;
/*
* Default 1.5 factor used to decide if a jet matches a parton
* in merging: if DR(parton,jet)<rclusfactor*rclus the parton
* and jet are said to have been matched.
*/
double rclusfactor_;
/*
* Determines whether to detect the hard process or to manually determine which particles
* to include in the merging. If False, then the ihrd code below is used.
*/
bool hpdetect_;
/*
* The AlpGen hard process code. Relation to the AlpGen process names:
* 1: wqq, 2: zqq, 3: wjet, 4: zjet, 5: vbjet, 6: 2Q, 8: QQh, 9: Njet,
* 10: wcjet, 11: phjet, 12: hjet, 13: top, 14: wphjet, 15: wphqq,
* 16: 2Qph.
*/
int ihrd_;
/*
* The number of light jets in the AlpGen process (i.e. the 'extra' ones).
*/
int njets_;
/*
* Mimimum parton-parton R-sep used for generation (used for hvq merging).
*/
double drjmin_;
/*
* This flags that the highest multiplicity ME-level process is
* being processed.
*/
mutable bool highestMultiplicity_;
/*
* This flags whether the etclus_ (merging scale) should be fixed or variable according to a prob. distribution around the mean
*/
bool etclusfixed_;
/*
* The forwards rapidity span of the calorimeter.
*/
double ycmax_;
/*
* The backwards rapidity span of the calorimeter.
*/
double ycmin_;
/*
* The jet algorithm used for parton-jet matching in the MLM procedure.
*/
int jetAlgorithm_;
/*
* The merging mode (FxFx vs tree-level) used.
*/
int mergemode_;
/*
* Allows the vetoing to be turned off completely - just for convenience.
*/
bool vetoIsTurnedOff_;
/*
* Allows the vetoing on heavy quark decay products to be turned off.
*/
bool vetoHeavyQ_;
/*
* Allows vetoing of heavy flavour
*/
bool vetoHeavyFlavour_;
/*
* Veto if there exist softer unmatched jets than matched
*/
bool vetoSoftThanMatched_;
/*
* Cosine of phi values of calorimeter cell centres.
* Goes phi~=0 to phi~=2*pi (index = 0 ---> ncphi).
* ==> Cosine goes from +1 ---> +1 (index = 0 ---> ncphi).
*/
vector<double> cphcal_;
/*
* Sine of phi values of calorimeter cell centres.
* Goes phi~=0 to phi~=2*pi (index = 0 ---> ncphi).
* ==> Sine goes 0 -> 1 -> 0 -> -1 -> 0 (index = 0 ---> ncphi).
*/
vector<double> sphcal_;
/*
* Cosine of theta values of calorimeter cell centres in Y.
* Goes bwds th~=pi to fwds th~=0 (index = 0 ---> ncy).
* ==> Cosine goes from -1 ---> +1 (index = 0 ---> ncy).
*/
vector<double> cthcal_;
/*
* Sine of theta values of calorimeter cell centres in Y.
* Goes bwds th~=pi to fwds th~=0 (index = 0 ---> ncy).
* ==> Sine goes from 0 ---> +1 ---> 0 (index = 0 ---> ncy).
*/
vector<double> sthcal_;
/*
* Transverse energy deposit in a given calorimeter cell.
* First array index corresponds to rapidity index of cell,
* second array index corresponds to phi cell index.
*/
vector<vector<Energy> > et_;
/*
* For a given calorimeter cell this holds the index of the jet
* that the cell was clustered into.
*/
vector<vector<int> > jetIdx_;
/*
* Vector holding the Lorentz 5 momenta of each jet.
*/
mutable vector<Lorentz5Momentum> pjet_;
/*
* Vector holding the Lorentz 5 momenta of each jet from ME partons
*/
mutable vector<Lorentz5Momentum> pjetME_;
/*
* Vector holding the list of FS particles resulting from
* the particle input to getDescendents.
*/
mutable ParticleVector tmpList_;
/*
* Variables for the C++ translation of the calini_m(), calsim_m(),
* getjet_m(...) and caldel_m() functions
*/
mutable vector<Energy> etjet_;
vector<Energy> etjetME_;
mutable double dely_, delphi_;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FxFxHandler. */
template <>
struct BaseClassTrait<Herwig::FxFxHandler,1> {
/** Typedef of the first base class of FxFxHandler. */
typedef Herwig::QTildeShowerHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FxFxHandler class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::FxFxHandler>
: public ClassTraitsBase<Herwig::FxFxHandler> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FxFxHandler"; }
/**
* The name of a file containing the dynamic library where the class
* FxFxHandler is implemented. It may also include several, space-separated,
* libraries if the class FxFxHandler depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "FxFxHandler.so"; }
};
/** @endcond */
}
#endif /* HERWIG_FxFxHandler_H */
diff --git a/MatrixElement/Hadron/MEDiffraction.h b/MatrixElement/Hadron/MEDiffraction.h
--- a/MatrixElement/Hadron/MEDiffraction.h
+++ b/MatrixElement/Hadron/MEDiffraction.h
@@ -1,303 +1,302 @@
// -*- C++ -*-
#ifndef HERWIG_MEDiffraction_H
#define HERWIG_MEDiffraction_H
//
// This is the declaration of the MEDiffraction class.
//
#include "Herwig/MatrixElement/HwMEBase.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEDiffraction class provides a simple colour singlet exchange matrix element
* to be used in the soft component of the multiple scattering model of the
* underlying event
*
* @see \ref MEDiffractionInterfaces "The interfaces"
* defined for MEDiffraction.
*/
class MEDiffraction: public HwMEBase {
public:
MEDiffraction();
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
/**
* Return the scale associated with the last set phase space point.
*/
virtual Energy2 scale() const;
/**
* Set the typed and momenta of the incoming and outgoing partons to
* be used in subsequent calls to me() and colourGeometries()
* according to the associated XComb object. If the function is
* overridden in a sub class the new function must call the base
* class one first.
*/
virtual void setKinematics();
/**
* The number of internal degrees of freedom used in the matrix
* element.
*/
virtual int nDim() const;
/**
* Generate internal degrees of freedom given nDim() uniform
* random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
* generator, the dSigHatDR should be a smooth function of these
* numbers, although this is not strictly necessary.
* @param r a pointer to the first of nDim() consecutive random numbers.
* @return true if the generation succeeded, otherwise false.
*/
virtual bool generateKinematics(const double * r);
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR() const;
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Get diagram selector. With the information previously supplied with the
* setKinematics method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
* @param dv the diagrams to be weighted.
* @return a Selector relating the given diagrams to their weights.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
* @param diag the diagram chosen.
* @return the possible colour geometries weighted by their
* relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
//@}
/**
* Expect the incoming partons in the laboratory frame
*/
/* virtual bool wantCMS() const { return false; } */
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before a run begins.
*/
virtual void doinitrun();
//@}
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/* The matrix element squared */
double theme2;
/* Use only delta as excited state */
bool deltaOnly;
/* Direction of the excited proton */
unsigned int diffDirection;
/* Number of clusters the dissociated proton decays into */
unsigned int dissociationDecay;
/* The mass of the consitutent quark */
Energy mq() const {return Energy(0.325*GeV);}
/* The mass of the constituent diquark */
Energy mqq() const {return Energy(0.650*GeV);}
/* The proton-pomeron slope */
double theprotonPomeronSlope;
/* The soft pomeron intercept */
double thesoftPomeronIntercept;
/* The soft pomeron slope */
double thesoftPomeronSlope;
/**
* Sample the diffractive mass squared M2 and the momentum transfer t
*/
pair<pair<Energy2,Energy2>,Energy2> diffractiveMassAndMomentumTransfer() const;
/**
* Random value for the diffractive mass squared M2 according to (M2/s0)^(-intercept)
*/
Energy2 randomM2() const;
/**
* Random value for t according to exp(diffSlope*t)
*/
Energy2 randomt(Energy2 M2) const;
/**
* Random value for t according to exp(diffSlope*t) for double diffraction
*/
Energy2 doublediffrandomt(Energy2 M12, Energy2 M22) const;
/**
* Returns the momenta of the two-body decay of momentum pp
*/
pair<Lorentz5Momentum,Lorentz5Momentum> twoBodyDecayMomenta(Lorentz5Momentum pp) const;
/**
* Returns the proton-pomeron slope
*/
InvEnergy2 protonPomeronSlope() const;
/**
* Returns the soft pomeron intercept
*/
double softPomeronIntercept() const;
//M12 and M22 are masses squared of
//outgoing particles
/**
* Returns the minimal possible value of momentum transfer t given the center
* of mass energy and diffractive masses
*/
Energy2 tminfun(Energy2 s, Energy2 M12, Energy2 M22) const;
/**
* Returns the maximal possible value of momentum transfer t given the center
* of mass energy and diffractive masses
*/
Energy2 tmaxfun(Energy2 s , Energy2 M12, Energy2 M22) const;
/**
* Returns the minimal possible value of diffractive mass
*/
//lowest possible mass given the constituent masses of quark and diquark
Energy2 M2min() const{return sqr(getParticleData(2212)->mass()+mq()+mqq());}
/**
* Returns the maximal possible value of diffractive mass
*/
Energy2 M2max() const{
return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass());
}//TODO:modify to get proper parameters
InvEnergy2 softPomeronSlope() const;
/* Kallen function */
- template<int L, int E, int Q, int DL, int DE, int DQ>
- Qty<2*L,2*E,2*Q,DL,DE,DQ> kallen(Qty<L,E,Q,DL,DE,DQ> a,
- Qty<L,E,Q,DL,DE,DQ> b,
- Qty<L,E,Q,DL,DE,DQ> c) const {
+ template<typename A, typename B, typename C>
+ auto kallen(A a, B b, C c) const -> decltype(a*a)
+ {
return a*a + b*b + c*c - 2.0*(a*b + b*c + c*a);
}
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEDiffraction & operator=(const MEDiffraction &);
bool isInRunPhase;
/* The proton mass */
Energy theProtonMass;
};
}
#endif /* HERWIG_MEDiffraction_H */
diff --git a/MatrixElement/Matchbox/Utility/SpinorHelicity.h b/MatrixElement/Matchbox/Utility/SpinorHelicity.h
--- a/MatrixElement/Matchbox/Utility/SpinorHelicity.h
+++ b/MatrixElement/Matchbox/Utility/SpinorHelicity.h
@@ -1,434 +1,435 @@
// -*- C++ -*-
//
// SpinorHelicity.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SpinorHelicity_H
#define HERWIG_SpinorHelicity_H
#include "ThePEG/Config/Complex.h"
#include "ThePEG/Vectors/LorentzVector.h"
#include <boost/operators.hpp>
namespace Herwig {
using namespace ThePEG;
namespace SpinorHelicity {
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Tag for |p>
*
*/
struct PlusSpinorTag {};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Tag for |p]
*
*/
struct MinusSpinorTag {};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Tag for <p|
*
*/
struct PlusConjugateSpinorTag {};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Tag for [p|
*
*/
struct MinusConjugateSpinorTag {};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Helpers for commonly encountered types.
*
*/
template<class Value>
struct SpinorMultiplicationTraits {
- typedef typename BinaryOpTraits<Value,Value>::MulT ResultType;
+ typedef decltype(sqr(std::declval<Value>())) ResultType;
typedef complex<ResultType> ComplexResultType;
typedef LorentzVector<ComplexResultType> ComplexVectorResultType;
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Helpers for Weyl spinors
*
*/
template<class Type>
struct WeylSpinorTraits;
// specialize for |p>
template<>
struct WeylSpinorTraits<PlusSpinorTag> {
template<class Value, class MValue>
static pair<complex<Value>,complex<Value> >
components(const LorentzVector<MValue>& p) {
if ( p.t() < ZERO ) {
pair<complex<Value>,complex<Value> > res =
components<Value,MValue>(-p);
// do not revert to *=, breaks with XCode 5.1
res.first = res.first * Complex(0.,1.);
res.second = res.second * Complex(0.,1.);
return res;
}
Energy pPlus = p.t() + p.x();
if ( abs(pPlus) < 1.e-10 * GeV ) {
return make_pair(complex<Value>(ZERO),
complex<Value>(sqrt(2.*p.t())));
}
return make_pair(complex<Value>(sqrt(pPlus)),
complex<Value>(p.z()/sqrt(pPlus),p.y()/sqrt(pPlus)));
}
};
// specialize for |p]
template<>
struct WeylSpinorTraits<MinusSpinorTag> {
template<class Value, class MValue>
static pair<complex<Value>,complex<Value> >
components(const LorentzVector<MValue>& p) {
if ( p.t() < ZERO ) {
pair<complex<Value>,complex<Value> > res =
components<Value,MValue>(-p);
// do not revert to *=, breaks with XCode 5.1
res.first = res.first * Complex(0.,1.);
res.second = res.second * Complex(0.,1.);
return res;
}
Energy pPlus = p.t() + p.x();
if ( abs(pPlus) < 1.e-10 * GeV ) {
return make_pair(complex<Value>(sqrt(2.*p.t())),
complex<Value>(ZERO));
}
return make_pair(complex<Value>(p.z()/sqrt(pPlus),-p.y()/sqrt(pPlus)),
-complex<Value>(sqrt(pPlus)));
}
};
// specialize for <p|
template<>
struct WeylSpinorTraits<PlusConjugateSpinorTag> {
typedef PlusSpinorTag ConjugateSpinorTag;
typedef MinusSpinorTag BarSpinorTag;
template<class Value, class MValue>
static pair<complex<Value>,complex<Value> >
components(const LorentzVector<MValue>& p) {
pair<complex<Value>,complex<Value> > res =
WeylSpinorTraits<PlusSpinorTag>::template components<Value>(p);
res.first = -res.first;
swap(res.first,res.second);
return res;
}
};
// specialize for [p|
template<>
struct WeylSpinorTraits<MinusConjugateSpinorTag> {
typedef MinusSpinorTag ConjugateSpinorTag;
typedef PlusSpinorTag BarSpinorTag;
template<class Value, class MValue>
static pair<complex<Value>,complex<Value> >
components(const LorentzVector<MValue>& p) {
pair<complex<Value>,complex<Value> > res =
WeylSpinorTraits<MinusSpinorTag>::template components<Value>(p);
res.second = -res.second;
swap(res.first,res.second);
return res;
}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Base class for Weyl spinors
*
*/
template<class Type, class Value>
class WeylSpinor {
public:
typedef complex<Value> ComplexType;
typedef pair<ComplexType,ComplexType> ComponentsType;
typedef Type Tag;
typedef WeylSpinorTraits<Tag> Traits;
typedef Value ValueType;
private:
/**
* The components
*/
ComponentsType theComponents;
public:
/**
* Construct from components
*/
explicit WeylSpinor(const ComponentsType& c = ComponentsType())
: theComponents(c) {}
/**
* Construct from momentum
*/
template<class MValue>
explicit WeylSpinor(const LorentzVector<MValue>& p)
: theComponents(Traits::template components<Value>(p)) {}
/**
* Return the components.
*/
const ComponentsType& components() const { return theComponents; }
/**
* Return the first component
*/
const ComplexType& s1() const { return theComponents.first; }
/**
* Return the second component
*/
const ComplexType& s2() const { return theComponents.second; }
};
/** Define |p> */
typedef WeylSpinor<PlusSpinorTag,SqrtEnergy> PlusSpinor;
/** Define |p] */
typedef WeylSpinor<MinusSpinorTag,SqrtEnergy> MinusSpinor;
/** Define <p| */
typedef WeylSpinor<PlusConjugateSpinorTag,SqrtEnergy> PlusConjugateSpinor;
/** Define [p| */
typedef WeylSpinor<MinusConjugateSpinorTag,SqrtEnergy> MinusConjugateSpinor;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Weyl spinor product
*
*/
template<class Type, class Value>
class SpinorProduct
: public boost::addable<SpinorProduct<Type,Value> >,
public boost::subtractable<SpinorProduct<Type,Value> >,
public boost::multipliable<SpinorProduct<Type,Value>, double>,
public boost::multipliable<SpinorProduct<Type,Value>, complex<double> > {
public:
typedef typename SpinorMultiplicationTraits<Value>::ComplexResultType ResultType;
typedef WeylSpinor<Type,Value> LeftSpinorType;
typedef typename WeylSpinorTraits<Type>::ConjugateSpinorTag RightSpinorTag;
typedef WeylSpinor<RightSpinorTag,Value> RightSpinorType;
private:
/**
* The result.
*/
ResultType theResult;
public:
/**
* Construct from two spinors; note that the
* spinor metric is included, when constructing spinors.
* Typedefs break zero products like <p|q]
*/
explicit SpinorProduct(const LeftSpinorType& left,
const RightSpinorType& right)
: theResult(left.s1()*right.s1()+left.s2()*right.s2()) {}
/**
* Implicitly convert to complex value
*/
operator ResultType() const { return theResult; }
/**
* Return result
*/
ResultType eval() const { return theResult; }
public:
SpinorProduct& operator+= (const SpinorProduct& other) {
theResult += other.theResult;
return *this;
}
SpinorProduct& operator-= (const SpinorProduct& other) {
theResult -= other.theResult;
return *this;
}
SpinorProduct& operator*= (double x) {
theResult *= x;
return *this;
}
SpinorProduct& operator*= (complex<double> x) {
theResult *= x;
return *this;
}
};
/** Define <pq> */
typedef SpinorProduct<PlusConjugateSpinorTag,SqrtEnergy> PlusSpinorProduct;
/** Define [pq] */
typedef SpinorProduct<MinusConjugateSpinorTag,SqrtEnergy> MinusSpinorProduct;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Weyl spinor current.
*
*/
template<class Type, class Value>
class SpinorCurrent
: public boost::addable<SpinorCurrent<Type,Value> >,
public boost::subtractable<SpinorCurrent<Type,Value> >,
public boost::multipliable<SpinorCurrent<Type,Value>, double>,
public boost::multipliable<SpinorCurrent<Type,Value>, complex<double> > {
public:
typedef typename SpinorMultiplicationTraits<Value>::ComplexVectorResultType ResultType;
typedef WeylSpinor<Type,Value> LeftSpinorType;
typedef typename WeylSpinorTraits<Type>::BarSpinorTag RightSpinorTag;
typedef WeylSpinor<RightSpinorTag,Value> RightSpinorType;
private:
ResultType theResult;
/**
* Calculate [p|\gamma^\mu|q>
*/
ResultType evaluate(const WeylSpinor<MinusConjugateSpinorTag,Value>& left,
const WeylSpinor<PlusSpinorTag,Value>& right) {
return
ResultType(right.s1()*left.s1()-right.s2()*left.s2(),
complex<double>(0.,1.)*(right.s1()*left.s2()-right.s2()*left.s1()),
right.s1()*left.s2()+right.s2()*left.s1(),
right.s1()*left.s1()+right.s2()*left.s2());
}
/**
* Calculate <p|\gamma^\mu|q]
*/
ResultType evaluate(const WeylSpinor<PlusConjugateSpinorTag,Value>& left,
const WeylSpinor<MinusSpinorTag,Value>& right) {
return
ResultType(-right.s1()*left.s1()+right.s2()*left.s2(),
-complex<double>(0.,1.)*(right.s1()*left.s2()-right.s2()*left.s1()),
-right.s1()*left.s2()-right.s2()*left.s1(),
right.s1()*left.s1()+right.s2()*left.s2());
}
public:
/**
* Construct from two spinors.
* Typedefs break zero products like <p|\gamma^\mu|q>
*/
explicit SpinorCurrent(const LeftSpinorType& left,
const RightSpinorType& right)
: theResult(evaluate(left,right)) {}
/**
* Implicitly convert to complex Lorentz vector
*/
operator ResultType() const { return theResult; }
/**
* Return result
*/
ResultType eval() const { return theResult; }
public:
SpinorCurrent& operator+= (const SpinorCurrent& other) {
theResult += other.theResult;
return *this;
}
SpinorCurrent& operator-= (const SpinorCurrent& other) {
theResult -= other.theResult;
return *this;
}
SpinorCurrent& operator*= (double x) {
theResult *= x;
return *this;
}
SpinorCurrent& operator*= (complex<double> x) {
theResult *= x;
return *this;
}
};
/** Define <p|\gamma^\mu|q] */
typedef SpinorCurrent<PlusConjugateSpinorTag,SqrtEnergy> PlusSpinorCurrent;
/** Define [p|\gamma^\mu|q> */
typedef SpinorCurrent<MinusConjugateSpinorTag,SqrtEnergy> MinusSpinorCurrent;
/**
* Return |c|^2
*/
template<class T>
- typename BinaryOpTraits<T,T>::MulT abs2(const complex<T>& x) {
+ auto abs2(const complex<T>& x) -> decltype((x*conj(x)).real())
+ {
return (x*conj(x)).real();
}
}
}
#endif // HERWIG_SpinorHelicity_H
diff --git a/MatrixElement/Reweighters/ReweightEW.h b/MatrixElement/Reweighters/ReweightEW.h
--- a/MatrixElement/Reweighters/ReweightEW.h
+++ b/MatrixElement/Reweighters/ReweightEW.h
@@ -1,169 +1,164 @@
// -*- C++ -*-
#ifndef Herwig_ReweightEW_H
#define Herwig_ReweightEW_H
//
// This is the declaration of the ReweightEW class.
//
#include "ThePEG/MatrixElement/ReweightBase.h"
#include <array>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the ReweightEW class.
*
* @see \ref ReweightEWInterfaces "The interfaces"
* defined for ReweightEW.
*/
class ReweightEW: public ReweightBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ReweightEW();
/**
* The destructor.
*/
virtual ~ReweightEW();
//@}
public:
/**
* Return the weight for the kinematical configuation provided by
* the assigned XComb object (in the LastXCombInfo base class).
*/
virtual double weight() const;
/**
* Return values of the last evaluation (double/doubles in GeV2)
*/
double lastS() const {return thelasts;}
double lastT() const {return thelastt;}
double lastK() const {return thelastk;}
void setSTK(double s, double t, double K);
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The last s
*/
mutable double thelasts;
/**
* The last t
*/
mutable double thelastt;
/**
* The last K-factor
*/
mutable double thelastk;
/**
* The table of K factors to be read from file
*/
std::array<std::array<double,6>,40001> tab;
/**
* EW K factor filename
*/
string filename;
public:
/**
* Computation of K factors from table (s and t in GeV)
*/
double EWKFac(unsigned int f, double s, double t) const;
private:
/**
* initialize tables
*/
void inittable();
private:
/**
-
- /** @name Standard Interfaced functions. */
- //@{
-
- /**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ReweightEW & operator=(const ReweightEW &);
};
}
#endif /* Herwig_ReweightEW_H */
diff --git a/Models/Feynrules/python/ufo2peg/particles.py b/Models/Feynrules/python/ufo2peg/particles.py
--- a/Models/Feynrules/python/ufo2peg/particles.py
+++ b/Models/Feynrules/python/ufo2peg/particles.py
@@ -1,206 +1,243 @@
from string import Template
# ignore these, they're in Hw++ already # TODO reset Hw++ settings instead
SMPARTICLES = {
1:'d',
2:'u',
3:'s',
4:'c',
5:'b',
6:'t', # think about this one later
11:'e-',
12:'nu_e',
13:'mu-',
14:'nu_mu',
15:'tau-',
16:'nu_tau',
21:'g',
22:'gamma',
23:'Z0',
24:'W+',
-1:'dbar',
-2:'ubar',
-3:'sbar',
-4:'cbar',
-5:'bbar',
-6:'tbar',
-11:'e+',
-12:'nu_ebar',
-13:'mu+',
-14:'nu_mubar',
-15:'tau+',
-16:'nu_taubar',
-24:'W-',
}
particleT = Template(
"""
create ThePEG::ParticleData $name
# values set to 999999 are recalculated later from other model parameters
setup $name $pdg_code $name $mass $width $wcut $ctau $charge $color $spin 0
"""
)
class ParticleConverter:
'Convert a FR particle to extract the information ThePEG needs.'
def __init__(self,p,parmsubs,modelparameters):
self.name = p.name
self.pdg_code = p.pdg_code
self.spin = p.spin
self.color = p.color
if self.color == 1:
self.color = 0
self.selfconjugate = 0
self.mass = parmsubs[str(p.mass)]
if type(self.mass) == str:
value = modelparameters[self.mass]
try:
value = value.real
except:
pass
newname = '%s_ABS' % self.mass
self.mass = '${%s}' % newname
modelparameters[newname] = abs(value)
else:
try:
self.mass = self.mass.real
except:
pass
self.mass = 999999. # abs(self.mass)
hbarc = 197.3269631e-15 # GeV mm (I hope ;-) )
self.width = parmsubs[str(p.width)]
if type(self.width) == str:
width = modelparameters[self.width]
ctau = (hbarc / width) if width != 0 else 0
newname = '%s_CTAU' % self.width
self.ctau = '${%s}' % newname
modelparameters[newname] = ctau
wcut = 10 * width
newname = '%s_WCUT' % self.width
self.wcut = '${%s}' % newname
modelparameters[newname] = wcut
self.width = '${%s}' % self.width
else:
self.ctau = 999999. # (hbarc / self.width) if self.width != 0 else 0
self.wcut = 999999. #10.0 * self.width
self.width = 999999. # was blank line before
self.charge = int(3 * p.charge)
def subs(self):
return self.__dict__
-
+def check_effective_vertex(FR,p,ig) :
+ for vertex in FR.all_vertices:
+ if(len(vertex.particles) != 3) : continue
+ if(p not in vertex.particles ) : continue
+ ng=0
+ for part in vertex.particles :
+ if(part.pdg_code==ig) : ng+=1
+ if(ng==2) :
+ return False
+ return True
def thepeg_particles(FR,parameters,modelname,modelparameters,forbidden_names):
plist = []
antis = {}
names = []
splittings = []
- done_splitting=[]
+ done_splitting_QCD=[]
+ done_splitting_QED=[]
for p in FR.all_particles:
if p.spin == -1:
continue
gsnames = ['goldstone',
'goldstoneboson',
'GoldstoneBoson']
def gstest(name):
try:
return getattr(p,name)
except AttributeError:
return False
if any(map(gstest, gsnames)):
continue
if p.pdg_code in SMPARTICLES:
continue
if p.pdg_code == 25:
plist.append(
"""
set /Herwig/Particles/h0:Mass_generator NULL
set /Herwig/Particles/h0:Width_generator NULL
rm /Herwig/Masses/HiggsMass
rm /Herwig/Widths/hWidth
"""
)
if p.name in forbidden_names:
print 'RENAMING PARTICLE',p.name,'as ',p.name+'_UFO'
p.name +="_UFO"
subs = ParticleConverter(p,parameters,modelparameters).subs()
plist.append( particleT.substitute(subs) )
pdg, name = subs['pdg_code'], subs['name']
names.append(name)
if -pdg in antis:
plist.append( 'makeanti %s %s\n' % (antis[-pdg], name) )
else:
plist.append( 'insert /Herwig/NewPhysics/NewModel:DecayParticles 0 %s\n' % name )
plist.append( 'insert /Herwig/Shower/ShowerHandler:DecayInShower 0 %s # %s' % (abs(pdg), name) )
antis[pdg] = name
selfconjugate = 1
class SkipMe(Exception):
pass
def spin_name(s):
spins = { 1 : 'Zero',
2 : 'Half',
3 : 'One' }
if s not in spins:
raise SkipMe()
else:
return spins[s]
def col_name(c):
cols = { 3 : 'Triplet',
6 : 'Sextet',
8 : 'Octet' }
return cols[c]
try:
- if p.color in [3,6,8] and abs(pdg) not in done_splitting: # which colors?
- done_splitting.append(abs(pdg))
- splitname = '{name}SplitFn'.format(name=p.name)
- sudname = '{name}Sudakov'.format(name=p.name)
+ # QCD splitting functions
+ if p.color in [3,6,8] and abs(pdg) not in done_splitting_QCD: # which colors?
+ done_splitting_QCD.append(abs(pdg))
+ splitname = '{name}SplitFnQCD'.format(name=p.name)
+ sudname = '{name}SudakovQCD'.format(name=p.name)
splittings.append(
"""
create Herwig::{s}{s}OneSplitFn {name}
set {name}:InteractionType QCD
set {name}:ColourStructure {c}{c}Octet
cp /Herwig/Shower/SudakovCommon {sudname}
set {sudname}:SplittingFunction {name}
do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {pname}->{pname},g; {sudname}
""".format(s=spin_name(p.spin), name=splitname,
c=col_name(p.color), pname=p.name, sudname=sudname)
)
except SkipMe:
pass
+ # QED splitting functions
+ try:
+ if p.charge != 0 and abs(pdg) not in done_splitting_QED:
+ done_splitting_QED.append(abs(pdg))
+ splitname = '{name}SplitFnQED'.format(name=p.name)
+ sudname = '{name}SudakovQED'.format(name=p.name)
+ splittings.append(
+"""
+create Herwig::{s}{s}OneSplitFn {name}
+set {name}:InteractionType QED
+set {name}:ColourStructure ChargedChargedNeutral
+cp /Herwig/Shower/SudakovCommon {sudname}
+set {sudname}:SplittingFunction {name}
+set {sudname}:Alpha /Herwig/Shower/AlphaQED
+do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {pname}->{pname},gamma; {sudname}
+""".format(s=spin_name(p.spin), name=splitname, pname=p.name, sudname=sudname)
+ )
+ except SkipMe:
+ pass
if p.charge == 0 and p.color == 1 and p.spin == 1:
- plist.append(
+ if(check_effective_vertex(FR,p,21)) :
+ plist.append(
+"""
+insert /Herwig/{ModelName}/V_GenericHGG:Bosons 0 {pname}
+""".format(pname=p.name, ModelName=modelname)
+ )
+ if(check_effective_vertex(FR,p,22)) :
+ plist.append(
"""
insert /Herwig/{ModelName}/V_GenericHPP:Bosons 0 {pname}
-insert /Herwig/{ModelName}/V_GenericHGG:Bosons 0 {pname}
""".format(pname=p.name, ModelName=modelname)
- )
+ )
+
return ''.join(plist)+''.join(splittings), names
diff --git a/Tests/Makefile.am b/Tests/Makefile.am
--- a/Tests/Makefile.am
+++ b/Tests/Makefile.am
@@ -1,365 +1,365 @@
AM_LDFLAGS += -module -avoid-version -rpath /dummy/path/not/used
EXTRA_DIST = Inputs python Rivet
EXTRA_LTLIBRARIES = LeptonTest.la GammaTest.la HadronTest.la DISTest.la
if WANT_LIBFASTJET
EXTRA_LTLIBRARIES += HadronJetTest.la LeptonJetTest.la
HadronJetTest_la_SOURCES = \
Hadron/VHTest.h Hadron/VHTest.cc\
Hadron/VTest.h Hadron/VTest.cc\
Hadron/HTest.h Hadron/HTest.cc
HadronJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
HadronJetTest_la_LIBADD = $(FASTJETLIBS)
LeptonJetTest_la_SOURCES = \
Lepton/TopDecay.h Lepton/TopDecay.cc
LeptonJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
-I$(FASTJETPATH)
LeptonJetTest_la_LIBADD = $(FASTJETLIBS)
endif
LeptonTest_la_SOURCES = \
Lepton/VVTest.h Lepton/VVTest.cc \
Lepton/VBFTest.h Lepton/VBFTest.cc \
Lepton/VHTest.h Lepton/VHTest.cc \
Lepton/FermionTest.h Lepton/FermionTest.cc
GammaTest_la_SOURCES = \
Gamma/GammaMETest.h Gamma/GammaMETest.cc \
Gamma/GammaPMETest.h Gamma/GammaPMETest.cc
DISTest_la_SOURCES = \
DIS/DISTest.h DIS/DISTest.cc
HadronTest_la_SOURCES = \
Hadron/HadronVVTest.h Hadron/HadronVVTest.cc\
Hadron/HadronVBFTest.h Hadron/HadronVBFTest.cc\
Hadron/WHTest.h Hadron/WHTest.cc\
Hadron/ZHTest.h Hadron/ZHTest.cc\
Hadron/VGammaTest.h Hadron/VGammaTest.cc\
Hadron/ZJetTest.h Hadron/ZJetTest.cc\
Hadron/WJetTest.h Hadron/WJetTest.cc\
Hadron/QQHTest.h Hadron/QQHTest.cc
REPO = $(top_builddir)/src/HerwigDefaults.rpo
HERWIG = $(top_builddir)/src/Herwig
HWREAD = $(HERWIG) read --repo $(REPO) -L $(builddir)/.libs -i $(top_builddir)/src
HWRUN = $(HERWIG) run
tests : tests-LEP tests-DIS tests-LHC tests-Gamma
if WANT_LIBFASTJET
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-Quarks test-LEP-Leptons \
test-LEP-TopDecay
else
tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons
endif
tests-DIS : test-DIS-Charged test-DIS-Neutral
if WANT_LIBFASTJET
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top test-LHC-Bottom \
test-LHC-WHJet test-LHC-ZHJet test-LHC-HJet test-LHC-ZShower test-LHC-WShower\
test-LHC-WHJet-Powheg test-LHC-ZHJet-Powheg test-LHC-HJet-Powheg \
test-LHC-ZShower-Powheg test-LHC-WShower-Powheg
else
tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top
endif
tests-Gamma : test-Gamma-FF test-Gamma-WW test-Gamma-P
if WANT_LIBFASTJET
test-LEP-% : Inputs/LEP-%.in LeptonTest.la LeptonJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LEP-% : Inputs/LEP-%.in LeptonTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
Rivet-LHC-Matchbox-% : Rivet/LHC-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-TVT-Matchbox-% : Rivet/TVT-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-TVT-Dipole-% : Rivet/TVT-Dipole-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-LHC-Dipole-% : Rivet/LHC-Dipole-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet/LEP-%.in :
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/DIS-%.in :
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/BFactory-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/TVT-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/LHC-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/Star-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/SppS-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet/ISR-%.in:
python/make_input_files.py $(notdir $(subst .in,,$@))
Rivet-LEP-Matchbox-% : Rivet/LEP-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-LEP-Dipole-% : Rivet/LEP-Dipole-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-BFactory-Matchbox-% : Rivet/BFactory-Matchbox-%.in
if [ ! -d Rivet-$(notdir $(subst .in,,$<)) ]; then mkdir Rivet-$(notdir $(subst .in,,$<)); fi;
cd Rivet-$(notdir $(subst .in,,$<)); echo `pwd`; \
../$(HERWIG) read --repo ../$(REPO) -L ../$(top_builddir)/lib -i ../$(top_builddir)/src ../$<; \
../$(HERWIG) run $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}; \
mv $(notdir $(subst .in,.yoda,$<)) ..; \
cd ..
Rivet-LEP-% : Rivet/LEP-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-BFactory-% : Rivet/BFactory-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-TVT-% : Rivet/TVT-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-DIS-% : Rivet/DIS-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-LHC-% : Rivet/LHC-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-Star-% : Rivet/Star-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-SppS-% : Rivet/SppS-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-ISR-% : Rivet/ISR-%.in
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
Rivet-inputfiles: $(shell echo Rivet/LEP{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{9.4,12,13,17,27.6,29,30.2,30.7,30.75,30,31.3,34.8,43.6,50,52,55,56,57,60.8,60,61.4,10,12.8,22,26.8,35,44,48.0,91,93.0,130,133,136,161,172,177,183,189,192,196,197,200,202,206,91-nopi}.in) \
$(shell echo Rivet/LEP{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Powheg,-Matchbox-Powheg}-14.in) \
$(shell echo Rivet/LEP{,-Dipole}-{10.5,11.96,12.8,13.96,16.86,21.84,26.8,28.48,35.44,48.0,97.0}-gg.in) \
$(shell echo Rivet/BFactory{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg}-{10.52,10.52-sym,10.54,10.45}.in) \
$(shell echo Rivet/BFactory{,-Dipole}-{Upsilon,Upsilon2,Upsilon4,Tau,10.58-res}.in) \
$(shell echo Rivet/DIS{,-NoME,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{e--LowQ2,e+-LowQ2,e+-HighQ2}.in) \
$(shell echo Rivet/TVT{,-Powheg,-Matchbox,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox-Powheg,-Merging}-{Run-I-Z,Run-I-W,Run-I-WZ,Run-II-Z-e,Run-II-Z-{,LowMass-,HighMass-}mu,Run-II-W}.in) \
$(shell echo Rivet/TVT{,-Dipole}-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet,PromptPhoton}.in) \
$(shell echo Rivet/TVT-Powheg-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.in) \
$(shell echo Rivet/TVT{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{Run-II-Jets-{0..11},Run-I-Jets-{1..8}}.in ) \
$(shell echo Rivet/TVT{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{630-Jets-{1..3},300-Jets-1,900-Jets-1}.in ) \
$(shell echo Rivet/TVT{,-Dipole}-{Run-I,Run-II,300,630,900}-UE.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-DiJets-{1..7}-{A,B,C}.in ) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{7,8,13}-Jets-{0..10}.in ) \
$(shell echo Rivet/LHC{,-Dipole}-{900,2360,2760,7,8,13}-UE.in ) \
$(shell echo Rivet/LHC{,-Dipole}-{900,7,13}-UE-Long.in ) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Charm-{1..5}.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Bottom-{0..8}.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-7-Top-{L,SL}.in) \
$(shell echo Rivet/LHC{,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Matchbox,-Matchbox-Powheg,-Merging}-{7,8,13}-Top-All.in) \
$(shell echo Rivet/Star{,-Dipole}-{UE,Jets-{1..4}}.in ) \
$(shell echo Rivet/SppS{,-Dipole}-{200,500,900,546}-UE.in ) \
$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{W-{e,mu},13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},Z-{e,mu,mu-SOPHTY},Z-LowMass-{e,mu},Z-MedMass-e,WZ,WW-{emu,ll},ZZ-{ll,lv},8-ZZ-lv,8-WW-ll,W-Z-{e,mu}}.in) \
$(shell echo Rivet/LHC{,-Dipole}-7-{W,Z}Gamma-{e,mu}.in) \
$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{7-W-Jet-{1..3}-e,7-Z-Jet-{0..3}-e,7-Z-Jet-0-mu}.in) \
$(shell echo Rivet/LHC{-Matchbox,-Matchbox-Powheg,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{Z-b,Z-bb,W-b,8-Z-jj}.in) \
$(shell echo Rivet/LHC{,-Dipole}-{7,8}-PromptPhoton-{1..4}.in) Rivet/LHC-GammaGamma-7.in \
$(shell echo Rivet/LHC{,-Powheg}-{7,8}-{DiPhoton-GammaGamma,DiPhoton-GammaJet}.in) \
$(shell echo Rivet/LHC{,-Powheg,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-{ggH,VBF,WH,ZH}.in) \
$(shell echo Rivet/LHC{,-Powheg,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-8-{{ggH,VBF,WH,ZH}{,-GammaGamma},ggH-WW}.in) \
$(shell echo Rivet/LHC{,-Matchbox,-Matchbox-Powheg,-Dipole,-Dipole-MCatNLO,-Dipole-Matchbox-Powheg,-Merging}-ggHJet.in)
# $(shell echo Rivet/ISR-{30,44,53,62}-UE.in ) $(shell echo Rivet/SppS-{53,63}-UE.in )
Rivet-LEP: $(shell echo Rivet-LEP{,-Powheg,-Matchbox,-Dipole}-{9.4,12,13,17,27.6,29,30.2,30.7,30.75,30,31.3,34.8,43.6,50,52,55,56,57,60.8,60,61.4,10,12.8,22,26.8,35,44,48.0,91,93.0,130,133,136,161,172,177,183,189,192,196,197,200,202,206,91-nopi}) \
$(shell echo Rivet-LEP-{10.5,11.96,12.8,13.96,16.86,21.84,26.8,28.48,35.44,48.0,97.0}-gg) \
- $(shell echo Rivet-LEP{,-Powheg,-Matchbox-Powheg}-14.in)
+ $(shell echo Rivet-LEP{,-Powheg,-Matchbox-Powheg}-14)
rm -rf Rivet-LEP
python/merge-LEP --with-gg LEP
python/merge-LEP LEP-Powheg
python/merge-LEP LEP-Matchbox
python/merge-LEP LEP-Dipole
rivet-mkhtml -o Rivet-LEP LEP.yoda:Hw LEP-Powheg.yoda:Hw-Powheg LEP-Matchbox.yoda:Hw-Matchbox LEP-Dipole.yoda:Hw-Dipole
Rivet-BFactory: $(shell echo Rivet-BFactory{,-Powheg,-Matchbox,-Dipole}-{10.52,10.52-sym,10.54,10.45}) \
$(shell echo Rivet-BFactory-{Upsilon,Upsilon2,Upsilon4,Tau,10.58-res,10.58})
rm -rf Rivet-BFactory
python/merge-BFactory BFactory
python/merge-BFactory BFactory-Powheg
python/merge-BFactory BFactory-Matchbox
python/merge-BFactory BFactory-Dipole
rivet-mkhtml -o Rivet-BFactory BFactory.yoda:Hw BFactory-Powheg.yoda:Hw-Powheg BFactory-Matchbox.yoda:Hw-Matchbox BFactory-Dipole.yoda:Hw-Dipole
Rivet-DIS: $(shell echo Rivet-DIS{,-NoME,-Powheg,-Matchbox,-Dipole}-{e--LowQ2,e+-LowQ2,e+-HighQ2})
rm -rf Rivet-DIS
python/merge-DIS DIS
python/merge-DIS DIS-Powheg
python/merge-DIS DIS-NoME
python/merge-DIS DIS-Matchbox
python/merge-DIS DIS-Dipole
rivet-mkhtml -o Rivet-DIS DIS.yoda:Hw DIS-Powheg.yoda:Hw-Powheg DIS-NoME.yoda:Hw-NoME DIS-Matchbox.yoda:Hw-Matchbox DIS-Dipole.yoda:Hw-Dipole
Rivet-TVT-WZ: $(shell echo Rivet-TVT{,-Powheg,-Matchbox,-Dipole}-{Run-I-Z,Run-I-W,Run-I-WZ,Run-II-Z-{e,{,LowMass-,HighMass-}mu},Run-II-W})
rm -rf Rivet-TVT-WZ
python/merge-TVT-EW TVT-Run-II-W.yoda TVT-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Run-I-{W,Z,WZ}.yoda -o TVT-WZ.yoda
python/merge-TVT-EW TVT-Powheg-Run-II-W.yoda TVT-Powheg-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Powheg-Run-I-{W,Z,WZ}.yoda -o TVT-Powheg-WZ.yoda
python/merge-TVT-EW TVT-Matchbox-Run-II-W.yoda TVT-Matchbox-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Matchbox-Run-I-{W,Z,WZ}.yoda -o TVT-Matchbox-WZ.yoda
python/merge-TVT-EW TVT-Dipole-Run-II-W.yoda TVT-Dipole-Run-II-Z-{e,{,LowMass-,HighMass-}mu}.yoda\
TVT-Dipole-Run-I-{W,Z,WZ}.yoda -o TVT-Dipole-WZ.yoda
rivet-mkhtml -o Rivet-TVT-WZ TVT-WZ.yoda:Hw TVT-Powheg-WZ.yoda:Hw-Powheg TVT-Matchbox-WZ.yoda:Hw-Matchbox TVT-Dipole-WZ.yoda:Hw-Dipole
Rivet-TVT-Photon: $(shell echo Rivet-TVT-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet,PromptPhoton}) \
$(shell echo Rivet-TVT-Powheg-Run-II-{DiPhoton-GammaGamma,DiPhoton-GammaJet})
rm -rf Rivet-TVT-Photon
python/merge-TVT-Photon TVT -o TVT-Photon.yoda
python/merge-TVT-Photon TVT-Powheg -o TVT-Powheg-Photon.yoda
rivet-mkhtml -o Rivet-TVT-Photon TVT-Photon.yoda:Hw TVT-Powheg-Photon.yoda:Hw-Powheg
Rivet-TVT-Jets: $(shell echo Rivet-TVT-{Run-II-Jets-{0..11},Run-I-Jets-{1..8}} ) \
$(shell echo Rivet-TVT-{630-Jets-{1..3},300-Jets-1,900-Jets-1} ) \
$(shell echo Rivet-TVT-{Run-I,Run-II,300,630,900}-UE)
rm -rf Rivet-TVT-Jets
python/merge-TVT-Jets TVT
rivet-mkhtml -o Rivet-TVT-Jets TVT-Jets.yoda:Hw
#python/merge-TVT-Energy TVT
#rivet-merge-CDF_2012_NOTE10874 TVT-300-Energy.yoda TVT-900-Energy.yoda TVT-1960-Energy.yoda
#flat2yoda RatioPlots.dat -o TVT-RatioPlots.yoda
Rivet-LHC-Jets: $(shell echo Rivet-LHC-7-DiJets-{1..7}-{A,B,C} ) \
$(shell echo Rivet-LHC-{7,8,13}-Jets-{0..10} ) \
$(shell echo Rivet-LHC-{900,2360,2760,7,8,13}-UE ) \
$(shell echo Rivet-LHC-{900,7,13}-UE-Long ) \
$(shell echo Rivet-LHC-7-Charm-{1..5}) \
$(shell echo Rivet-LHC-7-Bottom-{0..8}) \
$(shell echo Rivet-LHC-7-Top-{L,SL})\
$(shell echo Rivet-LHC-{7,8,13}-Top-All)
rm -rf Rivet-LHC-Jets
python/merge-LHC-Jets LHC
rivet-mkhtml -o Rivet-LHC-Jets LHC-Jets.yoda:Hw
Rivet-Star: $(shell echo Rivet-Star-{UE,Jets-{1..4}} )
rm -rf Rivet-Star
python/merge-Star Star
rivet-mkhtml -o Rivet-Star Star.yoda
Rivet-SppS: $(shell echo Rivet-ISR-{30,44,53,62}-UE ) \
$(shell echo Rivet-SppS-{53,63,200,500,900,546}-UE )
rm -rf Rivet-SppS
python/merge-SppS SppS
rivet-mkhtml -o Rivet-SppS SppS.yoda
Rivet-LHC-EW: $(shell echo Rivet-LHC{,-Matchbox,-Powheg,-Dipole}-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-SOPHTY},Z-LowMass-{e,mu},Z-MedMass-e,WZ,WW-{emu,ll},ZZ-{ll,lv},8-ZZ-lv,,8-WW-llW-Z-{e,mu}}) \
$(shell echo Rivet-LHC{,-Matchbox,-Dipole}-{7-W-Jet-{1..3}-e,7-Z-Jet-{0..3}-e,7-Z-Jet-0-mu}) \
$(shell echo Rivet-LHC{-Matchbox,-Dipole}-{Z-b,Z-bb,W-b,8-Z-jj}) \
$(shell echo Rivet-LHC-7-{W,Z}Gamma-{e,mu}) \
rm -rf Rivet-LHC-EW;
python/merge-LHC-EW LHC-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-Short},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda LHC-7-{W,Z}-Jet-{1,2,3}-e.yoda LHC-7-{W,Z}Gamma-{e,mu}.yoda -o LHC-EW.yoda;
python/merge-LHC-EW LHC-Matchbox-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-Short},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda LHC-Matchbox-7-{W,Z}-Jet-{1,2,3}-e.yoda -o LHC-Matchbox-EW.yoda;
python/merge-LHC-EW LHC-Dipole-{13-Z-{e,mu},8-Z-Mass{1..4}-{e,mu},W-{e,mu},Z-{e,mu,mu-Short},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda LHC-Dipole-7-{W,Z}-Jet-{1,2,3}-e.yoda -o LHC-Dipole-EW.yoda;
python/merge-LHC-EW LHC-Powheg-{W-{e,mu},Z-{e,mu},Z-LowMass-{e,mu},Z-MedMass-e,W-Z-{e,mu},WW-{emu,ll},WZ,ZZ-{ll,lv},8-ZZ-lv,8-WW-ll}.yoda -o LHC-Powheg-EW.yoda;
rivet-mkhtml -o Rivet-LHC-EW LHC-EW.yoda:Hw LHC-Powheg-EW.yoda:Hw-Powheg LHC-Matchbox-EW.yoda:Hw-Matchbox LHC-Matchbox-Z-b.yoda:Hw-Matchbox-Zb \
LHC-Matchbox-Z-bb.yoda:Hw-Matchbox-Zbb LHC-Matchbox-W-b.yoda:Hw-Matchbox-W-bb LHC-Dipole-EW.yoda:Hw-Dipole \
LHC-Dipole-Z-b.yoda:Hw-Dipole-Zb LHC-Dipole-Z-bb.yoda:Hw-Dipole-Zbb LHC-Dipole-W-b.yoda:Hw-Dipole-W-bb \
LHC-Z-mu-SOPHTY.yoda:Hw LHC-Powheg-Z-mu-SOPHTY.yoda:Hw-Powheg LHC-Matchbox-Z-mu-SOPHTY.yoda:Hw-Matchbox
Rivet-LHC-Photon: $(shell echo Rivet-LHC-{7,8}-PromptPhoton-{1..4}) Rivet-LHC-GammaGamma-7 \
$(shell echo Rivet-LHC{,-Powheg}-{7,8}-{DiPhoton-GammaGamma,DiPhoton-GammaJet})
rm -rf Rivet-LHC-Photon
python/merge-LHC-Photon LHC -o LHC-Photon.yoda
python/merge-LHC-Photon LHC-Powheg -o LHC-Powheg-Photon.yoda
rivet-mkhtml -o Rivet-LHC-Photon LHC-Photon.yoda:Hw LHC-Powheg-Photon.yoda:Hw-Powheg
Rivet-LHC-Higgs: $(shell echo Rivet-LHC{,-Powheg}-{ggH,VBF,WH,ZH})\
$(shell echo Rivet-LHC{,-Powheg}-8-{{ggH,VBF,WH,ZH}{,-GammaGamma},ggH-WW}) Rivet-LHC-ggHJet
rm -rf Rivet-LHC-Higgs
rivet-mkhtml -o Rivet-LHC-Higgs LHC-Powheg-ggH.yoda:gg-Powheg LHC-ggH.yoda:gg LHC-ggHJet.yoda:HJet \
LHC-Powheg-VBF.yoda:VBF-Powheg LHC-VBF.yoda:VBF LHC-WH.yoda:WH LHC-ZH.yoda:ZH \
LHC-Powheg-WH.yoda:WH-Powheg LHC-Powheg-ZH.yoda:ZH-Powheg
tests-Rivet : Rivet-LEP Rivet-BFactory Rivet-DIS Rivet-TVT-WZ Rivet-TVT-Photon Rivet-TVT-Jets Rivet-LHC-Jets Rivet-Star Rivet-SppS Rivet-LHC-EW Rivet-LHC-Photon Rivet-LHC-Higgs
test-Gamma-% : Inputs/Gamma-%.in GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
test-DIS-% : Inputs/DIS-%.in DISTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
if WANT_LIBFASTJET
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la HadronJetTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
else
test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la
$(HWREAD) $<
$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
endif
clean-local:
rm -f *.out *.log *.tex *.top *.run *.dump *.mult *.Bmult *.yoda
diff --git a/UnderlyingEvent/MPIHandler.h b/UnderlyingEvent/MPIHandler.h
--- a/UnderlyingEvent/MPIHandler.h
+++ b/UnderlyingEvent/MPIHandler.h
@@ -1,878 +1,878 @@
// -*- C++ -*-
//
// MPIHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MPIHandler_H
#define HERWIG_MPIHandler_H
//
// This is the declaration of the MPIHandler class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig/PDT/StandardMatchers.h"
#include "Herwig/Utilities/GSLBisection.h"
//#include "Herwig/Utilities/GSLMultiRoot.h"
#include "Herwig/Utilities/GSLIntegrator.h"
#include "Herwig/Shower/UEBase.h"
#include <cassert>
#include "ProcessHandler.h"
#include "MPIHandler.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup UnderlyingEvent
* \class MPIHandler
* This class is responsible for generating additional
* semi hard partonic interactions.
*
* \author Manuel B\"ahr
*
* @see \ref MPIHandlerInterfaces "The interfaces"
* defined for MPIHandler.
* @see ProcessHandler
* @see ShowerHandler
* @see HwRemDecayer
*/
class MPIHandler: public UEBase {
/**
* Maximum number of scatters
*/
static const unsigned int maxScatters_ = 99;
/**
* Class for the integration is a friend to access private members
*/
friend struct Eikonalization;
friend struct TotalXSecBisection;
friend struct slopeAndTotalXSec;
friend struct slopeInt;
friend struct slopeBisection;
public:
/** A vector of <code>SubProcessHandler</code>s. */
typedef vector<SubHdlPtr> SubHandlerList;
/** A vector of <code>Cut</code>s. */
typedef vector<CutsPtr> CutsList;
/** A vector of <code>ProcessHandler</code>s. */
typedef vector<ProHdlPtr> ProcessHandlerList;
/** A vector of cross sections. */
typedef vector<CrossSection> XSVector;
/** A pair of multiplicities: hard, soft. */
typedef pair<unsigned int, unsigned int> MPair;
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MPIHandler(): softMult_(0), identicalToUE_(-1),
PtOfQCDProc_(-1.0*GeV), Ptmin_(-1.0*GeV),
hardXSec_(0*millibarn), softXSec_(0*millibarn),
totalXSecExp_(0*millibarn),
softMu2_(ZERO), beta_(100.0/GeV2),
algorithm_(2), numSubProcs_(0),
colourDisrupt_(0.0), softInt_(true), twoComp_(true),
DLmode_(2), avgNhard_(0.0), avgNsoft_(0.0),
energyExtrapolation_(2), EEparamA_(0.6*GeV),
EEparamB_(37.5*GeV), refScale_(7000.*GeV),
pT0_(3.11*GeV), b_(0.21) {}
/**
* The destructor.
*/
virtual ~MPIHandler(){}
//@}
public:
/** @name Methods for the MPI generation. */
//@{
/*
* @return true if for this beam setup MPI can be generated
*/
virtual bool beamOK() const;
/**
* Return true or false depending on whether soft interactions are enabled.
*/
virtual bool softInt() const {return softInt_;}
/**
* Get the soft multiplicity from the pretabulated multiplicity
* distribution. Generated in multiplicity in the first place.
* @return the number of extra soft events in this collision
*/
virtual unsigned int softMultiplicity() const {return softMult_;}
/**
* Sample from the pretabulated multiplicity distribution.
* @return the number of extra events in this collision
*/
virtual unsigned int multiplicity(unsigned int sel=0);
/**
* Select a StandardXComb according to it's weight
* @return that StandardXComb Object
* @param sel is the subprocess that should be returned,
* if more than one is specified.
*/
virtual tStdXCombPtr generate(unsigned int sel=0);
//@}
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
/**
* Initialize this Multiple Interaction handler and all related objects needed to
* generate additional events.
*/
virtual void initialize();
/**
* Finalize this Multiple Interaction handler and all related objects needed to
* generate additional events.
*/
virtual void finalize();
/**
* Clean up the XCombs from our subprocesses after each event.
* ThePEG cannot see them, so the usual cleaning misses these.
*/
virtual void clean();
/**
* Write out accumulated statistics about integrated cross sections.
*/
void statistics() const;
/**
* The level of statistics. Controlls the amount of statistics
* written out after each run to the <code>EventGenerator</code>s
* <code>.out</code> file. Simply the EventHandler method is called here.
*/
int statLevel() const {return eventHandler()->statLevel();}
/**
* Return the hard cross section above ptmin
*/
CrossSection hardXSec() const { return hardXSec_; }
/**
* Return the soft cross section below ptmin
*/
CrossSection softXSec() const { return softXSec_; }
/**
* Return the inelastic cross section
*/
CrossSection inelasticXSec() const { return inelXSec_; }
/** @name Simple access functions. */
//@{
/**
* Return the ThePEG::EventHandler assigned to this handler.
* This methods shadows ThePEG::StepHandler::eventHandler(), because
* it is not virtual in ThePEG::StepHandler. This is ok, because this
* method would give a null-pointer at some stages, whereas this method
* gives access to the explicitely copied pointer (in initialize())
* to the ThePEG::EventHandler.
*/
tEHPtr eventHandler() const {return theHandler;}
/**
* Return the current handler
*/
static const MPIHandler * currentHandler() {
return currentHandler_;
}
/**
* Return theAlgorithm.
*/
virtual int Algorithm() const {return algorithm_;}
/**
* Return the ptmin parameter of the model
*/
virtual Energy Ptmin() const {
if(Ptmin_ > ZERO)
return Ptmin_;
else
throw Exception() << "MPIHandler::Ptmin called without initialize before"
<< Exception::runerror;
}
/**
* Return the slope of the soft pt spectrum as calculated.
*/
virtual InvEnergy2 beta() const {
if(beta_ != 100.0/GeV2)
return beta_;
else
throw Exception() << "MPIHandler::beta called without initialization"
<< Exception::runerror;
}
/**
* Return the pt Cutoff of the Interaction that is identical to the UE
* one.
*/
virtual Energy PtForVeto() const {return PtOfQCDProc_;}
/**
* Return the number of additional "hard" processes ( = multiple
* parton scattering)
*/
virtual unsigned int additionalHardProcs() const {return numSubProcs_-1;}
/**
* Return the fraction of colour disrupted connections to the
* suprocesses.
*/
virtual double colourDisrupt() const {return colourDisrupt_;}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* Access the list of sub-process handlers.
*/
const SubHandlerList & subProcesses()
const {return theSubProcesses;}
/**
* Access the list of sub-process handlers.
*/
SubHandlerList & subProcesses() {return theSubProcesses;}
/**
* Access the list of cuts.
*/
const CutsList & cuts() const {return theCuts;}
/**
* Access the list of cuts.
*/
CutsList & cuts() {return theCuts;}
/**
* Access the list of sub-process handlers.
*/
const ProcessHandlerList & processHandlers()
const {return theProcessHandlers;}
/**
* Access the list of sub-process handlers.
*/
ProcessHandlerList & processHandlers() {return theProcessHandlers;}
/**
* Method to calculate the individual probabilities for N scatters in the event.
* @param UEXSecs is(are) the inclusiv cross section(s) for the UE process(es).
*/
void Probs(XSVector UEXSecs);
/**
* Debug method to check the individual probabilities.
* @param filename is the file the output gets written to
*/
void MultDistribution(string filename) const;
/**
* Return the value of the Overlap function A(b) for a given impact
* parameter \a b.
* @param b impact parameter
* @param mu2 = inv hadron radius squared. 0 will use the value of
* invRadius_
* @return inverse area.
*/
InvArea OverlapFunction(Length b, Energy2 mu2=ZERO) const;
/**
* Method to calculate the poisson probability for expectation value
* \f$<n> = A(b)\sigma\f$, and multiplicity N.
*/
double poisson(Length b, CrossSection sigma,
unsigned int N, Energy2 mu2=ZERO) const;
/**
* Return n!
*/
double factorial (unsigned int n) const;
/**
* Returns the total cross section for the current CMenergy. The
* decision which parametrization will be used is steered by a
* external parameter of this class.
*/
CrossSection totalXSecExp() const;
/**
* Difference of the calculated total cross section and the
* experimental one from totalXSecExp.
* @param softXSec = the soft cross section that is used
* @param softMu2 = the soft radius, if 0 the hard radius will be used
*/
CrossSection totalXSecDiff(CrossSection softXSec,
Energy2 softMu2=ZERO) const;
/**
* Difference of the calculated elastic slope and the
* experimental one from slopeExp.
* @param softXSec = the soft cross section that is used
* @param softMu2 = the soft radius, if 0 the hard radius will be used
*/
InvEnergy2 slopeDiff(CrossSection softXSec,
Energy2 softMu2=ZERO) const;
/**
* Returns the value of the elastic slope for the current CMenergy.
* The decision which parametrization will be used is steered by a
* external parameter of this class.
*/
InvEnergy2 slopeExp() const;
/**
* Calculate the minimal transverse momentum from the extrapolation
*/
void overrideUECuts();
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MPIHandler & operator=(const MPIHandler &);
/**
* A pointer to the EventHandler that calls us. Has to be saved, because the
* method eventHandler() inherited from ThePEG::StepHandler returns a null-pointer
* sometimes. Leif changed that in r1053 so that a valid pointer is present, when
* calling doinitrun().
*/
tEHPtr theHandler;
/**
* The list of <code>SubProcessHandler</code>s.
*/
SubHandlerList theSubProcesses;
/**
* The kinematical cuts used for this collision handler.
*/
CutsList theCuts;
/**
* List of ProcessHandler used to sample different processes independently
*/
ProcessHandlerList theProcessHandlers;
/**
* A ThePEG::Selector where the individual Probabilities P_N are stored
* and the actual Multiplicities can be selected.
*/
Selector<MPair> theMultiplicities;
/**
* Variable to store the soft multiplicity generated for a event. This
* has to be stored as it is generated at the time of the hard
* additional interactions but used later on.
*/
unsigned int softMult_;
/**
* Variable to store the multiplicity of the second hard process
*/
vector<int> additionalMultiplicities_;
/**
* Variable to store the information, which process is identical to
* the UE one (QCD dijets).
* 0 means "real" hard one
* n>0 means the nth additional hard scatter
* -1 means no one!
*/
int identicalToUE_;
/**
* Variable to store the minimal pt of the process that is identical
* to the UE one. This only has to be set, if it can't be determined
* automatically (i.e. when reading QCD LesHouches files in).
*/
Energy PtOfQCDProc_;
/**
* Variable to store the parameter ptmin
*/
Energy Ptmin_;
/**
* Variable to store the hard cross section above ptmin
*/
CrossSection hardXSec_;
/**
* Variable to store the final soft cross section below ptmin
*/
CrossSection softXSec_;
/**
* Variable to store the inelastic cross section
*/
CrossSection inelXSec_;
/**
* Variable to store the total pp cross section (assuming rho=0!) as
* measured at LHC. If this variable is set, this value is used in the
* subsequent run instead of any of the Donnachie-Landshoff
* parametrizations.
*/
CrossSection totalXSecExp_;
/**
* Variable to store the soft radius, that is calculated during
* initialization for the two-component model.
*/
Energy2 softMu2_;
/**
* slope to the non-perturbative pt spectrum: \f$d\sigma/dp_T^2 = A \exp
* (- beta p_T^2)\f$. Its value is determined durint initialization.
*/
InvEnergy2 beta_;
/**
* Switch to be set from outside to determine the algorithm used for
* UE activity.
*/
int algorithm_;
/**
* Inverse hadron Radius squared \f$ (\mu^2) \f$. Used inside the overlap function.
*/
Energy2 invRadius_;
/**
* Member variable to store the actual number of separate SubProcesses
*/
unsigned int numSubProcs_;
/**
* Variable to store the relative number of colour disrupted
* connections to additional subprocesses. This variable is used in
* Herwig::HwRemDecayer but store here, to have access to all
* parameters through one Object.
*/
double colourDisrupt_;
/**
* Flag to store whether soft interactions, i.e. pt < ptmin should be
* simulated.
*/
bool softInt_;
/**
* Flag to steer wheather the soft part has a different radius, that
* will be dynamically fixed.
*/
bool twoComp_;
/**
* Switch to determine which Donnachie & Landshoff parametrization
* should be used.
*/
unsigned int DLmode_;
/**
* Variable to store the average hard multiplicity.
*/
double avgNhard_;
/**
* Variable to store the average soft multiplicity.
*/
double avgNsoft_;
/**
* The current handler
*/
static MPIHandler * currentHandler_;
/**
* Flag to store whether to calculate the minimal UE pt according to an
* extrapolation formula or whether to use MPIHandler:Cuts[0]:OneCuts[0]:MinKT
*/
unsigned int energyExtrapolation_;
/**
* Parameters for the energy extrapolation formula
*/
Energy EEparamA_;
Energy EEparamB_;
Energy refScale_;
Energy pT0_;
double b_;
protected:
/** @cond EXCEPTIONCLASSES */
/**
* Exception class used by the MultipleInteractionHandler, when something
* during initialization went wrong.
* \todo understand!!!
*/
class InitError: public Exception {};
/** @endcond */
};
}
namespace Herwig {
/**
* A struct for the 2D root finding that is necessary to determine the
* soft cross section and the soft radius that is needed to describe
* the total cross section correctly.
* NOT IN USE CURRENTLY
*/
struct slopeAndTotalXSec : public GSLHelper<CrossSection, CrossSection> {
public:
/**
* Constructor
*/
slopeAndTotalXSec(tcMPIHPtr handler): handler_(handler) {}
/** second argument type */
typedef Energy2 ArgType2;
/** second value type */
typedef InvEnergy2 ValType2;
/** first element of the vector like function to find root for
* @param softXSec soft cross-section
* @param softMu2 \f$\mu^2\f$
*/
CrossSection f1(ArgType softXSec, ArgType2 softMu2) const {
return handler_->totalXSecDiff(softXSec, softMu2);
}
/** second element of the vector like function to find root for
* @param softXSec soft cross-section
* @param softMu2 \f$\mu^2\f$
*/
InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const {
return handler_->slopeDiff(softXSec, softMu2);
}
/** provide the actual units of use */
virtual ValType vUnit() const {return 1.0*millibarn;}
/** otherwise rounding errors may get significant */
virtual ArgType aUnit() const {return 1.0*millibarn;}
/** provide the actual units of use */
ValType2 vUnit2() const {return 1.0/GeV2;}
/** otherwise rounding errors may get significant */
ArgType2 aUnit2() const {return GeV2;}
private:
/**
* Pointer to the handler
*/
tcMPIHPtr handler_;
};
/**
* A struct for the root finding that is necessary to determine the
* slope of the soft pt spectrum to match the soft cross section
*/
struct betaBisection : public GSLHelper<Energy2, InvEnergy2>{
public:
/**
* Constructor.
* @param soft = soft cross section, i.e. the integral of the soft
* pt spectrum f(u=p_T^2) = dsig exp(-beta*u/u_min)
* @param dsig = dsigma_hard/dp_T^2 at the p_T cutoff
* @param ptmin = p_T cutoff
*/
betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
: softXSec_(soft), dsig_(dsig), ptmin_(ptmin) {}
/**
* Operator that is used inside the GSLBisection class
*/
virtual Energy2 operator ()(InvEnergy2 beta) const
{
if( fabs(beta*GeV2) < 1.E-4 )
beta = (beta > ZERO) ? 1.E-4/GeV2 : -1.E-4/GeV2;
return (exp(beta*sqr(ptmin_)) - 1.0)/beta - softXSec_/dsig_;
}
/** provide the actual units of use */
virtual ValType vUnit() const {return 1.0*GeV2;}
/** provide the actual units of use */
virtual ArgType aUnit() const {return 1.0/GeV2;}
private:
/** soft cross section */
CrossSection softXSec_;
/** dsigma/dp_T^2 at ptmin */
DiffXSec dsig_;
/** pt cutoff */
Energy ptmin_;
};
/**
* A struct for the root finding that is necessary to determine the
* soft cross section and soft mu2 that are needed to describe the
* total cross section AND elastic slope correctly.
*/
struct slopeBisection : public GSLHelper<InvEnergy2, Energy2> {
public:
/** Constructor */
slopeBisection(tcMPIHPtr handler) : handler_(handler) {}
/**
* Return the difference of the calculated elastic slope to the
* experimental one for a given value of the soft mu2. During that,
* the soft cross section get fixed.
*/
InvEnergy2 operator ()(Energy2 arg) const;
/** Return the soft cross section that has been calculated */
CrossSection softXSec() const {return softXSec_;}
private:
/** const pointer to the MPIHandler to give access to member functions.*/
tcMPIHPtr handler_;
/** soft cross section that is determined on the fly.*/
mutable CrossSection softXSec_;
};
/**
* A struct for the root finding that is necessary to determine the
* soft cross section that is needed to describe the total cross
* section correctly.
*/
struct TotalXSecBisection : public GSLHelper<CrossSection, CrossSection> {
public:
/**
* Constructor
* @param handler The handler
* @param softMu2 \f$\mu^2\f$
*/
TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO):
handler_(handler), softMu2_(softMu2) {}
/**
* operator to return the cross section
* @param argument input cross section
*/
CrossSection operator ()(CrossSection argument) const {
return handler_->totalXSecDiff(argument, softMu2_);
}
/** provide the actual units of use */
virtual ValType vUnit() const {return 1.0*millibarn;}
/** otherwise rounding errors may get significant */
virtual ArgType aUnit() const {return 1.0*millibarn;}
private:
/**
* The handler
*/
tcMPIHPtr handler_;
/**
* \f$\mu^2\f$
*/
Energy2 softMu2_;
};
/**
* Typedef for derivative of the length
*/
- typedef Qty<1,-2,0> LengthDiff;
+ typedef decltype(mm/GeV2) LengthDiff;
/**
* A struct for the integrand for the slope
*/
struct slopeInt : public GSLHelper<LengthDiff, Length>{
public:
/** Constructor
* @param handler The handler
* @param hard The hard cross section
* @param soft The soft cross section
* @param softMu2 \f$\mu^2\f$
*/
slopeInt(tcMPIHPtr handler, CrossSection hard,
CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
: handler_(handler), hardXSec_(hard),
softXSec_(soft), softMu2_(softMu2) {}
/**
* Operator to return the answer
* @param arg The argument
*/
ValType operator ()(ArgType arg) const;
private:
/**
* Pointer to the Handler that calls this integrand
*/
tcMPIHPtr handler_;
/**
* The hard cross section to be eikonalized
*/
CrossSection hardXSec_;
/**
* The soft cross section to be eikonalized. Default is zero
*/
CrossSection softXSec_;
/**
* The inv radius^2 of the soft interactions.
*/
Energy2 softMu2_;
};
/**
* A struct for the eikonalization of the inclusive cross section.
*/
struct Eikonalization : public GSLHelper<Length, Length>{
/**
* The constructor
* @param handler is the pointer to the MPIHandler to get access to
* MPIHandler::OverlapFunction and member variables of the MPIHandler.
* @param option is a flag, whether the inelastic or the total
* @param handler The handler
* @param hard The hard cross section
* @param soft The soft cross section
* @param softMu2 \f$\mu^2\f$
* cross section should be returned (-2 or -1). For option = N > 0 the integrand
* is N*(A(b)*sigma)^N/N! exp(-A(b)*sigma) this is the P_N*sigma where
* P_N is the Probability of having exactly N interaction (including the hard one)
* This is equation 14 from "Jimmy4: Multiparton Interactions in HERWIG for the LHC"
*/
Eikonalization(tcMPIHPtr handler, int option, CrossSection hard,
CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
: theHandler(handler), theoption(option), hardXSec_(hard),
softXSec_(soft), softMu2_(softMu2) {}
/**
* Get the function value
*/
Length operator ()(Length argument) const;
private:
/**
* Pointer to the Handler that calls this integrand
*/
tcMPIHPtr theHandler;
/**
* A flag to switch between the calculation of total and inelastic cross section
* or calculations for the individual probabilities. See the constructor
*/
int theoption;
/**
* The hard cross section to be eikonalized
*/
CrossSection hardXSec_;
/**
* The soft cross section to be eikonalized. Default is zero
*/
CrossSection softXSec_;
/**
* The inv radius^2 of the soft interactions.
*/
Energy2 softMu2_;
};
}
#endif /* HERWIG_MPIHandler_H */
diff --git a/Utilities/GSLIntegrator.h b/Utilities/GSLIntegrator.h
--- a/Utilities/GSLIntegrator.h
+++ b/Utilities/GSLIntegrator.h
@@ -1,120 +1,122 @@
// -*- C++ -*-
//
// GSLIntegrator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_GSLIntegrator_H
#define HERWIG_GSLIntegrator_H
//
// This is the declaration of the GSLIntegrator class.
//
#include "ThePEG/Pointer/ReferenceCounted.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "gsl/gsl_integration.h"
#include "gsl/gsl_errno.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Utilities
* This class is designed to integrate a given function between
* 2 limits using the gsl QAGS integration subroutine.
*
* The function is supplied using a templated class that must define
* operator(argument). The units of the argument ArgType and return
* type ValType must be supplied in the integrand class using a typedef
* i.e. <br>
* <code> struct integrand { </code><br>
* <code> ... </code> <BR>
* <code>Energy operator(double arg) const;</code><BR>
* <code>typedef double ArgType</code><BR>
* <code>typedef Energy ValType</code><BR>
* <code> ... </code> <BR>
* <code>}</code> <BR>
*/
class GSLIntegrator : public Pointer::ReferenceCounted {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default Constructor uses values in GSL manual as parameters
**/
GSLIntegrator() : _abserr(1.0E-35), _relerr(1.0E-3), _nbins(1000) {}
/**
* Specify all the parameters.
* @param abserr Absolute error.
* @param relerr Relative error.
* @param nbins Number of bins
*/
GSLIntegrator(double abserr, double relerr, int nbins) :
_abserr(abserr), _relerr(relerr), _nbins(nbins) {}
//@}
+ /// helper type for the integration result
+ template <class T>
+ using ValT = decltype(std::declval<typename T::ValType>()
+ * std::declval<typename T::ArgType>());
+
/**
* The value of the integral
* @param function The integrand class that defines operator()
* @param lower The lower limit of integration.
* @param upper The upper limit of integration.
*/
template <class T>
- inline typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT
+ inline ValT<T>
value(const T & function,
const typename T::ArgType lower,
const typename T::ArgType upper) const;
/**
* The value of the integral
* @param function The integrand class that defines operator()
* @param lower The lower limit of integration.
* @param upper The upper limit of integration.
* @param error Returns the estimated error of the integral
*/
template <class T>
- inline typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT
+ inline ValT<T>
value(const T & function,
const typename T::ArgType lower,
const typename T::ArgType upper,
- typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT & error) const;
+ ValT<T> & error) const;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
GSLIntegrator & operator=(const GSLIntegrator &);
private:
/**
* The parameters controlling the absolute error.
*/
double _abserr;
/**
* The parameters controlling the relative error.
*/
double _relerr;
/**
* The maximum number of intervals to use.
*/
int _nbins;
};
}
#include "GSLIntegrator.tcc"
#endif /* HERWIG_GSLIntegrator_H */
diff --git a/Utilities/GSLIntegrator.tcc b/Utilities/GSLIntegrator.tcc
--- a/Utilities/GSLIntegrator.tcc
+++ b/Utilities/GSLIntegrator.tcc
@@ -1,118 +1,116 @@
// -*- C++ -*-
//
// GSLIntegrator.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined templated member
// functions of the GSLIntegrator class.
//
using namespace Herwig;
using namespace ThePEG;
namespace {
template <class T> struct param {
//The integrand function
const T & function;
};
template<class T> double integrand(double x , void * p) {
//Units of the argument and return type
const typename T::ValType ValUnit =
TypeTraits<typename T::ValType>::baseunit();
const typename T::ArgType ArgUnit =
TypeTraits<typename T::ArgType>::baseunit();
const T & f = ((struct param<T> *)p)->function;
return f(x * ArgUnit ) / ValUnit;
}
}
namespace Herwig {
using namespace ThePEG;
template <class T>
-inline typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT
+inline GSLIntegrator::ValT<T>
GSLIntegrator::value(const T & fn,
const typename T::ArgType lower,
- const typename T::ArgType upper) const {
- typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT error;
+ const typename T::ArgType upper) const
+{
+ GSLIntegrator::ValT<T> error;
return value(fn,lower,upper,error);
}
template <class T>
-inline typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT
+inline GSLIntegrator::ValT<T>
GSLIntegrator::value(const T & fn,
const typename T::ArgType lower,
const typename T::ArgType upper,
- typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT & error) const {
+ GSLIntegrator::ValT<T> & error) const
+{
typedef typename T::ValType ValType;
typedef typename T::ArgType ArgType;
const ValType ValUnit = TypeTraits<ValType>::baseunit();
const ArgType ArgUnit = TypeTraits<ArgType>::baseunit();
double result(0.), error2(0.);
param<T> parameters = { fn };
gsl_function integrationFunction;
integrationFunction.function = &integrand<T>;
integrationFunction.params = &parameters;
gsl_integration_workspace * workspace =
gsl_integration_workspace_alloc(_nbins);
//do integration
//Want to check error messages ourselves
gsl_error_handler_t * oldhandler = gsl_set_error_handler_off();
int status = gsl_integration_qags(&integrationFunction, lower/ArgUnit,
upper/ArgUnit, _abserr, _relerr, _nbins,
workspace, &result, &error2);
if( status > 0 ) {
CurrentGenerator::log() << "An error occurred in the GSL "
"integration subroutine:\n";
switch( status ) {
case GSL_EMAXITER:
CurrentGenerator::log() << "The maximum number of subdivisions "
"was exceeded.\n";
break;
case GSL_EROUND:
CurrentGenerator::log() << "Cannot reach tolerance because of "
"roundoff error, or roundoff error was detected in the "
"extrapolation table.\n";
break;
case GSL_ESING:
CurrentGenerator::log() << "A non-integrable singularity or "
"other bad integrand behavior was found in the integration "
"interval.\n";
break;
case GSL_EDIVERGE:
CurrentGenerator::log() << "The integral is divergent, "
"or too slowly convergent to be integrated numerically.\n";
break;
default:
CurrentGenerator::log() << "A general error occurred with code "
<< status << '\n';
}
result = 0.;
}
gsl_set_error_handler(oldhandler);
gsl_integration_workspace_free(workspace);
//fix units and return
error = error2* ValUnit * ArgUnit;
return result * ValUnit * ArgUnit;
}
}
diff --git a/Utilities/GaussianIntegrator.h b/Utilities/GaussianIntegrator.h
--- a/Utilities/GaussianIntegrator.h
+++ b/Utilities/GaussianIntegrator.h
@@ -1,128 +1,132 @@
// -*- C++ -*-
//
// GaussianIntegrator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_GaussianIntegrator_H
#define HERWIG_GaussianIntegrator_H
//
// This is the declaration of the GaussianIntegrator class.
//
#include "ThePEG/Pointer/ReferenceCounted.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include <vector>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Utilities
* \author Peter Richardson
* This class is designed to perform the integral of a function
* using Gaussian quadrature.The method is adaptive based on using 6th,12th,
* 24th,48th, or 96th order Gaussian quadrature combined with
* subdivision of the integral if this is insufficient.
*
* The class is templated on a simple class which should provide a
* T::operator () (double) const which provides the integrand for the function.
*/
class GaussianIntegrator : public Pointer::ReferenceCounted {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default Constructor
*/
GaussianIntegrator()
: _abserr(1.E-35), _relerr(5.E-5), _binwidth(1.E-5),
_maxint(100), _maxeval(100000) {
// setup the weights and abscissae
Init();
}
/**
* Specify all the parameters.
* @param abserr Absolute error.
* @param relerr Relative error.
* @param binwidth Width of the bin as a fraction of the integration region.
* @param maxint Maximum number of intervals
* @param maxeval Maximum number of function evaluations
*/
GaussianIntegrator(double abserr, double relerr, double binwidth,
int maxint, int maxeval)
- : _abserr(abserr), _relerr(relerr), _binwidth(binwidth), _maxint(maxint),
+ : _abserr(abserr), _relerr(relerr),
+ _binwidth(binwidth), _maxint(maxint),
_maxeval(maxeval) {
// setup the weights and abscissae
Init();
}
+ /// helper type for the integration result
+ template <class T>
+ using ValT = decltype(std::declval<typename T::ValType>()
+ * std::declval<typename T::ArgType>());
+
/**
* The value of the integral
* @param lower The lower limit of integration.
* @param upper The upper limit of integration.
*/
template <class T>
- inline typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT
- value(const T &,
- const typename T::ArgType lower,
- const typename T::ArgType upper) const;
+ inline ValT<T> value(const T &,
+ const typename T::ArgType lower,
+ const typename T::ArgType upper) const;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
GaussianIntegrator & operator=(const GaussianIntegrator &);
/**
* Initialise the weights and abscissae.
*/
void Init();
private:
/**
* The weights for the gaussian quadrature.
*/
std::vector< std::vector<double> > _weights;
/**
* The abscissae.
*/
std::vector< std::vector <double> > _abscissae;
/**
* The parameters controlling the error.
*/
double _abserr,_relerr;
/**
* The minimum width of a bin as a fraction of the integration region.
*/
double _binwidth;
/**
* Maximum number of bins.
*/
int _maxint;
/**
* Maximum number of function evaluations.
*/
int _maxeval;
};
}
#include "GaussianIntegrator.tcc"
#endif /* HERWIG_GaussianIntegrator_H */
diff --git a/Utilities/GaussianIntegrator.tcc b/Utilities/GaussianIntegrator.tcc
--- a/Utilities/GaussianIntegrator.tcc
+++ b/Utilities/GaussianIntegrator.tcc
@@ -1,114 +1,113 @@
// -*- C++ -*-
//
// GaussianIntegrator.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined templated member
// functions of the GaussianIntegrator class.
//
namespace Herwig {
using namespace ThePEG;
template <class T>
-inline typename BinaryOpTraits<typename T::ValType,
- typename T::ArgType>::MulT
+inline GaussianIntegrator::ValT<T>
GaussianIntegrator::value(const T & function,
const typename T::ArgType lower,
const typename T::ArgType upper) const {
typedef typename T::ValType ValType;
typedef typename T::ArgType ArgType;
const ValType ValUnit = TypeTraits<ValType>::baseunit();
const ArgType ArgUnit = TypeTraits<ArgType>::baseunit();
// vector for the limits of the bin
vector<double> lowerlim,upperlim;
// start with the whole interval as 1 bin
lowerlim.push_back(lower/ArgUnit);upperlim.push_back(upper/ArgUnit);
// set the minimum bin width
double xmin=_binwidth*abs(upper-lower)/ArgUnit;
// counters for the number of function evals
int neval=0;
// and number of bad intervals
int nbad=0;
// the output value
double output=0.;
// the loop for the evaluation
double mid,wid; unsigned int ibin,ix=0,iorder;
double testvalue,value,tolerance;
do {
// the bin we are doing (always the last one in the list)
ibin = lowerlim.size()-1;
// midpoint and width of the bin
mid=0.5*(upperlim[ibin]+lowerlim[ibin]);
wid=0.5*(upperlim[ibin]-lowerlim[ibin]);
value=0.;
iorder=0;
// compute a trail value using sixth order GQ
for(ix=0;ix<_weights[0].size();++ix) {
value+=_weights[0][ix]
*( function((mid+wid*_abscissae[0][ix])*ArgUnit)
+function((mid-wid*_abscissae[0][ix])*ArgUnit)
)/ValUnit;
++neval;
if(neval>_maxeval)
CurrentGenerator::log() << "Error in Gaussian Integrator: Setting to zero"
<< endl;
}
value *=wid;
// compute more accurate answers using higher order GQ
do {
// use the next order of quadrature
testvalue=value;
++iorder;
value=0.;
for(ix=0;ix<_weights[iorder].size();++ix) {
value+=_weights[iorder][ix]*
( function((mid+wid*_abscissae[iorder][ix])*ArgUnit)
+function((mid-wid*_abscissae[iorder][ix])*ArgUnit)
)/ValUnit;
++neval;
if(neval>_maxeval)
CurrentGenerator::log() << "Error in Gaussian Integrator: Setting to zero"
<< endl;
}
value *=wid;
tolerance=max(_abserr,_relerr*abs(value));
}
// keep going if possible and not accurate enough
while(iorder<_weights.size()-1&&abs(testvalue-value)>tolerance);
// now decide what to do
// accept this value
if(abs(testvalue-value)<tolerance) {
output+=value;
lowerlim.pop_back();upperlim.pop_back();
}
// bin too small to redivide contribution set to zero
else if(wid<xmin) {
++nbad;
lowerlim.pop_back(); upperlim.pop_back();
}
// otherwise split the bin into two
else {
// reset the limits for the bin
upperlim[ibin]=mid;
// set up a new bin
lowerlim.push_back(mid);
upperlim.push_back(mid+wid);
}
}
// keep going if there's still some bins to evaluate
while(lowerlim.size()>0);
// output an error message if needed
if(nbad!=0)
CurrentGenerator::log() << "Error in GaussianIntegrator: Bad Convergence for "
<< nbad << "intervals" << endl;
// return the answer
return output * ValUnit * ArgUnit;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:12 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805985
Default Alt Text
(107 KB)

Event Timeline