1 | // StringFragmentation.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 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 = 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-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 | 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 | // Iterate junction rest frame equation until convergence or too many tries. |
---|
308 | const double StringFragmentation::CONVJRFEQ = 1e-12; |
---|
309 | const int StringFragmentation::NTRYJRFEQ = 40; |
---|
310 | |
---|
311 | //-------------------------------------------------------------------------- |
---|
312 | |
---|
313 | // Initialize and save pointers. |
---|
314 | |
---|
315 | void StringFragmentation::init(Info* infoPtrIn, Settings& settings, |
---|
316 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, StringFlav* flavSelPtrIn, |
---|
317 | StringPT* pTSelPtrIn, StringZ* zSelPtrIn) { |
---|
318 | |
---|
319 | // Save pointers. |
---|
320 | infoPtr = infoPtrIn; |
---|
321 | particleDataPtr = particleDataPtrIn; |
---|
322 | rndmPtr = rndmPtrIn; |
---|
323 | flavSelPtr = flavSelPtrIn; |
---|
324 | pTSelPtr = pTSelPtrIn; |
---|
325 | zSelPtr = zSelPtrIn; |
---|
326 | |
---|
327 | // Initialize the StringFragmentation class. |
---|
328 | stopMass = zSelPtr->stopMass(); |
---|
329 | stopNewFlav = zSelPtr->stopNewFlav(); |
---|
330 | stopSmear = zSelPtr->stopSmear(); |
---|
331 | eNormJunction = settings.parm("StringFragmentation:eNormJunction"); |
---|
332 | eBothLeftJunction |
---|
333 | = settings.parm("StringFragmentation:eBothLeftJunction"); |
---|
334 | eMaxLeftJunction |
---|
335 | = settings.parm("StringFragmentation:eMaxLeftJunction"); |
---|
336 | eMinLeftJunction |
---|
337 | = settings.parm("StringFragmentation:eMinLeftJunction"); |
---|
338 | |
---|
339 | // Joining of nearby partons along the string. |
---|
340 | mJoin = settings.parm("FragmentationSystems:mJoin"); |
---|
341 | |
---|
342 | // Initialize the b parameter of the z spectrum, used when joining jets. |
---|
343 | bLund = zSelPtr->bAreaLund(); |
---|
344 | |
---|
345 | // Initialize the hadrons instance of an event record. |
---|
346 | hadrons.init( "(string fragmentation)", particleDataPtr); |
---|
347 | |
---|
348 | // Send on pointers to the two StringEnd instances. |
---|
349 | posEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr); |
---|
350 | negEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr); |
---|
351 | |
---|
352 | } |
---|
353 | |
---|
354 | //-------------------------------------------------------------------------- |
---|
355 | |
---|
356 | // Perform the fragmentation. |
---|
357 | |
---|
358 | bool StringFragmentation::fragment( int iSub, ColConfig& colConfig, |
---|
359 | Event& event) { |
---|
360 | |
---|
361 | // Find partons and their total four-momentum. |
---|
362 | iParton = colConfig[iSub].iParton; |
---|
363 | iPos = iParton[0]; |
---|
364 | if (iPos < 0) iPos = iParton[1]; |
---|
365 | int idPos = event[iPos].id(); |
---|
366 | iNeg = iParton.back(); |
---|
367 | int idNeg = event[iNeg].id(); |
---|
368 | pSum = colConfig[iSub].pSum; |
---|
369 | |
---|
370 | // Reset the local event record. |
---|
371 | hadrons.clear(); |
---|
372 | |
---|
373 | // For closed gluon string: pick first breakup region. |
---|
374 | isClosed = colConfig[iSub].isClosed; |
---|
375 | if (isClosed) iParton = findFirstRegion(iParton, event); |
---|
376 | |
---|
377 | // For junction topology: fragment off two of the string legs. |
---|
378 | // Then iParton overwritten to remaining leg + leftover diquark. |
---|
379 | pJunctionHadrons = 0.; |
---|
380 | hasJunction = colConfig[iSub].hasJunction; |
---|
381 | if (hasJunction && !fragmentToJunction(event)) return false; |
---|
382 | int junctionHadrons = hadrons.size(); |
---|
383 | if (hasJunction) { |
---|
384 | idPos = event[ iParton[0] ].id(); |
---|
385 | idNeg = event.back().id(); |
---|
386 | pSum -= pJunctionHadrons; |
---|
387 | } |
---|
388 | |
---|
389 | // Set up kinematics of string evolution ( = motion). |
---|
390 | system.setUp(iParton, event); |
---|
391 | stopMassNow = stopMass; |
---|
392 | int nExtraJoin = 0; |
---|
393 | |
---|
394 | // Fallback loop, when joining in the middle fails. Bailout if stuck. |
---|
395 | for ( int iTry = 0; ; ++iTry) { |
---|
396 | if (iTry > NTRYJOIN) { |
---|
397 | infoPtr->errorMsg("Error in StringFragmentation::fragment: " |
---|
398 | "stuck in joining"); |
---|
399 | if (hasJunction) ++nExtraJoin; |
---|
400 | if (nExtraJoin > 0) event.popBack(nExtraJoin); |
---|
401 | return false; |
---|
402 | } |
---|
403 | |
---|
404 | // After several failed tries join some (extra) nearby partons. |
---|
405 | if (iTry == NTRYJOIN / 3) nExtraJoin = extraJoin( 2., event); |
---|
406 | if (iTry == 2 * NTRYJOIN / 3) nExtraJoin += extraJoin( 4., event); |
---|
407 | |
---|
408 | // After several failed tries gradually allow larger stop mass. |
---|
409 | if (iTry > NTRYJOIN - NSTOPMASS) stopMassNow *= FACSTOPMASS; |
---|
410 | |
---|
411 | // Set up flavours of two string ends, and reset other info. |
---|
412 | setStartEnds(idPos, idNeg, system); |
---|
413 | pRem = pSum; |
---|
414 | |
---|
415 | // Begin fragmentation loop, interleaved from the two ends. |
---|
416 | bool fromPos; |
---|
417 | for ( ; ; ) { |
---|
418 | |
---|
419 | // Take a step either from the positive or the negative end. |
---|
420 | fromPos = (rndmPtr->flat() < 0.5); |
---|
421 | StringEnd& nowEnd = (fromPos) ? posEnd : negEnd; |
---|
422 | |
---|
423 | // Construct trial hadron and check that energy remains. |
---|
424 | nowEnd.newHadron(); |
---|
425 | if ( energyUsedUp(fromPos) ) break; |
---|
426 | |
---|
427 | // Construct kinematics of the new hadron and store it. |
---|
428 | Vec4 pHad = nowEnd.kinematicsHadron(system); |
---|
429 | int statusHad = (fromPos) ? 83 : 84; |
---|
430 | hadrons.append( nowEnd.idHad, statusHad, iPos, iNeg, |
---|
431 | 0, 0, 0, 0, pHad, nowEnd.mHad); |
---|
432 | if (pHad.e() < 0.) break; |
---|
433 | |
---|
434 | // Update string end and remaining momentum. |
---|
435 | nowEnd.update(); |
---|
436 | pRem -= pHad; |
---|
437 | |
---|
438 | // End of fragmentation loop. |
---|
439 | } |
---|
440 | |
---|
441 | // When done, join in the middle. If this works, then really done. |
---|
442 | if ( finalTwo(fromPos) ) break; |
---|
443 | |
---|
444 | // Else remove produced particles (except from first two junction legs) |
---|
445 | // and start all over. |
---|
446 | int newHadrons = hadrons.size() - junctionHadrons; |
---|
447 | hadrons.popBack(newHadrons); |
---|
448 | } |
---|
449 | |
---|
450 | // Junctions & extra joins: remove fictitious end, restore original partons. |
---|
451 | if (hasJunction) ++nExtraJoin; |
---|
452 | if (nExtraJoin > 0) { |
---|
453 | event.popBack(nExtraJoin); |
---|
454 | iParton = colConfig[iSub].iParton; |
---|
455 | } |
---|
456 | |
---|
457 | // Store the hadrons in the normal event record, ordered from one end. |
---|
458 | store(event); |
---|
459 | |
---|
460 | // Done. |
---|
461 | return true; |
---|
462 | |
---|
463 | } |
---|
464 | |
---|
465 | //-------------------------------------------------------------------------- |
---|
466 | |
---|
467 | // Find region where to put first string break for closed gluon loop. |
---|
468 | |
---|
469 | vector<int> StringFragmentation::findFirstRegion(vector<int>& iPartonIn, |
---|
470 | Event& event) { |
---|
471 | |
---|
472 | // Evaluate mass-squared for all adjacent gluon pairs. |
---|
473 | vector<double> m2Pair; |
---|
474 | double m2Sum = 0.; |
---|
475 | int size = iPartonIn.size(); |
---|
476 | for (int i = 0; i < size; ++i) { |
---|
477 | double m2Now = 0.5 * event[ iPartonIn[i] ].p() |
---|
478 | * event[ iPartonIn[(i + 1)%size] ].p(); |
---|
479 | m2Pair.push_back(m2Now); |
---|
480 | m2Sum += m2Now; |
---|
481 | } |
---|
482 | |
---|
483 | // Pick breakup region with probability proportional to mass-squared. |
---|
484 | double m2Reg = m2Sum * rndmPtr->flat(); |
---|
485 | int iReg = -1; |
---|
486 | do m2Reg -= m2Pair[++iReg]; |
---|
487 | while (m2Reg > 0. && iReg < size - 1); |
---|
488 | |
---|
489 | // Create reordered parton list, with breakup string region duplicated. |
---|
490 | vector<int> iPartonOut; |
---|
491 | for (int i = 0; i < size + 2; ++i) |
---|
492 | iPartonOut.push_back( iPartonIn[(i + iReg)%size] ); |
---|
493 | |
---|
494 | // Done. |
---|
495 | return iPartonOut; |
---|
496 | |
---|
497 | } |
---|
498 | |
---|
499 | //-------------------------------------------------------------------------- |
---|
500 | |
---|
501 | // Set flavours and momentum position for initial string endpoints. |
---|
502 | |
---|
503 | void StringFragmentation::setStartEnds( int idPos, int idNeg, |
---|
504 | StringSystem systemNow) { |
---|
505 | |
---|
506 | // Variables characterizing string endpoints: defaults for open string. |
---|
507 | double px = 0.; |
---|
508 | double py = 0.; |
---|
509 | double Gamma = 0.; |
---|
510 | double xPosFromPos = 1.; |
---|
511 | double xNegFromPos = 0.; |
---|
512 | double xPosFromNeg = 0.; |
---|
513 | double xNegFromNeg = 1.; |
---|
514 | |
---|
515 | // For closed gluon loop need to pick an initial flavour. |
---|
516 | if (isClosed) { |
---|
517 | do { |
---|
518 | int idTry = flavSelPtr->pickLightQ(); |
---|
519 | FlavContainer flavTry(idTry, 1); |
---|
520 | flavTry = flavSelPtr->pick( flavTry); |
---|
521 | flavTry = flavSelPtr->pick( flavTry); |
---|
522 | idPos = flavTry.id; |
---|
523 | idNeg = -idPos; |
---|
524 | } while (idPos == 0); |
---|
525 | |
---|
526 | // Also need pT and breakup vertex position in region. |
---|
527 | pair<double, double> pxy = pTSelPtr->pxy(); |
---|
528 | px = pxy.first; |
---|
529 | py = pxy.second; |
---|
530 | double m2Region = systemNow.regionLowPos(0).w2; |
---|
531 | double m2Temp = min( CLOSEDM2MAX, CLOSEDM2FRAC * m2Region); |
---|
532 | do { |
---|
533 | double zTemp = zSelPtr->zFrag( idPos, idNeg, m2Temp); |
---|
534 | xPosFromPos = 1. - zTemp; |
---|
535 | xNegFromPos = m2Temp / (zTemp * m2Region); |
---|
536 | } while (xNegFromPos > 1.); |
---|
537 | Gamma = xPosFromPos * xNegFromPos * m2Region; |
---|
538 | xPosFromNeg = xPosFromPos; |
---|
539 | xNegFromNeg = xNegFromPos; |
---|
540 | } |
---|
541 | |
---|
542 | // Initialize two string endpoints. |
---|
543 | posEnd.setUp( true, iPos, idPos, systemNow.iMax, px, py, |
---|
544 | Gamma, xPosFromPos, xNegFromPos); |
---|
545 | negEnd.setUp( false, iNeg, idNeg, systemNow.iMax, -px, -py, |
---|
546 | Gamma, xPosFromNeg, xNegFromNeg); |
---|
547 | |
---|
548 | // For closed gluon loop can allow popcorn on one side but not both. |
---|
549 | if (isClosed) { |
---|
550 | flavSelPtr->assignPopQ(posEnd.flavOld); |
---|
551 | flavSelPtr->assignPopQ(negEnd.flavOld); |
---|
552 | if (rndmPtr->flat() < 0.5) posEnd.flavOld.nPop = 0; |
---|
553 | else negEnd.flavOld.nPop = 0; |
---|
554 | posEnd.flavOld.rank = 1; |
---|
555 | negEnd.flavOld.rank = 1; |
---|
556 | } |
---|
557 | |
---|
558 | // Done. |
---|
559 | |
---|
560 | } |
---|
561 | |
---|
562 | //-------------------------------------------------------------------------- |
---|
563 | |
---|
564 | // Check remaining energy-momentum whether it is OK to continue. |
---|
565 | |
---|
566 | bool StringFragmentation::energyUsedUp(bool fromPos) { |
---|
567 | |
---|
568 | // If remaining negative energy then abort right away. |
---|
569 | if (pRem.e() < 0.) return true; |
---|
570 | |
---|
571 | // Calculate W2_minimum and done if remaining W2 is below it. |
---|
572 | double wMin = stopMassNow |
---|
573 | + particleDataPtr->constituentMass(posEnd.flavOld.id) |
---|
574 | + particleDataPtr->constituentMass(negEnd.flavOld.id); |
---|
575 | if (fromPos) wMin += stopNewFlav |
---|
576 | * particleDataPtr->constituentMass(posEnd.flavNew.id); |
---|
577 | else wMin += stopNewFlav |
---|
578 | * particleDataPtr->constituentMass(negEnd.flavNew.id); |
---|
579 | wMin *= 1. + (2. * rndmPtr->flat() - 1.) * stopSmear; |
---|
580 | w2Rem = pRem.m2Calc(); |
---|
581 | if (w2Rem < pow2(wMin)) return true; |
---|
582 | |
---|
583 | // Else still enough energy left to continue iteration. |
---|
584 | return false; |
---|
585 | |
---|
586 | } |
---|
587 | |
---|
588 | //-------------------------------------------------------------------------- |
---|
589 | |
---|
590 | // Produce the final two partons to complete the system. |
---|
591 | |
---|
592 | bool StringFragmentation::finalTwo(bool fromPos) { |
---|
593 | |
---|
594 | // Check whether we went too far in p+-. |
---|
595 | if (pRem.e() < 0. || w2Rem < 0. || (hadrons.size() > 0 |
---|
596 | && hadrons.back().e() < 0.) ) return false; |
---|
597 | if ( posEnd.iPosOld > negEnd.iPosOld || negEnd.iNegOld > posEnd.iNegOld) |
---|
598 | return false; |
---|
599 | if ( posEnd.iPosOld == negEnd.iPosOld && posEnd.xPosOld < negEnd.xPosOld) |
---|
600 | return false; |
---|
601 | if ( posEnd.iNegOld == negEnd.iNegOld && posEnd.xNegOld > negEnd.xNegOld) |
---|
602 | return false; |
---|
603 | |
---|
604 | // Construct the final hadron from the leftover flavours. |
---|
605 | // Impossible to join two diquarks. Also break if stuck for other reason. |
---|
606 | FlavContainer flav1 = (fromPos) ? posEnd.flavNew.anti() : posEnd.flavOld; |
---|
607 | FlavContainer flav2 = (fromPos) ? negEnd.flavOld : negEnd.flavNew.anti(); |
---|
608 | if (flav1.isDiquark() && flav2.isDiquark()) return false; |
---|
609 | int idHad; |
---|
610 | for (int iTry = 0; iTry < NTRYFLAV; ++iTry) { |
---|
611 | idHad = flavSelPtr->combine( flav1, flav2); |
---|
612 | if (idHad != 0) break; |
---|
613 | } |
---|
614 | if (idHad == 0) return false; |
---|
615 | |
---|
616 | // Store the final particle and its new pT, and construct its mass. |
---|
617 | if (fromPos) { |
---|
618 | negEnd.idHad = idHad; |
---|
619 | negEnd.pxNew = -posEnd.pxNew; |
---|
620 | negEnd.pyNew = -posEnd.pyNew; |
---|
621 | negEnd.mHad = particleDataPtr->mass(idHad); |
---|
622 | } else { |
---|
623 | posEnd.idHad = idHad; |
---|
624 | posEnd.pxNew = -negEnd.pxNew; |
---|
625 | posEnd.pyNew = -negEnd.pyNew; |
---|
626 | posEnd.mHad = particleDataPtr->mass(idHad); |
---|
627 | } |
---|
628 | |
---|
629 | // String region in which to do the joining. |
---|
630 | StringRegion region = finalRegion(); |
---|
631 | if (region.isEmpty) return false; |
---|
632 | |
---|
633 | // Project remaining momentum along longitudinal and transverse directions. |
---|
634 | region.project( pRem); |
---|
635 | double pxRem = region.px() - posEnd.pxOld - negEnd.pxOld; |
---|
636 | double pyRem = region.py() - posEnd.pyOld - negEnd.pyOld; |
---|
637 | double xPosRem = region.xPos(); |
---|
638 | double xNegRem = region.xNeg(); |
---|
639 | |
---|
640 | // Share extra pT kick evenly between final two hadrons. |
---|
641 | posEnd.pxOld += 0.5 * pxRem; |
---|
642 | posEnd.pyOld += 0.5 * pyRem; |
---|
643 | negEnd.pxOld += 0.5 * pxRem; |
---|
644 | negEnd.pyOld += 0.5 * pyRem; |
---|
645 | |
---|
646 | // Construct new pT and mT of the final two particles. |
---|
647 | posEnd.pxHad = posEnd.pxOld + posEnd.pxNew; |
---|
648 | posEnd.pyHad = posEnd.pyOld + posEnd.pyNew; |
---|
649 | posEnd.mT2Had = pow2(posEnd.mHad) + pow2(posEnd.pxHad) |
---|
650 | + pow2(posEnd.pyHad); |
---|
651 | negEnd.pxHad = negEnd.pxOld + negEnd.pxNew; |
---|
652 | negEnd.pyHad = negEnd.pyOld + negEnd.pyNew; |
---|
653 | negEnd.mT2Had = pow2(negEnd.mHad) + pow2(negEnd.pxHad) |
---|
654 | + pow2(negEnd.pyHad); |
---|
655 | |
---|
656 | // Construct remaining system transverse mass. |
---|
657 | double wT2Rem = w2Rem + pow2( posEnd.pxHad + negEnd.pxHad) |
---|
658 | + pow2( posEnd.pyHad + negEnd.pyHad); |
---|
659 | |
---|
660 | // Check that kinematics possible. |
---|
661 | if ( sqrt(wT2Rem) < sqrt(posEnd.mT2Had) + sqrt(negEnd.mT2Had) ) |
---|
662 | return false; |
---|
663 | double lambda2 = pow2( wT2Rem - posEnd.mT2Had - negEnd.mT2Had) |
---|
664 | - 4. * posEnd.mT2Had * negEnd.mT2Had; |
---|
665 | if (lambda2 <= 0.) return false; |
---|
666 | |
---|
667 | // Construct kinematics, as viewed in the transverse rest frame. |
---|
668 | double lambda = sqrt(lambda2); |
---|
669 | double probReverse = 1. / (1. + exp( min( EXPMAX, bLund * lambda) ) ); |
---|
670 | double xpzPos = 0.5 * lambda/ wT2Rem; |
---|
671 | if (probReverse > rndmPtr->flat()) xpzPos = -xpzPos; |
---|
672 | double xmDiff = (posEnd.mT2Had - negEnd.mT2Had) / wT2Rem; |
---|
673 | double xePos = 0.5 * (1. + xmDiff); |
---|
674 | double xeNeg = 0.5 * (1. - xmDiff ); |
---|
675 | |
---|
676 | // Translate this into kinematics in the string frame. |
---|
677 | Vec4 pHadPos = region.pHad( (xePos + xpzPos) * xPosRem, |
---|
678 | (xePos - xpzPos) * xNegRem, posEnd.pxHad, posEnd.pyHad); |
---|
679 | Vec4 pHadNeg = region.pHad( (xeNeg - xpzPos) * xPosRem, |
---|
680 | (xeNeg + xpzPos) * xNegRem, negEnd.pxHad, negEnd.pyHad); |
---|
681 | |
---|
682 | // Add produced particles to the event record. |
---|
683 | hadrons.append( posEnd.idHad, 83, posEnd.iEnd, negEnd.iEnd, |
---|
684 | 0, 0, 0, 0, pHadPos, posEnd.mHad); |
---|
685 | hadrons.append( negEnd.idHad, 84, posEnd.iEnd, negEnd.iEnd, |
---|
686 | 0, 0, 0, 0, pHadNeg, negEnd.mHad); |
---|
687 | |
---|
688 | // It worked. |
---|
689 | return true; |
---|
690 | |
---|
691 | } |
---|
692 | |
---|
693 | //-------------------------------------------------------------------------- |
---|
694 | |
---|
695 | // Construct a special joining region for the final two hadrons. |
---|
696 | |
---|
697 | StringRegion StringFragmentation::finalRegion() { |
---|
698 | |
---|
699 | // Simple case when both string ends are in the same region. |
---|
700 | if (posEnd.iPosOld == negEnd.iPosOld && posEnd.iNegOld == negEnd.iNegOld) |
---|
701 | return system.region( posEnd.iPosOld, posEnd.iNegOld); |
---|
702 | |
---|
703 | // Start out with empty region. (Empty used for error returns.) |
---|
704 | StringRegion region; |
---|
705 | |
---|
706 | // Add up all remaining p+. |
---|
707 | Vec4 pPosJoin; |
---|
708 | if ( posEnd.iPosOld == negEnd.iPosOld) { |
---|
709 | double xPosJoin = posEnd.xPosOld - negEnd.xPosOld; |
---|
710 | if (xPosJoin < 0.) return region; |
---|
711 | pPosJoin = system.regionLowPos(posEnd.iPosOld).pHad( xPosJoin, 0., 0., 0.); |
---|
712 | } else { |
---|
713 | for (int iPosNow = posEnd.iPosOld; iPosNow <= negEnd.iPosOld; ++iPosNow) { |
---|
714 | if (iPosNow == posEnd.iPosOld) pPosJoin |
---|
715 | += system.regionLowPos(iPosNow).pHad( posEnd.xPosOld, 0., 0., 0.); |
---|
716 | else if (iPosNow == negEnd.iPosOld) pPosJoin |
---|
717 | += system.regionLowPos(iPosNow).pHad( 1. - negEnd.xPosOld, 0., 0., 0.); |
---|
718 | else pPosJoin += system.regionLowPos(iPosNow).pHad( 1., 0., 0., 0.); |
---|
719 | } |
---|
720 | } |
---|
721 | |
---|
722 | // Add up all remaining p-. |
---|
723 | Vec4 pNegJoin; |
---|
724 | if ( negEnd.iNegOld == posEnd.iNegOld) { |
---|
725 | double xNegJoin = negEnd.xNegOld - posEnd.xNegOld; |
---|
726 | if (xNegJoin < 0.) return region; |
---|
727 | pNegJoin = system.regionLowNeg(negEnd.iNegOld).pHad( 0., xNegJoin, 0., 0.); |
---|
728 | } else { |
---|
729 | for (int iNegNow = negEnd.iNegOld; iNegNow <= posEnd.iNegOld; ++iNegNow) { |
---|
730 | if (iNegNow == negEnd.iNegOld) pNegJoin |
---|
731 | += system.regionLowNeg(iNegNow).pHad( 0., negEnd.xNegOld, 0., 0.); |
---|
732 | else if (iNegNow == posEnd.iNegOld) pNegJoin |
---|
733 | += system.regionLowNeg(iNegNow).pHad( 0., 1. - posEnd.xNegOld, 0., 0.); |
---|
734 | else pNegJoin += system.regionLowNeg(iNegNow).pHad( 0., 1., 0., 0.); |
---|
735 | } |
---|
736 | } |
---|
737 | |
---|
738 | // For a closed gluon loop pPosJoin == pNegJoin and the above does not work. |
---|
739 | // So reshuffle; "perfect" for g g systems, OK in general. |
---|
740 | Vec4 pTest = pPosJoin - pNegJoin; |
---|
741 | if ( abs(pTest.px()) + abs(pTest.py()) + abs(pTest.pz()) + abs(pTest.e()) |
---|
742 | < MATCHPOSNEG * (pPosJoin.e() + pNegJoin.e()) ) { |
---|
743 | Vec4 delta |
---|
744 | = system.regionLowPos(posEnd.iPosOld + 1).pHad( 1., 0., 0., 0.) |
---|
745 | - system.regionLowNeg(negEnd.iNegOld + 1).pHad( 0., 1., 0., 0.); |
---|
746 | pPosJoin -= delta; |
---|
747 | pNegJoin += delta; |
---|
748 | } |
---|
749 | |
---|
750 | // Construct a new region from remaining p+ and p-. |
---|
751 | region.setUp( pPosJoin, pNegJoin); |
---|
752 | if (region.isEmpty) return region; |
---|
753 | |
---|
754 | // Project the existing pTold vectors onto the new directions. |
---|
755 | Vec4 pTposOld = system.region( posEnd.iPosOld, posEnd.iNegOld).pHad( |
---|
756 | 0., 0., posEnd.pxOld, posEnd.pyOld); |
---|
757 | region.project( pTposOld); |
---|
758 | posEnd.pxOld = region.px(); |
---|
759 | posEnd.pyOld = region.py(); |
---|
760 | Vec4 pTnegOld = system.region( negEnd.iPosOld, negEnd.iNegOld).pHad( |
---|
761 | 0., 0., negEnd.pxOld, negEnd.pyOld); |
---|
762 | region.project( pTnegOld); |
---|
763 | negEnd.pxOld = region.px(); |
---|
764 | negEnd.pyOld = region.py(); |
---|
765 | |
---|
766 | // Done. |
---|
767 | return region; |
---|
768 | |
---|
769 | } |
---|
770 | |
---|
771 | //-------------------------------------------------------------------------- |
---|
772 | |
---|
773 | // Store the hadrons in the normal event record, ordered from one end. |
---|
774 | |
---|
775 | void StringFragmentation::store(Event& event) { |
---|
776 | |
---|
777 | // Starting position. |
---|
778 | int iFirst = event.size(); |
---|
779 | |
---|
780 | // Copy straight over from first two junction legs. |
---|
781 | if (hasJunction) { |
---|
782 | for (int i = 0; i < hadrons.size(); ++i) |
---|
783 | if (hadrons[i].status() == 85 || hadrons[i].status() == 86) |
---|
784 | event.append( hadrons[i] ); |
---|
785 | } |
---|
786 | |
---|
787 | // Loop downwards, copying all from the positive end. |
---|
788 | for (int i = 0; i < hadrons.size(); ++i) |
---|
789 | if (hadrons[i].status() == 83) event.append( hadrons[i] ); |
---|
790 | |
---|
791 | // Loop upwards, copying all from the negative end. |
---|
792 | for (int i = hadrons.size() - 1; i >= 0 ; --i) |
---|
793 | if (hadrons[i].status() == 84) event.append( hadrons[i] ); |
---|
794 | int iLast = event.size() - 1; |
---|
795 | |
---|
796 | // Set decay vertex when this is displaced. |
---|
797 | if (event[posEnd.iEnd].hasVertex()) { |
---|
798 | Vec4 vDec = event[posEnd.iEnd].vDec(); |
---|
799 | for (int i = iFirst; i <= iLast; ++i) event[i].vProd( vDec ); |
---|
800 | } |
---|
801 | |
---|
802 | // Set lifetime of hadrons. |
---|
803 | for (int i = iFirst; i <= iLast; ++i) |
---|
804 | event[i].tau( event[i].tau0() * rndmPtr->exp() ); |
---|
805 | |
---|
806 | // Mark original partons as hadronized and set their daughter range. |
---|
807 | for (int i = 0; i < int(iParton.size()); ++i) |
---|
808 | if (iParton[i] >= 0) { |
---|
809 | event[ iParton[i] ].statusNeg(); |
---|
810 | event[ iParton[i] ].daughters(iFirst, iLast); |
---|
811 | } |
---|
812 | |
---|
813 | } |
---|
814 | |
---|
815 | //-------------------------------------------------------------------------- |
---|
816 | |
---|
817 | // Fragment off two of the string legs in to a junction. |
---|
818 | |
---|
819 | bool StringFragmentation::fragmentToJunction(Event& event) { |
---|
820 | |
---|
821 | // Identify range of partons on the three legs. |
---|
822 | // (Each leg begins with an iParton[i] = -(10 + 10*junctionNumber + leg), |
---|
823 | // and partons then appear ordered from the junction outwards.) |
---|
824 | int legBeg[3] = { 0, 0, 0}; |
---|
825 | int legEnd[3] = { 0, 0, 0}; |
---|
826 | int leg = -1; |
---|
827 | // PS (4/10/2011) Protect against invalid systems |
---|
828 | if (iParton[0] > 0) { |
---|
829 | infoPtr->errorMsg("Error in StringFragmentation::fragment" |
---|
830 | "ToJunction: iParton[0] not a valid junctionNumber"); |
---|
831 | return false; |
---|
832 | } |
---|
833 | for (int i = 0; i < int(iParton.size()); ++i) { |
---|
834 | if (iParton[i] < 0) { |
---|
835 | if (leg == 2) { |
---|
836 | infoPtr->errorMsg("Error in StringFragmentation::fragment" |
---|
837 | "ToJunction: unprocessed multi-junction system"); |
---|
838 | return false; |
---|
839 | } |
---|
840 | legBeg[++leg] = i + 1; |
---|
841 | } |
---|
842 | else legEnd[leg] = i; |
---|
843 | } |
---|
844 | |
---|
845 | // Iterate from system rest frame towards the junction rest frame (JRF). |
---|
846 | RotBstMatrix MtoJRF, Mstep; |
---|
847 | MtoJRF.bstback(pSum); |
---|
848 | Vec4 pWTinJRF[3]; |
---|
849 | int iter = 0; |
---|
850 | double errInCM = 0.; |
---|
851 | do { |
---|
852 | ++iter; |
---|
853 | |
---|
854 | // Find weighted sum of momenta on the three sides of the junction. |
---|
855 | for (leg = 0; leg < 3; ++ leg) { |
---|
856 | pWTinJRF[leg] = 0.; |
---|
857 | double eWeight = 0.; |
---|
858 | for (int i = legBeg[leg]; i <= legEnd[leg]; ++i) { |
---|
859 | Vec4 pTemp = event[ iParton[i] ].p(); |
---|
860 | pTemp.rotbst(MtoJRF); |
---|
861 | pWTinJRF[leg] += pTemp * exp(-eWeight); |
---|
862 | eWeight += pTemp.e() / eNormJunction; |
---|
863 | if (eWeight > EJNWEIGHTMAX) break; |
---|
864 | } |
---|
865 | } |
---|
866 | |
---|
867 | // Store original deviation from 120 degree topology. |
---|
868 | if (iter == 1) errInCM = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5) |
---|
869 | + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5) |
---|
870 | + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5); |
---|
871 | |
---|
872 | // Find new JRF from the set of weighted momenta. |
---|
873 | Mstep = junctionRestFrame( pWTinJRF[0], pWTinJRF[1], pWTinJRF[2]); |
---|
874 | // Fortran code will not take full step after the first few |
---|
875 | // iterations. How implement this in terms of an M matrix?? |
---|
876 | MtoJRF.rotbst( Mstep ); |
---|
877 | } while (iter < 3 || (Mstep.deviation() > CONVJNREST && iter < NTRYJNREST) ); |
---|
878 | |
---|
879 | // If final deviation from 120 degrees is bigger than in CM then revert. |
---|
880 | double errInJRF = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5) |
---|
881 | + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5) |
---|
882 | + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5); |
---|
883 | if (errInJRF > errInCM) { |
---|
884 | infoPtr->errorMsg("Warning in StringFragmentation::fragmentTo" |
---|
885 | "Junction: bad convergence junction rest frame"); |
---|
886 | MtoJRF.reset(); |
---|
887 | MtoJRF.bstback(pSum); |
---|
888 | } |
---|
889 | |
---|
890 | // Opposite operation: boost from JRF to original system. |
---|
891 | RotBstMatrix MfromJRF = MtoJRF; |
---|
892 | MfromJRF.invert(); |
---|
893 | |
---|
894 | // Sum leg four-momenta in original frame and in JRF. |
---|
895 | Vec4 pInLeg[3], pInJRF[3]; |
---|
896 | for (leg = 0; leg < 3; ++leg) { |
---|
897 | pInLeg[leg] = 0.; |
---|
898 | for (int i = legBeg[leg]; i <= legEnd[leg]; ++i) |
---|
899 | pInLeg[leg] += event[ iParton[i] ].p(); |
---|
900 | pInJRF[leg] = pInLeg[leg]; |
---|
901 | pInJRF[leg].rotbst(MtoJRF); |
---|
902 | } |
---|
903 | |
---|
904 | // Pick the two legs with lowest energy in JRF. |
---|
905 | int legMin = (pInJRF[0].e() < pInJRF[1].e()) ? 0 : 1; |
---|
906 | int legMax = 1 - legMin; |
---|
907 | if (pInJRF[2].e() < min(pInJRF[0].e(), pInJRF[1].e()) ) legMin = 2; |
---|
908 | else if (pInJRF[2].e() > max(pInJRF[0].e(), pInJRF[1].e()) ) legMax = 2; |
---|
909 | int legMid = 3 - legMin - legMax; |
---|
910 | |
---|
911 | // Save info on which status codes belong with the three legs. |
---|
912 | int iJunction = (-iParton[0]) / 10 - 1; |
---|
913 | event.statusJunction( iJunction, legMin, 85); |
---|
914 | event.statusJunction( iJunction, legMid, 86); |
---|
915 | event.statusJunction( iJunction, legMax, 83); |
---|
916 | |
---|
917 | // Temporarily copy the partons on the low-energy legs, into the JRF, |
---|
918 | // in reverse order, so (anti)quark leg end first. |
---|
919 | vector<int> iPartonMin; |
---|
920 | for (int i = legEnd[legMin]; i >= legBeg[legMin]; --i) { |
---|
921 | int iNew = event.append( event[ iParton[i] ] ); |
---|
922 | event[iNew].rotbst(MtoJRF); |
---|
923 | iPartonMin.push_back( iNew ); |
---|
924 | } |
---|
925 | vector<int> iPartonMid; |
---|
926 | for (int i = legEnd[legMid]; i >= legBeg[legMid]; --i) { |
---|
927 | int iNew = event.append( event[ iParton[i] ] ); |
---|
928 | event[iNew].rotbst(MtoJRF); |
---|
929 | iPartonMid.push_back( iNew ); |
---|
930 | } |
---|
931 | |
---|
932 | // Find final weighted sum of momenta on each of the two legs. |
---|
933 | double eWeight = 0.; |
---|
934 | pWTinJRF[legMin] = 0.; |
---|
935 | for (int i = iPartonMin.size() - 1; i >= 0; --i) { |
---|
936 | pWTinJRF[legMin] += event[ iPartonMin[i] ].p() * exp(-eWeight); |
---|
937 | eWeight += event[ iPartonMin[i] ].e() / eNormJunction; |
---|
938 | if (eWeight > EJNWEIGHTMAX) break; |
---|
939 | } |
---|
940 | eWeight = 0.; |
---|
941 | pWTinJRF[legMid] = 0.; |
---|
942 | for (int i = iPartonMid.size() - 1; i >= 0; --i) { |
---|
943 | pWTinJRF[legMid] += event[ iPartonMid[i] ].p() * exp(-eWeight); |
---|
944 | eWeight += event[ iPartonMid[i] ].e() / eNormJunction; |
---|
945 | if (eWeight > EJNWEIGHTMAX) break; |
---|
946 | } |
---|
947 | |
---|
948 | // Define fictitious opposing partons in JRF and store as string ends. |
---|
949 | Vec4 pOppose = pWTinJRF[legMin]; |
---|
950 | pOppose.flip3(); |
---|
951 | int idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1; |
---|
952 | if (event[ iPartonMin[0] ].col() > 0) idOppose = -idOppose; |
---|
953 | int iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0, |
---|
954 | pOppose, 0.); |
---|
955 | iPartonMin.push_back( iOppose); |
---|
956 | pOppose = pWTinJRF[legMid]; |
---|
957 | pOppose.flip3(); |
---|
958 | idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1; |
---|
959 | if (event[ iPartonMid[0] ].col() > 0) idOppose = -idOppose; |
---|
960 | iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0, |
---|
961 | pOppose, 0.); |
---|
962 | iPartonMid.push_back( iOppose); |
---|
963 | |
---|
964 | // Set up kinematics of string evolution in low-energy temporary systems. |
---|
965 | systemMin.setUp(iPartonMin, event); |
---|
966 | systemMid.setUp(iPartonMid, event); |
---|
967 | |
---|
968 | // Outer fallback loop, when too little energy left for third leg. |
---|
969 | int idMin = 0; |
---|
970 | int idMid = 0; |
---|
971 | Vec4 pDiquark; |
---|
972 | for ( int iTryOuter = 0; ; ++iTryOuter) { |
---|
973 | |
---|
974 | // Middle fallback loop, when much unused energy in leg remnants. |
---|
975 | double eLeftMin = 0.; |
---|
976 | double eLeftMid = 0.; |
---|
977 | for ( int iTryMiddle = 0; ; ++iTryMiddle) { |
---|
978 | |
---|
979 | // Loop over the two lowest-energy legs. |
---|
980 | for (int legLoop = 0; legLoop < 2; ++ legLoop) { |
---|
981 | int legNow = (legLoop == 0) ? legMin : legMid; |
---|
982 | |
---|
983 | // Read in properties specific to this leg. |
---|
984 | StringSystem& systemNow = (legLoop == 0) ? systemMin : systemMid; |
---|
985 | int idPos = (legLoop == 0) ? event[ iPartonMin[0] ].id() |
---|
986 | : event[ iPartonMid[0] ].id(); |
---|
987 | idOppose = (legLoop == 0) ? event[ iPartonMin.back() ].id() |
---|
988 | : event[ iPartonMid.back() ].id(); |
---|
989 | double eInJRF = pInJRF[legNow].e(); |
---|
990 | int statusHad = (legLoop == 0) ? 85 : 86; |
---|
991 | |
---|
992 | // Inner fallback loop, when a diquark comes in to junction. |
---|
993 | double eUsed = 0.; |
---|
994 | for ( int iTryInner = 0; ; ++iTryInner) { |
---|
995 | if (iTryInner > 2 * NTRYJNMATCH) { |
---|
996 | infoPtr->errorMsg("Error in StringFragmentation::fragment" |
---|
997 | "ToJunction: caught in junction flavour loop"); |
---|
998 | event.popBack( iPartonMin.size() + iPartonMid.size() ); |
---|
999 | return false; |
---|
1000 | } |
---|
1001 | bool needBaryon = (abs(idPos) > 10 && iTryInner > NTRYJNMATCH); |
---|
1002 | double eExtra = (iTryInner > NTRYJNMATCH) ? EEXTRAJNMATCH : 0.; |
---|
1003 | |
---|
1004 | // Set up two string ends, and begin fragmentation loop. |
---|
1005 | setStartEnds(idPos, idOppose, systemNow); |
---|
1006 | eUsed = 0.; |
---|
1007 | int nHadrons = 0; |
---|
1008 | bool noNegE = true; |
---|
1009 | for ( ; ; ++nHadrons) { |
---|
1010 | |
---|
1011 | // Construct trial hadron from positive end. |
---|
1012 | posEnd.newHadron(); |
---|
1013 | Vec4 pHad = posEnd.kinematicsHadron(systemNow); |
---|
1014 | |
---|
1015 | // Negative energy signals failure in construction. |
---|
1016 | if (pHad.e() < 0. ) { noNegE = false; break; } |
---|
1017 | |
---|
1018 | // Break if passed system midpoint ( = junction) in energy. |
---|
1019 | // Exceptions: small systems, and/or with diquark end. |
---|
1020 | bool delayedBreak = false; |
---|
1021 | if (eUsed + pHad.e() + eExtra > eInJRF) { |
---|
1022 | if (nHadrons > 0 || !needBaryon) break; |
---|
1023 | delayedBreak = true; |
---|
1024 | } |
---|
1025 | |
---|
1026 | // Else construct kinematics of the new hadron and store it. |
---|
1027 | hadrons.append( posEnd.idHad, statusHad, iPos, iNeg, |
---|
1028 | 0, 0, 0, 0, pHad, posEnd.mHad); |
---|
1029 | |
---|
1030 | // Update string end and remaining momentum. |
---|
1031 | posEnd.update(); |
---|
1032 | eUsed += pHad.e(); |
---|
1033 | |
---|
1034 | // Delayed break in small systems, and/or with diquark end. |
---|
1035 | if (delayedBreak) { |
---|
1036 | ++nHadrons; |
---|
1037 | break; |
---|
1038 | } |
---|
1039 | } |
---|
1040 | |
---|
1041 | // End of fragmentation loop. Inner loopback if ends on a diquark. |
---|
1042 | if ( noNegE && abs(posEnd.flavOld.id) < 10 ) break; |
---|
1043 | hadrons.popBack(nHadrons); |
---|
1044 | } |
---|
1045 | |
---|
1046 | // End of one-leg fragmentation. Store end quark and remnant energy. |
---|
1047 | if (legNow == legMin) { |
---|
1048 | idMin = posEnd.flavOld.id; |
---|
1049 | eLeftMin = eInJRF - eUsed; |
---|
1050 | } else { |
---|
1051 | idMid = posEnd.flavOld.id; |
---|
1052 | eLeftMid = eInJRF - eUsed; |
---|
1053 | } |
---|
1054 | } |
---|
1055 | |
---|
1056 | // End of both-leg fragmentation. |
---|
1057 | // Middle loopback if too much energy left. |
---|
1058 | double eTrial = eBothLeftJunction + rndmPtr->flat() * eMaxLeftJunction; |
---|
1059 | if (iTryMiddle > NTRYJNMATCH |
---|
1060 | || ( min( eLeftMin, eLeftMid) < eBothLeftJunction |
---|
1061 | && max( eLeftMin, eLeftMid) < eTrial ) ) break; |
---|
1062 | hadrons.clear(); |
---|
1063 | } |
---|
1064 | |
---|
1065 | // Boost hadrons away from the JRF to the original frame. |
---|
1066 | for (int i = 0; i < hadrons.size(); ++i) { |
---|
1067 | hadrons[i].rotbst(MfromJRF); |
---|
1068 | // Recalculate energy to compensate for numerical precision loss |
---|
1069 | // in iterative calculation of MfromJRF. |
---|
1070 | hadrons[i].e( hadrons[i].eCalc() ); |
---|
1071 | pJunctionHadrons += hadrons[i].p(); |
---|
1072 | } |
---|
1073 | |
---|
1074 | // Outer loopback if combined diquark mass too negative |
---|
1075 | // or too little energy left in third leg. |
---|
1076 | pDiquark = pInLeg[legMin] + pInLeg[legMid] - pJunctionHadrons; |
---|
1077 | double m2Left = m2( pInLeg[legMax], pDiquark); |
---|
1078 | if (iTryOuter > NTRYJNMATCH || (pDiquark.mCalc() > MDIQUARKMIN |
---|
1079 | && m2Left > eMinLeftJunction * pInLeg[legMax].e()) ) break; |
---|
1080 | hadrons.clear(); |
---|
1081 | pJunctionHadrons = 0.; |
---|
1082 | } |
---|
1083 | |
---|
1084 | // Now found solution; no more loopback. Remove temporary parton copies. |
---|
1085 | event.popBack( iPartonMin.size() + iPartonMid.size() ); |
---|
1086 | |
---|
1087 | // Construct and store an effective diquark string end from the |
---|
1088 | // two remnant quark ends, for temporary usage. |
---|
1089 | int idDiquark = flavSelPtr->makeDiquark( idMin, idMid); |
---|
1090 | double mDiquark = pDiquark.mCalc(); |
---|
1091 | int iDiquark = event.append( idDiquark, 78, 0, 0, 0, 0, 0, 0, |
---|
1092 | pDiquark, mDiquark); |
---|
1093 | |
---|
1094 | // Find the partons on the last leg, again in reverse order. |
---|
1095 | vector<int> iPartonMax; |
---|
1096 | for (int i = legEnd[legMax]; i >= legBeg[legMax]; --i) |
---|
1097 | iPartonMax.push_back( iParton[i] ); |
---|
1098 | iPartonMax.push_back( iDiquark ); |
---|
1099 | |
---|
1100 | // Recluster gluons nearby to diquark end when taken too much energy. |
---|
1101 | int iPsize = iPartonMax.size(); |
---|
1102 | double m0Diquark = event[iDiquark].m0(); |
---|
1103 | while (iPsize > 2) { |
---|
1104 | Vec4 pGluNear = event[ iPartonMax[iPsize - 2] ].p(); |
---|
1105 | if ( pDiquark.mCalc() > 0. |
---|
1106 | && (pDiquark + 0.5 * pGluNear).mCalc() > m0Diquark + mJoin ) break; |
---|
1107 | pDiquark += pGluNear; |
---|
1108 | event[iDiquark].p( pDiquark ); |
---|
1109 | event[iDiquark].m( pDiquark.mCalc() ); |
---|
1110 | iPartonMax.pop_back(); |
---|
1111 | --iPsize; |
---|
1112 | iPartonMax[iPsize - 1] = iDiquark; |
---|
1113 | } |
---|
1114 | |
---|
1115 | // Modify parton list to remaining leg + remnant of the first two. |
---|
1116 | iParton = iPartonMax; |
---|
1117 | |
---|
1118 | // Done. |
---|
1119 | return true; |
---|
1120 | } |
---|
1121 | |
---|
1122 | //-------------------------------------------------------------------------- |
---|
1123 | |
---|
1124 | // Find the boost matrix to the rest frame of a junction, |
---|
1125 | // given the three respective endpoint four-momenta. |
---|
1126 | |
---|
1127 | RotBstMatrix StringFragmentation::junctionRestFrame(Vec4& p0, Vec4& p1, |
---|
1128 | Vec4& p2) { |
---|
1129 | |
---|
1130 | // Calculate masses and other invariants. |
---|
1131 | Vec4 pSumJun = p0 + p1 + p2; |
---|
1132 | double sHat = pSumJun.m2Calc(); |
---|
1133 | double pp[3][3]; |
---|
1134 | pp[0][0] = p0.m2Calc(); |
---|
1135 | pp[1][1] = p1.m2Calc(); |
---|
1136 | pp[2][2] = p2.m2Calc(); |
---|
1137 | pp[0][1] = pp[1][0] = p0 * p1; |
---|
1138 | pp[0][2] = pp[2][0] = p0 * p2; |
---|
1139 | pp[1][2] = pp[2][1] = p1 * p2; |
---|
1140 | |
---|
1141 | // Requirement (eiMax)_j = pi*pj/mj < (eiMax)_k = pi*pk/mk, used below, |
---|
1142 | // here rewritten as pi*pj * mk < pi*pk * mj and squared. |
---|
1143 | double eMax01 = pow2(pp[0][1]) * pp[2][2]; |
---|
1144 | double eMax02 = pow2(pp[0][2]) * pp[1][1]; |
---|
1145 | double eMax12 = pow2(pp[1][2]) * pp[0][0]; |
---|
1146 | |
---|
1147 | // Initially pick i to be the most massive parton. but allow other tries. |
---|
1148 | int i = (pp[1][1] > pp[0][0]) ? 1 : 0; |
---|
1149 | if (pp[2][2] > max(pp[0][0], pp[1][1])) i = 2; |
---|
1150 | int j, k; |
---|
1151 | double ei = 0.; |
---|
1152 | double ej = 0.; |
---|
1153 | double ek = 0.; |
---|
1154 | for (int iTry = 0; iTry < 3; ++iTry) { |
---|
1155 | |
---|
1156 | // Pick j to give minimal eiMax, and k the third vector. |
---|
1157 | if (i == 0) j = (eMax02 < eMax01) ? 2 : 1; |
---|
1158 | else if (i == 1) j = (eMax12 < eMax01) ? 2 : 0; |
---|
1159 | else j = (eMax12 < eMax02) ? 1 : 0; |
---|
1160 | k = 3 - i - j; |
---|
1161 | |
---|
1162 | // Alternative names according to i, j, k conventions. |
---|
1163 | double m2i = pp[i][i]; |
---|
1164 | double m2j = pp[j][j]; |
---|
1165 | double m2k = pp[k][k]; |
---|
1166 | double pipj = pp[i][j]; |
---|
1167 | double pipk = pp[i][k]; |
---|
1168 | double pjpk = pp[j][k]; |
---|
1169 | |
---|
1170 | // Trivial to find new parton energies if all three partons are massless. |
---|
1171 | if (m2i < M2MAXJRF) { |
---|
1172 | ei = sqrt( 2. * pipk * pipj / (3. * pjpk) ); |
---|
1173 | ej = sqrt( 2. * pjpk * pipj / (3. * pipk) ); |
---|
1174 | ek = sqrt( 2. * pipk * pjpk / (3. * pipj) ); |
---|
1175 | |
---|
1176 | // Else find three-momentum range for parton i and values at extremes. |
---|
1177 | } else { |
---|
1178 | // Minimum when i is at rest. |
---|
1179 | double piMin = 0.; |
---|
1180 | double eiMin = sqrt(m2i); |
---|
1181 | double ejMin = pipj / eiMin; |
---|
1182 | double ekMin = pipk / eiMin; |
---|
1183 | double pjMin = sqrtpos( ejMin*ejMin - m2j ); |
---|
1184 | double pkMin = sqrtpos( ekMin*ekMin - m2k ); |
---|
1185 | double fMin = ejMin * ekMin + 0.5 * pjMin * pkMin - pjpk; |
---|
1186 | // Maximum estimated when j + k is at rest, alternatively j at rest. |
---|
1187 | double eiMax = (pipj + pipk) / sqrt(m2j + m2k + 2. * pjpk); |
---|
1188 | if (m2j > M2MAXJRF) eiMax = min( eiMax, pipj / sqrt(m2j) ); |
---|
1189 | double piMax = sqrtpos( eiMax*eiMax - m2i ); |
---|
1190 | double temp = eiMax*eiMax - 0.25 *piMax*piMax; |
---|
1191 | double pjMax = (eiMax * sqrtpos( pipj*pipj - m2j * temp ) |
---|
1192 | - 0.5 * piMax * pipj) / temp; |
---|
1193 | double pkMax = (eiMax * sqrtpos( pipk*pipk - m2k * temp ) |
---|
1194 | - 0.5 * piMax * pipk) / temp; |
---|
1195 | double ejMax = sqrt(pjMax*pjMax + m2j); |
---|
1196 | double ekMax = sqrt(pkMax*pkMax + m2k); |
---|
1197 | double fMax = ejMax * ekMax + 0.5 * pjMax * pkMax - pjpk; |
---|
1198 | |
---|
1199 | // If unexpected values at upper endpoint then pick another parton. |
---|
1200 | if (fMax > 0.) { |
---|
1201 | int iPrel = (i + 1)%3; |
---|
1202 | if (pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;} |
---|
1203 | ++iTry; |
---|
1204 | iPrel = (i + 2)%3; |
---|
1205 | if (iTry < 3 && pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;} |
---|
1206 | } |
---|
1207 | |
---|
1208 | // Start binary + linear search to find solution inside range. |
---|
1209 | int iterMin = 0; |
---|
1210 | int iterMax = 0; |
---|
1211 | double pi = 0.5 * (piMin + piMax); |
---|
1212 | for (int iter = 0; iter < NTRYJRFEQ; ++iter) { |
---|
1213 | |
---|
1214 | // Derive momentum of other two partons and distance to root. |
---|
1215 | ei = sqrt(pi*pi + m2i); |
---|
1216 | temp = ei*ei - 0.25 * pi*pi; |
---|
1217 | double pj = (ei * sqrtpos( pipj*pipj - m2j * temp ) |
---|
1218 | - 0.5 * pi * pipj) / temp; |
---|
1219 | double pk = (ei * sqrtpos( pipk*pipk - m2k * temp ) |
---|
1220 | - 0.5 * pi * pipk) / temp; |
---|
1221 | ej = sqrt(pj*pj + m2j); |
---|
1222 | ek = sqrt(pk*pk + m2k); |
---|
1223 | double fNow = ej * ek + 0.5 * pj * pk - pjpk; |
---|
1224 | |
---|
1225 | // Replace lower or upper bound by new value. |
---|
1226 | if (fNow > 0.) { ++iterMin; piMin = pi; fMin = fNow;} |
---|
1227 | else {++iterMax; piMax = pi; fMax = fNow;} |
---|
1228 | |
---|
1229 | // Pick next i momentum to explore, hopefully closer to root. |
---|
1230 | if (2 * iter < NTRYJRFEQ |
---|
1231 | && (iterMin < 2 || iterMax < 2 || 4 * iter < NTRYJRFEQ)) |
---|
1232 | { pi = 0.5 * (piMin + piMax); continue;} |
---|
1233 | if (fMin < 0. || fMax > 0. || abs(fNow) < CONVJRFEQ * sHat) break; |
---|
1234 | pi = piMin + (piMax - piMin) * fMin / (fMin - fMax); |
---|
1235 | } |
---|
1236 | |
---|
1237 | // If arrived here then either succeeded or exhausted possibilities. |
---|
1238 | } break; |
---|
1239 | } |
---|
1240 | |
---|
1241 | // Now we know the energies in the junction rest frame. |
---|
1242 | double eNew[3] = { 0., 0., 0.}; |
---|
1243 | eNew[i] = ei; |
---|
1244 | eNew[j] = ej; |
---|
1245 | eNew[k] = ek; |
---|
1246 | |
---|
1247 | // Boost (copy of) partons to their rest frame. |
---|
1248 | RotBstMatrix Mmove; |
---|
1249 | Vec4 p0cm = p0; |
---|
1250 | Vec4 p1cm = p1; |
---|
1251 | Vec4 p2cm = p2; |
---|
1252 | Mmove.bstback(pSumJun); |
---|
1253 | p0cm.rotbst(Mmove); |
---|
1254 | p1cm.rotbst(Mmove); |
---|
1255 | p2cm.rotbst(Mmove); |
---|
1256 | |
---|
1257 | // Construct difference vectors and the boost to junction rest frame. |
---|
1258 | Vec4 pDir01 = p0cm / p0cm.e() - p1cm / p1cm.e(); |
---|
1259 | Vec4 pDir02 = p0cm / p0cm.e() - p2cm / p2cm.e(); |
---|
1260 | double pDiff01 = pDir01.pAbs2(); |
---|
1261 | double pDiff02 = pDir02.pAbs2(); |
---|
1262 | double pDiff0102 = dot3(pDir01, pDir02); |
---|
1263 | double eDiff01 = eNew[0] / p0cm.e() - eNew[1] / p1cm.e(); |
---|
1264 | double eDiff02 = eNew[0] / p0cm.e() - eNew[2] / p2cm.e(); |
---|
1265 | double denom = pDiff01 * pDiff02 - pDiff0102*pDiff0102; |
---|
1266 | double coef01 = (eDiff01 * pDiff02 - eDiff02 * pDiff0102) / denom; |
---|
1267 | double coef02 = (eDiff02 * pDiff01 - eDiff01 * pDiff0102) / denom; |
---|
1268 | Vec4 vJunction = coef01 * pDir01 + coef02 * pDir02; |
---|
1269 | vJunction.e( sqrt(1. + vJunction.pAbs2()) ); |
---|
1270 | |
---|
1271 | // Add two boosts, giving final result. |
---|
1272 | Mmove.bst(vJunction); |
---|
1273 | return Mmove; |
---|
1274 | |
---|
1275 | } |
---|
1276 | |
---|
1277 | //-------------------------------------------------------------------------- |
---|
1278 | |
---|
1279 | // When string fragmentation has failed several times, |
---|
1280 | // try to join some more nearby partons. |
---|
1281 | |
---|
1282 | int StringFragmentation::extraJoin(double facExtra, Event& event) { |
---|
1283 | |
---|
1284 | // Keep on looping while pairs found below joining threshold. |
---|
1285 | int nJoin = 0; |
---|
1286 | int iPsize = iParton.size(); |
---|
1287 | while (iPsize > 2) { |
---|
1288 | |
---|
1289 | // Look for the pair of neighbour partons (along string) with |
---|
1290 | // the smallest invariant mass (subtracting quark masses). |
---|
1291 | int iJoinMin = -1; |
---|
1292 | double mJoinMin = 2. * facExtra * mJoin; |
---|
1293 | for (int i = 0; i < iPsize - 1; ++i) { |
---|
1294 | Particle& parton1 = event[ iParton[i] ]; |
---|
1295 | Particle& parton2 = event[ iParton[i + 1] ]; |
---|
1296 | Vec4 pSumNow; |
---|
1297 | pSumNow += (parton2.isGluon()) ? 0.5 * parton1.p() : parton1.p(); |
---|
1298 | pSumNow += (parton2.isGluon()) ? 0.5 * parton2.p() : parton2.p(); |
---|
1299 | double mJoinNow = pSumNow.mCalc(); |
---|
1300 | if (!parton1.isGluon()) mJoinNow -= parton1.m0(); |
---|
1301 | if (!parton2.isGluon()) mJoinNow -= parton2.m0(); |
---|
1302 | if (mJoinNow < mJoinMin) { iJoinMin = i; mJoinMin = mJoinNow; } |
---|
1303 | } |
---|
1304 | |
---|
1305 | // Decide whether to join, if not finished. |
---|
1306 | if (iJoinMin < 0 || mJoinMin > facExtra * mJoin) return nJoin; |
---|
1307 | ++nJoin; |
---|
1308 | |
---|
1309 | // Create new joined parton. |
---|
1310 | int iJoin1 = iParton[iJoinMin]; |
---|
1311 | int iJoin2 = iParton[iJoinMin + 1]; |
---|
1312 | int idNew = (event[iJoin1].isGluon()) ? event[iJoin2].id() |
---|
1313 | : event[iJoin1].id(); |
---|
1314 | int colNew = event[iJoin1].col(); |
---|
1315 | int acolNew = event[iJoin2].acol(); |
---|
1316 | if (colNew == acolNew) { |
---|
1317 | colNew = event[iJoin2].col(); |
---|
1318 | acolNew = event[iJoin1].acol(); |
---|
1319 | } |
---|
1320 | Vec4 pNew = event[iJoin1].p() + event[iJoin2].p(); |
---|
1321 | |
---|
1322 | // Append joined parton to event record and reduce parton list. |
---|
1323 | int iNew = event.append( idNew, 73, min(iJoin1, iJoin2), |
---|
1324 | max(iJoin1, iJoin2), 0, 0, colNew, acolNew, pNew, pNew.mCalc() ); |
---|
1325 | iParton[iJoinMin] = iNew; |
---|
1326 | for (int i = iJoinMin + 1; i < iPsize - 1; ++i) |
---|
1327 | iParton[i] = iParton[i + 1]; |
---|
1328 | iParton.pop_back(); |
---|
1329 | --iPsize; |
---|
1330 | |
---|
1331 | // Done. |
---|
1332 | } |
---|
1333 | return nJoin; |
---|
1334 | } |
---|
1335 | |
---|
1336 | //========================================================================== |
---|
1337 | |
---|
1338 | } // end namespace Pythia8 |
---|