StRoot  1
Ropewalk.cc
1 // Ropewalk.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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
7 // RopeRandState, RopeDipoleEnd, OverlappingRopeDipole, RopeDipole, Ropewalk,
8 // RopeFragPars and FlavourRope classes.
9
10 #include "Pythia8/Ropewalk.h"
11
12 namespace Pythia8 {
13
14 //==========================================================================
15
16 // OverlappingRopeDipole class.
17 // This class describes dipoles overlapping with a given dipole.
18
19 //--------------------------------------------------------------------------
20
21 // Constructor sets up coordinates in other dipole's rest frame.
22
23 OverlappingRopeDipole::OverlappingRopeDipole(RopeDipole* d, double m0,
24  RotBstMatrix& r) : dipole(d), dir(1) {
25
26  // Coordinates in other dipole's rest frame
27  b1 = d->d1Ptr()->getParticlePtr()->vProd();
28  b1.rotbst(r);
29  b2 = d->d2Ptr()->getParticlePtr()->vProd();
30  b2.rotbst(r);
31  y1 = d->d1Ptr()->rap(m0,r);
32  y2 = d->d2Ptr()->rap(m0,r);
33  if (y1 < y2) dir = -1;
34
35 }
36
37 //--------------------------------------------------------------------------
38
39 // Calculate the overlap at given y and b.
40
41 bool OverlappingRopeDipole::overlap(double y, Vec4 ba, double r0) {
42
43  if (y < min(y1, y2) || y > max(y1, y2)) return false;
44  Vec4 bb = b1 + (b2 - b1) * (y - y1) / (y2 - y1);
45  Vec4 tmp = ba - bb;
46  return (tmp.pT() <= 2 * r0);
47
48 }
49
50 //--------------------------------------------------------------------------
51
52 // Has the dipole been hadronized?
53
55
57
58 }
59
60 //==========================================================================
61
62 // RopeDipole class.
63 // This class describes a dipoles in impact parameter space.
64
65 //--------------------------------------------------------------------------
66
67 // The RopeDipole constructor makes sure that d1 is always the colored
68 // end and d2 the anti-colored.
69
70 RopeDipole::RopeDipole(RopeDipoleEnd d1In, RopeDipoleEnd d2In, int iSubIn,
71  Info* infoPtrIn)
72  : d1(d1In), d2(d2In), iSub(iSubIn), hasRotFrom(false), hasRotTo(false),
74
75  // Test if d1 is colored end and d2 anti-colored.
76  if (d1In.getParticlePtr()->col() == d2In.getParticlePtr()->acol()
77  && d1In.getParticlePtr()->col() != 0) { }
78  else { d2 = d1In, d1 = d2In; }
79
80 }
81
82 //--------------------------------------------------------------------------
83
84 // Insert an excitation on dipole, if not already there.
85
86 void RopeDipole::addExcitation(double ylab, Particle* ex) {
87
88  pair<map<double, Particle*>::iterator, map<double, Particle*>::iterator>
89  ret = excitations.equal_range(ylab);
90  for (map<double, Particle*>::iterator itr = ret.first; itr != ret.second;
91  ++itr)
92  if (ex == itr->second) return;
93  excitations.insert( make_pair(ylab,ex) );
94
95 }
96
97 //--------------------------------------------------------------------------
98
99 // Propagate the dipole itself.
100
101 void RopeDipole::propagateInit(double deltat) {
102
103  // Dipole end momenta.
104  Vec4 pcm = d1.getParticlePtr()->p();
105  Vec4 pam = d2.getParticlePtr()->p();
106  double mTc = sqrt(pcm.pT2() + pcm.m2Calc());
107  double mTa = sqrt(pam.pT2() + pam.m2Calc());
108  if (mTc == 0 || mTa == 0)
109  infoPtr->errorMsg("Error in RopeDipole::propagateInit: Tried to"
110  "propagate a RopeDipoleEnd with mT = 0");
111
112  // New vertices in the lab frame.
113  Vec4 newv1 = Vec4(d1.getParticlePtr()->xProd() + deltat * pcm.px() / mTc,
114  d1.getParticlePtr()->yProd() + deltat * pcm.py() / mTc, 0, 0);
115  Vec4 newv2 = Vec4(d2.getParticlePtr()->xProd() + deltat * pam.px() / mTa,
116  d2.getParticlePtr()->yProd() + deltat * pam.py() / mTa, 0, 0);
117  // Set the new vertices deep.
118  d1.getParticlePtr()->vProd(newv1);
119  d2.getParticlePtr()->vProd(newv2);
120
121 }
122
123 //--------------------------------------------------------------------------
124
125 // Propagate both dipole ends as well as all excitations.
126
127 void RopeDipole::propagate(double deltat, double m0) {
128
129  // First propagate the dipole ends.
130  propagateInit(deltat);
131  for (map<double, Particle*>::iterator eItr = excitations.begin();
132  eItr != excitations.end(); ++eItr) {
133
134  Vec4 em = eItr->second->p();
135  em.rotbst(getDipoleLabFrame());
136  // Propagate excitations.
137
138  if (em.pT() > 0.0){
139  Vec4 newVert = Vec4(eItr->second->xProd() + deltat * em.px() / em.pT(),
140  eItr->second->yProd() + deltat * em.py() / em.pT(), 0, 0);
141  eItr->second->vProd(newVert);
142  }
143  else eItr->second->vProd(bInterpolateLab(eItr->first,m0));
144  }
145
146 }
147
148 //--------------------------------------------------------------------------
149
150 // Put gluon excitations on the dipole.
151
152 void RopeDipole::excitationsToString(double m0, Event& event) {
153
154  // Erase excitations below cut-off.
155  map<double, Particle*>::iterator pItr = excitations.begin();
156  while (pItr != excitations.end() ) {
157  if (pItr->second->pAbs() < 1e-6) {
158  map<double, Particle*>::iterator eraseMe = pItr;
159  ++pItr;
160  excitations.erase(eraseMe);
161  }
162  else ++pItr;
163  }
164
165  // We now colour connect the excitations to the dipole.
166  // The dipole is connected like this sketch:
167  // acol (d1) col ------ acol (d2) col.
168  int oldcol = d1.getParticlePtr()->col();
169  if (oldcol != d2.getParticlePtr()->acol()) {
170  infoPtr->errorMsg("Error in Ropewalk::RopeDipole::excitationsToString: "
171  "color indices do not match.");
172  return;
173  }
174  vector<int> daughters;
175
176  // Two cases depending on which end we should start at.
177  // We always start at min rapidity and connect from there.
178  if (d1.rap(m0) == minRapidity(m0)) {
179  int acol = oldcol;
180  for (map<double, Particle*>::iterator itr = excitations.begin();
181  itr != excitations.end(); ++itr) {
182  int col = event.nextColTag();
183  itr->second->status(51);
184  itr->second->mothers(d1.getNe(),d1.getNe());
185  itr->second->cols(col,acol);
186  daughters.push_back(event.append(Particle(*(itr->second))));
187  acol = col;
188  }
189  d2.getParticlePtr()->acol(acol);
190  event[d2.getNe()].acol(acol);
191  }
192  else {
193  int acol = oldcol;
194  for (map<double, Particle*>::reverse_iterator itr = excitations.rbegin();
195  itr != excitations.rend(); ++itr) {
196  int col = event.nextColTag();
197  itr->second->status(51);
198  itr->second->mothers(d1.getNe(),d1.getNe());
199  itr->second->cols(col,acol);
200  daughters.push_back(event.append(Particle(*(itr->second))));
201  acol = col;
202  }
203  d2.getParticlePtr()->acol(acol);
204  event[d2.getNe()].acol(acol);
205  }
206  bool stringEnd = false;
207  if (d2.getParticlePtr()->col() == 0) stringEnd = true;
208
209  // Update status codes and mother/daughter indexing.
210  event[d1.getNe()].statusNeg();
211  Particle cc1 = *d1.getParticlePtr();
212  cc1.statusPos();
213  cc1.mothers(d1.getNe(),d1.getNe());
214  daughters.push_back(event.append(cc1));
215  event[d1.getNe()].daughters( daughters[0], daughters[daughters.size() -1 ] );
216  if (stringEnd) {
217  event[d2.getNe()].statusNeg();
218  Particle cc2 = *d2.getParticlePtr();
219  cc2.statusPos();
220  cc2.mothers(d2.getNe(),d2.getNe());
221  int did = event.append(cc2);
222  event[d2.getNe()].daughters(did,did);
223  }
224
225 }
226
227 //--------------------------------------------------------------------------
228
229 // Redistribute momentum to two particles.
230
231 void RopeDipole::splitMomentum(Vec4 mom, Particle* p1, Particle* p2,
232  double frac) {
233
234  Vec4 p1new = p1->p() + frac * mom;
235  Vec4 p2new = p2->p() + (1. - frac) * mom;
236  p1->p(p1new);
237  p2->p(p2new);
238  return;
239
240 }
241
242 //--------------------------------------------------------------------------
243
244 // Recoil the dipole from adding a gluon. If the "dummy" option is set,
245 // the recoil will not be added, but only checked.
246 // Note: the gluon will not actually be added, only the recoil (if possible).
247
248 bool RopeDipole::recoil(Vec4& pg, bool dummy) {
249
250  // Keep track of direction.
251  int sign = 1;
252  if (d1.rap(1.0) > d2.rap(1.0)) sign = -1;
253
254  // Lightcone momenta after inserting the gluon.
255  Particle* epaPtr = d1.getParticlePtr();
256  Particle* epcPtr = d2.getParticlePtr();
257  double pplus = epcPtr->pPos() + epaPtr->pPos() - pg.pPos();
258  double pminus = epcPtr->pNeg() + epaPtr->pNeg() - pg.pNeg();
259
260  // The new lightcone momenta of the dipole ends.
261  double ppa = 0.0;
262  double ppc = 0.0;
263  double pma = 0.0;
264  double pmc = 0.0;
265  double mta2 = epaPtr->mT2();
266  double mtc2 = epcPtr->mT2();
267  double mta = sqrt(mta2);
268  double mtc = sqrt(mtc2);
269  if ( pplus * pminus <= pow2(mta + mtc)
270  || pplus <= 0.0 || pminus <= 0.0 ) return false;
271
272  // Calculate the new momenta.
273  double sqarg = pow2(pplus * pminus - mta2 - mtc2) - 4. * mta2 * mtc2;
274  if (sqarg <= 0.0 ) return false;
275  if (sign > 0) {
276  ppa = 0.5 * (pplus * pminus + mta2 - mtc2 + sqrt(sqarg)) / pminus;
277  pma = mta2 / ppa;
278  pmc = pminus - pma;
279  ppc = mtc2 / pmc;
280  // Check rapidity after recoil.
281  if ( ppa * mtc < ppc * mta ) return false;
282  }
283  else {
284  pma = 0.5 * (pplus * pminus + mta2 - mtc2 + sqrt(sqarg)) / pplus;
285  ppa = mta2 / pma;
286  ppc = pplus - ppa;
287  pmc = mtc2 / ppc;
288  // Check rapidity after recoil.
289  if ( ppa*mtc > ppc*mta ) return false;
290  }
291
292  // Set up and store the new momenta.
293  Vec4 shifta = Vec4( epaPtr->px(), epaPtr->py(),
294  0.5 * (ppa - pma), 0.5 * (ppa + pma));
295  Vec4 shiftc = Vec4( epcPtr->px(), epcPtr->py(),
296  0.5 * (ppc - pmc), 0.5 * (ppc + pmc));
297  if (!dummy) {
298  epaPtr->p(shifta);
299  epcPtr->p(shiftc);
300  }
301  return true;
302
303 }
304
305 //--------------------------------------------------------------------------
306
307 // Get the Lorentz matrix to go to the dipole rest frame.
308
309 RotBstMatrix RopeDipole::getDipoleRestFrame() {
310
311  if (hasRotTo) return rotTo;
312
313  RotBstMatrix r;
314  r.toCMframe(d1.getParticlePtr()->p(),d2.getParticlePtr()->p());
315  rotTo = r;
316  hasRotTo = true;
317  return rotTo;
318
319 }
320
321 //--------------------------------------------------------------------------
322
323 // Get the Lorentz matrix to go from the dipole rest frame to lab frame.
324
325 RotBstMatrix RopeDipole::getDipoleLabFrame() {
326
327  if(hasRotFrom) return rotFrom;
328
329  RotBstMatrix r;
330  r.fromCMframe(d1.getParticlePtr()->p(),d2.getParticlePtr()->p());
331  rotFrom = r;
332  hasRotFrom = true;
333  return rotFrom;
334
335 }
336 //--------------------------------------------------------------------------
337
338 // The dipole four-momentum.
339
340 Vec4 RopeDipole::dipoleMomentum() {
341
342  Vec4 ret = d1.getParticlePtr()->p() + d2.getParticlePtr()->p();
343  return ret;
344
345 }
346
347 //--------------------------------------------------------------------------
348
349 // Interpolate (linear) between dipole ends to get b position at given y.
350 // Here y must be given in dipole rest frame, and the resulting b-position
351 // will also be in the dipole rest frame.
352
353 Vec4 RopeDipole::bInterpolateDip(double y, double m0) {
354  if(!hasRotTo) getDipoleRestFrame();
355  Vec4 bb1 = d1.getParticlePtr()->vProd();
356  bb1.rotbst(rotTo);
357  Vec4 bb2 = d2.getParticlePtr()->vProd();
358  bb2.rotbst(rotTo);
359  double y1 = d1.rap(m0,rotTo);
360  double y2 = d2.rap(m0,rotTo);
361  return bb1 + y * (bb2 - bb1) / (y2 - y1);
362
363 }
364
365 //--------------------------------------------------------------------------
366
367 // Interpolate (linear) between dipole ends to get b position at given y.
368
369 Vec4 RopeDipole::bInterpolateLab(double y, double m0) {
370
371  Vec4 bb1 = d1.getParticlePtr()->vProd();
372  Vec4 bb2 = d2.getParticlePtr()->vProd();
373  double y1 = d1.rap(m0);
374  double y2 = d2.rap(m0);
375  return bb1 + y * (bb2 - bb1) / (y2 - y1);
376
377 }
378
379 //--------------------------------------------------------------------------
380
381 // Interpolate (linear) between dipole ends to get b position in
382 // a given frame at given y.
383
384 Vec4 RopeDipole::bInterpolate(double y, RotBstMatrix rb, double m0) {
385
386  Vec4 bb1 = d1.getParticlePtr()->vProd();
387  Vec4 bb2 = d2.getParticlePtr()->vProd();
388  bb1.rotbst(rb);
389  bb2.rotbst(rb);
390  double y1 = d1.rap(m0);
391  double y2 = d2.rap(m0);
392  return bb1 + y * (bb2 - bb1) / (y2 - y1);
393
394 }
395
396 //--------------------------------------------------------------------------
397
398 // Calculate the amount of overlapping dipoles at a given rapidity.
399 // Return the number of overlapping dipoles to the "left" and "right".
400
401 pair<int, int> RopeDipole::getOverlaps(double yfrac, double m0, double r0) {
402  // Transform yfrac to y in dipole rest frame
403  if (!hasRotTo) getDipoleRestFrame();
404  double yL = d1.rap(m0,rotTo);
405  double yS = d2.rap(m0,rotTo);
406  double yH = yS + (yL - yS) * yfrac;
407  int m = 0, n = 0;
408  for (size_t i = 0; i < overlaps.size(); ++i) {
409  if (overlaps[i].overlap( yfrac, bInterpolateDip(yH,m0), r0)
411  if (overlaps[i].dir > 0) ++m;
412  else ++n;
413  }
414  }
415  return make_pair(m,n);
416
417 }
418
419 //==========================================================================
420
421 // Exc class.
422 // It is a helper class to Ropewalk, used to describe a pair of excitations
423 // needed for shoving. It is kept away from users, as there are a lot of
424 // raw pointers floating around.
425
426 //--------------------------------------------------------------------------
427
428 struct Exc {
429
430 // The constructor.
431 Exc(double yIn, double m0In, int iIn, int jIn, int kIn, RopeDipole* dip1In,
432  RopeDipole* dip2In) : y(yIn), m0(m0In), i(iIn), j(jIn), k(kIn), pp1(NULL),
433  pp2(NULL), dip1(dip1In), dip2(dip2In) { }
434
435 // Set pointers to the two excitation particles.
436 void setParticlePtrs(Particle* p1, Particle* p2) {
437  pp1 = p1;
438  pp2 = p2;
439  // Set a pointer to the excitation in the dipole.
442 }
443
444 // Give the excitation a kick in the x and y direction,
445 void shove(double dpx, double dpy) {
446  // The old momenta.
447  Vec4 p2 = pp2->p();
448  Vec4 p1 = pp1->p();
449  // The new momenta, after the shove.
450  double mt2new = sqrt(pow2(p2.px() - dpx) + pow2(p2.py() - dpy));
451  double e2new = mt2new * cosh(y);
452  double p2znew = mt2new * sinh(y);
453  Vec4 p2new(p2.px() - dpx, p2.py() - dpy, p2znew, e2new);
454  double mt1new = sqrt(pow2(p1.px() + dpx) + pow2(p1.py() + dpy));
455  double e1new = mt1new * cosh(y);
456  double p1znew = mt1new * sinh(y);
457  Vec4 p1new(p1.px() + dpx, p1.py() + dpy, p1znew, e1new);
458  // The differences between the two.
459  Vec4 deltap1 = p1new - p1;
460  Vec4 deltap2 = p2new - p2;
461  // Now check if we can add these two excitations to the dipoles.
462  if ( dip2->recoil(deltap2) ) {
463  if ( dip1->recoil(deltap1) ) {
464  pp1->p(p1new);
465  pp2->p(p2new);
466  } else {
467  Vec4 dp2 = -deltap2;
468  dip2->recoil(dp2);
469  }
470  }
471 }
472
473 // The push direction as a four vector.
474 Vec4 direction() {
475  return dip1->bInterpolateDip(y,m0) -
476  dip2->bInterpolateDip(y,m0);
477 }
478
479 // Member variables, slice rapidity and small cut-off mass.
480 double y;
481 double m0;
482
483 // Local particle indices.
484 int i, j, k;
485 Particle* pp1;
486 Particle* pp2;
487 RopeDipole* dip1;
488 RopeDipole* dip2;
489 };
490
491 //==========================================================================
492
493 // Ropewalk class.
494 // This class keeps track of all the strings making up ropes for shoving
495 // as well as flavour enhancement.
496
497 //--------------------------------------------------------------------------
498
499 // The Ropewalk init function sets parameters and pointers.
500
501 bool Ropewalk::init(Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn) {
502
503  // Save pointers.
504  infoPtr = infoPtrIn;
505  rndmPtr = rndmPtrIn;
506
507  // Parameters of the ropewalk.
508  doShoving = settings.flag("Ropewalk:doShoving");
509  shoveMiniStrings = settings.flag("Ropewalk:shoveMiniStrings");
510  shoveJunctionStrings = settings.flag("Ropewalk:shoveJunctionStrings");
511  shoveGluonLoops = settings.flag("Ropewalk:shoveGluonLoops");
512  limitMom = settings.flag("Ropewalk:limitMom");
514  r0 = settings.parm("Ropewalk:r0");
515  m0 = settings.parm("Ropewalk:m0");
516  pTcut = settings.parm("Ropewalk:pTcut");
517  rCutOff = settings.parm("Ropewalk:rCutOff");
518  gAmplitude = settings.parm("Ropewalk:gAmplitude");
519  gExponent = settings.parm("Ropewalk:gExponent");
520  deltay = settings.parm("Ropewalk:deltay");
521  deltat = settings.parm("Ropewalk:deltat");
522  tShove = settings.parm("Ropewalk:tShove");
523  tInit = settings.parm("Ropewalk:tInit");
524  showerCut = settings.parm("TimeShower:pTmin");
525  alwaysHighest = settings.flag("Ropewalk:alwaysHighest");
526
527  // Check consistency.
528  if (deltat > tShove) {
529  infoPtr->errorMsg("Error in Ropewalk::init: "
530  "deltat cannot be larger than tShove");
531  return false;
532  }
533  return true;
534
535 }
536
537 //--------------------------------------------------------------------------
538
539 // Calculate the average string tension of the event, in units of the default
540 // string tension (ie. 1 GeV/fm), using random walk in colour space.
541
542 double Ropewalk::averageKappa() {
543
544  double kap = 0.;
545  double nd = 0;
546  for (DMap::iterator itr = dipoles.begin(); itr != dipoles.end(); ++itr) {
547
548  // Getting the overlaps: m, n.
549  pair<int,int> overlap = itr->second.getOverlaps( rndmPtr->flat(), m0, r0);
550
551  // Overlaps define the number of steps taken in the random walk.
552  // We define the present dipole to always point in the p-direction.
553  pair<double, double> pq = select( overlap.first + 1, overlap.second,
554  rndmPtr);
555  double enh = 0.25 * (2. + 2. * pq.first + pq.second);
556  kap += (enh > 1.0 ? enh : 1.0);
557  nd += 1.0;
558  }
559  return kap / nd;
560
561 }
562
563 //--------------------------------------------------------------------------
564
565 // Calculate the effective string tension a fraction yfrac in on the dipole,
566 // given by indices e1 and e2.
567
568 double Ropewalk::getKappaHere(int e1, int e2, double yfrac) {
569
570  multimap< pair<int,int>, RopeDipole >::iterator
571  itr = dipoles.find( make_pair(e1,e2) );
572  if (itr == dipoles.end()) itr = dipoles.find( make_pair(e2,e1) );
573  if (itr == dipoles.end()) return 1.0;
574  RopeDipole* d = &(itr->second);
576
577  // Get quantum numbers m and n.
578  pair<int, int> overlap = d->getOverlaps(yfrac, m0, r0);
579  pair<double, double> pq;
580  // If we are always in the highest multiplet, we need not do
581  // a random walk
582  if (alwaysHighest) {
583  pq = make_pair(overlap.first + 1, overlap.second);
584  }
585  // Invoke default random walk procedure.
586  else {
587  pq = select(overlap.first + 1, overlap.second, rndmPtr);
588  }
589  // Calculate enhancement factor.
590  double enh = 0.25 * (2. + 2. * pq.first + pq.second);
591  return (enh > 1.0 ? enh : 1.0);
592
593 }
594
595 //--------------------------------------------------------------------------
596
597 // Calculate all overlaps of all dipoles and store as OverlappingRopeDipoles.
598
599 bool Ropewalk::calculateOverlaps() {
600
601  // Go through all dipoles.
602  for (multimap< pair<int,int>, RopeDipole>::iterator itr = dipoles.begin();
603  itr != dipoles.end(); ++itr ) {
604  RopeDipole* d1 = &(itr->second);
605  if (d1->dipoleMomentum().m2Calc() < pow2(m0)) continue;
606
607  // RopeDipoles rapidities in dipole rest frame.
608  RotBstMatrix dipoleRestFrame = d1->getDipoleRestFrame();
609  double yc1 = d1->d1Ptr()->rap(m0, dipoleRestFrame);
610  double ya1 = d1->d2Ptr()->rap(m0, dipoleRestFrame);
611  if (yc1 <= ya1) continue;
612
613  // Go through all possible overlapping dipoles.
614  for (multimap< pair<int,int>, RopeDipole>::iterator itr2 = dipoles.begin();
615  itr2 != dipoles.end(); ++itr2) {
616  RopeDipole* d2 = &(itr2->second);
617
618  // Skip self and overlaps with miniscule dipoles.
619  if (d1 == d2) continue;
620  if (d2->dipoleMomentum().m2Calc() < pow2(m0)) continue;
621
622  // Ignore if not overlapping in rapidity.
623  OverlappingRopeDipole od(d2, m0, dipoleRestFrame);
624  if (min(od.y1, od.y2) > yc1 || max(od.y1, od.y2) < ya1 || od.y1 == od.y2)
625  continue;
626
628
629  }
630  }
631  return true;
632
633 }
634 //--------------------------------------------------------------------------
635
636 // Invoke the random walk of colour states.
637
638 pair<int, int> Ropewalk::select(int m, int n, Rndm* rndm) {
639
640  // Initial valuesM mm and nn are step counters.
641  int p = 0, q = 0;
642  int mm = m, nn = n;
643
644  // We must take all steps before terminating.
645  while (mm + nn > 0) {
646
647  // Take randomly a step in one direction.
648  if (rndm->flat() < 0.5 && mm > 0) {
649  --mm;
650
651  // Calculate the step probabilities.
652  double p1 = multiplicity(p + 1, q);
653  double p2 = multiplicity(p, q - 1);
654  double p3 = multiplicity(p - 1, q + 1);
655
656  // Normalization.
657  double sum = p1 + p2 + p3;
658  p1 /= sum, p2 /= sum, p3 /= sum;
659
660  // Select a state.
661  double pick = rndm->flat();
662  if (pick < p1) ++p;
663  else if (pick < p1 + p2) --q;
664  else --p, ++q;
665  }
666
667  // As above, but for nn.
668  else if (nn > 0) {
669  --nn;
670  double p1 = multiplicity(p, q + 1);
671  double p2 = multiplicity(p -1, q);
672  double p3 = multiplicity(p + 1, q - 1);
673  double sum = p1 + p2 + p3;
674  p1 /= sum, p2 /= sum, p3 /= sum;
675  double pick = rndm->flat();
676  if (pick < p1) ++q;
677  else if (pick < p1 + p2) --p;
678  else --q, ++p;
679  }
680  }
681
682  // Done.
683  return make_pair( (p < 0 ? 0 : p), (q < 0 ? 0 : q) );
684
685 }
686
687 //--------------------------------------------------------------------------
688
689 // Shove all dipoles in the event.
690
691 void Ropewalk::shoveTheDipoles(Event& event) {
692
693  // Possibility of some initial propagation.
694  if ( tInit > 0.0) {
695  for (DMap::iterator dItr = dipoles.begin(); dItr != dipoles.end(); ++dItr)
696  dItr->second.propagateInit(tInit);
697  }
698
699  // The rapidity slices.
700  multimap<double, RopeDipole *> rapDipoles;
701
702  // Order the dipoles in max rapidity.
703  double ymin = 0;
704  double ymax = 0;
705  for (DMap::iterator dItr = dipoles.begin(); dItr != dipoles.end(); ++dItr) {
706  RopeDipole* dip = &(dItr->second);
707  // Order the dipoles in max rapidity.
708  rapDipoles.insert( make_pair(dip->maxRapidity(m0), dip) );
709  // Find maximal and minimal rapidity to sample.
710  if (dip->minRapidity(m0) < ymin) ymin = dip->minRapidity(m0);
711  if (dip->maxRapidity(m0) > ymax) ymax = dip->maxRapidity(m0);
712  }
713
714  // Do the sampling from flat distribution.
715  vector<double> rapidities;
716  for (double y = ymin; y < ymax; y += deltay) rapidities.push_back(y);
717
718  // For each value of ySample, we have a vector of excitation pairs.
719  map<double, vector<Exc> > exPairs;
720  for (int i = 0, N = eParticles.size(); i < N; ++i) eParticles[i].clear();
721  eParticles.clear();
722  for (int i = 0, N = rapidities.size(); i < N; ++i) {
723  // Construct an empty vector of excitation particles.
724  eParticles.push_back( vector<Particle>() );
725
726  // Find dipoles sampled in this slice, and store them temporarily.
727  double ySample = rapidities[i];
728  vector<RopeDipole*> tmp;
729  for (multimap<double, RopeDipole*>::iterator
730  rItr = rapDipoles.lower_bound(ySample); rItr != rapDipoles.end(); ++rItr) {
731  if (rItr->second->minRapidity(m0) < ySample)
732  tmp.push_back(rItr->second);
733  }
734
735  // Construct excitation particles, one for each sampled dipole in this slice.
736  vector<int> eraseDipoles;
737
738  for (int j = 0, M = tmp.size(); j < M; ++j) {
739  Vec4 ex;
740  // Test if the dipole can bear the excitation.
741  if (!tmp[j]->recoil(ex,true) ) {
742  eraseDipoles.push_back(j);
743  }
744  }
745  // Erase dipoles which could not bear an excitation.
746  for (int j = 0, M = eraseDipoles.size(); j < M; ++j) {
747  tmp.erase( tmp.begin() + (eraseDipoles[j]-j) );
748  }
749  // Add the actual excitations, but only if we can create pairs.
750  if( int(tmp.size()) > 1)
751  for (int j = 0, M = tmp.size(); j < M; ++j) {
752  Vec4 ex;
753  // We boost the excitation back from dipole rest frame.
754  tmp[j]->recoil(ex,false);
755  Particle pp = Particle(21, 22, 0, 0, 0, 0, 0, 0, ex);
756  pp.vProd( tmp[j]->bInterpolateLab(ySample,m0) );
757  eParticles[i].push_back(pp);
758  }
759  // Construct all pairs of possible excitations in this slice.
760  exPairs[ySample] = vector<Exc>();
761  for (int j = 0, M = tmp.size(); j < M; ++j)
762  for (int k = 0; k < M; ++k) {
763  // Don't allow a string to shove itself.
764  if(j != k && tmp[j]->index() != tmp[k]->index() )
765  exPairs[ySample].push_back( Exc(ySample, m0, i, j, k, tmp[j],
766  tmp[k]) );
767  }
768  }
769
770  // Give the excitations pointers to the excitation particles.
771  for (map<double, vector<Exc> >::iterator slItr = exPairs.begin();
772  slItr != exPairs.end(); ++slItr) {
773  for (int i = 0, N = slItr->second.size(); i < N; ++i) {
774  Exc& ep = slItr->second[i];
775  ep.setParticlePtrs( &eParticles[ep.i][ep.j], &eParticles[ep.i][ep.k] );
776  }
777  }
778
779  // Shoving loop.
780  for (double t = tInit; t < tShove + tInit; t += deltat) {
781  // For all slices.
782  for (map<double, vector<Exc> >::iterator slItr = exPairs.begin();
783  slItr != exPairs.end(); ++slItr)
784  // For all excitation pairs.
785  for (int i = 0, N = slItr->second.size(); i < N; ++i) {
786  Exc& ep = slItr->second[i];
787  // The direction vector is a space-time four-vector.
788  Vec4 direction = ep.direction();
789  // The string radius is time dependent,
790  // growing with the speed of light.
791  // Minimal string size is 1 / shower cut-off
792  // converted to fm.
793  double rt = max(t, 1. / showerCut / 5.068);
794  rt = min(rt, r0 * gExponent);
795  double dist = direction.pT();
796  // Calculate the push, its direction and do the shoving.
797  if (dist < rCutOff * rt) {
798  // Gain function.
799  double gain = 0.5 * deltay * deltat * gAmplitude * dist / rt / rt
800  * exp( -0.25 * dist * dist / rt / rt);
801  double dpx = dist > 0.0 ? gain * direction.px() / dist: 0.0;
802  double dpy = dist > 0.0 ? gain * direction.py() / dist: 0.0;
803  ep.shove(dpx, dpy);
804  }
805  }
806
807  // Propagate the dipoles.
808  for (DMap::iterator dItr = dipoles.begin(); dItr != dipoles.end(); ++dItr)
809  dItr->second.propagate(deltat, m0);
810  }
811
812  // Add the excitations to the dipoles.
813  for (DMap::iterator dItr = dipoles.begin(); dItr != dipoles.end(); ++dItr) {
814  RopeDipole* dip = &(dItr->second);
815  if (dip->nExcitations() > 0) dip->excitationsToString( m0, event);
816  }
817 }
818
819 //--------------------------------------------------------------------------
820
821 // Extract all dipoles from an event.
822
823 bool Ropewalk::extractDipoles(Event& event, ColConfig& colConfig) {
824
825  // Go through all strings in the event.
826  dipoles.clear();
827  for (int iSub = 0; iSub < colConfig.size(); ++iSub) {
828
829  // We can exclude loops, junctions and ministrings from the Ropewalk.
830  if (colConfig[iSub].hasJunction && !shoveJunctionStrings) continue;
831  if (colConfig[iSub].isClosed && !shoveGluonLoops) continue;
832  if (colConfig[iSub].massExcess <= mStringMin && !shoveMiniStrings)
833  continue;
834
835  colConfig.collect(iSub,event);
836  vector<int> stringPartons = colConfig[iSub].iParton;
837  vector<RopeDipole> stringDipole;
838  bool stringStart = true;
839  RopeDipoleEnd previous;
840  for (int iPar = int(stringPartons.size() - 1); iPar > -1; --iPar) {
841  if (stringPartons[iPar] > 0) {
842  // Ordinary particle.
843  RopeDipoleEnd next( &event, stringPartons[iPar]);
844  // If we are at the first parton, no dipole.
845  if ( !stringStart) {
846  // Get the parton placement in Event Record.
847  pair<int,int> dipoleER = make_pair( stringPartons[iPar + 1],
848  stringPartons[iPar] );
849  RopeDipole test(previous, next, iSub, infoPtr);
850  if ( limitMom && test.dipoleMomentum().pT() < pTcut)
851  dipoles.insert( pair< pair<int, int>, RopeDipole>(dipoleER,
852  RopeDipole( previous, next, iSub, infoPtr)) );
853  else if (!limitMom)
854  dipoles.insert( pair< pair<int, int>, RopeDipole>(dipoleER,
855  RopeDipole( previous, next, iSub, infoPtr)));
856  }
857  previous = next;
858  stringStart = false;
859  }
860  else continue;
861  }
862  // We have created all dipoles.
863  }
864  return true;
865
866 }
867
868 //==========================================================================
869
870 // RopeFragPars recalculates fragmentation parameters according to a
871 // changed string tension. Helper class to FlavourRope.
872
873 //--------------------------------------------------------------------------
874
875 // Constants: could be changed here if desired, but normally should not.
876
877 // Initial step size for a calculation.
878 const double RopeFragPars::DELTAA = 0.1;
879
880 // Convergence criterion for a calculation.
881 const double RopeFragPars::ACONV = 0.001;
882
883 // Low z cut-off in fragmentation function.
884 const double RopeFragPars::ZCUT = 1.0e-4;
885
886 //--------------------------------------------------------------------------
887
888 // The init function sets up initial parameters from settings.
889
890 void RopeFragPars::init(Info* infoPtrIn, Settings& settings) {
891
892  // Info pointer.
893  infoPtr = infoPtrIn;
894
895  // The junction parameter.
896  beta = settings.parm("Ropewalk:beta");
897
898  // Initialize default values from input settings.
899  const int len = 9;
900  string params [len] = {"StringPT:sigma", "StringZ:aLund",
902  "StringFlav:probSQtoQQ", "StringFlav:probQQ1toQQ0", "StringFlav:probQQtoQ",
903  "StringFlav:kappa"};
904  double* variables[len] = {&sigmaIn, &aIn, &adiqIn, &bIn, &rhoIn, &xIn,
905  &yIn, &xiIn, &kappaIn};
906  for (int i = 0; i < len; ++i) *variables[i] = settings.parm(params[i]);
907
908  // Insert the h = 1 case immediately.
909  sigmaEff = sigmaIn, aEff = aIn, adiqEff = adiqIn, bEff = bIn,
910  rhoEff = rhoIn, xEff = xIn, yEff = yIn, xiEff = xiIn, kappaEff = kappaIn;
911  if (!insertEffectiveParameters(1.0)) infoPtr->errorMsg(
912  "Error in RopeFragPars::init: failed to insert defaults.");
913
914 }
915
916 //--------------------------------------------------------------------------
917
918 // Return parameters at given string tension, ordered by their
919 // name for easy insertion in settings.
920
921 map<string,double> RopeFragPars::getEffectiveParameters(double h) {
922
923  map<double, map<string, double> >::iterator parItr = parameters.find(h);
924
925  // If the parameters are already calculated, return them.
926  if ( parItr != parameters.end()) return parItr->second;
927
928  // Otherwise calculate them.
929  if (!calculateEffectiveParameters(h))
930  infoPtr->errorMsg("Error in RopeFragPars::getEffectiveParameters:"
931  " calculating effective parameters.");
932
933  // And insert them.
934  if (!insertEffectiveParameters(h))
935  infoPtr->errorMsg("Error in RopeFragPars::getEffectiveParameters:"
936  " inserting effective parameters.");
937
938  // And recurse.
939  return getEffectiveParameters(h);
940
941 }
942
943 //--------------------------------------------------------------------------
944
945 // Get the Fragmentation function a parameter from cache or calculate it.
946
947 double RopeFragPars::getEffectiveA(double thisb, double mT2, bool isDiquark) {
948
949  // Check for the trivial case.
950  if (thisb == bIn) return (isDiquark ? aIn + adiqIn : aIn);
951
952  // We order by b*mT2
953  map<double, double>* aMPtr = (isDiquark ? &aDiqMap : &aMap);
954  double bmT2 = mT2 * thisb;
955
956  // Check if we have already calculated this a value before.
957  map<double,double>::iterator aItr = aMPtr->find(bmT2);
958  if (aItr != aMPtr->end()) return aItr->second;
959
960  // Otherwise calculate it.
961  double ae = ( isDiquark ? aEffective(aIn + adiqIn, thisb, mT2)
962  : aEffective(aIn, thisb, mT2) );
963  if (isDiquark) {
964  double suba = getEffectiveA(thisb, mT2, false);
965  aMPtr->insert( make_pair(bmT2, ae - suba) );
966  }
967  else aMPtr->insert(make_pair(bmT2, ae));
968  return ae;
969
970 }
971
972 //--------------------------------------------------------------------------
973
974 // Calculate the effective parameters.
975
976 bool RopeFragPars::calculateEffectiveParameters(double h) {
977
978  if (h <= 0) return false;
979  double hinv = 1.0 / h;
980
982  // The string tension kappa.
983  kappaEff = kappaIn * h;
984  // Strangeness.
985  rhoEff = pow(rhoIn, hinv);
986  // Strange diquarks.
987  xEff = pow(xIn, hinv);
988  // Spin.
989  yEff = pow(yIn, hinv);
990  // pT distribution width.
991  sigmaEff = sigmaIn * sqrt(h);
992
993  // Derived quantity alpha.
994  double alpha = (1 + 2. * xIn * rhoIn + 9. * yIn + 6. * xIn * rhoIn * yIn
995  + 3. * yIn * xIn * xIn * rhoIn * rhoIn) / (2. + rhoIn);
996  double alphaEff = (1. + 2. * xEff * rhoEff + 9. * yEff
997  + 6. * xEff * rhoEff * yEff + 3. * yEff * xEff * xEff * rhoEff * rhoEff)
998  / (2. + rhoEff);
999
1000  // Baryons.
1001  xiEff = alphaEff * beta * pow( xiIn / alpha / beta, hinv);
1002  if (xiEff > 1.0) xiEff = 1.0;
1003  if (xiEff < xiIn) xiEff = xiIn;
1004
1005  // Fragmentation function b.
1006  bEff = (2. + rhoEff) / (2. + rhoIn) * bIn;
1007  if (bEff < bIn) bEff = bIn;
1008  if (bEff > 2.0) bEff = 2.0;
1009
1010  // We calculate a for a typical particle with mT2 = 1 GeV^2.
1011  aEff = getEffectiveA( bEff, 1.0, false);
1012  adiqEff = getEffectiveA( bEff, 1.0, true) - aEff;
1013
1014  return true;
1015
1016 }
1017
1018 //--------------------------------------------------------------------------
1019
1020 // Insert calculated parameters in cache for later (re-)use.
1021
1022 bool RopeFragPars::insertEffectiveParameters(double h) {
1023
1024  map<string,double> p;
1025  p["StringPT:sigma"] = sigmaEff;
1026  p["StringZ:bLund"] = bEff;
1027  p["StringFlav:probStoUD"] = rhoEff;
1028  p["StringFlav:probSQtoQQ"] = xEff;
1029  p["StringFlav:probQQ1toQQ0"] = yEff;
1030  p["StringFlav:probQQtoQ"] = xiEff;
1031  p["StringZ:aLund"] = aEff;
1033  p["StringFlav:kappa"] = kappaEff;
1034
1035  return (parameters.insert( make_pair(h,p)).second );
1036
1037 }
1038
1039 //--------------------------------------------------------------------------
1040
1041 // Calculate the a parameter.
1042
1043 double RopeFragPars::aEffective(double aOrig, double thisb, double mT2) {
1044
1045  // Calculate initial normalization constants.
1046  double N = integrateFragFun(aOrig, bIn, mT2);
1047  double NEff = integrateFragFun(aOrig, thisb, mT2);
1048  int s = (N < NEff) ? -1 : 1;
1049  double da = DELTAA;
1050  double aNew = aOrig - s * da;
1051
1052  // Iterate until we meet preset convergence criterion.
1053  do {
1054  // Calculate normalization with current a.
1055  NEff = integrateFragFun(aNew, thisb, mT2);
1056  if ( ((N < NEff) ? -1 : 1) != s ) {
1057  s = (N < NEff) ? -1 : 1;
1058  // If we have crossed over the solution, decrease
1059  // the step size and turn around.
1060  da /= 10.0;
1061  }
1062  aNew -= s * da;
1063  if (aNew < 0.0) {aNew = 0.1; break;}
1064  if (aNew > 2.0) {aNew = 2.0; break;}
1065  } while (da > ACONV);
1066  return aNew;
1067
1068 }
1069
1070 //--------------------------------------------------------------------------
1071
1072 // The Lund fragmentation function.
1073
1074 double RopeFragPars::fragf(double z, double a, double b, double mT2) {
1075
1076  if (z < ZCUT) return 0.0;
1077  return pow(1 - z, a) * exp(-b * mT2 / z) / z;
1078
1079 }
1080
1081 //--------------------------------------------------------------------------
1082
1083 // Integral of the Lund fragmentation function with given parameter values.
1084
1085 double RopeFragPars::integrateFragFun(double a, double b, double mT2) {
1086
1087  // Using Simpson's rule to integrate the Lund fragmentation function.
1088  double nextIter, nextComb;
1089  double thisComb = 0.0, thisIter = 0.0;
1090  // The target error on the integral should never be changed.
1091  double error = 1.0e-2;
1092
1093  // 20 is the max number of iterations, 3 is min. Should not be changed.
1094  for (int i = 1; i < 20; ++i) {
1095  nextIter = trapIntegrate( a, b, mT2, thisIter, i);
1096  nextComb = (4.0 * nextIter - thisIter) / 3.0;
1097  if (i > 3 && abs(nextComb - thisComb) < error * abs(nextComb))
1098  return nextComb;
1099  thisIter = nextIter;
1100  thisComb = nextComb;
1101  }
1102  infoPtr->errorMsg("RopeFragPars::integrateFragFun:"
1103  "No convergence of frag fun integral.");
1104  return 0.0;
1105
1106 }
1107
1108 //--------------------------------------------------------------------------
1109
1110 // Helper function for integration.
1111
1112 double RopeFragPars::trapIntegrate( double a, double b, double mT2,
1113  double sOld, int n) {
1114
1115  // Compute the nth correction to the integral of fragfunc between 0 and 1
1116  // using extended trapezoidal rule.
1117  if (n == 1) return 0.5 * (fragf(0.0, a, b, mT2) + fragf(1.0, a, b, mT2));
1118  // We want 2^(n-2) interior points (intp). Use bitwise shift to speed up.
1119  int intp = 1;
1120  intp <<= n - 2;
1121  double deltaz = 1.0 / double(intp);
1122  double z = 0.5 * deltaz;
1123  double sum = 0.0;
1124  // Do the integral.
1125  for (int i = 0; i < intp; ++i, z += deltaz) sum += fragf( z, a, b, mT2);
1126  return 0.5 * (sOld + sum / double(intp));
1127
1128 }
1129
1130 //==========================================================================
1131
1132 // The FlavourRope class takes care of placing a string breakup in
1133 // the event, and assigning the string breakup effective parameters.
1134
1135 //--------------------------------------------------------------------------
1136
1137 // Change the fragmentation parameters.
1138
1139 bool FlavourRope::doChangeFragPar(StringFlav* flavPtr, StringZ* zPtr,
1140  StringPT * pTPtr, double m2Had, vector<int> iParton, int endId) {
1141
1142  // The new parameters.
1143  map<string, double> newPar;
1144  if (doBuffon)
1145  newPar = fetchParametersBuffon(m2Had, iParton, endId);
1146  else
1147  newPar = fetchParameters(m2Had, iParton, endId);
1148  // Change settings to new settings.
1149  for (map<string, double>::iterator itr = newPar.begin(); itr!=newPar.end();
1150  ++itr) settingsPtr->parm( itr->first, itr->second);
1151  // Re-initialize flavour, z, and pT selection with new settings.
1152  flavPtr->init( *settingsPtr, particleDataPtr, rndmPtr, infoPtr);
1153  zPtr->init( *settingsPtr, *particleDataPtr, rndmPtr, infoPtr);
1154  pTPtr->init( *settingsPtr, particleDataPtr, rndmPtr, infoPtr);
1155  return true;
1156
1157 }
1158
1159 //--------------------------------------------------------------------------
1160
1161 // Find breakup placement and fetch effective parameters using Buffon.
1162
1164  vector<int> iParton, int endId) {
1165  // If effective string tension is set manually, use that.
1166  if (fixedKappa) return fp.getEffectiveParameters(h);
1167  if (!ePtr) {
1168  infoPtr->errorMsg("Error in FlavourRope::fetchParametersBuffon:"
1169  " Event pointer not set in FlavourRope");
1170  return fp.getEffectiveParameters(1.0);
1171  }
1176  }
1177  // Quark string ends, default mode
1178  if (endId != 21){
1179  // Test consistency
1180  if(ePtr->at(*(iParton.begin())).id() != endId &&
1181  ePtr->at(*(iParton.end() - 1)).id() != endId) {
1182  infoPtr->errorMsg("Error in FlavourRope::fetchParametersBuffon:"
1183  " Quark end inconsistency.");
1184  return fp.getEffectiveParameters(1.0);
1185  }
1186
1187  // First we must let the string vector point in the right direction
1188  if(ePtr->at(*(iParton.begin())).id() != endId)
1189  reverse(iParton.begin(),iParton.end());
1190
1191  // Initialize a bit
1193  double enh = 1.0;
1194  double dipFrac;
1195  vector<int>::iterator dipItr;
1196  // Find out when invariant mass exceeds m2Had
1197  for(dipItr = iParton.begin(); dipItr != iParton.end(); ++dipItr){
1200  // Approximate the fraction we are in on the dipole, this goes
1201  // in three cases.
1202  // We are at the beginning.
1204  dipFrac = 0;
1205  }
1206  // We are somewhere in the first dipole
1207  else if(dipItr - 1 == iParton.begin()){
1209  }
1210  else{
1211  if(ePtr->at(*(dipItr - 1)).id() != 21) {
1212  infoPtr->errorMsg("Error in FlavourRope::fetchParametersBuffon:"
1213  " Connecting partons should always be gluons.");
1214  return fp.getEffectiveParameters(1.0);
1215  }
1216
1219
1220  dipFrac = (sqrt(m2Had) - sqrt(m2Small)) /
1221  (sqrt(m2Big) - sqrt(m2Small));
1222  }
1223  break;
1224  }
1225  hadronic4Momentum += ePtr->at(*dipItr).id() == 21 ?
1226  0.5*ePtr->at(*dipItr).p() : ePtr->at(*dipItr).p();
1227  }
1228  // If we reached the end
1229  // we are in a small string that should just be collapsed
1230  if(dipItr == iParton.end())
1231  return fp.getEffectiveParameters(1.0);
1232  // Sanity check
1233  if(dipFrac < 0 || dipFrac > 1) {
1234  infoPtr->errorMsg("Error in FlavourRope::fetchParametersBuffon:"
1235  " Dipole exceed with fraction less than 0 or greater than 1.");
1236  return fp.getEffectiveParameters(1.0);
1237  }
1238  // We now figure out at what rapidity value,
1239  // in the lab system, the string is breaking
1240  double yBreak;
1241  // Trivial case, just inherit
1242  if(dipFrac == 0)
1243  yBreak = ePtr->at(*dipItr).y();
1244  else{
1245  // Sanity check
1246  if(dipItr == iParton.begin()) {
1247  infoPtr->errorMsg("Error in FlavourRope::fetchParametersBuffon:"
1248  " We are somehow before the first dipole on a string.");
1249  return fp.getEffectiveParameters(1.0);
1250  }
1251  double dy = ePtr->at(*dipItr).y() - ePtr->at(*(dipItr - 1)).y();
1252  yBreak = ePtr->at(*(dipItr - 1)).y() + dipFrac*dy;
1253  }
1254  // Count the number of partons in the whole
1255  // event within deltay of breaking point
1256  double p = 1;
1257  double q = 0;
1258
1259  for(int i = 0; i < ePtr->size(); ++i){
1260  // Don't double count partons from this
1261  // string (no self-overlap)
1262  if(find(iParton.begin(),iParton.end(), i) != iParton.end())
1263  continue;
1266  continue;
1267  double pRap = ePtr->at(i).y();
1268  if(pRap > yBreak - rapiditySpan && pRap < yBreak + rapiditySpan ){
1269  // Do a "Buffon" selection to decide whether
1270  // two strings overlap in impact parameter space
1271  // given ratio of string diameter to collision diameter
1272  double r1 = rndmPtr->flat();
1273  double r2 = rndmPtr->flat();
1274  double theta1 = 2*M_PI*rndmPtr->flat();
1275  double theta2 = 2*M_PI*rndmPtr->flat();
1276  // Overlap?
1277  if(4*pow2(stringProtonRatio) > pow2(sqrt(r1)*cos(theta1) -
1278  sqrt(r2)*cos(theta2)) + pow2(sqrt(r1)*sin(theta1) -
1279  sqrt(r2)*sin(theta2))) {
1280  if(rndmPtr->flat() < 0.5) p += 0.5;
1281  else q += 0.5;
1282  }
1283  }
1284  }
1285  enh = 0.25*(2.0*p+q+2.0);
1286
1287  return fp.getEffectiveParameters(enh);
1288  }
1289  // For closed gluon loops we cannot distinguish the ends.
1291  else{
1292  return fp.getEffectiveParameters(1.0);
1293  }
1294  return fp.getEffectiveParameters(1.0);
1295 }
1296
1297 //--------------------------------------------------------------------------
1298 // Find breakup placement and fetch effective parameters using Ropewalk.
1299
1301  vector<int> iParton, int endId) {
1302  // If effective string tension is set manually, use that.
1303  if (fixedKappa) return fp.getEffectiveParameters(h);
1304  if (!ePtr) {
1305  infoPtr->errorMsg("Error in FlavourRope::fetchParameters:"
1306  " Event pointer not set in FlavourRope");
1307  return fp.getEffectiveParameters(1.0);
1308  }
1309  Vec4 mom;
1310  int eventIndex = -1;
1311  // Set direction
1312  bool dirPos;
1313  if( ePtr->at(iParton[0]).id() == endId) dirPos = true;
1314  else if( ePtr->at(iParton[iParton.size() - 1]).id() == endId) dirPos = false;
1315  else {
1316  infoPtr->errorMsg("Error in FlavourRope::fetchParameters:"
1317  " Could not get string direction");
1318  return fp.getEffectiveParameters(1.0);
1319  }
1320
1321  for (int i = 0, N = iParton.size(); i < N; ++i) {
1322  // Change to right direction
1323  int j = (dirPos ? i : N - 1 - i);
1324  // Skip the junction entry.
1325  if ( iParton[j] < 0) continue;
1326  mom += ePtr->at(iParton[j]).p();
1327  if ( mom.m2Calc() > m2Had) {
1328  eventIndex = j;
1329  break;
1330  }
1331  }
1332
1333  // We figure out where we are on the dipole.
1334  // Set some values.
1335  double m2Here = mom.m2Calc();
1336  // The dipFrac signifies a fraction.
1337  double dipFrac = 0;
1338  // We are in the first dipole.
1339  if (eventIndex == -1 || eventIndex == 0) {
1340  eventIndex = 0;
1341  dipFrac = sqrt(m2Had / m2Here);
1342  }
1343  else {
1344  mom -= ePtr->at(iParton[eventIndex]).p();
1345  double m2Small = mom.m2Calc();
1346  dipFrac = (sqrt(m2Had) - sqrt(m2Small)) / (sqrt(m2Here) - sqrt(m2Small));
1347  }
1348  double enh = rwPtr->getKappaHere( iParton[eventIndex],
1349  iParton[eventIndex + 1], dipFrac);
1350  return fp.getEffectiveParameters(enh);
1351
1352 }
1353
1354
1355 //==========================================================================
1356
1357 } // End namespace Pythia8
void ae(int tracks=-1, int hits=-1)
This function is to search for the next non-empty event and draw it by looping over StBFChain (readin...
Definition: Ed.C:102
Definition: rb.hh:21
Definition: AgUStep.h:26