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/LHC-FR.in.template b/Models/Feynrules/python/ufo2peg/LHC-FR.in.template
--- a/Models/Feynrules/python/ufo2peg/LHC-FR.in.template
+++ b/Models/Feynrules/python/ufo2peg/LHC-FR.in.template
@@ -1,74 +1,92 @@
read snippets/PPCollider.in
read ${ModelName}.model
cd /Herwig/NewPhysics
####################################
#
# Modify the required process here
#
####################################
+insert HPConstructor:Incoming 0 /Herwig/Particles/d
+insert HPConstructor:Incoming 0 /Herwig/Particles/dbar
insert HPConstructor:Incoming 0 /Herwig/Particles/u
-insert HPConstructor:Incoming 1 /Herwig/Particles/ubar
+insert HPConstructor:Incoming 0 /Herwig/Particles/ubar
+insert HPConstructor:Incoming 0 /Herwig/Particles/s
+insert HPConstructor:Incoming 0 /Herwig/Particles/sbar
+insert HPConstructor:Incoming 0 /Herwig/Particles/c
+insert HPConstructor:Incoming 0 /Herwig/Particles/cbar
+insert HPConstructor:Incoming 0 /Herwig/Particles/b
+insert HPConstructor:Incoming 0 /Herwig/Particles/bbar
+insert HPConstructor:Incoming 0 /Herwig/Particles/g
${Particles}
set HPConstructor:Processes SingleParticleInclusive
#############################################################
## Additionally, you can use new particles as intermediates
## with the ResConstructor:
#############################################################
+# insert ResConstructor:Incoming 0 /Herwig/Particles/d
+# insert ResConstructor:Incoming 0 /Herwig/Particles/dbar
# insert ResConstructor:Incoming 0 /Herwig/Particles/u
# insert ResConstructor:Incoming 0 /Herwig/Particles/ubar
+# insert ResConstructor:Incoming 0 /Herwig/Particles/s
+# insert ResConstructor:Incoming 0 /Herwig/Particles/sbar
+# insert ResConstructor:Incoming 0 /Herwig/Particles/c
+# insert ResConstructor:Incoming 0 /Herwig/Particles/cbar
+# insert ResConstructor:Incoming 0 /Herwig/Particles/b
+# insert ResConstructor:Incoming 0 /Herwig/Particles/bbar
+# insert ResConstructor:Incoming 0 /Herwig/Particles/g
#
# insert ResConstructor:Intermediates 0 [PARTICLE_NAME]
#
# insert ResConstructor:Outgoing 0 /Herwig/Particles/e+
# insert ResConstructor:Outgoing 1 /Herwig/Particles/W+
# insert ResConstructor:Outgoing 2 /Herwig/Particles/Z0
# insert ResConstructor:Outgoing 3 /Herwig/Particles/gamma
###########################################################
# Specialized 2->3 higgs constructors are also available,
# where incoming lines don't need to be set.
###########################################################
## Higgs + t tbar
# set /Herwig/NewPhysics/QQHiggsConstructor:QuarkType Top
# insert /Herwig/NewPhysics/QQHiggsConstructor:HiggsBoson 0 [HIGGS_NAME]
#
## Higgs VBF
# insert /Herwig/NewPhysics/HiggsVBFConstructor:HiggsBoson 0 [HIGGS_NAME]
#
## Higgs + W/Z, with full 2->3 ME
# set /Herwig/NewPhysics/HVConstructor:CollisionType Hadron
# insert /Herwig/NewPhysics/HVConstructor:VectorBoson 0 /Herwig/Particles/Z0
# insert /Herwig/NewPhysics/HVConstructor:HiggsBoson 0 [HIGGS_NAME]
####################################
####################################
####################################
# Intrinsic pT tune extrapolated to LHC energy
set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 2.2*GeV
# disable default cuts if required
# cd /Herwig/EventHandlers
# create ThePEG::Cuts /Herwig/Cuts/NoCuts
# set EventHandler:Cuts /Herwig/Cuts/NoCuts
# Other parameters for run
cd /Herwig/Generators
set EventGenerator:EventHandler:LuminosityFunction:Energy 13000.0
set EventGenerator:NumberOfEvents 10000000
set EventGenerator:RandomNumberGenerator:Seed 31122001
set EventGenerator:DebugLevel 0
set EventGenerator:EventHandler:StatLevel Full
set EventGenerator:PrintEvent 100
set EventGenerator:MaxErrors 10000
saverun LHC-${ModelName} EventGenerator
diff --git a/Shower/Core/Base/ShowerParticle.cc b/Shower/Core/Base/ShowerParticle.cc
--- a/Shower/Core/Base/ShowerParticle.cc
+++ b/Shower/Core/Base/ShowerParticle.cc
@@ -1,595 +1,584 @@
// -*- C++ -*-
//
// ShowerParticle.cc 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, non-templated member
// functions of the ShowerParticle class.
//
#include "ShowerParticle.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include <ostream>
using namespace Herwig;
PPtr ShowerParticle::clone() const {
return new_ptr(*this);
}
PPtr ShowerParticle::fullclone() const {
return new_ptr(*this);
}
ClassDescription<ShowerParticle> ShowerParticle::initShowerParticle;
// Definition of the static class description member.
void ShowerParticle::vetoEmission(ShowerPartnerType, Energy scale) {
scales_.QED = min(scale,scales_.QED );
scales_.QED_noAO = min(scale,scales_.QED_noAO );
scales_.QCD_c = min(scale,scales_.QCD_c );
scales_.QCD_c_noAO = min(scale,scales_.QCD_c_noAO );
scales_.QCD_ac = min(scale,scales_.QCD_ac );
scales_.QCD_ac_noAO = min(scale,scales_.QCD_ac_noAO);
}
void ShowerParticle::addPartner(EvolutionPartner in ) {
partners_.push_back(in);
}
namespace {
LorentzRotation boostToShower(Lorentz5Momentum & porig, tShowerBasisPtr basis) {
LorentzRotation output;
assert(basis);
if(basis->frame()==ShowerBasis::BackToBack) {
// we are doing the evolution in the back-to-back frame
// work out the boostvector
Boost boostv(-(basis->pVector()+basis->nVector()).boostVector());
// momentum of the parton
Lorentz5Momentum ptest(basis->pVector());
// construct the Lorentz boost
output = LorentzRotation(boostv);
ptest *= output;
Axis axis(ptest.vect().unit());
// now rotate so along the z axis as needed for the splitting functions
if(axis.perp2()>1e-10) {
double sinth(sqrt(1.-sqr(axis.z())));
output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
else if(axis.z()<0.) {
output.rotate(Constants::pi,Axis(1.,0.,0.));
}
porig = output*basis->pVector();
porig.setX(ZERO);
porig.setY(ZERO);
}
else {
output = LorentzRotation(-basis->pVector().boostVector());
porig = output*basis->pVector();
porig.setX(ZERO);
porig.setY(ZERO);
porig.setZ(ZERO);
}
return output;
}
RhoDMatrix bosonMapping(ShowerParticle & particle,
const Lorentz5Momentum & porig,
VectorSpinPtr vspin,
const LorentzRotation & rot) {
// rotate the original basis
vector<LorentzPolarizationVector> sbasis;
for(unsigned int ix=0;ix<3;++ix) {
sbasis.push_back(vspin->getProductionBasisState(ix));
sbasis.back().transform(rot);
}
// splitting basis
vector<LorentzPolarizationVector> fbasis;
bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
- VectorWaveFunction wave(porig,particle.dataPtr(),outgoing);
+ VectorWaveFunction wave(porig,particle.dataPtr(),incoming);
for(unsigned int ix=0;ix<3;++ix) {
if(massless&&ix==1) {
fbasis.push_back(LorentzPolarizationVector());
}
else {
- wave.reset(ix);
+ wave.reset(ix,vector_phase);
fbasis.push_back(wave.wave());
}
}
// work out the mapping
RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
- mapping(ix,iy)= sbasis[iy].dot(fbasis[ix].conjugate());
+ mapping(ix,iy)= -sbasis[iy].conjugate().dot(fbasis[ix]);
if(particle.id()<0)
mapping(ix,iy)=conj(mapping(ix,iy));
}
}
- // \todo need to fix this
- mapping = RhoDMatrix(PDT::Spin1,false);
- if(massless) {
- mapping(0,0) = 1.;
- mapping(2,2) = 1.;
- }
- else {
- mapping(0,0) = 1.;
- mapping(1,1) = 1.;
- mapping(2,2) = 1.;
- }
return mapping;
}
RhoDMatrix fermionMapping(ShowerParticle & particle,
const Lorentz5Momentum & porig,
FermionSpinPtr fspin,
const LorentzRotation & rot) {
// extract the original basis states
vector<LorentzSpinor<SqrtEnergy> > sbasis;
for(unsigned int ix=0;ix<2;++ix) {
sbasis.push_back(fspin->getProductionBasisState(ix));
sbasis.back().transform(rot);
}
// calculate the states in the splitting basis
vector<LorentzSpinor<SqrtEnergy> > fbasis;
SpinorWaveFunction wave(porig,particle.dataPtr(),
particle.id()>0 ? incoming : outgoing);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
fbasis.push_back(wave.dimensionedWave());
}
RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
for(unsigned int ix=0;ix<2;++ix) {
if(fbasis[0].s2()==complex<SqrtEnergy>()) {
mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
}
else {
mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
}
}
return mapping;
}
FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
const Lorentz5Momentum & porig,
const LorentzRotation & rot,
Helicity::Direction dir) {
// calculate the splitting basis for the branching
// and rotate back to construct the basis states
LorentzRotation rinv = rot.inverse();
SpinorWaveFunction wave;
if(particle.id()>0)
wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
else
wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
basis.transform(rinv);
fspin->setBasisState(ix,basis);
fspin->setDecayState(ix,basis);
}
particle.spinInfo(fspin);
return fspin;
}
VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
const Lorentz5Momentum & porig,
const LorentzRotation & rot,
Helicity::Direction dir) {
// calculate the splitting basis for the branching
// and rotate back to construct the basis states
LorentzRotation rinv = rot.inverse();
bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
VectorWaveFunction wave(porig,particle.dataPtr(),dir);
VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
for(unsigned int ix=0;ix<3;++ix) {
LorentzPolarizationVector basis;
if(massless&&ix==1) {
basis = LorentzPolarizationVector();
}
else {
- wave.reset(ix);
+ wave.reset(ix,vector_phase);
basis = wave.wave();
}
basis *= rinv;
vspin->setBasisState(ix,basis);
vspin->setDecayState(ix,basis);
}
particle.spinInfo(vspin);
vspin-> DMatrix() = RhoDMatrix(PDT::Spin1);
vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
if(massless) {
vspin-> DMatrix()(0,0) = 0.5;
vspin->rhoMatrix()(0,0) = 0.5;
vspin-> DMatrix()(2,2) = 0.5;
vspin->rhoMatrix()(2,2) = 0.5;
}
return vspin;
}
}
RhoDMatrix ShowerParticle::extractRhoMatrix(bool forward) {
// get the spin density matrix and the mapping
RhoDMatrix mapping;
SpinPtr inspin;
bool needMapping = getMapping(inspin,mapping);
// set the decayed flag
inspin->decay();
// get the spin density matrix
RhoDMatrix rho = forward ? inspin->rhoMatrix() : inspin->DMatrix();
// map to the shower basis if needed
if(needMapping) {
RhoDMatrix rhop(rho.iSpin(),false);
for(int ixb=0;ixb<rho.iSpin();++ixb) {
for(int iyb=0;iyb<rho.iSpin();++iyb) {
if(mapping(iyb,ixb)==0.)continue;
for(int iya=0;iya<rho.iSpin();++iya) {
if(rho(iya,iyb)==0.)continue;
for(int ixa=0;ixa<rho.iSpin();++ixa) {
rhop(ixa,ixb) += rho(iya,iyb)*mapping(iya,ixa)*conj(mapping(iyb,ixb));
}
}
}
}
rhop.normalize();
rho = rhop;
}
return rho;
}
bool ShowerParticle::getMapping(SpinPtr & output, RhoDMatrix & mapping) {
// if the particle is not from the hard process
if(!this->perturbative()) {
// mapping is the identity
output=this->spinInfo();
mapping=RhoDMatrix(this->dataPtr()->iSpin());
if(output) {
return false;
}
else {
Lorentz5Momentum porig;
LorentzRotation rot = boostToShower(porig,showerBasis());
Helicity::Direction dir = this->isFinalState() ? outgoing : incoming;
if(this->dataPtr()->iSpin()==PDT::Spin0) {
assert(false);
}
else if(this->dataPtr()->iSpin()==PDT::Spin1Half) {
output = createFermionSpinInfo(*this,porig,rot,dir);
}
else if(this->dataPtr()->iSpin()==PDT::Spin1) {
output = createVectorSpinInfo(*this,porig,rot,dir);
}
else {
assert(false);
}
return false;
}
}
// if particle is final-state and is from the hard process
else if(this->isFinalState()) {
assert(this->perturbative()==1 || this->perturbative()==2);
// get transform to shower frame
Lorentz5Momentum porig;
LorentzRotation rot = boostToShower(porig,showerBasis());
// the rest depends on the spin of the particle
PDT::Spin spin(this->dataPtr()->iSpin());
mapping=RhoDMatrix(spin,false);
// do the spin dependent bit
if(spin==PDT::Spin0) {
ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(this->spinInfo());
if(!sspin) {
ScalarWaveFunction::constructSpinInfo(this,outgoing,true);
}
output=this->spinInfo();
return false;
}
else if(spin==PDT::Spin1Half) {
FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
// spin info exists get information from it
if(fspin) {
output=fspin;
mapping = fermionMapping(*this,porig,fspin,rot);
return true;
}
// spin info does not exist create it
else {
output = createFermionSpinInfo(*this,porig,rot,outgoing);
return false;
}
}
else if(spin==PDT::Spin1) {
VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
// spin info exists get information from it
if(vspin) {
output=vspin;
mapping = bosonMapping(*this,porig,vspin,rot);
return true;
}
else {
output = createVectorSpinInfo(*this,porig,rot,outgoing);
return false;
}
}
// not scalar/fermion/vector
else
assert(false);
}
// incoming to hard process
else if(this->perturbative()==1 && !this->isFinalState()) {
// get the basis vectors
// get transform to shower frame
Lorentz5Momentum porig;
LorentzRotation rot = boostToShower(porig,showerBasis());
porig *= this->x();
// the rest depends on the spin of the particle
PDT::Spin spin(this->dataPtr()->iSpin());
mapping=RhoDMatrix(spin);
// do the spin dependent bit
if(spin==PDT::Spin0) {
cerr << "testing spin 0 not yet implemented " << endl;
assert(false);
}
// spin-1/2
else if(spin==PDT::Spin1Half) {
FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
// spin info exists get information from it
if(fspin) {
output=fspin;
mapping = fermionMapping(*this,porig,fspin,rot);
return true;
}
// spin info does not exist create it
else {
output = createFermionSpinInfo(*this,porig,rot,incoming);
return false;
}
}
// spin-1
else if(spin==PDT::Spin1) {
VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
// spinInfo exists map it
if(vspin) {
output=vspin;
mapping = bosonMapping(*this,porig,vspin,rot);
return true;
}
// create the spininfo
else {
output = createVectorSpinInfo(*this,porig,rot,incoming);
return false;
}
}
assert(false);
}
// incoming to decay
else if(this->perturbative() == 2 && !this->isFinalState()) {
// get the basis vectors
Lorentz5Momentum porig;
LorentzRotation rot=boostToShower(porig,showerBasis());
// the rest depends on the spin of the particle
PDT::Spin spin(this->dataPtr()->iSpin());
mapping=RhoDMatrix(spin);
// do the spin dependent bit
if(spin==PDT::Spin0) {
cerr << "testing spin 0 not yet implemented " << endl;
assert(false);
}
// spin-1/2
else if(spin==PDT::Spin1Half) {
// FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
// // spin info exists get information from it
// if(fspin) {
// output=fspin;
// mapping = fermionMapping(*this,porig,fspin,rot);
// return true;
// // spin info does not exist create it
// else {
// output = createFermionSpinInfo(*this,porig,rot,incoming);
// return false;
// }
// }
assert(false);
}
// // spin-1
// else if(spin==PDT::Spin1) {
// VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
// // spinInfo exists map it
// if(vspin) {
// output=vspin;
// mapping = bosonMapping(*this,porig,vspin,rot);
// return true;
// }
// // create the spininfo
// else {
// output = createVectorSpinInfo(*this,porig,rot,incoming);
// return false;
// }
// }
// assert(false);
assert(false);
}
else
assert(false);
return true;
}
void ShowerParticle::constructSpinInfo(bool timeLike) {
// now construct the required spininfo and calculate the basis states
PDT::Spin spin(dataPtr()->iSpin());
if(spin==PDT::Spin0) {
ScalarWaveFunction::constructSpinInfo(this,outgoing,timeLike);
}
// calculate the basis states and construct the SpinInfo for a spin-1/2 particle
else if(spin==PDT::Spin1Half) {
// outgoing particle
if(id()>0) {
vector<LorentzSpinorBar<SqrtEnergy> > stemp;
SpinorBarWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
SpinorBarWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
}
// outgoing antiparticle
else {
vector<LorentzSpinor<SqrtEnergy> > stemp;
SpinorWaveFunction::calculateWaveFunctions(stemp,this,outgoing);
SpinorWaveFunction::constructSpinInfo(stemp,this,outgoing,timeLike);
}
}
// calculate the basis states and construct the SpinInfo for a spin-1 particle
else if(spin==PDT::Spin1) {
bool massless(id()==ParticleID::g||id()==ParticleID::gamma);
vector<Helicity::LorentzPolarizationVector> vtemp;
- VectorWaveFunction::calculateWaveFunctions(vtemp,this,outgoing,massless);
+ VectorWaveFunction::calculateWaveFunctions(vtemp,this,outgoing,massless,vector_phase);
VectorWaveFunction::constructSpinInfo(vtemp,this,outgoing,timeLike,massless);
}
else {
throw Exception() << "Spins higher than 1 are not yet implemented in "
<< "FS_QtildaShowerKinematics1to2::constructVertex() "
<< Exception::runerror;
}
}
void ShowerParticle::initializeDecay() {
Lorentz5Momentum p, n, ppartner, pcm;
assert(perturbative()!=1);
// this is for the initial decaying particle
if(perturbative()==2) {
ShowerBasisPtr newBasis;
p = momentum();
Lorentz5Momentum ppartner(partner()->momentum());
// removed to make inverse recon work properly
//if(partner->thePEGBase()) ppartner=partner->thePEGBase()->momentum();
pcm=ppartner;
Boost boost(p.findBoostToCM());
pcm.boost(boost);
n = Lorentz5Momentum( ZERO,0.5*p.mass()*pcm.vect().unit());
n.boost( -boost);
newBasis = new_ptr(ShowerBasis());
newBasis->setBasis(p,n,ShowerBasis::Rest);
showerBasis(newBasis,false);
}
else {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
}
}
void ShowerParticle::initializeInitialState(PPtr parent) {
// For the time being we are considering only 1->2 branching
Lorentz5Momentum p, n, pthis, pcm;
assert(perturbative()!=2);
if(perturbative()==1) {
ShowerBasisPtr newBasis;
// find the partner and its momentum
if(!partner()) return;
if(partner()->isFinalState()) {
Lorentz5Momentum pa = -momentum()+partner()->momentum();
Lorentz5Momentum pb = momentum();
Energy scale=parent->momentum().t();
Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
Axis axis(pa.vect().unit());
LorentzRotation rot;
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
if(axis.perp2()>1e-20) {
rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
rot.rotateX(Constants::pi);
}
if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
pb *= rot;
if(pb.perp2()/GeV2>1e-20) {
Boost trans = -1./pb.e()*pb.vect();
trans.setZ(0.);
rot.boost(trans);
}
pbasis *=rot;
rot.invert();
n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
}
else {
pcm = parent->momentum();
p = Lorentz5Momentum(ZERO, pcm.vect());
n = Lorentz5Momentum(ZERO, -pcm.vect());
}
newBasis = new_ptr(ShowerBasis());
newBasis->setBasis(p,n,ShowerBasis::BackToBack);
showerBasis(newBasis,false);
}
else {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(children()[0])->showerBasis(),true);
}
}
void ShowerParticle::initializeFinalState() {
// set the basis vectors
Lorentz5Momentum p,n;
if(perturbative()!=0) {
ShowerBasisPtr newBasis;
// find the partner() and its momentum
if(!partner()) return;
Lorentz5Momentum ppartner(partner()->momentum());
// momentum of the emitting particle
p = momentum();
Lorentz5Momentum pcm;
// if the partner is a final-state particle then the reference
// vector is along the partner in the rest frame of the pair
if(partner()->isFinalState()) {
Boost boost=(p + ppartner).findBoostToCM();
pcm = ppartner;
pcm.boost(boost);
n = Lorentz5Momentum(ZERO,pcm.vect());
n.boost( -boost);
}
else if(!partner()->isFinalState()) {
// if the partner is an initial-state particle then the reference
// vector is along the partner which should be massless
if(perturbative()==1) {
n = Lorentz5Momentum(ZERO,ppartner.vect());
}
// if the partner is an initial-state decaying particle then the reference
// vector is along the backwards direction in rest frame of decaying particle
else {
Boost boost=ppartner.findBoostToCM();
pcm = p;
pcm.boost(boost);
n = Lorentz5Momentum( ZERO, -pcm.vect());
n.boost( -boost);
}
}
newBasis = new_ptr(ShowerBasis());
newBasis->setBasis(p,n,ShowerBasis::BackToBack);
showerBasis(newBasis,false);
}
else if(initiatesTLS()) {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0]->children()[0])->showerBasis(),true);
}
else {
showerBasis(dynamic_ptr_cast<ShowerParticlePtr>(parents()[0])->showerBasis(),true);
}
}
void ShowerParticle::setShowerMomentum(bool timeLike) {
Energy m = this->mass() > ZERO ? this->mass() : this->data().mass();
// calculate the momentum of the assuming on-shell
Energy2 pt2 = sqr(this->showerParameters().pt);
double alpha = timeLike ? this->showerParameters().alpha : this->x();
double beta = 0.5*(sqr(m) + pt2 - sqr(alpha)*showerBasis()->pVector().m2())/(alpha*showerBasis()->p_dot_n());
Lorentz5Momentum porig=showerBasis()->sudakov2Momentum(alpha,beta,
this->showerParameters().ptx,
this->showerParameters().pty);
porig.setMass(m);
this->set5Momentum(porig);
}
diff --git a/Shower/Core/Base/ShowerParticle.h b/Shower/Core/Base/ShowerParticle.h
--- a/Shower/Core/Base/ShowerParticle.h
+++ b/Shower/Core/Base/ShowerParticle.h
@@ -1,530 +1,528 @@
// -*- C++ -*-
//
// ShowerParticle.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_ShowerParticle_H
#define HERWIG_ShowerParticle_H
//
// This is the declaration of the ShowerParticle class.
//
#include "ThePEG/EventRecord/Particle.h"
#include "Herwig/Shower/Core/SplittingFunctions/SplittingFunction.fh"
#include "Herwig/Shower/Core/ShowerConfig.h"
#include "ShowerBasis.h"
#include "ShowerKinematics.h"
#include "ShowerParticle.fh"
#include <iosfwd>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* This class represents a particle in the showering process.
* It inherits from the Particle class of ThePEG and has some
* specifics information useful only during the showering process.
*
* Notice that:
* - for forward evolution, it is clear what is meant by parent/child;
* for backward evolution, however, it depends whether we want
* to keep a physical picture or a Monte-Carlo effective one.
* In the former case, an incoming particle (emitting particle)
* splits into an emitted particle and the emitting particle after
* the emission: the latter two are then children of the
* emitting particle, the parent. In the Monte-Carlo effective
* picture, we have that the particle close to the hard subprocess,
* with higher (space-like) virtuality, splits into an emitted particle
* and the emitting particle at lower virtuality: the latter two are,
* in this case, the children of the first one, the parent. However we
* choose a more physical picture where the new emitting particle is the
* parented of the emitted final-state particle and the original emitting
* particle.
* - the pointer to a SplitFun object is set only in the case
* that the particle has undergone a shower emission. This is similar to
* the case of the decay of a normal Particle where
* the pointer to a Decayer object is set only in the case
* that the particle has undergone to a decay.
* In the case of particle connected directly to the hard subprocess,
* there is no pointer to the hard subprocess, but there is a method
* isFromHardSubprocess() which returns true only in this case.
*
* @see Particle
* @see ShowerConfig
* @see ShowerKinematics
*/
class ShowerParticle: public Particle {
public:
/**
* Struct for all the info on an evolution partner
*/
struct EvolutionPartner {
/**
* Constructor
*/
EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType t,
Energy s) : partner(p), weight(w), type(t), scale(s)
{}
/**
* The partner
*/
tShowerParticlePtr partner;
/**
* Weight
*/
double weight;
/**
* Type
*/
ShowerPartnerType type;
/**
* The assoicated evolution scale
*/
Energy scale;
};
/**
* Struct to store the evolution scales
*/
struct EvolutionScales {
/**
* Constructor
*/
EvolutionScales() : QED(),QCD_c(),QCD_ac(),
QED_noAO(),QCD_c_noAO(),QCD_ac_noAO(),
Max_Q2(Constants::MaxEnergy2)
{}
/**
* QED scale
*/
Energy QED;
/**
* QCD colour scale
*/
Energy QCD_c;
/**
* QCD anticolour scale
*/
Energy QCD_ac;
/**
* QED scale
*/
Energy QED_noAO;
/**
* QCD colour scale
*/
Energy QCD_c_noAO;
/**
* QCD anticolour scale
*/
Energy QCD_ac_noAO;
/**
* Maximum allowed virtuality of the particle
*/
Energy2 Max_Q2;
};
/** @name Construction and descruction functions. */
//@{
/**
* Standard Constructor. Note that the default constructor is
* private - there is no particle without a pointer to a
* ParticleData object.
* @param x the ParticleData object
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase() {}
/**
* Copy constructor from a ThePEG Particle
* @param x ThePEG particle
* @param pert Where the particle came from
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase(&x) {}
//@}
public:
/**
* Set a preliminary momentum for the particle
*/
void setShowerMomentum(bool timelike);
/**
* Construct the spin info object for a shower particle
*/
void constructSpinInfo(bool timelike);
/**
* Perform any initial calculations needed after the branching has been selected
*/
void initializeDecay();
/**
* Perform any initial calculations needed after the branching has been selected
* @param parent The beam particle
*/
void initializeInitialState(PPtr parent);
/**
* Perform any initial calculations needed after the branching has been selected
*/
void initializeFinalState();
/**
* Access/Set various flags about the state of the particle
*/
//@{
/**
* Access the flag that tells if the particle is final state
* or initial state.
*/
bool isFinalState() const { return _isFinalState; }
/**
* Access the flag that tells if the particle is initiating a
* time like shower when it has been emitted in an initial state shower.
*/
bool initiatesTLS() const { return _initiatesTLS; }
/**
* Access the flag which tells us where the particle came from
* This is 0 for a particle produced in the shower, 1 if the particle came
* from the hard sub-process and 2 is it came from a decay.
*/
unsigned int perturbative() const { return _perturbative; }
//@}
/**
* Set/Get the momentum fraction for initial-state particles
*/
//@{
/**
* For an initial state particle get the fraction of the beam momentum
*/
void x(double x) { _x = x; }
/**
* For an initial state particle set the fraction of the beam momentum
*/
double x() const { return _x; }
//@}
/**
* Set/Get methods for the ShowerKinematics objects
*/
//@{
/**
* Access/ the ShowerKinematics object.
*/
const ShoKinPtr & showerKinematics() const { return _showerKinematics; }
/**
* Set the ShowerKinematics object.
*/
void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; }
//@}
/**
* Set/Get methods for the ShowerBasis objects
*/
//@{
/**
* Access/ the ShowerBasis object.
*/
const ShowerBasisPtr & showerBasis() const { return _showerBasis; }
/**
* Set the ShowerBasis object.
*/
void showerBasis(const ShowerBasisPtr in, bool copy) {
if(!copy)
_showerBasis = in;
else {
_showerBasis = new_ptr(ShowerBasis());
_showerBasis->setBasis(in->pVector(),in->nVector(),in->frame());
}
}
//@}
/**
* Members relating to the initial evolution scale and partner for the particle
*/
//@{
/**
* Veto emission at a given scale
*/
void vetoEmission(ShowerPartnerType type, Energy scale);
/**
* Access to the evolution scales
*/
const EvolutionScales & scales() const {return scales_;}
/**
* Access to the evolution scales
*/
EvolutionScales & scales() {return scales_;}
/**
* Return the virtual mass\f$
*/
Energy virtualMass() const { return _vMass; }
/**
* Set the virtual mass
*/
void virtualMass(Energy mass) { _vMass = mass; }
/**
* Return the partner
*/
tShowerParticlePtr partner() const { return _partner; }
/**
* Set the partner
*/
void partner(const tShowerParticlePtr partner) { _partner = partner; }
/**
* Get the possible partners
*/
vector<EvolutionPartner> & partners() { return partners_; }
/**
* Add a possible partners
*/
void addPartner(EvolutionPartner in );
/**
* Clear the evolution partners
*/
void clearPartners() { partners_.clear(); }
/**
* Return the progenitor of the shower
*/
tShowerParticlePtr progenitor() const { return _progenitor; }
/**
* Set the progenitor of the shower
*/
void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; }
//@}
/**
* Members to store and provide access to variables for a specific
* shower evolution scheme
*/
//@{
struct Parameters {
Parameters() : alpha(1.), beta(), ptx(), pty(), pt() {}
double alpha;
double beta;
Energy ptx;
Energy pty;
Energy pt;
};
/**
* Set the vector containing dimensionless variables
*/
Parameters & showerParameters() { return _parameters; }
//@}
/**
* If this particle came from the hard process get a pointer to ThePEG particle
* it came from
*/
const tcPPtr thePEGBase() const { return _thePEGBase; }
public:
/**
* Extract the rho matrix including mapping needed in the shower
*/
RhoDMatrix extractRhoMatrix(bool forward);
-protected:
-
/**
* For a particle which came from the hard process get the spin density and
* the mapping required to the basis used in the Shower
* @param rho The \f$\rho\f$ matrix
* @param mapping The mapping
* @param showerkin The ShowerKinematics object
*/
bool getMapping(SpinPtr &, RhoDMatrix & map);
protected:
/**
* Standard clone function.
*/
virtual PPtr clone() const;
/**
* Standard clone function.
*/
virtual PPtr fullclone() 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<ShowerParticle> initShowerParticle;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerParticle & operator=(const ShowerParticle &);
private:
/**
* Whether the particle is in the final or initial state
*/
bool _isFinalState;
/**
* Whether the particle came from
*/
unsigned int _perturbative;
/**
* Does a particle produced in the backward shower initiate a time-like shower
*/
bool _initiatesTLS;
/**
* Dimensionless parameters
*/
Parameters _parameters;
/**
* The beam energy fraction for particle's in the initial state
*/
double _x;
/**
* The shower kinematics for the particle
*/
ShoKinPtr _showerKinematics;
/**
* The shower basis for the particle
*/
ShowerBasisPtr _showerBasis;
/**
* Storage of the evolution scales
*/
EvolutionScales scales_;
/**
* Virtual mass
*/
Energy _vMass;
/**
* Partners
*/
tShowerParticlePtr _partner;
/**
* Pointer to ThePEG Particle this ShowerParticle was created from
*/
const tcPPtr _thePEGBase;
/**
* Progenitor
*/
tShowerParticlePtr _progenitor;
/**
* Partners
*/
vector<EvolutionPartner> partners_;
};
inline ostream & operator<<(ostream & os, const ShowerParticle::EvolutionScales & es) {
os << "Scales: QED=" << es.QED / GeV
<< " QCD_c=" << es.QCD_c / GeV
<< " QCD_ac=" << es.QCD_ac / GeV
<< " QED_noAO=" << es.QED_noAO / GeV
<< " QCD_c_noAO=" << es.QCD_c_noAO / GeV
<< " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV
<< '\n';
return os;
}
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ShowerParticle. */
template <>
struct BaseClassTrait<Herwig::ShowerParticle,1> {
/** Typedef of the first base class of ShowerParticle. */
typedef Particle NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ShowerParticle class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ShowerParticle>
: public ClassTraitsBase<Herwig::ShowerParticle> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ShowerParticle"; }
/** Create a Event object. */
static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); }
};
/** @endcond */
}
#endif /* HERWIG_ShowerParticle_H */
diff --git a/Shower/Core/Base/ShowerVertex.cc b/Shower/Core/Base/ShowerVertex.cc
--- a/Shower/Core/Base/ShowerVertex.cc
+++ b/Shower/Core/Base/ShowerVertex.cc
@@ -1,76 +1,95 @@
// -*- C++ -*-
//
// ShowerVertex.cc 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, non-templated member
// functions of the ShowerVertex class.
//
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ShowerVertex.h"
using namespace Herwig;
using namespace Herwig::Helicity;
using namespace ThePEG;
DescribeNoPIOClass<ShowerVertex,HelicityVertex>
describeShowerVertex ("Herwig::ShowerVertex","");
void ShowerVertex::Init() {
static ClassDocumentation<ShowerVertex> documentation
("The ShowerVertex class is the implementation of a "
"vertex for a shower for the Herwig spin correlation algorithm");
}
// method to get the rho matrix for a given outgoing particle
RhoDMatrix ShowerVertex::getRhoMatrix(int i, bool) const {
assert(matrixElement_->nOut()==2);
// calculate the incoming spin density matrix
RhoDMatrix input=incoming()[0]->rhoMatrix();
if(convertIn_) input = mapIncoming(input);
// get the rho matrices for the outgoing particles
vector<RhoDMatrix> rhoout;
for(unsigned int ix=0,N=outgoing().size();ix<N;++ix) {
if(int(ix)!=i)
rhoout.push_back(outgoing()[ix]->DMatrix());
}
// calculate the spin density matrix
return matrixElement_->calculateRhoMatrix(i,input,rhoout);
}
// method to get the D matrix for an incoming particle
RhoDMatrix ShowerVertex::getDMatrix(int) const {
assert(matrixElement_->nOut()==2);
// get the decay matrices for the outgoing particles
vector<RhoDMatrix> Dout;
for(unsigned int ix=0,N=outgoing().size();ix<N;++ix) {
Dout.push_back(outgoing()[ix]->DMatrix());
}
- // calculate the spin density matrix and return the answer
- return matrixElement_->calculateDMatrix(Dout);
+ // calculate the spin density matrix
+ RhoDMatrix rho = matrixElement_->calculateDMatrix(Dout);
+ // map if needed
+ if(convertIn_) {
+ RhoDMatrix rhop(rho.iSpin(),false);
+ for(int ixb=0;ixb<rho.iSpin();++ixb) {
+ for(int iyb=0;iyb<rho.iSpin();++iyb) {
+ if(inMatrix_(iyb,ixb)==0.)continue;
+ for(int iya=0;iya<rho.iSpin();++iya) {
+ if(rho(iya,iyb)==0.)continue;
+ for(int ixa=0;ixa<rho.iSpin();++ixa) {
+ rhop(ixa,ixb) += rho(iya,iyb)*inMatrix_(ixa,iya)*conj(inMatrix_(ixb,iyb));
+ }
+ }
+ }
+ }
+ rhop.normalize();
+ rho = rhop;
+ }
+ // return the answer
+ return rho;
}
RhoDMatrix ShowerVertex::mapIncoming(RhoDMatrix rho) const {
RhoDMatrix output(rho.iSpin());
for(int ixa=0;ixa<rho.iSpin();++ixa) {
for(int ixb=0;ixb<rho.iSpin();++ixb) {
for(int iya=0;iya<rho.iSpin();++iya) {
for(int iyb=0;iyb<rho.iSpin();++iyb) {
output(ixa,ixb) += rho(iya,iyb)*inMatrix_(iya,ixa)*conj(inMatrix_(iyb,ixb));
}
}
}
}
output.normalize();
return output;
}
diff --git a/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc b/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc
--- a/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc
+++ b/Shower/QTilde/Default/FS_QTildeShowerKinematics1to2.cc
@@ -1,200 +1,204 @@
// -*- C++ -*-
//
// FS_QTildeShowerKinematics1to2.cc 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, non-templated member
// functions of the FS_QTildeShowerKinematics1to2 class.
//
#include "FS_QTildeShowerKinematics1to2.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig/Shower/Core/SplittingFunctions/SplittingFunction.h"
#include "Herwig/Shower/Core/Base/ShowerParticle.h"
#include "ThePEG/Utilities/Debug.h"
#include "Herwig/Shower/QTilde/QTildeShowerHandler.h"
#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
#include "Herwig/Shower/QTilde/Base/ShowerModel.h"
#include "Herwig/Shower/QTilde/Base/KinematicsReconstructor.h"
#include "Herwig/Shower/Core/Base/ShowerVertex.h"
using namespace Herwig;
void FS_QTildeShowerKinematics1to2::
updateParameters(tShowerParticlePtr theParent,
tShowerParticlePtr theChild0,
tShowerParticlePtr theChild1,
bool setAlpha) const {
const ShowerParticle::Parameters & parent = theParent->showerParameters();
ShowerParticle::Parameters & child0 = theChild0->showerParameters();
ShowerParticle::Parameters & child1 = theChild1->showerParameters();
// determine alphas of children according to interpretation of z
if ( setAlpha ) {
child0.alpha = z() * parent.alpha;
child1.alpha = (1.-z()) * parent.alpha;
}
// set the values
double cphi = cos(phi());
double sphi = sin(phi());
child0.ptx = pT() * cphi + z() * parent.ptx;
child0.pty = pT() * sphi + z() * parent.pty;
child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) );
child1.ptx = -pT() * cphi + (1.-z())* parent.ptx;
child1.pty = -pT() * sphi + (1.-z())* parent.pty;
child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );
}
void FS_QTildeShowerKinematics1to2::
updateChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType partnerType,
bool massVeto) const {
assert(children.size()==2);
// calculate the scales
splittingFn()->evaluateFinalStateScales(partnerType,scale(),z(),parent,
children[0],children[1]);
// set the maximum virtual masses
if(massVeto) {
Energy2 q2 = z()*(1.-z())*sqr(scale());
IdList ids(3);
ids[0] = parent->dataPtr();
ids[1] = children[0]->dataPtr();
ids[2] = children[1]->dataPtr();
const vector<Energy> & virtualMasses = SudakovFormFactor()->virtualMasses(ids);
if(ids[0]->id()!=ParticleID::g && ids[0]->id()!=ParticleID::gamma ) {
q2 += sqr(virtualMasses[0]);
}
// limits on further evolution
children[0]->scales().Max_Q2 = z() *(q2-sqr(virtualMasses[2])/(1.-z()));
children[1]->scales().Max_Q2 = (1.-z())*(q2-sqr(virtualMasses[1])/ z() );
}
// update the parameters
updateParameters(parent, children[0], children[1], true);
// set up the colour connections
splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false);
// make the products children of the parent
parent->addChild(children[0]);
parent->addChild(children[1]);
// set the momenta of the children
for(ShowerParticleVector::const_iterator pit=children.begin();
pit!=children.end();++pit) {
(**pit).showerBasis(parent->showerBasis(),true);
(**pit).setShowerMomentum(true);
}
// sort out the helicity stuff
if(! dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->correlations()) return;
SpinPtr pspin(parent->spinInfo());
if(!pspin || !dynamic_ptr_cast<tcQTildeShowerHandlerPtr>(ShowerHandler::currentHandler())->spinCorrelations() ) return;
Energy2 t = sqr(scale())*z()*(1.-z());
IdList ids;
ids.push_back(parent->dataPtr());
ids.push_back(children[0]->dataPtr());
ids.push_back(children[1]->dataPtr());
// create the vertex
SVertexPtr vertex(new_ptr(ShowerVertex()));
// set the matrix element
vertex->ME(splittingFn()->matrixElement(z(),t,ids,phi(),true));
+ RhoDMatrix mapping;
+ SpinPtr inspin;
+ bool needMapping = parent->getMapping(inspin,mapping);
+ if(needMapping) vertex->incomingBasisTransform(mapping);
// set the incoming particle for the vertex
parent->spinInfo()->decayVertex(vertex);
for(ShowerParticleVector::const_iterator pit=children.begin();
pit!=children.end();++pit) {
// construct the spin info for the children
(**pit).constructSpinInfo(true);
// connect the spinInfo object to the vertex
(*pit)->spinInfo()->productionVertex(vertex);
}
}
void FS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr parent,
const ParticleVector & children ) const {
assert(children.size() == 2);
ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(children[0]);
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(children[1]);
parent->showerParameters().beta=
c1->showerParameters().beta + c2->showerParameters().beta;
Lorentz5Momentum pnew = c1->momentum() + c2->momentum();
Energy2 m2 = sqr(pT())/z()/(1.-z()) + sqr(c1->mass())/z()
+ sqr(c2->mass())/(1.-z());
pnew.setMass(sqrt(m2));
parent->set5Momentum( pnew );
}
void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr last,
Energy mass) const {
// set beta component and consequently all missing data from that,
// using the nominal (i.e. PDT) mass.
Energy theMass = mass > ZERO ? mass : last->data().constituentMass();
Lorentz5Momentum pVector = last->showerBasis()->pVector();
ShowerParticle::Parameters & lastParam = last->showerParameters();
Energy2 denom = 2. * lastParam.alpha * last->showerBasis()->p_dot_n();
if(abs(denom)/(sqr(pVector.e())+pVector.rho2())<1e-10) {
throw KinematicsReconstructionVeto();
}
lastParam.beta = ( sqr(theMass) + sqr(lastParam.pt)
- sqr(lastParam.alpha) * pVector.m2() )
/ denom;
// set that new momentum
Lorentz5Momentum newMomentum = last->showerBasis()->
sudakov2Momentum( lastParam.alpha, lastParam.beta,
lastParam.ptx , lastParam.pty);
newMomentum.setMass(theMass);
newMomentum.rescaleEnergy();
if(last->data().stable()) {
last->set5Momentum( newMomentum );
}
else {
last->boost(last->momentum().findBoostToCM());
last->boost(newMomentum.boostVector());
}
}
void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType) const {
IdList ids(3);
ids[0] = parent->dataPtr();
ids[1] = children[0]->dataPtr();
ids[2] = children[1]->dataPtr();
const vector<Energy> & virtualMasses = SudakovFormFactor()->virtualMasses(ids);
if(children[0]->children().empty()) children[0]->virtualMass(virtualMasses[1]);
if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]);
// compute the new pT of the branching
Energy2 pt2=sqr(z()*(1.-z()))*sqr(scale())
- sqr(children[0]->virtualMass())*(1.-z())
- sqr(children[1]->virtualMass())* z() ;
if(ids[0]->id()!=ParticleID::g) pt2 += z()*(1.-z())*sqr(virtualMasses[0]);
if(pt2>ZERO) {
pT(sqrt(pt2));
}
else {
pt2=ZERO;
pT(ZERO);
}
Energy2 q2 =
sqr(children[0]->virtualMass())/z() +
sqr(children[1]->virtualMass())/(1.-z()) +
pt2/z()/(1.-z());
parent->virtualMass(sqrt(q2));
}
void FS_QTildeShowerKinematics1to2::
resetChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children) const {
updateParameters(parent, children[0], children[1], false);
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->children().empty()) continue;
ShowerParticleVector newChildren;
for(unsigned int iy=0;iy<children[ix]->children().size();++iy)
newChildren.push_back(dynamic_ptr_cast<ShowerParticlePtr>
(children[ix]->children()[iy]));
children[ix]->showerKinematics()->resetChildren(children[ix],newChildren);
}
}
diff --git a/Shower/QTilde/Default/QTildeShowerKinematics1to2.cc b/Shower/QTilde/Default/QTildeShowerKinematics1to2.cc
--- a/Shower/QTilde/Default/QTildeShowerKinematics1to2.cc
+++ b/Shower/QTilde/Default/QTildeShowerKinematics1to2.cc
@@ -1,133 +1,133 @@
// -*- C++ -*-
//
// QTildeShowerKinematics1to2.cc 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, non-templated member
// functions of the QTildeShowerKinematics1to2 class.
//
#include "QTildeShowerKinematics1to2.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "Herwig/Shower/Core/Base/ShowerParticle.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/LorentzSpinorBar.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
vector<Lorentz5Momentum> QTildeShowerKinematics1to2::getBasis() const {
vector<Lorentz5Momentum> dum;
dum.push_back( _pVector );
dum.push_back( _nVector );
return dum;
}
void QTildeShowerKinematics1to2::setBasis(const Lorentz5Momentum &p,
const Lorentz5Momentum & n,
Frame inframe) {
_pVector=p;
_nVector=n;
frame(inframe);
Boost beta_bb;
if(frame()==BackToBack) {
beta_bb = -(_pVector + _nVector).boostVector();
}
else if(frame()==Rest) {
beta_bb = -pVector().boostVector();
}
else
assert(false);
Lorentz5Momentum p_bb = pVector();
Lorentz5Momentum n_bb = nVector();
p_bb.boost( beta_bb );
n_bb.boost( beta_bb );
// rotate to have z-axis parallel to p/n
Axis axis;
if(frame()==BackToBack) {
axis = p_bb.vect().unit();
}
else if(frame()==Rest) {
axis = n_bb.vect().unit();
}
else
assert(false);
LorentzRotation rot;
if(axis.perp2()>1e-10) {
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
else if(axis.z()<0.) {
rot.rotate(Constants::pi,Axis(1.,0.,0.));
}
_xPerp=LorentzVector<double>(1.,0.,0.,0.);
_yPerp=LorentzVector<double>(0.,1.,0.,0.);
_xPerp.transform(rot);
_yPerp.transform(rot);
// boost back
_xPerp.boost( -beta_bb );
_yPerp.boost( -beta_bb );
}
void QTildeShowerKinematics1to2::setMomentum(tShowerParticlePtr particle,
bool timeLike) const {
Energy mass = particle->mass() > ZERO ? particle->mass() : particle->data().mass();
// calculate the momentum of the assuming on-shell
Energy2 pt2 = sqr(particle->showerParameters().pt);
double alpha = timeLike ? particle->showerParameters().alpha : particle->x();
double beta = 0.5*(sqr(mass) + pt2 - sqr(alpha)*pVector().m2())/(alpha*p_dot_n());
Lorentz5Momentum porig=sudakov2Momentum(alpha,beta,
particle->showerParameters().ptx,
particle->showerParameters().pty);
porig.setMass(mass);
particle->set5Momentum(porig);
}
void QTildeShowerKinematics1to2::constructSpinInfo(tShowerParticlePtr particle,
bool timeLike) const {
// now construct the required spininfo and calculate the basis states
PDT::Spin spin(particle->dataPtr()->iSpin());
if(spin==PDT::Spin0) {
ScalarWaveFunction::constructSpinInfo(particle,outgoing,timeLike);
}
// calculate the basis states and construct the SpinInfo for a spin-1/2 particle
else if(spin==PDT::Spin1Half) {
// outgoing particle
if(particle->id()>0) {
vector<LorentzSpinorBar<SqrtEnergy> > stemp;
SpinorBarWaveFunction::calculateWaveFunctions(stemp,particle,outgoing);
SpinorBarWaveFunction::constructSpinInfo(stemp,particle,outgoing,timeLike);
}
// outgoing antiparticle
else {
vector<LorentzSpinor<SqrtEnergy> > stemp;
SpinorWaveFunction::calculateWaveFunctions(stemp,particle,outgoing);
SpinorWaveFunction::constructSpinInfo(stemp,particle,outgoing,timeLike);
}
}
// calculate the basis states and construct the SpinInfo for a spin-1 particle
else if(spin==PDT::Spin1) {
bool massless(particle->id()==ParticleID::g||particle->id()==ParticleID::gamma);
vector<Helicity::LorentzPolarizationVector> vtemp;
VectorWaveFunction::calculateWaveFunctions(vtemp,particle,outgoing,massless);
- VectorWaveFunction::constructSpinInfo(vtemp,particle,outgoing,timeLike,massless);
+ VectorWaveFunction::constructSpinInfo(vtemp,particle,outgoing,timeLike,massless,vector_phase);
}
else {
throw Exception() << "Spins higher than 1 are not yet implemented in "
<< "FS_QtildaShowerKinematics1to2::constructVertex() "
<< Exception::runerror;
}
}
void QTildeShowerKinematics1to2::transform(const LorentzRotation & r) {
_pVector *= r;
_nVector *= r;
_xPerp *= r;
_yPerp *= r;
}
diff --git a/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.cc
@@ -1,143 +1,143 @@
// -*- C++ -*-
//
// HalfOneHalfSplitFn.cc 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, non-templated member
// functions of the HalfOneHalfSplitFn class.
//
#include "HalfOneHalfSplitFn.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
using namespace Herwig;
DescribeNoPIOClass<HalfOneHalfSplitFn,Herwig::SplittingFunction>
describeHalfOneHalfSplitFn ("Herwig::HalfOneHalfSplitFn","HwShower.so");
void HalfOneHalfSplitFn::Init() {
static ClassDocumentation<HalfOneHalfSplitFn> documentation
("The HalfOneHalfSplitFn class implements the splitting "
"function for q -> g q");
}
double HalfOneHalfSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & ) const {
double val=(2.*(1.-z)+sqr(z))/z;
if(mass) {
Energy m = ids[0]->mass();
val-=2.*sqr(m)/t;
}
return colourFactor(ids)*val;
}
double HalfOneHalfSplitFn::overestimateP(const double z,
const IdList &ids) const {
return 2.*colourFactor(ids)/z;
}
double HalfOneHalfSplitFn::ratioP(const double z, const Energy2 t,
const IdList &ids,const bool mass, const RhoDMatrix & ) const {
double val=2.*(1.-z)+sqr(z);
if(mass) {
Energy m=ids[0]->mass();
val -=2.*sqr(m)*z/t;
}
return 0.5*val;
}
double HalfOneHalfSplitFn::integOverP(const double z, const IdList & ids,
unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return 2.*colourFactor(ids)*log(z);
case 1:
return -2.*colourFactor(ids)/z;
case 2:
return 2.*colourFactor(ids)*log(z/(1.-z));
case 3:
default:
throw Exception() << "HalfOneHalfSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double HalfOneHalfSplitFn::invIntegOverP(const double r,
const IdList & ids,
unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return exp(0.5*r/colourFactor(ids));
case 1:
return -2.*colourFactor(ids)/r;
case 2:
return 1./(1.+exp(-0.5*r/colourFactor(ids)));
case 3:
default:
throw Exception() << "HalfOneHalfSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool HalfOneHalfSplitFn::accept(const IdList &ids) const {
// 3 particles and in and out fermion same
if(ids.size()!=3 || ids[0]!=ids[2]) return false;
if(ids[0]->iSpin()!=PDT::Spin1Half ||
ids[1]->iSpin()!=PDT::Spin1) return false;
return checkColours(ids);
}
vector<pair<int, Complex> >
HalfOneHalfSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
vector<pair<int, Complex> >
HalfOneHalfSplitFn::generatePhiBackward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix & rho) {
assert(rho.iSpin()==PDT::Spin1);
double mt = sqr(ids[0]->mass())/t;
double diag = (1.+sqr(1.-z))/z - 2.*mt;
double off = 2.*(1.-z)/z*(1.-mt*z/(1.-z));
double max = diag+2.*abs(rho(0,2))*off;
vector<pair<int, Complex> > output;
output.push_back(make_pair( 0, (rho(0,0)+rho(2,2))*diag/max));
output.push_back(make_pair( 2, -rho(0,2) * off/max));
- output.push_back(make_pair(-2, -rho(2,0) * off/max));
+ output.push_back(make_pair(-2, -rho(2,0) * off/max));
return output;
}
DecayMEPtr HalfOneHalfSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool timeLike) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1,PDT::Spin1Half)));
Energy m = !timeLike ? ZERO : ids[0]->mass();
double mt = m/sqrt(t);
double root = sqrt(1.-z*sqr(m)/(1.-z)/t);
double romz = sqrt(1.-z);
double rz = sqrt(z);
- Complex phase = exp(Complex(0.,1.)*phi);
+ Complex phase = exp(-Complex(0.,1.)*phi);
(*kernal)(0,0,0) = -root/rz/phase;
(*kernal)(1,2,1) = -conj((*kernal)(0,0,0));
(*kernal)(0,2,0) = root/rz*(1.-z)*phase;
(*kernal)(1,0,1) = -conj((*kernal)(0,2,0));
(*kernal)(1,2,0) = mt*z/romz;
(*kernal)(0,0,1) = conj((*kernal)(1,2,0));
(*kernal)(0,2,1) = 0.;
(*kernal)(1,0,0) = 0.;
return kernal;
}
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 = ¶meters;
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;
}
}