herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
HwRemDecayer.h
1// -*- C++ -*-
2//
3// HwRemDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4// Copyright (C) 2002-2019 The Herwig Collaboration
5//
6// Herwig is licenced under version 3 of the GPL, see COPYING for details.
7// Please respect the MCnet academic guidelines, see GUIDELINES for details.
8//
9#ifndef HERWIG_HwRemDecayer_H
10#define HERWIG_HwRemDecayer_H
11//
12// This is the declaration of the HwRemDecayer class.
13//
14
15#include "ThePEG/PDT/RemnantDecayer.h"
16#include "ThePEG/Handlers/EventHandler.h"
17#include "ThePEG/Repository/EventGenerator.h"
18#include "ThePEG/EventRecord/SubProcess.h"
19#include "ThePEG/PDF/BeamParticleData.h"
20#include "Herwig/Shower/ShowerAlpha.h"
21#include "Herwig/PDT/StandardMatchers.h"
23#include "HwRemDecayer.fh"
24
25namespace Herwig {
26using namespace ThePEG;
57
58public:
59
61 typedef vector<pair<tPPtr, tPPtr> > PartnerMap;
62
63public:
64
69 multiPeriph_(true), quarkPair_(false),
70 ptmin_(-1.*GeV), beta_(ZERO),
71 maxtrySoft_(10),
72 colourDisrupt_(1.0),
73 ladderbFactor_(0.0),
74 ladderPower_(-0.08),
75 ladderNorm_(1.0),
76 ladderMult_(1.0),
77 gaussWidth_(0.1),
78 valOfN_(0),
79 initTotRap_(0),
80 _kinCutoff(0.75*GeV),
81 _forcedSplitScale(2.5*GeV),
82 _range(1.1), _zbin(0.05),_ybin(0.),
83 _nbinmax(100), DISRemnantOpt_(0),
86
94 virtual bool accept(const DecayMode &) const {
95 return true;
96 }
97
102 virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const;
103
108 virtual bool multiCapable() const {
109 return true;
110 }
111
119 virtual ParticleVector decay(const DecayMode & dm, const Particle & p, Step & step) const;
121
122public:
123
129
137
143 void persistentInput(PersistentIStream & is, int version);
145
152 static void Init();
153
157 void initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
158 Energy forcedSplitScale);
159
165 void initSoftInteractions(Energy ptmin, InvEnergy2 beta);
166
174 void doSplit(pair<tPPtr, tPPtr> partons, pair<tcPDFPtr, tcPDFPtr> pdfs, bool first);
175
183 void finalize(double colourDisrupt=0.0, unsigned int softInt=0);
184
188 void findChildren(tPPtr,vector<PPtr> &) const;
189
190protected:
191
198 virtual IBPtr clone() const {return new_ptr(*this);}
199
204 virtual IBPtr fullclone() const {return new_ptr(*this);}
206
207protected:
208
216 virtual void doinit() {
217 Interfaced::doinit();
218 _ybin=0.25/_zbin;
219 mg_ = getParticleData(ParticleID::g)->constituentMass();
220 }
222
223private:
224
230
231public:
232
238
242 inline void extract(int id) {
243 for(unsigned int i=0; i<flav.size(); i++) {
244 if(id == sign*flav[i]){
245 if(hadron->id() == ParticleID::gamma ||
247 hadron->id() == ParticleID::reggeon) {
248 flav[0] = id;
249 flav[1] = -id;
250 extracted = 0;
251 flav.resize(2);
252 }
253 else if (hadron->id() == ParticleID::pomeron && pomeronStructure==0) {
254 extracted = 0;
255 }
256 else {
257 extracted = i;
258 }
259 break;
260 }
261 }
262 }
263
268 long RemID() const;
269
274 bool isSeaQuark(tcPPtr parton) const {
275 return ((parton->id() != ParticleID::g) && ( !isValenceQuark(parton) ) );
276 }
277
281 bool isValenceQuark(tcPPtr parton) const {
282 return isValenceQuark(parton->id());
283 }
284
289 bool isSeaQuarkData(tcPDPtr partonData) const {
290 return ((partonData->id() != ParticleID::g) && ( !isValenceQuarkData(partonData) ) );
291 }
292
296 bool isValenceQuarkData(tcPDPtr partonData) const {
297 int id(sign*partonData->id());
298 return find(flav.begin(),flav.end(),id) != flav.end();
299 }
300
304 bool isValenceQuark(int id) const {
305 return find(flav.begin(),flav.end(),sign*id) != flav.end();
306 }
307
309 vector<int> flav;
310
313
315 int sign;
316
319
321 unsigned int pomeronStructure;
322 };
323
327 const pair<HadronContent, HadronContent>& content() const {
328 return theContent;
329 }
330
335
340 theContent.first = getHadronContent(beam.first);
341 theContent.second = getHadronContent(beam.second);
342 }
343
344private:
345
358 void split(tPPtr parton, HadronContent & content, tRemPPtr rem,
359 Lorentz5Momentum & used, PartnerMap & partners, tcPDFPtr pdf, bool first);
360
367 void mergeColour(tPPtr p1, tPPtr p2, bool anti) const;
368
375 void fixColours(PartnerMap partners, bool anti, double disrupt) const;
376
380 void setRemMasses(PPair diquarks) const;
381
387 PPtr finalSplit(const tRemPPtr rem, long remID,
388 Lorentz5Momentum usedMomentum) const {
389 // Create the remnant and set its momentum, also reset all of the decay
390 // products from the hadron
391 PPtr remnant = new_ptr(Particle(getParticleData(remID)));
392 Lorentz5Momentum prem(rem->momentum()-usedMomentum);
393 prem.setMass(getParticleData(remID)->constituentMass());
394 prem.rescaleEnergy();
395 remnant->set5Momentum(prem);
396 // Add the remnant to the step, but don't do colour connections
397 thestep->addDecayProduct(rem,remnant,false);
398 return remnant;
399 }
400
401
414 PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx,
415 Lorentz5Momentum &pf, Lorentz5Momentum &p,
416 HadronContent & content) const;
417
422 bool isPartonic(tPPtr parton) const;
423
426
430 Energy softPt() const;
431
436 void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2,
437 Lorentz5Momentum &g1, Lorentz5Momentum &g2) const;
438
442 void doSoftInteractions(unsigned int N){
443 if(!multiPeriph_){
444 doSoftInteractions_old(N);} //outdated model for soft interactions
445 else{
446 doSoftInteractions_multiPeriph(N); // Multiperipheral model
447 }
448 }
449
453 void doSoftInteractions_old(unsigned int N);
454
458 void doSoftInteractions_multiPeriph(unsigned int N);
459
463 bool doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy energy, unsigned int &its)
464 const;
465
469 LorentzRotation rotate(const LorentzMomentum &p) const;
470
475 template <typename T>
476 inline T gaussDistribution(T mean, T stdev) const{
477 double x = rnd();
478 x = sqrt(-2.*log(x));
479 double y;
480 randAzm(x,x,y);
481 return mean + stdev*x;
482 }
483
484
493 inline double randUng(double A, double B) const{
494 double prun;
495 if(A == 0.) prun = 0.;
496 else prun = 1./(1.+B*1.2533/A);
497 if(rnd() < prun) return 2.*(rnd()-0.5)*A;
498 else {
499 double temp = gaussDistribution(0.,B);
500 if(temp < 0) return temp - abs(A);
501 else return temp + abs(A);
502 }
503 }
504 template <typename T>
505 inline void randAzm(T pt, T &px, T &py) const{
506 double c,s,cs;
507 while(true) {
508 c = 2.*rnd()-1.;
509 s = 2.*rnd()-1.;
510 cs = c*c+s*s;
511 if(cs <= 1.&&cs!=0.) break;
512 }
513 T qt = pt/cs;
514 px = (c*c-s*s)*qt;
515 py = 2.*c*s*qt;
516 }
517
518 inline Energy randExt(Energy AM0,InvEnergy B) const{
519 double r = rnd();
520 // Starting value
521 Energy am = AM0-log(r)/B;
522 for(int i = 1; i<20; ++i) {
523 double a = exp(-B*(am-AM0))/(1.+B*AM0);
524 double f = (1.+B*am)*a-r;
525 InvEnergy df = -B*B*am*a;
526 Energy dam = -f/df;
527 am += dam;
528 if(am<AM0) am = AM0 + .001*MeV;
529 if(abs(dam) < .001*MeV) break;
530 }
531 return am;
532 }
533
540 tPPtr addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const;
542
547 pair<bool, bool> theanti;
548
552 pair<double, double> theX;
553
555 pair<HadronContent, HadronContent> theContent;
556
558 pair<Lorentz5Momentum, Lorentz5Momentum> theUsed;
559
564 pair<PartnerMap, PartnerMap> theMaps;
565
572
577 pair<RemPPtr, RemPPtr> theRems;
578
583
587 mutable Ptr<BeamParticleData>::const_pointer theBeamData;
588
592 mutable tcPDFPtr _pdf;
593
594private:
595
600
605
610
615
618
623
627 Energy ptmin_;
628
632 InvEnergy2 beta_;
633
638 unsigned int maxtrySoft_;
639
645
651
656
661
662 double ladderMult_;
669
674 double valOfN_;
675
681
683
686
691
696
700 double _range;
701
705 double _zbin;
706
710 double _ybin;
711
716
720 ShowerAlphaPtr _alphaS;
721
725 ShowerAlphaPtr _alphaEM;
726
730 unsigned int DISRemnantOpt_;
731
735 unsigned int PtDistribution_;
736
740 unsigned int pomeronStructure_;
742
746 Energy mg_;
747
748};
749
750
751}
752
753#endif /* HERWIG_HwRemDecayer_H */
The HwRemDecayer class is responsible for the decay of the remnants.
Definition: HwRemDecayer.h:56
double valOfN_
Variable to store the current total multiplicity of a ladder.
Definition: HwRemDecayer.h:674
Energy softPt() const
Produce pt values according to dN/dp_T = N p_T exp(-beta_*p_T^2)
void split(tPPtr parton, HadronContent &content, tRemPPtr rem, Lorentz5Momentum &used, PartnerMap &partners, tcPDFPtr pdf, bool first)
Do the forced Splitting of the Remnant with respect to the extracted parton parton.
double colourDisrupt_
Variable to store the relative number of colour disrupted connections to additional soft subprocesses...
Definition: HwRemDecayer.h:644
double initTotRap_
Variable to store the initial total rapidity between of the remnants.
Definition: HwRemDecayer.h:680
virtual bool accept(const DecayMode &) const
Check if this decayer can perfom the decay specified by the given decay mode.
Definition: HwRemDecayer.h:94
double randUng(double A, double B) const
This returns a random number with a flat distribution [-A,A] plus gaussian tail with stdev B TODO: Sh...
Definition: HwRemDecayer.h:493
unsigned int DISRemnantOpt_
Option for the DIS remnant.
Definition: HwRemDecayer.h:730
pair< RemPPtr, RemPPtr > theRems
Pair of Remnant pointers.
Definition: HwRemDecayer.h:577
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
virtual bool multiCapable() const
Return true if this decayed can extract more than one parton from a particle.
Definition: HwRemDecayer.h:108
Energy _forcedSplitScale
The PDF freezing scale as set in ShowerHandler.
Definition: HwRemDecayer.h:695
tPPair softRems_
Pair of soft Remnant pointers, i.e.
Definition: HwRemDecayer.h:622
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
void doSoftInteractions_old(unsigned int N)
Create N soft gluon interactions (old version)
tcPDFPtr _pdf
The PDF for the current initial-state shower.
Definition: HwRemDecayer.h:592
vector< pair< tPPtr, tPPtr > > PartnerMap
Typedef to store information about colour partners.
Definition: HwRemDecayer.h:61
void setRemMasses(PPair diquarks) const
Set the momenta of the Remnants properly and boost the decay particles.
Energy mg_
The gluon constituent mass.
Definition: HwRemDecayer.h:746
tcPPtr theBeam
The beam particle data for the current incoming hadron.
Definition: HwRemDecayer.h:582
double ladderPower_
Variable of the parameterization of the ladder multiplicity.
Definition: HwRemDecayer.h:655
pair< double, double > theX
variable to sum up the x values of the extracted particles
Definition: HwRemDecayer.h:552
tPPtr addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const
Method to add a particle to the step.
bool isPartonic(tPPtr parton) const
Check if a particle is a parton from a hadron or not.
HadronContent getHadronContent(tcPPtr hadron) const
Return a HadronContent struct from a PPtr to a hadron.
double _ybin
Size of the bins in y for the interpolation.
Definition: HwRemDecayer.h:710
double ladderNorm_
Variable of the parameterization of the ladder multiplicity.
Definition: HwRemDecayer.h:660
T gaussDistribution(T mean, T stdev) const
Methods to generate random distributions also all stolen form UA5Handler.
Definition: HwRemDecayer.h:476
HwRemDecayer & operator=(const HwRemDecayer &)=delete
The assignment operator is private and must never be called.
StepPtr thestep
Variable to hold a pointer to the current step.
Definition: HwRemDecayer.h:571
int _nbinmax
Maximum number of bins for the z interpolation.
Definition: HwRemDecayer.h:715
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
Definition: HwRemDecayer.h:204
HwRemDecayer()
The default constructor.
Definition: HwRemDecayer.h:68
PPtr finalSplit(const tRemPPtr rem, long remID, Lorentz5Momentum usedMomentum) const
This creates a parton from the remaining flavours of the hadron.
Definition: HwRemDecayer.h:387
InvEnergy2 beta_
slope of the soft pt-spectrum: dN/dp_T = N p_T exp(-beta*p_T^2)
Definition: HwRemDecayer.h:632
const pair< HadronContent, HadronContent > & content() const
Return the hadron content objects for the incoming particles.
Definition: HwRemDecayer.h:327
ShowerAlphaPtr _alphaS
Pointer to the object calculating the QCD coupling.
Definition: HwRemDecayer.h:720
LorentzRotation rotate(const LorentzMomentum &p) const
This returns the rotation matrix needed to rotate p into the z axis.
void setHadronContent(tPPair beam)
Set the hadron contents.
Definition: HwRemDecayer.h:339
unsigned int pomeronStructure_
Option for the treatment of the pomeron structure.
Definition: HwRemDecayer.h:740
void doSoftInteractions(unsigned int N)
Create N soft gluon interactions.
Definition: HwRemDecayer.h:442
virtual IBPtr clone() const
Make a simple clone of this object.
Definition: HwRemDecayer.h:198
void initSoftInteractions(Energy ptmin, InvEnergy2 beta)
Initialize the soft scattering machinery.
virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const
Return true if this decayer can handle the extraction of the extracted parton from the given parti...
virtual ParticleVector decay(const DecayMode &dm, const Particle &p, Step &step) const
Perform a decay for a given DecayMode and a given Particle instance.
virtual void doinit()
Initialize this object after the setup phase before saving an EventGenerator to disk.
Definition: HwRemDecayer.h:216
Energy ptmin_
ptcut of the UE model
Definition: HwRemDecayer.h:627
pair< HadronContent, HadronContent > theContent
Pair of HadronContent structs to know about the quark content of the beams.
Definition: HwRemDecayer.h:555
void fixColours(PartnerMap partners, bool anti, double disrupt) const
Set the colour connections.
double _zbin
Size of the bins in z for the interpolation.
Definition: HwRemDecayer.h:705
void finalize(double colourDisrupt=0.0, unsigned int softInt=0)
Perform the final creation of the diquarks.
void mergeColour(tPPtr p1, tPPtr p2, bool anti) const
Merge the colour lines of two particles.
PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx, Lorentz5Momentum &pf, Lorentz5Momentum &p, HadronContent &content) const
This takes the particle and find a splitting for np -> p + child and creates the correct kinematics a...
double _range
Range for emission.
Definition: HwRemDecayer.h:700
void doSoftInteractions_multiPeriph(unsigned int N)
Create N soft gluon interactions with multiperhpheral kinematics.
void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2, Lorentz5Momentum &g1, Lorentz5Momentum &g2) const
Get the 2 pairs of 5Momenta for the scattering.
static void Init()
The standard Init function used to initialize the interfaces.
double ladderbFactor_
Variable to store the additive factor of the multiperipheral ladder multiplicity.
Definition: HwRemDecayer.h:650
bool allowLeptons_
Switch to control handling of charged leptons in proton.
Definition: HwRemDecayer.h:604
pair< bool, bool > theanti
A flag which indicates, whether the extracted valence quark was a anti particle.
Definition: HwRemDecayer.h:547
Ptr< BeamParticleData >::const_pointer theBeamData
the beam data
Definition: HwRemDecayer.h:587
unsigned int maxtrySoft_
Maximum number of attempts for the regeneration of an additional soft scattering, before the number o...
Definition: HwRemDecayer.h:638
pair< Lorentz5Momentum, Lorentz5Momentum > theUsed
Pair of Lorentz5Momentum to keep track of the forced splitting product momenta.
Definition: HwRemDecayer.h:558
void initialize(pair< tRemPPtr, tRemPPtr > rems, tPPair beam, Step &step, Energy forcedSplitScale)
Do several checks and initialization, for remnantdecay inside ShowerHandler.
pair< PartnerMap, PartnerMap > theMaps
Pair of PartnerMap's to store the particles, which will be colour connected in the end.
Definition: HwRemDecayer.h:564
bool doPhaseSpaceGenerationGluons(vector< Lorentz5Momentum > &softGluons, Energy energy, unsigned int &its) const
Phase space generation for the ladder partons.
bool quarkPair_
True if kinematics is to be calculated for quarks.
Definition: HwRemDecayer.h:614
double gaussWidth_
Variable to store the gaussian width of the fluctuation of the longitudinal momentum fraction.
Definition: HwRemDecayer.h:668
bool multiPeriph_
Switch to control using multiperipheral kinemaics.
Definition: HwRemDecayer.h:609
ShowerAlphaPtr _alphaEM
Pointer to the object calculating the QED coupling.
Definition: HwRemDecayer.h:725
unsigned int PtDistribution_
Option for the pT generation.
Definition: HwRemDecayer.h:735
Energy _kinCutoff
The kinematic cut-off.
Definition: HwRemDecayer.h:690
bool allowTop_
Switch to control handling of top quarks in proton.
Definition: HwRemDecayer.h:599
void findChildren(tPPtr, vector< PPtr > &) const
Find the children.
void doSplit(pair< tPPtr, tPPtr > partons, pair< tcPDFPtr, tcPDFPtr > pdfs, bool first)
Perform the acual forced splitting.
bool used() const
PDPtr getParticleData(PID) const
tcPPtr parent(tcRemPPtr remnant) const
-*- C++ -*-
Qty< 0, 1, 0 > Energy
ThePEG::Ptr< Particle >::transient_pointer tPPtr
pair< PPtr, PPtr > PPair
double sqrt(int x)
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
ThePEG::Ptr< Particle >::transient_const_pointer tcPPtr
ThePEG::Ptr< PDFBase >::transient_const_pointer tcPDFPtr
ThePEG::Ptr< Step >::pointer StepPtr
pair< tPPtr, tPPtr > tPPair
ThePEG::Ptr< Particle >::pointer PPtr
vector< PPtr > ParticleVector
constexpr ZeroUnit ZERO
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
Iterator find(IteratorRange< Iterator > r, const T &t)
struct that is used to catch exceptions which are thrown due to energy conservation issues of additio...
Definition: HwRemDecayer.h:128
Simple struct to store info about baryon quark and di-quark constituents.
Definition: HwRemDecayer.h:237
bool isSeaQuark(tcPPtr parton) const
Method to determine whether parton is a quark from the sea.
Definition: HwRemDecayer.h:274
int sign
-1 if the particle is an anti-particle.
Definition: HwRemDecayer.h:315
long RemID() const
Return a proper particle ID assuming that id has been removed from the hadron.
tcPDPtr hadron
The ParticleData objects of the hadron.
Definition: HwRemDecayer.h:318
int extracted
The array index of the extracted particle.
Definition: HwRemDecayer.h:312
bool isSeaQuarkData(tcPDPtr partonData) const
Method to determine whether parton is a quark from the sea.
Definition: HwRemDecayer.h:289
void extract(int id)
manually extract the valence flavour id.
Definition: HwRemDecayer.h:242
unsigned int pomeronStructure
Pomeron treatment.
Definition: HwRemDecayer.h:321
bool isValenceQuark(tcPPtr parton) const
Method to determine whether parton is a valence quark.
Definition: HwRemDecayer.h:281
bool isValenceQuarkData(tcPDPtr partonData) const
Method to determine whether parton is a valence quark.
Definition: HwRemDecayer.h:296
bool isValenceQuark(int id) const
Method to determine whether parton is a valence quark.
Definition: HwRemDecayer.h:304
vector< int > flav
The valence flavours of the corresponding baryon.
Definition: HwRemDecayer.h:309