herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
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 
25 namespace Herwig {
26 using namespace ThePEG;
57 
58 public:
59 
61  typedef vector<pair<tPPtr, tPPtr> > PartnerMap;
62 
63 public:
64 
68  HwRemDecayer() : allowTop_(false), multiPeriph_(true), quarkPair_(false),
69  ptmin_(-1.*GeV), beta_(ZERO),
70  maxtrySoft_(10),
71  colourDisrupt_(1.0),
72  ladderbFactor_(0.0),
73  ladderPower_(-0.08),
74  ladderNorm_(1.0),
75  ladderMult_(1.0),
76  gaussWidth_(0.1),
77  valOfN_(0),
78  initTotRap_(0),
79  _kinCutoff(0.75*GeV),
80  _forcedSplitScale(2.5*GeV),
81  _range(1.1), _zbin(0.05),_ybin(0.),
82  _nbinmax(100), DISRemnantOpt_(0),
83  PtDistribution_(0),
84  pomeronStructure_(0), mg_(ZERO) {}
85 
93  virtual bool accept(const DecayMode &) const {
94  return true;
95  }
96 
101  virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const;
102 
107  virtual bool multiCapable() const {
108  return true;
109  }
110 
118  virtual ParticleVector decay(const DecayMode & dm, const Particle & p, Step & step) const;
120 
121 public:
122 
128 
135  void persistentOutput(PersistentOStream & os) const;
136 
142  void persistentInput(PersistentIStream & is, int version);
144 
151  static void Init();
152 
156  void initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
157  Energy forcedSplitScale);
158 
164  void initSoftInteractions(Energy ptmin, InvEnergy2 beta);
165 
173  void doSplit(pair<tPPtr, tPPtr> partons, pair<tcPDFPtr, tcPDFPtr> pdfs, bool first);
174 
182  void finalize(double colourDisrupt=0.0, unsigned int softInt=0);
183 
187  void findChildren(tPPtr,vector<PPtr> &) const;
188 
189 protected:
190 
197  virtual IBPtr clone() const {return new_ptr(*this);}
198 
203  virtual IBPtr fullclone() const {return new_ptr(*this);}
205 
206 protected:
207 
215  virtual void doinit() {
217  _ybin=0.25/_zbin;
218  mg_ = getParticleData(ParticleID::g)->constituentMass();
219  }
221 
222 private:
223 
228  HwRemDecayer & operator=(const HwRemDecayer &);
229 
230 public:
231 
236  struct HadronContent {
237 
241  inline void extract(int id) {
242  for(unsigned int i=0; i<flav.size(); i++) {
243  if(id == sign*flav[i]){
244  if(hadron->id() == ParticleID::gamma ||
245  (hadron->id() == ParticleID::pomeron && pomeronStructure==1) ||
246  hadron->id() == ParticleID::reggeon) {
247  flav[0] = id;
248  flav[1] = -id;
249  extracted = 0;
250  flav.resize(2);
251  }
252  else if (hadron->id() == ParticleID::pomeron && pomeronStructure==0) {
253  extracted = 0;
254  }
255  else {
256  extracted = i;
257  }
258  break;
259  }
260  }
261  }
262 
267  long RemID() const;
268 
273  bool isSeaQuark(tcPPtr parton) const {
274  return ((parton->id() != ParticleID::g) && ( !isValenceQuark(parton) ) );
275  }
276 
280  bool isValenceQuark(tcPPtr parton) const {
281  return isValenceQuark(parton->id());
282  }
283 
288  bool isSeaQuarkData(tcPDPtr partonData) const {
289  return ((partonData->id() != ParticleID::g) && ( !isValenceQuarkData(partonData) ) );
290  }
291 
295  bool isValenceQuarkData(tcPDPtr partonData) const {
296  int id(sign*partonData->id());
297  return find(flav.begin(),flav.end(),id) != flav.end();
298  }
299 
303  bool isValenceQuark(int id) const {
304  return find(flav.begin(),flav.end(),sign*id) != flav.end();
305  }
306 
308  vector<int> flav;
309 
312 
314  int sign;
315 
318 
320  unsigned int pomeronStructure;
321  };
322 
326  const pair<HadronContent, HadronContent>& content() const {
327  return theContent;
328  }
329 
333  HadronContent getHadronContent(tcPPtr hadron) const;
334 
339  theContent.first = getHadronContent(beam.first);
340  theContent.second = getHadronContent(beam.second);
341  }
342 
343 private:
344 
357  void split(tPPtr parton, HadronContent & content, tRemPPtr rem,
358  Lorentz5Momentum & used, PartnerMap & partners, tcPDFPtr pdf, bool first);
359 
366  void mergeColour(tPPtr p1, tPPtr p2, bool anti) const;
367 
374  void fixColours(PartnerMap partners, bool anti, double disrupt) const;
375 
379  void setRemMasses() const;
380 
386  PPtr finalSplit(const tRemPPtr rem, long remID,
387  Lorentz5Momentum usedMomentum) const {
388  // Create the remnant and set its momentum, also reset all of the decay
389  // products from the hadron
390  PPtr remnant = new_ptr(Particle(getParticleData(remID)));
391  Lorentz5Momentum prem(rem->momentum()-usedMomentum);
392  prem.setMass(getParticleData(remID)->constituentMass());
393  prem.rescaleEnergy();
394  remnant->set5Momentum(prem);
395  // Add the remnant to the step, but don't do colour connections
396  thestep->addDecayProduct(rem,remnant,false);
397  return remnant;
398  }
399 
400 
413  PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx,
414  Lorentz5Momentum &pf, Lorentz5Momentum &p,
415  HadronContent & content) const;
416 
421  bool isPartonic(tPPtr parton) const;
422 
425 
429  Energy softPt() const;
430 
435  void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2,
436  Lorentz5Momentum &g1, Lorentz5Momentum &g2) const;
437 
441  void doSoftInteractions(unsigned int N){
442  if(!multiPeriph_){
443  doSoftInteractions_old(N);} //outdated model for soft interactions
444  else{
445  doSoftInteractions_multiPeriph(N); // Multiperipheral model
446  }
447  }
448 
452  void doSoftInteractions_old(unsigned int N);
453 
457  void doSoftInteractions_multiPeriph(unsigned int N);
458 
462  bool doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy energy, unsigned int &its)
463  const;
464 
468  LorentzRotation rotate(const LorentzMomentum &p) const;
469 
474  template <typename T>
475  inline T gaussDistribution(T mean, T stdev) const{
476  double x = rnd();
477  x = sqrt(-2.*log(x));
478  double y;
479  randAzm(x,x,y);
480  return mean + stdev*x;
481  }
482 
483 
492  inline double randUng(double A, double B) const{
493  double prun;
494  if(A == 0.) prun = 0.;
495  else prun = 1./(1.+B*1.2533/A);
496  if(rnd() < prun) return 2.*(rnd()-0.5)*A;
497  else {
498  double temp = gaussDistribution(0.,B);
499  if(temp < 0) return temp - abs(A);
500  else return temp + abs(A);
501  }
502  }
503  template <typename T>
504  inline void randAzm(T pt, T &px, T &py) const{
505  double c,s,cs;
506  while(true) {
507  c = 2.*rnd()-1.;
508  s = 2.*rnd()-1.;
509  cs = c*c+s*s;
510  if(cs <= 1.&&cs!=0.) break;
511  }
512  T qt = pt/cs;
513  px = (c*c-s*s)*qt;
514  py = 2.*c*s*qt;
515  }
516 
517  inline Energy randExt(Energy AM0,InvEnergy B) const{
518  double r = rnd();
519  // Starting value
520  Energy am = AM0-log(r)/B;
521  for(int i = 1; i<20; ++i) {
522  double a = exp(-B*(am-AM0))/(1.+B*AM0);
523  double f = (1.+B*am)*a-r;
524  InvEnergy df = -B*B*am*a;
525  Energy dam = -f/df;
526  am += dam;
527  if(am<AM0) am = AM0 + .001*MeV;
528  if(abs(dam) < .001*MeV) break;
529  }
530  return am;
531  }
532 
539  tPPtr addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const;
541 
546  pair<bool, bool> theanti;
547 
551  pair<double, double> theX;
552 
554  pair<HadronContent, HadronContent> theContent;
555 
557  pair<Lorentz5Momentum, Lorentz5Momentum> theUsed;
558 
563  pair<PartnerMap, PartnerMap> theMaps;
564 
571 
576  pair<RemPPtr, RemPPtr> theRems;
577 
581  mutable tcPPtr theBeam;
582 
586  mutable Ptr<BeamParticleData>::const_pointer theBeamData;
587 
591  mutable tcPDFPtr _pdf;
592 
593 private:
594 
598  bool allowTop_;
599 
604 
609 
612 
617 
621  Energy ptmin_;
622 
626  InvEnergy2 beta_;
627 
632  unsigned int maxtrySoft_;
633 
639 
645 
649  double ladderPower_;
650 
654  double ladderNorm_;
655 
656  double ladderMult_;
662  double gaussWidth_;
663 
668  double valOfN_;
669 
674  double initTotRap_;
675 
677 
680 
684  Energy _kinCutoff;
685 
690 
694  double _range;
695 
699  double _zbin;
700 
704  double _ybin;
705 
709  int _nbinmax;
710 
714  ShowerAlphaPtr _alphaS;
715 
719  ShowerAlphaPtr _alphaEM;
720 
724  unsigned int DISRemnantOpt_;
725 
729  unsigned int PtDistribution_;
730 
734  unsigned int pomeronStructure_;
736 
740  Energy mg_;
741 
742 };
743 
744 
745 }
746 
747 #endif /* HERWIG_HwRemDecayer_H */
virtual void doinit()
Initialize this object after the setup phase before saving an EventGenerator to disk.
Definition: HwRemDecayer.h:215
ThePEG::Ptr< Particle >::pointer PPtr
Energy _kinCutoff
The kinematic cut-off.
Definition: HwRemDecayer.h:684
vector< pair< tPPtr, tPPtr > > PartnerMap
Typedef to store information about colour partners.
Definition: HwRemDecayer.h:61
const pair< HadronContent, HadronContent > & content() const
Return the hadron content objects for the incoming particles.
Definition: HwRemDecayer.h:326
bool multiPeriph_
Switch to control using multiperipheral kinemaics.
Definition: HwRemDecayer.h:603
double sqrt(int x)
tcPDFPtr _pdf
The PDF for the current initial-state shower.
Definition: HwRemDecayer.h:591
int extracted
The array index of the extracted particle.
Definition: HwRemDecayer.h:311
pair< double, double > theX
variable to sum up the x values of the extracted particles
Definition: HwRemDecayer.h:551
unsigned int DISRemnantOpt_
Option for the DIS remnant.
Definition: HwRemDecayer.h:724
int sign
-1 if the particle is an anti-particle.
Definition: HwRemDecayer.h:314
StepPtr thestep
Variable to hold a pointer to the current step.
Definition: HwRemDecayer.h:570
HwRemDecayer()
The default constructor.
Definition: HwRemDecayer.h:68
ThePEG::Ptr< PDFBase >::transient_const_pointer tcPDFPtr
void extract(int id)
manually extract the valence flavour id.
Definition: HwRemDecayer.h:241
vector< int > flav
The valence flavours of the corresponding baryon.
Definition: HwRemDecayer.h:308
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
Energy _forcedSplitScale
The PDF freezing scale as set in ShowerHandler.
Definition: HwRemDecayer.h:689
pair< PartnerMap, PartnerMap > theMaps
Pair of PartnerMap&#39;s to store the particles, which will be colour connected in the end...
Definition: HwRemDecayer.h:563
bool isValenceQuark(int id) const
Method to determine whether parton is a valence quark.
Definition: HwRemDecayer.h:303
unsigned int pomeronStructure
Pomeron treatment.
Definition: HwRemDecayer.h:320
pair< HadronContent, HadronContent > theContent
Pair of HadronContent structs to know about the quark content of the beams.
Definition: HwRemDecayer.h:554
InvEnergy2 beta_
slope of the soft pt-spectrum: dN/dp_T = N p_T exp(-beta*p_T^2)
Definition: HwRemDecayer.h:626
void doSoftInteractions(unsigned int N)
Create N soft gluon interactions.
Definition: HwRemDecayer.h:441
virtual bool accept(const DecayMode &) const
Check if this decayer can perfom the decay specified by the given decay mode.
Definition: HwRemDecayer.h:93
unsigned int pomeronStructure_
Option for the treatment of the pomeron structure.
Definition: HwRemDecayer.h:734
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
Definition: HwRemDecayer.h:203
pair< tPPtr, tPPtr > tPPair
pair< bool, bool > theanti
A flag which indicates, whether the extracted valence quark was a anti particle.
Definition: HwRemDecayer.h:546
Iterator find(IteratorRange< Iterator > r, const T &t)
virtual bool multiCapable() const
Return true if this decayed can extract more than one parton from a particle.
Definition: HwRemDecayer.h:107
virtual IBPtr clone() const
Make a simple clone of this object.
Definition: HwRemDecayer.h:197
bool isValenceQuarkData(tcPDPtr partonData) const
Method to determine whether parton is a valence quark.
Definition: HwRemDecayer.h:295
void setMass(Energy a)
bool quarkPair_
True if kinematics is to be calculated for quarks.
Definition: HwRemDecayer.h:608
PPtr finalSplit(const tRemPPtr rem, long remID, Lorentz5Momentum usedMomentum) const
This creates a parton from the remaining flavours of the hadron.
Definition: HwRemDecayer.h:386
Simple struct to store info about baryon quark and di-quark constituents.
Definition: HwRemDecayer.h:236
tcPDPtr hadron
The ParticleData objects of the hadron.
Definition: HwRemDecayer.h:317
bool isSeaQuark(tcPPtr parton) const
Method to determine whether parton is a quark from the sea.
Definition: HwRemDecayer.h:273
double ladderbFactor_
Variable to store the additive factor of the multiperipheral ladder multiplicity. ...
Definition: HwRemDecayer.h:644
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
double gaussWidth_
Variable to store the gaussian width of the fluctuation of the longitudinal momentum fraction...
Definition: HwRemDecayer.h:662
bool allowTop_
Switch to control handling of top quarks in proton.
Definition: HwRemDecayer.h:598
double ladderNorm_
Variable of the parameterization of the ladder multiplicity.
Definition: HwRemDecayer.h:654
ThePEG::Ptr< Particle >::transient_pointer tPPtr
T gaussDistribution(T mean, T stdev) const
Methods to generate random distributions also all stolen form UA5Handler.
Definition: HwRemDecayer.h:475
double initTotRap_
Variable to store the initial total rapidity between of the remnants.
Definition: HwRemDecayer.h:674
ThePEG::Ptr< Particle >::transient_const_pointer tcPPtr
Ptr< BeamParticleData >::const_pointer theBeamData
the beam data
Definition: HwRemDecayer.h:586
tPPair softRems_
Pair of soft Remnant pointers, i.e.
Definition: HwRemDecayer.h:616
struct that is used to catch exceptions which are thrown due to energy conservation issues of additio...
Definition: HwRemDecayer.h:127
double colourDisrupt_
Variable to store the relative number of colour disrupted connections to additional soft subprocesses...
Definition: HwRemDecayer.h:638
Energy ptmin_
ptcut of the UE model
Definition: HwRemDecayer.h:621
tcPPtr theBeam
The beam particle data for the current incoming hadron.
Definition: HwRemDecayer.h:581
double _zbin
Size of the bins in z for the interpolation.
Definition: HwRemDecayer.h:699
bool isSeaQuarkData(tcPDPtr partonData) const
Method to determine whether parton is a quark from the sea.
Definition: HwRemDecayer.h:288
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:492
-*- C++ -*-
virtual void doinit()
int _nbinmax
Maximum number of bins for the z interpolation.
Definition: HwRemDecayer.h:709
unsigned int PtDistribution_
Option for the pT generation.
Definition: HwRemDecayer.h:729
pair< Lorentz5Momentum, Lorentz5Momentum > theUsed
Pair of Lorentz5Momentum to keep track of the forced splitting product momenta.
Definition: HwRemDecayer.h:557
double sign(double x)
Small helper.
Definition: RandomHelpers.h:30
double ladderPower_
Variable of the parameterization of the ladder multiplicity.
Definition: HwRemDecayer.h:649
unsigned int maxtrySoft_
Maximum number of attempts for the regeneration of an additional soft scattering, before the number o...
Definition: HwRemDecayer.h:632
Energy mg_
The gluon constituent mass.
Definition: HwRemDecayer.h:740
void setHadronContent(tPPair beam)
Set the hadron contents.
Definition: HwRemDecayer.h:338
The HwRemDecayer class is responsible for the decay of the remnants.
Definition: HwRemDecayer.h:56
double _ybin
Size of the bins in y for the interpolation.
Definition: HwRemDecayer.h:704
double valOfN_
Variable to store the current total multiplicity of a ladder.
Definition: HwRemDecayer.h:668
vector< PPtr > ParticleVector
bool isValenceQuark(tcPPtr parton) const
Method to determine whether parton is a valence quark.
Definition: HwRemDecayer.h:280
ThePEG::Ptr< Step >::pointer StepPtr
constexpr ZeroUnit ZERO
pair< RemPPtr, RemPPtr > theRems
Pair of Remnant pointers.
Definition: HwRemDecayer.h:576
ShowerAlphaPtr _alphaEM
Pointer to the object calculating the QED coupling.
Definition: HwRemDecayer.h:719
double _range
Range for emission.
Definition: HwRemDecayer.h:694
ShowerAlphaPtr _alphaS
Pointer to the object calculating the QCD coupling.
Definition: HwRemDecayer.h:714