Line data Source code
1 : // StringFragmentation.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 :
6 : // Function definitions (not found in the header) for the StringEnd and
7 : // StringFragmentation classes.
8 :
9 : #include "Pythia8/StringFragmentation.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // The StringEnd class.
16 :
17 : //--------------------------------------------------------------------------
18 :
19 : // Constants: could be changed here if desired, but normally should not.
20 :
21 : // Avoid unphysical solutions to equation system.
22 : const double StringEnd::TINY = 1e-6;
23 :
24 : // Assume two (eX, eY) regions are related if pT2 differs by less.
25 : const double StringEnd::PT2SAME = 0.01;
26 :
27 : //--------------------------------------------------------------------------
28 :
29 : // Set up initial endpoint values from input.
30 :
31 : void StringEnd::setUp(bool fromPosIn, int iEndIn, int idOldIn, int iMaxIn,
32 : double pxIn, double pyIn, double GammaIn, double xPosIn, double xNegIn) {
33 :
34 : // Simple transcription from input.
35 0 : fromPos = fromPosIn;
36 0 : iEnd = iEndIn;
37 0 : iMax = iMaxIn;
38 0 : flavOld = FlavContainer(idOldIn);
39 0 : pxOld = pxIn;
40 0 : pyOld = pyIn;
41 0 : GammaOld = GammaIn;
42 0 : iPosOld = (fromPos) ? 0 : iMax;
43 0 : iNegOld = (fromPos) ? iMax : 0;
44 0 : xPosOld = xPosIn;
45 0 : xNegOld = xNegIn;
46 :
47 0 : }
48 :
49 : //--------------------------------------------------------------------------
50 :
51 : // Fragment off one hadron from the string system, in flavour and pT.
52 :
53 : void StringEnd::newHadron() {
54 :
55 : // Pick new flavour and form a new hadron.
56 0 : do {
57 0 : flavNew = flavSelPtr->pick( flavOld);
58 0 : idHad = flavSelPtr->combine( flavOld, flavNew);
59 0 : } while (idHad == 0);
60 :
61 : // Pick its transverse momentum.
62 0 : pair<double, double> pxy = pTSelPtr->pxy();
63 0 : pxNew = pxy.first;
64 0 : pyNew = pxy.second;
65 0 : pxHad = pxOld + pxNew;
66 0 : pyHad = pyOld + pyNew;
67 :
68 : // Pick its mass and thereby define its transverse mass.
69 0 : mHad = particleDataPtr->mSel(idHad);
70 0 : mT2Had = pow2(mHad) + pow2(pxHad) + pow2(pyHad);
71 :
72 0 : }
73 :
74 : //--------------------------------------------------------------------------
75 :
76 : // Fragment off one hadron from the string system, in momentum space,
77 : // by taking steps from positive end.
78 :
79 : Vec4 StringEnd::kinematicsHadron( StringSystem& system) {
80 :
81 : // Pick fragmentation step z and calculate new Gamma.
82 0 : zHad = zSelPtr->zFrag( flavOld.id, flavNew.id, mT2Had);
83 0 : GammaNew = (1. - zHad) * (GammaOld + mT2Had / zHad);
84 :
85 : // Set up references that are direction-neutral;
86 : // ...Dir for direction of iteration and ...Inv for its inverse.
87 0 : int& iDirOld = (fromPos) ? iPosOld : iNegOld;
88 0 : int& iInvOld = (fromPos) ? iNegOld : iPosOld;
89 0 : int& iDirNew = (fromPos) ? iPosNew : iNegNew;
90 0 : int& iInvNew = (fromPos) ? iNegNew : iPosNew;
91 0 : double& xDirOld = (fromPos) ? xPosOld : xNegOld;
92 0 : double& xInvOld = (fromPos) ? xNegOld : xPosOld;
93 0 : double& xDirNew = (fromPos) ? xPosNew : xNegNew;
94 0 : double& xInvNew = (fromPos) ? xNegNew : xPosNew;
95 0 : double& xDirHad = (fromPos) ? xPosHad : xNegHad;
96 0 : double& xInvHad = (fromPos) ? xNegHad : xPosHad;
97 :
98 : // Start search for new breakup in the old region.
99 0 : iDirNew = iDirOld;
100 0 : iInvNew = iInvOld;
101 0 : Vec4 pTNew;
102 :
103 : // Each step corresponds to trying a new string region.
104 0 : for (int iStep = 0; ; ++iStep) {
105 :
106 : // Referance to current string region.
107 0 : StringRegion& region = system.region( iPosNew, iNegNew);
108 :
109 : // Now begin special section for rapid processing of low region.
110 0 : if (iStep == 0 && iPosOld + iNegOld == iMax) {
111 :
112 : // A first step within a low region is easy.
113 0 : if (mT2Had < zHad * xDirOld * (1. - xInvOld) * region.w2) {
114 :
115 : // Translate into x coordinates.
116 0 : xDirHad = zHad * xDirOld;
117 0 : xInvHad = mT2Had / (xDirHad * region.w2);
118 0 : xDirNew = xDirOld - xDirHad;
119 0 : xInvNew = xInvOld + xInvHad;
120 :
121 : // Find and return four-momentum of the produced particle.
122 0 : return region.pHad( xPosHad, xNegHad, pxHad, pyHad);
123 :
124 : // A first step out of a low region also OK, if there are more regions.
125 : // Negative energy signals failure, i.e. in last region.
126 : } else {
127 0 : --iInvNew;
128 0 : if (iInvNew < 0) return Vec4(0., 0., 0., -1.);
129 :
130 : // Momentum taken by stepping out of region. Continue to next region.
131 0 : xInvHad = 1. - xInvOld;
132 0 : xDirHad = 0.;
133 0 : pSoFar = region.pHad( xPosHad, xNegHad, pxOld, pyOld);
134 0 : continue;
135 : }
136 :
137 : // Else, for first step, take into account starting pT.
138 0 : } else if (iStep == 0) {
139 0 : pSoFar = region.pHad( 0., 0., pxOld, pyOld);
140 0 : pTNew = region.pHad( 0., 0., pxNew, pyNew);
141 0 : }
142 :
143 : // Now begin normal treatment of nontrivial regions.
144 : // Set up four-vectors in a region not visited before.
145 0 : if (!region.isSetUp) region.setUp(
146 0 : system.regionLowPos(iPosNew).pPos,
147 0 : system.regionLowNeg(iNegNew).pNeg, true);
148 :
149 : // If new region is vanishingly small, continue immediately to next.
150 : // Negative energy signals failure to do this, i.e. moved too low.
151 0 : if (region.isEmpty) {
152 0 : xDirHad = (iDirNew == iDirOld) ? xDirOld : 1.;
153 0 : xInvHad = 0.;
154 0 : pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.);
155 0 : ++iDirNew;
156 0 : if (iDirNew + iInvNew > iMax) return Vec4(0., 0., 0., -1.);
157 0 : continue;
158 : }
159 :
160 : // Reexpress pTNew w.r.t. base vectors in new region, if possible.
161 : // Recall minus sign from normalization e_x^2 = e_y^2 = -1.
162 0 : double pxNewTemp = -pTNew * region.eX;
163 0 : double pyNewTemp = -pTNew * region.eY;
164 0 : if (abs( pxNewTemp * pxNewTemp + pyNewTemp * pyNewTemp
165 0 : - pxNew * pxNew - pyNew * pyNew) < PT2SAME) {
166 0 : pxNew = pxNewTemp;
167 0 : pyNew = pyNewTemp;
168 0 : }
169 :
170 : // Four-momentum taken so far, including new pT.
171 0 : Vec4 pTemp = pSoFar + region.pHad( 0., 0., pxNew, pyNew);
172 :
173 : // Derive coefficients for m2 expression.
174 : // cM2 * x+ + cM3 * x- + cM4 * x+ * x- = m^2 - cM1;
175 0 : double cM1 = pTemp.m2Calc();
176 0 : double cM2 = 2. * (pTemp * region.pPos);
177 0 : double cM3 = 2. * (pTemp * region.pNeg);
178 0 : double cM4 = region.w2;
179 0 : if (!fromPos) swap( cM2, cM3);
180 :
181 : // Derive coefficients for Gamma expression.
182 : // cGam2 * x+ + cGam3 * x- + cGam4 * x+ * x- = Gamma_new - cGam1;
183 : double cGam1 = 0.;
184 : double cGam2 = 0.;
185 : double cGam3 = 0.;
186 : double cGam4 = 0.;
187 0 : for (int iInv = iInvNew; iInv <= iMax - iDirNew; ++iInv) {
188 : double xInv = 1.;
189 0 : if (iInv == iInvNew) xInv = (iInvNew == iInvOld) ? xInvOld : 0.;
190 0 : for (int iDir = iDirNew; iDir <= iMax - iInv; ++iDir) {
191 0 : double xDir = (iDir == iDirOld) ? xDirOld : 1.;
192 0 : int iPos = (fromPos) ? iDir : iInv;
193 0 : int iNeg = (fromPos) ? iInv : iDir;
194 0 : StringRegion& regionGam = system.region( iPos, iNeg);
195 0 : if (!regionGam.isSetUp) regionGam.setUp(
196 0 : system.regionLowPos(iPos).pPos,
197 0 : system.regionLowNeg(iNeg).pNeg, true);
198 0 : double w2 = regionGam.w2;
199 0 : cGam1 += xDir * xInv * w2;
200 0 : if (iDir == iDirNew) cGam2 -= xInv * w2;
201 0 : if (iInv == iInvNew) cGam3 += xDir * w2;
202 0 : if (iDir == iDirNew && iInv == iInvNew) cGam4 -= w2;
203 : }
204 : }
205 :
206 : // Solve (m2, Gamma) equation system => r2 * x-^2 + r1 * x- + r0 = 0.
207 0 : double cM0 = pow2(mHad) - cM1;
208 0 : double cGam0 = GammaNew - cGam1;
209 0 : double r2 = cM3 * cGam4 - cM4 * cGam3;
210 0 : double r1 = cM4 * cGam0 - cM0 * cGam4 + cM3 * cGam2 - cM2 * cGam3;
211 0 : double r0 = cM2 * cGam0 - cM0 * cGam2;
212 0 : double root = sqrtpos( r1*r1 - 4. * r2 * r0 );
213 0 : if (abs(r2) < TINY || root < TINY) return Vec4(0., 0., 0., -1.);
214 0 : xInvHad = 0.5 * (root / abs(r2) - r1 / r2);
215 0 : xDirHad = (cM0 - cM3 * xInvHad) / (cM2 + cM4 * xInvHad);
216 :
217 : // Define position of new trial vertex.
218 0 : xDirNew = (iDirNew == iDirOld) ? xDirOld - xDirHad : 1. - xDirHad;
219 0 : xInvNew = (iInvNew == iInvOld) ? xInvOld + xInvHad : xInvHad;
220 :
221 : // Step up to new region if new x- > 1.
222 0 : if (xInvNew > 1.) {
223 0 : xInvHad = (iInvNew == iInvOld) ? 1. - xInvOld : 1.;
224 0 : xDirHad = 0.;
225 0 : pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.);
226 0 : --iInvNew;
227 0 : if (iInvNew < 0) return Vec4(0., 0., 0., -1.);
228 0 : continue;
229 :
230 : // Step down to new region if new x+ < 0.
231 0 : } else if (xDirNew < 0.) {
232 0 : xDirHad = (iDirNew == iDirOld) ? xDirOld : 1.;
233 0 : xInvHad = 0.;
234 0 : pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.);
235 0 : ++iDirNew;
236 0 : if (iDirNew + iInvNew > iMax) return Vec4(0., 0., 0., -1.);
237 0 : continue;
238 : }
239 :
240 : // Else we have found the correct region, and can return four-momentum.
241 0 : return pSoFar + region.pHad( xPosHad, xNegHad, pxNew, pyNew);
242 :
243 : // End of "infinite" loop of stepping to new region.
244 0 : }
245 :
246 0 : }
247 :
248 : //--------------------------------------------------------------------------
249 :
250 : // Update string end information after a hadron has been removed.
251 :
252 : void StringEnd::update() {
253 :
254 0 : flavOld.anti(flavNew);
255 0 : iPosOld = iPosNew;
256 0 : iNegOld = iNegNew;
257 0 : pxOld = -pxNew;
258 0 : pyOld = -pyNew;
259 0 : GammaOld = GammaNew;
260 0 : xPosOld = xPosNew;
261 0 : xNegOld = xNegNew;
262 :
263 0 : }
264 :
265 : //==========================================================================
266 :
267 : // The StringFragmentation class.
268 :
269 : //--------------------------------------------------------------------------
270 :
271 : // Constants: could be changed here if desired, but normally should not.
272 : // These are of technical nature, as described for each.
273 :
274 : // Maximum number of tries to (flavour-, energy) join the two string ends.
275 : const int StringFragmentation::NTRYFLAV = 10;
276 : const int StringFragmentation::NTRYJOIN = 30;
277 :
278 : // The last few times gradually increase the stop mass to make it easier.
279 : const int StringFragmentation::NSTOPMASS = 15;
280 : const double StringFragmentation::FACSTOPMASS = 1.05;
281 :
282 : // For closed string, pick a Gamma by taking a step with fictitious mass.
283 : const double StringFragmentation::CLOSEDM2MAX = 25.;
284 : const double StringFragmentation::CLOSEDM2FRAC = 0.1;
285 :
286 : // Do not allow too large argument to exp function.
287 : const double StringFragmentation::EXPMAX = 50.;
288 :
289 : // Matching criterion that p+ and p- not the same (can happen in gg loop).
290 : const double StringFragmentation::MATCHPOSNEG = 1e-4;
291 :
292 : // For pull on junction, do not trace too far down each leg.
293 : const double StringFragmentation::EJNWEIGHTMAX = 10.;
294 :
295 : // Iterate junction rest frame boost until convergence or too many tries.
296 : const double StringFragmentation::CONVJNREST = 1e-5;
297 : const int StringFragmentation::NTRYJNREST = 20;
298 :
299 : // Fail and try again when two legs combined to diquark (3 loops).
300 : const int StringFragmentation::NTRYJNMATCH = 20;
301 : const double StringFragmentation::EEXTRAJNMATCH = 0.5;
302 : const double StringFragmentation::MDIQUARKMIN = -2.0;
303 :
304 : // Consider junction-leg parton as massless if m2 tiny.
305 : const double StringFragmentation::M2MAXJRF = 1e-4;
306 :
307 : // Protect against numerical precision giving zero or negative m2.
308 : const double StringFragmentation::M2MINJRF = 1e-4;
309 :
310 : // Iterate junction rest frame equation until convergence or too many tries.
311 : const double StringFragmentation::CONVJRFEQ = 1e-12;
312 : const int StringFragmentation::NTRYJRFEQ = 40;
313 :
314 : //--------------------------------------------------------------------------
315 :
316 : // Initialize and save pointers.
317 :
318 : void StringFragmentation::init(Info* infoPtrIn, Settings& settings,
319 : ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, StringFlav* flavSelPtrIn,
320 : StringPT* pTSelPtrIn, StringZ* zSelPtrIn) {
321 :
322 : // Save pointers.
323 0 : infoPtr = infoPtrIn;
324 0 : particleDataPtr = particleDataPtrIn;
325 0 : rndmPtr = rndmPtrIn;
326 0 : flavSelPtr = flavSelPtrIn;
327 0 : pTSelPtr = pTSelPtrIn;
328 0 : zSelPtr = zSelPtrIn;
329 :
330 : // Initialize the StringFragmentation class.
331 0 : stopMass = zSelPtr->stopMass();
332 0 : stopNewFlav = zSelPtr->stopNewFlav();
333 0 : stopSmear = zSelPtr->stopSmear();
334 0 : eNormJunction = settings.parm("StringFragmentation:eNormJunction");
335 0 : eBothLeftJunction
336 0 : = settings.parm("StringFragmentation:eBothLeftJunction");
337 0 : eMaxLeftJunction
338 0 : = settings.parm("StringFragmentation:eMaxLeftJunction");
339 0 : eMinLeftJunction
340 0 : = settings.parm("StringFragmentation:eMinLeftJunction");
341 :
342 : // Joining of nearby partons along the string.
343 0 : mJoin = settings.parm("FragmentationSystems:mJoin");
344 :
345 : // Initialize the b parameter of the z spectrum, used when joining jets.
346 0 : bLund = zSelPtr->bAreaLund();
347 :
348 : // Initialize the hadrons instance of an event record.
349 0 : hadrons.init( "(string fragmentation)", particleDataPtr);
350 :
351 : // Send on pointers to the two StringEnd instances.
352 0 : posEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr);
353 0 : negEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr);
354 :
355 0 : }
356 :
357 : //--------------------------------------------------------------------------
358 :
359 : // Perform the fragmentation.
360 :
361 : bool StringFragmentation::fragment( int iSub, ColConfig& colConfig,
362 : Event& event) {
363 :
364 : // Find partons and their total four-momentum.
365 0 : iParton = colConfig[iSub].iParton;
366 0 : iPos = iParton[0];
367 0 : if (iPos < 0) iPos = iParton[1];
368 0 : int idPos = event[iPos].id();
369 0 : iNeg = iParton.back();
370 0 : int idNeg = event[iNeg].id();
371 0 : pSum = colConfig[iSub].pSum;
372 :
373 : // Reset the local event record.
374 0 : hadrons.clear();
375 :
376 : // For closed gluon string: pick first breakup region.
377 0 : isClosed = colConfig[iSub].isClosed;
378 0 : if (isClosed) iParton = findFirstRegion(iParton, event);
379 :
380 : // For junction topology: fragment off two of the string legs.
381 : // Then iParton overwritten to remaining leg + leftover diquark.
382 0 : pJunctionHadrons = 0.;
383 0 : hasJunction = colConfig[iSub].hasJunction;
384 0 : if (hasJunction && !fragmentToJunction(event)) return false;
385 0 : int junctionHadrons = hadrons.size();
386 0 : if (hasJunction) {
387 0 : idPos = event[ iParton[0] ].id();
388 0 : idNeg = event.back().id();
389 0 : pSum -= pJunctionHadrons;
390 0 : }
391 :
392 : // Set up kinematics of string evolution ( = motion).
393 0 : system.setUp(iParton, event);
394 0 : stopMassNow = stopMass;
395 : int nExtraJoin = 0;
396 :
397 : // Fallback loop, when joining in the middle fails. Bailout if stuck.
398 0 : for ( int iTry = 0; ; ++iTry) {
399 0 : if (iTry > NTRYJOIN) {
400 0 : infoPtr->errorMsg("Error in StringFragmentation::fragment: "
401 : "stuck in joining");
402 0 : if (hasJunction) ++nExtraJoin;
403 0 : if (nExtraJoin > 0) event.popBack(nExtraJoin);
404 0 : return false;
405 : }
406 :
407 : // After several failed tries join some (extra) nearby partons.
408 0 : if (iTry == NTRYJOIN / 3) nExtraJoin = extraJoin( 2., event);
409 0 : if (iTry == 2 * NTRYJOIN / 3) nExtraJoin += extraJoin( 4., event);
410 :
411 : // After several failed tries gradually allow larger stop mass.
412 0 : if (iTry > NTRYJOIN - NSTOPMASS) stopMassNow *= FACSTOPMASS;
413 :
414 : // Set up flavours of two string ends, and reset other info.
415 0 : setStartEnds(idPos, idNeg, system);
416 0 : pRem = pSum;
417 :
418 : // Begin fragmentation loop, interleaved from the two ends.
419 : bool fromPos;
420 :
421 : // Variables used to help identifying baryons from junction splittings.
422 : bool usedPosJun = false, usedNegJun = false;
423 :
424 0 : for ( ; ; ) {
425 :
426 : // Take a step either from the positive or the negative end.
427 0 : fromPos = (rndmPtr->flat() < 0.5);
428 0 : StringEnd& nowEnd = (fromPos) ? posEnd : negEnd;
429 :
430 : // Construct trial hadron and check that energy remains.
431 0 : nowEnd.newHadron();
432 0 : if ( energyUsedUp(fromPos) ) break;
433 :
434 : // Construct kinematics of the new hadron and store it.
435 0 : Vec4 pHad = nowEnd.kinematicsHadron(system);
436 0 : int statusHad = (fromPos) ? 83 : 84;
437 :
438 : // Change status code if hadron from junction.
439 0 : if (abs(nowEnd.idHad) > 1000 && abs(nowEnd.idHad) < 10000) {
440 0 : if (fromPos && event[iPos].statusAbs() == 74 && !usedPosJun) {
441 : statusHad = 87;
442 : usedPosJun = true;
443 0 : }
444 0 : if (!fromPos && event[iNeg].statusAbs() == 74 && !usedNegJun) {
445 : statusHad = 88;
446 : usedNegJun = true;
447 0 : }
448 0 : if (!fromPos && hasJunction && !usedNegJun) {
449 : statusHad = 88;
450 : usedNegJun = true;
451 0 : }
452 : }
453 0 : hadrons.append( nowEnd.idHad, statusHad, iPos, iNeg,
454 0 : 0, 0, 0, 0, pHad, nowEnd.mHad);
455 0 : if (pHad.e() < 0.) break;
456 :
457 : // Update string end and remaining momentum.
458 0 : nowEnd.update();
459 0 : pRem -= pHad;
460 :
461 : // End of fragmentation loop.
462 0 : }
463 :
464 : // When done, join in the middle. If this works, then really done.
465 0 : if ( finalTwo(fromPos, event, usedPosJun, usedNegJun) ) break;
466 :
467 : // Else remove produced particles (except from first two junction legs)
468 : // and start all over.
469 0 : int newHadrons = hadrons.size() - junctionHadrons;
470 0 : hadrons.popBack(newHadrons);
471 0 : }
472 :
473 :
474 : // Junctions & extra joins: remove fictitious end, restore original partons.
475 0 : if (hasJunction) ++nExtraJoin;
476 0 : if (nExtraJoin > 0) {
477 0 : event.popBack(nExtraJoin);
478 0 : iParton = colConfig[iSub].iParton;
479 0 : }
480 :
481 : // Store the hadrons in the normal event record, ordered from one end.
482 0 : store(event);
483 :
484 : // Done.
485 0 : return true;
486 :
487 0 : }
488 :
489 : //--------------------------------------------------------------------------
490 :
491 : // Find region where to put first string break for closed gluon loop.
492 :
493 : vector<int> StringFragmentation::findFirstRegion(vector<int>& iPartonIn,
494 : Event& event) {
495 :
496 : // Evaluate mass-squared for all adjacent gluon pairs.
497 0 : vector<double> m2Pair;
498 : double m2Sum = 0.;
499 0 : int size = iPartonIn.size();
500 0 : for (int i = 0; i < size; ++i) {
501 0 : double m2Now = 0.5 * event[ iPartonIn[i] ].p()
502 0 : * event[ iPartonIn[(i + 1)%size] ].p();
503 0 : m2Pair.push_back(m2Now);
504 0 : m2Sum += m2Now;
505 0 : }
506 :
507 : // Pick breakup region with probability proportional to mass-squared.
508 0 : double m2Reg = m2Sum * rndmPtr->flat();
509 : int iReg = -1;
510 0 : do m2Reg -= m2Pair[++iReg];
511 0 : while (m2Reg > 0. && iReg < size - 1);
512 :
513 : // Create reordered parton list, with breakup string region duplicated.
514 0 : vector<int> iPartonOut;
515 0 : for (int i = 0; i < size + 2; ++i)
516 0 : iPartonOut.push_back( iPartonIn[(i + iReg)%size] );
517 :
518 : // Done.
519 : return iPartonOut;
520 :
521 0 : }
522 :
523 : //--------------------------------------------------------------------------
524 :
525 : // Set flavours and momentum position for initial string endpoints.
526 :
527 : void StringFragmentation::setStartEnds( int idPos, int idNeg,
528 : StringSystem systemNow) {
529 :
530 : // Variables characterizing string endpoints: defaults for open string.
531 : double px = 0.;
532 : double py = 0.;
533 : double Gamma = 0.;
534 : double xPosFromPos = 1.;
535 : double xNegFromPos = 0.;
536 : double xPosFromNeg = 0.;
537 : double xNegFromNeg = 1.;
538 :
539 : // For closed gluon loop need to pick an initial flavour.
540 0 : if (isClosed) {
541 : do {
542 0 : int idTry = flavSelPtr->pickLightQ();
543 0 : FlavContainer flavTry(idTry, 1);
544 0 : flavTry = flavSelPtr->pick( flavTry);
545 0 : flavTry = flavSelPtr->pick( flavTry);
546 0 : idPos = flavTry.id;
547 0 : idNeg = -idPos;
548 0 : } while (idPos == 0);
549 :
550 : // Also need pT and breakup vertex position in region.
551 0 : pair<double, double> pxy = pTSelPtr->pxy();
552 : px = pxy.first;
553 : py = pxy.second;
554 0 : double m2Region = systemNow.regionLowPos(0).w2;
555 0 : double m2Temp = min( CLOSEDM2MAX, CLOSEDM2FRAC * m2Region);
556 0 : do {
557 0 : double zTemp = zSelPtr->zFrag( idPos, idNeg, m2Temp);
558 0 : xPosFromPos = 1. - zTemp;
559 0 : xNegFromPos = m2Temp / (zTemp * m2Region);
560 0 : } while (xNegFromPos > 1.);
561 0 : Gamma = xPosFromPos * xNegFromPos * m2Region;
562 : xPosFromNeg = xPosFromPos;
563 : xNegFromNeg = xNegFromPos;
564 0 : }
565 :
566 : // Initialize two string endpoints.
567 0 : posEnd.setUp( true, iPos, idPos, systemNow.iMax, px, py,
568 : Gamma, xPosFromPos, xNegFromPos);
569 0 : negEnd.setUp( false, iNeg, idNeg, systemNow.iMax, -px, -py,
570 : Gamma, xPosFromNeg, xNegFromNeg);
571 :
572 : // For closed gluon loop can allow popcorn on one side but not both.
573 0 : if (isClosed) {
574 0 : flavSelPtr->assignPopQ(posEnd.flavOld);
575 0 : flavSelPtr->assignPopQ(negEnd.flavOld);
576 0 : if (rndmPtr->flat() < 0.5) posEnd.flavOld.nPop = 0;
577 0 : else negEnd.flavOld.nPop = 0;
578 0 : posEnd.flavOld.rank = 1;
579 0 : negEnd.flavOld.rank = 1;
580 0 : }
581 :
582 : // Done.
583 :
584 0 : }
585 :
586 : //--------------------------------------------------------------------------
587 :
588 : // Check remaining energy-momentum whether it is OK to continue.
589 :
590 : bool StringFragmentation::energyUsedUp(bool fromPos) {
591 :
592 : // If remaining negative energy then abort right away.
593 0 : if (pRem.e() < 0.) return true;
594 :
595 : // Calculate W2_minimum and done if remaining W2 is below it.
596 0 : double wMin = stopMassNow
597 0 : + particleDataPtr->constituentMass(posEnd.flavOld.id)
598 0 : + particleDataPtr->constituentMass(negEnd.flavOld.id);
599 0 : if (fromPos) wMin += stopNewFlav
600 0 : * particleDataPtr->constituentMass(posEnd.flavNew.id);
601 0 : else wMin += stopNewFlav
602 0 : * particleDataPtr->constituentMass(negEnd.flavNew.id);
603 0 : wMin *= 1. + (2. * rndmPtr->flat() - 1.) * stopSmear;
604 0 : w2Rem = pRem.m2Calc();
605 0 : if (w2Rem < pow2(wMin)) return true;
606 :
607 : // Else still enough energy left to continue iteration.
608 0 : return false;
609 :
610 0 : }
611 :
612 : //--------------------------------------------------------------------------
613 :
614 : // Produce the final two partons to complete the system.
615 :
616 : bool StringFragmentation::finalTwo(bool fromPos, Event& event, bool usedPosJun,
617 : bool usedNegJun) {
618 :
619 : // Check whether we went too far in p+-.
620 0 : if (pRem.e() < 0. || w2Rem < 0. || (hadrons.size() > 0
621 0 : && hadrons.back().e() < 0.) ) return false;
622 0 : if ( posEnd.iPosOld > negEnd.iPosOld || negEnd.iNegOld > posEnd.iNegOld)
623 0 : return false;
624 0 : if ( posEnd.iPosOld == negEnd.iPosOld && posEnd.xPosOld < negEnd.xPosOld)
625 0 : return false;
626 0 : if ( posEnd.iNegOld == negEnd.iNegOld && posEnd.xNegOld > negEnd.xNegOld)
627 0 : return false;
628 :
629 : // Construct the final hadron from the leftover flavours.
630 : // Impossible to join two diquarks. Also break if stuck for other reason.
631 0 : FlavContainer flav1 = (fromPos) ? posEnd.flavNew.anti() : posEnd.flavOld;
632 0 : FlavContainer flav2 = (fromPos) ? negEnd.flavOld : negEnd.flavNew.anti();
633 0 : if (flav1.isDiquark() && flav2.isDiquark()) return false;
634 : int idHad;
635 0 : for (int iTry = 0; iTry < NTRYFLAV; ++iTry) {
636 0 : idHad = flavSelPtr->combine( flav1, flav2);
637 0 : if (idHad != 0) break;
638 : }
639 0 : if (idHad == 0) return false;
640 :
641 : // Store the final particle and its new pT, and construct its mass.
642 0 : if (fromPos) {
643 0 : negEnd.idHad = idHad;
644 0 : negEnd.pxNew = -posEnd.pxNew;
645 0 : negEnd.pyNew = -posEnd.pyNew;
646 0 : negEnd.mHad = particleDataPtr->mSel(idHad);
647 0 : } else {
648 0 : posEnd.idHad = idHad;
649 0 : posEnd.pxNew = -negEnd.pxNew;
650 0 : posEnd.pyNew = -negEnd.pyNew;
651 0 : posEnd.mHad = particleDataPtr->mSel(idHad);
652 : }
653 :
654 : // String region in which to do the joining.
655 0 : StringRegion region = finalRegion();
656 0 : if (region.isEmpty) return false;
657 :
658 : // Project remaining momentum along longitudinal and transverse directions.
659 0 : region.project( pRem);
660 0 : double pxRem = region.px() - posEnd.pxOld - negEnd.pxOld;
661 0 : double pyRem = region.py() - posEnd.pyOld - negEnd.pyOld;
662 0 : double xPosRem = region.xPos();
663 0 : double xNegRem = region.xNeg();
664 :
665 : // Share extra pT kick evenly between final two hadrons.
666 0 : posEnd.pxOld += 0.5 * pxRem;
667 0 : posEnd.pyOld += 0.5 * pyRem;
668 0 : negEnd.pxOld += 0.5 * pxRem;
669 0 : negEnd.pyOld += 0.5 * pyRem;
670 :
671 : // Construct new pT and mT of the final two particles.
672 0 : posEnd.pxHad = posEnd.pxOld + posEnd.pxNew;
673 0 : posEnd.pyHad = posEnd.pyOld + posEnd.pyNew;
674 0 : posEnd.mT2Had = pow2(posEnd.mHad) + pow2(posEnd.pxHad)
675 0 : + pow2(posEnd.pyHad);
676 0 : negEnd.pxHad = negEnd.pxOld + negEnd.pxNew;
677 0 : negEnd.pyHad = negEnd.pyOld + negEnd.pyNew;
678 0 : negEnd.mT2Had = pow2(negEnd.mHad) + pow2(negEnd.pxHad)
679 0 : + pow2(negEnd.pyHad);
680 :
681 : // Construct remaining system transverse mass.
682 0 : double wT2Rem = w2Rem + pow2( posEnd.pxHad + negEnd.pxHad)
683 0 : + pow2( posEnd.pyHad + negEnd.pyHad);
684 :
685 : // Check that kinematics possible.
686 0 : if ( sqrt(wT2Rem) < sqrt(posEnd.mT2Had) + sqrt(negEnd.mT2Had) )
687 0 : return false;
688 0 : double lambda2 = pow2( wT2Rem - posEnd.mT2Had - negEnd.mT2Had)
689 0 : - 4. * posEnd.mT2Had * negEnd.mT2Had;
690 0 : if (lambda2 <= 0.) return false;
691 :
692 : // Construct kinematics, as viewed in the transverse rest frame.
693 0 : double lambda = sqrt(lambda2);
694 0 : double probReverse = 1. / (1. + exp( min( EXPMAX, bLund * lambda) ) );
695 0 : double xpzPos = 0.5 * lambda/ wT2Rem;
696 0 : if (probReverse > rndmPtr->flat()) xpzPos = -xpzPos;
697 0 : double xmDiff = (posEnd.mT2Had - negEnd.mT2Had) / wT2Rem;
698 0 : double xePos = 0.5 * (1. + xmDiff);
699 0 : double xeNeg = 0.5 * (1. - xmDiff );
700 :
701 : // Translate this into kinematics in the string frame.
702 0 : Vec4 pHadPos = region.pHad( (xePos + xpzPos) * xPosRem,
703 0 : (xePos - xpzPos) * xNegRem, posEnd.pxHad, posEnd.pyHad);
704 0 : Vec4 pHadNeg = region.pHad( (xeNeg - xpzPos) * xPosRem,
705 0 : (xeNeg + xpzPos) * xNegRem, negEnd.pxHad, negEnd.pyHad);
706 :
707 : // Update status codes for junction baryons.
708 : int statusHadPos = 83;
709 : int statusHadNeg = 84;
710 :
711 0 : if (fromPos) {
712 0 : if (abs(posEnd.idHad) > 1000 && abs(posEnd.idHad) < 10000) {
713 0 : if (event[iPos].statusAbs() == 74 && !usedPosJun) {
714 : statusHadPos = 87;
715 : usedPosJun = true;
716 0 : }
717 : }
718 0 : if (abs(idHad) > 1000 && abs(idHad) < 10000) {
719 0 : if ((!usedNegJun && (event[iNeg].statusAbs() == 74 || hasJunction))
720 0 : || (!usedPosJun && event[iPos].statusAbs() == 74)) {
721 : statusHadNeg = 88;
722 0 : }
723 : }
724 : } else {
725 0 : if (abs(negEnd.idHad) > 1000 && abs(negEnd.idHad) < 10000) {
726 0 : if (!usedNegJun && (event[iNeg].statusAbs() == 74 || hasJunction)) {
727 : statusHadNeg = 87;
728 : usedNegJun = true;
729 0 : }
730 : }
731 0 : if (abs(idHad) > 1000 && abs(idHad) < 10000) {
732 0 : if ((!usedNegJun && (event[iNeg].statusAbs() == 74 || hasJunction))
733 0 : || (!usedPosJun && event[iPos].statusAbs() == 74)) {
734 : statusHadPos = 88;
735 0 : }
736 : }
737 : }
738 :
739 : // Add produced particles to the event record.
740 0 : hadrons.append( posEnd.idHad, statusHadPos, posEnd.iEnd, negEnd.iEnd,
741 0 : 0, 0, 0, 0, pHadPos, posEnd.mHad);
742 0 : hadrons.append( negEnd.idHad, statusHadNeg, posEnd.iEnd, negEnd.iEnd,
743 0 : 0, 0, 0, 0, pHadNeg, negEnd.mHad);
744 :
745 : // It worked.
746 : return true;
747 :
748 0 : }
749 :
750 : //--------------------------------------------------------------------------
751 :
752 : // Construct a special joining region for the final two hadrons.
753 :
754 : StringRegion StringFragmentation::finalRegion() {
755 :
756 : // Simple case when both string ends are in the same region.
757 0 : if (posEnd.iPosOld == negEnd.iPosOld && posEnd.iNegOld == negEnd.iNegOld)
758 0 : return system.region( posEnd.iPosOld, posEnd.iNegOld);
759 :
760 : // Start out with empty region. (Empty used for error returns.)
761 0 : StringRegion region;
762 :
763 : // Add up all remaining p+.
764 0 : Vec4 pPosJoin;
765 0 : if ( posEnd.iPosOld == negEnd.iPosOld) {
766 0 : double xPosJoin = posEnd.xPosOld - negEnd.xPosOld;
767 0 : if (xPosJoin < 0.) return region;
768 0 : pPosJoin = system.regionLowPos(posEnd.iPosOld).pHad( xPosJoin, 0., 0., 0.);
769 0 : } else {
770 0 : for (int iPosNow = posEnd.iPosOld; iPosNow <= negEnd.iPosOld; ++iPosNow) {
771 0 : if (iPosNow == posEnd.iPosOld) pPosJoin
772 0 : += system.regionLowPos(iPosNow).pHad( posEnd.xPosOld, 0., 0., 0.);
773 0 : else if (iPosNow == negEnd.iPosOld) pPosJoin
774 0 : += system.regionLowPos(iPosNow).pHad( 1. - negEnd.xPosOld, 0., 0., 0.);
775 0 : else pPosJoin += system.regionLowPos(iPosNow).pHad( 1., 0., 0., 0.);
776 : }
777 : }
778 :
779 : // Add up all remaining p-.
780 0 : Vec4 pNegJoin;
781 0 : if ( negEnd.iNegOld == posEnd.iNegOld) {
782 0 : double xNegJoin = negEnd.xNegOld - posEnd.xNegOld;
783 0 : if (xNegJoin < 0.) return region;
784 0 : pNegJoin = system.regionLowNeg(negEnd.iNegOld).pHad( 0., xNegJoin, 0., 0.);
785 0 : } else {
786 0 : for (int iNegNow = negEnd.iNegOld; iNegNow <= posEnd.iNegOld; ++iNegNow) {
787 0 : if (iNegNow == negEnd.iNegOld) pNegJoin
788 0 : += system.regionLowNeg(iNegNow).pHad( 0., negEnd.xNegOld, 0., 0.);
789 0 : else if (iNegNow == posEnd.iNegOld) pNegJoin
790 0 : += system.regionLowNeg(iNegNow).pHad( 0., 1. - posEnd.xNegOld, 0., 0.);
791 0 : else pNegJoin += system.regionLowNeg(iNegNow).pHad( 0., 1., 0., 0.);
792 : }
793 : }
794 :
795 : // For a closed gluon loop pPosJoin == pNegJoin and the above does not work.
796 : // So reshuffle; "perfect" for g g systems, OK in general.
797 0 : Vec4 pTest = pPosJoin - pNegJoin;
798 0 : if ( abs(pTest.px()) + abs(pTest.py()) + abs(pTest.pz()) + abs(pTest.e())
799 0 : < MATCHPOSNEG * (pPosJoin.e() + pNegJoin.e()) ) {
800 0 : Vec4 delta
801 0 : = system.regionLowPos(posEnd.iPosOld + 1).pHad( 1., 0., 0., 0.)
802 0 : - system.regionLowNeg(negEnd.iNegOld + 1).pHad( 0., 1., 0., 0.);
803 : // If reshuffle did not help then pick random axis to break tie.
804 : // (Needed for low-mass q-g-qbar with q-qbar perfectly parallel.)
805 0 : if ( abs(delta.px()) + abs(delta.py()) + abs(delta.pz()) + abs(delta.e())
806 0 : < MATCHPOSNEG * (pPosJoin.e() + pNegJoin.e()) ) {
807 0 : double cthe = 2. * rndmPtr->flat() - 1.;
808 0 : double sthe = sqrtpos(1. - cthe * cthe);
809 0 : double phi = 2. * M_PI * rndmPtr->flat();
810 0 : delta = 0.5 * min( pPosJoin.e(), pNegJoin.e())
811 0 : * Vec4( sthe * sin(phi), sthe * cos(phi), cthe, 0.);
812 0 : infoPtr->errorMsg("Warning in StringFragmentation::finalRegion: "
813 : "random axis needed to break tie");
814 0 : }
815 0 : pPosJoin -= delta;
816 0 : pNegJoin += delta;
817 0 : }
818 :
819 : // Construct a new region from remaining p+ and p-.
820 0 : region.setUp( pPosJoin, pNegJoin);
821 0 : if (region.isEmpty) return region;
822 :
823 : // Project the existing pTold vectors onto the new directions.
824 0 : Vec4 pTposOld = system.region( posEnd.iPosOld, posEnd.iNegOld).pHad(
825 0 : 0., 0., posEnd.pxOld, posEnd.pyOld);
826 0 : region.project( pTposOld);
827 0 : posEnd.pxOld = region.px();
828 0 : posEnd.pyOld = region.py();
829 0 : Vec4 pTnegOld = system.region( negEnd.iPosOld, negEnd.iNegOld).pHad(
830 0 : 0., 0., negEnd.pxOld, negEnd.pyOld);
831 0 : region.project( pTnegOld);
832 0 : negEnd.pxOld = region.px();
833 0 : negEnd.pyOld = region.py();
834 :
835 : // Done.
836 0 : return region;
837 :
838 0 : }
839 :
840 : //--------------------------------------------------------------------------
841 :
842 : // Store the hadrons in the normal event record, ordered from one end.
843 :
844 : void StringFragmentation::store(Event& event) {
845 :
846 : // Starting position.
847 0 : int iFirst = event.size();
848 :
849 : // Copy straight over from first two junction legs.
850 0 : if (hasJunction) {
851 0 : for (int i = 0; i < hadrons.size(); ++i)
852 0 : if (hadrons[i].status() == 85 || hadrons[i].status() == 86)
853 0 : event.append( hadrons[i] );
854 0 : }
855 :
856 : // Loop downwards, copying all from the positive end.
857 0 : for (int i = 0; i < hadrons.size(); ++i)
858 0 : if (hadrons[i].status() == 83 || hadrons[i].status() == 87)
859 0 : event.append( hadrons[i] );
860 :
861 : // Loop upwards, copying all from the negative end.
862 0 : for (int i = hadrons.size() - 1; i >= 0 ; --i)
863 0 : if (hadrons[i].status() == 84 || hadrons[i].status() == 88)
864 0 : event.append( hadrons[i] );
865 :
866 0 : int iLast = event.size() - 1;
867 :
868 : // Set decay vertex when this is displaced.
869 0 : if (event[posEnd.iEnd].hasVertex()) {
870 0 : Vec4 vDec = event[posEnd.iEnd].vDec();
871 0 : for (int i = iFirst; i <= iLast; ++i) event[i].vProd( vDec );
872 0 : }
873 :
874 : // Set lifetime of hadrons.
875 0 : for (int i = iFirst; i <= iLast; ++i)
876 0 : event[i].tau( event[i].tau0() * rndmPtr->exp() );
877 :
878 : // Mark original partons as hadronized and set their daughter range.
879 0 : for (int i = 0; i < int(iParton.size()); ++i)
880 0 : if (iParton[i] >= 0) {
881 0 : event[ iParton[i] ].statusNeg();
882 0 : event[ iParton[i] ].daughters(iFirst, iLast);
883 0 : }
884 :
885 0 : }
886 :
887 : //--------------------------------------------------------------------------
888 :
889 : // Fragment off two of the string legs in to a junction.
890 :
891 : bool StringFragmentation::fragmentToJunction(Event& event) {
892 :
893 : // Identify range of partons on the three legs.
894 : // (Each leg begins with an iParton[i] = -(10 + 10*junctionNumber + leg),
895 : // and partons then appear ordered from the junction outwards.)
896 0 : int legBeg[3] = { 0, 0, 0};
897 0 : int legEnd[3] = { 0, 0, 0};
898 : int leg = -1;
899 : // PS (4/10/2011) Protect against invalid systems
900 0 : if (iParton[0] > 0) {
901 0 : infoPtr->errorMsg("Error in StringFragmentation::fragment"
902 : "ToJunction: iParton[0] not a valid junctionNumber");
903 0 : return false;
904 : }
905 0 : for (int i = 0; i < int(iParton.size()); ++i) {
906 0 : if (iParton[i] < 0) {
907 0 : if (leg == 2) {
908 0 : infoPtr->errorMsg("Error in StringFragmentation::fragment"
909 : "ToJunction: unprocessed multi-junction system");
910 0 : return false;
911 : }
912 0 : legBeg[++leg] = i + 1;
913 0 : }
914 0 : else legEnd[leg] = i;
915 : }
916 :
917 : // Iterate from system rest frame towards the junction rest frame (JRF).
918 0 : RotBstMatrix MtoJRF, Mstep;
919 0 : MtoJRF.bstback(pSum);
920 0 : Vec4 pWTinJRF[3];
921 : int iter = 0;
922 : double errInCM = 0.;
923 0 : do {
924 0 : ++iter;
925 :
926 : // Find weighted sum of momenta on the three sides of the junction.
927 0 : for (leg = 0; leg < 3; ++ leg) {
928 0 : pWTinJRF[leg] = 0.;
929 : double eWeight = 0.;
930 0 : for (int i = legBeg[leg]; i <= legEnd[leg]; ++i) {
931 0 : Vec4 pTemp = event[ iParton[i] ].p();
932 0 : pTemp.rotbst(MtoJRF);
933 0 : pWTinJRF[leg] += pTemp * exp(-eWeight);
934 0 : eWeight += pTemp.e() / eNormJunction;
935 0 : if (eWeight > EJNWEIGHTMAX) break;
936 0 : }
937 : }
938 :
939 : // Store original deviation from 120 degree topology.
940 0 : if (iter == 1) errInCM = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5)
941 0 : + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5)
942 0 : + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5);
943 :
944 : // Check that no two legs has unreasonably small invariant mass.
945 0 : if ( (pWTinJRF[0] + pWTinJRF[1]).m2Calc() < M2MINJRF
946 0 : || (pWTinJRF[0] + pWTinJRF[2]).m2Calc() < M2MINJRF
947 0 : || (pWTinJRF[1] + pWTinJRF[2]).m2Calc() < M2MINJRF ) {
948 0 : infoPtr->errorMsg("Error in StringFragmentation::fragmentToJunction: "
949 : "too small leg-leg invariant mass");
950 0 : return false;
951 : }
952 :
953 : // Find new JRF from the set of weighted momenta.
954 0 : Mstep = junctionRestFrame( pWTinJRF[0], pWTinJRF[1], pWTinJRF[2]);
955 : // Fortran code will not take full step after the first few
956 : // iterations. How implement this in terms of an M matrix??
957 0 : MtoJRF.rotbst( Mstep );
958 0 : } while (iter < 3 || (Mstep.deviation() > CONVJNREST && iter < NTRYJNREST) );
959 :
960 : // If final deviation from 120 degrees is bigger than in CM then revert.
961 0 : double errInJRF = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5)
962 0 : + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5)
963 0 : + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5);
964 0 : if (errInJRF > errInCM + CONVJNREST) {
965 0 : infoPtr->errorMsg("Warning in StringFragmentation::fragmentTo"
966 : "Junction: bad convergence junction rest frame");
967 0 : MtoJRF.reset();
968 0 : MtoJRF.bstback(pSum);
969 0 : }
970 :
971 : // Opposite operation: boost from JRF to original system.
972 0 : RotBstMatrix MfromJRF = MtoJRF;
973 0 : MfromJRF.invert();
974 :
975 : // Sum leg four-momenta in original frame and in JRF.
976 0 : Vec4 pInLeg[3], pInJRF[3];
977 0 : for (leg = 0; leg < 3; ++leg) {
978 0 : pInLeg[leg] = 0.;
979 0 : for (int i = legBeg[leg]; i <= legEnd[leg]; ++i)
980 0 : pInLeg[leg] += event[ iParton[i] ].p();
981 0 : pInJRF[leg] = pInLeg[leg];
982 0 : pInJRF[leg].rotbst(MtoJRF);
983 : }
984 :
985 : // Pick the two legs with lowest energy in JRF.
986 0 : int legMin = (pInJRF[0].e() < pInJRF[1].e()) ? 0 : 1;
987 0 : int legMax = 1 - legMin;
988 0 : if (pInJRF[2].e() < min(pInJRF[0].e(), pInJRF[1].e()) ) legMin = 2;
989 0 : else if (pInJRF[2].e() > max(pInJRF[0].e(), pInJRF[1].e()) ) legMax = 2;
990 0 : int legMid = 3 - legMin - legMax;
991 :
992 : // Save info on which status codes belong with the three legs.
993 0 : int iJunction = (-iParton[0]) / 10 - 1;
994 0 : event.statusJunction( iJunction, legMin, 85);
995 0 : event.statusJunction( iJunction, legMid, 86);
996 0 : event.statusJunction( iJunction, legMax, 83);
997 :
998 : // Temporarily copy the partons on the low-energy legs, into the JRF,
999 : // in reverse order, so (anti)quark leg end first.
1000 0 : vector<int> iPartonMin;
1001 0 : for (int i = legEnd[legMin]; i >= legBeg[legMin]; --i) {
1002 0 : int iNew = event.append( event[ iParton[i] ] );
1003 0 : event[iNew].rotbst(MtoJRF);
1004 0 : iPartonMin.push_back( iNew );
1005 0 : }
1006 0 : vector<int> iPartonMid;
1007 0 : for (int i = legEnd[legMid]; i >= legBeg[legMid]; --i) {
1008 0 : int iNew = event.append( event[ iParton[i] ] );
1009 0 : event[iNew].rotbst(MtoJRF);
1010 0 : iPartonMid.push_back( iNew );
1011 0 : }
1012 :
1013 : // Find final weighted sum of momenta on each of the two legs.
1014 : double eWeight = 0.;
1015 0 : pWTinJRF[legMin] = 0.;
1016 0 : for (int i = iPartonMin.size() - 1; i >= 0; --i) {
1017 0 : pWTinJRF[legMin] += event[ iPartonMin[i] ].p() * exp(-eWeight);
1018 0 : eWeight += event[ iPartonMin[i] ].e() / eNormJunction;
1019 0 : if (eWeight > EJNWEIGHTMAX) break;
1020 : }
1021 : eWeight = 0.;
1022 0 : pWTinJRF[legMid] = 0.;
1023 0 : for (int i = iPartonMid.size() - 1; i >= 0; --i) {
1024 0 : pWTinJRF[legMid] += event[ iPartonMid[i] ].p() * exp(-eWeight);
1025 0 : eWeight += event[ iPartonMid[i] ].e() / eNormJunction;
1026 0 : if (eWeight > EJNWEIGHTMAX) break;
1027 : }
1028 :
1029 : // Define fictitious opposing partons in JRF and store as string ends.
1030 0 : Vec4 pOppose = pWTinJRF[legMin];
1031 0 : pOppose.flip3();
1032 0 : int idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1;
1033 0 : if (event[ iPartonMin[0] ].col() > 0) idOppose = -idOppose;
1034 0 : int iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0,
1035 0 : pOppose, 0.);
1036 0 : iPartonMin.push_back( iOppose);
1037 0 : pOppose = pWTinJRF[legMid];
1038 0 : pOppose.flip3();
1039 0 : idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1;
1040 0 : if (event[ iPartonMid[0] ].col() > 0) idOppose = -idOppose;
1041 0 : iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0,
1042 0 : pOppose, 0.);
1043 0 : iPartonMid.push_back( iOppose);
1044 :
1045 : // Set up kinematics of string evolution in low-energy temporary systems.
1046 0 : systemMin.setUp(iPartonMin, event);
1047 0 : systemMid.setUp(iPartonMid, event);
1048 :
1049 : // Outer fallback loop, when too little energy left for third leg.
1050 : int idMin = 0;
1051 : int idMid = 0;
1052 0 : Vec4 pDiquark;
1053 0 : for ( int iTryOuter = 0; ; ++iTryOuter) {
1054 :
1055 : // Middle fallback loop, when much unused energy in leg remnants.
1056 0 : double eLeftMin = 0.;
1057 0 : double eLeftMid = 0.;
1058 0 : for ( int iTryMiddle = 0; ; ++iTryMiddle) {
1059 :
1060 : // Loop over the two lowest-energy legs.
1061 0 : for (int legLoop = 0; legLoop < 2; ++ legLoop) {
1062 0 : int legNow = (legLoop == 0) ? legMin : legMid;
1063 :
1064 : // Read in properties specific to this leg.
1065 0 : StringSystem& systemNow = (legLoop == 0) ? systemMin : systemMid;
1066 0 : int idPos = (legLoop == 0) ? event[ iPartonMin[0] ].id()
1067 0 : : event[ iPartonMid[0] ].id();
1068 0 : idOppose = (legLoop == 0) ? event[ iPartonMin.back() ].id()
1069 0 : : event[ iPartonMid.back() ].id();
1070 0 : double eInJRF = pInJRF[legNow].e();
1071 0 : int statusHad = (legLoop == 0) ? 85 : 86;
1072 :
1073 : // Inner fallback loop, when a diquark comes in to junction.
1074 : double eUsed = 0.;
1075 0 : for ( int iTryInner = 0; ; ++iTryInner) {
1076 0 : if (iTryInner > 2 * NTRYJNMATCH) {
1077 0 : infoPtr->errorMsg("Error in StringFragmentation::fragment"
1078 : "ToJunction: caught in junction flavour loop");
1079 0 : event.popBack( iPartonMin.size() + iPartonMid.size() );
1080 0 : return false;
1081 : }
1082 0 : bool needBaryon = (abs(idPos) > 10 && iTryInner > NTRYJNMATCH);
1083 0 : double eExtra = (iTryInner > NTRYJNMATCH) ? EEXTRAJNMATCH : 0.;
1084 :
1085 : // Set up two string ends, and begin fragmentation loop.
1086 0 : setStartEnds(idPos, idOppose, systemNow);
1087 : eUsed = 0.;
1088 : int nHadrons = 0;
1089 : bool noNegE = true;
1090 0 : for ( ; ; ++nHadrons) {
1091 :
1092 : // Construct trial hadron from positive end.
1093 0 : posEnd.newHadron();
1094 0 : Vec4 pHad = posEnd.kinematicsHadron(systemNow);
1095 :
1096 : // Negative energy signals failure in construction.
1097 0 : if (pHad.e() < 0. ) { noNegE = false; break; }
1098 :
1099 : // Break if passed system midpoint ( = junction) in energy.
1100 : // Exceptions: small systems, and/or with diquark end.
1101 : bool delayedBreak = false;
1102 0 : if (eUsed + pHad.e() + eExtra > eInJRF) {
1103 0 : if (nHadrons > 0 || !needBaryon) break;
1104 : delayedBreak = true;
1105 0 : }
1106 :
1107 : // Else construct kinematics of the new hadron and store it.
1108 0 : hadrons.append( posEnd.idHad, statusHad, iPos, iNeg,
1109 0 : 0, 0, 0, 0, pHad, posEnd.mHad);
1110 :
1111 : // Update string end and remaining momentum.
1112 0 : posEnd.update();
1113 0 : eUsed += pHad.e();
1114 :
1115 : // Delayed break in small systems, and/or with diquark end.
1116 0 : if (delayedBreak) {
1117 0 : ++nHadrons;
1118 0 : break;
1119 : }
1120 0 : }
1121 :
1122 : // Possible to produce zero hadrons if the endpoint is not a diquark.
1123 0 : if (iTryInner > NTRYJNMATCH && !noNegE && nHadrons == 0 &&
1124 0 : abs(idPos) < 10) break;
1125 :
1126 : // End of fragmentation loop. Inner loopback if ends on a diquark.
1127 0 : if ( noNegE && abs(posEnd.flavOld.id) < 10 ) break;
1128 0 : hadrons.popBack(nHadrons);
1129 0 : }
1130 :
1131 : // End of one-leg fragmentation. Store end quark and remnant energy.
1132 0 : if (legNow == legMin) {
1133 0 : idMin = posEnd.flavOld.id;
1134 0 : eLeftMin = eInJRF - eUsed;
1135 0 : } else {
1136 : idMid = posEnd.flavOld.id;
1137 0 : eLeftMid = eInJRF - eUsed;
1138 : }
1139 0 : }
1140 :
1141 : // End of both-leg fragmentation.
1142 : // Middle loopback if too much energy left.
1143 0 : double eTrial = eBothLeftJunction + rndmPtr->flat() * eMaxLeftJunction;
1144 0 : if (iTryMiddle > NTRYJNMATCH
1145 0 : || ( min( eLeftMin, eLeftMid) < eBothLeftJunction
1146 0 : && max( eLeftMin, eLeftMid) < eTrial ) ) break;
1147 0 : hadrons.clear();
1148 0 : }
1149 :
1150 : // Boost hadrons away from the JRF to the original frame.
1151 0 : for (int i = 0; i < hadrons.size(); ++i) {
1152 0 : hadrons[i].rotbst(MfromJRF);
1153 : // Recalculate energy to compensate for numerical precision loss
1154 : // in iterative calculation of MfromJRF.
1155 0 : hadrons[i].e( hadrons[i].eCalc() );
1156 0 : pJunctionHadrons += hadrons[i].p();
1157 : }
1158 :
1159 : // Outer loopback if combined diquark mass too negative
1160 : // or too little energy left in third leg.
1161 0 : pDiquark = pInLeg[legMin] + pInLeg[legMid] - pJunctionHadrons;
1162 0 : double m2Left = m2( pInLeg[legMax], pDiquark);
1163 0 : if (iTryOuter > NTRYJNMATCH || (pDiquark.mCalc() > MDIQUARKMIN
1164 0 : && m2Left > eMinLeftJunction * pInLeg[legMax].e()) ) break;
1165 0 : hadrons.clear();
1166 0 : pJunctionHadrons = 0.;
1167 0 : }
1168 :
1169 : // Now found solution; no more loopback. Remove temporary parton copies.
1170 0 : event.popBack( iPartonMin.size() + iPartonMid.size() );
1171 :
1172 : // Construct and store an effective diquark string end from the
1173 : // two remnant quark ends, for temporary usage.
1174 0 : int idDiquark = flavSelPtr->makeDiquark( idMin, idMid);
1175 0 : double mDiquark = pDiquark.mCalc();
1176 0 : int iDiquark = event.append( idDiquark, 78, 0, 0, 0, 0, 0, 0,
1177 0 : pDiquark, mDiquark);
1178 :
1179 : // Find the partons on the last leg, again in reverse order.
1180 0 : vector<int> iPartonMax;
1181 0 : for (int i = legEnd[legMax]; i >= legBeg[legMax]; --i)
1182 0 : iPartonMax.push_back( iParton[i] );
1183 0 : iPartonMax.push_back( iDiquark );
1184 :
1185 : // Recluster gluons nearby to diquark end when taken too much energy.
1186 0 : int iPsize = iPartonMax.size();
1187 0 : double m0Diquark = event[iDiquark].m0();
1188 0 : while (iPsize > 2) {
1189 0 : Vec4 pGluNear = event[ iPartonMax[iPsize - 2] ].p();
1190 0 : if ( (pDiquark + 0.5 * pGluNear).mCalc() > m0Diquark + mJoin ) break;
1191 0 : pDiquark += pGluNear;
1192 0 : event[iDiquark].p( pDiquark );
1193 0 : event[iDiquark].m( pDiquark.mCalc() );
1194 0 : iPartonMax.pop_back();
1195 0 : --iPsize;
1196 0 : iPartonMax[iPsize - 1] = iDiquark;
1197 0 : }
1198 :
1199 : // Modify parton list to remaining leg + remnant of the first two.
1200 0 : iParton = iPartonMax;
1201 :
1202 : // Done.
1203 : return true;
1204 0 : }
1205 :
1206 : //--------------------------------------------------------------------------
1207 :
1208 : // Find the boost matrix to the rest frame of a junction,
1209 : // given the three respective endpoint four-momenta.
1210 :
1211 : RotBstMatrix StringFragmentation::junctionRestFrame(Vec4& p0, Vec4& p1,
1212 : Vec4& p2) {
1213 :
1214 : // Calculate masses and other invariants.
1215 0 : Vec4 pSumJun = p0 + p1 + p2;
1216 0 : double sHat = pSumJun.m2Calc();
1217 0 : double pp[3][3];
1218 0 : pp[0][0] = p0.m2Calc();
1219 0 : pp[1][1] = p1.m2Calc();
1220 0 : pp[2][2] = p2.m2Calc();
1221 0 : pp[0][1] = pp[1][0] = p0 * p1;
1222 0 : pp[0][2] = pp[2][0] = p0 * p2;
1223 0 : pp[1][2] = pp[2][1] = p1 * p2;
1224 :
1225 : // Requirement (eiMax)_j = pi*pj/mj < (eiMax)_k = pi*pk/mk, used below,
1226 : // here rewritten as pi*pj * mk < pi*pk * mj and squared.
1227 0 : double eMax01 = pow2(pp[0][1]) * pp[2][2];
1228 0 : double eMax02 = pow2(pp[0][2]) * pp[1][1];
1229 0 : double eMax12 = pow2(pp[1][2]) * pp[0][0];
1230 :
1231 : // Initially pick i to be the most massive parton. but allow other tries.
1232 0 : int i = (pp[1][1] > pp[0][0]) ? 1 : 0;
1233 0 : if (pp[2][2] > max(pp[0][0], pp[1][1])) i = 2;
1234 : int j, k;
1235 : double ei = 0.;
1236 : double ej = 0.;
1237 : double ek = 0.;
1238 0 : for (int iTry = 0; iTry < 3; ++iTry) {
1239 :
1240 : // Pick j to give minimal eiMax, and k the third vector.
1241 0 : if (i == 0) j = (eMax02 < eMax01) ? 2 : 1;
1242 0 : else if (i == 1) j = (eMax12 < eMax01) ? 2 : 0;
1243 0 : else j = (eMax12 < eMax02) ? 1 : 0;
1244 0 : k = 3 - i - j;
1245 :
1246 : // Alternative names according to i, j, k conventions.
1247 0 : double m2i = pp[i][i];
1248 0 : double m2j = pp[j][j];
1249 0 : double m2k = pp[k][k];
1250 0 : double pipj = pp[i][j];
1251 0 : double pipk = pp[i][k];
1252 0 : double pjpk = pp[j][k];
1253 :
1254 : // avoid too small or negative values
1255 0 : if (pipj < 1e-6)
1256 0 : pipj = 1e-6;
1257 0 : if (pipk < 1e-6)
1258 0 : pipk = 1e-6;
1259 0 : if (pjpk < 1e-6)
1260 0 : pjpk = 1e-6;
1261 :
1262 : // Trivial to find new parton energies if all three partons are massless.
1263 0 : if (m2i < M2MAXJRF) {
1264 0 : ei = sqrt( 2. * pipk * pipj / (3. * pjpk) );
1265 0 : ej = sqrt( 2. * pjpk * pipj / (3. * pipk) );
1266 0 : ek = sqrt( 2. * pipk * pjpk / (3. * pipj) );
1267 :
1268 : // Else find three-momentum range for parton i and values at extremes.
1269 0 : } else {
1270 : // Minimum when i is at rest.
1271 : double piMin = 0.;
1272 0 : double eiMin = sqrt(m2i);
1273 0 : double ejMin = pipj / eiMin;
1274 0 : double ekMin = pipk / eiMin;
1275 0 : double pjMin = sqrtpos( ejMin*ejMin - m2j );
1276 0 : double pkMin = sqrtpos( ekMin*ekMin - m2k );
1277 0 : double fMin = ejMin * ekMin + 0.5 * pjMin * pkMin - pjpk;
1278 : // Maximum estimated when j + k is at rest, alternatively j at rest.
1279 0 : double eiMax = (pipj + pipk) / sqrt(max(0., m2j) + max(0., m2k) + 2. * pjpk);
1280 0 : if (m2j > M2MAXJRF) eiMax = min( eiMax, pipj / sqrt(m2j) );
1281 0 : double piMax = sqrtpos( eiMax*eiMax - m2i );
1282 0 : double temp = eiMax*eiMax - 0.25 *piMax*piMax;
1283 0 : double pjMax = (eiMax * sqrtpos( pipj*pipj - m2j * temp )
1284 0 : - 0.5 * piMax * pipj) / temp;
1285 0 : double pkMax = (eiMax * sqrtpos( pipk*pipk - m2k * temp )
1286 0 : - 0.5 * piMax * pipk) / temp;
1287 : double ejMax = 1e-06;
1288 0 : if(pjMax*pjMax + m2j > 1e-06) ejMax = sqrt(pjMax*pjMax + m2j);
1289 : double ekMax = 1e-06;
1290 0 : if(pkMax*pkMax + m2k > 1e-06) ekMax = sqrt(pkMax*pkMax + m2k);
1291 0 : double fMax = ejMax * ekMax + 0.5 * pjMax * pkMax - pjpk;
1292 :
1293 : // If unexpected values at upper endpoint then pick another parton.
1294 0 : if (fMax > 0.) {
1295 0 : int iPrel = (i + 1)%3;
1296 0 : if (pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;}
1297 0 : ++iTry;
1298 0 : iPrel = (i + 2)%3;
1299 0 : if (iTry < 3 && pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;}
1300 0 : }
1301 :
1302 : // Start binary + linear search to find solution inside range.
1303 : int iterMin = 0;
1304 : int iterMax = 0;
1305 0 : double pi = 0.5 * (piMin + piMax);
1306 0 : for (int iter = 0; iter < NTRYJRFEQ; ++iter) {
1307 :
1308 : // Derive momentum of other two partons and distance to root.
1309 0 : ei = sqrt(pi*pi + m2i);
1310 0 : temp = ei*ei - 0.25 * pi*pi;
1311 0 : double pj = (ei * sqrtpos( pipj*pipj - m2j * temp )
1312 0 : - 0.5 * pi * pipj) / temp;
1313 0 : double pk = (ei * sqrtpos( pipk*pipk - m2k * temp )
1314 0 : - 0.5 * pi * pipk) / temp;
1315 0 : ej = sqrt(pj*pj + max(0., m2j));
1316 0 : ek = sqrt(pk*pk + max(0., m2k));
1317 0 : double fNow = ej * ek + 0.5 * pj * pk - pjpk;
1318 :
1319 : // Replace lower or upper bound by new value.
1320 0 : if (fNow > 0.) { ++iterMin; piMin = pi; fMin = fNow;}
1321 0 : else {++iterMax; piMax = pi; fMax = fNow;}
1322 :
1323 : // Pick next i momentum to explore, hopefully closer to root.
1324 0 : if (2 * iter < NTRYJRFEQ
1325 0 : && (iterMin < 2 || iterMax < 2 || 4 * iter < NTRYJRFEQ))
1326 0 : { pi = 0.5 * (piMin + piMax); continue;}
1327 0 : if (fMin < 0. || fMax > 0. || abs(fNow) < CONVJRFEQ * sHat) break;
1328 0 : pi = piMin + (piMax - piMin) * fMin / (fMin - fMax);
1329 0 : }
1330 :
1331 : // If arrived here then either succeeded or exhausted possibilities.
1332 0 : } break;
1333 0 : }
1334 :
1335 : // Now we know the energies in the junction rest frame.
1336 0 : double eNew[3] = { 0., 0., 0.};
1337 0 : eNew[i] = ei;
1338 0 : eNew[j] = ej;
1339 0 : eNew[k] = ek;
1340 :
1341 : // Boost (copy of) partons to their rest frame.
1342 0 : RotBstMatrix Mmove;
1343 0 : Vec4 p0cm = p0;
1344 0 : Vec4 p1cm = p1;
1345 0 : Vec4 p2cm = p2;
1346 0 : Mmove.bstback(pSumJun);
1347 0 : p0cm.rotbst(Mmove);
1348 0 : p1cm.rotbst(Mmove);
1349 0 : p2cm.rotbst(Mmove);
1350 :
1351 : // Construct difference vectors and the boost to junction rest frame.
1352 0 : Vec4 pDir01 = p0cm / p0cm.e() - p1cm / p1cm.e();
1353 0 : Vec4 pDir02 = p0cm / p0cm.e() - p2cm / p2cm.e();
1354 0 : double pDiff01 = pDir01.pAbs2();
1355 0 : double pDiff02 = pDir02.pAbs2();
1356 0 : double pDiff0102 = dot3(pDir01, pDir02);
1357 0 : double eDiff01 = eNew[0] / p0cm.e() - eNew[1] / p1cm.e();
1358 0 : double eDiff02 = eNew[0] / p0cm.e() - eNew[2] / p2cm.e();
1359 0 : double denom = max(1e-6, pDiff01 * pDiff02 - pDiff0102 * pDiff0102);
1360 0 : double coef01 = (eDiff01 * pDiff02 - eDiff02 * pDiff0102) / denom;
1361 0 : double coef02 = (eDiff02 * pDiff01 - eDiff01 * pDiff0102) / denom;
1362 0 : Vec4 vJunction = coef01 * pDir01 + coef02 * pDir02;
1363 0 : vJunction.e( sqrt(1. + vJunction.pAbs2()) );
1364 :
1365 : // Add two boosts, giving final result.
1366 0 : Mmove.bst(vJunction);
1367 : return Mmove;
1368 :
1369 0 : }
1370 :
1371 : //--------------------------------------------------------------------------
1372 :
1373 : // When string fragmentation has failed several times,
1374 : // try to join some more nearby partons.
1375 :
1376 : int StringFragmentation::extraJoin(double facExtra, Event& event) {
1377 :
1378 : // Keep on looping while pairs found below joining threshold.
1379 : int nJoin = 0;
1380 0 : int iPsize = iParton.size();
1381 0 : while (iPsize > 2) {
1382 :
1383 : // Look for the pair of neighbour partons (along string) with
1384 : // the smallest invariant mass (subtracting quark masses).
1385 : int iJoinMin = -1;
1386 0 : double mJoinMin = 2. * facExtra * mJoin;
1387 0 : for (int i = 0; i < iPsize - 1; ++i) {
1388 0 : Particle& parton1 = event[ iParton[i] ];
1389 0 : Particle& parton2 = event[ iParton[i + 1] ];
1390 0 : Vec4 pSumNow;
1391 0 : pSumNow += (parton2.isGluon()) ? 0.5 * parton1.p() : parton1.p();
1392 0 : pSumNow += (parton2.isGluon()) ? 0.5 * parton2.p() : parton2.p();
1393 0 : double mJoinNow = pSumNow.mCalc();
1394 0 : if (!parton1.isGluon()) mJoinNow -= parton1.m0();
1395 0 : if (!parton2.isGluon()) mJoinNow -= parton2.m0();
1396 0 : if (mJoinNow < mJoinMin) { iJoinMin = i; mJoinMin = mJoinNow; }
1397 0 : }
1398 :
1399 : // Decide whether to join, if not finished.
1400 0 : if (iJoinMin < 0 || mJoinMin > facExtra * mJoin) return nJoin;
1401 0 : ++nJoin;
1402 :
1403 : // Create new joined parton.
1404 0 : int iJoin1 = iParton[iJoinMin];
1405 0 : int iJoin2 = iParton[iJoinMin + 1];
1406 0 : int idNew = (event[iJoin1].isGluon()) ? event[iJoin2].id()
1407 0 : : event[iJoin1].id();
1408 0 : int colNew = event[iJoin1].col();
1409 0 : int acolNew = event[iJoin2].acol();
1410 0 : if (colNew == acolNew) {
1411 0 : colNew = event[iJoin2].col();
1412 0 : acolNew = event[iJoin1].acol();
1413 0 : }
1414 0 : Vec4 pNew = event[iJoin1].p() + event[iJoin2].p();
1415 :
1416 : // Append joined parton to event record and reduce parton list.
1417 0 : int iNew = event.append( idNew, 73, min(iJoin1, iJoin2),
1418 0 : max(iJoin1, iJoin2), 0, 0, colNew, acolNew, pNew, pNew.mCalc() );
1419 0 : iParton[iJoinMin] = iNew;
1420 0 : for (int i = iJoinMin + 1; i < iPsize - 1; ++i)
1421 0 : iParton[i] = iParton[i + 1];
1422 0 : iParton.pop_back();
1423 0 : --iPsize;
1424 :
1425 : // Done.
1426 0 : }
1427 0 : return nJoin;
1428 0 : }
1429 :
1430 : //==========================================================================
1431 :
1432 : } // end namespace Pythia8
|
