diff --git a/.hgtags b/.hgtags --- a/.hgtags +++ b/.hgtags @@ -1,40 +1,41 @@ 168ae2110e964d62fbc1331a1c2e095952a67748 release-2-5-2 3abb4fa42e20e332796c2572334c2d77204cd0e0 release-2-4-2 4796ca080aafd5daa3b7349b015cb1df944428a2 release-2-5-0 76da042f056eb153981b4d005d5474ffb90a5e88 release-2-4-1 81a684a558413c69df314365eabf09893ffd43d8 release-2-6-0 bd75cd00d99f4bdbaed992daf98f0a73c0f91e9b release-2-4-0 ff6ecc8d49ce10299303b050394bd5cb5837f1c3 release-2-5-1 d0389f5453b2c210923e1adc7b872b18269de668 release-2-6-1 f8998033021185942533b824607285feb3fbd2dc release-2-6-1a cead23e428b9aacaf2d709e722624e54f844498b release-2-6-1b 191db4655439045f912cb21bd905e729d59ec7bc release-2-6-2 edb538156e9c3d64bb842934b4cebf0126aeb9ea release-2-6-3 eb4a104591859ecac18746b1ad54d6aa0c2a5d1a release-2-7-0 568971ac5b3c1d044c9259f2280a8304fc5a62e9 trunk-before-QED 6e3edb6cfeb4ee48687eb4eb3d016026fc59d602 trunk-after-QED 633abb80b571aa23088957df60e9b0000bbb8a22 release-2-7-1 1bdde095d2346c15ee548e5406a96f0fc6d6e0f1 beforeHQ a0f9fb821396092bdbeee532bcb0bd624f58335b before_MB_merge 270c1e6b34aa7f758f1d9868c4d3e1ec4bf4e709 herwig-7-0-0 6e0f198c1c2603ecd1a0b6cfe40105cda4bd58c5 herwig-7-0-1 566c1de845a8070559cda45b1bdb40afa18cb2cc herwig-7-0-2 f5c4aa956880f2def763ebd57de7b5bfa55cb1db herwig-7-0-3 65282dedfc2e4bec184e68678dbf4c553c968f38 herwig-7-0-4 541e7790b65ed423c86780bf66ec30e6b99b5a18 herwig-7-1-0 dd35a1c12d57c047169e8c5fb18644972d49c6ac herwig-7-1-1 0d651b079756b63713e32a1341d81e4dfc7eeb7b herwig-7-1-2 4b97934bc41c861c4be04f563ffa68a94a982560 herwig-7-1-3 97aca5398cfa1f3273804f03fa96fa0fa23eca61 herwig-7-1-4 3d69fbe18c682c98891c5f9204947f2eb7a72686 herwig-7-1-5 392e0bdc94f11d067dc24792a4b470b09eb8fdf7 herwig-7-2-0 392e0bdc94f11d067dc24792a4b470b09eb8fdf7 herwig-7-2-0 af22cb052ed5e3fd5323ec9a24693962efe8144d herwig-7-2-0 f3047b8819217a3ea264adf423ab183c9ca9b3a1 herwig-7-1-6 f3047b8819217a3ea264adf423ab183c9ca9b3a1 herwig-7-1-6 c9519d355aeab14d43e4655faf78e2c7324b34a6 herwig-7-1-6 af22cb052ed5e3fd5323ec9a24693962efe8144d herwig-7-2-0 51f480cba67c56c2b6d5a65383254d6aafcd5cbd herwig-7-2-0 b76dedb1f1f1bc42d676bca00019852b32003502 herwig-7-2-1 474898bfc7f68305c654d6323f844c6cbd0921fc herwig-7.2.2 6178a87dde19c2547c33d44c66f17e429977cab2 herwig-7-2-3 +9c12a19ed7fd602882ce332fe35ec0ef55d18d6e herwig-7-3-0 diff --git a/Decay/FormFactors/CzyzNucleonFormFactor.cc b/Decay/FormFactors/CzyzNucleonFormFactor.cc --- a/Decay/FormFactors/CzyzNucleonFormFactor.cc +++ b/Decay/FormFactors/CzyzNucleonFormFactor.cc @@ -1,339 +1,342 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the CzyzNucleonFormFactor class. // #include "CzyzNucleonFormFactor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/ResonanceHelpers.h" using namespace Herwig; CzyzNucleonFormFactor::CzyzNucleonFormFactor() { // masses and widths of rho resonances rhoMasses_ = {775.49*MeV, 1465*MeV, 1720*MeV, 2.12*GeV,2.32647*GeV}; rhoWidths_ = {149.10*MeV, 400*MeV, 250*MeV, 0.3 *GeV,0.4473 *GeV}; // masses and width of omega resonances omegaMasses_ = {782.65*MeV, 1425*MeV, 1670*MeV, 2.0707 *GeV, 2.34795 *GeV}; omegaWidths_ = {8.49 *MeV, 215*MeV, 315*MeV, 1.03535*GeV, 1.173975*GeV}; // c_1 couplings c1Re_ = {1.,-0.467,-0.177, 0.301}; c1Im_ = {0.,-0.385, 0.149, 0.264}; // c_2 couplings c2Re_ = {1., 0.0521, -0.00308,-0.348}; c2Im_ = {0.,-3.04, 2.38,-0.104}; // c_3 couplings c3Re_ = {1.,-7.88,10.2}; c3Im_ = {0., 5.67, -1.94}; // c_4 couplings c4Re_ = {1.,-0.832, 0.405}; c4Im_ = {0., 0.308,-0.25}; // Magnetic moments mup_ = 2.792847356; mun_ = -1.9130427; + // initialise a_ and b_ + a_ = mup_-mun_-1.; + b_ = -mup_-mun_+1.; // set up the form factors addFormFactor(2212,2212,2,2,2,1,1,1); addFormFactor(2112,2112,2,2,2,1,1,1); initialModes(numberOfFactors()); } IBPtr CzyzNucleonFormFactor::clone() const { return new_ptr(*this); } IBPtr CzyzNucleonFormFactor::fullclone() const { return new_ptr(*this); } void CzyzNucleonFormFactor::doinit() { BaryonFormFactor::doinit(); static const Complex ii(0.,1.); // calculate c_1 c1_.clear(); assert(c1Re_.size()==4 && c1Im_.size()==4); // c_1 1 -> 4 complex fact(ZERO); for(unsigned int ix=0;ix<4;++ix) { c1_.push_back(c1Re_[ix]+ii*c1Im_[ix]); fact += c1_[ix]*sqr(omegaMasses_[ix]); } c1_.push_back(-fact/sqr(omegaMasses_[4])); // calculate c_2 c2_.clear(); assert(c2Re_.size()==4 && c2Im_.size()==4); // c_2 1 -> 4 fact = ZERO; for(unsigned int ix=0;ix<4;++ix) { c2_.push_back(c2Re_[ix]+ii*c2Im_[ix]); fact += c2_[ix]*sqr(rhoMasses_[ix]); } c2_.push_back(-fact/sqr(rhoMasses_[4])); // calculate c_3 c3_.clear(); assert(c3Re_.size()==3 && c3Im_.size()==3); // c_3 1 -> 4 fact = ZERO; complex fact2(ZERO); for(unsigned int ix=0;ix<3;++ix) { c3_.push_back(c3Re_[ix]+ii*c3Im_[ix]); fact += c3_[ix]*sqr(omegaMasses_[ix]); fact2 += c3_[ix]*sqr(omegaMasses_[ix])* (sqr(omegaMasses_[ix])-sqr(omegaMasses_[4]) +ii*(omegaMasses_[4]*omegaWidths_[4]-omegaMasses_[ix]*omegaWidths_[ix])); } c3_.push_back(fact2/sqr(omegaMasses_[3])/ (sqr(omegaMasses_[4])-sqr(omegaMasses_[3]) +ii*(omegaMasses_[3]*omegaWidths_[3]-omegaMasses_[4]*omegaWidths_[4])) ); fact += c3_[3]*sqr(omegaMasses_[3]); c3_.push_back(-fact/sqr(omegaMasses_[4])); // c_4 1 -> 4 fact = ZERO; fact2 = ZERO; for(unsigned int ix=0;ix<3;++ix) { c4_.push_back(c4Re_[ix]+ii*c4Im_[ix]); fact += c4_[ix]*sqr(rhoMasses_[ix]); fact2 += c4_[ix]*sqr(rhoMasses_[ix])* (sqr(rhoMasses_[ix])-sqr(rhoMasses_[4]) +ii*(rhoMasses_[4]*rhoWidths_[4]-rhoMasses_[ix]*rhoWidths_[ix])); } c4_.push_back(fact2/sqr(rhoMasses_[3])/ (sqr(rhoMasses_[4])-sqr(rhoMasses_[3]) +ii*(rhoMasses_[3]*rhoWidths_[3]-rhoMasses_[4]*rhoWidths_[4])) ); fact += c4_[3]*sqr(rhoMasses_[3]); c4_.push_back(-fact/sqr(rhoMasses_[4])); // a and b parameters - a_ = mup_-mun_-1; + a_ = mup_-mun_-1.; b_ = -mup_-mun_+1.; } void CzyzNucleonFormFactor::persistentOutput(PersistentOStream & os) const { os << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV) << ounit(omegaMasses_,GeV) << ounit(omegaWidths_,GeV) << c1Re_ << c1Im_ << c2Re_ << c2Im_ << c3Re_ << c3Im_ << c4Re_ << c4Im_ << c1_ << c2_ << c3_ << c4_ << mup_ << mun_ << a_ << b_; } void CzyzNucleonFormFactor::persistentInput(PersistentIStream & is, int) { is >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV) >> iunit(omegaMasses_,GeV) >> iunit(omegaWidths_,GeV) >> c1Re_ >> c1Im_ >> c2Re_ >> c2Im_ >> c3Re_ >> c3Im_ >> c4Re_ >> c4Im_ >> c1_ >> c2_ >> c3_ >> c4_ >> mup_ >> mun_ >> a_ >> b_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigCzyzNucleonFormFactor("Herwig::CzyzNucleonFormFactor", "HwFormFactors.so"); void CzyzNucleonFormFactor::Init() { static ClassDocumentation documentation ("The CzyzNucleonFormFactor class implements the model of " "Phys.Rev. D90 (2014) no.11, 114021 for the nucleon form factor", "The nucleon form factor model of \\cite{Czyz:2014sha} was used", "\\bibitem{Czyz:2014sha}\n" "H.~Czyż, J.~H.~Kühn and S.~Tracz,\n" "%``Nucleon form factors and final state radiative corrections to $e^+e^- \\to p\\bar{p}γ$,''\n" "Phys.\\ Rev.\\ D {\\bf 90} (2014) no.11, 114021\n" "doi:10.1103/PhysRevD.90.114021\n" "[arXiv:1407.7995 [hep-ph]].\n" "%%CITATION = doi:10.1103/PhysRevD.90.114021;%%\n" "%5 citations counted in INSPIRE as of 25 Aug 2018\n"); - + static ParVector interfaceRhoMasses ("RhoMasses", "The masses of the rho mesons for the form factor", &CzyzNucleonFormFactor::rhoMasses_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); - + static ParVector interfaceRhoWidths ("RhoWidths", "The widths of the rho mesons for the form factor", &CzyzNucleonFormFactor::rhoWidths_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); - + static ParVector interfaceOmegaMasses ("OmegaMasses", "The masses of the omega mesons for the form factor", &CzyzNucleonFormFactor::omegaMasses_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); - + static ParVector interfaceOmegaWidths ("OmegaWidths", "The widths of the omega mesons for the form factor", &CzyzNucleonFormFactor::omegaWidths_, GeV, 5, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static ParVector interfacec1Real ("c1Real", "The real part of the c_1 coupling", &CzyzNucleonFormFactor::c1Re_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec1Imag ("c1Imag", "The imaginary part of the c_1 coupling", &CzyzNucleonFormFactor::c1Im_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec2Real ("c2Real", "The real part of the c_2 coupling", &CzyzNucleonFormFactor::c2Re_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec2Imag ("c2Imag", "The imaginary part of the c_2 coupling", &CzyzNucleonFormFactor::c2Im_, 4, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec3Real ("c3Real", "The real part of the c_3 coupling", &CzyzNucleonFormFactor::c3Re_, 3, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec3Imag ("c3Imag", "The imaginary part of the c_3 coupling", &CzyzNucleonFormFactor::c3Im_, 3, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec4Real ("c4Real", "The real part of the c_4 coupling", &CzyzNucleonFormFactor::c4Re_, 3, 1., -100., 100.0, false, false, Interface::limited); static ParVector interfacec4Imag ("c4Imag", "The imaginary part of the c_4 coupling", &CzyzNucleonFormFactor::c4Im_, 3, 1., -100., 100.0, false, false, Interface::limited); static Parameter interfaceMuProton ("MuProton", "The proton magnetic moment", &CzyzNucleonFormFactor::mup_, 2.792, 0.0, 10.0, false, false, Interface::limited); static Parameter interfaceMuNeutron ("MuNeutron", "The proton magnetic moment", &CzyzNucleonFormFactor::mun_, -1.913, 0.0, 10.0, false, false, Interface::limited); } void CzyzNucleonFormFactor:: SpinHalfSpinHalfFormFactor(Energy2 q2,int iloc, int id0,int id1,Energy,Energy, Complex & f1v,Complex & f2v,Complex & f3v, Complex & f1a,Complex & f2a,Complex & f3a, FlavourInfo flavour, Virtuality virt) { f1a = f2a = f3a = f1v = f2v = f3v = 0.; assert(abs(id0)==abs(id1)); if(iloc==0) assert(abs(id0)==2212); else assert(abs(id0)==2112); assert(virt==TimeLike); if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return; if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return; if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return; bool I0 = true, I1 = true; if(flavour.I!=IsoSpin::IUnknown) { if(flavour.I==IsoSpin::IZero) { I1=false; if(flavour.I3!=IsoSpin::I3Unknown and flavour.I3!=IsoSpin::I3Zero) return; } else if(flavour.I==IsoSpin::IOne) { I0=false; if(flavour.I3!=IsoSpin::I3Unknown and flavour.I3!=IsoSpin::I3Zero) return; } } if(flavour.I3!=IsoSpin::I3Unknown and flavour.I3!=IsoSpin::I3Zero) return; // calculate the form factors Complex F1S(0.),F1V(0.),F2S(0.),F2V(0.); Complex n1(0.),n2(0.),n3(0.),n4(0.); for(unsigned int ix=0;ix<5;++ix) { if(I0) { F1S += c1_[ix]*Resonance::BreitWignerFW(q2,omegaMasses_[ix],omegaWidths_[ix]); F2S += c3_[ix]*Resonance::BreitWignerFW(q2,omegaMasses_[ix],omegaWidths_[ix]); } if(I1) { F1V += c2_[ix]*Resonance::BreitWignerFW(q2, rhoMasses_[ix], rhoWidths_[ix]); F2V += c4_[ix]*Resonance::BreitWignerFW(q2, rhoMasses_[ix], rhoWidths_[ix]); } n1 += c1_[ix]; n2 += c2_[ix]; n3 += c3_[ix]; n4 += c4_[ix]; } F1S *= 0.5 /n1; F1V *= 0.5 /n2; F2S *= -0.5*b_/n3; F2V *= 0.5*a_/n4; if(iloc==0) { f1v = F1S + F1V; f2v = F2S + F2V; } else { f1v = F1S - F1V; f2v = F2S - F2V; } } void CzyzNucleonFormFactor:: dataBaseOutput(ofstream& output,bool header, bool create) const { if(header) output << "update decayers set parameters=\""; - if(create) output << "create Herwig::CzyzNucleonFormFactor " + if(create) output << "create Herwig::CzyzNucleonFormFactor " << name() << " \n"; for(unsigned int ix=0;ix<5;++ix) { output << "newdef " << name() << ":RhoMasses " << ix << " " << rhoMasses_[ix]/GeV << "\n"; output << "newdef " << name() << ":RhoWidths " << ix << " " << rhoWidths_[ix]/GeV << "\n"; output << "newdef " << name() << ":OmegaMasses " << ix << " " << omegaMasses_[ix]/GeV << "\n"; output << "newdef " << name() << ":OmegaWidths " << ix << " " << omegaWidths_[ix]/GeV << "\n"; } for(unsigned int ix=0;ix<4;++ix) { output << "newdef " << name() << ":c1Real " << ix << " " << c1Re_[ix] << "\n"; output << "newdef " << name() << ":c1Imag " << ix << " " << c1Im_[ix] << "\n"; output << "newdef " << name() << ":c2Real " << ix << " " << c2Re_[ix] << "\n"; output << "newdef " << name() << ":c2Imag " << ix << " " << c2Im_[ix] << "\n"; } for(unsigned int ix=0;ix<3;++ix) { output << "newdef " << name() << ":c3Real " << ix << " " << c3Re_[ix] << "\n"; output << "newdef " << name() << ":c3Imag " << ix << " " << c3Im_[ix] << "\n"; output << "newdef " << name() << ":c4Real " << ix << " " << c4Re_[ix] << "\n"; output << "newdef " << name() << ":c4Imag " << ix << " " << c4Im_[ix] << "\n"; } output << "newdef " << name() << ":MuProton " << mup_ << "\n"; output << "newdef " << name() << ":MuNeutron " << mun_ << "\n"; BaryonFormFactor::dataBaseOutput(output,false,false); - if(header) output << "\n\" where BINARY ThePEGName=\"" + if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,1187 +1,1187 @@ // -*- 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 #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" #include "Cluster.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include using namespace Herwig; DescribeClass describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so"); ClusterFissioner::ClusterFissioner() : _clMaxLight(3.35*GeV), _clMaxBottom(3.35*GeV), _clMaxCharm(3.35*GeV), _clMaxExotic(3.35*GeV), _clPowLight(2.0), _clPowBottom(2.0), _clPowCharm(2.0), _clPowExotic(2.0), _pSplitLight(1.0), _pSplitBottom(1.0), _pSplitCharm(1.0), _pSplitExotic(1.0), _fissionPwtUquark(1), _fissionPwtDquark(1), _fissionPwtSquark(0.5), _fissionCluster(0), _kinematicThresholdChoice(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)) {} IBPtr ClusterFissioner::clone() const { return new_ptr(*this); } IBPtr ClusterFissioner::fullclone() const { return new_ptr(*this); } void ClusterFissioner::persistentOutput(PersistentOStream & os) const { os << _hadronSelector << ounit(_clMaxLight,GeV) << ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom << _clPowCharm << _clPowExotic << _pSplitLight << _pSplitBottom << _pSplitCharm << _pSplitExotic << _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark << ounit(_btClM,GeV) << _iopRem << _kinematicThresholdChoice << ounit(_kappa, GeV/meter) << _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure << _probPowFactor << _probShift << ounit(_kinThresholdShift,sqr(GeV)); } void ClusterFissioner::persistentInput(PersistentIStream & is, int) { is >> _hadronSelector >> iunit(_clMaxLight,GeV) >> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom >> _clPowCharm >> _clPowExotic >> _pSplitLight >> _pSplitBottom >> _pSplitCharm >> _pSplitExotic >> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark >> iunit(_btClM,GeV) >> _iopRem >> _kinematicThresholdChoice >> iunit(_kappa, GeV/meter) >> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure >> _probPowFactor >> _probShift >> iunit(_kinThresholdShift,sqr(GeV)); } void ClusterFissioner::Init() { static ClassDocumentation documentation ("Class responsibles for chopping up the clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::ClusterFissioner::_hadronSelector, false, false, true, false); // ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])", &ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])", &ClusterFissioner::_clMaxBottom, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])", &ClusterFissioner::_clMaxCharm, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])", &ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); // ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks", &ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks", &ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks", &ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter 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 interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks", &ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks", &ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks", &ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks", &ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false); static Switch 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 selector class.", 0); static SwitchOption interfaceFissionNew (interfaceFission, "new", "Alternative cluster fission which does not depend on the hadron selector class", 1); static Parameter interfaceFissionPwtUquark ("FissionPwtUquark", "Weight for fission in U quarks", &ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtDquark ("FissionPwtDquark", "Weight for fission in D quarks", &ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtSquark ("FissionPwtSquark", "Weight for fission in S quarks", &ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch 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 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 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 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 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 interfaceFissionMassScale ("FissionMassScale", "Cluster fission mass scale", &ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); static Parameter interfaceProbPowFactor ("ProbablityPowerFactor", "Power factor in ClausterFissioner bell probablity function", &ClusterFissioner::_probPowFactor, 2.0, 1.0, 20.0, false, false, Interface::limited); static Parameter interfaceProbShift ("ProbablityShift", "Shifts from the center in ClausterFissioner bell probablity function", &ClusterFissioner::_probShift, 0.0, -10.0, 10.0, false, false, Interface::limited); - // static Parameter interfaceKineticThresholdShift - // ("KineticThresholdShift", - // "Shifts from the kinetic threshold in ClausterFissioner", - // &ClusterFissioner::_kinThresholdShift, sqr(GeV), 0.*sqr(GeV), -10.0*sqr(GeV), 10.0*sqr(GeV), - // false, false, Interface::limited); + static Parameter interfaceKineticThresholdShift + ("KineticThresholdShift", + "Shifts from the kinetic threshold in ClausterFissioner", + &ClusterFissioner::_kinThresholdShift, sqr(GeV), 0.*sqr(GeV), -10.0*sqr(GeV), 10.0*sqr(GeV), + false, false, Interface::limited); static Switch 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); } 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 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 & 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(ct.first.first); ClusterPtr two = dynamic_ptr_cast(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' 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); } } } namespace { /** * Check if can't make a hadron from the partons */ bool cantMakeHadron(tcPPtr p1, tcPPtr p2) { return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr()); } /** * Check if can't make a diquark from the partons */ bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) { long id1 = p1->id(), id2 = p2->id(); return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0); } } 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; // pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2); 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); pair res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2, ptrQ1, pQ1, newPtr1, pQone, newPtr2, pQtwo, ptrQ2, pQ2); // derive the masses of the children Mc1 = res.first; Mc2 = res.second; // static kinematic threshold if(_kinematicThresholdChoice == 0) { if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) 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; } /************************** * 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 = _hadronSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); } toHadron2 = _hadronSelector->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 = _hadronSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2); if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); } } else if(!toHadron2) { toHadron2 = _hadronSelector->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) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { swap(newPtr1,newPtr2); } // check again if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(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); pair res = drawNewMasses(mmax, soft1, soft2, pClu1, pClu2, ptrQ[iq1], pQ1, newPtr1, pQone, newPtr2, pQtwo, ptrQ[iq1], pQ2); Mc1 = res.first; Mc2 = res.second; if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue; // check if need to force meson clster to hadron toHadron = _hadronSelector->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 = _hadronSelector->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 = _hadronSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2); if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); } } else if(!toDiQuark) { Mc2 = _hadronSelector->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 = (_hadronSelector->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; } void ClusterFissioner::drawNewFlavour(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. double prob_d; double prob_u; double prob_s; switch(_fissionCluster){ case 0: prob_d = _hadronSelector->pwt(1); prob_u = _hadronSelector->pwt(2); prob_s = _hadronSelector->pwt(3); break; case 1: prob_d = _fissionPwtDquark; prob_u = _fissionPwtUquark; prob_s = _fissionPwtSquark; 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()); } void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) 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. double prob_d; double prob_u; double prob_s = 0.; double scale = abs(double(sqr(_m0Fission)/mass2)); // Choose which splitting weights you wish to use switch(_fissionCluster){ // 0: ClusterFissioner and ClusterDecayer use the same weights case 0: prob_d = _hadronSelector->pwt(1); prob_u = _hadronSelector->pwt(2); /* Strangeness enhancement: Case 1: probability scaling Case 2: Exponential scaling */ if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_hadronSelector->pwt(3),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 = _fissionPwtDquark; prob_u = _fissionPwtUquark; if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,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::ProbablityFunction(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::isHeavy(tcClusterPtr clu) { // default double clpow = _clPowLight; Energy clmax = _clMaxLight; // particle data for constituents tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()}; for(unsigned int ix=0;ixnumComponents(),3);++ix) { cptr[ix]=clu->particle(ix)->dataPtr(); } // different parameters for exotic, bottom and charm clusters if(CheckId::isExotic(cptr[0],cptr[1],cptr[1])) { clpow = _clPowExotic; clmax = _clMaxExotic; } else if(CheckId::hasBottom(cptr[0],cptr[1],cptr[1])) { clpow = _clPowBottom; clmax = _clMaxBottom; } else if(CheckId::hasCharm(cptr[0],cptr[1],cptr[1])) { clpow = _clPowCharm; clmax = _clMaxCharm; } //regular checks // 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 = ProbablityFunction(scale,threshold); scale = clu->mass()/GeV; threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV; canSplitMinimally = ProbablityFunction(scale,threshold); } return aboveCutoff && canSplitMinimally; } diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in --- a/src/defaults/Hadronization.in +++ b/src/defaults/Hadronization.in @@ -1,133 +1,133 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default hadronization # # There are no user servicable parts inside. # # Anything that follows below should only be touched if you # know what you're doing. ############################################################# cd /Herwig/Particles create ThePEG::ParticleData Cluster setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1 create ThePEG::ParticleData Remnant setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1 mkdir /Herwig/Hadronization cd /Herwig/Hadronization create Herwig::ClusterHadronizationHandler ClusterHadHandler create Herwig::PartonSplitter PartonSplitter create Herwig::ClusterFinder ClusterFinder create Herwig::ColourReconnector ColourReconnector create Herwig::ClusterFissioner ClusterFissioner create Herwig::LightClusterDecayer LightClusterDecayer create Herwig::ClusterDecayer ClusterDecayer create Herwig::HwppSelector HadronSelector newdef ClusterHadHandler:PartonSplitter PartonSplitter newdef ClusterHadHandler:ClusterFinder ClusterFinder newdef ClusterHadHandler:ColourReconnector ColourReconnector newdef ClusterHadHandler:ClusterFissioner ClusterFissioner newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer newdef ClusterHadHandler:ClusterDecayer ClusterDecayer newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2 newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter newdef ClusterHadHandler:UnderlyingEventHandler NULL newdef ClusterFissioner:HadronSelector HadronSelector newdef LightClusterDecayer:HadronSelector HadronSelector newdef ClusterDecayer:HadronSelector HadronSelector newdef ClusterFinder:HadronSelector HadronSelector # ColourReconnector Default Parameters newdef ColourReconnector:ColourReconnection Yes newdef ColourReconnector:Algorithm Baryonic # Statistical CR Parameters: newdef ColourReconnector:AnnealingFactor 0.9 newdef ColourReconnector:AnnealingSteps 50 newdef ColourReconnector:TriesPerStepFactor 5.0 newdef ColourReconnector:InitialTemperature 0.1 # Plain and Baryonic CR Paramters newdef ColourReconnector:ReconnectionProbability 0.95 newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 # BaryonicMesonic and BaryonicMesonic CR Paramters newdef ColourReconnector:ReconnectionProbability3Mto3M 0.5 newdef ColourReconnector:ReconnectionProbability3MtoBBbar 0.5 newdef ColourReconnector:ReconnectionProbabilityBbarBto3M 0.5 newdef ColourReconnector:ReconnectionProbability2Bto2B 0.05 newdef ColourReconnector:ReconnectionProbabilityMBtoMB 0.5 newdef ColourReconnector:StepFactor 1.0 newdef ColourReconnector:MesonToBaryonFactor 1.333 # General Parameters and switches newdef ColourReconnector:MaxDistance 1.0e50 newdef ColourReconnector:OctetTreatment All newdef ColourReconnector:CR2BeamClusters No newdef ColourReconnector:Junction Yes newdef ColourReconnector:PrePlainCR No newdef ColourReconnector:LocalCR No newdef ColourReconnector:CausalCR No # Debugging newdef ColourReconnector:Debug No # set ClusterFissioner parameters set /Herwig/Hadronization/ClusterFissioner:KinematicThreshold Dynamic -set /Herwig/Hadronization/ClusterFissioner:ProbablityPowerFactor 3.3462 -set /Herwig/Hadronization/ClusterFissioner:ProbablityShift 6.0442 -# set /Herwig/Hadronization/ClusterFissioner:KineticThresholdShift -1.5321 +set /Herwig/Hadronization/ClusterFissioner:ProbablityPowerFactor 6.486 +set /Herwig/Hadronization/ClusterFissioner:ProbablityShift -0.87875 +# set /Herwig/Hadronization/ClusterFissioner:KineticThresholdShift 0.08844 # Clustering parameters for light quarks -newdef ClusterFissioner:ClMaxLight 3.2344 -newdef ClusterFissioner:ClPowLight 2.6464 -newdef ClusterFissioner:PSplitLight 0.7235 +newdef ClusterFissioner:ClMaxLight 3.528693 +newdef ClusterFissioner:ClPowLight 1.849375 +newdef ClusterFissioner:PSplitLight 0.914156 newdef ClusterDecayer:ClDirLight 1 newdef ClusterDecayer:ClSmrLight 0.78 # Clustering parameters for b-quarks newdef ClusterFissioner:ClMaxBottom 3.757 newdef ClusterFissioner:ClPowBottom 0.547 newdef ClusterFissioner:PSplitBottom 0.625 newdef ClusterDecayer:ClDirBottom 1 newdef ClusterDecayer:ClSmrBottom 0.078 newdef HadronSelector:SingleHadronLimitBottom 0.000 # Clustering parameters for c-quarks newdef ClusterFissioner:ClMaxCharm 3.950 newdef ClusterFissioner:ClPowCharm 2.559 newdef ClusterFissioner:PSplitCharm 0.994 newdef ClusterDecayer:ClDirCharm 1 newdef ClusterDecayer:ClSmrCharm 0.163 newdef HadronSelector:SingleHadronLimitCharm 0.000 # Clustering parameters for exotic quarks # (e.g. hadronizing Susy particles) newdef ClusterFissioner:ClMaxExotic 2.7*GeV newdef ClusterFissioner:ClPowExotic 1.46 newdef ClusterFissioner:PSplitExotic 1.00 newdef ClusterDecayer:ClDirExotic 1 newdef ClusterDecayer:ClSmrExotic 0. newdef HadronSelector:SingleHadronLimitExotic 0. # newdef PartonSplitter:SplitPwtSquark 0.824135 newdef PartonSplitter:Split uds # newdef HadronSelector:PwtDquark 1.0 newdef HadronSelector:PwtUquark 1.0 -newdef HadronSelector:PwtSquark 0.35743 +newdef HadronSelector:PwtSquark 0.374094 newdef HadronSelector:PwtCquark 0.0 newdef HadronSelector:PwtBquark 0.0 -newdef HadronSelector:PwtDIquark 0.36452 -newdef HadronSelector:SngWt 0.88022 -newdef HadronSelector:DecWt 0.34622 +newdef HadronSelector:PwtDIquark 0.33107 +newdef HadronSelector:SngWt 0.89050 +newdef HadronSelector:DecWt 0.41628 newdef HadronSelector:Mode 1 newdef HadronSelector:BelowThreshold All create Herwig::SpinHadronizer SpinHadronizer diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in --- a/src/defaults/Shower.in +++ b/src/defaults/Shower.in @@ -1,464 +1,464 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default parton shower # # Useful switches for users are marked near the top of # this file. # # Don't edit this file directly, but reset the switches # in your own input files! ############################################################ library HwMPI.so library HwShower.so library HwMatching.so mkdir /Herwig/Shower cd /Herwig/Shower create Herwig::QTildeShowerHandler ShowerHandler newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer # use LO PDFs for Shower, can be changed later newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF ##################################### # initial setup, don't change these! ##################################### create Herwig::SplittingGenerator SplittingGenerator create Herwig::ShowerAlphaQCD AlphaQCD create Herwig::ShowerAlphaQCD AlphaQCDFSR create Herwig::ShowerAlphaQED AlphaQED set AlphaQED:CouplingSource Thompson create Herwig::ShowerAlphaQED AlphaEW set AlphaEW:CouplingSource MZ create Herwig::PartnerFinder PartnerFinder newdef PartnerFinder:PartnerMethod 0 newdef PartnerFinder:ScaleChoice 0 create Herwig::KinematicsReconstructor KinematicsReconstructor newdef KinematicsReconstructor:ReconstructionOption Colour3 newdef KinematicsReconstructor:InitialStateReconOption SofterFraction newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost newdef KinematicsReconstructor:FinalFinalWeight Yes newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED newdef ShowerHandler:PartnerFinder PartnerFinder newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef ShowerHandler:SplittingGenerator SplittingGenerator newdef ShowerHandler:Interactions QEDQCD newdef ShowerHandler:SpinCorrelations Yes newdef ShowerHandler:SoftCorrelations Singular ################################################################## # Intrinsic pT # # Recommended: # 1.9 GeV for Tevatron W/Z production. # 2.1 GeV for LHC W/Z production at 10 TeV # 2.2 GeV for LHC W/Z production at 14 TeV # # Set all parameters to 0 to disable ################################################################## newdef ShowerHandler:IntrinsicPtGaussian 2.2*GeV newdef ShowerHandler:IntrinsicPtBeta 0 newdef ShowerHandler:IntrinsicPtGamma 0*GeV newdef ShowerHandler:IntrinsicPtIptmax 0*GeV ############################################################# # Set up truncated shower handler. ############################################################# create Herwig::PowhegShowerHandler PowhegShowerHandler set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PartnerFinder PartnerFinder newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator newdef PowhegShowerHandler:Interactions QEDQCD newdef PowhegShowerHandler:SpinCorrelations Yes newdef PowhegShowerHandler:SoftCorrelations Singular newdef PowhegShowerHandler:IntrinsicPtGaussian 2.2*GeV newdef PowhegShowerHandler:IntrinsicPtBeta 0 newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV newdef PowhegShowerHandler:EvolutionScheme DotProduct ############################################################# # End of interesting user servicable section. # # Anything that follows below should only be touched if you # know what you're doing. # # Really. ############################################################# # # a few default values newdef ShowerHandler:MECorrMode 1 newdef ShowerHandler:EvolutionScheme DotProduct newdef AlphaQCD:ScaleFactor 1.0 newdef AlphaQCD:NPAlphaS 2 newdef AlphaQCD:Qmin 0.935 newdef AlphaQCD:NumberOfLoops 2 newdef AlphaQCD:AlphaIn 0.1185 newdef AlphaQCDFSR:ScaleFactor 1.0 newdef AlphaQCDFSR:NPAlphaS 2 newdef AlphaQCDFSR:Qmin 0.935 newdef AlphaQCDFSR:NumberOfLoops 2 -newdef AlphaQCDFSR:AlphaIn 0.1186 +newdef AlphaQCDFSR:AlphaIn 0.102337 # # # Lets set up all the splittings create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn set QtoQGammaSplitFn:InteractionType QED set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral set QtoQGammaSplitFn:AngularOrdered Yes set QtoQGammaSplitFn:StrictAO Yes create Herwig::HalfHalfOneSplitFn QtoQGSplitFn newdef QtoQGSplitFn:InteractionType QCD newdef QtoQGSplitFn:ColourStructure TripletTripletOctet set QtoQGSplitFn:AngularOrdered Yes set QtoQGSplitFn:StrictAO Yes create Herwig::OneOneOneSplitFn GtoGGSplitFn newdef GtoGGSplitFn:InteractionType QCD newdef GtoGGSplitFn:ColourStructure OctetOctetOctet set GtoGGSplitFn:AngularOrdered Yes set GtoGGSplitFn:StrictAO Yes create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn newdef WtoWGammaSplitFn:InteractionType QED newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral set WtoWGammaSplitFn:AngularOrdered Yes set WtoWGammaSplitFn:StrictAO Yes create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn newdef GtoQQbarSplitFn:InteractionType QCD newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet set GtoQQbarSplitFn:AngularOrdered Yes set GtoQQbarSplitFn:StrictAO Yes create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn newdef GammatoQQbarSplitFn:InteractionType QED newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged set GammatoQQbarSplitFn:AngularOrdered Yes set GammatoQQbarSplitFn:StrictAO Yes create Herwig::HalfOneHalfSplitFn QtoGQSplitFn newdef QtoGQSplitFn:InteractionType QCD newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet set QtoGQSplitFn:AngularOrdered Yes set QtoGQSplitFn:StrictAO Yes create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn newdef QtoGammaQSplitFn:InteractionType QED newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged set QtoGammaQSplitFn:AngularOrdered Yes set QtoGammaQSplitFn:StrictAO Yes create Herwig::HalfHalfOneEWSplitFn QtoQWZSplitFn newdef QtoQWZSplitFn:InteractionType EW newdef QtoQWZSplitFn:ColourStructure EW create Herwig::HalfHalfOneEWSplitFn LtoLWZSplitFn newdef LtoLWZSplitFn:InteractionType EW newdef LtoLWZSplitFn:ColourStructure EW create Herwig::HalfHalfZeroEWSplitFn QtoQHSplitFn newdef QtoQHSplitFn:InteractionType EW newdef QtoQHSplitFn:ColourStructure EW create Herwig::OneOneOneEWSplitFn VtoVVSplitFn newdef VtoVVSplitFn:InteractionType EW newdef VtoVVSplitFn:ColourStructure EW create Herwig::OneOneOneQEDSplitFn GammatoWWSplitFn newdef GammatoWWSplitFn:InteractionType QED newdef GammatoWWSplitFn:ColourStructure ChargedNeutralCharged create Herwig::OneOneZeroEWSplitFn VtoVHSplitFn newdef VtoVHSplitFn:InteractionType EW newdef VtoVHSplitFn:ColourStructure EW # # Now the Sudakovs create Herwig::PTCutOff PTCutOff -newdef PTCutOff:pTmin 0.7473*GeV +newdef PTCutOff:pTmin 0.654714*GeV create Herwig::SudakovFormFactor SudakovCommon newdef SudakovCommon:Alpha AlphaQCD newdef SudakovCommon:Cutoff PTCutOff newdef SudakovCommon:PDFmax 1.0 cp SudakovCommon QtoQGSudakov newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn newdef QtoQGSudakov:PDFmax 1.9 cp SudakovCommon QtoQGammaSudakov set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn set QtoQGammaSudakov:Alpha AlphaQED set QtoQGammaSudakov:PDFmax 1.9 cp QtoQGammaSudakov LtoLGammaSudakov cp PTCutOff LtoLGammaPTCutOff # Technical parameter to stop evolution. set LtoLGammaPTCutOff:pTmin 0.000001 set LtoLGammaSudakov:Cutoff LtoLGammaPTCutOff cp SudakovCommon QtoQWZSudakov set QtoQWZSudakov:SplittingFunction QtoQWZSplitFn set QtoQWZSudakov:Alpha AlphaEW set QtoQWZSudakov:PDFmax 1.9 cp SudakovCommon LtoLWZSudakov set LtoLWZSudakov:SplittingFunction LtoLWZSplitFn set LtoLWZSudakov:Alpha AlphaEW set LtoLWZSudakov:PDFmax 1.9 cp SudakovCommon QtoQHSudakov set QtoQHSudakov:SplittingFunction QtoQHSplitFn set QtoQHSudakov:Alpha AlphaEW set QtoQHSudakov:PDFmax 1.9 cp QtoQHSudakov LtoLHSudakov cp SudakovCommon VtoVVSudakov set VtoVVSudakov:SplittingFunction VtoVVSplitFn set VtoVVSudakov:Alpha AlphaQED cp SudakovCommon GammatoWWSudakov set GammatoWWSudakov:SplittingFunction GammatoWWSplitFn set GammatoWWSudakov:Alpha AlphaEW set GammatoWWSudakov:PDFmax 1.9 cp SudakovCommon VtoVHSudakov set VtoVHSudakov:SplittingFunction VtoVHSplitFn set VtoVHSudakov:Alpha AlphaEW set VtoVHSudakov:PDFmax 1.9 cp SudakovCommon GtoGGSudakov newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn newdef GtoGGSudakov:PDFmax 2.0 cp SudakovCommon WtoWGammaSudakov newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn set WtoWGammaSudakov:Alpha AlphaQED cp SudakovCommon GtoQQbarSudakov newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GammatoQQbarSudakov newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn set GammatoQQbarSudakov:Alpha AlphaQED newdef GammatoQQbarSudakov:PDFmax 10000.0 cp SudakovCommon GtobbbarSudakov newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtobbbarSudakov:PDFmax 40000.0 cp SudakovCommon GtoccbarSudakov newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoccbarSudakov:PDFmax 2000.0 cp SudakovCommon QtoGQSudakov newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn cp SudakovCommon QtoGammaQSudakov newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn set QtoGammaQSudakov:Alpha AlphaQED newdef QtoGammaQSudakov:PDFmax 2000.0 cp SudakovCommon utoGuSudakov newdef utoGuSudakov:SplittingFunction QtoGQSplitFn newdef utoGuSudakov:PDFFactor OverOneMinusZ newdef utoGuSudakov:PDFmax 5.0 cp SudakovCommon dtoGdSudakov newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn newdef dtoGdSudakov:PDFFactor OverOneMinusZ #Use a different Sudakov for FS QCD splittings in order to use a different value of alphaS cp QtoQGSudakov QtoQGSudakovFSR cp GtoGGSudakov GtoGGSudakovFSR cp GtoQQbarSudakov GtoQQbarSudakovFSR cp GtobbbarSudakov GtobbbarSudakovFSR cp GtobbbarSudakov GtoccbarSudakovFSR set QtoQGSudakovFSR:Alpha AlphaQCDFSR set GtoGGSudakovFSR:Alpha AlphaQCDFSR set GtoQQbarSudakovFSR:Alpha AlphaQCDFSR set GtobbbarSudakovFSR:Alpha AlphaQCDFSR set GtoccbarSudakovFSR:Alpha AlphaQCDFSR # # Now add the final splittings # do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakovFSR do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakovFSR do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakovFSR do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakovFSR do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakovFSR do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakovFSR # do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakovFSR # do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakovFSR do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakovFSR do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakovFSR do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakovFSR do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakovFSR do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakovFSR # do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov # do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov # # Now lets add the initial splittings. Remember the form a->b,c; means # that the current particle b is given and we backward branch to new # particle a which is initial state and new particle c which is final state # do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov # #do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov #do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov #do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov #do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov #do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->e-,e+; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->mu-,mu+; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->tau-,tau+; GammatoQQbarSudakov # do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov # do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov # # Electroweak # do SplittingGenerator:AddFinalSplitting u->u,Z0; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting d->d,Z0; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting s->s,Z0; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting c->c,Z0; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting b->b,Z0; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting t->t,Z0; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting u->u,Z0; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting d->d,Z0; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting s->s,Z0; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting c->c,Z0; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting b->b,Z0; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting t->t,Z0; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting u->d,W+; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting c->s,W+; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting d->u,W-; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting s->c,W-; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting u->d,W+; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting c->s,W+; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting d->u,W-; QtoQWZSudakov do SplittingGenerator:AddInitialSplitting s->c,W-; QtoQWZSudakov do SplittingGenerator:AddFinalSplitting c->c,h0; QtoQHSudakov do SplittingGenerator:AddFinalSplitting b->b,h0; QtoQHSudakov do SplittingGenerator:AddFinalSplitting t->t,h0; QtoQHSudakov do SplittingGenerator:AddInitialSplitting c->c,h0; QtoQHSudakov do SplittingGenerator:AddInitialSplitting b->b,h0; QtoQHSudakov do SplittingGenerator:AddInitialSplitting t->t,h0; QtoQHSudakov do SplittingGenerator:AddFinalSplitting gamma->W+,W-; GammatoWWSudakov do SplittingGenerator:AddFinalSplitting Z0->W+,W-; VtoVVSudakov do SplittingGenerator:AddFinalSplitting W+->W+,gamma; VtoVVSudakov do SplittingGenerator:AddFinalSplitting W+->W+,Z0; VtoVVSudakov do SplittingGenerator:AddFinalSplitting W+->W+,h0; VtoVHSudakov do SplittingGenerator:AddFinalSplitting Z0->Z0,h0; VtoVHSudakov # # Electroweak l -> l V # #do SplittingGenerator:AddFinalSplitting e-->e-,Z0; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting mu-->mu-,Z0; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting tau-->tau-,Z0; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting nu_e->nu_e,Z0; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting nu_mu->nu_mu,Z0; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting nu_tau->nu_tau,Z0; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting e-->nu_e,W-; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting mu-->nu_mu,W-; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting tau-->nu_tau,W-; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting nu_e->e-,W+; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting nu_mu->mu-,W+; LtoLWZSudakov #do SplittingGenerator:AddFinalSplitting nu_tau->tau-,W+; LtoLWZSudakov