herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
HadronSelector.h
1 // -*- C++ -*-
2 //
3 // HadronSelector.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_HadronSelector_H
10 #define HERWIG_HadronSelector_H
11 //
12 // This is the declaration of the HadronSelector class.
13 //
14 
15 #include "ThePEG/Interface/Interfaced.h"
16 #include <ThePEG/Persistency/PersistentOStream.h>
17 #include <ThePEG/Persistency/PersistentIStream.h>
18 #include <ThePEG/PDT/ParticleData.h>
20 #include <ThePEG/Repository/EventGenerator.h>
21 #include "HadronSelector.fh"
22 
23 namespace Herwig {
24 
25 using namespace ThePEG;
26 
48 class HadronSelector: public Interfaced {
49 
50 public:
51 
73  class HadronInfo {
74 
75  public:
76 
84  HadronInfo(long idin=0, tPDPtr datain=tPDPtr(),
85  double swtin=1., Energy massin=ZERO)
86  : id(idin), ptrData(datain), swtef(swtin), wt(1.0), overallWeight(0.0),
87  mass(massin)
88  {}
89 
93  bool operator<(const HadronInfo &x) const {
94  if(mass!=x.mass) return mass < x.mass;
95  else return id < x.id;
96  }
97 
101  long id;
102 
107 
111  double swtef;
112 
116  double wt;
117 
122 
126  Energy mass;
127 
131  void rescale(double x) const {
132  const_cast<HadronInfo*>(this)->overallWeight *= x;
133  }
134 
139  const HadronInfo & hi ) {
140  os << hi.id << hi.ptrData << hi.swtef << hi.wt
141  << hi.overallWeight << ounit(hi.mass,GeV);
142  return os;
143  }
144 
148  friend ostream & operator<< (ostream & os, const HadronInfo & hi );
149 
154  HadronInfo & hi ) {
155  is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt
156  >> hi.overallWeight >> iunit(hi.mass,GeV);
157  return is;
158  }
159  };
160 
169  class Kupco {
170 
171  public:
172 
180  Kupco(tcPDPtr inidQ,tcPDPtr inhad1,tcPDPtr inhad2, Energy inwgt)
181  : idQ(inidQ),hadron1(inhad1),hadron2(inhad2),weight(inwgt)
182  {}
183 
188 
193 
198 
202  Energy weight;
203  };
204 
205 public:
206 
214  typedef set<HadronInfo> KupcoData;
216 
220  typedef map<pair<long,long>,KupcoData> HadronTable;
221 
222 public:
223 
227  HadronSelector(unsigned int);
228 
237  virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass, tcPDPtr par1,
238  tcPDPtr par2,tcPDPtr par3 = PDPtr()) const
239  = 0;
240 
248  tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const;
249 
279  pair<tcPDPtr,tcPDPtr> lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
280  tcPDPtr ptr3 = PDPtr ()) const;
281 
289  tcPDPtr ptr3 = PDPtr ()) const {
290  pair<tcPDPtr,tcPDPtr> pairData = lightestHadronPair(ptr1, ptr2, ptr3);
291  if ( ! pairData.first || ! pairData.second ) return ZERO;
292  return ( pairData.first->mass() + pairData.second->mass() );
293  }
294 
305  tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
306  tcPDPtr ptr3 = PDPtr ()) const;
307 
320  vector<pair<tcPDPtr,double> >
321  hadronsBelowThreshold(Energy threshold,
322  tcPDPtr ptr1, tcPDPtr ptr2,
323  tcPDPtr ptr3 = PDPtr ()) const;
324 
331  Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
332 #ifndef NDEBUG
333  tcPDPtr ptr3 = PDPtr ()) const {
334 #else
335  tcPDPtr = PDPtr ()) const {
336 #endif
337  // The method assumes ptr3 == empty
338  assert(!(ptr3));
339  // find entry in the table
340  pair<long,long> ids(abs(ptr1->id()),abs(ptr2->id()));
341  HadronTable::const_iterator tit=_table.find(ids);
342  // throw exception if flavours wrong
343  if(tit==_table.end()||tit->second.empty())
344  throw Exception() << "HadronSelector::massLightestHadron "
345  << "failed for particle" << ptr1->id() << " "
346  << ptr2->id()
348  // return the mass
349  return tit->second.begin()->mass;
350  }
351 
357  Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const;
358 
366  double pwtDquark() const {
367  return _pwtDquark;
368  }
369 
373  double pwtUquark() const {
374  return _pwtUquark;
375  }
376 
380  double pwtSquark() const {
381  return _pwtSquark;
382  }
383 
387  double pwtCquark() const {
388  return _pwtCquark;
389  }
390 
394  double pwtBquark() const {
395  return _pwtBquark;
396  }
397 
401  double pwtDIquark() const {
402  return _pwtDIquark;
403  }
405 
406 public:
407 
414  void persistentOutput(PersistentOStream & os) const;
415 
421  void persistentInput(PersistentIStream & is, int version);
423 
430  static void Init();
431 
432 protected:
433 
446  virtual void doinit();
448 
449 protected:
450 
489  virtual void constructHadronTable();
490 
494  const HadronTable & table() const {
495  return _table;
496  }
497 
501  const vector<PDPtr> & partons() const {
502  return _partons;
503  }
504 
508  double pwt(long pid) const {
509  map<long,double>::const_iterator it = _pwt.find(abs(pid));
510  assert( it != _pwt.end() );
511  return it->second;
512  }
513 
540  double probabilityMixing(const double angleMix,
541  const int order) const {
542  static double convert=Constants::pi/180.0;
543  if (order == 1)
544  return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
545  else if (order == 2)
546  return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
547  else
548  return 1.;
549  }
550 
555  virtual double mixingStateWeight(long id) const;
557 
562  double specialWeight(long id) const;
563 
569  int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const;
570 
571 private:
572 
577  HadronSelector & operator=(const HadronSelector &) = delete;
578 
579 private:
580 
584  vector<PDPtr> _partons;
585 
589  vector<PDPtr> _forbidden;
590 
598  double _pwtDquark;
599 
603  double _pwtUquark;
604 
608  double _pwtSquark;
609 
613  double _pwtCquark;
614 
618  double _pwtBquark;
619 
623  double _pwtDIquark;
624 
628  map<long,double> _pwt;
630 
638  double _etamix;
639 
643  double _phimix;
644 
648  double _h1mix;
649 
653  double _f0mix;
654 
658  double _f1mix;
659 
663  double _f2mix;
664 
668  double _eta2mix;
669 
673  double _omhmix;
674 
678  double _ph3mix;
679 
683  double _eta2Smix;
684 
688  double _phi2Smix;
690 
699  vector<double> _weight1S0;
700 
704  vector<double> _weight3S1;
705 
709  vector<double> _weight1P1;
710 
714  vector<double> _weight3P0;
715 
719  vector<double> _weight3P1;
720 
724  vector<double> _weight3P2;
725 
729  vector<double> _weight1D2;
730 
734  vector<double> _weight3D1;
735 
739  vector<double> _weight3D2;
740 
744  vector<double> _weight3D3;
746 
750  vector<vector<vector<double> > > _repwt;
751 
759  double _sngWt;
760 
764  double _decWt;
766 
770  HadronTable _table;
771 
779  enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4};
781 
785  unsigned int _topt;
786 
790  unsigned int _trial;
791 
800  double _limBottom;
801  double _limCharm;
802  double _limExotic;
804 
808  unsigned int belowThreshold_;
809 };
810 
811 
812 }
813 
814 #endif /* HERWIG_HadronSelector_H */
double pwtUquark() const
The up quark weight.
double _f2mix
The mixing angle.
void rescale(double x) const
Rescale the weight for a given hadron.
Kupco(tcPDPtr inidQ, tcPDPtr inhad1, tcPDPtr inhad2, Energy inwgt)
Constructor.
double sqrt(int x)
This class selects the hadron flavours of a cluster decay.
vector< double > _weight3P1
The weights for the multiplets.
double _pwtBquark
The probability of producting a bottom quark.
ThePEG::Ptr< ParticleData >::transient_pointer tPDPtr
double _pwtDquark
The weights for the different quarks and diquarks.
double _phi2Smix
The mixing angle.
bool operator<(const HadronInfo &x) const
Comparision operator on mass.
vector< double > _weight1P1
The weights for the multiplets.
constexpr auto sqr(const T &x) -> decltype(x *x)
Energy weight
Weight factor of this componation.
double pwtBquark() const
The bottom quark weight.
double pwtCquark() const
The charm quark weight.
double _sngWt
Singlet and Decuplet weights.
tcPDPtr hadron1
The ParticleData object for the first hadron produced.
vector< vector< vector< double > > > _repwt
The weights for the excited meson multiplets.
double _decWt
The decuplet weight.
map< pair< long, long >, KupcoData > HadronTable
The hadron table type.
vector< double > _weight1D2
The weights for the multiplets.
map< long, double > _pwt
Weights for quarks and diquarks.
double _pwtDIquark
The probability of producting a diquark.
tcPDPtr idQ
id of the quark drawn from the vacuum.
ostream & operator<<(ostream &, const DecayIntegrator &)
Output information on the DecayIntegrator for debugging purposes.
HadronTable _table
The table of hadron data.
double pwtDIquark() const
The diquark weight.
double _f1mix
The mixing angle.
double probabilityMixing(const double angleMix, const int order) const
Methods for the mixing of mesons.
MesonMultiplets
Enums so arrays can be statically allocated.
double _h1mix
The mixing angle.
tPDPtr ptrData
pointer to ParticleData, to get the spin, etc...
unsigned int _trial
Which particles to produce for debugging purposes.
vector< double > _weight3P0
The weights for the multiplets.
vector< double > _weight3D1
The weights for the multiplets.
double _ph3mix
The mixing angle.
PersistentIStream & operator>>(PersistentIStream &is, map< ShowerInteraction, T, Cmp, A > &m)
Read a map with ShowerInteraction as the key.
vector< double > _weight3D3
The weights for the multiplets.
OUnit< T, UT > ounit(const T &t, const UT &ut)
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
ThePEG::Ptr< ParticleData >::pointer PDPtr
HadronInfo(long idin=0, tPDPtr datain=tPDPtr(), double swtin=1., Energy massin=ZERO)
Constructor.
double _pwtCquark
The probability of producting a charm quark.
double pwt(long pid) const
Access the parton weights.
Energy mass
The hadrons mass.
Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3=PDPtr()) const
Returns the mass of the lightest pair of hadrons with the given particles.
vector< double > _weight3P2
The weights for the multiplets.
vector< double > _weight3D2
The weights for the multiplets.
double _f0mix
The mixing angle.
double _eta2Smix
The mixing angle.
double swtef
singlet/decuplet/orbital factor
double _etamix
The mixing angles for the mesons containing light quarks.
double _pwtUquark
The probability of producting an up quark.
double pwtDquark() const
Return the weights for the different quarks and diquarks.
tcPDPtr hadron2
The ParticleData object for the second hadron produced.
vector< PDPtr > _forbidden
The PDG codes of the hadrons which cannot be produced in the hadronization.
vector< double > _weight3S1
The weights for the multiplets.
double _phimix
The mixing angle.
double overallWeight
(2*J+1)*wt*swtef
-*- C++ -*-
unsigned int belowThreshold_
Option for the selection of hadrons below the pair threshold.
double _eta2mix
The mixing angle.
vector< PDPtr > _partons
The PDG codes of the constituent particles allowed.
double _omhmix
The mixing angle.
set< HadronInfo > KupcoData
The helper classes.
const vector< PDPtr > & partons() const
Access to the list of partons.
Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3=PDPtr()) const
Return the nominal mass of the hadron returned by lightestHadron()
double pwtSquark() const
The strange quark weight.
double _pwtSquark
The probability of producting a strange quark.
Class used to store all the hadron information for easy access.
Class designed to make STL routines easy to use.
constexpr ZeroUnit ZERO
unsigned int _topt
Option for the construction of the tables.
vector< double > _weight1S0
The weights for the various meson multiplets to be used to supress the production of particular state...
IUnit< T, UT > iunit(T &t, const UT &ut)
const HadronTable & table() const
Access to the table of hadrons.