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