1 | // FragmentationSystems.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 |
---|
7 | // ColConfig, StringRegion and StringSystem classes. |
---|
8 | |
---|
9 | #include "FragmentationSystems.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // The ColConfig class. |
---|
16 | |
---|
17 | //-------------------------------------------------------------------------- |
---|
18 | |
---|
19 | // Constants: could be changed here if desired, but normally should not. |
---|
20 | // These are of technical nature, as described for each. |
---|
21 | |
---|
22 | // A typical u/d constituent mass. |
---|
23 | const double ColConfig::CONSTITUENTMASS = 0.325; |
---|
24 | |
---|
25 | //-------------------------------------------------------------------------- |
---|
26 | |
---|
27 | // Initialize and save pointers. |
---|
28 | |
---|
29 | void ColConfig::init(Info* infoPtrIn, Settings& settings, |
---|
30 | StringFlav* flavSelPtrIn) { |
---|
31 | |
---|
32 | // Save pointers. |
---|
33 | infoPtr = infoPtrIn; |
---|
34 | flavSelPtr = flavSelPtrIn; |
---|
35 | |
---|
36 | // Joining of nearby partons along the string. |
---|
37 | mJoin = settings.parm("FragmentationSystems:mJoin"); |
---|
38 | |
---|
39 | // For consistency ensure that mJoin is bigger than in StringRegion. |
---|
40 | mJoin = max( mJoin, 2. * StringRegion::MJOIN); |
---|
41 | |
---|
42 | // Simplification of q q q junction topology to quark - diquark one. |
---|
43 | mJoinJunction = settings.parm("FragmentationSystems:mJoinJunction"); |
---|
44 | mStringMin = settings.parm("HadronLevel:mStringMin"); |
---|
45 | |
---|
46 | } |
---|
47 | |
---|
48 | //-------------------------------------------------------------------------- |
---|
49 | |
---|
50 | // Insert a new colour singlet system in ascending mass order. |
---|
51 | // Calculate its properties. Join nearby partons. |
---|
52 | |
---|
53 | bool ColConfig::insert( vector<int>& iPartonIn, Event& event) { |
---|
54 | |
---|
55 | // Find momentum and invariant mass of system, minus endpoint masses. |
---|
56 | Vec4 pSumIn; |
---|
57 | double mSumIn = 0.; |
---|
58 | bool hasJunctionIn = false; |
---|
59 | int nJunctionLegs = 0; |
---|
60 | for (int i = 0; i < int(iPartonIn.size()); ++i) { |
---|
61 | if (iPartonIn[i] < 0) { |
---|
62 | hasJunctionIn = true; |
---|
63 | ++nJunctionLegs; |
---|
64 | } else { |
---|
65 | pSumIn += event[ iPartonIn[i] ].p(); |
---|
66 | if (!event[ iPartonIn[i] ].isGluon()) |
---|
67 | mSumIn += event[ iPartonIn[i] ].constituentMass(); |
---|
68 | } |
---|
69 | } |
---|
70 | double massIn = pSumIn.mCalc(); |
---|
71 | double massExcessIn = massIn - mSumIn; |
---|
72 | |
---|
73 | // Check for rare triple- and higher junction systems (like J-Jbar-J) |
---|
74 | if (nJunctionLegs >= 5) { |
---|
75 | infoPtr->errorMsg("Error in ColConfig::insert: " |
---|
76 | "junction topology too complicated; too many junction legs"); |
---|
77 | return false; |
---|
78 | } |
---|
79 | // Check that junction systems have at least three legs. |
---|
80 | else if (nJunctionLegs > 0 && nJunctionLegs <= 2) { |
---|
81 | infoPtr->errorMsg("Error in ColConfig::insert: " |
---|
82 | "junction topology inconsistent; too few junction legs"); |
---|
83 | return false; |
---|
84 | } |
---|
85 | |
---|
86 | // Check that momenta do not contain not-a-number. |
---|
87 | if (abs(massExcessIn) >= 0.); |
---|
88 | else { |
---|
89 | infoPtr->errorMsg("Error in ColConfig::insert: " |
---|
90 | "not-a-number system mass"); |
---|
91 | return false; |
---|
92 | } |
---|
93 | |
---|
94 | // Identify closed gluon loop. Assign "endpoint" masses as light quarks. |
---|
95 | bool isClosedIn = (iPartonIn[0] >= 0 && event[ iPartonIn[0] ].col() != 0 |
---|
96 | && event[ iPartonIn[0] ].acol() != 0 ); |
---|
97 | if (isClosedIn) massExcessIn -= 2. * CONSTITUENTMASS; |
---|
98 | |
---|
99 | // For junction topology: join two nearby legs into a diquark. |
---|
100 | if (hasJunctionIn && joinJunction( iPartonIn, event, massExcessIn)) |
---|
101 | hasJunctionIn = false; |
---|
102 | |
---|
103 | // Loop while > 2 partons left and hope of finding joining pair. |
---|
104 | bool hasJoined = true; |
---|
105 | while (hasJoined && iPartonIn.size() > 2) { |
---|
106 | |
---|
107 | // Look for the pair of neighbour partons (along string) with |
---|
108 | // the smallest invariant mass (subtracting quark masses). |
---|
109 | int iJoinMin = -1; |
---|
110 | double mJoinMin = 2. * mJoin; |
---|
111 | int nSize = iPartonIn.size(); |
---|
112 | int nPair = (isClosedIn) ? nSize : nSize - 1; |
---|
113 | for (int i = 0; i < nPair; ++i) { |
---|
114 | // Keep three legs of junction separate. |
---|
115 | if (iPartonIn[i] < 0 || iPartonIn[(i + 1)%nSize] < 0) continue; |
---|
116 | Particle& parton1 = event[ iPartonIn[i] ]; |
---|
117 | Particle& parton2 = event[ iPartonIn[(i + 1)%nSize] ]; |
---|
118 | // Avoid joining non-partons, e.g. gluino/squark for R-hadron. |
---|
119 | if (!parton1.isParton() || !parton2.isParton()) continue; |
---|
120 | Vec4 pSumNow; |
---|
121 | pSumNow += (parton1.isGluon()) ? 0.5 * parton1.p() : parton1.p(); |
---|
122 | pSumNow += (parton2.isGluon()) ? 0.5 * parton2.p() : parton2.p(); |
---|
123 | double mJoinNow = pSumNow.mCalc(); |
---|
124 | if (!parton1.isGluon()) mJoinNow -= parton1.m(); |
---|
125 | if (!parton2.isGluon()) mJoinNow -= parton2.m(); |
---|
126 | if (mJoinNow < mJoinMin) { iJoinMin = i; mJoinMin = mJoinNow; } |
---|
127 | } |
---|
128 | |
---|
129 | // If sufficiently nearby then join into one new parton. |
---|
130 | // Note: error sensitivity to mJoin indicates unstable precedure?? |
---|
131 | hasJoined = false; |
---|
132 | if (mJoinMin < mJoin) { |
---|
133 | int iJoin1 = iPartonIn[iJoinMin]; |
---|
134 | int iJoin2 = iPartonIn[(iJoinMin + 1)%nSize]; |
---|
135 | int idNew = (event[iJoin1].isGluon()) ? event[iJoin2].id() |
---|
136 | : event[iJoin1].id(); |
---|
137 | int colNew = event[iJoin1].col(); |
---|
138 | int acolNew = event[iJoin2].acol(); |
---|
139 | if (colNew == acolNew) { |
---|
140 | colNew = event[iJoin2].col(); |
---|
141 | acolNew = event[iJoin1].acol(); |
---|
142 | } |
---|
143 | Vec4 pNew = event[iJoin1].p() + event[iJoin2].p(); |
---|
144 | |
---|
145 | // Append joined parton to event record. |
---|
146 | int iNew = event.append( idNew, 73, min(iJoin1, iJoin2), |
---|
147 | max(iJoin1, iJoin2), 0, 0, colNew, acolNew, pNew, pNew.mCalc() ); |
---|
148 | |
---|
149 | // Displaced lifetime/vertex; mothers should be same but prefer quark. |
---|
150 | int iVtx = (event[iJoin1].isGluon()) ? iJoin2 : iJoin1; |
---|
151 | event[iNew].tau( event[iVtx].tau() ); |
---|
152 | if (event[iVtx].hasVertex()) event[iNew].vProd( event[iVtx].vProd() ); |
---|
153 | |
---|
154 | // Mark joined partons and reduce remaining system. |
---|
155 | event[iJoin1].statusNeg(); |
---|
156 | event[iJoin2].statusNeg(); |
---|
157 | event[iJoin1].daughter1(iNew); |
---|
158 | event[iJoin2].daughter1(iNew); |
---|
159 | if (iJoinMin == nSize - 1) iPartonIn[0] = iNew; |
---|
160 | else { |
---|
161 | iPartonIn[iJoinMin] = iNew; |
---|
162 | for (int i = iJoinMin + 1; i < nSize - 1; ++i) |
---|
163 | iPartonIn[i] = iPartonIn[i + 1]; |
---|
164 | } |
---|
165 | iPartonIn.pop_back(); |
---|
166 | |
---|
167 | // If joined,then loopback to look for more. |
---|
168 | hasJoined = true; |
---|
169 | } |
---|
170 | } |
---|
171 | |
---|
172 | // Store new colour singlet system at the end. |
---|
173 | singlets.push_back( ColSinglet(iPartonIn, pSumIn, massIn, |
---|
174 | massExcessIn, hasJunctionIn, isClosedIn) ); |
---|
175 | |
---|
176 | // Now move around, so that smallest mass excesses come first. |
---|
177 | int iInsert = singlets.size() - 1; |
---|
178 | for (int iSub = singlets.size() - 2; iSub >= 0; --iSub) { |
---|
179 | if (massExcessIn > singlets[iSub].massExcess) break; |
---|
180 | singlets[iSub + 1] = singlets[iSub]; |
---|
181 | iInsert = iSub; |
---|
182 | } |
---|
183 | if (iInsert < int(singlets.size()) - 1) singlets[iInsert] = |
---|
184 | ColSinglet(iPartonIn, pSumIn, massIn, massExcessIn, |
---|
185 | hasJunctionIn, isClosedIn); |
---|
186 | |
---|
187 | // Done. |
---|
188 | return true; |
---|
189 | } |
---|
190 | |
---|
191 | //-------------------------------------------------------------------------- |
---|
192 | |
---|
193 | // Join two legs of junction to a diquark for small invariant masses. |
---|
194 | // Note: for junction system, iPartonIn points to structure |
---|
195 | // (-code0) g...g.q0 (-code1) g...g.q1 (-code2) g...g.q2 |
---|
196 | |
---|
197 | bool ColConfig::joinJunction( vector<int>& iPartonIn, Event& event, |
---|
198 | double massExcessIn) { |
---|
199 | |
---|
200 | // Find four-momentum and endpoint quarks and masses on the three legs. |
---|
201 | Vec4 pLeg[3]; |
---|
202 | double mLeg[3] = { 0., 0., 0.}; |
---|
203 | int idAbsLeg[3]; |
---|
204 | int leg = -1; |
---|
205 | for (int i = 0; i < int(iPartonIn.size()); ++ i) { |
---|
206 | if (iPartonIn[i] < 0) ++leg; |
---|
207 | else { |
---|
208 | pLeg[leg] += event[ iPartonIn[i] ].p(); |
---|
209 | mLeg[leg] = event[ iPartonIn[i] ].m(); |
---|
210 | idAbsLeg[leg] = event[ iPartonIn[i] ].idAbs(); |
---|
211 | } |
---|
212 | } |
---|
213 | |
---|
214 | // Calculate invariant mass of three pairs, minus endpoint masses. |
---|
215 | double m01 = (pLeg[0] + pLeg[1]).mCalc() - mLeg[0] - mLeg[1]; |
---|
216 | double m02 = (pLeg[0] + pLeg[2]).mCalc() - mLeg[0] - mLeg[2]; |
---|
217 | double m12 = (pLeg[1] + pLeg[2]).mCalc() - mLeg[1] - mLeg[2]; |
---|
218 | |
---|
219 | // Find lowest-mass pair not involving diquark. |
---|
220 | double mMin = mJoinJunction + 1.; |
---|
221 | int legA = -1; |
---|
222 | int legB = -1; |
---|
223 | if (m01 < mMin && idAbsLeg[0] < 9 && idAbsLeg[1] < 9) { |
---|
224 | mMin = m01; |
---|
225 | legA = 0; |
---|
226 | legB = 1; |
---|
227 | } |
---|
228 | if (m02 < mMin && idAbsLeg[0] < 9 && idAbsLeg[2] < 9) { |
---|
229 | mMin = m02; |
---|
230 | legA = 0; |
---|
231 | legB = 2; |
---|
232 | } |
---|
233 | if (m12 < mMin && idAbsLeg[1] < 9 && idAbsLeg[2] < 9) { |
---|
234 | mMin = m12; |
---|
235 | legA = 1; |
---|
236 | legB = 2; |
---|
237 | } |
---|
238 | int legC = 3 - legA - legB; |
---|
239 | |
---|
240 | // Nothing to do if no two legs have small invariant mass, and |
---|
241 | // system as a whole is above MiniStringFragmentation threshold. |
---|
242 | if (legA == -1 || (mMin > mJoinJunction && massExcessIn > mStringMin)) |
---|
243 | return false; |
---|
244 | |
---|
245 | // Construct separate index arrays for the three legs. |
---|
246 | vector<int> iLegA, iLegB, iLegC; |
---|
247 | leg = -1; |
---|
248 | for (int i = 0; i < int(iPartonIn.size()); ++ i) { |
---|
249 | if (iPartonIn[i] < 0) ++leg; |
---|
250 | else if( leg == legA) iLegA.push_back( iPartonIn[i] ); |
---|
251 | else if( leg == legB) iLegB.push_back( iPartonIn[i] ); |
---|
252 | else if( leg == legC) iLegC.push_back( iPartonIn[i] ); |
---|
253 | } |
---|
254 | |
---|
255 | // First step: successively combine any gluons on the two legs. |
---|
256 | // (Presumably overkill; not likely to be (m)any extra gluons.) |
---|
257 | // (Do as successive binary joinings, so only need two mothers.) |
---|
258 | for (leg = 0; leg < 2; ++leg) { |
---|
259 | vector<int>& iLegNow = (leg == 0) ? iLegA : iLegB; |
---|
260 | int sizeNow = iLegNow.size(); |
---|
261 | for (int i = sizeNow - 2; i >= 0; --i) { |
---|
262 | int iQ = iLegNow.back(); |
---|
263 | int iG = iLegNow[i]; |
---|
264 | int colNew = (event[iQ].id() > 0) ? event[iG].col() : 0; |
---|
265 | int acolNew = (event[iQ].id() < 0) ? event[iG].acol() : 0; |
---|
266 | Vec4 pNew = event[iQ].p() + event[iG].p(); |
---|
267 | int iNew = event.append( event[iQ].id(), 74, min(iQ, iG), |
---|
268 | max(iQ, iG), 0, 0, colNew, acolNew, pNew, pNew.mCalc() ); |
---|
269 | |
---|
270 | // Mark joined partons and update iLeg end. |
---|
271 | event[iQ].statusNeg(); |
---|
272 | event[iG].statusNeg(); |
---|
273 | event[iQ].daughter1(iNew); |
---|
274 | event[iG].daughter1(iNew); |
---|
275 | iLegNow.back() = iNew; |
---|
276 | } |
---|
277 | } |
---|
278 | |
---|
279 | // Second step: combine two quarks into a diquark. |
---|
280 | int iQA = iLegA.back(); |
---|
281 | int iQB = iLegB.back(); |
---|
282 | int idQA = event[iQA].id(); |
---|
283 | int idQB = event[iQB].id(); |
---|
284 | int idNew = flavSelPtr->makeDiquark( idQA, idQB ); |
---|
285 | // Diquark colour is opposite to parton closest to junction on third leg. |
---|
286 | int colNew = (idNew > 0) ? 0 : event[ iLegC[0] ].acol(); |
---|
287 | int acolNew = (idNew > 0) ? event[ iLegC[0] ].col() : 0; |
---|
288 | Vec4 pNew = pLeg[legA] + pLeg[legB]; |
---|
289 | int iNew = event.append( idNew, 74, min(iQA, iQB), max( iQA, iQB), |
---|
290 | 0, 0, colNew, acolNew, pNew, pNew.mCalc() ); |
---|
291 | |
---|
292 | // Mark joined partons and reduce remaining system. |
---|
293 | event[iQA].statusNeg(); |
---|
294 | event[iQB].statusNeg(); |
---|
295 | event[iQA].daughter1(iNew); |
---|
296 | event[iQB].daughter1(iNew); |
---|
297 | iPartonIn.resize(0); |
---|
298 | iPartonIn.push_back( iNew); |
---|
299 | for (int i = 0; i < int(iLegC.size()) ; ++i) |
---|
300 | iPartonIn.push_back( iLegC[i]); |
---|
301 | |
---|
302 | // Remove junction from event record list, identifying by colour. |
---|
303 | int iJun = -1; |
---|
304 | for (int i = 0; i < event.sizeJunction(); ++i) |
---|
305 | for (int j = 0; j < 3; ++ j) |
---|
306 | if ( event.colJunction(i,j) == max(colNew, acolNew) ) iJun = i; |
---|
307 | if (iJun >= 0) event.eraseJunction(iJun); |
---|
308 | |
---|
309 | // Done, having eliminated junction. |
---|
310 | return true; |
---|
311 | |
---|
312 | } |
---|
313 | |
---|
314 | //-------------------------------------------------------------------------- |
---|
315 | |
---|
316 | // Collect all partons of singlet to be consecutively ordered. |
---|
317 | |
---|
318 | void ColConfig::collect(int iSub, Event& event, bool skipTrivial) { |
---|
319 | |
---|
320 | // Check that all partons have positive energy. |
---|
321 | for (int i = 0; i < singlets[iSub].size(); ++i) { |
---|
322 | int iNow = singlets[iSub].iParton[i]; |
---|
323 | if (iNow > 0 && event[iNow].e() < 0.) |
---|
324 | infoPtr->errorMsg("Warning in ColConfig::collect: " |
---|
325 | "negative-energy parton encountered"); |
---|
326 | } |
---|
327 | |
---|
328 | // Partons may already have been collected, e.g. at ministring collapse. |
---|
329 | if (singlets[iSub].isCollected) return; |
---|
330 | singlets[iSub].isCollected = true; |
---|
331 | |
---|
332 | // Check if partons already "by chance" happen to be ordered. |
---|
333 | bool inOrder = true; |
---|
334 | for (int i = 0; i < singlets[iSub].size() - 1; ++i) { |
---|
335 | int iFirst = singlets[iSub].iParton[i]; |
---|
336 | if (iFirst < 0) continue; |
---|
337 | int iSecond = singlets[iSub].iParton[i + 1]; |
---|
338 | if (iSecond < 0) iSecond = singlets[iSub].iParton[i + 2]; |
---|
339 | if (iSecond != iFirst + 1) { inOrder = false; break;} |
---|
340 | } |
---|
341 | |
---|
342 | // Normally done if in order, but sometimes may need to copy anyway. |
---|
343 | if (inOrder && skipTrivial) return; |
---|
344 | |
---|
345 | // Copy down system. Update current partons. |
---|
346 | for (int i = 0; i < singlets[iSub].size(); ++i) { |
---|
347 | int iOld = singlets[iSub].iParton[i]; |
---|
348 | if (iOld < 0) continue; |
---|
349 | int iNew = event.copy(iOld, 71); |
---|
350 | singlets[iSub].iParton[i] = iNew; |
---|
351 | } |
---|
352 | |
---|
353 | // Done. |
---|
354 | } |
---|
355 | |
---|
356 | //-------------------------------------------------------------------------- |
---|
357 | |
---|
358 | // Find to which singlet system a particle belongs. |
---|
359 | |
---|
360 | int ColConfig::findSinglet(int i) { |
---|
361 | |
---|
362 | // Loop through all systems and all members in them. |
---|
363 | for (int iSub = 0; iSub < int(singlets.size()); ++iSub) |
---|
364 | for (int iMem = 0; iMem < singlets[iSub].size(); ++iMem) |
---|
365 | if (singlets[iSub].iParton[iMem] == i) return iSub; |
---|
366 | |
---|
367 | // Done without having found particle; return -1 = error code. |
---|
368 | return -1; |
---|
369 | } |
---|
370 | |
---|
371 | //-------------------------------------------------------------------------- |
---|
372 | |
---|
373 | // List all currently identified singlets. |
---|
374 | |
---|
375 | void ColConfig::list(ostream& os) const { |
---|
376 | |
---|
377 | // Header. Loop over all individual singlets. |
---|
378 | os << "\n -------- Colour Singlet Systems Listing -------------------\n"; |
---|
379 | for (int iSub = 0; iSub < int(singlets.size()); ++iSub) { |
---|
380 | |
---|
381 | // List all partons belonging to each singlet. |
---|
382 | os << " singlet " << iSub << " contains " ; |
---|
383 | for (int i = 0; i < singlets[iSub].size(); ++i) |
---|
384 | os << singlets[iSub].iParton[i] << " "; |
---|
385 | os << "\n"; |
---|
386 | |
---|
387 | // Done. |
---|
388 | } |
---|
389 | } |
---|
390 | |
---|
391 | //========================================================================== |
---|
392 | |
---|
393 | // The StringRegion class. |
---|
394 | |
---|
395 | // Currently a number of simplifications, in particular ?? |
---|
396 | // 1) No popcorn baryon production. |
---|
397 | // 2) Simplified treatment of pT in stepping and joining. |
---|
398 | |
---|
399 | //-------------------------------------------------------------------------- |
---|
400 | |
---|
401 | // Constants: could be changed here if desired, but normally should not. |
---|
402 | // These are of technical nature, as described for each. |
---|
403 | |
---|
404 | // If a string region is smaller thsan this it is assumed empty. |
---|
405 | const double StringRegion::MJOIN = 0.1; |
---|
406 | |
---|
407 | // Avoid division by zero. |
---|
408 | const double StringRegion::TINY = 1e-20; |
---|
409 | |
---|
410 | //-------------------------------------------------------------------------- |
---|
411 | |
---|
412 | // Set up four-vectors for longitudinal and transverse directions. |
---|
413 | |
---|
414 | void StringRegion::setUp(Vec4 p1, Vec4 p2, bool isMassless) { |
---|
415 | |
---|
416 | // Simple case: the two incoming four-vectors guaranteed massless. |
---|
417 | if (isMassless) { |
---|
418 | |
---|
419 | // Calculate w2, minimum value. Lightcone directions = input. |
---|
420 | w2 = 2. * (p1 * p2); |
---|
421 | if (w2 < MJOIN*MJOIN) {isSetUp = true; isEmpty = true; return;} |
---|
422 | pPos = p1; |
---|
423 | pNeg = p2; |
---|
424 | |
---|
425 | // Else allow possibility of masses for incoming partons (also gluons!). |
---|
426 | } else { |
---|
427 | |
---|
428 | // Generic four-momentum combinations. |
---|
429 | double m1Sq = p1 * p1; |
---|
430 | double m2Sq = p2 * p2; |
---|
431 | double p1p2 = p1 * p2; |
---|
432 | w2 = m1Sq + 2. * p1p2 + m2Sq; |
---|
433 | double rootSq = pow2(p1p2) - m1Sq * m2Sq; |
---|
434 | |
---|
435 | // If crazy kinematics (should not happen!) modify energies. |
---|
436 | if (w2 <= 0. || rootSq <= 0.) { |
---|
437 | if (m1Sq < 0.) m1Sq = 0.; |
---|
438 | p1.e( sqrt(m1Sq + p1.pAbs2()) ); |
---|
439 | if (m2Sq < 0.) m2Sq = 0.; |
---|
440 | p2.e( sqrt(m2Sq + p2.pAbs2()) ); |
---|
441 | p1p2 = p1 * p2; |
---|
442 | w2 = m1Sq + 2. * p1p2 + m2Sq; |
---|
443 | rootSq = pow2(p1p2) - m1Sq * m2Sq; |
---|
444 | } |
---|
445 | |
---|
446 | // If still small invariant mass then empty region (e.g. in gg system). |
---|
447 | if (w2 < MJOIN*MJOIN) {isSetUp = true; isEmpty = true; return;} |
---|
448 | |
---|
449 | // Find two lightconelike longitudinal four-vector directions. |
---|
450 | double root = sqrt( max(TINY, rootSq) ); |
---|
451 | double k1 = 0.5 * ( (m2Sq + p1p2) / root - 1.); |
---|
452 | double k2 = 0.5 * ( (m1Sq + p1p2) / root - 1.); |
---|
453 | pPos = (1. + k1) * p1 - k2 * p2; |
---|
454 | pNeg = (1. + k2) * p2 - k1 * p1; |
---|
455 | } |
---|
456 | |
---|
457 | // Find two spacelike transverse four-vector directions. |
---|
458 | // Begin by picking two sensible trial directions. |
---|
459 | Vec4 eDiff = pPos / pPos.e() - pNeg / pNeg.e(); |
---|
460 | double eDx = pow2( eDiff.px() ); |
---|
461 | double eDy = pow2( eDiff.py() ); |
---|
462 | double eDz = pow2( eDiff.pz() ); |
---|
463 | if (eDx < min(eDy, eDz)) { |
---|
464 | eX = Vec4( 1., 0., 0., 0.); |
---|
465 | eY = (eDy < eDz) ? Vec4( 0., 1., 0., 0.) : Vec4( 0., 0., 1., 0.); |
---|
466 | } else if (eDy < eDz) { |
---|
467 | eX = Vec4( 0., 1., 0., 0.); |
---|
468 | eY = (eDx < eDz) ? Vec4( 1., 0., 0., 0.) : Vec4( 0., 0., 1., 0.); |
---|
469 | } else { |
---|
470 | eX = Vec4( 0., 0., 1., 0.); |
---|
471 | eY = (eDx < eDy) ? Vec4( 1., 0., 0., 0.) : Vec4( 0., 1., 0., 0.); |
---|
472 | } |
---|
473 | |
---|
474 | // Then construct orthogonal linear combinations. |
---|
475 | double pPosNeg = pPos * pNeg; |
---|
476 | double kXPos = eX * pPos / pPosNeg; |
---|
477 | double kXNeg = eX * pNeg / pPosNeg; |
---|
478 | double kXX = 1. / sqrt( 1. + 2. * kXPos * kXNeg * pPosNeg ); |
---|
479 | double kYPos = eY * pPos / pPosNeg; |
---|
480 | double kYNeg = eY * pNeg / pPosNeg; |
---|
481 | double kYX = kXX * (kXPos * kYNeg + kXNeg * kYPos) * pPosNeg; |
---|
482 | double kYY = 1. / sqrt(1. + 2. * kYPos * kYNeg * pPosNeg - pow2(kYX)); |
---|
483 | eX = kXX * (eX - kXNeg * pPos - kXPos * pNeg); |
---|
484 | eY = kYY * (eY - kYNeg * pPos - kYPos * pNeg - kYX * eX); |
---|
485 | |
---|
486 | // Done. |
---|
487 | isSetUp = true; |
---|
488 | isEmpty = false; |
---|
489 | |
---|
490 | } |
---|
491 | |
---|
492 | //-------------------------------------------------------------------------- |
---|
493 | |
---|
494 | // Project a four-momentum onto (x+, x-, px, py). |
---|
495 | |
---|
496 | void StringRegion::project(Vec4 pIn) { |
---|
497 | |
---|
498 | // Perform projections by four-vector multiplication. |
---|
499 | xPosProj = 2. * (pIn * pNeg) / w2; |
---|
500 | xNegProj = 2. * (pIn * pPos) / w2; |
---|
501 | pxProj = - (pIn * eX); |
---|
502 | pyProj = - (pIn * eY); |
---|
503 | |
---|
504 | } |
---|
505 | |
---|
506 | //========================================================================== |
---|
507 | |
---|
508 | // The StringSystem class. |
---|
509 | |
---|
510 | //-------------------------------------------------------------------------- |
---|
511 | |
---|
512 | // Set up system from parton list. |
---|
513 | |
---|
514 | void StringSystem::setUp(vector<int>& iSys, Event& event) { |
---|
515 | |
---|
516 | // Figure out how big the system is. (Closed gluon loops?) |
---|
517 | sizePartons = iSys.size(); |
---|
518 | sizeStrings = sizePartons - 1; |
---|
519 | sizeRegions = (sizeStrings * (sizeStrings + 1)) / 2; |
---|
520 | indxReg = 2 * sizeStrings + 1; |
---|
521 | iMax = sizeStrings - 1; |
---|
522 | |
---|
523 | // Reserve space for the required number of regions. |
---|
524 | system.clear(); |
---|
525 | system.resize(sizeRegions); |
---|
526 | |
---|
527 | // Set up the lowest-lying regions. |
---|
528 | for (int i = 0; i < sizeStrings; ++i) { |
---|
529 | Vec4 p1 = event[ iSys[i] ].p(); |
---|
530 | if ( event[ iSys[i] ].isGluon() ) p1 *= 0.5; |
---|
531 | Vec4 p2 = event[ iSys[i+1] ].p(); |
---|
532 | if ( event[ iSys[i+1] ].isGluon() ) p2 *= 0.5; |
---|
533 | system[ iReg(i, iMax - i) ].setUp( p1, p2, false); |
---|
534 | } |
---|
535 | |
---|
536 | } |
---|
537 | |
---|
538 | //========================================================================== |
---|
539 | |
---|
540 | } // end namespace Pythia8 |
---|