StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ropewalk.h
1 // Ropewalk.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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 // Header file for Rope Hadronization. The Ropewalk takes care of setting
7 // up string geometry and calculating overlaps etc. The RopeDipole classes
8 // take care of the dynamics of string shoving. Flavour-composition-changing
9 // ropes are handled by FlavourRope which changes parameters, and RopeFragPars
10 // which calculates parameters for the rope.
11 //
12 // The file contains the following classes: RopeDipoleEnd,
13 // OverlappingRopeDipole, RopeDipole, Ropewalk, RopeFragPars and FlavourRope.
14 
15 #ifndef Pythia8_Ropewalk_H
16 #define Pythia8_Ropewalk_H
17 
18 #include "Pythia8/Basics.h"
19 #include "Pythia8/Event.h"
20 #include "Pythia8/FragmentationSystems.h"
21 #include "Pythia8/Info.h"
22 #include "Pythia8/ParticleData.h"
23 #include "Pythia8/Settings.h"
24 #include "Pythia8/PythiaStdlib.h"
25 #include "Pythia8/StringInteractions.h"
26 
27 namespace Pythia8 {
28 
29 //==================================================================
30 
31 // Define the end of a dipole, containing a pointer to the particle,
32 // and its index in the event record.
33 // Includes some methods for kinematics output.
34 
35 class RopeDipoleEnd {
36 
37 public:
38 
39  // Constructors sets event pointer and event record index.
40  RopeDipoleEnd() : e(NULL), ne(-1) { }
41  RopeDipoleEnd(Event* eIn, int neIn) : e(eIn), ne(neIn) { }
42 
43  // Get a pointer to the particle.
44  Particle* getParticlePtr() { if (!e) return NULL; return &(*e)[ne]; }
45 
46  // Get the particle index in event record.
47  int getNe() {return ne;}
48 
49  // Output methods for (modified) rapidity.
50  double labrap() { return getParticlePtr()->y(); }
51  double rap(double m0){ return getParticlePtr()->y(m0); }
52  double rap(double m0, RotBstMatrix& r) { return getParticlePtr()->y(m0, r); }
53 
54 private:
55 
56  // Pointer to the event and the particle index in event record.
57  Event* e;
58  int ne;
59 
60 };
61 
62 //==================================================================
63 
64 // A dipole has many dipoles overlapping with it. The OverlappingRopeDipole
65 // class does bookkeeping of this. Holds a pointer to the original dipole.
66 
67 // Forward declaration of RopeDipole class.
68 class RopeDipole;
69 
70 class OverlappingRopeDipole {
71 
72 public:
73 
74  // Constructor sets up coordinates in the rest frame of other dipole.
75  OverlappingRopeDipole(RopeDipole* d, double m0, RotBstMatrix& r);
76 
77  // Calculate the overlap at given y and b.
78  bool overlap(double y, Vec4 ba, double r0);
79 
80  // Has the dipole been hadronized?
81  bool hadronized();
82 
83 private:
84 
85  // Pointer to the original dipole.
86  RopeDipole* dipole;
87 
88 // Data members are made public rather than making
89 // the RopeDipole class a friend.
90 public:
91 
92  int dir;
93  double y1, y2;
94  Vec4 b1, b2;
95 
96 };
97 
98 //==================================================================
99 
100 // The RopeDipole class holds information about a colour dipole, as well as
101 // functionality to do shoving and to calculate effective string tension.
102 
103 class RopeDipole {
104 
105 public:
106 
107  // The RopeDipole constructor makes sure that d1 is always the colored
108  // end and d2 the anti-colored.
109  RopeDipole(RopeDipoleEnd d1In, RopeDipoleEnd d2In, int iSubIn,
110  Info* infoPtrIn);
111 
112  // Insert an excitation on dipole, if not already there.
113  void addExcitation(double ylab, Particle* ex);
114 
115  // Deep access to the RopeDipoleEnds (needed by OverlappingRopeDipole).
116  RopeDipoleEnd* d1Ptr() { return &d1; }
117  RopeDipoleEnd* d2Ptr() { return &d2; }
118 
119  // Get the rotation matrix to go to dipole rest frame.
120  RotBstMatrix getDipoleRestFrame();
121  RotBstMatrix getDipoleLabFrame();
122 
123  // Get the dipole momentum four-vector.
124  Vec4 dipoleMomentum();
125  // Get the spatial point interpolated to given rapidity.
126  // In the dipole rest frame.
127  Vec4 bInterpolateDip(double y, double m0);
128  // In the lab frame.
129  Vec4 bInterpolateLab(double y, double m0);
130  // Given a Lorentz matrix.
131  Vec4 bInterpolate(double y, RotBstMatrix rb, double m0);
132 
133  // Get the quantum numbers m,n characterizing all dipole overlaps
134  // at a given rapidity value.
135  pair<int, int> getOverlaps(double yfrac, double m0, double r0);
136 
137  // Add an overlapping dipole.
138  void addOverlappingDipole(OverlappingRopeDipole& d) {
139  overlaps.push_back(d); }
140 
141  // Get the maximal and minimal rapidity of the dipole.
142  double maxRapidity(double m0) { return (max(d1.rap(m0), d2.rap(m0))); }
143  double minRapidity(double m0) { return (min(d1.rap(m0), d2.rap(m0))); }
144 
145  // Get the maximal and minimal boosted rapidity of the dipole.
146  double maxRapidity(double m0, RotBstMatrix& r) { return (max(d1.rap(m0,r),
147  d2.rap(m0,r))); }
148  double minRapidity(double m0, RotBstMatrix& r) { return (min(d1.rap(m0,r),
149  d2.rap(m0,r))); }
150 
151  // Propagate the dipole itself.
152  void propagateInit(double deltat);
153 
154  // Propagate both dipole ends as well as all excitations.
155  void propagate(double deltat, double m0);
156 
157  // Redistribute momentum to two particles.
158  void splitMomentum(Vec4 mom, Particle* p1, Particle* p2, double frac = 0.5);
159 
160  // Put gluon excitations on the dipole.
161  void excitationsToString(double m0, Event& event);
162 
163  // Test if the dipole is hadronized.
164  bool hadronized() { return isHadronized; }
165 
166  // Get the (event colconfig) index.
167  int index() { return iSub; }
168 
169  // Recoil the dipole from adding a gluon. If the "dummy" option is set,
170  // the recoil will not be added, but only checked.
171  // Note: the gluon will not actually be added, only the recoil (if possible).
172  bool recoil(Vec4& pg, bool dummy = false);
173 
174  // Set dipole hadronized flag.
175  void hadronized(bool h) { isHadronized = h; }
176 
177  // The number of excitations on the dipole.
178  int nExcitations() { return int(excitations.size()); }
179 
180 private:
181 
182  // The ends (ie. particles) of the dipole.
183  RopeDipoleEnd d1, d2;
184 
185  // The propagated positions in the lab frame.
186  Vec4 b1, b2;
187 
188  // The string index (internal to the event).
189  int iSub;
190 
191  // Lorentz matrices to go to and from dipole rest frame.
192  RotBstMatrix rotFrom, rotTo;
193  bool hasRotFrom, hasRotTo;
194 
195  // The dipoles overlapping with this one.
196  vector<OverlappingRopeDipole> overlaps;
197 
198  // All excitations belonging to this dipole ordered in rapidity in lab frame.
199  map<double, Particle*> excitations;
200 
201  bool isHadronized;
202  Info* infoPtr;
203 
204 };
205 
206 //==================================================================
207 
208 // The Ropewalk class keeps track of all the strings making up ropes
209 // for shoving as well as flavour enhancement.
210 
211 class Ropewalk : public StringInteractions {
212 
213 public:
214 
215  // Constructor.
216  Ropewalk() : r0(), m0(), pTcut(),
217  shoveJunctionStrings(),
218  shoveMiniStrings(), shoveGluonLoops(), mStringMin(), limitMom(), rCutOff(),
219  gAmplitude(), gExponent(), deltay(), deltat(), tShove(), tInit(),
220  showerCut(), alwaysHighest() {}
221 
222  // The Ropewalk init function sets parameters and pointers.
223  virtual bool init();
224 
225  // Extract all dipoles from an event.
226  bool extractDipoles(Event& event, ColConfig& colConfig);
227 
228  // Calculate all overlaps of all dipoles and store as OverlappingRopeDipoles.
229  bool calculateOverlaps();
230 
231  // Calculate the effective string tension a fraction yfrac in on the dipole
232  // given by indices e1 and e2.
233  double getKappaHere(int e1, int e2, double yfrac);
234 
235  // The multiplicity of a colour state given its quantum numbers.
236  double multiplicity(double p, double q) {
237  return ( p < 0 || q < 0 || p + q == 0 )
238  ? 0.0 : 0.5 * (p + 1) * (q + 1) * (p + q + 2);
239  }
240 
241  // Calculate the average string tension of the event, in units of the default
242  // string tension (ie. 1 GeV/fm), using random walk in colour space.
243  double averageKappa();
244 
245  // Invoke the random walk and select a state.
246  pair<int, int> select(int m, int n, Rndm* rndm);
247 
248  // Shove all dipoles in the event.
249  void shoveTheDipoles(Event& event);
250 
251 private:
252 
253  // Parameters of the ropewalk.
254  double r0, m0, pTcut;
255  // Include junction strings in shoving.
256  bool shoveJunctionStrings;
257  // Include ministrings in shoving.
258  bool shoveMiniStrings;
259  // Include gluon loops in shoving.
260  bool shoveGluonLoops;
261  // Limit between string frag and ministring.
262  double mStringMin;
263  // Limit participating dipoles by their momentum.
264  bool limitMom;
265  // Radius cutoff in multiples of r0.
266  double rCutOff;
267  // Parameters of the shoving.
268  double gAmplitude, gExponent;
269  // Rapidity slicing.
270  double deltay;
271  // Time steps.
272  double deltat;
273  // Shove time.
274  double tShove;
275  // Intial propagation time.
276  double tInit;
277  // Final state shower cut-off.
278  double showerCut;
279  // Assume we are always in highest multiplet.
280  bool alwaysHighest;
281 
282  // All dipoles in the event sorted by event record.
283  // Index of the two partons.
284  typedef multimap<pair<int,int>, RopeDipole> DMap;
285  DMap dipoles;
286 
287  // All excitations.
288  vector< vector<Particle> > eParticles;
289 
290  // Random walk states and their weights.
291  vector<pair<int, int> > states;
292  vector<double> weights;
293 
294  // Private assignment operator.
295  Ropewalk& operator=(const Ropewalk) = delete;
296 
297 };
298 
299 //==================================================================
300 
301 // RopeFragPars recalculates fragmentation parameters according to a
302 // changed string tension. Helper class to FlavourRope.
303 
304 class RopeFragPars : public PhysicsBase {
305 
306 public:
307 
308  // Constructor.
309  RopeFragPars() : aIn(), adiqIn(), bIn(), rhoIn(), xIn(),
310  yIn(), xiIn(), sigmaIn(), kappaIn(), aEff(), adiqEff(), bEff(),
311  rhoEff(), xEff(), yEff(), xiEff(), sigmaEff(), kappaEff(),
312  beta() {}
313 
314  // The init function sets up initial parameters from settings.
315  bool init();
316 
317  // Return parameters at given string tension, ordered by their
318  // name for easy insertion in settings.
319  map<string,double> getEffectiveParameters(double h);
320 
321 private:
322 
323  // Constants: can only be changed in the code itself.
324  static const double DELTAA, ACONV, ZCUT;
325 
326  // Get the Fragmentation function a parameter from cache or calculate it.
327  double getEffectiveA(double thisb, double mT2, bool isDiquark);
328 
329  // Calculate the effective parameters.
330  bool calculateEffectiveParameters(double h);
331 
332  // Insert calculated parameters in cache for later (re-)use.
333  bool insertEffectiveParameters(double h);
334 
335  // Calculate the a parameter.
336  double aEffective(double aOrig, double thisb, double mT2);
337 
338  // The Lund fragmentation function.
339  double fragf(double z, double a, double b, double mT2);
340 
341  // Integral of the Lund fragmentation function with given parameter values.
342  double integrateFragFun(double a, double b, double mT2);
343 
344  // Helper function for integration.
345  double trapIntegrate(double a, double b, double mT2, double sOld, int n);
346 
347  // Parameter caches to re-use calculations. Sets of parameters, ordered in h.
348  map<double, map<string, double> > parameters;
349 
350  // Values of the a-parameter ordered in b*mT2 grid.
351  map<double, double> aMap;
352 
353  // The same for aDiquark.
354  map<double, double> aDiqMap;
355 
356  // Initial values of parameters.
357  double aIn, adiqIn, bIn, rhoIn, xIn, yIn, xiIn, sigmaIn, kappaIn;
358 
359  // Effective values of parameters.
360  double aEff, adiqEff, bEff, rhoEff, xEff, yEff, xiEff, sigmaEff, kappaEff;
361 
362  // Junction parameter.
363  double beta;
364 
365 };
366 
367 //==================================================================
368 
369 // The FlavourRope class takes care of placing a string breakup in
370 // the event, and assigning the string breakup effective parameters.
371 // It is a UserHooks derived class, and one must make sure to add it
372 // to the UserHooksVector in the main program or somewhere else.
373 
374 class FlavourRope : public FragmentationModifierBase {
375 
376 public:
377 
378  // Constructor.
379  FlavourRope(Ropewalk & rwIn) : rwPtr(&rwIn), ePtr(), doBuffon(),
380  rapiditySpan(), stringProtonRatio(), fixedKappa(), h() {}
381 
382  // Initialize. Set pointers.
383  virtual bool init() override {
384 
385  // Initialize event pointer such that it can be tested.
386  ePtr = NULL;
387  h = parm("Ropewalk:presetKappa");
388  fixedKappa = flag("Ropewalk:setFixedKappa");
389  doBuffon = flag("Ropewalk:doBuffon");
390  rapiditySpan = parm("Ropewalk:rapiditySpan");
391  stringProtonRatio = parm("Ropewalk:stringProtonRatio");
392  // Initialize FragPar.
393  fp.init();
394  return true;
395  }
396 
397  // Change the fragmentation parameters.
398  bool doChangeFragPar(StringFlav* flavPtr, StringZ* zPtr,
399  StringPT * pTPtr, double m2Had, vector<int> iParton, int endId) override;
400 
401  // Set enhancement manually.
402  void setEnhancement(double hIn) { h = hIn;}
403 
404  // Set pointer to the event.
405  void setEventPtr(Event& event) { ePtr = &event;}
406 
407  // Initialise the current event just before string fragmentation
408  // starts.
409  virtual bool initEvent(Event& event, ColConfig& colConfig) override;
410 
411 protected:
412 
413  virtual void onInitInfoPtr() override {
414  registerSubObject(fp);
415  }
416 
417 private:
418 
419  // Find breakup placement and fetch effective parameters.
420  // For model depending on vertex information.
421  map<string, double> fetchParameters(double m2Had, vector<int> iParton,
422  int endId);
423  // For simple Buffon model.
424  map<string, double> fetchParametersBuffon(double m2Had, vector<int> iParton,
425  int endId);
426 
427  // Pointer to the ropewalk object.
428  Ropewalk* rwPtr;
429 
430  // Pointer to the event.
431  Event* ePtr;
432 
433  // The object which handles change in parameters.
434  RopeFragPars fp;
435 
436  // Use Buffon prescription without vertex information.
437  bool doBuffon;
438 
439  // Parameters of Buffon prescription
440  double rapiditySpan, stringProtonRatio;
441 
442  // Keep track of hadronized dipoles in Buffon prescription.
443  vector<int> hadronized;
444 
445  // Use preset kappa from settings.
446  bool fixedKappa;
447 
448  // Locally stored string tension.
449  double h;
450 
451 };
452 
453 //==========================================================================
454 
455 // Interface to RopeWalk via an ShoverBase object.
456 
458 
459 public:
460 
461  // RopeWalk is a friend.
462  friend class RopeWalk;
463 
464  // The costructor needs a RopeWalk object be consistent.
465  RopewalkShover(Ropewalk & rwIn) : rwPtr(&rwIn) {}
466 
467  // Empty virtual destructor.
468  virtual ~RopewalkShover() {};
469  // Main virtual function, called before the hadronisation.
470  virtual bool stringRepulsion(Event & event, ColConfig & colConfig);
471 
472 private:
473 
474  // The Ropewalk object doing all the work.
475  Ropewalk * rwPtr;
476 
477 };
478 
479 } // end namespace Pythia8
480 
481 #endif // Pythia8_Ropewalk_H
Definition: rb.hh:21
Definition: AgUStep.h:26