PYTHIA  8.312
FragmentationFlavZpT.h
1 // FragmentationFlavZpT.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file contains helper classes for fragmentation.
7 // StringFlav is used to select quark and hadron flavours.
8 // StringPT is used to select transverse momenta.
9 // StringZ is used to sample the fragmentation function f(z).
10 
11 #ifndef Pythia8_FragmentationFlavZpT_H
12 #define Pythia8_FragmentationFlavZpT_H
13 
14 #include "Pythia8/Basics.h"
15 #include "Pythia8/MathTools.h"
16 #include "Pythia8/ParticleData.h"
17 #include "Pythia8/PhysicsBase.h"
18 #include "Pythia8/PythiaStdlib.h"
19 #include "Pythia8/Settings.h"
20 
21 namespace Pythia8 {
22 
23 //==========================================================================
24 
25 // Functions for unnormalised and average Lund FF.
26 
27 double LundFFRaw(double z, double a, double b, double c, double mT2);
28 
29 double LundFFAvg(double a, double b, double c, double mT2, double tol);
30 
31 //==========================================================================
32 
33 // The FlavContainer class is a simple container for flavour,
34 // including the extra properties needed for popcorn baryon handling.
35 // id = current flavour.
36 // rank = current rank; 0 for endpoint flavour and then increase by 1.
37 // nPop = number of popcorn mesons yet to be produced (1 or 0).
38 // idPop = (absolute sign of) popcorn quark, shared between B and Bbar.
39 // idVtx = (absolute sign of) vertex (= non-shared) quark in diquark.
40 
42 
43 public:
44 
45  // Constructor.
46  FlavContainer(int idIn = 0, int rankIn = 0, int nPopIn = 0,
47  int idPopIn = 0, int idVtxIn = 0) : id(idIn), rank(rankIn),
48  nPop(nPopIn), idPop(idPopIn), idVtx(idVtxIn) {}
49 
50  // Copy constructor.
52  id = flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
53  idVtx = flav.idVtx;}
54 
55  // Overloaded equal operator.
56  FlavContainer& operator=(const FlavContainer& flav) { if (this != &flav) {
57  id = flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
58  idVtx = flav.idVtx; } return *this; }
59 
60  // Invert flavour.
61  FlavContainer& anti() {id = -id; return *this;}
62 
63  // Read in a container into another, without/with id sign flip.
64  FlavContainer& copy(const FlavContainer& flav) { if (this != &flav) {
65  id = flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
66  idVtx = flav.idVtx; } return *this; }
67  FlavContainer& anti(const FlavContainer& flav) { if (this != &flav) {
68  id = -flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
69  idVtx = flav.idVtx; } return *this; }
70 
71  // Check whether is diquark.
72  bool isDiquark() {int idAbs = abs(id);
73  return (idAbs > 1000 && idAbs < 10000 && (idAbs/10)%10 == 0);}
74 
75  // Stored properties.
76  int id, rank, nPop, idPop, idVtx;
77 
78 };
79 
80 //==========================================================================
81 
82 // The StringFlav class is used to select quark and hadron flavours.
83 
84 class StringFlav : public PhysicsBase {
85 
86 public:
87 
88  // Constructor.
90  suppressLeadingB(),
91  mT2suppression(), useWidthPre(), probQQtoQ(), probStoUD(), probSQtoQQ(),
92  probQQ1toQQ0(), probQandQQ(), probQandS(), probQandSinQQ(), probQQ1corr(),
93  probQQ1corrInv(), probQQ1norm(), probQQ1join(), mesonRate(),
94  mesonRateSum(), mesonMix1(), mesonMix2(), etaSup(), etaPrimeSup(),
95  decupletSup(), baryonCGSum(), baryonCGMax(), popcornRate(), popcornSpair(),
96  popcornSmeson(), barCGMax(), scbBM(), popFrac(), popS(), dWT(),
97  lightLeadingBSup(), heavyLeadingBSup(), qqKappa(), closePackingFacPT2(),
98  closePackingFacQQ2(), probStoUDSav(), probQQtoQSav(), probSQtoQQSav(),
99  probQQ1toQQ0Sav(), alphaQQSav(), sigmaHad(), widthPreStrange(),
100  widthPreDiquark(), thermalModel(), mesonNonetL1(), temperature(),
101  tempPreFactor(), nNewQuark(), mesMixRate1(), mesMixRate2(), mesMixRate3(),
102  baryonOctWeight(), baryonDecWeight(), closePacking(), exponentMPI(),
103  exponentNSP(), hadronIDwin(0), idNewWin(0), hadronMassWin(-1.0) {}
104 
105  // Destructor.
106  virtual ~StringFlav() {}
107 
108  // Initialize data members.
109  virtual void init();
110 
111  // Initialise parameters when using close packing.
112  virtual void init(double kappaRatio, double strangeFac, double probQQmod);
113 
114  // Pick a light d, u or s quark according to fixed ratios.
115  int pickLightQ() { double rndmFlav = probQandS * rndmPtr->flat();
116  if (rndmFlav < 1.) return 1;
117  if (rndmFlav < 2.) return 2;
118  return 3; }
119 
120  // Pick a new flavour (including diquarks) given an incoming one,
121  // either by old standard Gaussian or new alternative exponential.
122  virtual FlavContainer pick(FlavContainer& flavOld, double pT = -1.0,
123  double kappaRatio = 0.0, bool allowPop = true) {
124  hadronIDwin = 0; idNewWin = 0; hadronMassWin = -1.0;
125  if ( (thermalModel || mT2suppression) && (pT >= 0.0) )
126  return pickThermal(flavOld, pT, kappaRatio);
127  return pickGauss(flavOld, allowPop); }
128  virtual FlavContainer pickGauss(FlavContainer& flavOld,
129  bool allowPop = true);
130  virtual FlavContainer pickThermal(FlavContainer& flavOld,
131  double pT, double kappaRatio);
132 
133  // Combine two flavours (including diquarks) to produce a hadron.
134  virtual int combine(FlavContainer& flav1, FlavContainer& flav2);
135 
136  // Ditto, simplified input argument for simple configurations.
137  virtual int combineId( int id1, int id2, bool keepTrying = true) {
138  FlavContainer flag1(id1); FlavContainer flag2(id2);
139  for (int i = 0; i < 100; ++i) { int idNew = combine( flag1, flag2);
140  if (idNew != 0 || !keepTrying) return idNew;} return 0;}
141 
142  // Combine three (di-) quark flavours into two hadrons.
143  virtual pair<int,int> combineDiquarkJunction(int id1, int id2, int id3);
144 
145  // Combine two flavours to produce a hadron with lowest possible mass.
146  virtual int combineToLightest( int id1, int id2);
147 
148  // Lightest flavour-neutral meson.
149  virtual int idLightestNeutralMeson() { return 111; }
150 
151  // Return chosen hadron in case of thermal model.
152  virtual int getHadronIDwin() { return hadronIDwin; }
153 
154  // Combine two flavours into hadron for last two remaining flavours
155  // for thermal model.
156  virtual int combineLastThermal(FlavContainer& flav1, FlavContainer& flav2,
157  double pT, double kappaRatio);
158 
159  // General function, decides whether to just return the hadron id
160  // if thermal model was use or whether to combine the two flavours.
161  virtual int getHadronID(FlavContainer& flav1, FlavContainer& flav2,
162  double pT = -1.0, double kappaRatio = 0, bool finalTwo = false) {
163  if (finalTwo) return ((thermalModel || mT2suppression) ?
164  combineLastThermal(flav1, flav2, pT, kappaRatio)
165  : combine(flav1, flav2));
166  if ((thermalModel || mT2suppression)&& (hadronIDwin != 0)
167  && (idNewWin != 0)) return getHadronIDwin();
168  return combine(flav1, flav2); }
169 
170  // Return hadron mass. Used one if present, pick otherwise.
171  virtual double getHadronMassWin(int idHad) { return
172  ((hadronMassWin < 0.0) ? particleDataPtr->mSel(idHad) : hadronMassWin); }
173 
174  // Assign popcorn quark inside an original (= rank 0) diquark.
175  void assignPopQ(FlavContainer& flav);
176 
177  // Combine two quarks to produce a diquark.
178  int makeDiquark(int id1, int id2, int idHad = 0);
179 
180  // Check if quark-diquark combination should be added. If so add.
181  void addQuarkDiquark(vector< pair<int,int> >& quarkCombis,
182  int qID, int diqID, int hadronID) {
183  bool allowed = true;
184  for (int iCombi = 0; iCombi < int(quarkCombis.size()); iCombi++)
185  if ( (qID == quarkCombis[iCombi].first ) &&
186  (diqID == quarkCombis[iCombi].second) ) allowed = false;
187  if (allowed) quarkCombis.push_back( (hadronID > 0) ?
188  make_pair( qID, diqID) : make_pair(-qID, -diqID) ); }
189 
190  // Get spin counter for mesons.
191  int getMesonSpinCounter(int hadronID) { hadronID = abs(hadronID);
192  int j = (hadronID % 10);
193  if (hadronID < 1000) return ((j==1) ? 0 : ( (j==3) ? 1 : 5 ));
194  if (hadronID < 20000) return ((j==1) ? 3 : 2);
195  if (hadronID > 20000) return 4;
196  return -1; }
197 
198  // Get the flavour and spin ratios calculated from the diquark weights.
199  // i: (0) q -> B B, (1) q -> B M B, (2) qq -> M B
200  // j: (0) s/u popcorn ratio, (1/2) s/u ratio for vertex quark if popcorn
201  // quark is u/d or s, (3) q/q' vertex quark ratio if popcorn quark is
202  // light and = q, (4/5/6) (spin 1)/(spin 0) ratio for su, us and ud
203  double getFlavourSpinRatios(int i, int j) {
204  return (i < 3 && j < 7) ? dWT[i][j] : -1.0;}
205 
206  // Calculate the flavor variations.
207  void variations(int idIn, bool early, bool noChoice);
208 
209 protected:
210 
211  // Initialise derived parameters.
212  virtual void initDerived();
213 
214  // Constants: could only be changed in the code itself.
215  static const int mesonMultipletCode[6];
216  static const double baryonCGOct[6], baryonCGDec[6];
217 
218  // Settings for Gaussian model.
219  bool suppressLeadingB, mT2suppression, useWidthPre;
220  double probQQtoQ, probStoUD, probSQtoQQ, probQQ1toQQ0, probQandQQ,
221  probQandS, probQandSinQQ, probQQ1corr, probQQ1corrInv, probQQ1norm,
222  probQQ1join[4], mesonRate[4][6], mesonRateSum[4], mesonMix1[2][6],
223  mesonMix2[2][6], etaSup, etaPrimeSup, decupletSup, baryonCGSum[6],
224  baryonCGMax[6], popcornRate, popcornSpair, popcornSmeson, barCGMax[8],
225  scbBM[3], popFrac, popS[3], dWT[3][7], lightLeadingBSup,
226  heavyLeadingBSup;
227  bool qqKappa;
228  double closePackingFacPT2, closePackingFacQQ2, probStoUDSav, probQQtoQSav,
229  probSQtoQQSav, probQQ1toQQ0Sav, alphaQQSav;
230  double sigmaHad, widthPreStrange, widthPreDiquark;
231 
232  // Settings for thermal model.
233  bool thermalModel, mesonNonetL1;
234  double temperature, tempPreFactor;
235  int nNewQuark;
236  double mesMixRate1[2][6], mesMixRate2[2][6], mesMixRate3[2][6];
237  double baryonOctWeight[6][6][6][2], baryonDecWeight[6][6][6][2];
238 
239  // Settings used by both models.
241  double exponentMPI, exponentNSP;
242 
243  // Key = hadron id, value = list of constituent ids.
244  map< int, vector< pair<int,int> > > hadronConstIDs;
245  // Key = initial (di)quark id, value = list of possible hadron ids
246  // + nr in hadronConstIDs.
247  map< int, vector< pair<int,int> > > possibleHadrons;
248  // Key = initial (di)quark id, value = prefactor to multiply rate.
249  map< int, vector<double> > possibleRatePrefacs;
250  // Similar, but for combining the last two (di)quarks. Key = (di)quark pair.
251  map< pair<int,int>, vector< pair<int,int> > > possibleHadronsLast;
252  map< pair<int,int>, vector<double> > possibleRatePrefacsLast;
253 
254  // Selection in thermal model.
255  int hadronIDwin, idNewWin;
256  double hadronMassWin;
257 
258 };
259 
260 //==========================================================================
261 
262 // The StringZ class is used to sample the fragmentation function f(z).
263 
264 class StringZ : public PhysicsBase {
265 
266 public:
267 
268  // Constructor.
269  StringZ() : useNonStandC(), useNonStandB(), useNonStandH(), usePetersonC(),
270  usePetersonB(), usePetersonH(), mc2(), mb2(), aLund(), bLund(),
271  aExtraSQuark(), aExtraDiquark(), rFactC(), rFactB(), rFactH(), aNonC(),
272  aNonB(), aNonH(), bNonC(), bNonB(), bNonH(), epsilonC(), epsilonB(),
273  epsilonH(), stopM(), stopNF(), stopS() {}
274 
275  // Destructor.
276  virtual ~StringZ() {}
277 
278  // Initialize data members.
279  virtual void init();
280 
281  // Fragmentation function: top-level to determine parameters.
282  virtual double zFrag( int idOld, int idNew = 0, double mT2 = 1.);
283 
284  // Fragmentation function: select z according to provided parameters.
285  virtual double zLund( double a, double b, double c = 1.,
286  double head = 1., double bNow = 0., int idFrag = 0,
287  bool isOldSQuark = false, bool isNewSQuark = false,
288  bool isOldDiquark = false, bool isNewDiquark = false);
289  virtual double zPeterson( double epsilon);
290  virtual double zLundMax( double a, double b, double c = 1.);
291 
292  // Parameters for stopping in the middle; overloaded for Hidden Valley.
293  virtual double stopMass() {return stopM;}
294  virtual double stopNewFlav() {return stopNF;}
295  virtual double stopSmear() {return stopS;}
296 
297  // a and b fragmentation parameters needed in some operations.
298  virtual double aAreaLund() {return aLund;}
299  virtual double bAreaLund() {return bLund;}
300 
301  // Method to derive bLund from <z> (for fixed a and reference mT2).
302  bool deriveBLund();
303 
304 protected:
305 
306  // Constants: could only be changed in the code itself.
307  static const double CFROMUNITY, AFROMZERO, AFROMC, EXPMAX;
308 
309  // Initialization data, to be read from Settings.
310  bool useNonStandC, useNonStandB, useNonStandH,
311  usePetersonC, usePetersonB, usePetersonH;
312  double mc2, mb2, aLund, bLund, aExtraSQuark, aExtraDiquark, rFactC,
313  rFactB, rFactH, aNonC, aNonB, aNonH, bNonC, bNonB, bNonH,
314  epsilonC, epsilonB, epsilonH, stopM, stopNF, stopS;
315 
316 };
317 
318 //==========================================================================
319 
320 // The StringPT class is used to select select transverse momenta.
321 
322 class StringPT : public PhysicsBase {
323 
324 public:
325 
326  // Constructor.
327  StringPT() : useWidthPre(), sigmaQ(), enhancedFraction(), enhancedWidth(),
328  sigma2Had(), widthPreStrange(), widthPreDiquark(), closePackingFacPT2(),
329  thermalModel(), temperature(), tempPreFactor(), fracSmallX(),
330  closePacking(), exponentMPI(), exponentNSP() {}
331 
332  // Destructor.
333  virtual ~StringPT() {}
334 
335  // Initialize data members.
336  virtual void init();
337 
338  // General function, return px and py as a pair in the same call
339  // in either model.
340  pair<double, double> pxy(int idIn, double kappaRatio = 0.0) {
341  return (thermalModel ? pxyThermal(idIn, kappaRatio) :
342  pxyGauss(idIn, kappaRatio)); }
343  pair<double, double> pxyGauss(int idIn = 0, double kappaRatio = 0.0);
344  pair<double, double> pxyThermal(int idIn, double kappaRatio = 0.0);
345 
346  // Gaussian suppression of given pT2; used in MiniStringFragmentation.
347  double suppressPT2(double pT2) { return (thermalModel ?
348  exp(-sqrt(pT2)/temperature) : exp(-pT2/sigma2Had)); }
349 
350 protected:
351 
352  // Constants: could only be changed in the code itself.
353  static const double SIGMAMIN;
354 
355  // Initialization data, to be read from Settings.
356  // Gaussian model.
358  double sigmaQ, enhancedFraction, enhancedWidth, sigma2Had,
359  widthPreStrange, widthPreDiquark, closePackingFacPT2;
360  // Thermal model.
362  double temperature, tempPreFactor, fracSmallX;
363  // Both.
365  double exponentMPI, exponentNSP;
366 
367 private:
368 
369  // Evaluate Bessel function K_{1/4}(x).
370  double BesselK14(double x);
371 
372 };
373 
374 //==========================================================================
375 
376 } // end namespace Pythia8
377 
378 #endif // Pythia8_FragmentationFlavZpT_H
int hadronIDwin
Selection in thermal model.
Definition: FragmentationFlavZpT.h:255
FlavContainer & operator=(const FlavContainer &flav)
Overloaded equal operator.
Definition: FragmentationFlavZpT.h:56
virtual double aAreaLund()
a and b fragmentation parameters needed in some operations.
Definition: FragmentationFlavZpT.h:298
Definition: PhysicsBase.h:27
virtual double getHadronMassWin(int idHad)
Return hadron mass. Used one if present, pick otherwise.
Definition: FragmentationFlavZpT.h:171
int getMesonSpinCounter(int hadronID)
Get spin counter for mesons.
Definition: FragmentationFlavZpT.h:191
map< int, vector< pair< int, int > > > possibleHadrons
Definition: FragmentationFlavZpT.h:247
static const double SIGMAMIN
Constants: could only be changed in the code itself.
Definition: FragmentationFlavZpT.h:353
bool thermalModel
Settings for thermal model.
Definition: FragmentationFlavZpT.h:233
The StringPT class is used to select select transverse momenta.
Definition: FragmentationFlavZpT.h:322
virtual ~StringPT()
Destructor.
Definition: FragmentationFlavZpT.h:333
virtual int getHadronID(FlavContainer &flav1, FlavContainer &flav2, double pT=-1.0, double kappaRatio=0, bool finalTwo=false)
Definition: FragmentationFlavZpT.h:161
bool thermalModel
Thermal model.
Definition: FragmentationFlavZpT.h:361
bool isDiquark()
Check whether is diquark.
Definition: FragmentationFlavZpT.h:72
Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3)
Permutation operator.
Definition: HelicityBasics.cc:56
virtual ~StringZ()
Destructor.
Definition: FragmentationFlavZpT.h:276
virtual ~StringFlav()
Destructor.
Definition: FragmentationFlavZpT.h:106
void addQuarkDiquark(vector< pair< int, int > > &quarkCombis, int qID, int diqID, int hadronID)
Check if quark-diquark combination should be added. If so add.
Definition: FragmentationFlavZpT.h:181
pair< double, double > pxy(int idIn, double kappaRatio=0.0)
Definition: FragmentationFlavZpT.h:340
The StringZ class is used to sample the fragmentation function f(z).
Definition: FragmentationFlavZpT.h:264
FlavContainer & anti()
Invert flavour.
Definition: FragmentationFlavZpT.h:61
FlavContainer(int idIn=0, int rankIn=0, int nPopIn=0, int idPopIn=0, int idVtxIn=0)
Constructor.
Definition: FragmentationFlavZpT.h:46
virtual FlavContainer pick(FlavContainer &flavOld, double pT=-1.0, double kappaRatio=0.0, bool allowPop=true)
Definition: FragmentationFlavZpT.h:122
int id
Stored properties.
Definition: FragmentationFlavZpT.h:76
double getFlavourSpinRatios(int i, int j)
Definition: FragmentationFlavZpT.h:203
bool suppressLeadingB
Settings for Gaussian model.
Definition: FragmentationFlavZpT.h:219
bool useWidthPre
Definition: FragmentationFlavZpT.h:357
virtual int combineId(int id1, int id2, bool keepTrying=true)
Ditto, simplified input argument for simple configurations.
Definition: FragmentationFlavZpT.h:137
double suppressPT2(double pT2)
Gaussian suppression of given pT2; used in MiniStringFragmentation.
Definition: FragmentationFlavZpT.h:347
FlavContainer & copy(const FlavContainer &flav)
Read in a container into another, without/with id sign flip.
Definition: FragmentationFlavZpT.h:64
virtual int idLightestNeutralMeson()
Lightest flavour-neutral meson.
Definition: FragmentationFlavZpT.h:149
FlavContainer(const FlavContainer &flav)
Copy constructor.
Definition: FragmentationFlavZpT.h:51
double LundFFRaw(double z, double a, double b, double c, double mT2)
Functions for unnormalised and average Lund FF.
Definition: FragmentationFlavZpT.cc:21
StringZ()
Constructor.
Definition: FragmentationFlavZpT.h:269
map< int, vector< pair< int, int > > > hadronConstIDs
Key = hadron id, value = list of constituent ids.
Definition: FragmentationFlavZpT.h:244
virtual double stopMass()
Parameters for stopping in the middle; overloaded for Hidden Valley.
Definition: FragmentationFlavZpT.h:293
Definition: FragmentationFlavZpT.h:41
map< pair< int, int >, vector< pair< int, int > > > possibleHadronsLast
Similar, but for combining the last two (di)quarks. Key = (di)quark pair.
Definition: FragmentationFlavZpT.h:251
static const double EXPMAX
Do not take exponent of too large or small number.
Definition: FragmentationFlavZpT.h:307
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
bool closePacking
Both.
Definition: FragmentationFlavZpT.h:364
The StringFlav class is used to select quark and hadron flavours.
Definition: FragmentationFlavZpT.h:84
double LundFFAvg(double a, double b, double c, double mT2, double tol)
Average, <z>, of Lund FF.
Definition: FragmentationFlavZpT.cc:30
bool closePacking
Settings used by both models.
Definition: FragmentationFlavZpT.h:240
bool useNonStandC
Initialization data, to be read from Settings.
Definition: FragmentationFlavZpT.h:310
int pickLightQ()
Pick a light d, u or s quark according to fixed ratios.
Definition: FragmentationFlavZpT.h:115
StringPT()
Constructor.
Definition: FragmentationFlavZpT.h:327
StringFlav()
Constructor.
Definition: FragmentationFlavZpT.h:89
map< int, vector< double > > possibleRatePrefacs
Key = initial (di)quark id, value = prefactor to multiply rate.
Definition: FragmentationFlavZpT.h:249
virtual int getHadronIDwin()
Return chosen hadron in case of thermal model.
Definition: FragmentationFlavZpT.h:152