PYTHIA  8.316
FragmentationFlavZpT.h
1 // FragmentationFlavZpT.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2025 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 // The FlavContainer class is a simple container for flavour,
26 // including the extra properties needed for popcorn baryon handling.
27 // id = current flavour.
28 // rank = current rank; 0 for endpoint flavour and then increase by 1.
29 // nPop = number of popcorn mesons yet to be produced (1 or 0).
30 // idPop = (absolute sign of) popcorn quark, shared between B and Bbar.
31 // idVtx = (absolute sign of) vertex (= non-shared) quark in diquark.
32 
34 
35 public:
36 
37  // Constructor.
38  FlavContainer(int idIn = 0, int rankIn = 0, int nPopIn = 0,
39  int idPopIn = 0, int idVtxIn = 0) : id(idIn), rank(rankIn),
40  nPop(nPopIn), idPop(idPopIn), idVtx(idVtxIn) {}
41 
42  // Copy constructor.
44  id = flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
45  idVtx = flav.idVtx;}
46 
47  // Overloaded equal operator.
48  FlavContainer& operator=(const FlavContainer& flav) { if (this != &flav) {
49  id = flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
50  idVtx = flav.idVtx; } return *this; }
51 
52  // Invert flavour.
53  FlavContainer& anti() {id = -id; return *this;}
54 
55  // Read in a container into another, without/with id sign flip.
56  FlavContainer& copy(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  FlavContainer& anti(const FlavContainer& flav) { if (this != &flav) {
60  id = -flav.id; rank = flav.rank; nPop = flav.nPop; idPop = flav.idPop;
61  idVtx = flav.idVtx; } return *this; }
62 
63  // Check whether is diquark.
64  bool isDiquark() {int idAbs = abs(id);
65  return (idAbs > 1000 && idAbs < 10000 && (idAbs/10)%10 == 0);}
66 
67  // Stored properties.
68  int id, rank, nPop, idPop, idVtx;
69 
70 };
71 
72 //==========================================================================
73 
74 // The StringFlav class is used to select quark and hadron flavours.
75 
76 class StringFlav : public PhysicsBase {
77 
78 public:
79 
80  // Constructor.
81  StringFlav() : suppressLeadingB(), probQQtoQ(), probStoUD(), probSQtoQQ(),
82  probQQ1toQQ0(), probQandQQ(), probQandS(), probQandSinQQ(), probQQ1corr(),
83  probQQ1corrInv(), probQQ1norm(), probQQ1join(), mesonRate(),
84  mesonRateSum(), mesonMix1(), mesonMix2(), etaSup(), etaPrimeSup(),
85  decupletSup(), baryonCGSum(), baryonCGMax(), popcornRate(), popcornSpair(),
86  popcornSmeson(), barCGMax(), scbBM(), popFrac(), popS(), dWT(),
87  lightLeadingBSup(), heavyLeadingBSup(), probStoUDSav(), probQQtoQSav(),
88  probSQtoQQSav(), probQQ1toQQ0Sav(), alphaQQSav(), closePacking(),
89  doEnhanceDiquark(), enhanceStrange(), enhancePT(), enhanceDiquark(),
90  exponentMPI(), exponentNSP() {}
91 
92  // Destructor.
93  virtual ~StringFlav() {}
94 
95  // Initialize data members.
96  virtual void init();
97 
98  // Initialise parameters when using close packing.
99  virtual void init(double kappaModifier, double strangeJunc,
100  double probQQmod);
101 
102  // Pick a light d, u or s quark according to fixed ratios.
103  int pickLightQ() { double rndmFlav = probQandS * rndmPtr->flat();
104  if (rndmFlav < 1.) return 1;
105  if (rndmFlav < 2.) return 2;
106  return 3; }
107 
108  // Pick a new flavour (including diquarks) given an incoming one.
109  // Optional arguments: pT, kappaModifier, allowPop.
110  virtual FlavContainer pick(FlavContainer& flavOld,
111  double = -1., double = -1., bool allowPop = true);
112 
113  // Combine two flavours (including diquarks) to produce a hadron.
114  virtual int combine(FlavContainer& flav1, FlavContainer& flav2);
115 
116  // Ditto, simplified input argument for simple configurations.
117  virtual int combineId( int id1, int id2, bool keepTrying = true) {
118  FlavContainer flag1(id1); FlavContainer flag2(id2);
119  for (int i = 0; i < 100; ++i) { int idNew = combine( flag1, flag2);
120  if (idNew != 0 || !keepTrying) return idNew;} return 0;}
121 
122  // Combine three (di-) quark flavours into two hadrons.
123  virtual pair<int,int> combineDiquarkJunction(int id1, int id2, int id3);
124 
125  // Combine two flavours to produce a hadron with lowest possible mass.
126  virtual int combineToLightest( int id1, int id2);
127 
128  // Lightest flavour-neutral meson.
129  virtual int idLightestNeutralMeson() { return 111; }
130 
131  // General function, decides whether to just return the hadron id
132  // if thermal model was use or whether to combine the two flavours.
133  virtual int getHadronID(FlavContainer& flav1, FlavContainer& flav2,
134  double = -1.0, double = -1.0, bool = false) {
135  return combine(flav1, flav2); }
136 
137  // Return hadron mass.
138  virtual double getHadronMassWin(int idHad) {
139  return particleDataPtr->mSel(idHad); }
140 
141  // Assign popcorn quark inside an original (= rank 0) diquark.
142  void assignPopQ(FlavContainer& flav);
143 
144  // Combine two quarks to produce a diquark.
145  int makeDiquark(int id1, int id2, int idHad = 0);
146 
147  // Check if quark-diquark combination should be added. If so add.
148  void addQuarkDiquark(vector< pair<int,int> >& quarkCombis,
149  int qID, int diqID, int hadronID) {
150  bool allowed = true;
151  for (int iCombi = 0; iCombi < int(quarkCombis.size()); iCombi++)
152  if ( (qID == quarkCombis[iCombi].first ) &&
153  (diqID == quarkCombis[iCombi].second) ) allowed = false;
154  if (allowed) quarkCombis.push_back( (hadronID > 0) ?
155  make_pair( qID, diqID) : make_pair(-qID, -diqID) ); }
156 
157  // Get spin counter for mesons.
158  int getMesonSpinCounter(int hadronID) { hadronID = abs(hadronID);
159  int j = (hadronID % 10);
160  if (hadronID < 1000) return ((j==1) ? 0 : ( (j==3) ? 1 : 5 ));
161  if (hadronID < 20000) return ((j==1) ? 3 : 2);
162  if (hadronID > 20000) return 4;
163  return -1; }
164 
165  // Get the flavour and spin ratios calculated from the diquark weights.
166  // i: (0) q -> B B, (1) q -> B M B, (2) qq -> M B
167  // j: (0) s/u popcorn ratio, (1/2) s/u ratio for vertex quark if popcorn
168  // quark is u/d or s, (3) q/q' vertex quark ratio if popcorn quark is
169  // light and = q, (4/5/6) (spin 1)/(spin 0) ratio for su, us and ud
170  double getFlavourSpinRatios(int i, int j) {
171  return (i < 3 && j < 7) ? dWT[i][j] : -1.0;}
172 
173 protected:
174 
175  // Initialise derived parameters.
176  virtual void initDerived();
177 
178  // Constants: could only be changed in the code itself.
179  static const int mesonMultipletCode[6];
180  static const double baryonCGOct[6], baryonCGDec[6];
181 
182  // Settings for default Gaussian model.
184  double probQQtoQ, probStoUD, probSQtoQQ, probQQ1toQQ0, probQandQQ,
185  probQandS, probQandSinQQ, probQQ1corr, probQQ1corrInv, probQQ1norm,
186  probQQ1join[4], mesonRate[4][6], mesonRateSum[4], mesonMix1[2][6],
187  mesonMix2[2][6], etaSup, etaPrimeSup, decupletSup, baryonCGSum[6],
188  baryonCGMax[6], popcornRate, popcornSpair, popcornSmeson, barCGMax[8],
189  scbBM[3], popFrac, popS[3], dWT[3][7], lightLeadingBSup,
190  heavyLeadingBSup;
191  bool qqKappa;
192  double probStoUDSav, probQQtoQSav, probSQtoQQSav, probQQ1toQQ0Sav,
193  alphaQQSav;
194 
195  // Settings for closepacking.
196  bool closePacking, doEnhanceDiquark;
197  double enhanceStrange, enhancePT, enhanceDiquark, exponentMPI, exponentNSP;
198 
199  // Fragmentation weights container.
201 
202 };
203 
204 //==========================================================================
205 
206 // Functions for unnormalised, <z>, and RMSD(z) of Lund FF. The two latter
207 // return negative values in case of failure.
208 
209 double LundFFRaw(double z, double a, double b, double c, double mT2);
210 
211 double LundFFAvg(double a, double b, double mT2, double tol);
212 
213 double LundFFRms(double a, double b, double mT2, double tol);
214 
215 //==========================================================================
216 
217 // The StringZ class is used to sample the fragmentation function f(z).
218 
219 class StringZ : public PhysicsBase {
220 
221 public:
222 
223  // Constructor.
224  StringZ() : useNonStandC(), useNonStandB(), useNonStandH(), usePetersonC(),
225  usePetersonB(), usePetersonH(), useOldAExtra(), mc2(), mb2(),
226  aLund(), bLund(), aExtraSQuark(), aExtraDiquark(), rFactC(),
227  rFactB(), rFactH(), aNonC(), aNonB(), aNonH(), bNonC(), bNonB(),
228  bNonH(), epsilonC(), epsilonB(), epsilonH(), stopM(), stopNF(),
229  stopS() {}
230 
231  // Destructor.
232  virtual ~StringZ() {}
233 
234  // Initialize data members.
235  virtual bool init();
236 
237  // Fragmentation function: top-level to determine parameters.
238  virtual double zFrag( int idOld, int idNew = 0, double mT2 = 1.);
239 
240  // Fragmentation function: select z according to provided parameters.
241  virtual double zLund( double a, double b, double c = 1.,
242  double head = 1., double bNow = 0., int idFrag = 0,
243  bool isOldSQuark = false, bool isNewSQuark = false,
244  bool isOldDiquark = false, bool isNewDiquark = false);
245  virtual double zPeterson( double epsilon);
246  virtual double zLundMax( double a, double b, double c = 1.);
247 
248  // Parameters for stopping in the middle; overloaded for Hidden Valley.
249  virtual double stopMass() {return stopM;}
250  virtual double stopNewFlav() {return stopNF;}
251  virtual double stopSmear() {return stopS;}
252 
253  // a and b fragmentation parameters needed in some operations.
254  virtual double aAreaLund() {return aLund;}
255  virtual double bAreaLund() {return bLund;}
256 
257  // Method to derive both a and b parameters (from <z> and RMSD(z)).
258  bool deriveABLund( bool derivaA = false, bool deriveAExtraDiquark = false,
259  bool deriveAExtraSQuark = false);
260  // Method to derive bLund from <z> (for fixed a and reference mT2).
261  double deriveBLund( double avgZ, double a, double mT2ref);
262 
263  protected:
264 
265  // Constants: could only be changed in the code itself.
266  static const double CFROMUNITY, AFROMZERO, AFROMC, EXPMAX;
267 
268  // Initialization data, to be read from Settings.
269  bool useNonStandC, useNonStandB, useNonStandH,
270  usePetersonC, usePetersonB, usePetersonH, useOldAExtra;
271  double mc2, mb2, aLund, bLund, aExtraSQuark, aExtraDiquark, rFactC,
272  rFactB, rFactH, aNonC, aNonB, aNonH, bNonC, bNonB, bNonH,
273  epsilonC, epsilonB, epsilonH, stopM, stopNF, stopS;
274 
275  // Fragmentation weights container.
277 
278 };
279 
280 //==========================================================================
281 
282 // The StringPT class is used to select select transverse momenta.
283 
284 class StringPT : public PhysicsBase {
285 
286 public:
287 
288  // Constructor.
289  StringPT() : sigmaQ(), enhancedFraction(), enhancedWidth(), sigma2Had(),
290  closePacking(), enhancePT(), exponentMPI(), exponentNSP() {}
291 
292  // Destructor.
293  virtual ~StringPT() {}
294 
295  // Initialize data members.
296  virtual void init();
297 
298  // Return px and py as a pair.
299  virtual pair<double, double> pxy(int idIn = 0, double kappaModifier = -1.0);
300 
301  // Gaussian suppression of given pT2; used in MiniStringFragmentation.
302  virtual double suppressPT2(double pT2) { return exp(-pT2/sigma2Had); }
303 
304 protected:
305 
306  // Constants: could only be changed in the code itself.
307  static const double SIGMAMIN;
308 
309  // Initialization data, to be read from Settings.
310  double sigmaQ, enhancedFraction, enhancedWidth, sigma2Had;
311  // Optional enhanced pT widths for strange and/or diquarks.
312  bool useWidthPre{false};
313  double widthPreStrange{1.}, widthPreQQ0{1.}, widthPreQQ1{1.};
314  // Special for closepacking.
316  double enhancePT, exponentMPI, exponentNSP;
317 
318 private:
319 
320  // Fragmentation weights container.
321  WeightsFragmentation* wgtsPtr{};
322 
323 };
324 
325 //==========================================================================
326 
327 } // end namespace Pythia8
328 
329 #endif // Pythia8_FragmentationFlavZpT_H
virtual int getHadronID(FlavContainer &flav1, FlavContainer &flav2, double=-1.0, double=-1.0, bool=false)
Definition: FragmentationFlavZpT.h:133
double LundFFRms(double a, double b, double mT2, double tol)
Definition: FragmentationFlavZpT.cc:738
FlavContainer & operator=(const FlavContainer &flav)
Overloaded equal operator.
Definition: FragmentationFlavZpT.h:48
virtual double aAreaLund()
a and b fragmentation parameters needed in some operations.
Definition: FragmentationFlavZpT.h:254
Definition: PhysicsBase.h:26
virtual double getHadronMassWin(int idHad)
Return hadron mass.
Definition: FragmentationFlavZpT.h:138
int getMesonSpinCounter(int hadronID)
Get spin counter for mesons.
Definition: FragmentationFlavZpT.h:158
static const double SIGMAMIN
Constants: could only be changed in the code itself.
Definition: FragmentationFlavZpT.h:307
The StringPT class is used to select select transverse momenta.
Definition: FragmentationFlavZpT.h:284
virtual ~StringPT()
Destructor.
Definition: FragmentationFlavZpT.h:293
double LundFFAvg(double a, double b, double mT2, double tol)
Definition: FragmentationFlavZpT.cc:704
bool isDiquark()
Check whether is diquark.
Definition: FragmentationFlavZpT.h:64
Wave4 epsilon(Wave4 w1, Wave4 w2, Wave4 w3)
Permutation operator.
Definition: HelicityBasics.cc:56
virtual ~StringZ()
Destructor.
Definition: FragmentationFlavZpT.h:232
virtual ~StringFlav()
Destructor.
Definition: FragmentationFlavZpT.h:93
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:148
virtual double suppressPT2(double pT2)
Gaussian suppression of given pT2; used in MiniStringFragmentation.
Definition: FragmentationFlavZpT.h:302
The StringZ class is used to sample the fragmentation function f(z).
Definition: FragmentationFlavZpT.h:219
FlavContainer & anti()
Invert flavour.
Definition: FragmentationFlavZpT.h:53
FlavContainer(int idIn=0, int rankIn=0, int nPopIn=0, int idPopIn=0, int idVtxIn=0)
Constructor.
Definition: FragmentationFlavZpT.h:38
int id
Stored properties.
Definition: FragmentationFlavZpT.h:68
double getFlavourSpinRatios(int i, int j)
Definition: FragmentationFlavZpT.h:170
double sigmaQ
Initialization data, to be read from Settings.
Definition: FragmentationFlavZpT.h:310
bool suppressLeadingB
Settings for default Gaussian model.
Definition: FragmentationFlavZpT.h:183
virtual int combineId(int id1, int id2, bool keepTrying=true)
Ditto, simplified input argument for simple configurations.
Definition: FragmentationFlavZpT.h:117
FlavContainer & copy(const FlavContainer &flav)
Read in a container into another, without/with id sign flip.
Definition: FragmentationFlavZpT.h:56
virtual int idLightestNeutralMeson()
Lightest flavour-neutral meson.
Definition: FragmentationFlavZpT.h:129
FlavContainer(const FlavContainer &flav)
Copy constructor.
Definition: FragmentationFlavZpT.h:43
double LundFFRaw(double z, double a, double b, double c, double mT2)
Functions for the Lund symmetric FF: unnormalised, average, and RMSD.
Definition: FragmentationFlavZpT.cc:689
StringZ()
Constructor.
Definition: FragmentationFlavZpT.h:224
virtual double stopMass()
Parameters for stopping in the middle; overloaded for Hidden Valley.
Definition: FragmentationFlavZpT.h:249
Definition: FragmentationFlavZpT.h:33
Definition: Weights.h:354
static const double EXPMAX
Do not take exponent of too large or small number.
Definition: FragmentationFlavZpT.h:266
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
bool closePacking
Special for closepacking.
Definition: FragmentationFlavZpT.h:315
The StringFlav class is used to select quark and hadron flavours.
Definition: FragmentationFlavZpT.h:76
bool closePacking
Settings for closepacking.
Definition: FragmentationFlavZpT.h:196
bool useNonStandC
Initialization data, to be read from Settings.
Definition: FragmentationFlavZpT.h:269
int pickLightQ()
Pick a light d, u or s quark according to fixed ratios.
Definition: FragmentationFlavZpT.h:103
StringPT()
Constructor.
Definition: FragmentationFlavZpT.h:289
StringFlav()
Constructor.
Definition: FragmentationFlavZpT.h:81