Page MenuHomeHEPForge
Files F9501674

No OneTemporary
Actions

View Options

diff --git a/Decay/Dalitz/DalitzKMatrix.h b/Decay/Dalitz/DalitzKMatrix.h
--- a/Decay/Dalitz/DalitzKMatrix.h
+++ b/Decay/Dalitz/DalitzKMatrix.h
@@ -1,145 +1,145 @@
// -*- C++ -*-
#ifndef Herwig_DalitzKMatrix_H
#define Herwig_DalitzKMatrix_H
//
// This is the declaration of the DalitzKMatrix class.
//
#include "DalitzResonance.h"
#include "Herwig/Decay/FormFactors/KMatrix.h"
namespace Herwig {
using namespace ThePEG;
/**
* The DalitzKMatrix class allows the use of \f$K\f$-matrices in Dalitz decays
*/
class DalitzKMatrix: public DalitzResonance {
public:
/**
* The default constructor.
*/
DalitzKMatrix()
{}
/**
* Constructor specifiying the parameters
*/
DalitzKMatrix(long pid, ResonanceType::Type rtype, Energy m, Energy w,
unsigned int d1, unsigned int d2, unsigned int s,
double mag, double phi, InvEnergy rr,
unsigned int imat, unsigned int chan,
Energy2 sc, unsigned int itype,
vector<pair<double,double> > beta,
vector<pair<double,vector<double > > > coeffs)
: DalitzResonance(pid,rtype,m,w,d1,d2,s,mag,phi,rr),
imat_(imat), channel_(chan), sc_(sc), expType_(itype), coeffs_(coeffs) {
beta_.clear();
for(unsigned int ix=0;ix<beta.size();++ix) {
beta_.push_back(beta[ix].first*exp(Complex(0.,beta[ix].second)));
}
}
public:
/**
* Return the Breit-Wigner times the form factor
*/
virtual Complex BreitWigner(const Energy & mAB, const Energy & mA, const Energy & mB) const;
/**
* Output the parameters
*/
virtual void dataBaseOutput(ofstream & output);
public:
/**
* Set the K-matrix
*/
void setKMatrix(KMatrixPtr mat) {kMatrix_ = mat;}
/**
* Location of the matrix
*/
unsigned int imatrix() const {return imat_;}
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();
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DalitzKMatrix & operator=(const DalitzKMatrix &) = delete;
private:
/**
* The K-matrix for the channel
*/
KMatrixPtr kMatrix_;
/**
* Which \f$K-matrix\f$ to do
*/
- unsigned int imat_;
+ unsigned int imat_ = 0;
/**
* Which channel to use from the K-matrix
*/
- unsigned int channel_;
+ unsigned int channel_ = 0;
/**
* Expansion point for the constant terms
*/
- Energy2 sc_;
+ Energy2 sc_ = ZERO;
/**
* Coefficients of the poles
*/
vector<Complex> beta_;
/**
* Type of expansion
*/
- unsigned int expType_;
+ unsigned int expType_ = 0;
/**
* Coefficients for the series expansion
*/
vector<pair<double,vector<double > > > coeffs_;
};
}
#endif /* Herwig_DalitzKMatrix_H */
diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
--- a/Hadronization/ClusterFissioner.cc
+++ b/Hadronization/ClusterFissioner.cc
@@ -1,1677 +1,1677 @@
// -*- C++ -*-
//
// ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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.
//
//
// Thisk is the implementation of the non-inlined, non-templated member
// functions of the ClusterFissioner class.
//
#include "ClusterFissioner.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "Herwig/Utilities/Kinematics.h"
#include "Cluster.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include "ThePEG/Interface/ParMap.h"
#include "Herwig/Utilities/AlphaS.h"
using namespace Herwig;
DescribeClass<ClusterFissioner,Interfaced>
describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so");
ClusterFissioner::ClusterFissioner() :
_clMaxLight(3.35*GeV),
_clMaxDiquark(3.35*GeV),
_clMaxExotic(3.35*GeV),
_clPowLight(2.0),
_clPowDiquark(2.0),
_clPowExotic(2.0),
_pSplitLight(1.0),
_pSplitExotic(1.0),
_phaseSpaceWeights(0),
_dim(4),
_fissionCluster(0),
_kinematicThresholdChoice(0),
_pwtDIquark(0.0),
_diquarkClusterFission(0),
_btClM(1.0*GeV),
_iopRem(1),
_kappa(1.0e15*GeV/meter),
_enhanceSProb(0),
_m0Fission(2.*GeV),
_massMeasure(0),
_probPowFactor(4.0),
_probShift(0.0),
_kinThresholdShift(1.0*sqr(GeV)),
_strictDiquarkKinematics(0),
_covariantBoost(false),
_hadronizingStrangeDiquarks(2),
_writeOut(0)
{
}
ClusterFissioner::~ClusterFissioner(){
}
IBPtr ClusterFissioner::clone() const {
return new_ptr(*this);
}
IBPtr ClusterFissioner::fullclone() const {
return new_ptr(*this);
}
void ClusterFissioner::persistentOutput(PersistentOStream & os) const {
os << ounit(_clMaxLight,GeV) << ounit(_clMaxHeavy,GeV) << ounit(_clMaxDiquark,GeV) << ounit(_clMaxExotic,GeV)
<< _clPowLight << _clPowHeavy << _clPowDiquark << _clPowExotic
<< _pSplitLight << _pSplitHeavy << _pSplitExotic
<< _fissionCluster << _fissionPwt
<< _pwtDIquark
<< _diquarkClusterFission
<< ounit(_btClM,GeV)
<< _iopRem << ounit(_kappa, GeV/meter)
<< _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure
<< _dim << _phaseSpaceWeights
<< _hadronSpectrum << _kinematicThresholdChoice
<< _probPowFactor << _probShift << ounit(_kinThresholdShift,sqr(GeV))
<< _strictDiquarkKinematics
<< _covariantBoost
<< _hadronizingStrangeDiquarks
<< _writeOut
;
}
void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
is >> iunit(_clMaxLight,GeV) >> iunit(_clMaxHeavy,GeV) >> iunit(_clMaxDiquark,GeV) >> iunit(_clMaxExotic,GeV)
>> _clPowLight >> _clPowHeavy >> _clPowDiquark >> _clPowExotic
>> _pSplitLight >> _pSplitHeavy >> _pSplitExotic
>> _fissionCluster >> _fissionPwt
>> _pwtDIquark
>> _diquarkClusterFission
>> iunit(_btClM,GeV)
>> _iopRem >> iunit(_kappa, GeV/meter)
>> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure
>> _dim >> _phaseSpaceWeights
>> _hadronSpectrum >> _kinematicThresholdChoice
>> _probPowFactor >> _probShift >> iunit(_kinThresholdShift,sqr(GeV))
>> _strictDiquarkKinematics
>> _covariantBoost
>> _hadronizingStrangeDiquarks
>> _writeOut
;
}
void ClusterFissioner::doinit() {
Interfaced::doinit();
if (_writeOut){
std::ofstream out("data_CluFis.dat", std::ios::out);
out.close();
}
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( _pSplitHeavy.find(id) == _pSplitHeavy.end() ||
_clPowHeavy.find(id) == _clPowHeavy.end() ||
_clMaxHeavy.find(id) == _clMaxHeavy.end() ){
std::cout << "id = "<<id << std::endl;
throw InitException() << "not all parameters have been set for heavy quark cluster fission";
}
}
// for default Pwts not needed to initialize
if (_fissionCluster==0) return;
for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
if ( _fissionPwt.find(id) == _fissionPwt.end() )
// check that all relevant weights are set
throw InitException() << "fission weights for light quarks have not been set";
}
double pwtDquark=_fissionPwt.find(ParticleID::d)->second;
double pwtUquark=_fissionPwt.find(ParticleID::u)->second;
double pwtSquark=_fissionPwt.find(ParticleID::s)->second;
// TODO better solution for this magic number alternative
_fissionPwt[1103] = _pwtDIquark * pwtDquark * pwtDquark;
_fissionPwt[2101] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
_fissionPwt[2203] = _pwtDIquark * pwtUquark * pwtUquark;
if (_hadronizingStrangeDiquarks>0) {
_fissionPwt[3101] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
_fissionPwt[3201] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
if (_hadronizingStrangeDiquarks==2) {
_fissionPwt[3303] = _pwtDIquark* pwtSquark * pwtSquark;
}
}
}
void ClusterFissioner::Init() {
static ClassDocumentation<ClusterFissioner> documentation
("Class responsibles for chopping up the clusters");
static Reference<ClusterFissioner,HadronSpectrum> interfaceHadronSpectrum
("HadronSpectrum",
"Set the Hadron spectrum for this cluster fissioner.",
&ClusterFissioner::_hadronSpectrum, false, false, true, false);
// ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,Energy>
interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])",
&ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 100.0*GeV,
false,false,false);
static Parameter<ClusterFissioner,Energy>
interfaceClMaxDiquark ("ClMaxDiquark","cluster max mass for light hadronizing diquarks (unit [GeV])",
&ClusterFissioner::_clMaxDiquark, GeV, 3.35*GeV, ZERO, 100.0*GeV,
false,false,false);
static ParMap<ClusterFissioner,Energy> interfaceClMaxHeavy
("ClMaxHeavy",
"ClMax for heavy quarks",
&ClusterFissioner::_clMaxHeavy, GeV, -1, 3.35*GeV, ZERO, 100.0*GeV,
false, false, Interface::upperlim);
static Parameter<ClusterFissioner,Energy>
interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])",
&ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 100.0*GeV,
false,false,false);
// ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,double>
interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks",
&ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false);
static ParMap<ClusterFissioner,double> interfaceClPowHeavy
("ClPowHeavy",
"ClPow for heavy quarks",
&ClusterFissioner::_clPowHeavy, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
static Parameter<ClusterFissioner,double>
interfaceClPowDiquark ("ClPowDiquark","cluster mass exponent for light hadronizing diquarks",
&ClusterFissioner::_clPowDiquark, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks",
&ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false);
// PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,double>
interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks",
&ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false);
static ParMap<ClusterFissioner,double> interfacePSplitHeavy
("PSplitHeavy",
"PSplit for heavy quarks",
&ClusterFissioner::_pSplitHeavy, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
static Parameter<ClusterFissioner,double>
interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks",
&ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false);
static Switch<ClusterFissioner,int> interfaceFission
("Fission",
"Option for different Fission options",
&ClusterFissioner::_fissionCluster, 1, false, false);
static SwitchOption interfaceFissionDefault
(interfaceFission,
"Default",
"Normal cluster fission which depends on the hadron spectrum class.",
0);
static SwitchOption interfaceFissionNew
(interfaceFission,
"New",
"Alternative cluster fission which does not depend on the hadron spectrum class",
1);
static SwitchOption interfaceFissionNewDiquarkSuppression
(interfaceFission,
"NewDiquarkSuppression",
"Alternative cluster fission which does not depend on the hadron spectrum class"
" and includes a suppression of AlphaS^2(Mc) for Diquark Production during "
"Cluster Fission",
-1);
static Switch<ClusterFissioner,int> interfaceDiquarkClusterFission
("DiquarkClusterFission",
"Allow clusters to fission to 1 or 2 diquark Clusters or Turn off diquark fission completely",
&ClusterFissioner::_diquarkClusterFission, 0, false, false);
static SwitchOption interfaceDiquarkClusterFissionAll
(interfaceDiquarkClusterFission,
"All",
"Allow diquark clusters and baryon clusters to fission to new diquark Clusters",
2);
static SwitchOption interfaceDiquarkClusterFissionOnlyBaryonClusters
(interfaceDiquarkClusterFission,
"OnlyBaryonClusters",
"Allow only baryon clusters to fission to new diquark Clusters",
1);
static SwitchOption interfaceDiquarkClusterFissionNo
(interfaceDiquarkClusterFission,
"No",
"Don't allow clusters to fission to new diquark Clusters",
0);
static SwitchOption interfaceDiquarkClusterFissionOff
(interfaceDiquarkClusterFission,
"Off",
"Don't allow clusters fission to draw diquarks ",
-1);
static ParMap<ClusterFissioner,double> interfaceFissionPwt
("FissionPwt",
"The weights for quarks in the fission process.",
&ClusterFissioner::_fissionPwt, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
static Switch<ClusterFissioner,int> interfaceRemnantOption
("RemnantOption",
"Option for the treatment of remnant clusters",
&ClusterFissioner::_iopRem, 1, false, false);
static SwitchOption interfaceRemnantOptionSoft
(interfaceRemnantOption,
"Soft",
"Both clusters produced in the fission of the beam cluster"
" are treated as soft clusters.",
0);
static SwitchOption interfaceRemnantOptionHard
(interfaceRemnantOption,
"Hard",
"Only the cluster containing the remnant is treated as a soft cluster.",
1);
static SwitchOption interfaceRemnantOptionVeryHard
(interfaceRemnantOption,
"VeryHard",
"Even remnant clusters are treated as hard, i.e. all clusters the same",
2);
static Parameter<ClusterFissioner,Energy> interfaceBTCLM
("SoftClusterFactor",
"Parameter for the mass spectrum of remnant clusters",
&ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<ClusterFissioner,Tension> interfaceStringTension
("StringTension",
"String tension used in vertex displacement calculation",
&ClusterFissioner::_kappa, GeV/meter,
1.0e15*GeV/meter, ZERO, ZERO,
false, false, Interface::lowerlim);
static Switch<ClusterFissioner,int> interfaceEnhanceSProb
("EnhanceSProb",
"Option for enhancing strangeness",
&ClusterFissioner::_enhanceSProb, 0, false, false);
static SwitchOption interfaceEnhanceSProbNo
(interfaceEnhanceSProb,
"No",
"No strangeness enhancement.",
0);
static SwitchOption interfaceEnhanceSProbScaled
(interfaceEnhanceSProb,
"Scaled",
"Scaled strangeness enhancement",
1);
static SwitchOption interfaceEnhanceSProbExponential
(interfaceEnhanceSProb,
"Exponential",
"Exponential strangeness enhancement",
2);
static Switch<ClusterFissioner,int> interfaceMassMeasure
("MassMeasure",
"Option to use different mass measures",
&ClusterFissioner::_massMeasure,0,false,false);
static SwitchOption interfaceMassMeasureMass
(interfaceMassMeasure,
"Mass",
"Mass Measure",
0);
static SwitchOption interfaceMassMeasureLambda
(interfaceMassMeasure,
"Lambda",
"Lambda Measure",
1);
static Parameter<ClusterFissioner,Energy> interfaceFissionMassScale
("FissionMassScale",
"Cluster fission mass scale",
&ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV,
false, false, Interface::limited);
static Parameter<ClusterFissioner,double> interfaceProbPowFactor
("ProbabilityPowerFactor",
"Power factor in ClusterFissioner bell probablity function",
&ClusterFissioner::_probPowFactor, 2.0, 0.001, 20.0,
false, false, Interface::limited);
static Parameter<ClusterFissioner,double> interfaceProbShift
("ProbabilityShift",
"Shifts from the center in ClausterFissioner bell probablity function",
&ClusterFissioner::_probShift, 0.0, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<ClusterFissioner,Energy2> interfaceKineticThresholdShift
("KineticThresholdShift",
"Shifts from the kinetic threshold in ClusterFissioner",
&ClusterFissioner::_kinThresholdShift, sqr(GeV), 0.*sqr(GeV), -10.0*sqr(GeV), 10.0*sqr(GeV),
false, false, Interface::limited);
static Switch<ClusterFissioner,int> interfaceKinematicThreshold
("KinematicThreshold",
"Option for using static or dynamic kinematic thresholds in cluster splittings",
&ClusterFissioner::_kinematicThresholdChoice, 0, false, false);
static SwitchOption interfaceKinematicThresholdStatic
(interfaceKinematicThreshold,
"Static",
"Set static kinematic thresholds for cluster splittings.",
0);
static SwitchOption interfaceKinematicThresholdDynamic
(interfaceKinematicThreshold,
"Dynamic",
"Set dynamic kinematic thresholds for cluster splittings.",
1);
static Switch<ClusterFissioner,bool> interfaceCovariantBoost
("CovariantBoost",
"Use single Covariant Boost for Cluster Fission",
&ClusterFissioner::_covariantBoost, false, false, false);
static SwitchOption interfaceCovariantBoostYes
(interfaceCovariantBoost,
"Yes",
"Use Covariant boost",
true);
static SwitchOption interfaceCovariantBoostNo
(interfaceCovariantBoost,
"No",
"Do NOT use Covariant boost",
false);
static Switch<ClusterFissioner,int> interfaceStrictDiquarkKinematics
("StrictDiquarkKinematics",
"Option for selecting different selection criterions of diquarks for ClusterFission",
&ClusterFissioner::_strictDiquarkKinematics, 0, false, false);
static SwitchOption interfaceStrictDiquarkKinematicsLoose
(interfaceStrictDiquarkKinematics,
"Loose",
"No kinematic threshold for diquark selection except for Mass bigger than 2 baryons",
0);
static SwitchOption interfaceStrictDiquarkKinematicsStrict
(interfaceStrictDiquarkKinematics,
"Strict",
"Resulting clusters are at least as heavy as 2 lightest baryons",
1);
static Parameter<ClusterFissioner,double> interfacePwtDIquark
("PwtDIquark",
"specific probability for choosing a d diquark",
&ClusterFissioner::_pwtDIquark, 0.0, 0.0, 10.0,
false, false, Interface::limited);
static Switch<ClusterFissioner,int> interfacePhaseSpaceWeights
("PhaseSpaceWeights",
"Include phase space weights.",
&ClusterFissioner::_phaseSpaceWeights, 0, false, false);
static SwitchOption interfacePhaseSpaceWeightsNo
(interfacePhaseSpaceWeights,
"No",
"Do not include the effect of cluster phase space",
0);
static SwitchOption interfacePhaseSpaceWeightsYes
(interfacePhaseSpaceWeights,
"Yes",
"Do include the effect of cluster fission phase space "
"related to constituent masses."
"Note: Need static Threshold choice",
1);
static SwitchOption interfacePhaseSpaceWeightsUseHadronMasses
(interfacePhaseSpaceWeights,
"UseHadronMasses",
"Do include the effect of cluster fission phase space "
"related to hadron masses."
"Note: Need static Threshold choice",
2);
static SwitchOption interfacePhaseSpaceWeightsNoConstituentMasses
(interfacePhaseSpaceWeights,
"NoConstituentMasses",
"Do not include the effect of cluster fission phase space "
"related to constituent masses."
"Note: Need static Threshold choice",
3);
static Parameter<ClusterFissioner,double>
interfaceDim ("Dimension","Dimension in which phase space weights are calculated",
&ClusterFissioner::_dim, 0, 4.0, 0.0, 10.0,false,false,false);
// Allowing for strange diquarks in the ClusterFission
static Switch<ClusterFissioner,unsigned int> interfaceHadronizingStrangeDiquarks
("HadronizingStrangeDiquarks",
"Option for adding strange diquarks to Cluster Fission (if Fission = New or Hybrid is enabled)",
&ClusterFissioner::_hadronizingStrangeDiquarks, 0, false, false);
static SwitchOption interfaceHadronizingStrangeDiquarksNo
(interfaceHadronizingStrangeDiquarks,
"No",
"No strangeness containing diquarks during Cluster Fission",
0);
static SwitchOption interfaceHadronizingStrangeDiquarksOnlySingleStrange
(interfaceHadronizingStrangeDiquarks,
"OnlySingleStrange",
"Only one strangeness containing diquarks during Cluster Fission i.e. su,sd",
1);
static SwitchOption interfaceHadronizingStrangeDiquarksAll
(interfaceHadronizingStrangeDiquarks,
"All",
"All strangeness containing diquarks during Cluster Fission i.e. su,sd,ss",
2);
}
tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) {
// return if no clusters
if (clusters.empty()) return tPVector();
/*****************
* Loop over the (input) collection of cluster pointers, and store in
* the vector splitClusters all the clusters that need to be split
* (these are beam clusters, if soft underlying event is off, and
* heavy non-beam clusters).
********************/
stack<ClusterPtr> splitClusters;
for(ClusterVector::iterator it = clusters.begin() ;
it != clusters.end() ; ++it) {
/**************
* Skip 3-component clusters that have been redefined (as 2-component
* clusters) or not available clusters. The latter check is indeed
* redundant now, but it is used for possible future extensions in which,
* for some reasons, some of the clusters found by ClusterFinder are tagged
* straight away as not available.
**************/
if((*it)->isRedefined() || !(*it)->isAvailable()) continue;
// if the cluster is a beam cluster add it to the vector of clusters
// to be split or if it is heavy
if((*it)->isBeamCluster() || isHeavy(*it)) splitClusters.push(*it);
}
tPVector finalhadrons;
cut(splitClusters, clusters, finalhadrons, softUEisOn);
return finalhadrons;
}
void ClusterFissioner::cut(stack<ClusterPtr> & clusterStack,
ClusterVector &clusters, tPVector & finalhadrons,
bool softUEisOn) {
/**************************************************
* This method does the splitting of the cluster pointed by cluPtr
* and "recursively" by all of its cluster children, if heavy. All of these
* new children clusters are added (indeed the pointers to them) to the
* collection of cluster pointers collecCluPtr. The method works as follows.
* Initially the vector vecCluPtr contains just the input pointer to the
* cluster to be split. Then it will be filled "recursively" by all
* of the cluster's children that are heavy enough to require, in their turn,
* to be split. In each loop, the last element of the vector vecCluPtr is
* considered (only once because it is then removed from the vector).
* This approach is conceptually recursive, but avoid the overhead of
* a concrete recursive function. Furthermore it requires minimal changes
* in the case that the fission of an heavy cluster could produce more
* than two cluster children as assumed now.
*
* Draw the masses: for normal, non-beam clusters a power-like mass dist
* is used, whereas for beam clusters a fast-decreasing exponential mass
* dist is used instead (to avoid many iterative splitting which could
* produce an unphysical large transverse energy from a supposed soft beam
* remnant process).
****************************************/
// Here we recursively loop over clusters in the stack and cut them
while (!clusterStack.empty()) {
// take the last element of the vector
ClusterPtr iCluster = clusterStack.top(); clusterStack.pop();
// split it
cutType ct = iCluster->numComponents() == 2 ?
cutTwo(iCluster, finalhadrons, softUEisOn) :
cutThree(iCluster, finalhadrons, softUEisOn);
// There are cases when we don't want to split, even if it fails mass test
if(!ct.first.first || !ct.second.first) {
// if an unsplit beam cluster leave if for the underlying event
if(iCluster->isBeamCluster() && softUEisOn)
iCluster->isAvailable(false);
continue;
}
// check if clusters
ClusterPtr one = dynamic_ptr_cast<ClusterPtr>(ct.first.first);
ClusterPtr two = dynamic_ptr_cast<ClusterPtr>(ct.second.first);
// is a beam cluster must be split into two clusters
if(iCluster->isBeamCluster() && (!one||!two) && softUEisOn) {
iCluster->isAvailable(false);
continue;
}
// There should always be a intermediate quark(s) from the splitting
assert(ct.first.second && ct.second.second);
/// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors
iCluster->addChild(ct.first.first);
// iCluster->addChild(ct.first.second);
// ct.first.second->addChild(ct.first.first);
iCluster->addChild(ct.second.first);
// iCluster->addChild(ct.second.second);
// ct.second.second->addChild(ct.second.first);
// Sometimes the clusters decay C -> H + C' or C -> H + H' rather then C -> C' + C''
if(one) {
clusters.push_back(one);
if(one->isBeamCluster() && softUEisOn)
one->isAvailable(false);
if(isHeavy(one) && one->isAvailable())
clusterStack.push(one);
}
if(two) {
clusters.push_back(two);
if(two->isBeamCluster() && softUEisOn)
two->isAvailable(false);
if(isHeavy(two) && two->isAvailable())
clusterStack.push(two);
}
}
}
ClusterFissioner::cutType
ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons,
bool softUEisOn) {
// need to make sure only 2-cpt clusters get here
assert(cluster->numComponents() == 2);
tPPtr ptrQ1 = cluster->particle(0);
tPPtr ptrQ2 = cluster->particle(1);
Energy Mc = cluster->mass();
assert(ptrQ1);
assert(ptrQ2);
// And check if those particles are from a beam remnant
bool rem1 = cluster->isBeamRemnant(0);
bool rem2 = cluster->isBeamRemnant(1);
// workout which distribution to use
bool soft1(false),soft2(false);
switch (_iopRem) {
case 0:
soft1 = rem1 || rem2;
soft2 = rem2 || rem1;
break;
case 1:
soft1 = rem1;
soft2 = rem2;
break;
}
// Initialization for the exponential ("soft") mass distribution.
static const int max_loop = 1000;
int counter = 0;
Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO;
tcPDPtr toHadron1, toHadron2;
PPtr newPtr1 = PPtr ();
PPtr newPtr2 = PPtr ();
bool succeeded = false;
Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
do
{
succeeded = false;
++counter;
// get a flavour for the qqbar pair
drawNewFlavour(newPtr1,newPtr2,cluster);
// check for right ordering
assert (ptrQ2);
assert (newPtr2);
assert (ptrQ2->dataPtr());
assert (newPtr2->dataPtr());
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
swap(newPtr1, newPtr2);
// check again
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
throw Exception()
<< "ClusterFissioner cannot split the cluster ("
<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
<< ") into hadrons.\n" << Exception::runerror;
}
}
// Check that new clusters can produce particles and there is enough
// phase space to choose the drawn flavour
m1 = ptrQ1->data().constituentMass();
m2 = ptrQ2->data().constituentMass();
m = newPtr1->data().constituentMass();
// Do not split in the case there is no phase space available
if(Mc < m1+m + m2+m) continue;
pQ1.setMass(m1);
pQone.setMass(m);
pQtwo.setMass(m);
pQ2.setMass(m2);
double weightMasses = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
ptrQ1, pQ1, newPtr1, pQone,
newPtr2, pQtwo, ptrQ2, pQ2);
if (weightMasses==0.0)
continue;
// derive the masses of the children
Mc1 = pClu1.mass();
Mc2 = pClu2.mass();
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
if (Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue;
if (_phaseSpaceWeights==2 &&
( Mc1 < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),newPtr1->dataPtr())
|| Mc2 < spectrum()->massLightestHadronPair(ptrQ2->dataPtr(),newPtr2->dataPtr()) ))
continue;
// dynamic kinematic threshold
}
else if(_kinematicThresholdChoice == 1) {
bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false;
if( C1 || C2 || C3 ) continue;
}
if ( _phaseSpaceWeights && phaseSpaceVeto(Mc,Mc1,Mc2,m,m1,m2, ptrQ1, ptrQ2, newPtr1, 0.0) ) {
// reduce counter as it regards only the mass sampling
counter--;
continue;
}
/**************************
* New (not present in Fortran Herwig):
* check whether the fragment masses Mc1 and Mc2 are above the
* threshold for the production of the lightest pair of hadrons with the
* right flavours. If not, then set by hand the mass to the lightest
* single hadron with the right flavours, in order to solve correctly
* the kinematics, and (later in this method) create directly such hadron
* and add it to the children hadrons of the cluster that undergoes the
* fission (i.e. the one pointed by iCluPtr). Notice that in this special
* case, the heavy cluster that undergoes the fission has one single
* cluster child and one single hadron child. We prefer this approach,
* rather than to create a light cluster, with the mass set equal to
* the lightest hadron, and let then the class LightClusterDecayer to do
* the job to decay it to that single hadron, for two reasons:
* First, because the sum of the masses of the two constituents can be,
* in this case, greater than the mass of that hadron, hence it would
* be impossible to solve the kinematics for such two components, and
* therefore we would have a cluster whose components are undefined.
* Second, the algorithm is faster, because it avoids the reshuffling
* procedure that would be necessary if we used LightClusterDecayer
* to decay the light cluster to the lightest hadron.
****************************/
// override chosen masses if needed
toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
// if a beam cluster not allowed to decay to hadrons
if(cluster->isBeamCluster() && (toHadron1||toHadron2) && softUEisOn)
continue;
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > Mc) {
if(!toHadron1) {
toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
}
else if(!toHadron2) {
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
}
}
succeeded = (Mc >= Mc1+Mc2);
}
while (!succeeded && counter < max_loop);
if(counter >= max_loop) {
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
// Determined the (5-components) momenta (all in the LAB frame)
Lorentz5Momentum pClu = cluster->momentum(); // known
Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
/******************
* The previous methods have determined the kinematics and positions
* of C -> C1 + C2.
* In the case that one of the two product is light, that means either
* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
* of the components of that light product have not been determined,
* and a (light) cluster will not be created: the heavy father cluster
* decays, in this case, into a single (not-light) cluster and a
* single hadron. In the other, "normal", cases the father cluster
* decays into two clusters, each of which has well defined components.
* Notice that, in the case of components which point to particles, the
* momenta of the components is properly set to the new values, whereas
* we do not change the momenta of the pointed particles, because we
* want to keep all of the information (that is the new momentum of a
* component after the splitting, which is contained in the _momentum
* member of the Component class, and the (old) momentum of that component
* before the splitting, which is contained in the momentum of the
* pointed particle). Please not make confusion of this only apparent
* inconsistency!
********************/
LorentzPoint posC,pos1,pos2;
posC = cluster->vertex();
calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
cutType rval;
if(toHadron1) {
rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toHadron2) {
rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
finalhadrons.push_back(rval.second.first);
}
else {
rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
}
return rval;
}
ClusterFissioner::cutType
ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons,
bool softUEisOn) {
// need to make sure only 3-cpt clusters get here
assert(cluster->numComponents() == 3);
// extract quarks
tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)};
assert( ptrQ[0] && ptrQ[1] && ptrQ[2] );
// find maximum mass pair
Energy mmax(ZERO);
Lorentz5Momentum pDiQuark;
int iq1(-1),iq2(-1);
Lorentz5Momentum psum;
for(int q1=0;q1<3;++q1) {
psum+= ptrQ[q1]->momentum();
for(int q2=q1+1;q2<3;++q2) {
Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum();
ptest.rescaleMass();
Energy mass = ptest.m();
if(mass>mmax) {
mmax = mass;
pDiQuark = ptest;
iq1 = q1;
iq2 = q2;
}
}
}
// and the spectators
int iother(-1);
for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix;
assert(iq1>=0&&iq2>=0&&iother>=0);
// And check if those particles are from a beam remnant
bool rem1 = cluster->isBeamRemnant(iq1);
bool rem2 = cluster->isBeamRemnant(iq2);
// workout which distribution to use
bool soft1(false),soft2(false);
switch (_iopRem) {
case 0:
soft1 = rem1 || rem2;
soft2 = rem2 || rem1;
break;
case 1:
soft1 = rem1;
soft2 = rem2;
break;
}
// Initialization for the exponential ("soft") mass distribution.
static const int max_loop = 1000;
int counter = 0;
Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO;
tcPDPtr toHadron;
bool toDiQuark(false);
PPtr newPtr1 = PPtr(),newPtr2 = PPtr();
PDPtr diquark;
bool succeeded = false;
Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
do {
succeeded = false;
++counter;
// get a flavour for the qqbar pair
drawNewFlavour(newPtr1,newPtr2,cluster);
// randomly pick which will be (anti)diquark and which a mesonic cluster
if(UseRandom::rndbool()) {
swap(iq1,iq2);
swap(rem1,rem2);
}
// check first order
if(cantMakeHadron(ptrQ[iq1], newPtr1) || !spectrum()->canMakeDiQuark(ptrQ[iq2], newPtr2)) {
swap(newPtr1,newPtr2);
}
// check again
if(cantMakeHadron(ptrQ[iq1], newPtr1) || !spectrum()->canMakeDiQuark(ptrQ[iq2], newPtr2)) {
throw Exception()
<< "ClusterFissioner cannot split the cluster ("
<< ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName()
<< ") into a hadron and diquark.\n" << Exception::runerror;
}
// Check that new clusters can produce particles and there is enough
// phase space to choose the drawn flavour
m1 = ptrQ[iq1]->data().constituentMass();
m2 = ptrQ[iq2]->data().constituentMass();
m = newPtr1->data().constituentMass();
// Do not split in the case there is no phase space available
if(mmax < m1+m + m2+m) continue;
pQ1.setMass(m1);
pQone.setMass(m);
pQtwo.setMass(m);
pQ2.setMass(m2);
double weightMasses = drawNewMasses(mmax, soft1, soft2, pClu1, pClu2,
ptrQ[iq1], pQ1, newPtr1, pQone,
newPtr2, pQtwo, ptrQ[iq1], pQ2);
if (weightMasses == 0.0) continue;
Mc1 = pClu1.mass();
Mc2 = pClu2.mass();
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue;
if ( _phaseSpaceWeights && phaseSpaceVeto(mmax,Mc1,Mc2,m,m1,m2) ) {
// reduce counter as it regards only the mass sampling
counter--;
continue;
}
// check if need to force meson clster to hadron
toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
// check if need to force diquark cluster to be on-shell
toDiQuark = false;
diquark = spectrum()->makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
if(Mc2 < diquark->constituentMass()) {
Mc2 = diquark->constituentMass(); pClu2.setMass(Mc2);
toDiQuark = true;
}
// if a beam cluster not allowed to decay to hadrons
if(cluster->isBeamCluster() && toHadron && softUEisOn)
continue;
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > mmax) {
if(!toHadron) {
toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2);
if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
}
else if(!toDiQuark) {
Mc2 = _hadronSpectrum->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); pClu2.setMass(Mc2);
toDiQuark = true;
}
}
succeeded = (mmax >= Mc1+Mc2);
}
while (!succeeded && counter < max_loop);
// check no of tries
if(counter >= max_loop) return cutType();
// Determine the (5-components) momenta (all in the LAB frame)
Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum();
calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
// positions of the new clusters
LorentzPoint pos1,pos2;
Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum();
calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2);
// first the mesonic cluster/meson
cutType rval;
if(toHadron) {
rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toDiQuark) {
rem2 |= cluster->isBeamRemnant(iother);
PPtr newDiQuark = diquark->produceParticle(pClu2);
rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2,
ptrQ[iother]->momentum(), rem2);
}
else {
rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2,
ptrQ[iother],cluster->isBeamRemnant(iother));
}
cluster->isAvailable(false);
return rval;
}
ClusterFissioner::PPair
ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
const LorentzPoint &b) const {
PPair rval;
if(hadron->coloured()) {
rval.first = (_hadronSpectrum->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle();
}
else
rval.first = hadron->produceParticle();
rval.second = newPtr;
rval.first->set5Momentum(a);
rval.first->setVertex(b);
return rval;
}
ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr,
const Lorentz5Momentum & a,
const LorentzPoint & b,
const Lorentz5Momentum & c,
const Lorentz5Momentum & d,
bool isRem,
tPPtr spect, bool remSpect) const {
PPair rval;
rval.second = newPtr;
ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect));
rval.first = cluster;
cluster->set5Momentum(a);
cluster->setVertex(b);
assert(cluster->particle(0)->id() == ptrQ->id());
cluster->particle(0)->set5Momentum(c);
cluster->particle(1)->set5Momentum(d);
cluster->setBeamRemnant(0,isRem);
if(remSpect) cluster->setBeamRemnant(2,remSpect);
return rval;
}
/**
* Calculate the phase space weight for M1*M2*(2 body PhaseSpace) ignore constituent masses
*/
double ClusterFissioner::weightFlatPhaseSpaceNoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2) const {
double M_temp = Mc/GeV;
double M1_temp = Mc1/GeV;
double M2_temp = Mc2/GeV;
if (sqr(M_temp)<sqr(M1_temp+M2_temp)) {
// This should be checked before
throw Exception()
<< "ERROR in ClusterFissioner::weightFlatPhaseSpaceNoConstituentMasses\n"
<< "ClusterFissioner has not checked Masses properly\n"
<< "Mc = " << M_temp << "\n"
<< "Mc1 = " << M1_temp << "\n"
<< "Mc2 = " << M2_temp << "\n"
<< Exception::warning;
return 0.0;
}
double lam = Kinematics::kaellen(M_temp, M1_temp, M2_temp);
double ratio;
// new weight with the Jacobi factor M1*M2 of the Mass integration
double PSweight = M1_temp*M2_temp*pow(sqrt(lam),_dim-3.);
// overestimate only possible for dim>=3.0
assert(_dim>=3.0);
// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
double overEstimate = pow(6.0*sqrt(3.0), 3.0 - _dim)*pow(M_temp, 2.*(_dim-2.));
ratio = PSweight/overEstimate;
if (!(ratio>=0) || !(ratio<=1)) {
throw Exception()
<< "ERROR in ClusterFissioner::weightFlatPhaseSpaceNoConstituentMasses\n"
<< "ratio = " <<ratio
<<" M "<<M_temp
<<" M1 "<<M1_temp
<<" M2 "<<M2_temp <<"\t"<<_dim<<"\t" << lam
<<"\t"<< overEstimate<<"\n\n"
<< Exception::runerror;
}
return ratio;
}
/**
* Calculate the phase space weight for M1*M2*(2 body PhaseSpace)^3
*/
double ClusterFissioner::weightPhaseSpaceConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2, const double power) const {
double M_temp = Mc/GeV;
double M1_temp = Mc1/GeV;
double M2_temp = Mc2/GeV;
double m_temp = m/GeV;
double m1_temp = m1/GeV;
double m2_temp = m2/GeV;
if (sqr(M_temp)<sqr(M1_temp+M2_temp)
|| sqr(M1_temp)<sqr(m1_temp+m_temp)
|| sqr(M2_temp)<sqr(m2_temp+m_temp)
) {
// This should be checked before
throw Exception()
<< "ERROR in ClusterFissioner::weightPhaseSpaceConstituentMasses\n"
<< "ClusterFissioner has not checked Masses properly\n"
<< "Mc = " << M_temp << "\n"
<< "Mc1 = " << M1_temp << "\n"
<< "Mc2 = " << M2_temp << "\n"
<< "m1 = " << m1_temp << "\n"
<< "m2 = " << m2_temp << "\n"
<< "m = " << m_temp << "\n"
<< Exception::warning;
return 0.0;
}
double lam1 = Kinematics::kaellen(M1_temp, m1_temp, m_temp);
double lam2 = Kinematics::kaellen(M2_temp, m2_temp, m_temp);
double lam3 = Kinematics::kaellen(M_temp, M1_temp, M2_temp);
double ratio;
// new weight with the Jacobi factor M1*M2 of the Mass integration
double PSweight = pow(lam1*lam2*lam3,(_dim-3.)/2.0)*pow(M1_temp*M2_temp,3.-_dim);
// overestimate only possible for dim>=3.0
assert(_dim>=3.0);
// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
double overEstimate = pow(6.0*sqrt(3.0), 3.0 - _dim)*pow(M_temp, 4.*_dim-12.);
ratio = PSweight/overEstimate;
if (!(ratio>=0)) std::cout << "ratio = " <<ratio<<" M "<<M_temp<<" M1 "<<M1_temp<<" M2 "<<M2_temp<<" m1 "<<m1_temp<<" m2 "<<m2_temp<<" m "<<m_temp<<"\t"<<_dim<<"\t" << lam1<<"\t"<< lam2<<"\t" << lam3 <<"\t"<< overEstimate<<"\n\n";
if (!(ratio>=0) || !(ratio<=1)) {
throw Exception()
<< "ERROR in ClusterFissioner::weightPhaseSpaceConstituentMasses\n"
<< "ratio = " <<ratio
<<" M "<<M_temp
<<" M1 "<<M1_temp
<<" M2 "<<M2_temp
<<" m1 "<<m1_temp
<<" m2 "<<m2_temp
<<" m "<<m_temp <<"\t"<<_dim
<<"\t" << lam1<<"\t"<< lam2<<"\t" << lam3
<<"\t"<< overEstimate<<"\n\n"
<< Exception::runerror;
}
// multiply by overestimate of power of matrix element to modulate the phase space with (M1*M2)^power
if (power) {
double powerLawOver = power<0 ? pow(Mc1*Mc2/((m1+m)*(m2+m)),power):pow(Mc1*Mc2/((Mc-(m1+m))*(Mc-(m2+m))),power);
ratio*=powerLawOver;
}
return ratio;
}
/**
* Calculate the phase space weight for M1*M2*(2 body PhaseSpace)^3
* using Hadron Masses
*/
double ClusterFissioner::weightFlatPhaseSpaceHadronMasses(const Energy Mc, const Energy Mc1, const Energy Mc2, tcPPtr pQ, tcPPtr pQ1, tcPPtr pQ2) const {
auto LHP1 = spectrum()->lightestHadronPair(pQ1->dataPtr(),pQ->dataPtr());
auto LHP2 = spectrum()->lightestHadronPair(pQ2->dataPtr(),pQ->dataPtr());
if (sqr(Mc1)<sqr(LHP1.first->mass()+LHP1.second->mass()))
return true;
if (sqr(Mc2)<sqr(LHP2.first->mass()+LHP2.second->mass()))
return true;
double lam1 = sqrt(Kinematics::kaellen(Mc1/GeV, LHP1.first->mass()/GeV, LHP1.second->mass()/GeV));
double lam2 = sqrt(Kinematics::kaellen(Mc2/GeV, LHP2.first->mass()/GeV, LHP2.second->mass()/GeV));
double lam3 = sqrt(Kinematics::kaellen(Mc/GeV, Mc1/GeV, Mc2/GeV));
double ratio;
// new weight with the Jacobi factor M1*M2 of the Mass integration
double PSweight = pow(lam1*lam2*lam3,_dim-3.)*pow(Mc1*Mc2/GeV2,3.-_dim);
// overestimate only possible for dim>=3.0
assert(_dim>=3.0);
// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
double overEstimate = pow(6.0*sqrt(3.0), 3.0 - _dim)*pow(Mc/GeV, 4.*_dim-12.);
ratio = PSweight/overEstimate;
if (!(ratio>=0) || !(ratio<=1)) {
throw Exception()
<< "ERROR in ClusterFissioner::weightFlatPhaseSpaceHadronMasses\n"
<< "ratio = " <<ratio
<<" M "<<Mc/GeV
<<" M1 "<<Mc1/GeV
<<" M2 "<<Mc2/GeV <<"\t"<<_dim<<"\t" << lam1<<"\t" << lam2 <<"\t" << lam3
<<"\t"<< overEstimate<<"\n\n"
<< Exception::runerror;
}
return ratio;
}
/**
* Veto for the phase space weight
* returns true if proposed Masses are rejected
* else returns false
*/
bool ClusterFissioner::phaseSpaceVeto(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2, tcPPtr pQ1, tcPPtr pQ2, tcPPtr pQ, const double power) const {
switch (_phaseSpaceWeights)
{
case 1:
return phaseSpaceVetoConstituentMasses(Mc, Mc1, Mc2, m, m1, m2, power);
case 2:
return phaseSpaceVetoHadronPairs(Mc, Mc1, Mc2, pQ, pQ1, pQ2);
case 3:
return phaseSpaceVetoNoConstituentMasses(Mc, Mc1, Mc2);
default:
assert(false);
}
}
/**
* Veto for the phase space weight
* returns true if proposed Masses are rejected
* else returns false
*/
bool ClusterFissioner::phaseSpaceVetoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2, const double power) const {
return (UseRandom::rnd()>weightPhaseSpaceConstituentMasses(Mc, Mc1, Mc2, m, m1, m2, power));
}
bool ClusterFissioner::phaseSpaceVetoNoConstituentMasses(const Energy Mc, const Energy Mc1, const Energy Mc2) const {
return (UseRandom::rnd()>weightFlatPhaseSpaceNoConstituentMasses(Mc, Mc1, Mc2));
}
bool ClusterFissioner::phaseSpaceVetoHadronPairs(const Energy Mc, const Energy Mc1, const Energy Mc2, tcPPtr pQ, tcPPtr pQ1, tcPPtr pQ2) const {
return (UseRandom::rnd()>weightFlatPhaseSpaceHadronMasses(Mc, Mc1, Mc2, pQ, pQ1, pQ2));
}
/**
* Calculate the masses and possibly kinematics of the cluster
* fission at hand; if calculateKineamtics is perfomring non-trivial
* steps kinematics claulcated here will be overriden. Currentl;y resorts to the default
*/
double ClusterFissioner::drawNewMasses(const Energy Mc, const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
tcPPtr, const Lorentz5Momentum& pQone,
tcPPtr, const Lorentz5Momentum& pQtwo,
tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const {
// power for splitting
double exp1 = !spectrum()->isExotic(ptrQ1->dataPtr()) ? _pSplitLight : _pSplitExotic;
double exp2 = !spectrum()->isExotic(ptrQ2->dataPtr()) ? _pSplitLight : _pSplitExotic;
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
assert(_pSplitHeavy.find(id) != _pSplitHeavy.end());
if ( spectrum()->hasHeavy(id,ptrQ1->dataPtr()) ) exp1 = _pSplitHeavy.find(id)->second;
if ( spectrum()->hasHeavy(id,ptrQ2->dataPtr()) ) exp2 = _pSplitHeavy.find(id)->second;
}
Energy M1 = drawChildMass(Mc,pQ1.mass(),pQ2.mass(),pQone.mass(),exp1,soft1);
Energy M2 = drawChildMass(Mc,pQ2.mass(),pQ1.mass(),pQtwo.mass(),exp2,soft2);
pClu1.setMass(M1);
pClu2.setMass(M2);
return 1.0; // succeeds
}
void ClusterFissioner::drawNewFlavourDiquarks(PPtr& newPtrPos,PPtr& newPtrNeg,
const ClusterPtr & clu) const {
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
unsigned hasDiquarks=0;
assert(clu->numComponents()==2);
tcPDPtr pD1=clu->particle(0)->dataPtr();
tcPDPtr pD2=clu->particle(1)->dataPtr();
bool isDiq1=DiquarkMatcher::Check(pD1->id());
if (isDiq1) hasDiquarks++;
bool isDiq2=DiquarkMatcher::Check(pD2->id());
if (isDiq2) hasDiquarks++;
assert(hasDiquarks<=2);
Energy Mc=(clu->momentum().mass());
Energy minMass;
double weight;
Selector<long> choice;
// adding quark-antiquark pairs to the selection list
for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
minMass=spectrum()->massLightestHadronPair(pD1,pD2);
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (abs(_fissionCluster)==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
// adding diquark-antidiquark pairs to the selection list
switch (hasDiquarks)
{
case 0:
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
if (_strictDiquarkKinematics) {
tPDPtr cand = getParticleData(id);
Energy mH1=spectrum()->massLightestHadron(pD2,cand);
Energy mH2=spectrum()->massLightestHadron(cand,pD1);
minMass = mH1 + mH2;
}
else {
minMass = spectrum()->massLightestBaryonPair(pD1,pD2);
}
if (Mc < minMass) continue;
if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
else assert(false);
if (_fissionCluster==-1)
weight*=sqr(Herwig::Math::alphaS(Mc, 0.25*GeV,3, 2));
choice.insert(weight,id);
}
break;
case 1:
if (_diquarkClusterFission<1) break;
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
tPDPtr diq = getParticleData(id);
if (isDiq1)
minMass = spectrum()->massLightestHadron(pD2,diq)
+ spectrum()->massLightestBaryonPair(diq,pD1);
else
minMass = spectrum()->massLightestHadron(pD1,diq)
+ spectrum()->massLightestBaryonPair(diq,pD2);
if (Mc < minMass) continue;
if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
else assert(false);
if (_fissionCluster==-1)
weight*=sqr(Herwig::Math::alphaS(Mc, 0.25*GeV,3, 2));
choice.insert(weight,id);
}
break;
case 2:
if (_diquarkClusterFission<2) break;
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
tPDPtr diq = getParticleData(id);
if (Mc < spectrum()->massLightestBaryonPair(pD1,pD2)) {
- throw Exception() << "Found Diquark Cluster:\n" << *clu << "\nwith MassCluster = "
- << Mc/GeV <<" GeV MassLightestBaryonPair = "
- << spectrum()->massLightestBaryonPair(pD1,pD2)/GeV
- << " GeV cannot decay" << Exception::eventerror;
+ throw Exception() << "Found Diquark Cluster:\n" << *clu << "\nwith MassCluster = "
+ << Mc/GeV <<" GeV MassLightestBaryonPair = "
+ << spectrum()->massLightestBaryonPair(pD1,pD2)/GeV
+ << " GeV cannot decay" << Exception::eventerror;
}
minMass = spectrum()->massLightestBaryonPair(pD1,diq)
+ spectrum()->massLightestBaryonPair(diq,pD2);
if (Mc < minMass) continue;
if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
else assert(false);
if (_fissionCluster==-1)
weight*=sqr(Herwig::Math::alphaS(Mc, 0.25*GeV,3, 2));
choice.insert(weight,id);
}
break;
default:
assert(false);
}
assert(choice.size()>0);
long idNew = choice.select(UseRandom::rnd());
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert(newPtrPos);
assert(newPtrNeg);
assert(newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
void ClusterFissioner::drawNewFlavourQuarks(PPtr& newPtrPos,PPtr& newPtrNeg) const {
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
Selector<long> choice;
switch(abs(_fissionCluster)){
case 0:
for ( const long& id : spectrum()->lightHadronizingQuarks() )
choice.insert(_hadronSpectrum->pwtQuark(id),id);
break;
case 1:
for ( const long& id : spectrum()->lightHadronizingQuarks() )
choice.insert(_fissionPwt.find(id)->second,id);
break;
default :
assert(false);
}
long idNew = choice.select(UseRandom::rnd());
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg,
Energy2 mass2) const {
if ( spectrum()->gluonId() != ParticleID::g )
throw Exception() << "strange enhancement only working with Standard Model hadronization"
<< Exception::runerror;
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
double prob_d = 0.;
double prob_u = 0.;
double prob_s = 0.;
double scale = abs(double(sqr(_m0Fission)/mass2));
// Choose which splitting weights you wish to use
switch(abs(_fissionCluster)){
// 0: ClusterFissioner and ClusterDecayer use the same weights
case 0:
prob_d = _hadronSpectrum->pwtQuark(ParticleID::d);
prob_u = _hadronSpectrum->pwtQuark(ParticleID::u);
/* Strangeness enhancement:
Case 1: probability scaling
Case 2: Exponential scaling
*/
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_hadronSpectrum->pwtQuark(ParticleID::s),scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
/* 1: ClusterFissioner uses its own unique set of weights,
i.e. decoupled from ClusterDecayer */
case 1:
prob_d = _fissionPwt.find(ParticleID::d)->second;
prob_u = _fissionPwt.find(ParticleID::u)->second;
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwt.find(ParticleID::s)->second,scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
default:
assert(false);
}
int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
long idNew = 0;
switch (choice) {
case 0: idNew = ThePEG::ParticleID::u; break;
case 1: idNew = ThePEG::ParticleID::d; break;
case 2: idNew = ThePEG::ParticleID::s; break;
}
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster) const {
Lorentz5Momentum pIn = cluster->momentum();
Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() +
cluster->particle(1)->mass());
Energy2 singletm2 = pIn.m2();
// Return either the cluster mass, or the lambda measure
return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2;
}
Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1,
const Energy m2, const Energy m,
const double expt, const bool soft) const {
/***************************
* This method, given in input the cluster mass Mclu of an heavy cluster C,
* made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2
* of, respectively, the children cluster C1, made of constituent masses m1
* and m, and cluster C2, of mass Mclu2 and made of constituent masses m2
* and m. The mass is extracted from one of the two following mass
* distributions:
* --- power-like ("normal" distribution)
* d(Prob) / d(M^exponent) = const
* where the exponent can be different from the two children C1 (exp1)
* and C2 (exponent2).
* --- exponential ("soft" distribution)
* d(Prob) / d(M^2) = exp(-b*M)
* where b = 2.0 / average.
* Such distributions are limited below by the masses of
* the constituents quarks, and above from the mass of decaying cluster C.
* The choice of which of the two mass distributions to use for each of the
* two cluster children is dictated by iRemnant (see below).
* If the number of attempts to extract a pair of mass values that are
* kinematically acceptable is above some fixed number (max_loop, see below)
* the method gives up and returns false; otherwise, when it succeeds, it
* returns true.
*
* These distributions have been modified from HERWIG:
* Before these were:
* Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 );
* The new one coded here is a more efficient version, same density
* but taking into account 'in phase space from' beforehand
***************************/
// hard cluster
if(!soft) {
return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt),
pow(m*UnitRemoval::InvE, expt)), 1./expt
)*UnitRemoval::E + m1;
}
// Otherwise it uses a soft mass distribution
else {
static const InvEnergy b = 2.0 / _btClM;
Energy max = M-m1-m2-2.0*m;
double rmin = b*max;
rmin = ( rmin < 50 ) ? exp(-rmin) : 0.;
double r1;
do {
r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0);
}
while (r1 < rmin);
return m1 + m - log(r1)/b;
}
}
void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu,
const Lorentz5Momentum & p0Q1,
const bool toHadron1,
const bool toHadron2,
Lorentz5Momentum & pClu1,
Lorentz5Momentum & pClu2,
Lorentz5Momentum & pQ1,
Lorentz5Momentum & pQbar,
Lorentz5Momentum & pQ,
Lorentz5Momentum & pQ2bar) const {
/******************
* This method solves the kinematics of the two body cluster decay:
* C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar)
* In input we receive the momentum of C, pClu, and the momentum
* of the quark Q1 (constituent of C), p0Q1, both in the LAB frame.
* Furthermore, two boolean variables inform whether the two fission
* products (C1, C2) decay immediately into a single hadron (in which
* case the cluster itself is identify with that hadron) and we do
* not have to solve the kinematics of the components (Q1,Qbar) for
* C1 and (Q,Q2bar) for C2.
* The output is given by the following momenta (all 5-components,
* and all in the LAB frame):
* pClu1 , pClu2 respectively of C1 , C2
* pQ1 , pQbar respectively of Q1 , Qbar in C1
* pQ , pQ2bar respectively of Q , Q2 in C2
* The assumption, suggested from the string model, is that, in C frame,
* C1 and its constituents Q1 and Qbar are collinear, and collinear to
* the direction of Q1 in C (that is before cluster decay); similarly,
* (always in the C frame) C2 and its constituents Q and Q2bar are
* collinear (and therefore anti-collinear with C1,Q1,Qbar).
* The solution is then obtained by using Lorentz boosts, as follows.
* The kinematics of C1 and C2 is solved in their parent C frame,
* and then boosted back in the LAB. The kinematics of Q1 and Qbar
* is solved in their parent C1 frame and then boosted back in the LAB;
* similarly, the kinematics of Q and Q2bar is solved in their parent
* C2 frame and then boosted back in the LAB. In each of the three
* "two-body decay"-like cases, we use the fact that the direction
* of the motion of the decay products is known in the rest frame of
* their parent. This is obvious for the first case in which the
* parent rest frame is C; but it is also true in the other two cases
* where the rest frames are C1 and C2. This is because C1 and C2
* are boosted w.r.t. C in the same direction where their components,
* respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame
* respectively.
* Of course, although the notation used assumed that C = (Q1 Q2bar)
* where Q1 is a quark and Q2bar an antiquark, indeed everything remain
* unchanged also in all following cases:
* Q1 quark, Q2bar antiquark; --> Q quark;
* Q1 antiquark , Q2bar quark; --> Q antiquark;
* Q1 quark, Q2bar diquark; --> Q quark
* Q1 antiquark, Q2bar anti-diquark; --> Q antiquark
* Q1 diquark, Q2bar quark --> Q antiquark
* Q1 anti-diquark, Q2bar antiquark; --> Q quark
**************************/
// Calculate the unit three-vector, in the C frame, along which
// all of the constituents and children clusters move.
Lorentz5Momentum u(p0Q1);
u.boost( -pClu.boostVector() ); // boost from LAB to C
// the unit three-vector is then u.vect().unit()
// Calculate the momenta of C1 and C2 in the (parent) C frame first,
// where the direction of C1 is u.vect().unit(), and then boost back in the
// LAB frame.
if (pClu.m() < pClu1.mass() + pClu2.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),
u.vect().unit(), pClu1, pClu2);
// In the case that cluster1 does not decay immediately into a single hadron,
// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron1) {
if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(),
u.vect().unit(), pQ1, pQbar);
}
// In the case that cluster2 does not decay immediately into a single hadron,
// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron2) {
if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(),
u.vect().unit(), pQ, pQ2bar);
}
}
void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu,
const LorentzPoint & positionClu,
const Lorentz5Momentum & pClu1,
const Lorentz5Momentum & pClu2,
LorentzPoint & positionClu1,
LorentzPoint & positionClu2) const {
// Determine positions of cluster children.
// See Marc Smith's thesis, page 127, formulas (4.122) and (4.123).
Energy Mclu = pClu.m();
Energy Mclu1 = pClu1.m();
Energy Mclu2 = pClu2.m();
// Calculate the unit three-vector, in the C frame, along which
// children clusters move.
Lorentz5Momentum u(pClu1);
u.boost( -pClu.boostVector() ); // boost from LAB to C frame
// the unit three-vector is then u.vect().unit()
Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2);
// First, determine the relative positions of the children clusters
// in the parent cluster reference frame.
Energy2 mag2 = u.vect().mag2();
InvEnergy fact = mag2>ZERO ? 1./sqrt(mag2) : 1./GeV;
Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t1 = Mclu/_kappa - x1;
LorentzDistance distanceClu1( x1 * fact * u.vect(), t1 );
Length x2 = (-0.25*Mclu + 0.5*(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t2 = Mclu/_kappa + x2;
LorentzDistance distanceClu2( x2 * fact * u.vect(), t2 );
// Then, transform such relative positions from the parent cluster
// reference frame to the Lab frame.
distanceClu1.boost( pClu.boostVector() );
distanceClu2.boost( pClu.boostVector() );
// Finally, determine the absolute positions in the Lab frame.
positionClu1 = positionClu + distanceClu1;
positionClu2 = positionClu + distanceClu2;
}
bool ClusterFissioner::ProbabilityFunction(double scale, double threshold) {
double cut = UseRandom::rnd(0.0,1.0);
return 1./(1.+pow(abs((threshold-_probShift)/scale),_probPowFactor)) > cut ? true : false;
}
bool ClusterFissioner::ProbabilityFunctionPower(double Mass, double threshold) {
double cut = UseRandom::rnd(0.0,1.0);
if ((Mass-threshold)<=0)
return false;
return 1.0/(1.0 + _probPowFactor*pow(1.0/(Mass-threshold),_clPowLight)) > cut ? true : false;
}
bool ClusterFissioner::isHeavy(tcClusterPtr clu) {
// particle data for constituents
tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
bool hasDiquark=0;
for(size_t ix=0;ix<min(clu->numComponents(),3);++ix) {
cptr[ix]=clu->particle(ix)->dataPtr();
// Assuming diquark masses are ordered with larger id corresponding to larger masses
if (DiquarkMatcher::Check(*(cptr[ix]))) {
hasDiquark=true;
break;
}
}
// different parameters for exotic, bottom and charm clusters
double clpow = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clPowLight : _clPowExotic;
Energy clmax = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clMaxLight : _clMaxExotic;
// if no heavy quark is found in the cluster, but diquarks are present use
// different ClMax and ClPow
if ( hasDiquark) {
clpow = _clPowDiquark;
clmax = _clMaxDiquark;
}
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( spectrum()->hasHeavy(id,cptr[0],cptr[1],cptr[1])) {
clpow = _clPowHeavy[id];
clmax = _clMaxHeavy[id];
}
}
// required test for SUSY clusters, since aboveCutoff alone
// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
static const Energy minmass
= getParticleData(ParticleID::d)->constituentMass();
bool aboveCutoff = false, canSplitMinimally = false;
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
aboveCutoff = (
pow(clu->mass()*UnitRemoval::InvE , clpow)
>
pow(clmax*UnitRemoval::InvE, clpow)
+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
);
canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
}
// dynamic kinematic threshold
else if(_kinematicThresholdChoice == 1) {
//some smooth probablity function to create a dynamic thershold
double scale = pow(clu->mass()/GeV , clpow);
double threshold = pow(clmax/GeV, clpow)
+ pow(clu->sumConstituentMasses()/GeV, clpow);
aboveCutoff = ProbabilityFunction(scale,threshold);
scale = clu->mass()/GeV;
threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
canSplitMinimally = ProbabilityFunction(scale,threshold);
}
// probablistic kinematic threshold
else if(_kinematicThresholdChoice == 2) {
// Consistent power law for CF probability
double Mass = clu->mass()/GeV;
double threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
aboveCutoff = ProbabilityFunctionPower(Mass,threshold + clmax/GeV);
canSplitMinimally = Mass - threshold>ZERO;
}
return aboveCutoff && canSplitMinimally;
}
diff --git a/MatrixElement/FxFx/FxFxHandler.cc b/MatrixElement/FxFx/FxFxHandler.cc
--- a/MatrixElement/FxFx/FxFxHandler.cc
+++ b/MatrixElement/FxFx/FxFxHandler.cc
@@ -1,1655 +1,1665 @@
#include "FxFxHandler.h"
#include "FxFxReader.h"
#include "FxFxReader.fh"
#include "FxFxEventHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include <queue>
#include "ThePEG/Utilities/Throw.h"
#include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"
#include <boost/algorithm/string.hpp>
using namespace Herwig;
using namespace ThePEG;
bool recordEntry(PPtr i,PPtr j) {
return (i->number()<j->number());
}
bool pTsortFunction(PPtr i,PPtr j) {
return (i->momentum().perp2()>j->momentum().perp2());
}
bool ETsortFunction(pair<Energy, Lorentz5Momentum> i,pair<Energy, Lorentz5Momentum> j) {
return (i.first>j.first);
}
bool isMomLessThanEpsilon(Lorentz5Momentum p,Energy epsilon) {
return (abs(p.x())<epsilon&&abs(p.y())<epsilon&&
abs(p.z())<epsilon&&abs(p.t())<epsilon);
}
FxFxHandler::FxFxHandler()
: ncy_(100),ncphi_(60),ihvy_(-999),nph_(-999),nh_(-999),
mergemode_(0), vetoHeavyQ_(true), vetoHeavyFlavour_(false), etclusmean_(20*GeV), jetAlgorithm_(1),vetoIsTurnedOff_(false),vetoSoftThanMatched_(false),etclusfixed_(true),epsetclus_(2.5*GeV),rclus_(0.4),etaclmax_(5.0),rclusfactor_(1.5), hpdetect_(true), ihrd_(-999),njets_(-999),drjmin_(-999), highestMultiplicity_(false), ycmax_(5.4),ycmin_(-5.4)
{}
void FxFxHandler::doinitrun() {
QTildeShowerHandler::doinitrun();
// et_ holds the ET deposited in the (ncy_ x ncphi_) calorimeter cells.
et_.resize(ncy_);
for(unsigned int ixx=0; ixx<et_.size(); ixx++) et_[ixx].resize(ncphi_);
// jetIdx_ for a given calorimeter cell this holds the index of the jet
// that the cell was clustered into.
jetIdx_.resize(ncy_);
for(unsigned int ixx=0; ixx<jetIdx_.size(); ixx++) jetIdx_[ixx].resize(ncphi_);
if(mergemode_ == 0) {
cout << "Merging mode is FxFx." << endl;
} else if (mergemode_ == 1) {
cout << "Merging mode is Tree level." << endl;
} else if (mergemode_ == 2) {
cout << "Merging mode is Tree level, using MadGraph Les Houches information." << endl;
}
if(hpdetect_ && mergemode_ != 2) {
cout << "Automatic detection of hard process enabled." << endl;
ihrd_ = -999;
}
}
IBPtr FxFxHandler::clone() const {
return new_ptr(*this);
}
IBPtr FxFxHandler::fullclone() const {
return new_ptr(*this);
}
void FxFxHandler::persistentOutput(PersistentOStream & os) const {
os << alphaS_
<< ncy_ << ncphi_ << ihvy_ << nph_ << nh_
<< ounit(etclusmean_,GeV) << rclus_ << etaclmax_ << rclusfactor_
<< ihrd_ << njets_ << drjmin_ << highestMultiplicity_
<< ycmax_ << ycmin_ << jetAlgorithm_ << vetoIsTurnedOff_ << vetoSoftThanMatched_ << etclusfixed_
<< cphcal_ << sphcal_ << cthcal_ << sthcal_ << ounit(epsetclus_,GeV) << vetoHeavyQ_ << mergemode_ << hpdetect_ << vetoHeavyFlavour_;
}
void FxFxHandler::persistentInput(PersistentIStream & is, int) {
is >> alphaS_
>> ncy_ >> ncphi_ >> ihvy_ >> nph_ >> nh_
>> iunit(etclusmean_,GeV) >> rclus_ >> etaclmax_ >> rclusfactor_
>> ihrd_ >> njets_ >> drjmin_ >> highestMultiplicity_
>> ycmax_ >> ycmin_ >> jetAlgorithm_ >> vetoIsTurnedOff_ >> vetoSoftThanMatched_ >> etclusfixed_
>> cphcal_ >> sphcal_ >> cthcal_ >> sthcal_ >> iunit(epsetclus_,GeV) >> vetoHeavyQ_ >> mergemode_ >> hpdetect_ >> vetoHeavyFlavour_;
}
ClassDescription<FxFxHandler> FxFxHandler::initFxFxHandler;
// Definition of the static class description member.
void FxFxHandler::Init() {
static ClassDocumentation<FxFxHandler> documentation
("The FxFxHandler class performs MEPS merging "
"using the MLM procedure.");
static Reference<FxFxHandler,ShowerAlpha> interfaceShowerAlpha
("ShowerAlpha",
"The object calculating the strong coupling constant",
&FxFxHandler::alphaS_, false, false, true, false, false);
static Parameter<FxFxHandler,int> interfaceihvy
("ihvy",
"heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t)",
&FxFxHandler::ihvy_, -999, -999, 7,
false, false, Interface::limited);
static Parameter<FxFxHandler,int> interfacenph
("nph",
"Number of photons in the AlpGen process",
&FxFxHandler::nph_, -999, -999, 7,
false, false, Interface::limited);
static Parameter<FxFxHandler,int> interfacenh
("nh",
"Number of higgses in the AlpGen process",
&FxFxHandler::nph_, -999, -999, 7,
false, false, Interface::limited);
static Parameter<FxFxHandler,Energy> interfaceETClus
("ETClus",
"The ET threshold defining a jet in the merging procedure",
&FxFxHandler::etclusmean_, GeV, 20*GeV, 0*GeV, 14000*GeV,
false, false, Interface::limited);
static Parameter<FxFxHandler,double> interfaceRClus
("RClus",
"The cone size used to define a jet in the merging procedure",
&FxFxHandler::rclus_, 0.4, 0.0, 4.0,
false, false, Interface::limited);
static Parameter<FxFxHandler,double> interfaceEtaClusMax
("EtaClusMax",
"The maximum |eta| used to define a jet in the merging procedure",
&FxFxHandler::etaclmax_, 5.0, 0.0, 15.0,
false, false, Interface::limited);
static Parameter<FxFxHandler,double> interfaceRClusFactor
("RClusFactor",
"The prefactor for RClus used to define the jet-parton matching "
"distance",
&FxFxHandler::rclusfactor_, 1.5, 0.0, 4.0,
false, false, Interface::limited);
static Parameter<FxFxHandler,int> interfaceihrd
("ihrd",
"The hard process code",
&FxFxHandler::ihrd_, -999, 0, 10000,
false, false, Interface::limited);
static Parameter<FxFxHandler,int> interfacenjetsmax
("njetsmax",
"The number of light jets in the maximum-multiplicity process",
&FxFxHandler::njets_, -999, 0, 10000,
false, false, Interface::limited);
static Parameter<FxFxHandler,double> interfacedrjmin
("drjmin",
"Mimimum parton-parton R-sep used for generation.",
&FxFxHandler::drjmin_, 0.7, 0.0, 4.0,
false, false, Interface::limited);
- static Parameter<FxFxHandler,bool> interfacehighestMultiplicity
+ static Switch<FxFxHandler,bool> interfacehighestMultiplicity
("highestMultiplicity",
- "If true it indicates that this is the highest multiplicity input "
- "ME-level configuration to be processed.",
- &FxFxHandler::highestMultiplicity_, 0, 0, 1,
- false, false, Interface::limited);
+ "If true it indicates that this is the highest multiplicity input ME-level configuration to be processed.",
+ &FxFxHandler::highestMultiplicity_, false, false, false);
+ static SwitchOption interfacehighestMultiplicityYes
+ (interfacehighestMultiplicity,
+ "Yes",
+ "Highest Multiplicity",
+ true);
+ static SwitchOption interfacehighestMultiplicityNo
+ (interfacehighestMultiplicity,
+ "No",
+ "Not the highest multiplicity",
+ false);
- static Parameter<FxFxHandler,bool> interfaceETClusFixed
+ static Switch<FxFxHandler,bool> interfaceETClusFixed
("ETClusFixed",
- "If false, indicates that the jet merging scale, etclus_ is allowed to vary"
- "according to epsetclus_",
- &FxFxHandler::etclusfixed_, 1, 0, 1,
- false, false, Interface::limited);
-
- static Parameter<FxFxHandler,Energy> interfaceEpsilonETClus
- ("EpsilonETClus",
- "The ET threshold defining a jet in the merging procedure",
- &FxFxHandler::epsetclus_, GeV, 2.5*GeV, 0*GeV, 100.0*GeV,
- false, false, Interface::limited);
+ "If false, indicates that the jet merging scale, etclus_ is allowed to vary according to epsetclus_",
+ &FxFxHandler::etclusfixed_, true, false, false);
+ static SwitchOption interfaceETClusFixedYes
+ (interfaceETClusFixed,
+ "Yes",
+ "Fixed",
+ true);
+ static SwitchOption interfaceETClusFixedNo
+ (interfaceETClusFixed,
+ "No",
+ "Varying",
+ false);
static Switch<FxFxHandler,int> interfaceJetAlgorithm
("JetAlgorithm",
"Determines the jet algorithm for finding jets in parton-jet "
"matching in the MLM procedure.",
&FxFxHandler::jetAlgorithm_, 1, false, false);
static SwitchOption AntiKt
(interfaceJetAlgorithm,
"AntiKt",
"The anti-kt jet algorithm.",
-1);
static SwitchOption CambridgeAachen
(interfaceJetAlgorithm,
"CambridgeAachen",
"The Cambridge-Aachen jet algorithm.",
0);
static SwitchOption Kt
(interfaceJetAlgorithm,
"Kt",
"The Kt jet algorithm.",
1);
static Switch<FxFxHandler,int> interfaceMergeMode
("MergeMode",
"The choice of merging mode",
&FxFxHandler::mergemode_, 0, false, false);
static SwitchOption FxFx
(interfaceMergeMode,
"FxFx",
"FxFx merging.",
0);
static SwitchOption Tree
(interfaceMergeMode,
"Tree",
"Tree-level merging.",
1);
static SwitchOption TreeMG
(interfaceMergeMode,
"TreeMG5",
"Tree-level merging using the MadGraph pt clustering information.",
2);
static Switch<FxFxHandler,bool> interfaceHardProcessDetection
("HardProcessDetection",
"The choice of merging mode",
&FxFxHandler::hpdetect_, true, false, false);
static SwitchOption Automatic
(interfaceHardProcessDetection,
"Automatic",
"Automatically determine which particles to include in the merging.",
true);
static SwitchOption Manual
(interfaceHardProcessDetection,
"Manual",
"Use the ihrd code to determine which particles to include in the merging.",
false);
static Switch<FxFxHandler,bool> interfaceVetoIsTurnedOff
("VetoIsTurnedOff",
"Allows the vetoing mechanism to be switched off.",
&FxFxHandler::vetoIsTurnedOff_, false, false, false);
static SwitchOption VetoingIsOn
(interfaceVetoIsTurnedOff,
"VetoingIsOn",
"The MLM merging veto mechanism is switched ON.",
false);
static SwitchOption VetoingIsOff
(interfaceVetoIsTurnedOff,
"VetoingIsOff",
"The MLM merging veto mechanism is switched OFF.",
true);
static Switch<FxFxHandler,bool> interfaceVetoHeavyFlavour
("VetoHeavyFlavour",
"Allows the heavy flavour vetoing mechanism to be switched off.",
&FxFxHandler::vetoHeavyFlavour_, false, false, false);
static SwitchOption HeavyFVetoingIsOn
(interfaceVetoHeavyFlavour,
"Yes",
"The MLM merging veto mechanism for heavy flavour is switched ON.",
true);
static SwitchOption HeavyFVetoingIsOff
(interfaceVetoHeavyFlavour,
"No",
"The MLM merging veto mechanism for heavy flavour is switched OFF.",
false);
static Switch<FxFxHandler,bool> interfaceHeavyQVeto
("HeavyQVeto",
"Allows the vetoing mechanism on the heavy quark products to be switched off.",
&FxFxHandler::vetoHeavyQ_, false, false, false);
static SwitchOption HQVetoingIsOn
(interfaceHeavyQVeto,
"Yes",
"The MLM merging veto on Heavy quark decay produts mechanism is switched ON.",
true);
static SwitchOption HQVetoingIsOff
(interfaceHeavyQVeto,
"No",
"The MLM merging veto on Heavy quark decay products mechanism is switched OFF.",
false);
static Switch<FxFxHandler,bool> interfaceVetoSoftThanMatched
("VetoSoftThanMatched",
"Allows the vetoing mechanism to be switched off.",
&FxFxHandler::vetoSoftThanMatched_, false, false, false);
static SwitchOption VetoSoftIsOn
(interfaceVetoSoftThanMatched,
"VetoSoftIsOn",
"The vetoing of highest-mult. events with jets softer than matched ones is ON",
true);
static SwitchOption VetoSoftIsOff
(interfaceVetoSoftThanMatched,
"VetoSoftIsOff",
"The vetoing of highest-mult. events with jets softer than matched ones is OFF.",
false);
}
void FxFxHandler::dofinish() {
QTildeShowerHandler::dofinish();
}
void FxFxHandler::doinit() {
//print error if HardProcID is not set in input file
if(ihrd_ == -999 && !hpdetect_) { cout << "Error: FxFxHandler:ihrd not set and FxFx:HardProcessDetection set to Manual!" << endl; exit(1); }
QTildeShowerHandler::doinit();
// Compute calorimeter edges in rapidity for GetJet algorithm.
ycmax_=etaclmax_+rclus_;
ycmin_=-ycmax_;
}
// Throws a veto according to MLM strategy ... when we finish writing it.
bool FxFxHandler::showerHardProcessVeto() const {
int debug_mode = 0;
if(vetoIsTurnedOff_) {
// cout << "Vetoing is turned OFF." << endl;
return false;
}
//if(debug_mode) { cout << "debug_mode = " << 5 << endl; }
// Skip veto for processes in which merging is not implemented:
if(ihrd_==7||ihrd_==8||ihrd_==13) {
ostringstream wstring;
wstring << "FxFxHandler::showerHardProcessVeto() - warning."
<< "MLM merging not implemented "
<< "processes 4Q (ihrd=7), QQh (ihrd=8), "
<< "(single) top (ihrd=13) \n";
generator()->logWarning( Exception(wstring.str(),
Exception::warning) );
return false;
}
// Fill preshowerISPs_ pair and preshowerFSPs_ particle pointer vector.
getPreshowerParticles();
// Fill showeredISHs_, showeredISPs and showeredRems pairs, as well as
// showeredFSPs_ particle pointer vector.
getShoweredParticles();
// Turn on some screen output debugging: 0 = none ---> 5 = very verbose.
doSanityChecks(debug_mode);
// Dimensions of each calorimter cell in y and phi.
dely_ = (ycmax_-ycmin_)/double(ncy_);
delphi_ = 2*M_PI/double(ncphi_);
// Fill partonsToMatch_ with only those pre-shower partons intended to
// used in jet-parton matching and fill particlesToCluster_ using only
// those final state particles (post-shower) which are supposed to go
// in the jet clustering used to do merging.
partonsToMatch_ = preshowerFSPs_;
particlesToCluster_ = showeredFSPs_ ;
// Filter out all but the 'extra' light-parton progenitors and their
// associated final state particles.
if(mergemode_ == 0 || mergemode_ == 1) { caldel_m(); }
else if(mergemode_ == 2) { caldel_mg(); }
double prob(1);
//if etclusfixed_ then set the etclus_ to the fixed chosen value
if(etclusfixed_) {
etclus_ = etclusmean_;
} else {
//else, if we wish to vary etclus_, we use the probability distribution
//choose a probability between 0 and 1
prob = rnd();
etclus_ = etclusran_(prob);
}
// Cluster particlesToCluster_ into jets with FastJet.
getFastJets(rclus_,etclus_,etaclmax_);
if(mergemode_ == 0) {
// Get npLO_ and npNLO_ for FxFx matching
getnpFxFx();
// print the npXLO_ values obtained
// cout << "HANDLER:\t\t\t\t" << npLO_ << "\t\t" << npNLO_ << endl;
//FxFx modifications start here.
// Sort partonsToMatch_ from high to low pT.
sort(partonsToMatch_.begin(),partonsToMatch_.end(),pTsortFunction);
// Count the number of jets.
int njets_found(pjet_.size());
// If the number of jets found is not equal to the number of partons in the Born
// (i.e., the number of partons in the S-event, or one less than the number of partons in an H-event),
// the jets cannot be matched and the event has to be rejected. The number of partons in the Born is written in the event file with a name “npNLO”
// if there are no jets to match and no jets have been found, do not veto.
if(njets_found == 0 && npNLO_ == 0) { /*cout << "njets_found = " << njets_found << " and npNLO = " << npNLO_ << ", accepting" << endl;*/ return false; }
//if the number of jets is smaller than npNLO -> reject the event.
if(njets_found < npNLO_) { /*cout << "njets_found = " << njets_found << " and npNLO = " << npNLO_ << ", rejecting" << endl;*/ return true; }
// For the maximum-multiplicity sample, the number of jets obtained does not have to be exactly equal to npNLO, it may also be larger;
if(njets_found > npNLO_ && npNLO_ != njets_) { /*cout << "njets_found = " << njets_found << " and npNLO = " << npNLO_ << ", rejecting" << endl;*/ return true; }
// Create the matrix-element jets.
// Cluster also the partons at the hard-matrix element level into jets with the same algorithm as above,
// but without the requirement of a minimal pT on the jets (or set it very small).
// By construction, for S-events you should find exactly npNLO jets, while for the H-events it is either npNLO or npNLO+1.
// Cluster partonsToMatch_ into jets with FastJet.
getFastJetsToMatch(rclus_);
//int me_njets_found(pjetME_.size());
// cout << "number of ME jets found = " << me_njets_found << "partons to match: " << partonsToMatch_.size() << endl;
// Match light progenitors to jets.
vector<int> jetToPartonMap(pjetME_.size(),-999);
Energy etmin(777e100*GeV);
// Match the jets.
// Try to match the “npNLO” hardest jets created post-shower with any of the jets pre-shower. Two jets are matched if the distance between them is smaller than 1.5*DeltaR.
// If not all the npNLO hardest shower jets are matched the event has to be rejected.
// Note that if the current event does not belong to the maximum multiplicity sample, this means that all the shower jets need to be matched, because the requirement above already rejects
// events that do not have npNLO shower jets.
// For those events, at the level of the matrix elements there can either be npNLO or npNLO+1 matrix-element jets, depending on S- or H-events and the kinematics of those partons.
// Still only the shower jets need to be matched, so an event should not be rejected if a matrix-element jet cannot be matched.
// For each parton, starting with the hardest one ...
for(int ixx=0; ixx<npNLO_; ixx++) {
// ... loop over all jets not already matched.
double DRmin(777e100);
int jetIndexForDRmin(-999);
for(unsigned int jxx=0; jxx<pjetME_.size(); jxx++) {
// ... and tag closest of the remaining ones
double DRpartonJet(partonJetDeltaR(pjetME_[ixx],pjet_[jxx]));
if(jetToPartonMap[jxx]<0&&DRpartonJet<DRmin) {
DRmin=DRpartonJet;
jetIndexForDRmin=jxx;
}
}
// If the parton-jet distance is less than the matching
// distance, the parton and jet match.
if(DRmin<rclus_*rclusfactor_&&jetIndexForDRmin>=0) {
jetToPartonMap[jetIndexForDRmin]=ixx;
if(ixx==0||etjet_[jetIndexForDRmin]<etmin)
etmin=etjet_[jetIndexForDRmin];
// Otherwise this parton is not matched so veto the event.
} else return true;
}
// Veto events where matched jets are softer than non-matched ones,
// in the inclusive (highestMultiplicity_ = true) mode, unless we
// are dealing with NLO input events.
if(npNLO_ == njets_ && vetoSoftThanMatched_) {
//cout << "highest mult. event being tested for softer jets than matched..." << endl;
for(unsigned int iyy=0; iyy<pjet_.size(); iyy++) {
//cout << "etjet_[iyy] = " << etjet_[iyy]/GeV << endl;
if(jetToPartonMap[iyy]<0&&etmin<etjet_[iyy]) { /*cout << "VETO!" << endl;*/ return true; }
}
}
if(!vetoHeavyQ_) {
//cout << "no heavy quark decay product veto!" << endl;
return false;
}
//end of FxFx part
// **************************************************************** //
// * Now look to the non-light partons for heavy quark processes. * //
// **************************************************************** //
if( (ihrd_<=2||ihrd_==6||ihrd_==10||ihrd_==15||ihrd_==16) || (hpdetect_==true && hvqfound==true) ) {
// Extract heavy quark progenitors and the radiation they
// produce and put it in the calorimeter.
caldel_hvq();
// Cluster particlesToCluster_ into jets with FastJet.
getFastJets(rclus_,etclus_,etaclmax_);
// If the radiation from the heavy quarks does not give rise
// to any jets we accept event.
if(pjet_.size() == 0) return false;
// If extra jets emerge from the jet clustering we only
// accept events where the jets formed by radiation from
// b and c quarks lies within drjmin_ of the heavy quark
// progenitor.
int nmjet(pjet_.size());
for(unsigned int ixx=0; ixx<pjet_.size(); ixx++) {
for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
if(!(abs(partonsToMatch_[jxx]->id())==4||abs(partonsToMatch_[jxx]->id())==5)) continue;
if(partonJetDeltaR(partonsToMatch_[jxx],pjet_[ixx])<drjmin_) {
nmjet--; // Decrease the number of unmatched jets.
etjet_[ixx]=0*GeV; // Set jet ET to zero to indicate it is 'matched'.
}
}
}
// If every jet matched to _at_least_one_ progenitor accept the event.
if(nmjet<=0) return false;
}
}
/*****
*****
***** END of mergemode_ == 0 (i.e. FxFx MERGING BLOCK)
*****/
if(mergemode_ == 1 || mergemode_ == 2) {
// determine whether event is of highest multiplicity or not
unsigned long int njets_long = njets_;
if(partonsToMatch_.size()==njets_long) { highestMultiplicity_ = true; }
// cout << "jet finding gives pjet_.size() = " << pjet_.size() << endl;
// If there are less jets than partons then parton-jet matching is
// bound to fail: reject the event already. Also, if the input is
// an NLO event file it will 99.5% of the time contain a number of
// light partons in the F.S. equal to that in the real emission
// process in the NLO calculation, moreover, it has already
// effectively merged njets_-1 and njets jet events. So in that
// case we do not reject events on the grounds they have jet
// multiplicity less than partonsToMatch_.size() but rather less
// jets than partonsToMatch.size()-1; such events are better
// described by the lower-by-one-unit (partonsToMatch_.size()-1)
// of multiplicity NLO event file, or the lower-by-two-units
// (partonsToMatch_.size()-2) of multiplicity LO event file.
// If it is not jet production apply rejection criterion as above.
if(ihrd_!=9) {
if(pjet_.size() < partonsToMatch_.size()) return true;
}
// Otherwise, in the case of jet production allow the lowest
// contributing multiplicity process (just at NLO), namely,
// dijet production, to give rise to 1-jet and even 0-jet
// events, since these can contribute to, for example, the
// inclusive jet cross section i.e. in this case the rejection
// is only applied in the case of the next-to-lowest multiplicity
// processes (>2 parton events at LO and >3 parton events at NLO).
else {
// KH - March 5th
// Removed the following line giving special treatment
// also to the LO events, to maintain consistency with
// the fortran algorithm, at least for now. So now jet
// production at LO is being treated the same as all
// other processes.
// if(partonsToMatch_.size()==2 && pjet_.size()<2) return false;
if(pjet_.size() < partonsToMatch_.size()) return true;
}
// Sort partonsToMatch_ from high to low pT.
sort(partonsToMatch_.begin(),partonsToMatch_.end(),pTsortFunction);
// Match light progenitors to jets.
vector<int> jetToPartonMap(pjet_.size(),-999);
Energy etmin(777e100*GeV);
// For each parton, starting with the hardest one ...
for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
// ... loop over all jets not already matched.
double DRmin(777e100);
int jetIndexForDRmin(-999);
for(unsigned int jxx=0; jxx<pjet_.size(); jxx++) {
// ... and tag closest of the remaining ones
double DRpartonJet(partonJetDeltaR(partonsToMatch_[ixx],pjet_[jxx]));
if(jetToPartonMap[jxx]<0&&DRpartonJet<DRmin) {
DRmin=DRpartonJet;
jetIndexForDRmin=jxx;
}
}
// If the parton-jet distance is less than the matching
// distance, the parton and jet match.
if(DRmin<rclus_*rclusfactor_&&jetIndexForDRmin>=0) {
jetToPartonMap[jetIndexForDRmin]=ixx;
if(ixx==0||etjet_[jetIndexForDRmin]<etmin)
etmin=etjet_[jetIndexForDRmin];
// Otherwise this parton is not matched so veto the event.
} else return true;
}
// Veto events with larger jet multiplicity from exclusive sample.
if(!highestMultiplicity_&&pjet_.size()>partonsToMatch_.size()) return true;
// Veto events where matched jets are softer than non-matched ones,
// in the inclusive (highestMultiplicity_ = true) mode, unless we
// are dealing with NLO input events.
if(highestMultiplicity_) {
for(unsigned int ixx=0; ixx<pjet_.size(); ixx++)
if(jetToPartonMap[ixx]<0&&etmin<etjet_[ixx]) return true;
}
if(!vetoHeavyQ_) {
//cout << "no heavy quark decay product veto!" << endl;
return false;
}
// **************************************************************** //
// * Now look to the non-light partons for heavy quark processes. * //
// **************************************************************** //
if( (ihrd_<=2||ihrd_==6||ihrd_==10||ihrd_==15||ihrd_==16) || (hpdetect_==true && hvqfound==true)) {
// Extract heavy quark progenitors and the radiation they
// produce and put it in the calorimeter.
caldel_hvq();
// Cluster particlesToCluster_ into jets with FastJet.
getFastJets(rclus_,etclus_,etaclmax_);
// If the radiation from the heavy quarks does not give rise
// to any jets we accept event.
if(pjet_.size() == 0) return false;
// If extra jets emerge from the jet clustering we only
// accept events where the jets formed by radiation from
// b and c quarks lies within drjmin_ of the heavy quark
// progenitor.
int nmjet(pjet_.size());
for(unsigned int ixx=0; ixx<pjet_.size(); ixx++) {
for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
if(!(abs(partonsToMatch_[jxx]->id())==4||abs(partonsToMatch_[jxx]->id())==5)) continue;
if(partonJetDeltaR(partonsToMatch_[jxx],pjet_[ixx])<drjmin_) {
nmjet--; // Decrease the number of unmatched jets.
etjet_[ixx]=0*GeV; // Set jet ET to zero to indicate it is 'matched'.
}
}
}
// If every jet matched to _at_least_one_ progenitor accept the event.
if(nmjet<=0) return false;
else {
// If unmatched jets remain, reject the event if highestMultiplicity_!=1
if(!highestMultiplicity_) return true;
else {
// If unmatched jets remain and highestMultiplicity is true then check
// that these are softer than all the matched ones (from the light-parton
// matching round).
Energy etmax(0.*GeV);
for(unsigned int ixx=0; ixx<pjet_.size(); ixx++) etmax=max(etjet_[ixx],etmax);
if(etmax>etmin) return true;
}
}
}
}
// Otherwise we accept the event ...
return false;
}
/* Function that returns the R distance
between a particle and a jet. */
double FxFxHandler::partonJetDeltaR(ThePEG::tPPtr partonptr, LorentzMomentum jetmom) const {
LorentzMomentum partonmom(partonptr->momentum());
// Calculate DY, DPhi and then DR
double DY(partonmom.eta()-jetmom.eta());
double DPhi(partonmom.phi()-jetmom.phi());
if(DPhi>M_PI) DPhi=2*M_PI-DPhi;
double DR(sqrt(sqr(DY)+sqr(DPhi)));
return DR;
}
double FxFxHandler::partonJetDeltaR(LorentzMomentum jetmom1, LorentzMomentum jetmom2) const {
// Calculate DY, DPhi and then DR
double DY(jetmom1.eta()-jetmom2.eta());
double DPhi(jetmom1.phi()-jetmom2.phi());
if(DPhi>M_PI) DPhi=2*M_PI-DPhi;
double DR(sqrt(sqr(DY)+sqr(DPhi)));
return DR;
}
// Get FastJets
void FxFxHandler::getFastJets(double rjet, Energy ejcut, double etajcut) const {
vector<fastjet::PseudoJet> particlesToCluster;
for(unsigned int ipar=0; ipar<particlesToCluster_.size(); ipar++) {
double y(particlesToCluster_[ipar]->momentum().eta());
if(y>=ycmin_&&y<=ycmax_) {
int absId(abs(particlesToCluster_[ipar]->id()));
// If it's not a lepton / top / photon it may go in the jet finder.
if(!(absId>=11&&absId<=16) && absId!=6 && absId!=22) {
// input particles into fastjet pseudojet
fastjet::PseudoJet p(particlesToCluster_[ipar]->momentum().x()/GeV,
particlesToCluster_[ipar]->momentum().y()/GeV,
particlesToCluster_[ipar]->momentum().z()/GeV,
particlesToCluster_[ipar]->momentum().e()/GeV);
p.set_user_index(ipar);
particlesToCluster.push_back(p);
}
}
}
fastjet::RecombinationScheme recombinationScheme = fastjet::E_scheme;
fastjet::Strategy strategy = fastjet::Best;
double R(rjet);
fastjet::JetDefinition theJetDefinition;
switch (jetAlgorithm_) {
case -1: theJetDefinition=fastjet::JetDefinition(fastjet::antikt_algorithm,
R,
recombinationScheme,
strategy); break;
case 0: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
R,
recombinationScheme,
strategy); break;
case 1: theJetDefinition=fastjet::JetDefinition(fastjet::kt_algorithm,
R,
recombinationScheme,
strategy); break;
default: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
R,
recombinationScheme,
strategy); break;
}
fastjet::ClusterSequence fastjetEvent(particlesToCluster,theJetDefinition);
vector<fastjet::PseudoJet> inclusiveJets = fastjetEvent.inclusive_jets();
inclusiveJets = fastjet::sorted_by_pt(inclusiveJets);
// Fill the array of jet momenta for the rest of the veto procedure.
pjet_.clear();
pjet_.resize(inclusiveJets.size());
etjet_.clear();
etjet_.resize(inclusiveJets.size());
for(unsigned int ffj=0; ffj<pjet_.size();ffj++) {
pjet_[ffj] = Lorentz5Momentum(inclusiveJets[ffj].px()*GeV,
inclusiveJets[ffj].py()*GeV,
inclusiveJets[ffj].pz()*GeV,
inclusiveJets[ffj].e()*GeV);
pjet_[ffj].rescaleMass();
etjet_[ffj] = pjet_[ffj].et();
}
// Throw the jet away if it's outside required eta region or
// has transverse energy below ejcut.
for(unsigned int fj=0; fj<pjet_.size(); fj++)
if(etjet_[fj]<ejcut||fabs(pjet_[fj].eta())>etajcut) {
pjet_.erase(pjet_.begin()+fj);
etjet_.erase(etjet_.begin()+fj);
fj--;
}
// Sort jets from high to low ET.
vector<pair<Energy, Lorentz5Momentum> > etjet_pjet;
for(unsigned int ixx=0; ixx<etjet_.size(); ixx++)
etjet_pjet.push_back(make_pair(etjet_[ixx],pjet_[ixx]));
sort(etjet_pjet.begin(),etjet_pjet.end(),ETsortFunction);
for(unsigned int ixx=0; ixx<etjet_.size(); ixx++) {
etjet_[ixx]=etjet_pjet[ixx].first;
pjet_[ixx]=etjet_pjet[ixx].second;
}
return;
}
// Get FastJets from partonsToMatch_
void FxFxHandler::getFastJetsToMatch(double rjet) const {
vector<fastjet::PseudoJet> particlesToCluster;
for(unsigned int ipar=0; ipar<partonsToMatch_.size(); ipar++) {
double y(partonsToMatch_[ipar]->momentum().eta());
if(y>=ycmin_&&y<=ycmax_) {
int absId(abs(partonsToMatch_[ipar]->id()));
// If it's not a lepton / top / photon it may go in the jet finder.
if(!(absId>=11&&absId<=16) && absId!=6 && absId!=22) {
// input particles into fastjet pseudojet
fastjet::PseudoJet p(partonsToMatch_[ipar]->momentum().x()/GeV,
partonsToMatch_[ipar]->momentum().y()/GeV,
partonsToMatch_[ipar]->momentum().z()/GeV,
partonsToMatch_[ipar]->momentum().e()/GeV);
p.set_user_index(ipar);
particlesToCluster.push_back(p);
}
}
}
fastjet::RecombinationScheme recombinationScheme = fastjet::E_scheme;
fastjet::Strategy strategy = fastjet::Best;
double R(rjet);
fastjet::JetDefinition theJetDefinition;
switch (jetAlgorithm_) {
case -1: theJetDefinition=fastjet::JetDefinition(fastjet::antikt_algorithm,
R,
recombinationScheme,
strategy); break;
case 0: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
R,
recombinationScheme,
strategy); break;
case 1: theJetDefinition=fastjet::JetDefinition(fastjet::kt_algorithm,
R,
recombinationScheme,
strategy); break;
default: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
R,
recombinationScheme,
strategy); break;
}
fastjet::ClusterSequence fastjetEvent(particlesToCluster,theJetDefinition);
vector<fastjet::PseudoJet> inclusiveJets = fastjetEvent.inclusive_jets();
inclusiveJets = fastjet::sorted_by_pt(inclusiveJets);
// Fill the array of jet momenta for the rest of the veto procedure.
pjetME_.clear();
pjetME_.resize(inclusiveJets.size());
for(unsigned int ffj=0; ffj<pjetME_.size();ffj++) {
pjetME_[ffj] = Lorentz5Momentum(inclusiveJets[ffj].px()*GeV,
inclusiveJets[ffj].py()*GeV,
inclusiveJets[ffj].pz()*GeV,
inclusiveJets[ffj].e()*GeV);
pjetME_[ffj].rescaleMass();
}
return;
}
// Deletes particles from partonsToMatch_ and particlesToCluster_
// vectors so that these contain only the partons to match to the
// jets and the particles used to build jets respectively. By and
// large the candidates for deletion are: vector bosons and their
// decay products, Higgs bosons, photons as well as _primary_, i.e.
// present in the lowest multiplicity process, heavy quarks and
// any related decay products.
void FxFxHandler::getDescendents(PPtr theParticle) const {
ParticleVector theChildren(theParticle->children());
for (unsigned int ixx=0; ixx<theChildren.size(); ixx++)
if(theChildren[ixx]->children().size()==0)
tmpList_.push_back(theChildren[ixx]);
else
getDescendents(theChildren[ixx]);
return;
}
void FxFxHandler::caldel_m() const {
preshowerFSPsToDelete_.clear();
showeredFSPsToDelete_.clear();
hvqfound = false;
if(hpdetect_ && mergemode_!=2) {
for(unsigned int ixx=0; ixx<preshowerFSPs_.size(); ixx++) {
tmpList_.clear();
/* Exclude the top quarks and any children
they may have produced from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6) {
hvqfound = true;
// cout << "preshowerFSPs_[ixx]->id() = " << preshowerFSPs_[ixx]->id() << " " << preshowerFSPs_[ixx]->momentum().perp2() << endl;
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]->parents()[0]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++) {
// cout << "tmpList_[jxx]->id() = " << tmpList_[jxx]->id() << " " << tmpList_[jxx]->momentum().perp2() << endl;
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
continue;
}
/* Exclude the v.bosons and any children they may
have produced from the jet parton matching. */
if( (abs(preshowerFSPs_[ixx]->parents()[0]->id())==23||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24||
abs(preshowerFSPs_[ixx]->id())==22||
abs(preshowerFSPs_[ixx]->id())==25|| (abs(preshowerFSPs_[ixx]->id()) < 17 && abs(preshowerFSPs_[ixx]->id()) > 10)) && (abs(preshowerFSPs_[ixx]->parents()[0]->id())!=6)) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
continue;
}
/* Exclude the bottom quarks and any children
they may have produced from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->id())==5&&ixx<2) {
hvqfound = true;
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
continue;
}
}
}
if(!hpdetect_) {
for(unsigned int ixx=0; ixx<preshowerFSPs_.size(); ixx++) {
tmpList_.clear();
if(ihrd_<=2) {
/* wqq... , zqq... */
/* Exclude the heavy quarks and any children they may
have produced as well as the v.boson and any children
it may have produced from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==23||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
if(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=4) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "wqq / zqq process should have 4 particles to omit from"
<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
}
} else if(ihrd_<=4) {
/* zjet... */
/* Exclude the v.boson and any children
it may have produced from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==23||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
/*if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=2) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "zjet process should have 2 particles to omit from"
<< "jet-parton matching." << Exception::eventerror;
}*/
} else if(ihrd_==5) {
/* vbjet... */
/* Exclude the v.bosons and any children they may
have produced from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==23||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24||
abs(preshowerFSPs_[ixx]->id())==22||
abs(preshowerFSPs_[ixx]->id())==25) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
} else if(ihrd_==6) {
/* 2Q... */
/* Exclude the heavy quarks and any children
they may have produced from the jet parton matching. */
if(ihvy_==6) {
if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
// cout << "preshowerFSPs_[ixx]->id() = " << preshowerFSPs_[ixx]->id() << " " << preshowerFSPs_[ixx]->momentum().perp2() << endl;
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
}
if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6) {
getDescendents(preshowerFSPs_[ixx]->parents()[0]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++) {
// cout << "tmpList_[jxx]->id() = " << tmpList_[jxx]->id() << " " << tmpList_[jxx]->momentum().perp2() << endl;
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
}
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=6) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "2Q process should have 6 particles to omit from"
<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
}
} else {
if(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=2) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "2Q process should have 2 particles to omit from"
<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
}
}
} else if(ihrd_==7) {
/* 4Q... */
/* There are no light jets for this process, so nothing to match. */
} else if(ihrd_==9) {
/* Njet... */
} else if(ihrd_==10) {
/* wcjet... */
/* Exclude the charm quark and any children it may
have produced as well as the v.boson and any children
it may have produced from the jet parton matching. */
if((abs(preshowerFSPs_[ixx]->id())==4&&ixx<1)||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=3) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "wcjet process should have 3 particles to omit from"
<< "jet-parton matching." << Exception::eventerror;
}
} else if(ihrd_==11) {
/* phjet... */
/* Exclude the hard photons from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->id())==22) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
unsigned int tmpUnsignedInt(nph_);
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "phjet process should have " << nph_ << " particles to omit from"
<< "jet-parton matching." << Exception::eventerror;
}
} else if(ihrd_==12) {
/* hjet... */
/* Exclude the higgs and any children it may have
produced from the jet parton matching. */
if(abs(preshowerFSPs_[ixx]->id())==25) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=1) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "hjet process should have 1 particle to omit from"
<< "jet-parton matching." << Exception::eventerror;
}
} else if(ihrd_==14) {
/* wphjet... */
/* Exclude the v.boson and any children it may have
produced from the jet parton matching. */
// AND WHAT ABOUT THE PHOTON? <--- CHECK THIS WITH AlpGen GUYs.
if(abs(preshowerFSPs_[ixx]->id())==22||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
unsigned int tmpUnsignedInt(2+nph_);
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "wphjet process should have " << 2+nph_ << " particles to omit from"
<< "jet-parton matching." << Exception::eventerror;
}
} else if(ihrd_==15) {
/* wphqq... <--- N.B. if q = top, it is not decayed. */
/* Exclude the heavy quarks and any children they may
have produced as well as the v.boson and any children
it may have produced from the jet parton matching. */
// AND WHAT ABOUT THE PHOTON? <--- CHECK THIS WITH AlpGen GUYs.
if(abs(preshowerFSPs_[ixx]->id())==22||
(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx==(preshowerFSPs_.size()-(2+nph_+1)))||
(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx==(preshowerFSPs_.size()-(2+nph_+2)))||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
unsigned int tmpUnsignedInt(4+nph_);
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "wphjet process should have " << 4+nph_ << " particles to omit from"
<< "jet-parton matching." << Exception::eventerror;
}
} else if(ihrd_==16) {
/* 2Qph... <--- if Q is a top it will decay. */
/* Exclude the hard photons and any children they
may have produced from the jet parton matching
as well as the heavy quarks and any children it
may have produced. */
// AND WHAT ABOUT THE PHOTON?
if(ihvy_==6) {
if(abs(preshowerFSPs_[ixx]->id())==22||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==6||
abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
unsigned int tmpUnsignedInt(6+nph_);
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "wphjet process should have " << 6+nph_ << " particles to omit from"
<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
}
} else {
if(abs(preshowerFSPs_[ixx]->id())==22||
(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2)) {
preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
getDescendents(preshowerFSPs_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
unsigned int tmpUnsignedInt(2+nph_);
if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "wphjet process should have " << 2+nph_ << " particles to omit from"
<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
}
}
}
}
}
// cout << "partonsToMatch_.size()= " << partonsToMatch_.size() << " preshowerFSPsToDelete_.size() = " << preshowerFSPsToDelete_.size() << endl;
// cout << "deleting preshowerFSPs" << endl;
for(unsigned int ixx=0; ixx<preshowerFSPsToDelete_.size(); ixx++) {
for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
if(preshowerFSPsToDelete_[ixx]==partonsToMatch_[jxx]) {
// cout << "deleting " << preshowerFSPsToDelete_[ixx]->id() << " with parent " << preshowerFSPsToDelete_[ixx]->parents()[0]->id() << endl; // " energy = " << preshowerFSPsToDelete_[ixx]->momentum().e()*GeV << endl;
partonsToMatch_.erase(partonsToMatch_.begin()+jxx);
break;
}
}
}
//cout << "partonsToMatch_.size() (AFTER) = " << partonsToMatch_.size() << endl;
//cout << "deleting showeredFSPs" << endl;
for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++) {
for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
if(showeredFSPsToDelete_[ixx]==particlesToCluster_[jxx]) {
// cout << "deleting " << showeredFSPsToDelete_[ixx]->id() << " with parent " << showeredFSPsToDelete_[ixx]->parents()[0]->id() << endl; //" energy = " << preshowerFSPsToDelete_[ixx]->momentum().e()*GeV << endl;
particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
break;
}
}
}
// Sanity check!
if(partonsToMatch_.size()>particlesToCluster_.size()) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "No. of ME level partons to be matched to jets = "
<< partonsToMatch_.size() << "\n"
<< "No. of showered particles to build jets from = "
<< particlesToCluster_.size() << "\n"
<< "There should be at least as many partons to\n"
<< "cluster as there are partons to match to.\n"
<< Exception::eventerror;
}
// cout << "partonsToMatch_.size() (AFTER2) = " << partonsToMatch_.size() << endl;
// Acid test.
/* unsigned int tmpUnsignedInt(njets_);
if(!inputIsNLO_&&partonsToMatch_.size()!=tmpUnsignedInt) {
for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
if(abs(partonsToMatch_[ixx]->id())>=6&&
abs(partonsToMatch_[ixx]->id())!=21)
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "Found a parton to match to which is not a quark or gluon!"
<< *partonsToMatch_[ixx] << "\n"
<< Exception::eventerror;
}
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "No. of ME level partons to be matched to jets = "
<< partonsToMatch_.size() << "\n"
<< "No. of light jets (njets) in AlpGen process = "
<< njets_ << "\n"
<< "These should be equal." << "\n"
<< Exception::eventerror;
} */
//cout << "partonsToMatch_.size() (AFTER3) = " << partonsToMatch_.size() << endl;
return;
}
// This looks for all descendents of a top up to but not including
// the W and b children.
void FxFxHandler::getTopRadiation(PPtr theParticle) const {
ParticleVector theChildren(theParticle->children());
for (unsigned int ixx=0; ixx<theChildren.size(); ixx++)
if(theChildren[ixx]->children().size()==0)
tmpList_.push_back(theChildren[ixx]);
else if(abs(theChildren[ixx]->id())==5||abs(theChildren[ixx]->id())==24)
return;
else
getTopRadiation(theChildren[ixx]);
return;
}
void FxFxHandler::caldel_mg() const {
preshowerFSPsToDelete_.clear();
showeredFSPsToDelete_.clear();
/*
* Get the MadGraph clustering information
* from the Les Houches event tags
*/
ptclust_.clear();
getptclust();
//for(unsigned int izz = 0; izz<ptclust_.size(); izz++) cout << "ptclust_[" << izz << "] = " << ptclust_[izz] << endl;
/*
* Get the COM energy of the colliding hadrons
*
*/
getECOM();
// cout << "ECOM_ = " << ECOM_ << endl;
for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
if(ptclust_[ixx] == ECOM_){
preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
tmpList_.clear();
getDescendents(partonsToMatch_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
}
for(unsigned int ixx=0; ixx<preshowerFSPsToDelete_.size(); ixx++) {
for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
if(preshowerFSPsToDelete_[ixx]==partonsToMatch_[jxx]) {
//cout << "deleting " << preshowerFSPsToDelete_[ixx]->id() << endl;
partonsToMatch_.erase(partonsToMatch_.begin()+jxx);
break;
}
}
}
for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++) {
for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
if(showeredFSPsToDelete_[ixx]==particlesToCluster_[jxx]) {
// cout << "deleting " << showeredFSPsToDelete_[ixx]->id() << " with parent " << showeredFSPsToDelete_[ixx]->parents()[0]->id() << endl; //" energy = " << preshowerFSPsToDelete_[ixx]->momentum().e()*GeV << endl;
particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
break;
}
}
}
// Sanity check!
if(partonsToMatch_.size()>particlesToCluster_.size()) {
throw Exception()
<< "FxFxHandler::caldel_m() - ERROR!\n"
<< "No. of ME level partons to be matched to jets = "
<< partonsToMatch_.size() << "\n"
<< "No. of showered particles to build jets from = "
<< particlesToCluster_.size() << "\n"
<< "There should be at least as many partons to\n"
<< "cluster as there are partons to match to.\n"
<< Exception::eventerror;
}
}
void FxFxHandler::caldel_hvq() const {
// Fill partonsToMatch_ with only those pre-shower partons intended to
// be used in heavy-quark-jet matching and fill particlesToCluster_ using
// only those final state particles (post-shower) which are supposed
// in the heavy-quark-jet clustering used to do merging. To begin with
// these are made from the corresponding sets of particles that were
// omitted from the initial jet-parton matching run.
partonsToMatch_ = preshowerFSPsToDelete_;
particlesToCluster_.resize(showeredFSPsToDelete_.size());
for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++)
particlesToCluster_[ixx] = showeredFSPsToDelete_[ixx];
// Reset the arrays of particles to delete so that the partonsToMatch_
// and particlesToCluster_ vectors can undergo further filtering.
preshowerFSPsToDelete_.clear();
showeredFSPsToDelete_.clear();
// Determine further particles in partonsToMatch_ and particlesToCluster_
// for deletion.
for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
// If the progenitor particle is not a heavy
// quark we delete it and its descendents.
if(abs(partonsToMatch_[ixx]->id())<4||abs(partonsToMatch_[ixx]->id())>6) {
preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
tmpList_.clear();
getDescendents(partonsToMatch_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
// If the progenitor is a b quark from a top decay drop
// it & it's descendents too!
} else if(abs(partonsToMatch_[ixx]->id())==5&&
partonsToMatch_[ixx]->parents().size()>0&&
abs(partonsToMatch_[ixx]->parents()[0]->id())==6) {
preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
tmpList_.clear();
getDescendents(partonsToMatch_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
// If (it's a hvy quark not from a top decay and) it has a W/Z/H
// as a parent [ditto].
} else if(partonsToMatch_[ixx]->parents().size()>0&&
(abs(partonsToMatch_[ixx]->parents()[0]->id())==23||
abs(partonsToMatch_[ixx]->parents()[0]->id())==24||
abs(partonsToMatch_[ixx]->parents()[0]->id())==25)) {
preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
tmpList_.clear();
getDescendents(partonsToMatch_[ixx]);
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
showeredFSPsToDelete_.push_back(tmpList_[jxx]);
}
}
// Now do the necessary deleting from partonsToMatch_ and
// particlesToCluster_.
for(unsigned int ixx=0; ixx<preshowerFSPsToDelete_.size(); ixx++) {
for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
if(preshowerFSPsToDelete_[ixx]==partonsToMatch_[jxx]) {
partonsToMatch_.erase(partonsToMatch_.begin()+jxx);
break;
}
}
}
for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++) {
for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
if(showeredFSPsToDelete_[ixx]==particlesToCluster_[jxx]) {
particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
break;
}
}
}
// Now we return to get the decaying top quarks and any
// radiation they produced.
ParticleVector intermediates(lastXCombPtr()->subProcess()->intermediates());
for(unsigned int ixx=0; ixx<intermediates.size(); ixx++) {
if(abs(intermediates[ixx]->id())==6) {
partonsToMatch_.push_back(intermediates[ixx]);
tmpList_.clear();
getTopRadiation(partonsToMatch_.back());
for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
particlesToCluster_.push_back(tmpList_[jxx]);
}
}
// If there are any heavy quark progenitors we have to remove
// the final (showered) instance of them from particlesToCluster.
ParticleVector evolvedHeavyQuarks;
for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
if(abs(partonsToMatch_[ixx]->id())>=4&&abs(partonsToMatch_[ixx]->id())<=6) {
theProgenitor = partonsToMatch_[ixx];
// Follow the heavy quark line down to where it stops branching.
while(theProgenitor->children().size()>0) {
theLastProgenitor = theProgenitor;
for(unsigned int jxx=0; jxx<theProgenitor->children().size(); jxx++) {
if(theProgenitor->children()[jxx]->id()==theProgenitor->id())
theProgenitor=theProgenitor->children()[jxx];
}
// If the progenitor had children but none of them had
// the same particle id as it, then it must have undergone
// a decay rather than a branching, i.e. it is the end of
// the evolution line, so,
if(theProgenitor==theLastProgenitor) break;
}
evolvedHeavyQuarks.push_back(theProgenitor);
}
}
// Now delete the evolved heavy quark from the particlesToCluster.
for(unsigned int ixx=0; ixx<evolvedHeavyQuarks.size(); ixx++) {
for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
if(evolvedHeavyQuarks[ixx]==particlesToCluster_[jxx]) {
particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
break;
}
}
}
return;
}
// get npLO_ and npNLO_
void FxFxHandler::getnpFxFx() const {
split_vector_type SplitVec;
// pull the optional weights from the current event
map<string,double> optionalEventWeights = eventHandler()->currentEvent()->optionalWeights();
// loop over the optional weights and find np values
for (map<string,double>::const_iterator it=optionalEventWeights.begin(); it!=optionalEventWeights.end(); ++it){
// split the line
boost::split( SplitVec, it->first, boost::is_any_of(" ") );
// if np is found, store the information
if(SplitVec[0] == "np") {
npLO_ = atof(SplitVec[1].c_str());
npNLO_ = atof(SplitVec[2].c_str());
}
}
return;
}
// get hadron COM energy
void FxFxHandler::getECOM() const {
split_vector_type SplitVec;
// pull the optional weights from the current event
map<string,double> optionalEventWeights = eventHandler()->currentEvent()->optionalWeights();
// loop over the optional weights and find np values
for (map<string,double>::const_iterator it=optionalEventWeights.begin(); it!=optionalEventWeights.end(); ++it){
// split the line
boost::split( SplitVec, it->first, boost::is_any_of(" ") );
// if np is found, store the information
if(SplitVec[0] == "ecom") {
ECOM_ = it->second;
}
}
return;
}
// get pt_clust information
void FxFxHandler::getptclust() const {
split_vector_type SplitVec;
// pull the optional weights from the current event
map<string,double> optionalEventWeights = eventHandler()->currentEvent()->optionalWeights();
// loop over the optional weights and find np values
string str_eq = "=";
string str_quote = "\"";
int countptc_(0);
for (map<string,double>::const_iterator it=optionalEventWeights.begin(); it!=optionalEventWeights.end(); ++it){
// split the line
double wgtid = it->second;
//cout << "wgtdid = " << wgtid << endl;
if(wgtid == -333) {
// cout << "it->first for -333 = " << it->first << endl;
boost::split( SplitVec, it->first, boost::is_any_of(" ") );
// if np is found, store the information
double ptclust(0.);
string stringtohandle("");
for(unsigned int i_ = 0; i_ < SplitVec.size(); ++i_) {
if(SplitVec[i_].find("pt_clust_") != std::string::npos) {
//cout << "pt_clust_ found in " << SplitVec[i_] << endl;
countptc_++;
// cout << "SplitVec[i_] = " << SplitVec[i_] << endl;
stringtohandle = SplitVec[i_];
stringtohandle.erase(0,10);
erase_substr(stringtohandle,str_eq);
erase_substr(stringtohandle,str_quote);
// cout << "stringtohandle = " << stringtohandle << endl;
ptclust = atof(stringtohandle.c_str());
ptclust_.push_back(ptclust);
}
}
}
}
return;
}
void FxFxHandler::getPreshowerParticles() const {
// LH file initial-state partons:
preshowerISPs_ = lastXCombPtr()->subProcess()->incoming();
// LH file final-state partICLEs:
preshowerFSPs_ = lastXCombPtr()->subProcess()->outgoing();
return;
}
void FxFxHandler::getShoweredParticles() const {
// Post-shower initial-state hadrons:
showeredISHs_ = eventHandler()->currentEvent()->incoming();
// Post-shower initial-state partons:
for(unsigned int ixx=0; ixx<(showeredISHs_.first)->children().size(); ixx++)
if(((showeredISHs_.first)->children()[ixx]->id())<6||
((showeredISHs_.first)->children()[ixx]->id())==21)
showeredISPs_.first=(showeredISHs_.first)->children()[ixx];
for(unsigned int ixx=0; ixx<(showeredISHs_.second)->children().size(); ixx++)
if(((showeredISHs_.second)->children()[ixx]->id())<6||
((showeredISHs_.second)->children()[ixx]->id())==21)
showeredISPs_.second=(showeredISHs_.second)->children()[ixx];
// Post-shower final-state partICLEs plus remnants (to be removed later):
showeredFSPs_ = eventHandler()->currentEvent()->getFinalState();
// Post-shower final-state remnants:
for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++) {
if(showeredFSPs_[ixx]->PDGName()=="Rem:p+"||
showeredFSPs_[ixx]->PDGName()=="Rem:pbar-") {
if(showeredFSPs_[ixx]->parents()[0]->parents()[0]==
showeredISHs_.first)
showeredRems_.first=showeredFSPs_[ixx];
else if(showeredFSPs_[ixx]->parents()[0]->parents()[0]==
showeredISHs_.second)
showeredRems_.second=showeredFSPs_[ixx];
}
}
// Now delete found remnants from the showeredFSPs vector for consistency.
for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
if(showeredFSPs_[ixx]->PDGName()=="Rem:p+")
showeredFSPs_.erase(showeredFSPs_.begin()+ixx);
for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
if(showeredFSPs_[ixx]->PDGName()=="Rem:pbar-")
showeredFSPs_.erase(showeredFSPs_.begin()+ixx);
sort(showeredFSPs_.begin(),showeredFSPs_.end(),recordEntry);
return;
}
void FxFxHandler::doSanityChecks(int debugLevel) const {
// When checking momentum conservation in the form
// p_in - p_out, any momentum component bigger / less
// than + / - epsilon will result in the p_in - p_out
// vector being flagged as "non-null", triggering a
// warning that momentum conservation is violated.
Energy epsilon(0.5*GeV);
if(debugLevel>=5) epsilon=1e-9*GeV;
// Print out what was found for the incoming and outgoing
// partons in the lastXCombPtr regardless.
if(debugLevel>=5) {
cout << "\n\n\n\n";
cout << "****************************************************" << endl;
cout << " The following are the hard subprocess momenta from " << "\n"
<< " lastXCombPtr and should be basically identical to " << "\n"
<< " the input LH file momenta." << "\n\n";
cout << " Incoming particles:" << "\n"
<< *(preshowerISPs_.first) << "\n"
<< *(preshowerISPs_.second) << endl;
cout << " Outgoing particles:" << endl;
for(unsigned int ixx=0; ixx<preshowerFSPs_.size(); ixx++)
cout << *(preshowerFSPs_[ixx]) << endl;
}
// Print out what was found for the incoming and outgoing
// partons after the shower.
if(debugLevel>=5) {
cout << "\n\n";
cout << "****************************************************" << endl;
cout << " The following are the particles left at the end of" << "\n"
<< " the showering step." << "\n\n";
cout << " Incoming hadrons:" << "\n"
<< *(showeredISHs_.first) << "\n"
<< *(showeredISHs_.second) << endl;
cout << " Incoming partons:" << "\n"
<< *(showeredISPs_.first) << "\n"
<< *(showeredISPs_.second) << endl;
cout << " Outgoing partons:" << endl;
for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
cout << *(showeredFSPs_[ixx]) << endl;
cout << " Outgoing remnants:" << "\n"
<< *(showeredRems_.first) << "\n"
<< *(showeredRems_.second) << endl;
}
// Check if we correctly identified all initial and final-state
// particles by testing momentum is conserved.
if(debugLevel>=4) {
Lorentz5Momentum tmpMom;
tmpMom += showeredISPs_.first->momentum();
tmpMom += showeredISPs_.second->momentum();
for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
tmpMom -= showeredFSPs_[ixx]->momentum();
if(!isMomLessThanEpsilon(tmpMom,epsilon))
cout << "Total parton mom.in - total parton mom.out = "
<< tmpMom/GeV << endl;
tmpMom = showeredISHs_.first->momentum()
- showeredRems_.first->momentum() -showeredISPs_.first->momentum();
if(!isMomLessThanEpsilon(tmpMom,epsilon))
cout << "First p_hadron-p_remnant-p_incoming " << tmpMom/GeV << endl;
tmpMom = showeredISHs_.second->momentum()
- showeredRems_.second->momentum()-showeredISPs_.second->momentum();
if(!isMomLessThanEpsilon(tmpMom,epsilon))
cout << "Second p_hadron-p_remnant-p_incoming " << tmpMom/GeV << endl;
}
// Check if what we found to be the remnant is consistent with
// what we identified as the parent incoming hadron i.e. p+
// goes with Rem:p+ and pbar- goes with Rem:pbar-.
if(debugLevel>=0) {
string tmpString;
tmpString=showeredRems_.first->PDGName();
tmpString=tmpString.substr(tmpString.find_first_of(":")+1,
string::npos);
if(showeredISHs_.first->PDGName()!=tmpString) {
cout << "FxFxHandler::showerHardProcessVeto" << "\n"
<< "Fatal error in pairing of remnant and parent hadron." << "\n"
<< "Remnant = " << *(showeredRems_.first) << "\n"
<< "Parent hadron = " << *(showeredISHs_.first)
<< endl;
cout << showeredISHs_.first->PDGName() << endl;
cout << tmpString << endl;
}
tmpString=showeredRems_.second->PDGName();
tmpString=tmpString.substr(tmpString.find_first_of(":")+1,
string::npos);
if(showeredISHs_.second->PDGName()!=tmpString) {
cout << "FxFxHandler::showerHardProcessVeto" << "\n"
<< "Fatal error in pairing of remnant and parent hadron." << "\n"
<< "Remnant = " << *(showeredRems_.second) << "\n"
<< "Parent hadron = " << *(showeredISHs_.second)
<< endl;
cout << showeredISHs_.second->PDGName() << endl;
cout << tmpString << endl;
}
}
return;
}
void FxFxHandler::printMomVec(vector<Lorentz5Momentum> momVec) {
cout << "\n\n";
// Label columns.
printf("%5s %9s %9s %9s %9s %9s %9s %9s %9s %9s\n",
"jet #",
"px","py","pz","E",
"eta","phi","pt","et","mass");
// Print out the details for each jet
for (unsigned int ixx=0; ixx<momVec.size(); ixx++) {
printf("%5u %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f\n",
ixx,
double(momVec[ixx].x()/GeV),double(momVec[ixx].y()/GeV),
double(momVec[ixx].z()/GeV),double(momVec[ixx].t()/GeV),
double(momVec[ixx].eta()) ,double(momVec[ixx].phi()),
double(momVec[ixx].perp()/GeV),
double(momVec[ixx].et()/GeV),
double(momVec[ixx].m()/GeV));
}
}
Energy FxFxHandler::etclusran_(double petc) const {
return (((2 * epsetclus_)/M_PI) * asin(2 * petc - 1) + etclusmean_);
}
void FxFxHandler::erase_substr(std::string& subject, const std::string& search) const {
size_t pos = 0;
while((pos = subject.find(search, pos)) != std::string::npos) {
subject.erase( pos, search.length() );
}
}
diff --git a/MatrixElement/FxFx/FxFxReader.h b/MatrixElement/FxFx/FxFxReader.h
--- a/MatrixElement/FxFx/FxFxReader.h
+++ b/MatrixElement/FxFx/FxFxReader.h
@@ -1,1016 +1,1016 @@
// -*- C++ -*-
//
// FxFxReader.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_FxFxReader_H
#define THEPEG_FxFxReader_H
// This is the declaration of the FxFxReader class.
#include "FxFx.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Utilities/ObjectIndexer.h"
#include "ThePEG/Utilities/Exception.h"
#include "ThePEG/Utilities/XSecStat.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "ThePEG/PDF/PartonBin.fh"
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "FxFxEventHandler.fh"
#include "FxFxReader.fh"
#include "ThePEG/Utilities/CFile.h"
#include <cstdio>
#include <cstring>
namespace ThePEG {
/**
* FxFxReader is an abstract base class to be used for objects
* which reads event files or streams from matrix element
* generators. Derived classes must at least implement the open() and
* doReadEvent() methods to read in information about the whole run into
* the HEPRUP variable and next event into the HEPEUP variable
* respectively. Also the close() function to close the file or stream
* read must be implemented. Although these functions are named as if
* we are reading from event files, they could just as well implement
* the actual generation of events.
*
* After filling the HEPRUP and HEPEUP variables, which are protected
* and easily accesible from the sub-class, this base class will then
* be responsible for transforming this data to the ThePEG Event
* record in the getEvent() method. <code>FxFxReader</code>s can
* only be used inside FxFxEventHandler objects.
*
* In the initialization the virtual open() and scan() functions are
* called. Here the derived class must provide the information about
* the processes in the variables corresponding to the HEPRUP common
* block. Note that the IDWTUP is required to be +/- 1, and sub
* classes are required to change the information accordingly to
* ensure the correct corss section sampling. Note also that the
* controlling FxFxEventHandler may choose to generate weighted
* events even if IDWTUP is 1.
*
* Note that the information given per process in e.g. the XSECUP and
* XMAXUP vectors is not used by the FxFxEventHandler and by
* default the FxFxReader is not assumed to be able to actively
* choose between the sub-processes. Instead, the
* FxFxEventHandler can handle several FxFxReader objects
* and choose between them. However, a sub-class of FxFxReader
* may set the flag isActive, in which case it is assumed to be able
* to select between its sub-processes itself.
*
* The FxFxReader may be assigned a number ReweightBase objects
* which either completely reweights the events produced (in the
* reweights vector), or only biases the selection without influencing
* the cross section (in the preweights vector). Note that it is the
* responsibility of a sub-class to call the reweight() function and
* multiply the weight according to its return value (typically done
* in the readEvent() function).
*
* @see \ref FxFxReaderInterfaces "The interfaces"
* defined for FxFxReader.
* @see Event
* @see FxFxEventHandler
*/
class FxFxReader: public HandlerBase, public LastXCombInfo<> {
/**
* FxFxEventHandler should have access to our private parts.
*/
friend class FxFxEventHandler;
/**
* Map for accumulating statistics of cross sections per process
* number.
*/
typedef map<int,XSecStat> StatMap;
/**
* Map of XComb objects describing the incoming partons indexed by
* the corresponding PartonBin pair.
*/
typedef map<tcPBPair,XCombPtr> XCombMap;
/**
* A vector of pointers to ReweightBase objects.
*/
typedef vector<ReweightPtr> ReweightVector;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor. If the optional argument is true, the reader
* is assumed to be able to produce events on demand for a given
* process.
*/
FxFxReader(bool active = false);
/**
* Copy-constructor.
*/
FxFxReader(const FxFxReader &);
/**
* Destructor.
*/
virtual ~FxFxReader();
//@}
public:
/** @name Main virtual fuctions to be overridden in
* sub-classes. They are named as if we are reading from event
* files, but could equally well implement the actual generation of
* events. */
//@{
/**
* Open a file or stream with events and read in the run information
* into the heprup variable.
*/
virtual void open() = 0;
/**
* Read the next event from the file or stream into the
* corresponding protected variables. Return false if there is no
* more events.
*/
virtual bool doReadEvent() = 0;
/**
* Close the file or stream from which events have been read.
*/
virtual void close() = 0;
/**
* return the weight names
*/
// virtual vector<string> optWeightsNamesFunc();
virtual vector<string> optWeightsNamesFunc() = 0;
//virtual vector<string*> optWeightNamesFunc() = 0;
vector<string> optionalWeightsNames;
/**
* The ID (e.g. 100x, 2001) for the weight
*/
// vector<string> optionalWeightsNames;
//@}
/** @name Other important function which may be overridden in
* sub-classes which wants to bypass the basic HEPRUP or HEPEUP
* variables or otherwise facilitate the conversion to ThePEG
* objects. */
//@{
/**
* Initialize. This function is called by the FxFxEventHandler
* to which this object is assigned.
*/
virtual void initialize(FxFxEventHandler & eh);
/**
* Calls readEvent() or uncacheEvent() to read information into the
* FxFx common block variables. This function is called by the
* FxFxEventHandler if this reader has been selectod to
* produce an event.
*
* @return the weight asociated with this event. If negative weights
* are allowed it should be between -1 and 1, otherwise between 0
* and 1. If outside these limits the previously estimated maximum
* is violated. Note that the estimated maximum then should be
* updated from the outside.
*/
virtual double getEvent();
/**
* Calls doReadEvent() and performs pre-defined reweightings. A
* sub-class overrides this function it must make sure that the
* corresponding reweightings are done.
*/
virtual bool readEvent();
/**
* Skip \a n events. Used by FxFxEventHandler to make sure
* that a file is scanned an even number of times in case the events
* are not ramdomly distributed in the file.
*/
virtual void skip(long n);
/**
* Get an XComb object. Converts the information in the Les Houches
* common block variables to an XComb object describing the sub
* process. This is the way information is conveyed from the reader
* to the controlling FxFxEventHandler.
*/
tXCombPtr getXComb();
/**
* Get a SubProcess object corresponding to the information in the
* Les Houches common block variables.
*/
tSubProPtr getSubProcess();
/**
* Scan the file or stream to obtain information about cross section
* weights and particles etc. This function should fill the
* variables corresponding to the /HEPRUP/ common block. The
* function returns the number of events scanned.
*/
virtual long scan();
/**
* Take the information corresponding to the HEPRUP common block and
* initialize the statistics for this reader.
*/
virtual void initStat();
/**
* Reweights the current event using the reweights and preweights
* vectors. It is the responsibility of the sub-class to call this
* function after the HEPEUP information has been retrieved.
*/
double reweight();
/**
* Converts the information in the Les Houches common block
* variables into a Particle objects.
*/
virtual void fillEvent();
/**
* Removes the particles created in the last generated event,
* preparing to produce a new one.
*/
void reset();
/**
* Possibility for subclasses to recover from non-conformant
* settings of XMAXUP when an event file has been scanned with \a
* neve events. Should set weightScale so that the average XMAXUP
* times weightScale gives the cross section for a process. (This is
* needed for MadEvent).
*/
virtual void setWeightScale(long neve);
//@}
/** @name Access information about the current event. */
//@{
/**
* Return the size of this event in bytes. To be used for the cache
* file. \a npart is the number of particles. If \a npart is 0, the
* number is taken from NUP.
*/
static size_t eventSize(int N) {
return (N + 1)*sizeof(int) + // IDPRUP, ISTUP
(7*N + 4)*sizeof(double) + // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP,
// VTIMUP, SPINUP
N*sizeof(long) + // IDUP
2*N*sizeof(pair<int,int>) + // MOTHUP, ICOLUP
sizeof(pair<double,double>) + // XPDWUP.
2*sizeof(double); // lastweight and preweight
}
/**
* The current event weight given by XWGTUP times possible
* reweighting. Note that this is not necessarily the same as what
* is returned by getEvent(), which is scaled with the maximum
* weight.
*/
double eventWeight() const { return hepeup.XWGTUP*lastweight; }
/**
* Return the optional named weights associated to the current event.
*/
const map<string,double>& optionalEventWeights() const { return optionalWeights; }
/**
* Return the Les Houches event number associated with the current event
*/
const long& LHEEventNum() const { return LHEeventnum; }
/**
* Return the optional npLO and npNLO
*/
const int& optionalEventnpLO() const { return optionalnpLO; }
const int& optionalEventnpNLO() const { return optionalnpNLO; }
/**
* The pair of PartonBinInstance objects describing the current
* incoming partons in the event.
*/
const PBIPair & partonBinInstances() const { return thePartonBinInstances; }
/**
* Return the instances of the beam particles for the current event.
*/
const PPair & beams() const { return theBeams; }
/**
* Return the instances of the incoming particles to the sub process
* for the current event.
*/
const PPair & incoming() const { return theIncoming; }
/**
* Return the instances of the outgoing particles from the sub process
* for the current event.
*/
const PVector & outgoing() const { return theOutgoing; }
/**
* Return the instances of the intermediate particles in the sub
* process for the current event.
*/
const PVector & intermediates() const { return theIntermediates; }
/**
* If this reader is to be used (possibly together with others) for
* CKKW reweighting and veto, this should give the multiplicity of
* outgoing particles in the highest multiplicity matrix element in
* the group.
*/
int maxMultCKKW() const { return theMaxMultCKKW; }
/**
* If this reader is to be used (possibly together with others) for
* CKKW reweighting and veto, this should give the multiplicity of
* outgoing particles in the lowest multiplicity matrix element in
* the group.
*/
int minMultCKKW() const { return theMinMultCKKW; } //@}
/** @name Other inlined access functions. */
//@{
/**
* The number of events found in this reader. If less than zero the
* number of events are unlimited.
*/
long NEvents() const { return theNEvents; }
/**
* The number of events produced so far. Is reset to zero if an
* event file is reopened.
*/
long currentPosition() const { return position; }
/**
* The maximum number of events to scan to collect information about
* processes and cross sections. If less than 0, all events will be
* scanned.
*/
long maxScan() const { return theMaxScan; }
/**
* Return true if this reader is active.
*/
bool active() const { return isActive; }
/**
* True if negative weights may be produced.
*/
bool negativeWeights() const { return heprup.IDWTUP < 0; }
/**
* The collected cross section statistics for this reader.
*/
const XSecStat & xSecStats() const { return stats; }
/**
* Collected statistics about the individual processes.
*/
const StatMap & processStats() const { return statmap; }
/**
* Select the current event. It will later be rejected with a
* probability given by \a weight.
*/
void select(double weight) {
stats.select(weight);
statmap[hepeup.IDPRUP].select(weight);
}
/**
* Accept the current event assuming it was previously selcted.
*/
void accept() {
stats.accept();
statmap[hepeup.IDPRUP].accept();
}
/**
* Reject the current event assuming it was previously accepted.
*/
void reject(double w) {
stats.reject(w);
statmap[hepeup.IDPRUP].reject(w);
}
/**
* Increase the overestimated cross section for this reader.
*/
virtual void increaseMaxXSec(CrossSection maxxsec);
/**
* The PartonExtractor object used to construct remnants.
*/
tPExtrPtr partonExtractor() const { return thePartonExtractor; }
/**
* Return a possibly null pointer to a CascadeHandler to be used for
* CKKW-reweighting.
*/
tCascHdlPtr CKKWHandler() const { return theCKKW; }
/**
* The pairs of PartonBin objects describing the partons which can
* be extracted by the PartonExtractor object.
*/
const PartonPairVec & partonBins() const { return thePartonBins; }
/**
* The map of XComb objects indexed by the corresponding PartonBin
* pair.
*/
const XCombMap & xCombs() const { return theXCombs; }
/**
* The Cuts object to be used for this reader.
*/
const Cuts & cuts() const { return *theCuts; }
//@}
protected:
/** @name Functions for manipulating cache files. */
//@{
/**
* Name of file used to cache the events form the reader in a
* fast-readable form. If empty, no cache file will be generated.
*/
string cacheFileName() const { return theCacheFileName; }
/**
* Determines whether to apply cuts to events converting them to
* ThePEG format.
*/
bool cutEarly() const { return doCutEarly; }
/**
* File stream for the cache.
*/
CFile cacheFile() const { return theCacheFile;}
/**
* Open the cache file for reading.
*/
void openReadCacheFile();
/**
* Open the cache file for writing.
*/
void openWriteCacheFile();
/**
* Close the cache file;
*/
void closeCacheFile();
/**
* Write the current event to the cache file.
*/
void cacheEvent() const;
/**
* Read an event from the cache file. Return false if something went wrong.
*/
bool uncacheEvent();
/**
* Reopen a reader. If we have reached the end of an event file,
* reopen it and issue a warning if we have used up a large fraction
* of it.
*/
void reopen();
/**
* Helper function to write a variable to a memory location
*/
template <typename T>
static char * mwrite(char * pos, const T & t, size_t n = 1) {
std::memcpy(pos, &t, n*sizeof(T));
return pos + n*sizeof(T);
}
/**
* Helper function to read a variable from a memory location
*/
template <typename T>
static const char * mread(const char * pos, T & t, size_t n = 1) {
- std::memcpy(&t, pos, n*sizeof(T));
+ std::memcpy((void *)&t, pos, n*sizeof(T));
return pos + n*sizeof(T);
}
//@}
/** @name Auxilliary virtual methods which may be verridden by sub-classes. */
//@{
/**
* Check the existence of a pair of PartonBin objects corresponding
* to the current event.
*
* @return false if no pair of suitable PartonBin objects was found.
*/
virtual bool checkPartonBin();
/**
* Create instances of all particles in the event and store them
* in particleIndex.
*/
virtual void createParticles();
/**
* Using the already created particles create a pair of
* PartonBinInstance objects corresponding to the incoming
* partons. Return the corresponding PartonBin objects.
*/
virtual tcPBPair createPartonBinInstances();
/**
* Create instances of the incoming beams in the event and store
* them in particleIndex. If no beam particles are included in the
* event they are created from the run info.
*/
virtual void createBeams();
/**
* Go through the mother indices and connect up the Particles.
*/
virtual void connectMothers();
//@}
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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
protected:
/** @name Set functions for some variables not in the Les Houches accord. */
//@{
/**
* The number of events in this reader. If less than zero the number
* of events is unlimited.
*/
void NEvents(long x) { theNEvents = x; }
/**
* The map of XComb objects indexed by the corresponding PartonBin
* pair.
*/
XCombMap & xCombs() { return theXCombs; }
//@}
/** @name Standard (and non-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();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish() {
close();
HandlerBase::dofinish();
}
/**
* Return true if this object needs to be initialized before all
* other objects because it needs to extract PDFs from the event file.
*/
virtual bool preInitialize() const;
/**
* Called from doinit() to extract PDFs from the event file and add
* the corresponding objects to the current EventGenerator.
*/
virtual void initPDFs();
//@}
protected:
/**
* The HEPRUP common block.
*/
HEPRUP heprup;
/**
* The HEPEUP common block.
*/
HEPEUP hepeup;
/**
* The ParticleData objects corresponding to the incoming particles.
*/
tcPDPair inData;
/**
* The PDFBase objects which has been used for the beam particle
* when generating the events being read. Specified in the interface
* or derived from PDFGUP and PDFSUP.
*/
pair<PDFPtr,PDFPtr> inPDF;
/**
* The PDFBase object to be used in the subsequent generation.
*/
pair<cPDFPtr,cPDFPtr> outPDF;
/**
* The PartonExtractor object used to construct remnants.
*/
PExtrPtr thePartonExtractor;
/**
* A pointer to a CascadeHandler to be used for CKKW-reweighting.
*/
tCascHdlPtr theCKKW;
/**
* The pairs of PartonBin objects describing the partons which can
* be extracted by the PartonExtractor object.
*/
PartonPairVec thePartonBins;
/**
* The map of XComb objects indexed by the corresponding PartonBin
* pair.
*/
XCombMap theXCombs;
/**
* The Cuts object to be used for this reader.
*/
CutsPtr theCuts;
/**
* The number of events in this reader. If less than zero the number
* of events is unlimited.
*/
long theNEvents;
/**
* The number of events produced by this reader so far. Is reset
* every time an event file is reopened.
*/
long position;
/**
* The number of times this reader has been reopened.
*/
int reopened;
/**
* The maximum number of events to scan to collect information about
* processes and cross sections. If less than 0, all events will be
* scanned.
*/
long theMaxScan;
/**
* Flag to tell whether we are in the process of scanning.
*/
bool scanning;
/**
* True if this is an active reader.
*/
bool isActive;
/**
* Name of file used to cache the events form the reader in a
* fast-readable form. If empty, no cache file will be generated.
*/
string theCacheFileName;
/**
* Determines whether to apply cuts to events before converting them
* to ThePEG format.
*/
bool doCutEarly;
/**
* Collect statistics for this reader.
*/
XSecStat stats;
/**
* Collect statistics for each individual process.
*/
StatMap statmap;
/**
* The pair of PartonBinInstance objects describing the current
* incoming partons in the event.
*/
PBIPair thePartonBinInstances;
/**
* Association between ColourLines and colour indices in the current
* translation.
*/
ObjectIndexer<long,ColourLine> colourIndex;
/**
* Association between Particles and indices in the current
* translation.
*/
ObjectIndexer<long,Particle> particleIndex;
/**
* The instances of the beam particles for the current event.
*/
PPair theBeams;
/**
* The instances of the incoming particles to the sub process for
* the current event.
*/
PPair theIncoming;
/**
* The instances of the outgoing particles from the sub process for
* the current event.
*/
PVector theOutgoing;
/**
* The instances of the intermediate particles in the sub process for
* the current event.
*/
PVector theIntermediates;
/**
* File stream for the cache.
*/
CFile theCacheFile;
/**
* The reweight objects modifying the weights of this reader.
*/
ReweightVector reweights;
/**
* The preweight objects modifying the weights of this reader.
*/
ReweightVector preweights;
/**
* The factor with which this reader was last pre-weighted.
*/
double preweight;
/**
* Should the event be reweighted by PDFs used by the PartonExtractor?
*/
bool reweightPDF;
/**
* Should PDFBase objects be constructed from the information in the
* event file in the initialization?
*/
bool doInitPDFs;
/**
* If this reader is to be used (possibly together with others) for
* CKKW reweighting and veto, this should give the multiplicity of
* outgoing particles in the highest multiplicity matrix element in
* the group.
*/
int theMaxMultCKKW;
/**
* If this reader is to be used (possibly together with others) for
* CKKW reweighting and veto, this should give the multiplicity of
* outgoing particles in the lowest multiplicity matrix element in
* the group.
*/
int theMinMultCKKW;
/**
* The weight multiplying the last read event due to PDF
* reweighting, CKKW reweighting or assigned reweight and preweight
* objects.
*/
double lastweight;
/**
* The optional weights associated to the last read events.
*/
map<string,double> optionalWeights;
/**
* The event number
*/
long LHEeventnum;
/**
* If the maximum cross section of this reader has been increased
* with increaseMaxXSec(), this is the total factor with which it
* has been increased.
*/
double maxFactor;
/**
* npLO for FxFx merging
*/
int optionalnpLO;
/**
* npNLO for FxFx merging
*/
int optionalnpNLO;
/**
* The (reweighted) XWGTUP value should be scaled with this cross
* section when compared to the overestimated cross section.
*/
CrossSection weightScale;
/**
* Individual scales for different sub-processes if reweighted.
*/
vector<double> xSecWeights;
/**
* Individual maximum weights for individual (possibly reweighted)
* processes.
*/
map<int,double> maxWeights;
/**
* Is set to true when getEvent() is called from skip(int).
*/
bool skipping;
/**
* Option for the treatment of the momenta supplied
*/
unsigned int theMomentumTreatment;
/**
* Set to true if warnings about possible weight incompatibilities
* should be issued.
*/
bool useWeightWarnings;
/**
* Option to allow reopening of the file
*/
bool theReOpenAllowed;
/**
* Use the spin information
*/
bool theIncludeSpin;
private:
/** Access function for the interface. */
void setBeamA(long id);
/** Access function for the interface. */
long getBeamA() const;
/** Access function for the interface. */
void setBeamB(long id);
/** Access function for the interface. */
long getBeamB() const;
/** Access function for the interface. */
void setEBeamA(Energy e);
/** Access function for the interface. */
Energy getEBeamA() const;
/** Access function for the interface. */
void setEBeamB(Energy e);
/** Access function for the interface. */
Energy getEBeamB() const;
/** Access function for the interface. */
void setPDFA(PDFPtr);
/** Access function for the interface. */
PDFPtr getPDFA() const;
/** Access function for the interface. */
void setPDFB(PDFPtr);
/** Access function for the interface. */
PDFPtr getPDFB() const;
private:
/**
* Describe an abstract base class with persistent data.
*/
static AbstractClassDescription<FxFxReader> initFxFxReader;
/**
* Private and non-existent assignment operator.
*/
FxFxReader & operator=(const FxFxReader &) = delete;
public:
/** @cond EXCEPTIONCLASSES */
/** Exception class used by FxFxReader in case inconsistencies
* are encountered. */
class FxFxInconsistencyError: public Exception {};
/** Exception class used by FxFxReader in case more events
than available are requested. */
class FxFxReopenWarning: public Exception {};
/** Exception class used by FxFxReader in case reopening an
event file fails. */
class FxFxReopenError: public Exception {};
/** Exception class used by FxFxReader in case there is
information missing in the initialization phase. */
class FxFxInitError: public InitException {};
/** @endcond */
};
/// Stream output for HEPEUP
ostream & operator<<(ostream & os, const HEPEUP & h);
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the
* base class of FxFxReader.
*/
template <>
struct BaseClassTrait<FxFxReader,1>: public ClassTraitsType {
/** Typedef of the base class of FxFxReader. */
typedef HandlerBase NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* FxFxReader class and the shared object where it is
* defined.
*/
template <>
struct ClassTraits<FxFxReader>
: public ClassTraitsBase<FxFxReader> {
/**
* Return the class name.
*/
static string className() { return "Herwig::FxFxReader"; }
/**
* Return the name of the shared library to be loaded to get access
* to the FxFxReader class and every other class it uses
* (except the base class).
*/
static string library() { return "HwFxFx.so"; }
};
/** @endcond */
}
#endif /* THEPEG_FxFxReader_H */
diff --git a/Shower/QTilde/Dark/HalfHalfOneDarkSplitFn.cc b/Shower/QTilde/Dark/HalfHalfOneDarkSplitFn.cc
--- a/Shower/QTilde/Dark/HalfHalfOneDarkSplitFn.cc
+++ b/Shower/QTilde/Dark/HalfHalfOneDarkSplitFn.cc
@@ -1,171 +1,171 @@
// -*- C++ -*-
//
// HalfHalfOneDarkSplitFn.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 HalfHalfOneDarkSplitFn class.
//
#include "HalfHalfOneDarkSplitFn.h"
#include "HiddenValleyModel.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <cassert>
using namespace Herwig;
DescribeNoPIOClass<HalfHalfOneDarkSplitFn,Herwig::Sudakov1to2FormFactor>
describeHalfHalfOneDarkSplitFn ("Herwig::HalfHalfOneDarkSplitFn","HwShower.so");
void HalfHalfOneDarkSplitFn::Init() {
// documentation
static ClassDocumentation<HalfHalfOneDarkSplitFn> documentation
("The HalfHalfOneDarkSplitFn class implements the dark coloured q -> qg splitting function");
}
bool HalfHalfOneDarkSplitFn::checkColours(const IdList & ids) const {
if(colourStructure()==TripletTripletOctet) {
if(ids[0]!=ids[1]) return false;
if((ids[0]->iColour()==PDT::DarkColourFundamental||
ids[0]->iColour()==PDT::DarkColourAntiFundamental) &&
ids[2]->iColour()==PDT::DarkColourAdjoint) return true;
return false;
}
else if(colourStructure()==OctetOctetOctet) {
for(unsigned int ix=0;ix<3;++ix) {
if(ids[ix]->iColour()!=PDT::DarkColourAdjoint) return false;
}
return true;
}
else if(colourStructure()==OctetTripletTriplet) {
if(ids[0]->iColour()!=PDT::DarkColourAdjoint) return false;
if(ids[1]->iColour()==PDT::DarkColourFundamental&&
ids[2]->iColour()==PDT::DarkColourAntiFundamental)
return true;
if(ids[1]->iColour()==PDT::DarkColourAntiFundamental&&
ids[2]->iColour()==PDT::DarkColourFundamental)
return true;
return false;
}
else if(colourStructure()==TripletOctetTriplet) {
if(ids[0]!=ids[2]) return false;
if((ids[0]->iColour()==PDT::DarkColourFundamental||
ids[0]->iColour()==PDT::DarkColourAntiFundamental) &&
ids[1]->iColour()==PDT::DarkColourAdjoint) return true;
return false;
}
else {
assert(false);
}
return false;
}
double HalfHalfOneDarkSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & ) const {
double val = (1. + sqr(z))/(1.-z);
if(mass) {
Energy m = ids[0]->mass();
val -= 2.*sqr(m)/t;
}
return val;
}
double HalfHalfOneDarkSplitFn::overestimateP(const double z,
- const IdList & ids) const {
+ const IdList & ) const {
return 2./(1.-z);
}
double HalfHalfOneDarkSplitFn::ratioP(const double z, const Energy2 t,
- const IdList & ids, const bool mass, const RhoDMatrix & ) const {
+ const IdList & ids, const bool mass, const RhoDMatrix & ) const {
double val = 1. + sqr(z);
if(mass) {
Energy m = ids[0]->mass();
val -= 2.*sqr(m)*(1.-z)/t;
}
return 0.5*val;
}
double HalfHalfOneDarkSplitFn::integOverP(const double z,
- const IdList & ids,
- unsigned int PDFfactor) const {
+ const IdList & ,
+ unsigned int PDFfactor) const {
switch (PDFfactor) {
case 0:
return -2.*Math::log1m(z);
case 1:
return 2.*log(z/(1.-z));
case 2:
return 2./(1.-z);
case 3:
default:
throw Exception() << "HalfHalfOneDarkSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
-double HalfHalfOneDarkSplitFn::invIntegOverP(const double r, const IdList & ids,
- unsigned int PDFfactor) const {
+double HalfHalfOneDarkSplitFn::invIntegOverP(const double r, const IdList &,
+ unsigned int PDFfactor) const {
switch (PDFfactor) {
case 0:
return 1. - exp(- 0.5*r);
case 1:
return 1./(1.-exp(-0.5*r));
case 2:
return 1.-2./r;
case 3:
default:
throw Exception() << "HalfHalfOneDarkSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool HalfHalfOneDarkSplitFn::accept(const IdList &ids) const {
// 3 particles and in and out fermion same
if(ids.size()!=3 || ids[0]!=ids[1]) return false;
if(ids[0]->iSpin()!=PDT::Spin1Half ||
ids[2]->iSpin()!=PDT::Spin1) return false;
return checkColours(ids);
}
vector<pair<int, Complex> >
HalfHalfOneDarkSplitFn::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 {{ {0, 1.} }};
}
vector<pair<int, Complex> >
HalfHalfOneDarkSplitFn::generatePhiBackward(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 {{ {0, 1.} }};
}
DecayMEPtr HalfHalfOneDarkSplitFn::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::Spin1Half,PDT::Spin1)));
Energy m = !timeLike ? ZERO : ids[0]->mass();
double mt = m/sqrt(t);
double root = sqrt(1.-(1.-z)*sqr(m)/z/t);
double romz = sqrt(1.-z);
double rz = sqrt(z);
Complex phase = exp(Complex(0.,1.)*phi);
(*kernal)(0,0,0) = -root/romz*phase;
(*kernal)(1,1,2) = -conj((*kernal)(0,0,0));
(*kernal)(0,0,2) = root/romz*z/phase;
(*kernal)(1,1,0) = -conj((*kernal)(0,0,2));
(*kernal)(1,0,2) = mt*(1.-z)/rz;
(*kernal)(0,1,0) = conj((*kernal)(1,0,2));
(*kernal)(0,1,2) = 0.;
(*kernal)(1,0,0) = 0.;
return kernal;
}
diff --git a/Shower/QTilde/Dark/OneHalfHalfDarkSplitFn.cc b/Shower/QTilde/Dark/OneHalfHalfDarkSplitFn.cc
--- a/Shower/QTilde/Dark/OneHalfHalfDarkSplitFn.cc
+++ b/Shower/QTilde/Dark/OneHalfHalfDarkSplitFn.cc
@@ -1,187 +1,187 @@
// -*- C++ -*-
//
// OneHalfHalfDarkSplitFn.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 OneHalfHalfDarkSplitFn class.
//
#include "OneHalfHalfDarkSplitFn.h"
#include "HiddenValleyModel.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeNoPIOClass<OneHalfHalfDarkSplitFn,Herwig::Sudakov1to2FormFactor>
describeOneHalfHalfDarkSplitFn ("Herwig::OneHalfHalfDarkSplitFn","HwShower.so");
void OneHalfHalfDarkSplitFn::Init() {
static ClassDocumentation<OneHalfHalfDarkSplitFn> documentation
("The OneHalfHalfDarkSplitFn class implements the splitting function for g->q qbar");
}
bool OneHalfHalfDarkSplitFn::checkColours(const IdList & ids) const {
if(colourStructure()==TripletTripletOctet) {
if(ids[0]!=ids[1]) return false;
if((ids[0]->iColour()==PDT::DarkColourFundamental||
ids[0]->iColour()==PDT::DarkColourAntiFundamental) &&
ids[2]->iColour()==PDT::DarkColourAdjoint) return true;
return false;
}
else if(colourStructure()==OctetOctetOctet) {
for(unsigned int ix=0;ix<3;++ix) {
if(ids[ix]->iColour()!=PDT::DarkColourAdjoint) return false;
}
return true;
}
else if(colourStructure()==OctetTripletTriplet) {
if(ids[0]->iColour()!=PDT::DarkColourAdjoint) return false;
if(ids[1]->iColour()==PDT::DarkColourFundamental&&
ids[2]->iColour()==PDT::DarkColourAntiFundamental)
return true;
if(ids[1]->iColour()==PDT::DarkColourAntiFundamental&&
ids[2]->iColour()==PDT::DarkColourFundamental)
return true;
return false;
}
else if(colourStructure()==TripletOctetTriplet) {
if(ids[0]!=ids[2]) return false;
if((ids[0]->iColour()==PDT::DarkColourFundamental||
ids[0]->iColour()==PDT::DarkColourAntiFundamental) &&
ids[1]->iColour()==PDT::DarkColourAdjoint) return true;
return false;
}
else {
assert(false);
}
return false;
}
double OneHalfHalfDarkSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix &) const {
double zz = z*(1.-z);
double val=1.-2.*zz;
if(mass) {
Energy m = ids[1]->mass();
val +=2.*sqr(m)/t;
}
return val;
}
double OneHalfHalfDarkSplitFn::overestimateP(const double,
- const IdList &ids) const {
+ const IdList &) const {
return 1.;
}
double OneHalfHalfDarkSplitFn::ratioP(const double z, const Energy2 t,
- const IdList &ids, const bool mass, const RhoDMatrix &) const {
+ const IdList &ids, const bool mass, const RhoDMatrix &) const {
double zz = z*(1.-z);
double val = 1.-2.*zz;
if(mass) {
Energy m = ids[1]->mass();
val+= 2.*sqr(m)/t;
}
return val;
}
-double OneHalfHalfDarkSplitFn::integOverP(const double z, const IdList & ids,
- unsigned int PDFfactor) const {
+double OneHalfHalfDarkSplitFn::integOverP(const double z, const IdList & ,
+ unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return z;
case 1:
return log(z);
case 2:
return -log(1.-z);
case 3:
return log(z/(1.-z));
case 4:
return 2.*sqrt(z);
case 5:
return (2./3.)*z*sqrt(z);
default:
throw Exception() << "OneHalfHalfDarkSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double OneHalfHalfDarkSplitFn::invIntegOverP(const double r,
- const IdList & ids,
- unsigned int PDFfactor) const {
+ const IdList & ,
+ unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return r;
case 1:
return exp(r);
case 2:
return 1.-exp(-r);
case 3:
return 1./(1.+exp(-r));
case 4:
return 0.25*sqr(r);
case 5:
return pow(1.5*r,2./3.);
default:
throw Exception() << "OneHalfHalfDarkSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool OneHalfHalfDarkSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3) return false;
if(ids[1]!=ids[2]->CC()) return false;
if(ids[1]->iSpin()!=PDT::Spin1Half) return false;
if(ids[0]->iSpin()!=PDT::Spin1) return false;
return OneHalfHalfDarkSplitFn::checkColours(ids);
}
vector<pair<int, Complex> >
OneHalfHalfDarkSplitFn::generatePhiForward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix & rho) {
assert(rho.iSpin()==PDT::Spin1);
double modRho = abs(rho(0,2));
Energy mq = ids[1]->mass();
Energy2 mq2 = sqr(mq);
double fact = z*(1.-z)-mq2/t;
double max = 1.+2.*fact*(-1.+2.*modRho);
vector<pair<int, Complex> > output;
output.push_back(make_pair( 0,(rho(0,0)+rho(2,2))*(1.-2.*fact)/max));
output.push_back(make_pair(-2,2.*fact*rho(0,2)/max));
output.push_back(make_pair( 2,2.*fact*rho(2,0)/max));
return output;
}
vector<pair<int, Complex> >
OneHalfHalfDarkSplitFn::generatePhiBackward(const double, const Energy2, const IdList &,
const RhoDMatrix & ) {
// no dependance
return {{ {0, 1.} }};
}
DecayMEPtr OneHalfHalfDarkSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool timeLike) {
static const Complex ii(0.,1.);
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half)));
double mt = !timeLike ? ZERO : ids[1]->mass()/sqrt(t);
double root =sqrt(1.-sqr(mt)/z/(1.-z));
(*kernal)(0,0,0) = mt/sqrt(z*(1.-z));
(*kernal)(2,1,1) = (*kernal)(0,0,0);
(*kernal)(0,0,1) = -z*root*exp(-ii*phi);
(*kernal)(2,1,0) = -conj((*kernal)(0,0,1));
(*kernal)(0,1,0) = (1.-z)*exp(-ii*phi)*root;
(*kernal)(2,0,1) = -conj((*kernal)(0,1,0));
(*kernal)(0,1,1) = 0.;
(*kernal)(2,0,0) = 0.;
return kernal;
}
diff --git a/Shower/QTilde/Dark/OneOneOneDarkSplitFn.cc b/Shower/QTilde/Dark/OneOneOneDarkSplitFn.cc
--- a/Shower/QTilde/Dark/OneOneOneDarkSplitFn.cc
+++ b/Shower/QTilde/Dark/OneOneOneDarkSplitFn.cc
@@ -1,167 +1,167 @@
// -*- C++ -*-
//
// OneOneOneDarkSplitFn.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 OneOneOneDarkSplitFn class.
//
#include "OneOneOneDarkSplitFn.h"
#include "HiddenValleyModel.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeNoPIOClass<OneOneOneDarkSplitFn,Herwig::Sudakov1to2FormFactor>
describeOneOneOneDarkSplitFn ("Herwig::OneOneOneDarkSplitFn","HwShower.so");
void OneOneOneDarkSplitFn::Init() {
static ClassDocumentation<OneOneOneDarkSplitFn> documentation
("The OneOneOneDarkSplitFn class implements the g -> gg splitting function");
}
bool OneOneOneDarkSplitFn::checkColours(const IdList & ids) const {
if(colourStructure()==TripletTripletOctet) {
if(ids[0]!=ids[1]) return false;
if((ids[0]->iColour()==PDT::DarkColourFundamental||
ids[0]->iColour()==PDT::DarkColourAntiFundamental) &&
ids[2]->iColour()==PDT::DarkColourAdjoint) return true;
return false;
}
else if(colourStructure()==OctetOctetOctet) {
for(unsigned int ix=0;ix<3;++ix) {
if(ids[ix]->iColour()!=PDT::DarkColourAdjoint) return false;
}
return true;
}
else if(colourStructure()==OctetTripletTriplet) {
if(ids[0]->iColour()!=PDT::DarkColourAdjoint) return false;
if(ids[1]->iColour()==PDT::DarkColourFundamental&&
ids[2]->iColour()==PDT::DarkColourAntiFundamental)
return true;
if(ids[1]->iColour()==PDT::DarkColourAntiFundamental&&
ids[2]->iColour()==PDT::DarkColourFundamental)
return true;
return false;
}
else if(colourStructure()==TripletOctetTriplet) {
if(ids[0]!=ids[2]) return false;
if((ids[0]->iColour()==PDT::DarkColourFundamental||
ids[0]->iColour()==PDT::DarkColourAntiFundamental) &&
ids[1]->iColour()==PDT::DarkColourAdjoint) return true;
return false;
}
else {
assert(false);
}
return false;
}
double OneOneOneDarkSplitFn::P(const double z, const Energy2,
- const IdList & ids, const bool, const RhoDMatrix &)const {
+ const IdList & , const bool, const RhoDMatrix &)const {
return sqr(1.-z*(1.-z))/(z*(1.-z));
}
double OneOneOneDarkSplitFn::overestimateP(const double z,
- const IdList & ids) const {
+ const IdList & ) const {
return (1/z + 1/(1.-z));
}
double OneOneOneDarkSplitFn::ratioP(const double z, const Energy2,
- const IdList & , const bool, const RhoDMatrix &) const {
+ const IdList & , const bool, const RhoDMatrix &) const {
return sqr(1.-z*(1.-z));
}
double OneOneOneDarkSplitFn::invIntegOverP(const double r,
- const IdList & ids,
- unsigned int PDFfactor) const {
+ const IdList & ,
+ unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return 1./(1.+exp(-r));
case 1:
case 2:
case 3:
default:
throw Exception() << "OneOneOneDarkSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
-double OneOneOneDarkSplitFn::integOverP(const double z, const IdList & ids,
- unsigned int PDFfactor) const {
+double OneOneOneDarkSplitFn::integOverP(const double z, const IdList & ,
+ unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
assert(z>0.&&z<1.);
return log(z/(1.-z));
case 1:
case 2:
case 3:
default:
throw Exception() << "OneOneOneDarkSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool OneOneOneDarkSplitFn::accept(const IdList & ids) const {
if(ids.size()!=3) return false;
for(unsigned int ix=0;ix<ids.size();++ix) {
if(ids[0]->iSpin()!=PDT::Spin1) return false;
}
return OneOneOneDarkSplitFn::checkColours(ids);
}
vector<pair<int, Complex> >
OneOneOneDarkSplitFn::generatePhiForward(const double z, const Energy2, const IdList &,
const RhoDMatrix & rho) {
assert(rho.iSpin()==PDT::Spin1);
double modRho = abs(rho(0,2));
double max = 2.*z*modRho*(1.-z)+sqr(1.-(1.-z)*z)/(z*(1.-z));
vector<pair<int, Complex> > output;
output.push_back(make_pair( 0,(rho(0,0)+rho(2,2))*sqr(1.-(1.-z)*z)/(z*(1.-z))/max));
output.push_back(make_pair(-2,-rho(0,2)*z*(1.-z)/max));
output.push_back(make_pair( 2,-rho(2,0)*z*(1.-z)/max));
return output;
}
vector<pair<int, Complex> >
OneOneOneDarkSplitFn::generatePhiBackward(const double z, const Energy2, const IdList &,
const RhoDMatrix & rho) {
assert(rho.iSpin()==PDT::Spin1);
double diag = sqr(1 - (1 - z)*z)/(1 - z)/z;
double off = (1.-z)/z;
double max = 2.*abs(rho(0,2))*off+diag;
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));
return output;
}
DecayMEPtr OneOneOneDarkSplitFn::matrixElement(const double z, const Energy2,
const IdList &, const double phi,
bool) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin1)));
double omz = 1.-z;
double root = sqrt(z*omz);
Complex phase = exp(Complex(0.,1.)*phi);
(*kernal)(0,0,0) = phase/root;
(*kernal)(2,2,2) = -conj((*kernal)(0,0,0));
(*kernal)(0,0,2) = -sqr(z)/root/phase;
(*kernal)(2,2,0) = -conj((*kernal)(0,0,2));
(*kernal)(0,2,0) = -sqr(omz)/root/phase;
(*kernal)(2,0,2) = -conj((*kernal)(0,2,0));
(*kernal)(0,2,2) = 0.;
(*kernal)(2,0,0) = 0.;
return kernal;
}
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,117 +1,117 @@
// -*- C++ -*-
//
// HalfOneHalfSplitFn.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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::Sudakov1to2FormFactor>
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::integOverP(const double z, const IdList & ids,
+double HalfOneHalfSplitFn::integOverP(const double z, const IdList &,
unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return 2.*log(z);
case 1:
return -2./z;
case 2:
return 2.*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,
+ const IdList & ,
unsigned int PDFfactor) const {
switch(PDFfactor) {
case 0:
return exp(0.5*r);
case 1:
return -2./r;
case 2:
return 1./(1.+exp(-0.5*r));
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 {{ {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));
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);
(*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/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h b/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h
--- a/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h
+++ b/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h
@@ -1,180 +1,180 @@
// -*- C++ -*-
//
// HalfOneHalfSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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_HalfOneHalfSplitFn_H
#define HERWIG_HalfOneHalfSplitFn_H
//
// This is the declaration of the HalfOneHalfSplitFn class.
//
#include "Sudakov1to2FormFactor.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
*
* This classs provides the concrete implementation of the exact leading-order
* splitting function for \f$\frac12\to 1\frac12\f$.
*
* In this case the splitting function is given by
* \f[P(z,t) = C\left(\frac{2(1-z)+z^2}{z}-2\frac{m^2_q}t\right),\f]
* where \f$C\f$ is the corresponding colour factor.
* Our choice for the overestimate is
* \f[P_{\rm over}(z) = 2C\frac1z,\f]
* therefore the integral is
* \f[\int P_{\rm over}(z) {\rm d}z = 2C\ln z,\f]
* and its inverse is
* \f[\exp\left(\frac{r}{2C}\right).\f]
*
* @see SplittingFunction
*/
class HalfOneHalfSplitFn: public Sudakov1to2FormFactor {
public:
/**
* Concrete implementation of the method to determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual bool accept(const IdList & ids) const;
/**
* Methods to return the splitting function.
*/
//@{
/**
* The concrete implementation of the overestimate of the splitting function,
* \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
- virtual double overestimateP(const double z, const IdList & ids) const {
+ virtual double overestimateP(const double z, const IdList & ) const {
return 2./z;
}
/**
* The concrete implementation of the
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,t)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
* @param rho The spin density matrix
*/
virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
- const bool mass, const RhoDMatrix & rho) const {
+ 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;
}
/**
* The concrete implementation of the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double integOverP(const double z, const IdList & ids,
unsigned int PDFfactor=0) const;
/**
* The concrete implementation of the inverse of the indefinite integral.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor=0) const;
//@}
/**
* Method to calculate the azimuthal angle for forward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiForward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix &);
/**
* Method to calculate the azimuthal angle for backward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiBackward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix &);
/**
* Calculate the matrix element for the splitting
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
*/
virtual DecayMEPtr matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi, bool timeLike);
public:
/**
* 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 {return new_ptr(*this);}
/** 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 {return new_ptr(*this);}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HalfOneHalfSplitFn & operator=(const HalfOneHalfSplitFn &) = delete;
};
}
#endif /* HERWIG_HalfOneHalfSplitFn_H */
diff --git a/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h b/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h
--- a/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h
+++ b/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h
@@ -1,195 +1,195 @@
// -*- C++ -*-
//
// OneHalfHalfSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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_OneHalfHalfSplitFn_H
#define HERWIG_OneHalfHalfSplitFn_H
//
// This is the declaration of the OneHalfHalfSplitFn class.
//
#include "Sudakov1to2FormFactor.h"
namespace Herwig {
using namespace ThePEG;
/**\ingroup Shower
*
* This class provides the concrete implementation of the exact leading-order
* splitting function for \f$1\to \frac12\frac12\f$.
*
* In this case the splitting function is given by
* \f[P(z,t) =C\left(1-2z(1-z)+2\frac{m_q^2}{t}\right),\f]
* where \f$C\f$ is the corresponding colour factor
* Our choice for the overestimate is
* \f[P_{\rm over}(z) = C,\f]
* therefore the integral is
* \f[\int P_{\rm over}(z) {\rm d}z =Cz,\f]
* and its inverse is
* \f[\frac{r}{C}\f]
*
* @see \ref OneHalfHalfSplitFnInterfaces "The interfaces"
* defined for OneHalfHalfSplitFn.
*/
class OneHalfHalfSplitFn: public Sudakov1to2FormFactor {
public:
/**
* Concrete implementation of the method to determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual bool accept(const IdList & ids) const {
if(ids.size()!=3) return false;
if(ids[1]!=ids[2]->CC()) return false;
if(ids[1]->iSpin()!=PDT::Spin1Half) return false;
if(ids[0]->iSpin()!=PDT::Spin1) return false;
return true;
}
/**
* Methods to return the splitting function.
*/
//@{
/**
* The concrete implementation of the overestimate of the splitting function,
* \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
- virtual double overestimateP(const double z, const IdList & ids) const {
+ virtual double overestimateP(const double , const IdList & ) const {
return 1.;
}
/**
* The concrete implementation of the
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
* @param rho The spin density matrix
*/
virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
- const bool mass, const RhoDMatrix & rho) const {
+ const bool mass, const RhoDMatrix &) const {
double zz = z*(1.-z);
double val = 1.-2.*zz;
if(mass) {
Energy m = ids[1]->mass();
val+= 2.*sqr(m)/t;
}
return val;
}
/**
* The concrete implementation of the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double integOverP(const double z, const IdList & ids,
unsigned int PDFfactor=0) const;
/**
* The concrete implementation of the inverse of the indefinite integral.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor=0) const;
//@}
/**
* Method to calculate the azimuthal angle
* @param particle The particle which is branching
* @param showerkin The ShowerKinematics object
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiForward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix &);
/**
* Method to calculate the azimuthal angle
* @param particle The particle which is branching
* @param showerkin The ShowerKinematics object
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiBackward(const double, const Energy2, const IdList &,
const RhoDMatrix &) {
// no dependance
return {{ {0, 1.} }};
}
/**
* Calculate the matrix element for the splitting
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
*/
virtual DecayMEPtr matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi, bool timeLike);
public:
/**
* 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 {return new_ptr(*this);}
/** 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 {return new_ptr(*this);}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
OneHalfHalfSplitFn & operator=(const OneHalfHalfSplitFn &) = delete;
};
}
#endif /* HERWIG_OneHalfHalfSplitFn_H */
diff --git a/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h b/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h
--- a/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h
+++ b/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h
@@ -1,225 +1,225 @@
// -*- C++ -*-
//
// OneOneOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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_OneOneOneSplitFn_H
#define HERWIG_OneOneOneSplitFn_H
//
// This is the declaration of the OneOneOneSplitFn class.
//
#include "Sudakov1to2FormFactor.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
*
* This class provides the concrete implementation of the exact leading-order
* splitting function for \f$1\to 11\f$.
*
* In this case the splitting function is given by
* \f[P(z) =2C*\left(\frac{z}{1-z}+\frac{1-z}{z}+z(1-z)\right),\f]
* where \f$C=\f$ is the corresponding colour factor.
* Our choice for the overestimate is
* \f[P_{\rm over}(z) = 2C\left(\frac1z+\frac1{1-z}\right),\f]
* therefore the integral is
* \f[\int P_{\rm over}(z) {\rm d}z =2C\ln\left(\frac{z}{1-z}\right),\f]
* and its inverse is
* \f[\frac{\exp\left(\frac{r}{2C}\right)}{(1+\exp\left(\frac{r}{2C}\right)}\f]
*
*
* @see \ref OneOneOneSplitFnInterfaces "The interfaces"
* defined for OneOneOneSplitFn.
*/
class OneOneOneSplitFn: public Sudakov1to2FormFactor {
public:
/**
* Concrete implementation of the method to determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
bool accept(const IdList & ids) const {
if(ids.size()!=3) return false;
for(unsigned int ix=0;ix<ids.size();++ix) {
if(ids[ix]->iSpin()!=PDT::Spin1) return false;
}
return true;
}
/**
* Methods to return the splitting function.
*/
//@{
/**
* The concrete implementation of the overestimate of the splitting function,
* \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
- double overestimateP(const double z, const IdList & ids) const {
+ double overestimateP(const double z, const IdList &) const {
return 1/z + 1/(1.-z);
}
/**
* The concrete implementation of the
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,t)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
* @param rho The spin density matrix
*/
virtual double ratioP(const double z, const Energy2, const IdList &,
const bool, const RhoDMatrix &) const {
return sqr(1.-z*(1.-z));
}
/**
* The concrete implementation of the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double integOverP(const double z, const IdList &,
unsigned int PDFfactor=0) const {
assert(PDFfactor==0);
assert(z>0.&&z<1.);
return log(z/(1.-z));
}
/**
* The concrete implementation of the inverse of the indefinite integral.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList &,
unsigned int PDFfactor=0) const {
assert(PDFfactor==0);
return 1./(1.+exp(-r));
}
//@}
/**
* Method to calculate the azimuthal angle for forward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
vector<pair<int,Complex> >
generatePhiForward(const double z, const Energy2, const IdList &,
const RhoDMatrix & rho) {
assert(rho.iSpin()==PDT::Spin1);
double modRho = abs(rho(0,2));
double max = 2.*z*modRho*(1.-z)+sqr(1.-(1.-z)*z)/(z*(1.-z));
vector<pair<int, Complex> > output;
output.push_back(make_pair( 0,(rho(0,0)+rho(2,2))*sqr(1.-(1.-z)*z)/(z*(1.-z))/max));
output.push_back(make_pair(-2,-rho(0,2)*z*(1.-z)/max));
output.push_back(make_pair( 2,-rho(2,0)*z*(1.-z)/max));
return output;
}
/**
* Method to calculate the azimuthal angle for backward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
vector<pair<int,Complex> >
generatePhiBackward(const double z, const Energy2, const IdList &,
const RhoDMatrix & rho) {
assert(rho.iSpin()==PDT::Spin1);
double diag = sqr(1 - (1 - z)*z)/(1 - z)/z;
double off = (1.-z)/z;
double max = 2.*abs(rho(0,2))*off+diag;
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));
return output;
}
/**
* Calculate the matrix element for the splitting
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
*/
virtual DecayMEPtr matrixElement(const double z, const Energy2,
const IdList &, const double phi, bool) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin1)));
double omz = 1.-z;
double root = sqrt(z*omz);
Complex phase = exp(Complex(0.,1.)*phi);
(*kernal)(0,0,0) = phase/root;
(*kernal)(2,2,2) = -conj((*kernal)(0,0,0));
(*kernal)(0,0,2) = -sqr(z)/root/phase;
(*kernal)(2,2,0) = -conj((*kernal)(0,0,2));
(*kernal)(0,2,0) = -sqr(omz)/root/phase;
(*kernal)(2,0,2) = -conj((*kernal)(0,2,0));
(*kernal)(0,2,2) = 0.;
(*kernal)(2,0,0) = 0.;
return kernal;
}
public:
/**
* 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 {return new_ptr(*this);}
/** 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 {return new_ptr(*this);}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
OneOneOneSplitFn & operator=(const OneOneOneSplitFn &) = delete;
};
}
#endif /* HERWIG_OneOneOneSplitFn_H */
diff --git a/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.h b/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.h
--- a/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.h
+++ b/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.h
@@ -1,256 +1,256 @@
// -*- C++ -*-
#ifndef Herwig_ZeroZeroOneEWSplitFn_H
#define Herwig_ZeroZeroOneEWSplitFn_H
//
// This is the declaration of the ZeroZeroOneEWSplitFn class.
//
#include "Sudakov1to2FormFactor.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
/**
* The ZeroZeroOneEWSplitFn class implements the splitting function for
* \f$\1\to q\1 0\f$ where the spin-1 particles are the W / Z massive
* electroweak gauge bosons and the spin-0 particle is the massive Higgs
* boson.
*
* @see \ref ZeroZeroOneEWSplitFnInterfaces "The interfaces"
* defined for ZeroZeroOneEWSplitFn.
*/
class ZeroZeroOneEWSplitFn: public Sudakov1to2FormFactor {
public:
/**
* Concrete implementation of the method to determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual bool accept(const IdList & ids) const {
if(ids.size()!=3)
return false;
if(_couplingValueIm==0.&&_couplingValueRe==0.)
return false;
if(ids[0]->iCharge()!=ids[1]->iCharge()+ids[2]->iCharge())
return false;
if(ids[0]->iSpin()==PDT::Spin0 && ids[1]->iSpin()==PDT::Spin1 && ids[2]->iSpin()==PDT::Spin0)
return true;
else if(ids[0]->iSpin()==PDT::Spin0 && ids[1]->iSpin()==PDT::Spin0 && ids[2]->iSpin()==PDT::Spin1)
return true;
else
return false;
}
/**
* Methods to return the splitting function.
*/
//@{
/**
* The concrete implementation of the overestimate of the splitting function,
* \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual double overestimateP(const double z, const IdList & ids) const {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
return norm(ghhv)*2*z/(1.-z);
}
/**
* The concrete implementation of the
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,t)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
* @param rho The spin density matrix
*/
virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
- const bool mass, const RhoDMatrix & rho) const {
+ const bool mass, const RhoDMatrix & ) const {
// ratio in the massless limit
double val = 1.;
// the massive limit
if(mass){
// get the running mass
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
val += m0t2-(1./z)*m1t2+((1.-z)/(4.*z))*m2t2;
}
return val;
}
/**
* The concrete implementation of the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double integOverP(const double z, const IdList & ids,
unsigned int PDFfactor=0) const {
assert(PDFfactor==0);
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
double pre = norm(ghhv)*2.;
return -pre*(z+log(1.-z));
}
/**
* The concrete implementation of the inverse of the indefinite integral.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor=0) const {
assert(PDFfactor==0);
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
double pre = norm(ghhv)*2.;
return 1.-exp(-(1.+r/pre));
}
//@}
/**
* Method to calculate the azimuthal angle
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
- generatePhiForward(const double z, const Energy2 t, const IdList & ids,
- const RhoDMatrix &) {
+ 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.));
}
/**
* Method to calculate the azimuthal angle for backward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
- generatePhiBackward(const double z, const Energy2 t, const IdList & ids,
+ generatePhiBackward(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.));
}
/**
* Calculate the matrix element for the splitting
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
*/
virtual DecayMEPtr matrixElement(const double z, const Energy2 t,
- const IdList & ids, const double phi, bool timeLike) {
+ const IdList & ids, const double phi, bool ) {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1)));
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
double sqrtmass = sqrt(m0t2*z*(1.-z)-m1t2*(1.-z)-m2t2*z+z*(1.-z));
// assign kernel
(*kernal)(0,0,0) = -phase*ghhv*sqrtmass/(1.-z); // 111
(*kernal)(1,0,0) = -sqrt(m2t2)*(1.+z)/sqrt(2.*(1.-z)); // 211 -> 411
(*kernal)(2,0,0) = cphase*ghhv*sqrtmass/(1.-z); // 311
// return the answer
return kernal;
}
protected:
/**
* Get the couplings
*/
- void getCouplings(Complex & g, const IdList & ids) const {
+ void getCouplings(Complex & g, const IdList & ) const {
g = Complex(_couplingValueRe,_couplingValueIm);
}
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 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 assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ZeroZeroOneEWSplitFn & operator=(const ZeroZeroOneEWSplitFn &) = delete;
private:
/**
* numerical value of the splitting coupling to be imported for BSM splittings
*/
double _couplingValueIm = 0.;
double _couplingValueRe = 0.;
};
}
#endif /* Herwig_ZeroZeroOneEWSplitFn_H */
diff --git a/Shower/QTilde/SplittingFunctions/ZeroZeroZeroEWSplitFn.h b/Shower/QTilde/SplittingFunctions/ZeroZeroZeroEWSplitFn.h
--- a/Shower/QTilde/SplittingFunctions/ZeroZeroZeroEWSplitFn.h
+++ b/Shower/QTilde/SplittingFunctions/ZeroZeroZeroEWSplitFn.h
@@ -1,231 +1,231 @@
// -*- C++ -*-
#ifndef Herwig_ZeroZeroZeroEWSplitFn_H
#define Herwig_ZeroZeroZeroEWSplitFn_H
//
// This is the declaration of the ZeroZeroZeroEWSplitFn class.
//
#include "Sudakov1to2FormFactor.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
namespace Herwig {
using namespace ThePEG;
/**
* The ZeroZeroZeroEWSplitFn class implements the splitting function for
* \f$\frac12\to q\frac12 h\f$ where the spin-0 higgs particle is a massive scalar boson.
*
* @see \ref ZeroZeroZeroEWSplitFnInterfaces "The interfaces"
* defined for ZeroZeroZeroEWSplitFn.
*/
class ZeroZeroZeroEWSplitFn: public Sudakov1to2FormFactor {
public:
/**
* Concrete implementation of the method to determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual bool accept(const IdList & ids) const;
/**
* Methods to return the splitting function.
*/
//@{
/**
* The concrete implementation of the overestimate of the splitting function,
* \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual double overestimateP(const double z, const IdList & ids) const {
Complex ghhh(0.,0.);
getCouplings(ghhh,ids);
return norm(ghhh)/(2.*z*(1.-z));
}
/**
* The concrete implementation of the
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,t)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
* @param rho The spin density matrix
*/
- virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
- const bool mass, const RhoDMatrix & rho) const {
+ virtual double ratioP(const double z, const Energy2 t, const IdList & ,
+ const bool , const RhoDMatrix & ) const {
return z*(1.-z)/t*GeV2;
}
/**
* The concrete implementation of the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double integOverP(const double z, const IdList & ids,
unsigned int PDFfactor=0) const;
/**
* The concrete implementation of the inverse of the indefinite integral.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor=0) const;
//@}
/**
* Method to calculate the azimuthal angle
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiForward(const double, const Energy2, const IdList &,
const RhoDMatrix &) {
// scalar so no dependence
return {{ {0, 1.} }};
}
/**
* Method to calculate the azimuthal angle for backward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiBackward(const double, const Energy2, const IdList &,
const RhoDMatrix &) {
// scalar so no dependence
assert(false);
return {{ {0, 1.} }};
}
/**
* Calculate the matrix element for the splitting
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
*/
virtual DecayMEPtr matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi, bool timeLike);
protected:
/**
* Get the couplings
* @param g The couplings ( for the SM case, this is g_SM/(e*mH^2) )
* @param ids The PDG codes for the particles in the splitting.
*/
void getCouplings(Complex & g, const IdList & ids) const;
/**
* Get the couplings with running masses
* @param g The couplings ( for the SM case, this is g_SM/(e*mH^2) )
* @param ids The PDG codes for the particles in the splitting.
* @param t The scale \f$t=2p_j\cdot p_k\f$.
*/
void getCouplings(Complex & g, const IdList & ids, const Energy2 t) const;
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 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;
//@}
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();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ZeroZeroZeroEWSplitFn & operator=(const ZeroZeroZeroEWSplitFn &) = delete;
private:
/**
* 1 / sin theta_w
*/
double gw_;
/**
* numerical value of the splitting coupling to be imported for BSM splittings
*/
double _couplingValueIm = 0.;
double _couplingValueRe = 0.;
/**
* Pointer to the SM object.
*/
tcHwSMPtr _theSM;
};
}
#endif /* Herwig_ZeroZeroZeroEWSplitFn_H */
diff --git a/Utilities/tests/utilitiesTestsKinematics.h b/Utilities/tests/utilitiesTestsKinematics.h
--- a/Utilities/tests/utilitiesTestsKinematics.h
+++ b/Utilities/tests/utilitiesTestsKinematics.h
@@ -1,154 +1,154 @@
// -*- C++ -*-
//
// utilitiesTestKinematics.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration, 2015 Marco A. Harrendorf
//
// 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_Utilities_Test_Kinematics_H
#define HERWIG_Utilities_Test_Kinematics_H
#include <boost/test/unit_test.hpp>
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Config/Unitsystem.h"
using namespace Herwig::Kinematics;
using namespace ThePEG::Units;
struct FixKinematics1 {
FixKinematics1()
{BOOST_TEST_MESSAGE( "setup fixture for utilitiesKinematicsTestTest" ); }
~FixKinematics1() { BOOST_TEST_MESSAGE( "teardown fixture for utilitiesKinematicsTest" ); }
};
/*
* Start of boost unit tests for Kinematics.h
*
* @todo Implement unit test for threeBodyDecay
*/
BOOST_AUTO_TEST_SUITE(utilitiesKinematicsTest)
/*
* Boost unit tests
*
*/
BOOST_AUTO_TEST_CASE(LorentzTest)
{
Energy mi,mj,mq;
Energy Pi,Pj,Pq;
- double Phi_i,Phi_j,Phi_q;
- double The_i,The_j,The_q;
+ double Phi_i,Phi_j;
+ double The_i,The_j;
Energy mRes_q;
- Energy P_Res_q;
- double Phi_Res_q;
- double The_Res_q;
+ // Energy P_Res_q;
+ // double Phi_Res_q;
+ // double The_Res_q;
using namespace ThePEG;
mi = 0.0*GeV;
Pi = 0.0*GeV;
Phi_i = 0.0;
The_i = 0.0;
Lorentz5Momentum pi(mi,Momentum3(Pi*cos(Phi_i)*sin(The_i),Pi*sin(Phi_i)*sin(The_i),Pi*cos(The_i)));
mj = 0.0*GeV;
Pj = 0.0*GeV;
Phi_j = 0.0;
The_j = 0.0;
Lorentz5Momentum pj(mj,Momentum3(Pj*cos(Phi_j)*sin(The_j),Pj*sin(Phi_j)*sin(The_j),Pj*cos(The_j)));
// Lorentz5Momentum pq(mq,Momentum3(Pj*cos(Phi_j)*sin(The_j),Pj*sin(Phi_j)*sin(The_j),Pj*cos(The_j)));
double eps=1e-14;
BOOST_CHECK_CLOSE(pi.mass()/GeV, pi.mass()/GeV, eps);
BOOST_CHECK_CLOSE(pi.m()/GeV , pi.m()/GeV, eps);
BOOST_CHECK_CLOSE(pi.x()/GeV , pi.x()/GeV, eps);
BOOST_CHECK_CLOSE(pi.y()/GeV , pi.y()/GeV, eps);
BOOST_CHECK_CLOSE(pi.z()/GeV , pi.z()/GeV, eps);
}
BOOST_AUTO_TEST_CASE(generateAnglesTest)
{
double flatMinusPiToPlusPi, flatNullToTwoPi;
for(int i = 0; i < 100; ++i) {
generateAngles(flatMinusPiToPlusPi, flatNullToTwoPi);
BOOST_CHECK( -M_PI <= flatMinusPiToPlusPi );
BOOST_CHECK( flatMinusPiToPlusPi <= M_PI);
BOOST_CHECK( 0. <= flatNullToTwoPi);
BOOST_CHECK(flatNullToTwoPi <= 2*M_PI);
}
}
BOOST_AUTO_TEST_CASE(unitDirectionTest)
{
using namespace Herwig::Kinematics;
BOOST_CHECK_EQUAL( unitDirection(1.1, -1), Axis() );
BOOST_CHECK_EQUAL( unitDirection(-1.1, -1), Axis() );
BOOST_CHECK_EQUAL( unitDirection(1.1, 0), Axis() );
BOOST_CHECK_EQUAL( unitDirection(-1.1, -1), Axis() );
BOOST_CHECK_EQUAL( unitDirection(1, 0), Axis(0, 0, 1) );
BOOST_CHECK_EQUAL( unitDirection(1, M_PI/2.), Axis(0, 0, 1) );
BOOST_CHECK(unitDirection(0, M_PI/2).almostEqual(Axis(0, 1, 0), 0.001) );
BOOST_CHECK_EQUAL( unitDirection(0, 0), Axis(1, 0, 0) );
}
BOOST_AUTO_TEST_CASE(pstarTwoBodyDecayTest)
{
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(-100*GeV), Energy(60*GeV), Energy(60*GeV))/GeV, Energy(0*GeV)/GeV);
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(-40*GeV), Energy(40*GeV))/GeV, Energy(0*GeV)/GeV);
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(-40*GeV), Energy(-40*GeV))/GeV, Energy(0*GeV)/GeV);
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(60*GeV), Energy(60*GeV))/GeV, Energy(0*GeV)/GeV);
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(50*GeV), Energy(50*GeV))/GeV, Energy(0*GeV)/GeV);
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(10*GeV), Energy(6*GeV), Energy(3*GeV))/GeV, Energy(std::sqrt(19*91)/20.*GeV)/GeV);
BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(10*GeV), Energy(3*GeV), Energy(6*GeV))/GeV, Energy(std::sqrt(19*91)/20.*GeV)/GeV);
}
BOOST_AUTO_TEST_CASE(twoBodyDecayTest1)
{
Lorentz5Momentum decayProductOne(GeV);
Lorentz5Momentum decayProductTwo(GeV);
BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(40*GeV), Axis(), decayProductOne, decayProductTwo));
BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), Axis(), decayProductOne, decayProductTwo));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(60*GeV), Energy(60*GeV), Axis(), decayProductOne, decayProductTwo)));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(-100*GeV), Energy(40*GeV), Energy(40*GeV), Axis(), decayProductOne, decayProductTwo)));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(-40*GeV), Energy(40*GeV), Axis(), decayProductOne, decayProductTwo)));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(-40*GeV), Axis(), decayProductOne, decayProductTwo)));
twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), Axis(1,0,0), decayProductOne, decayProductTwo);
BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(50*GeV)/GeV);
BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(50*GeV)/GeV);
twoBodyDecay(Lorentz5Momentum(10*GeV), Energy(6*GeV), Energy(3*GeV), Axis(1,0,0), decayProductOne, decayProductTwo);
BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(6*GeV, Momentum3(std::sqrt(19*91)/20.*GeV, ThePEG::ZERO, ThePEG::ZERO))/GeV);
BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(3*GeV, Momentum3(-(std::sqrt(19*91)/20.*GeV), ThePEG::ZERO, ThePEG::ZERO))/GeV);
}
BOOST_AUTO_TEST_CASE(twoBodyDecayTest2)
{
Lorentz5Momentum decayProductOne(GeV);
Lorentz5Momentum decayProductTwo(GeV);
BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo));
BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), 1, M_PI/2., decayProductOne, decayProductTwo));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(60*GeV), Energy(60*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(-100*GeV), Energy(40*GeV), Energy(40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(-40*GeV), Energy(40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(-40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), 1, M_PI/2., decayProductOne, decayProductTwo);
BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(50*GeV)/GeV);
BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(50*GeV)/GeV);
twoBodyDecay(Lorentz5Momentum(10*GeV), Energy(6*GeV), Energy(3*GeV), 1, M_PI/2., decayProductOne, decayProductTwo);
BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(6*GeV, Momentum3(ThePEG::ZERO, ThePEG::ZERO, std::sqrt(19*91)/20.*GeV))/GeV);
BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(3*GeV, Momentum3(ThePEG::ZERO, ThePEG::ZERO, -(std::sqrt(19*91)/20.*GeV)))/GeV);
}
BOOST_AUTO_TEST_SUITE_END()
#endif /* HERWIG_Utilities_Test_Kinematics_H */

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:56 PM (5 h, 30 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4478760
Default Alt Text
(235 KB)

Event Timeline