]>
Commit | Line | Data |
---|---|---|
b584e2f5 | 1 | // StringFragmentation.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2010 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 "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 | fromPos = fromPosIn; | |
36 | iEnd = iEndIn; | |
37 | iMax = iMaxIn; | |
38 | flavOld = FlavContainer(idOldIn); | |
39 | pxOld = pxIn; | |
40 | pyOld = pyIn; | |
41 | GammaOld = GammaIn; | |
42 | iPosOld = (fromPos) ? 0 : iMax; | |
43 | iNegOld = (fromPos) ? iMax : 0; | |
44 | xPosOld = xPosIn; | |
45 | xNegOld = xNegIn; | |
46 | ||
47 | } | |
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 | do { | |
57 | flavNew = flavSelPtr->pick( flavOld); | |
58 | idHad = flavSelPtr->combine( flavOld, flavNew); | |
59 | } while (idHad == 0); | |
60 | ||
61 | // Pick its transverse momentum. | |
62 | pair<double, double> pxy = pTSelPtr->pxy(); | |
63 | pxNew = pxy.first; | |
64 | pyNew = pxy.second; | |
65 | pxHad = pxOld + pxNew; | |
66 | pyHad = pyOld + pyNew; | |
67 | ||
68 | // Pick its mass and thereby define its transverse mass. | |
69 | mHad = particleDataPtr->mass(idHad); | |
70 | mT2Had = pow2(mHad) + pow2(pxHad) + pow2(pyHad); | |
71 | ||
72 | } | |
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 | zHad = zSelPtr->zFrag( flavOld.id, flavNew.id, mT2Had); | |
83 | 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 | int& iDirOld = (fromPos) ? iPosOld : iNegOld; | |
88 | int& iInvOld = (fromPos) ? iNegOld : iPosOld; | |
89 | int& iDirNew = (fromPos) ? iPosNew : iNegNew; | |
90 | int& iInvNew = (fromPos) ? iNegNew : iPosNew; | |
91 | double& xDirOld = (fromPos) ? xPosOld : xNegOld; | |
92 | double& xInvOld = (fromPos) ? xNegOld : xPosOld; | |
93 | double& xDirNew = (fromPos) ? xPosNew : xNegNew; | |
94 | double& xInvNew = (fromPos) ? xNegNew : xPosNew; | |
95 | double& xDirHad = (fromPos) ? xPosHad : xNegHad; | |
96 | double& xInvHad = (fromPos) ? xNegHad : xPosHad; | |
97 | ||
98 | // Start search for new breakup in the old region. | |
99 | iDirNew = iDirOld; | |
100 | iInvNew = iInvOld; | |
101 | Vec4 pTNew; | |
102 | ||
103 | // Each step corresponds to trying a new string region. | |
104 | for (int iStep = 0; ; ++iStep) { | |
105 | ||
106 | // Referance to current string region. | |
107 | StringRegion& region = system.region( iPosNew, iNegNew); | |
108 | ||
109 | // Now begin special section for rapid processing of low region. | |
110 | if (iStep == 0 && iPosOld + iNegOld == iMax) { | |
111 | ||
112 | // A first step within a low region is easy. | |
113 | if (mT2Had < zHad * xDirOld * (1. - xInvOld) * region.w2) { | |
114 | ||
115 | // Translate into x coordinates. | |
116 | xDirHad = zHad * xDirOld; | |
117 | xInvHad = mT2Had / (xDirHad * region.w2); | |
118 | xDirNew = xDirOld - xDirHad; | |
119 | xInvNew = xInvOld + xInvHad; | |
120 | ||
121 | // Find and return four-momentum of the produced particle. | |
122 | 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 | --iInvNew; | |
128 | if (iInvNew < 0) return Vec4(0., 0., 0., -1.); | |
129 | ||
130 | // Momentum taken by stepping out of region. Continue to next region. | |
131 | xInvHad = 1. - xInvOld; | |
132 | xDirHad = 0.; | |
133 | pSoFar = region.pHad( xPosHad, xNegHad, pxOld, pyOld); | |
134 | continue; | |
135 | } | |
136 | ||
137 | // Else, for first step, take into account starting pT. | |
138 | } else if (iStep == 0) { | |
139 | pSoFar = region.pHad( 0., 0., pxOld, pyOld); | |
140 | pTNew = region.pHad( 0., 0., pxNew, pyNew); | |
141 | } | |
142 | ||
143 | // Now begin normal treatment of nontrivial regions. | |
144 | // Set up four-vectors in a region not visited before. | |
145 | if (!region.isSetUp) region.setUp( | |
146 | system.regionLowPos(iPosNew).pPos, | |
147 | 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 | if (region.isEmpty) { | |
152 | xDirHad = (iDirNew == iDirOld) ? xDirOld : 1.; | |
153 | xInvHad = 0.; | |
154 | pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.); | |
155 | ++iDirNew; | |
156 | if (iDirNew + iInvNew > iMax) return Vec4(0., 0., 0., -1.); | |
157 | 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 | double pxNewTemp = -pTNew * region.eX; | |
163 | double pyNewTemp = -pTNew * region.eY; | |
164 | if (abs( pxNewTemp * pxNewTemp + pyNewTemp * pyNewTemp | |
165 | - pxNew * pxNew - pyNew * pyNew) < PT2SAME) { | |
166 | pxNew = pxNewTemp; | |
167 | pyNew = pyNewTemp; | |
168 | } | |
169 | ||
170 | // Four-momentum taken so far, including new pT. | |
171 | 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 | double cM1 = pTemp.m2Calc(); | |
176 | double cM2 = 2. * (pTemp * region.pPos); | |
177 | double cM3 = 2. * (pTemp * region.pNeg); | |
178 | double cM4 = region.w2; | |
179 | 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 | for (int iInv = iInvNew; iInv <= iMax - iDirNew; ++iInv) { | |
188 | double xInv = 1.; | |
189 | if (iInv == iInvNew) xInv = (iInvNew == iInvOld) ? xInvOld : 0.; | |
190 | for (int iDir = iDirNew; iDir <= iMax - iInv; ++iDir) { | |
191 | double xDir = (iDir == iDirOld) ? xDirOld : 1.; | |
192 | int iPos = (fromPos) ? iDir : iInv; | |
193 | int iNeg = (fromPos) ? iInv : iDir; | |
194 | StringRegion& regionGam = system.region( iPos, iNeg); | |
195 | if (!regionGam.isSetUp) regionGam.setUp( | |
196 | system.regionLowPos(iPos).pPos, | |
197 | system.regionLowNeg(iNeg).pNeg, true); | |
198 | double w2 = regionGam.w2; | |
199 | cGam1 += xDir * xInv * w2; | |
200 | if (iDir == iDirNew) cGam2 -= xInv * w2; | |
201 | if (iInv == iInvNew) cGam3 += xDir * w2; | |
202 | if (iDir == iDirNew && iInv == iInvNew) cGam4 -= w2; | |
203 | } | |
204 | } | |
205 | ||
206 | // Solve (m2, Gamma) equation system => r2 * x-^2 + r1 * x- + r0 = 0. | |
207 | double cM0 = pow2(mHad) - cM1; | |
208 | double cGam0 = GammaNew - cGam1; | |
209 | double r2 = cM3 * cGam4 - cM4 * cGam3; | |
210 | double r1 = cM4 * cGam0 - cM0 * cGam4 + cM3 * cGam2 - cM2 * cGam3; | |
211 | double r0 = cM2 * cGam0 - cM0 * cGam2; | |
212 | double root = sqrtpos( r1*r1 - 4. * r2 * r0 ); | |
213 | if (abs(r2) < TINY || root < TINY) return Vec4(0., 0., 0., -1.); | |
214 | xInvHad = 0.5 * (root / abs(r2) - r1 / r2); | |
215 | xDirHad = (cM0 - cM3 * xInvHad) / (cM2 + cM4 * xInvHad); | |
216 | ||
217 | // Define position of new trial vertex. | |
218 | xDirNew = (iDirNew == iDirOld) ? xDirOld - xDirHad : 1. - xDirHad; | |
219 | xInvNew = (iInvNew == iInvOld) ? xInvOld + xInvHad : xInvHad; | |
220 | ||
221 | // Step up to new region if new x- > 1. | |
222 | if (xInvNew > 1.) { | |
223 | xInvHad = (iInvNew == iInvOld) ? 1. - xInvOld : 1.; | |
224 | xDirHad = 0.; | |
225 | pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.); | |
226 | --iInvNew; | |
227 | if (iInvNew < 0) return Vec4(0., 0., 0., -1.); | |
228 | continue; | |
229 | ||
230 | // Step down to new region if new x+ < 0. | |
231 | } else if (xDirNew < 0.) { | |
232 | xDirHad = (iDirNew == iDirOld) ? xDirOld : 1.; | |
233 | xInvHad = 0.; | |
234 | pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.); | |
235 | ++iDirNew; | |
236 | if (iDirNew + iInvNew > iMax) return Vec4(0., 0., 0., -1.); | |
237 | continue; | |
238 | } | |
239 | ||
240 | // Else we have found the correct region, and can return four-momentum. | |
241 | return pSoFar + region.pHad( xPosHad, xNegHad, pxNew, pyNew); | |
242 | ||
243 | // End of "infinite" loop of stepping to new region. | |
244 | } | |
245 | ||
246 | } | |
247 | ||
248 | //-------------------------------------------------------------------------- | |
249 | ||
250 | // Update string end information after a hadron has been removed. | |
251 | ||
252 | void StringEnd::update() { | |
253 | ||
254 | flavOld.anti(flavNew); | |
255 | iPosOld = iPosNew; | |
256 | iNegOld = iNegNew; | |
257 | pxOld = -pxNew; | |
258 | pyOld = -pyNew; | |
259 | GammaOld = GammaNew; | |
260 | xPosOld = xPosNew; | |
261 | xNegOld = xNegNew; | |
262 | ||
263 | } | |
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 = 25; | |
277 | ||
278 | // The last few times gradually increase the stop mass to make it easier. | |
279 | const int StringFragmentation::NSTOPMASS = 10; | |
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-6; | |
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 | ||
302 | // Consider junction-leg parton as massless if m2 tiny. | |
303 | const double StringFragmentation::M2MAXJRF = 1e-4; | |
304 | ||
305 | // Iterate junction rest frame equation until convergence or too many tries. | |
306 | const double StringFragmentation::CONVJRFEQ = 1e-12; | |
307 | const int StringFragmentation::NTRYJRFEQ = 40; | |
308 | ||
309 | //-------------------------------------------------------------------------- | |
310 | ||
311 | // Initialize and save pointers. | |
312 | ||
313 | void StringFragmentation::init(Info* infoPtrIn, Settings& settings, | |
314 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, StringFlav* flavSelPtrIn, | |
315 | StringPT* pTSelPtrIn, StringZ* zSelPtrIn) { | |
316 | ||
317 | // Save pointers. | |
318 | infoPtr = infoPtrIn; | |
319 | particleDataPtr = particleDataPtrIn; | |
320 | rndmPtr = rndmPtrIn; | |
321 | flavSelPtr = flavSelPtrIn; | |
322 | pTSelPtr = pTSelPtrIn; | |
323 | zSelPtr = zSelPtrIn; | |
324 | ||
325 | // Initialize the StringFragmentation class. | |
326 | stopMass = settings.parm("StringFragmentation:stopMass"); | |
327 | stopNewFlav = settings.parm("StringFragmentation:stopNewFlav"); | |
328 | stopSmear = settings.parm("StringFragmentation:stopSmear"); | |
329 | eNormJunction = settings.parm("StringFragmentation:eNormJunction"); | |
330 | eBothLeftJunction | |
331 | = settings.parm("StringFragmentation:eBothLeftJunction"); | |
332 | eMaxLeftJunction | |
333 | = settings.parm("StringFragmentation:eMaxLeftJunction"); | |
334 | eMinLeftJunction | |
335 | = settings.parm("StringFragmentation:eMinLeftJunction"); | |
336 | ||
337 | // Initialize the b parameter of the z spectrum, used when joining jets. | |
338 | bLund = settings.parm("StringZ:bLund"); | |
339 | ||
340 | // Initialize the hadrons instance of an event record. | |
341 | hadrons.init( "(string fragmentation)", particleDataPtr); | |
342 | ||
343 | // Send on pointers to the two StringEnd instances. | |
344 | posEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr); | |
345 | negEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr); | |
346 | ||
347 | } | |
348 | ||
349 | //-------------------------------------------------------------------------- | |
350 | ||
351 | // Perform the fragmentation. | |
352 | ||
353 | bool StringFragmentation::fragment( int iSub, ColConfig& colConfig, | |
354 | Event& event) { | |
355 | ||
356 | // Find partons and their total four-momentum. | |
357 | iParton = colConfig[iSub].iParton; | |
358 | iPos = iParton[0]; | |
359 | if (iPos < 0) iPos = iParton[1]; | |
360 | int idPos = event[iPos].id(); | |
361 | iNeg = iParton.back(); | |
362 | int idNeg = event[iNeg].id(); | |
363 | pSum = colConfig[iSub].pSum; | |
364 | ||
365 | // Reset the local event record. | |
366 | hadrons.clear(); | |
367 | ||
368 | // For closed gluon string: pick first breakup region. | |
369 | isClosed = colConfig[iSub].isClosed; | |
370 | if (isClosed) iParton = findFirstRegion(iParton, event); | |
371 | ||
372 | // For junction topology: fragment off two of the string legs. | |
373 | // Then iParton overwritten to remaining leg + leftover diquark. | |
374 | pJunctionHadrons = 0.; | |
375 | hasJunction = colConfig[iSub].hasJunction; | |
376 | if (hasJunction && !fragmentToJunction(event)) return false; | |
377 | int junctionHadrons = hadrons.size(); | |
378 | if (hasJunction) { | |
379 | idPos = event[ iParton[0] ].id(); | |
380 | idNeg = event.back().id(); | |
381 | pSum -= pJunctionHadrons; | |
382 | } | |
383 | ||
384 | // Set up kinematics of string evolution ( = motion). | |
385 | system.setUp(iParton, event); | |
386 | stopMassNow = stopMass; | |
387 | ||
388 | // Fallback loop, when joining in the middle fails. Bailout if stuck. | |
389 | for ( int iTry = 0; ; ++iTry) { | |
390 | if (iTry > NTRYJOIN) { | |
391 | infoPtr->errorMsg("Error in StringFragmentation::fragment: " | |
392 | "stuck in joining"); | |
393 | if (hasJunction) event.popBack(1); | |
394 | return false; | |
395 | } | |
396 | ||
397 | // After several failed tries gradually allow larger stop mass. | |
398 | if (iTry > NTRYJOIN - NSTOPMASS) stopMassNow *= FACSTOPMASS; | |
399 | ||
400 | // Set up flavours of two string ends, and reset other info. | |
401 | setStartEnds(idPos, idNeg, system); | |
402 | pRem = pSum; | |
403 | ||
404 | // Begin fragmentation loop, interleaved from the two ends. | |
405 | bool fromPos; | |
406 | for ( ; ; ) { | |
407 | ||
408 | // Take a step either from the positive or the negative end. | |
409 | fromPos = (rndmPtr->flat() < 0.5); | |
410 | StringEnd& nowEnd = (fromPos) ? posEnd : negEnd; | |
411 | ||
412 | // Construct trial hadron and check that energy remains. | |
413 | nowEnd.newHadron(); | |
414 | if ( energyUsedUp(fromPos) ) break; | |
415 | ||
416 | // Construct kinematics of the new hadron and store it. | |
417 | Vec4 pHad = nowEnd.kinematicsHadron(system); | |
418 | int statusHad = (fromPos) ? 83 : 84; | |
419 | hadrons.append( nowEnd.idHad, statusHad, iPos, iNeg, | |
420 | 0, 0, 0, 0, pHad, nowEnd.mHad); | |
421 | if (pHad.e() < 0.) break; | |
422 | ||
423 | // Update string end and remaining momentum. | |
424 | nowEnd.update(); | |
425 | pRem -= pHad; | |
426 | ||
427 | // End of fragmentation loop. | |
428 | } | |
429 | ||
430 | // When done, join in the middle. If this works, then really done. | |
431 | if ( finalTwo(fromPos) ) break; | |
432 | ||
433 | // Else remove produced particles (except from first two junction legs) | |
434 | // and start all over. | |
435 | int newHadrons = hadrons.size() - junctionHadrons; | |
436 | hadrons.popBack(newHadrons); | |
437 | } | |
438 | ||
439 | // Junctions: remove fictitious end, restore original parton list | |
440 | if (hasJunction) { | |
441 | event.popBack(1); | |
442 | iParton = colConfig[iSub].iParton; | |
443 | } | |
444 | ||
445 | // Store the hadrons in the normal event record, ordered from one end. | |
446 | store(event); | |
447 | ||
448 | // Done. | |
449 | return true; | |
450 | ||
451 | } | |
452 | ||
453 | //-------------------------------------------------------------------------- | |
454 | ||
455 | // Find region where to put first string break for closed gluon loop. | |
456 | ||
457 | vector<int> StringFragmentation::findFirstRegion(vector<int>& iPartonIn, | |
458 | Event& event) { | |
459 | ||
460 | // Evaluate mass-squared for all adjacent gluon pairs. | |
461 | vector<double> m2Pair; | |
462 | double m2Sum = 0.; | |
463 | int size = iPartonIn.size(); | |
464 | for (int i = 0; i < size; ++i) { | |
465 | double m2Now = 0.5 * event[ iPartonIn[i] ].p() | |
466 | * event[ iPartonIn[(i + 1)%size] ].p(); | |
467 | m2Pair.push_back(m2Now); | |
468 | m2Sum += m2Now; | |
469 | } | |
470 | ||
471 | // Pick breakup region with probability proportional to mass-squared. | |
472 | double m2Reg = m2Sum * rndmPtr->flat(); | |
473 | int iReg = -1; | |
474 | do m2Reg -= m2Pair[++iReg]; | |
475 | while (m2Reg > 0. && iReg < size - 1); | |
476 | ||
477 | // Create reordered parton list, with breakup string region duplicated. | |
478 | vector<int> iPartonOut; | |
479 | for (int i = 0; i < size + 2; ++i) | |
480 | iPartonOut.push_back( iPartonIn[(i + iReg)%size] ); | |
481 | ||
482 | // Done. | |
483 | return iPartonOut; | |
484 | ||
485 | } | |
486 | ||
487 | //-------------------------------------------------------------------------- | |
488 | ||
489 | // Set flavours and momentum position for initial string endpoints. | |
490 | ||
491 | void StringFragmentation::setStartEnds( int idPos, int idNeg, | |
492 | StringSystem systemNow) { | |
493 | ||
494 | // Variables characterizing string endpoints: defaults for open string. | |
495 | double px = 0.; | |
496 | double py = 0.; | |
497 | double Gamma = 0.; | |
498 | double xPosFromPos = 1.; | |
499 | double xNegFromPos = 0.; | |
500 | double xPosFromNeg = 0.; | |
501 | double xNegFromNeg = 1.; | |
502 | ||
503 | // For closed gluon loop need to pick an initial flavour. | |
504 | if (isClosed) { | |
505 | do { | |
506 | int idTry = flavSelPtr->pickLightQ(); | |
507 | FlavContainer flavTry(idTry, 1); | |
508 | flavTry = flavSelPtr->pick( flavTry); | |
509 | flavTry = flavSelPtr->pick( flavTry); | |
510 | idPos = flavTry.id; | |
511 | idNeg = -idPos; | |
512 | } while (idPos == 0); | |
513 | ||
514 | // Also need pT and breakup vertex position in region. | |
515 | pair<double, double> pxy = pTSelPtr->pxy(); | |
516 | px = pxy.first; | |
517 | py = pxy.second; | |
518 | double m2Region = systemNow.regionLowPos(0).w2; | |
519 | double m2Temp = min( CLOSEDM2MAX, CLOSEDM2FRAC * m2Region); | |
520 | do { | |
521 | double zTemp = zSelPtr->zFrag( idPos, idNeg, m2Temp); | |
522 | xPosFromPos = 1. - zTemp; | |
523 | xNegFromPos = m2Temp / (zTemp * m2Region); | |
524 | } while (xNegFromPos > 1.); | |
525 | Gamma = xPosFromPos * xNegFromPos * m2Region; | |
526 | xPosFromNeg = xPosFromPos; | |
527 | xNegFromNeg = xNegFromPos; | |
528 | } | |
529 | ||
530 | // Initialize two string endpoints. | |
531 | posEnd.setUp( true, iPos, idPos, systemNow.iMax, px, py, | |
532 | Gamma, xPosFromPos, xNegFromPos); | |
533 | negEnd.setUp( false, iNeg, idNeg, systemNow.iMax, -px, -py, | |
534 | Gamma, xPosFromNeg, xNegFromNeg); | |
535 | ||
536 | // For closed gluon loop can allow popcorn on one side but not both. | |
537 | if (isClosed) { | |
538 | flavSelPtr->assignPopQ(posEnd.flavOld); | |
539 | flavSelPtr->assignPopQ(negEnd.flavOld); | |
540 | if (rndmPtr->flat() < 0.5) posEnd.flavOld.nPop = 0; | |
541 | else negEnd.flavOld.nPop = 0; | |
542 | posEnd.flavOld.rank = 1; | |
543 | negEnd.flavOld.rank = 1; | |
544 | } | |
545 | ||
546 | // Done. | |
547 | ||
548 | } | |
549 | ||
550 | //-------------------------------------------------------------------------- | |
551 | ||
552 | // Check remaining energy-momentum whether it is OK to continue. | |
553 | ||
554 | bool StringFragmentation::energyUsedUp(bool fromPos) { | |
555 | ||
556 | // If remaining negative energy then abort right away. | |
557 | if (pRem.e() < 0.) return true; | |
558 | ||
559 | // Calculate W2_minimum and done if remaining W2 is below it. | |
560 | double wMin = stopMassNow | |
561 | + particleDataPtr->constituentMass(posEnd.flavOld.id) | |
562 | + particleDataPtr->constituentMass(negEnd.flavOld.id); | |
563 | if (fromPos) wMin += stopNewFlav | |
564 | * particleDataPtr->constituentMass(posEnd.flavNew.id); | |
565 | else wMin += stopNewFlav | |
566 | * particleDataPtr->constituentMass(negEnd.flavNew.id); | |
567 | wMin *= 1. + (2. * rndmPtr->flat() - 1.) * stopSmear; | |
568 | w2Rem = pRem.m2Calc(); | |
569 | if (w2Rem < pow2(wMin)) return true; | |
570 | ||
571 | // Else still enough energy left to continue iteration. | |
572 | return false; | |
573 | ||
574 | } | |
575 | ||
576 | //-------------------------------------------------------------------------- | |
577 | ||
578 | // Produce the final two partons to complete the system. | |
579 | ||
580 | bool StringFragmentation::finalTwo(bool fromPos) { | |
581 | ||
582 | // Check whether we went too far in p+-. | |
583 | if (pRem.e() < 0. || w2Rem < 0. || (hadrons.size() > 0 | |
584 | && hadrons.back().e() < 0.) ) return false; | |
585 | if ( posEnd.iPosOld > negEnd.iPosOld || negEnd.iNegOld > posEnd.iNegOld) | |
586 | return false; | |
587 | if ( posEnd.iPosOld == negEnd.iPosOld && posEnd.xPosOld < negEnd.xPosOld) | |
588 | return false; | |
589 | if ( posEnd.iNegOld == negEnd.iNegOld && posEnd.xNegOld > negEnd.xNegOld) | |
590 | return false; | |
591 | ||
592 | // Construct the final hadron from the leftover flavours. | |
593 | // Impossible to join two diquarks. Also break if stuck for other reason. | |
594 | FlavContainer flav1 = (fromPos) ? posEnd.flavNew.anti() : posEnd.flavOld; | |
595 | FlavContainer flav2 = (fromPos) ? negEnd.flavOld : negEnd.flavNew.anti(); | |
596 | if (abs(flav1.id) > 8 && abs(flav2.id) > 8) return false; | |
597 | int idHad; | |
598 | for (int iTry = 0; iTry < NTRYFLAV; ++iTry) { | |
599 | idHad = flavSelPtr->combine( flav1, flav2); | |
600 | if (idHad != 0) break; | |
601 | } | |
602 | if (idHad == 0) return false; | |
603 | ||
604 | // Store the final particle and its new pT, and construct its mass. | |
605 | if (fromPos) { | |
606 | negEnd.idHad = idHad; | |
607 | negEnd.pxNew = -posEnd.pxNew; | |
608 | negEnd.pyNew = -posEnd.pyNew; | |
609 | negEnd.mHad = particleDataPtr->mass(idHad); | |
610 | } else { | |
611 | posEnd.idHad = idHad; | |
612 | posEnd.pxNew = -negEnd.pxNew; | |
613 | posEnd.pyNew = -negEnd.pyNew; | |
614 | posEnd.mHad = particleDataPtr->mass(idHad); | |
615 | } | |
616 | ||
617 | // String region in which to do the joining. | |
618 | StringRegion region = finalRegion(); | |
619 | if (region.isEmpty) return false; | |
620 | ||
621 | // Project remaining momentum along longitudinal and transverse directions. | |
622 | region.project( pRem); | |
623 | double pxRem = region.px() - posEnd.pxOld - negEnd.pxOld; | |
624 | double pyRem = region.py() - posEnd.pyOld - negEnd.pyOld; | |
625 | double xPosRem = region.xPos(); | |
626 | double xNegRem = region.xNeg(); | |
627 | ||
628 | // Share extra pT kick evenly between final two hadrons. | |
629 | posEnd.pxOld += 0.5 * pxRem; | |
630 | posEnd.pyOld += 0.5 * pyRem; | |
631 | negEnd.pxOld += 0.5 * pxRem; | |
632 | negEnd.pyOld += 0.5 * pyRem; | |
633 | ||
634 | // Construct new pT and mT of the final two particles. | |
635 | posEnd.pxHad = posEnd.pxOld + posEnd.pxNew; | |
636 | posEnd.pyHad = posEnd.pyOld + posEnd.pyNew; | |
637 | posEnd.mT2Had = pow2(posEnd.mHad) + pow2(posEnd.pxHad) | |
638 | + pow2(posEnd.pyHad); | |
639 | negEnd.pxHad = negEnd.pxOld + negEnd.pxNew; | |
640 | negEnd.pyHad = negEnd.pyOld + negEnd.pyNew; | |
641 | negEnd.mT2Had = pow2(negEnd.mHad) + pow2(negEnd.pxHad) | |
642 | + pow2(negEnd.pyHad); | |
643 | ||
644 | // Construct remaining system transverse mass. | |
645 | double wT2Rem = w2Rem + pow2( posEnd.pxHad + negEnd.pxHad) | |
646 | + pow2( posEnd.pyHad + negEnd.pyHad); | |
647 | ||
648 | // Check that kinematics possible. | |
649 | if ( sqrt(wT2Rem) < sqrt(posEnd.mT2Had) + sqrt(negEnd.mT2Had) ) | |
650 | return false; | |
651 | double lambda2 = pow2( wT2Rem - posEnd.mT2Had - negEnd.mT2Had) | |
652 | - 4. * posEnd.mT2Had * negEnd.mT2Had; | |
653 | if (lambda2 <= 0.) return false; | |
654 | ||
655 | // Construct kinematics, as viewed in the transverse rest frame. | |
656 | double lambda = sqrt(lambda2); | |
657 | double probReverse = 1. / (1. + exp( min( EXPMAX, bLund * lambda) ) ); | |
658 | double xpzPos = 0.5 * lambda/ wT2Rem; | |
659 | if (probReverse > rndmPtr->flat()) xpzPos = -xpzPos; | |
660 | double xmDiff = (posEnd.mT2Had - negEnd.mT2Had) / wT2Rem; | |
661 | double xePos = 0.5 * (1. + xmDiff); | |
662 | double xeNeg = 0.5 * (1. - xmDiff ); | |
663 | ||
664 | // Translate this into kinematics in the string frame. | |
665 | Vec4 pHadPos = region.pHad( (xePos + xpzPos) * xPosRem, | |
666 | (xePos - xpzPos) * xNegRem, posEnd.pxHad, posEnd.pyHad); | |
667 | Vec4 pHadNeg = region.pHad( (xeNeg - xpzPos) * xPosRem, | |
668 | (xeNeg + xpzPos) * xNegRem, negEnd.pxHad, negEnd.pyHad); | |
669 | ||
670 | // Add produced particles to the event record. | |
671 | hadrons.append( posEnd.idHad, 83, posEnd.iEnd, negEnd.iEnd, | |
672 | 0, 0, 0, 0, pHadPos, posEnd.mHad); | |
673 | hadrons.append( negEnd.idHad, 84, posEnd.iEnd, negEnd.iEnd, | |
674 | 0, 0, 0, 0, pHadNeg, negEnd.mHad); | |
675 | ||
676 | // It worked. | |
677 | return true; | |
678 | ||
679 | } | |
680 | ||
681 | //-------------------------------------------------------------------------- | |
682 | ||
683 | // Construct a special joining region for the final two hadrons. | |
684 | ||
685 | StringRegion StringFragmentation::finalRegion() { | |
686 | ||
687 | // Simple case when both string ends are in the same region. | |
688 | if (posEnd.iPosOld == negEnd.iPosOld && posEnd.iNegOld == negEnd.iNegOld) | |
689 | return system.region( posEnd.iPosOld, posEnd.iNegOld); | |
690 | ||
691 | // Start out with empty region. (Empty used for error returns.) | |
692 | StringRegion region; | |
693 | ||
694 | // Add up all remaining p+. | |
695 | Vec4 pPosJoin; | |
696 | if ( posEnd.iPosOld == negEnd.iPosOld) { | |
697 | double xPosJoin = posEnd.xPosOld - negEnd.xPosOld; | |
698 | if (xPosJoin < 0.) return region; | |
699 | pPosJoin = system.regionLowPos(posEnd.iPosOld).pHad( xPosJoin, 0., 0., 0.); | |
700 | } else { | |
701 | for (int iPosNow = posEnd.iPosOld; iPosNow <= negEnd.iPosOld; ++iPosNow) { | |
702 | if (iPosNow == posEnd.iPosOld) pPosJoin | |
703 | += system.regionLowPos(iPosNow).pHad( posEnd.xPosOld, 0., 0., 0.); | |
704 | else if (iPosNow == negEnd.iPosOld) pPosJoin | |
705 | += system.regionLowPos(iPosNow).pHad( 1. - negEnd.xPosOld, 0., 0., 0.); | |
706 | else pPosJoin += system.regionLowPos(iPosNow).pHad( 1., 0., 0., 0.); | |
707 | } | |
708 | } | |
709 | ||
710 | // Add up all remaining p-. | |
711 | Vec4 pNegJoin; | |
712 | if ( negEnd.iNegOld == posEnd.iNegOld) { | |
713 | double xNegJoin = negEnd.xNegOld - posEnd.xNegOld; | |
714 | if (xNegJoin < 0.) return region; | |
715 | pNegJoin = system.regionLowNeg(negEnd.iNegOld).pHad( 0., xNegJoin, 0., 0.); | |
716 | } else { | |
717 | for (int iNegNow = negEnd.iNegOld; iNegNow <= posEnd.iNegOld; ++iNegNow) { | |
718 | if (iNegNow == negEnd.iNegOld) pNegJoin | |
719 | += system.regionLowNeg(iNegNow).pHad( 0., negEnd.xNegOld, 0., 0.); | |
720 | else if (iNegNow == posEnd.iNegOld) pNegJoin | |
721 | += system.regionLowNeg(iNegNow).pHad( 0., 1. - posEnd.xNegOld, 0., 0.); | |
722 | else pNegJoin += system.regionLowNeg(iNegNow).pHad( 0., 1., 0., 0.); | |
723 | } | |
724 | } | |
725 | ||
726 | // For a closed gluon loop pPosJoin == pNegJoin and the above does not work. | |
727 | // So reshuffle; "perfect" for g g systems, OK in general. | |
728 | Vec4 pTest = pPosJoin - pNegJoin; | |
729 | if ( abs(pTest.px()) + abs(pTest.py()) + abs(pTest.pz()) + abs(pTest.e()) | |
730 | < MATCHPOSNEG * (pPosJoin.e() + pNegJoin.e()) ) { | |
731 | Vec4 delta | |
732 | = system.regionLowPos(posEnd.iPosOld + 1).pHad( 1., 0., 0., 0.) | |
733 | - system.regionLowNeg(negEnd.iNegOld + 1).pHad( 0., 1., 0., 0.); | |
734 | pPosJoin -= delta; | |
735 | pNegJoin += delta; | |
736 | } | |
737 | ||
738 | // Construct a new region from remaining p+ and p-. | |
739 | region.setUp( pPosJoin, pNegJoin); | |
740 | if (region.isEmpty) return region; | |
741 | ||
742 | // Project the existing pTold vectors onto the new directions. | |
743 | Vec4 pTposOld = system.region( posEnd.iPosOld, posEnd.iNegOld).pHad( | |
744 | 0., 0., posEnd.pxOld, posEnd.pyOld); | |
745 | region.project( pTposOld); | |
746 | posEnd.pxOld = region.px(); | |
747 | posEnd.pyOld = region.py(); | |
748 | Vec4 pTnegOld = system.region( negEnd.iPosOld, negEnd.iNegOld).pHad( | |
749 | 0., 0., negEnd.pxOld, negEnd.pyOld); | |
750 | region.project( pTnegOld); | |
751 | negEnd.pxOld = region.px(); | |
752 | negEnd.pyOld = region.py(); | |
753 | ||
754 | // Done. | |
755 | return region; | |
756 | ||
757 | } | |
758 | ||
759 | //-------------------------------------------------------------------------- | |
760 | ||
761 | // Store the hadrons in the normal event record, ordered from one end. | |
762 | ||
763 | void StringFragmentation::store(Event& event) { | |
764 | ||
765 | // Starting position. | |
766 | int iFirst = event.size(); | |
767 | ||
768 | // Copy straight over from first two junction legs. | |
769 | if (hasJunction) { | |
770 | for (int i = 0; i < hadrons.size(); ++i) | |
771 | if (hadrons[i].status() == 85 || hadrons[i].status() == 86) | |
772 | event.append( hadrons[i] ); | |
773 | } | |
774 | ||
775 | // Loop downwards, copying all from the positive end. | |
776 | for (int i = 0; i < hadrons.size(); ++i) | |
777 | if (hadrons[i].status() == 83) event.append( hadrons[i] ); | |
778 | ||
779 | // Loop upwards, copying all from the negative end. | |
780 | for (int i = hadrons.size() - 1; i >= 0 ; --i) | |
781 | if (hadrons[i].status() == 84) event.append( hadrons[i] ); | |
782 | int iLast = event.size() - 1; | |
783 | ||
784 | // Set decay vertex when this is displaced. | |
785 | if (event[posEnd.iEnd].hasVertex()) { | |
786 | Vec4 vDec = event[posEnd.iEnd].vDec(); | |
787 | for (int i = iFirst; i <= iLast; ++i) event[i].vProd( vDec ); | |
788 | } | |
789 | ||
790 | // Set lifetime of hadrons. | |
791 | for (int i = iFirst; i <= iLast; ++i) | |
792 | event[i].tau( event[i].tau0() * rndmPtr->exp() ); | |
793 | ||
794 | // Mark original partons as hadronized and set their daughter range. | |
795 | for (int i = 0; i < int(iParton.size()); ++i) | |
796 | if (iParton[i] >= 0) { | |
797 | event[ iParton[i] ].statusNeg(); | |
798 | event[ iParton[i] ].daughters(iFirst, iLast); | |
799 | } | |
800 | ||
801 | } | |
802 | ||
803 | //-------------------------------------------------------------------------- | |
804 | ||
805 | // Fragment off two of the string legs in to a junction. | |
806 | ||
807 | bool StringFragmentation::fragmentToJunction(Event& event) { | |
808 | ||
809 | // Identify range of partons on the three legs. | |
810 | // (Each leg begins with an iParton[i] = -(10 + 10*junctionNuber + leg), | |
811 | // and partons then appear ordered from the junction outwards.) | |
812 | int legBeg[3], legEnd[3]; | |
813 | int leg = -1; | |
814 | for (int i = 0; i < int(iParton.size()); ++i) { | |
815 | if (iParton[i] < 0) legBeg[++leg] = i + 1; | |
816 | else legEnd[leg] = i; | |
817 | } | |
818 | ||
819 | // Iterate from system rest frame towards the junction rest frame (JRF). | |
820 | RotBstMatrix MtoJRF, Mstep; | |
821 | MtoJRF.bstback(pSum); | |
822 | Vec4 pWTinJRF[3]; | |
823 | int iter = 0; | |
824 | double errInCM = 0.; | |
825 | do { | |
826 | ++iter; | |
827 | ||
828 | // Find weighted sum of momenta on the three sides of the junction. | |
829 | for (leg = 0; leg < 3; ++ leg) { | |
830 | pWTinJRF[leg] = 0.; | |
831 | double eWeight = 0.; | |
832 | for (int i = legBeg[leg]; i <= legEnd[leg]; ++i) { | |
833 | Vec4 pTemp = event[ iParton[i] ].p(); | |
834 | pTemp.rotbst(MtoJRF); | |
835 | pWTinJRF[leg] += pTemp * exp(-eWeight); | |
836 | eWeight += pTemp.e() / eNormJunction; | |
837 | if (eWeight > EJNWEIGHTMAX) break; | |
838 | } | |
839 | } | |
840 | ||
841 | // Store original deviation from 120 degree topology. | |
842 | if (iter == 1) errInCM = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5) | |
843 | + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5) | |
844 | + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5); | |
845 | ||
846 | // Find new JRF from the set of weighted momenta. | |
847 | Mstep = junctionRestFrame( pWTinJRF[0], pWTinJRF[1], pWTinJRF[2]); | |
848 | // Fortran code will not take full step after the first few | |
849 | // iterations. How implement this in terms of an M matrix?? | |
850 | MtoJRF.rotbst( Mstep ); | |
851 | } while (iter < 3 || (Mstep.deviation() > CONVJNREST && iter < NTRYJNREST) ); | |
852 | ||
853 | // If final deviation from 120 degrees is bigger than in CM then revert. | |
854 | double errInJRF = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5) | |
855 | + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5) | |
856 | + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5); | |
857 | if (errInJRF > errInCM) { | |
858 | infoPtr->errorMsg("Warning in StringFragmentation::fragmentTo" | |
859 | "Junction: bad convergence junction rest frame"); | |
860 | MtoJRF.reset(); | |
861 | MtoJRF.bstback(pSum); | |
862 | } | |
863 | ||
864 | // Opposite operation: boost from JRF to original system. | |
865 | RotBstMatrix MfromJRF = MtoJRF; | |
866 | MfromJRF.invert(); | |
867 | ||
868 | // Sum leg four-momenta in original frame and in JRF. | |
869 | Vec4 pInLeg[3], pInJRF[3]; | |
870 | for (leg = 0; leg < 3; ++leg) { | |
871 | pInLeg[leg] = 0.; | |
872 | for (int i = legBeg[leg]; i <= legEnd[leg]; ++i) | |
873 | pInLeg[leg] += event[ iParton[i] ].p(); | |
874 | pInJRF[leg] = pInLeg[leg]; | |
875 | pInJRF[leg].rotbst(MtoJRF); | |
876 | } | |
877 | ||
878 | // Pick the two legs with lowest energy in JRF. | |
879 | int legMin = (pInJRF[0].e() < pInJRF[1].e()) ? 0 : 1; | |
880 | int legMax = 1 - legMin; | |
881 | if (pInJRF[2].e() < min(pInJRF[0].e(), pInJRF[1].e()) ) legMin = 2; | |
882 | else if (pInJRF[2].e() > max(pInJRF[0].e(), pInJRF[1].e()) ) legMax = 2; | |
883 | int legMid = 3 - legMin - legMax; | |
884 | ||
885 | // Save info on which status codes belong with the three legs. | |
886 | int iJunction = (-iParton[0]) / 10 - 1; | |
887 | event.statusJunction( iJunction, legMin, 85); | |
888 | event.statusJunction( iJunction, legMid, 86); | |
889 | event.statusJunction( iJunction, legMax, 83); | |
890 | ||
891 | // Temporarily copy the partons on the low-energy legs, into the JRF, | |
892 | // in reverse order, so (anti)quark leg end first. | |
893 | vector<int> iPartonMin; | |
894 | for (int i = legEnd[legMin]; i >= legBeg[legMin]; --i) { | |
895 | int iNew = event.append( event[ iParton[i] ] ); | |
896 | event[iNew].rotbst(MtoJRF); | |
897 | iPartonMin.push_back( iNew ); | |
898 | } | |
899 | vector<int> iPartonMid; | |
900 | for (int i = legEnd[legMid]; i >= legBeg[legMid]; --i) { | |
901 | int iNew = event.append( event[ iParton[i] ] ); | |
902 | event[iNew].rotbst(MtoJRF); | |
903 | iPartonMid.push_back( iNew ); | |
904 | } | |
905 | ||
906 | // Find final weighted sum of momenta on each of the two legs. | |
907 | double eWeight = 0.; | |
908 | pWTinJRF[legMin] = 0.; | |
909 | for (int i = iPartonMin.size() - 1; i >= 0; --i) { | |
910 | pWTinJRF[legMin] += event[ iPartonMin[i] ].p() * exp(-eWeight); | |
911 | eWeight += event[ iPartonMin[i] ].e() / eNormJunction; | |
912 | if (eWeight > EJNWEIGHTMAX) break; | |
913 | } | |
914 | eWeight = 0.; | |
915 | pWTinJRF[legMid] = 0.; | |
916 | for (int i = iPartonMid.size() - 1; i >= 0; --i) { | |
917 | pWTinJRF[legMid] += event[ iPartonMid[i] ].p() * exp(-eWeight); | |
918 | eWeight += event[ iPartonMid[i] ].e() / eNormJunction; | |
919 | if (eWeight > EJNWEIGHTMAX) break; | |
920 | } | |
921 | ||
922 | // Define fictitious opposing partons in JRF and store as string ends. | |
923 | Vec4 pOppose = pWTinJRF[legMin]; | |
924 | pOppose.flip3(); | |
925 | int idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1; | |
926 | if (event[ iPartonMin[0] ].col() > 0) idOppose = -idOppose; | |
927 | int iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0, | |
928 | pOppose, 0.); | |
929 | iPartonMin.push_back( iOppose); | |
930 | pOppose = pWTinJRF[legMid]; | |
931 | pOppose.flip3(); | |
932 | idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1; | |
933 | if (event[ iPartonMid[0] ].col() > 0) idOppose = -idOppose; | |
934 | iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0, | |
935 | pOppose, 0.); | |
936 | iPartonMid.push_back( iOppose); | |
937 | ||
938 | // Set up kinematics of string evolution in low-energy temporary systems. | |
939 | systemMin.setUp(iPartonMin, event); | |
940 | systemMid.setUp(iPartonMid, event); | |
941 | ||
942 | // Outer fallback loop, when too little energy left for third leg. | |
943 | int idMin = 0; | |
944 | int idMid = 0; | |
945 | Vec4 pDiquark; | |
946 | for ( int iTryOuter = 0; ; ++iTryOuter) { | |
947 | ||
948 | // Middle fallback loop, when much unused energy in leg remnants. | |
949 | double eLeftMin = 0.; | |
950 | double eLeftMid = 0.; | |
951 | for ( int iTryMiddle = 0; ; ++iTryMiddle) { | |
952 | ||
953 | // Loop over the two lowest-energy legs. | |
954 | for (int legLoop = 0; legLoop < 2; ++ legLoop) { | |
955 | int legNow = (legLoop == 0) ? legMin : legMid; | |
956 | ||
957 | // Read in properties specific to this leg. | |
958 | StringSystem& systemNow = (legLoop == 0) ? systemMin : systemMid; | |
959 | int idPos = (legLoop == 0) ? event[ iPartonMin[0] ].id() | |
960 | : event[ iPartonMid[0] ].id(); | |
961 | idOppose = (legLoop == 0) ? event[ iPartonMin.back() ].id() | |
962 | : event[ iPartonMid.back() ].id(); | |
963 | double eInJRF = pInJRF[legNow].e(); | |
964 | int statusHad = (legLoop == 0) ? 85 : 86; | |
965 | ||
966 | // Inner fallback loop, when a diquark comes in to junction. | |
967 | double eUsed = 0.; | |
968 | for ( int iTryInner = 0; ; ++iTryInner) { | |
969 | if (iTryInner > NTRYJNMATCH) { | |
970 | infoPtr->errorMsg("Error in StringFragmentation::fragment" | |
971 | "ToJunction: caught in junction flavour loop"); | |
972 | event.popBack( iPartonMin.size() + iPartonMid.size() ); | |
973 | return false; | |
974 | } | |
975 | ||
976 | // Set up two string ends, and begin fragmentation loop. | |
977 | setStartEnds(idPos, idOppose, systemNow); | |
978 | eUsed = 0.; | |
979 | int nHadrons = 0; | |
980 | bool noNegE = true; | |
981 | for ( ; ; ++nHadrons) { | |
982 | ||
983 | // Construct trial hadron from positive end. | |
984 | posEnd.newHadron(); | |
985 | Vec4 pHad = posEnd.kinematicsHadron(systemNow); | |
986 | ||
987 | // Negative energy signals failure in construction. | |
988 | if (pHad.e() < 0. ) { noNegE = false; break; } | |
989 | ||
990 | // Break if passed system midpoint ( = junction) in energy. | |
991 | if (eUsed + pHad.e() > eInJRF) break; | |
992 | ||
993 | // Else construct kinematics of the new hadron and store it. | |
994 | hadrons.append( posEnd.idHad, statusHad, iPos, iNeg, | |
995 | 0, 0, 0, 0, pHad, posEnd.mHad); | |
996 | ||
997 | // Update string end and remaining momentum. | |
998 | posEnd.update(); | |
999 | eUsed += pHad.e(); | |
1000 | } | |
1001 | ||
1002 | // End of fragmentation loop. Inner loopback if ends on a diquark. | |
1003 | if ( noNegE && abs(posEnd.flavOld.id) < 10 ) break; | |
1004 | hadrons.popBack(nHadrons); | |
1005 | } | |
1006 | ||
1007 | // End of one-leg fragmentation. Store end quark and remnant energy. | |
1008 | if (legNow == legMin) { | |
1009 | idMin = posEnd.flavOld.id; | |
1010 | eLeftMin = eInJRF - eUsed; | |
1011 | } else { | |
1012 | idMid = posEnd.flavOld.id; | |
1013 | eLeftMid = eInJRF - eUsed; | |
1014 | } | |
1015 | } | |
1016 | ||
1017 | // End of both-leg fragmentation. Middle loopback if too much energy left. | |
1018 | double eTrial = eBothLeftJunction + rndmPtr->flat() * eMaxLeftJunction; | |
1019 | if (iTryMiddle > NTRYJNMATCH | |
1020 | || ( min( eLeftMin, eLeftMid) < eBothLeftJunction | |
1021 | && max( eLeftMin, eLeftMid) < eTrial ) ) break; | |
1022 | hadrons.clear(); | |
1023 | } | |
1024 | ||
1025 | // Boost hadrons away from the JRF to the original frame. | |
1026 | for (int i = 0; i < hadrons.size(); ++i) { | |
1027 | hadrons[i].rotbst(MfromJRF); | |
1028 | // Recalculate energy to compensate for numerical precision loss | |
1029 | // in iterative calculation of MfromJRF. | |
1030 | hadrons[i].e( hadrons[i].eCalc() ); | |
1031 | pJunctionHadrons += hadrons[i].p(); | |
1032 | } | |
1033 | ||
1034 | // Outer loopback if too little energy left in third leg. | |
1035 | pDiquark = pInLeg[legMin] + pInLeg[legMid] - pJunctionHadrons; | |
1036 | double m2Left = m2( pInLeg[legMax], pDiquark); | |
1037 | if (iTryOuter > NTRYJNMATCH | |
1038 | || m2Left > eMinLeftJunction * pInLeg[legMax].e() ) break; | |
1039 | hadrons.clear(); | |
1040 | pJunctionHadrons = 0.; | |
1041 | } | |
1042 | ||
1043 | // Now found solution; no more loopback. Remove temporary parton copies. | |
1044 | event.popBack( iPartonMin.size() + iPartonMid.size() ); | |
1045 | ||
1046 | // Construct and store an effective diquark string end from the | |
1047 | // two remnant quark ends, for temporary usage. | |
1048 | int idDiquark = flavSelPtr->makeDiquark( idMin, idMid); | |
1049 | double mDiquark = pDiquark.mCalc(); | |
1050 | int iDiquark = event.append( idDiquark, 78, 0, 0, 0, 0, 0, 0, | |
1051 | pDiquark, mDiquark); | |
1052 | ||
1053 | // Find the partons on the last leg, again in reverse order. | |
1054 | vector<int> iPartonMax; | |
1055 | for (int i = legEnd[legMax]; i >= legBeg[legMax]; --i) | |
1056 | iPartonMax.push_back( iParton[i] ); | |
1057 | iPartonMax.push_back( iDiquark ); | |
1058 | ||
1059 | // Modify parton list to remaining leg + remnant of the first two. | |
1060 | iParton = iPartonMax; | |
1061 | ||
1062 | // Done. | |
1063 | return true; | |
1064 | } | |
1065 | ||
1066 | //-------------------------------------------------------------------------- | |
1067 | ||
1068 | // Find the boost matrix to the rest frame of a junction, | |
1069 | // given the three respective endpoint four-momenta. | |
1070 | ||
1071 | RotBstMatrix StringFragmentation::junctionRestFrame(Vec4& p0, Vec4& p1, | |
1072 | Vec4& p2) { | |
1073 | ||
1074 | // Calculate masses and other invariants. | |
1075 | Vec4 pSumJun = p0 + p1 + p2; | |
1076 | double sHat = pSumJun.m2Calc(); | |
1077 | double pp[3][3]; | |
1078 | pp[0][0] = p0.m2Calc(); | |
1079 | pp[1][1] = p1.m2Calc(); | |
1080 | pp[2][2] = p2.m2Calc(); | |
1081 | pp[0][1] = pp[1][0] = p0 * p1; | |
1082 | pp[0][2] = pp[2][0] = p0 * p2; | |
1083 | pp[1][2] = pp[2][1] = p1 * p2; | |
1084 | ||
1085 | // Requirement (eiMax)_j = pi*pj/mj < (eiMax)_k = pi*pk/mk, used below, | |
1086 | // here rewritten as pi*pj * mk < pi*pk * mj and squared. | |
1087 | double eMax01 = pow2(pp[0][1]) * pp[2][2]; | |
1088 | double eMax02 = pow2(pp[0][2]) * pp[1][1]; | |
1089 | double eMax12 = pow2(pp[1][2]) * pp[0][0]; | |
1090 | ||
1091 | // Initially pick i to be the most massive parton. but allow other tries. | |
1092 | int i = (pp[1][1] > pp[0][0]) ? 1 : 0; | |
1093 | if (pp[2][2] > max(pp[0][0], pp[1][1])) i = 2; | |
1094 | int j, k; | |
1095 | double ei = 0.; | |
1096 | double ej = 0.; | |
1097 | double ek = 0.; | |
1098 | for (int iTry = 0; iTry < 3; ++iTry) { | |
1099 | ||
1100 | // Pick j to give minimal eiMax, and k the third vector. | |
1101 | if (i == 0) j = (eMax02 < eMax01) ? 2 : 1; | |
1102 | else if (i == 1) j = (eMax12 < eMax01) ? 2 : 0; | |
1103 | else j = (eMax12 < eMax02) ? 1 : 0; | |
1104 | k = 3 - i - j; | |
1105 | ||
1106 | // Alternative names according to i, j, k conventions. | |
1107 | double m2i = pp[i][i]; | |
1108 | double m2j = pp[j][j]; | |
1109 | double m2k = pp[k][k]; | |
1110 | double pipj = pp[i][j]; | |
1111 | double pipk = pp[i][k]; | |
1112 | double pjpk = pp[j][k]; | |
1113 | ||
1114 | // Trivial to find new parton energies if all three partons are massless. | |
1115 | if (m2i < M2MAXJRF) { | |
1116 | ei = sqrt( 2. * pipk * pipj / (3. * pjpk) ); | |
1117 | ej = sqrt( 2. * pjpk * pipj / (3. * pipk) ); | |
1118 | ek = sqrt( 2. * pipk * pjpk / (3. * pipj) ); | |
1119 | ||
1120 | // Else find three-momentum range for parton i and values at extremes. | |
1121 | } else { | |
1122 | // Minimum when i is at rest. | |
1123 | double piMin = 0.; | |
1124 | double eiMin = sqrt(m2i); | |
1125 | double ejMin = pipj / eiMin; | |
1126 | double ekMin = pipk / eiMin; | |
1127 | double pjMin = sqrtpos( ejMin*ejMin - m2j ); | |
1128 | double pkMin = sqrtpos( ekMin*ekMin - m2k ); | |
1129 | double fMin = ejMin * ekMin + 0.5 * pjMin * pkMin - pjpk; | |
1130 | // Maximum estimated when j + k is at rest, alternatively j at rest. | |
1131 | double eiMax = (pipj + pipk) / sqrt(m2j + m2k + 2. * pjpk); | |
1132 | if (m2j > M2MAXJRF) eiMax = min( eiMax, pipj / sqrt(m2j) ); | |
1133 | double piMax = sqrtpos( eiMax*eiMax - m2i ); | |
1134 | double temp = eiMax*eiMax - 0.25 *piMax*piMax; | |
1135 | double pjMax = (eiMax * sqrtpos( pipj*pipj - m2j * temp ) | |
1136 | - 0.5 * piMax * pipj) / temp; | |
1137 | double pkMax = (eiMax * sqrtpos( pipk*pipk - m2k * temp ) | |
1138 | - 0.5 * piMax * pipk) / temp; | |
1139 | double ejMax = sqrt(pjMax*pjMax + m2j); | |
1140 | double ekMax = sqrt(pkMax*pkMax + m2k); | |
1141 | double fMax = ejMax * ekMax + 0.5 * pjMax * pkMax - pjpk; | |
1142 | ||
1143 | // If unexpected values at upper endpoint then pick another parton. | |
1144 | if (fMax > 0.) { | |
1145 | int iPrel = (i + 1)%3; | |
1146 | if (pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;} | |
1147 | ++iTry; | |
1148 | iPrel = (i + 2)%3; | |
1149 | if (iTry < 3 && pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;} | |
1150 | } | |
1151 | ||
1152 | // Start binary + linear search to find solution inside range. | |
1153 | int iterMin = 0; | |
1154 | int iterMax = 0; | |
1155 | double pi = 0.5 * (piMin + piMax); | |
1156 | for (int iter = 0; iter < NTRYJRFEQ; ++iter) { | |
1157 | ||
1158 | // Derive momentum of other two partons and distance to root. | |
1159 | ei = sqrt(pi*pi + m2i); | |
1160 | temp = ei*ei - 0.25 * pi*pi; | |
1161 | double pj = (ei * sqrtpos( pipj*pipj - m2j * temp ) | |
1162 | - 0.5 * pi * pipj) / temp; | |
1163 | double pk = (ei * sqrtpos( pipk*pipk - m2k * temp ) | |
1164 | - 0.5 * pi * pipk) / temp; | |
1165 | ej = sqrt(pj*pj + m2j); | |
1166 | ek = sqrt(pk*pk + m2k); | |
1167 | double fNow = ej * ek + 0.5 * pj * pk - pjpk; | |
1168 | ||
1169 | // Replace lower or upper bound by new value. | |
1170 | if (fNow > 0.) { ++iterMin; piMin = pi; fMin = fNow;} | |
1171 | else {++iterMax; piMax = pi; fMax = fNow;} | |
1172 | ||
1173 | // Pick next i momentum to explore, hopefully closer to root. | |
1174 | if (2 * iter < NTRYJRFEQ | |
1175 | && (iterMin < 2 || iterMax < 2 || 4 * iter < NTRYJRFEQ)) | |
1176 | { pi = 0.5 * (piMin + piMax); continue;} | |
1177 | if (fMin < 0. || fMax > 0. || abs(fNow) < CONVJRFEQ * sHat) break; | |
1178 | pi = piMin + (piMax - piMin) * fMin / (fMin - fMax); | |
1179 | } | |
1180 | ||
1181 | // If arrived here then either succeeded or exhausted possibilities. | |
1182 | } break; | |
1183 | } | |
1184 | ||
1185 | // Now we know the energies in the junction rest frame. | |
1186 | double eNew[3]; eNew[i] = ei; eNew[j] = ej; eNew[k] = ek; | |
1187 | ||
1188 | // Boost (copy of) partons to their rest frame. | |
1189 | RotBstMatrix Mmove; | |
1190 | Vec4 p0cm = p0; | |
1191 | Vec4 p1cm = p1; | |
1192 | Vec4 p2cm = p2; | |
1193 | Mmove.bstback(pSumJun); | |
1194 | p0cm.rotbst(Mmove); | |
1195 | p1cm.rotbst(Mmove); | |
1196 | p2cm.rotbst(Mmove); | |
1197 | ||
1198 | // Construct difference vectors and the boost to junction rest frame. | |
1199 | Vec4 pDir01 = p0cm / p0cm.e() - p1cm / p1cm.e(); | |
1200 | Vec4 pDir02 = p0cm / p0cm.e() - p2cm / p2cm.e(); | |
1201 | double pDiff01 = pDir01.pAbs2(); | |
1202 | double pDiff02 = pDir02.pAbs2(); | |
1203 | double pDiff0102 = dot3(pDir01, pDir02); | |
1204 | double eDiff01 = eNew[0] / p0cm.e() - eNew[1] / p1cm.e(); | |
1205 | double eDiff02 = eNew[0] / p0cm.e() - eNew[2] / p2cm.e(); | |
1206 | double denom = pDiff01 * pDiff02 - pDiff0102*pDiff0102; | |
1207 | double coef01 = (eDiff01 * pDiff02 - eDiff02 * pDiff0102) / denom; | |
1208 | double coef02 = (eDiff02 * pDiff01 - eDiff01 * pDiff0102) / denom; | |
1209 | Vec4 vJunction = coef01 * pDir01 + coef02 * pDir02; | |
1210 | vJunction.e( sqrt(1. + vJunction.pAbs2()) ); | |
1211 | ||
1212 | // Add two boosts, giving final result. | |
1213 | Mmove.bst(vJunction); | |
1214 | return Mmove; | |
1215 | ||
1216 | } | |
1217 | ||
1218 | //========================================================================== | |
1219 | ||
1220 | } // end namespace Pythia8 |