Line data Source code
1 : // JunctionSplitting.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 : // JunctionSplitting class.
8 :
9 : // Setup the list of colours, this is needed later for finding colour chains.
10 :
11 : #include "Pythia8/JunctionSplitting.h"
12 :
13 : namespace Pythia8 {
14 :
15 : //==========================================================================
16 :
17 : // The JunctionSplitting class.
18 :
19 : //--------------------------------------------------------------------------
20 :
21 : // Constants: could be changed here if desired, but normally should not.
22 : // These are of technical nature, as described for each.
23 :
24 : // For breaking J-J string, pick a Gamma by taking a step with fictitious mass.
25 : const double JunctionSplitting::JJSTRINGM2MAX = 25.;
26 : const double JunctionSplitting::JJSTRINGM2FRAC = 0.1;
27 :
28 : // Iterate junction rest frame boost until convergence or too many tries.
29 : const double JunctionSplitting::CONVJNREST = 1e-5;
30 : const int JunctionSplitting::NTRYJNREST = 20;
31 :
32 : // Typical average transvere primary hadron mass <mThad>.
33 : const double JunctionSplitting::MTHAD = 0.9;
34 :
35 : // Minimum delta R between two partons. This is to avoid problems
36 : // with infinities.
37 : const double JunctionSplitting::MINDELTAR = 1e-7;
38 :
39 : //--------------------------------------------------------------------------
40 :
41 : // Initialize the class and all the created classes.
42 :
43 : void JunctionSplitting::init( Info* infoPtrIn, Settings& settings,
44 : Rndm* rndmPtrIn, ParticleData* particleDataPtrIn) {
45 :
46 0 : infoPtr = infoPtrIn;
47 0 : rndmPtr = rndmPtrIn;
48 :
49 : // Initialize
50 0 : colTrace.init(infoPtrIn);
51 0 : stringLength.init(infoPtrIn, settings);
52 :
53 : // Initialize auxiliary fragmentation classes.
54 0 : flavSel.init(settings, rndmPtr);
55 0 : pTSel.init(settings, *particleDataPtrIn, rndmPtr);
56 0 : zSel.init(settings, *particleDataPtrIn, rndmPtr);
57 :
58 : // Initialize string and ministring fragmentation.
59 0 : stringFrag.init(infoPtr, settings, particleDataPtrIn, rndmPtr,
60 : &flavSel, &pTSel, &zSel);
61 :
62 : // For junction processing.
63 0 : eNormJunction = settings.parm("StringFragmentation:eNormJunction");
64 0 : allowDoubleJunRem = settings.flag("ColourReconnection:allowDoubleJunRem");
65 0 : }
66 :
67 : //--------------------------------------------------------------------------
68 :
69 : // Check that all colours are connected in physical way. Also split
70 : // junction pairs, such that the hadronization can handle the configuration.
71 :
72 : bool JunctionSplitting::checkColours( Event& event) {
73 :
74 : // Not really a colour check, but a check all numbers are valid.
75 0 : for (int i = 0; i < event.size(); ++i)
76 0 : if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0.
77 0 : && abs(event[i].pz()) >= 0. && abs(event[i].e()) >= 0.
78 0 : && abs(event[i].m()) >= 0.);
79 : else {
80 0 : infoPtr->errorMsg("Warning in JunctionSplitting::CheckColours: "
81 : "not-a-number energy/momentum/mass");
82 0 : return false;
83 : }
84 :
85 : // Check if any singlet gluons were made, and if so return false.
86 0 : for (int i = 0; i < event.size(); ++i) {
87 0 : if (event[i].isFinal() && event[i].col() != 0 &&
88 0 : event[i].col() == event[i].acol()) {
89 0 : infoPtr->errorMsg("Warning in JunctionSplitting::CheckColours:"
90 : "Made a gluon colour singlet. Redoing colour configuration");
91 0 : return false;
92 : }
93 : }
94 :
95 : // Need to try and split junction structures.
96 0 : colTrace.setupColList(event);
97 0 : vector<int> iParton;
98 0 : vector<vector <int > > iPartonJun, iPartonAntiJun;
99 0 : getPartonLists(event, iPartonJun, iPartonAntiJun);
100 :
101 : // Try to split up the junction chains by splitting gluons
102 0 : if (!splitJunGluons(event, iPartonJun, iPartonAntiJun) ) {
103 0 : infoPtr->errorMsg("Warning in JunctionSplitting::CheckColours:"
104 : "Not possible to split junctions. Making new colour configuration");
105 0 : return false;
106 : }
107 :
108 : // Remove junctions if more than 2 are connected.
109 0 : if (!splitJunChains(event) ) {
110 0 : infoPtr->errorMsg("Warning in JunctionSplitting::CheckColours:"
111 : "Not possible to split junctions. Making new colour configuration");
112 0 : return false;
113 : }
114 :
115 : // Split up junction pairs.
116 0 : getPartonLists(event, iPartonJun, iPartonAntiJun);
117 0 : if (!splitJunPairs(event, iPartonJun, iPartonAntiJun) ) {
118 0 : infoPtr->errorMsg("Warning in JunctionSplitting::CheckColours:"
119 : "Not possible to split junctions. Making new colour configuration");
120 0 : return false;
121 : }
122 :
123 : // Done checking.
124 0 : return true;
125 0 : }
126 :
127 : //--------------------------------------------------------------------------
128 :
129 : // Split connected junction chains into separated, mainly by splitting gluons
130 : // into q-qbar pairs. If the junctions are directly connected
131 : // other methods are applied.
132 :
133 : bool JunctionSplitting::splitJunGluons(Event& event,
134 : vector<vector< int > >& iPartonJun, vector<vector< int > >& iPartonAntiJun) {
135 :
136 : // Loop over all junctions and all junction legs.
137 0 : for (int iJun = 0; iJun < int(iPartonJun.size()); ++iJun) {
138 :
139 : // Fill in vector of the legs content.
140 0 : vector<vector <int> > iJunLegs;
141 0 : iJunLegs.resize(3);
142 : int leg = -1;
143 0 : for (int i = 0; i < int(iPartonJun[iJun].size()); ++i) {
144 0 : if ( iPartonJun[iJun][i]/10 == iPartonJun[iJun][0]/10) ++leg;
145 0 : iJunLegs[leg].push_back(iPartonJun[iJun][i]);
146 : }
147 :
148 : // Loop over legs.
149 0 : for (int i = 0;i < int(iJunLegs.size()); ++i) {
150 :
151 : // If it is not connected to another junction, no need to do anything.
152 0 : if (iJunLegs[i].back() > 0)
153 : continue;
154 0 : int identJun = iJunLegs[i][0];
155 : // If no gluons in between two junctions, not possible to do anything.
156 0 : if (iJunLegs[i].size() == 2)
157 0 : continue;
158 :
159 : int identAntiJun = 0, iAntiLeg = -1;
160 :
161 : // Pick a new quark at random; for simplicity no diquarks.
162 : int colQ = 0, acolQ = 0;
163 0 : int idQ = flavSel.pickLightQ();
164 :
165 : // If a single gluon in between the two junctions, change it to a
166 : // quark-anti quark system.
167 0 : if ( iJunLegs[i].size() == 3) {
168 :
169 : // Store the new q qbar pair, sharing gluon colour and momentum.
170 0 : colQ = event[ iJunLegs[i][1] ].col();
171 0 : acolQ = event[ iJunLegs[i][1] ].acol();
172 0 : Vec4 pQ = 0.5 * event[ iJunLegs[i][1] ].p();
173 0 : double mQ = 0.5 * event[ iJunLegs[i][1] ].m();
174 0 : int iQ = event.append( idQ, 75, iJunLegs[i][1], 0, 0, 0, colQ, 0,
175 0 : pQ, mQ );
176 0 : int iQbar = event.append( -idQ, 75, iJunLegs[i][1], 0, 0, 0, 0, acolQ,
177 0 : pQ, mQ );
178 :
179 : // Mark split gluon.
180 0 : event[ iJunLegs[i][1] ].statusNeg();
181 0 : event[ iJunLegs[i][1] ].daughters( iQ, iQbar);
182 :
183 : // Update junction and anti junction list.
184 0 : identAntiJun = iJunLegs[i].back();
185 0 : int iOld = iJunLegs[i][1];
186 : bool erasing = false;
187 0 : for (int j = 0; j < int(iPartonJun[iJun].size()); ++j) {
188 0 : if (iPartonJun[iJun][j] == iOld)
189 0 : erasing = true;
190 0 : if (iPartonJun[iJun][j] == identAntiJun) {
191 0 : iPartonJun[iJun][j] = iQ;
192 0 : break;
193 : }
194 0 : if (erasing) {
195 0 : iPartonJun[iJun].erase(iPartonJun[iJun].begin() + j);
196 0 : --j;
197 0 : }
198 : }
199 :
200 : // Find the connected anti junction from the list of anti junctions.
201 : int iAntiJun = -1;
202 0 : for (int j = 0; j < int(iPartonAntiJun.size()); j++)
203 0 : if ( iPartonAntiJun[j][0]/10 == identAntiJun/10) {
204 : iAntiJun = j;
205 0 : break;
206 : }
207 : // If no anti junction found, something went wrong earlier.
208 0 : if (iAntiJun == -1) {
209 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunChain:"
210 : "Something went wrong in finding anti junction");
211 0 : return false;
212 : }
213 :
214 : // Update the anti junction list.
215 : erasing = false;
216 0 : for (int j = 0; j < int(iPartonAntiJun[iAntiJun].size()); ++j) {
217 0 : if ( iPartonAntiJun[iAntiJun][j] / 10 == identAntiJun / 10)
218 0 : iAntiLeg++;
219 0 : if ( iPartonAntiJun[iAntiJun][j] == identJun) {
220 0 : iPartonAntiJun[iAntiJun][j] = iQbar;
221 0 : break;
222 : }
223 : }
224 0 : }
225 : // If more than a single gluon, decide depending on mass.
226 0 : else if (iJunLegs[i].size() > 3) {
227 :
228 : // Evaluate mass-squared for all adjacent gluon pairs.
229 0 : vector<double> m2Pair;
230 : double m2Sum = 0.;
231 0 : for (int j = 1; j < int(iJunLegs[i].size()) - 2; ++j) {
232 0 : double m2Now = 0.5 * event[ iJunLegs[i][j] ].p()
233 0 : * event[ iJunLegs[i][j + 1] ].p();
234 0 : m2Pair.push_back(m2Now);
235 0 : m2Sum += m2Now;
236 0 : }
237 :
238 : // Pick breakup region with probability proportional to mass-squared.
239 0 : double m2Reg = m2Sum * rndmPtr->flat();
240 : int iReg = -1;
241 0 : do m2Reg -= m2Pair[++iReg];
242 0 : while (m2Reg > 0. && iReg < int(m2Pair.size()) - 1);
243 0 : m2Reg = m2Pair[iReg];
244 :
245 : // increase iReg with one, since it should not point towards itself.
246 0 : iReg++;
247 :
248 : // Pick breaking point of string in chosen region (symmetrically).
249 0 : double m2Temp = min( JJSTRINGM2MAX, JJSTRINGM2FRAC * m2Reg);
250 0 : double xPos = 0.5;
251 0 : double xNeg = 0.5;
252 0 : do {
253 0 : double zTemp = zSel.zFrag( idQ, 0, m2Temp);
254 0 : xPos = 1. - zTemp;
255 0 : xNeg = m2Temp / (zTemp * m2Reg);
256 0 : } while (xNeg > 1.);
257 0 : if (rndmPtr->flat() > 0.5) swap(xPos, xNeg);
258 :
259 : // Pick up two "mother" gluons of breakup. Mark them decayed.
260 0 : Particle& gJun = event[ iJunLegs[i][iReg] ];
261 0 : Particle& gAnti = event[ iJunLegs[i][iReg + 1] ];
262 0 : gJun.statusNeg();
263 0 : gAnti.statusNeg();
264 0 : int dau1 = event.size();
265 0 : gJun.daughters(dau1, dau1 + 3);
266 0 : gAnti.daughters(dau1, dau1 + 3);
267 0 : int mother1 = min( iJunLegs[i][iReg], iJunLegs[i][iReg + 1]);
268 0 : int mother2 = max( iJunLegs[i][iReg], iJunLegs[i][iReg + 1]);
269 :
270 : // Need to store variables, since it is not safe to use references
271 : // with append.
272 0 : int gJunCol = gJun.col();
273 0 : int gJunAcol = gJun.acol();
274 0 : int gAntiAcol = gAnti.acol();
275 0 : Vec4 gJunP = gJun.p();
276 0 : Vec4 gAntiP = gAnti.p();
277 0 : double gJunM = gJun.m();
278 0 : double gAntiM = gAnti.m();
279 :
280 : // Can keep one of old colours but need one new so unambiguous.
281 : colQ = gJunAcol;
282 0 : acolQ = event.nextColTag();
283 :
284 : // Store copied gluons with reduced momenta.
285 0 : int iGjun = event.append( 21, 75, mother1, mother2, 0, 0,
286 0 : gJunCol, gJunAcol, (1. - 0.5 * xPos) * gJunP,
287 0 : (1. - 0.5 * xPos) * gJunM);
288 0 : event.append( 21, 75, mother1, mother2, 0, 0,
289 0 : acolQ, gAntiAcol, (1. - 0.5 * xNeg) * gAntiP,
290 0 : (1. - 0.5 * xNeg) * gAntiM);
291 :
292 : // Store the new q qbar pair with remaining momenta.
293 0 : int iQ = event.append( idQ, 75, mother1, mother2, 0, 0,
294 0 : colQ, 0, 0.5 * xNeg * gAntiP, 0.5 * xNeg * gAntiM );
295 0 : int iQbar = event.append( -idQ, 75, mother1, mother2, 0, 0,
296 0 : 0, acolQ, 0.5 * xPos * gJunP, 0.5 * xPos * gJunM );
297 :
298 : // Update the list of junctions to reflect the splitting.
299 0 : identAntiJun = iJunLegs[i].back();
300 : bool erasing = false;
301 0 : for (int j = 0; j < int(iPartonJun[iJun].size()); ++j) {
302 0 : if (iPartonJun[iJun][j] == mother1 ||
303 0 : iPartonJun[iJun][j] == mother2)
304 0 : erasing = true;
305 :
306 0 : if ( iPartonJun[iJun][j] == identAntiJun) {
307 0 : iPartonJun[iJun][j] = iQ;
308 0 : iPartonJun[iJun].insert(iPartonJun[iJun].begin() + j, iGjun);
309 0 : break;
310 : }
311 0 : if (erasing) {
312 0 : iPartonJun[iJun].erase(iPartonJun[iJun].begin() + j);
313 0 : j--;
314 0 : }
315 : }
316 :
317 : // Find the connected anti junction from the list of anti junctions.
318 : int iAntiJun = -1;
319 0 : for (int j = 0; j < int(iPartonAntiJun.size());j++)
320 0 : if ( iPartonAntiJun[j][0]/10 == identAntiJun/10) {
321 : iAntiJun = j;
322 0 : break;
323 : }
324 : // If no anti junction found, something went wrong earlier.
325 0 : if (iAntiJun == -1) {
326 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunChain:"
327 : "Something went wrong in finding anti junction");
328 0 : return false;
329 : }
330 :
331 : // Update the anti junction list to reflect the splitting
332 0 : for (int j = 0; j < int(iPartonAntiJun[iAntiJun].size()); ++j) {
333 0 : if ( iPartonAntiJun[iAntiJun][j] / 10 == identAntiJun / 10)
334 0 : iAntiLeg++;
335 0 : if (iPartonAntiJun[iAntiJun][j] == identJun) {
336 0 : iPartonAntiJun[iAntiJun][j] = iQbar;
337 0 : break;
338 : }
339 : }
340 0 : }
341 :
342 : // Update end colours for both g -> q qbar and g g -> g g q qbar.
343 0 : event.endColJunction((-identJun)/10 - 1, i, colQ);
344 0 : event.endColJunction((-identAntiJun)/10 - 1, iAntiLeg, acolQ);
345 0 : }
346 0 : }
347 :
348 : // Done.
349 0 : return true;
350 0 : }
351 :
352 : //--------------------------------------------------------------------------
353 :
354 : // Fix chains that contain more than two junctions.
355 : // This is done by removing the minimum needed amount of junctions.
356 : // Might need to make choice based on String length, now randomly chosen.
357 :
358 : bool JunctionSplitting::splitJunChains(Event& event) {
359 :
360 : // Get junction chains.
361 0 : event.saveJunctionSize();
362 0 : vector<vector<int> > junChains = colTrace.getJunChains(event);
363 :
364 : // Remove junctions.
365 0 : vector<int> junRem;
366 0 : for (int i = 0; i < int(junChains.size()); ++i) {
367 0 : if (junChains[i].size() < 3)
368 : continue;
369 :
370 0 : vector<int> cols, acols;
371 0 : for (int j = 0; j < int(junChains[i].size()); ++j) {
372 :
373 0 : junRem.push_back(junChains[i][j]);
374 0 : if (event.kindJunction(junChains[i][j]) % 2 == 0)
375 0 : for (int jLeg = 0; jLeg < 3; ++jLeg)
376 0 : acols.push_back(event.colJunction(junChains[i][j],jLeg));
377 : else
378 0 : for (int jLeg = 0; jLeg < 3; ++jLeg)
379 0 : cols.push_back(event.colJunction(junChains[i][j],jLeg));
380 : }
381 :
382 0 : for (int j = 0; j < int(cols.size()); ++j)
383 0 : for (int k = 0; k < int(acols.size()); ++k)
384 0 : if (cols[j] == acols[k]) {
385 0 : cols.erase(cols.begin() + j);
386 0 : acols.erase(acols.begin() + k);
387 0 : j--;
388 0 : break;
389 : }
390 :
391 : // Find junctions if we have more colours than anti colours
392 0 : while (cols.size() > acols.size()) {
393 0 : int i1 = int(rndmPtr->flat() *cols.size());
394 0 : int col1 = cols[i1];
395 0 : cols.erase(cols.begin() + i1);
396 0 : int i2 = int(rndmPtr->flat() *cols.size());
397 0 : int col2 = cols[i2];
398 0 : cols.erase(cols.begin() + i2);
399 0 : int i3 = int(rndmPtr->flat() *cols.size());
400 0 : int col3 = cols[i3];
401 0 : cols.erase(cols.begin() + i3);
402 0 : event.appendJunction(1, col1, col2, col3);
403 : }
404 :
405 : // Find junctions if we have more colours than anti colours
406 0 : while (acols.size() > cols.size()) {
407 0 : int i1 = int(rndmPtr->flat() *acols.size());
408 0 : int acol1 = acols[i1];
409 0 : acols.erase(acols.begin() + i1);
410 0 : int i2 = int(rndmPtr->flat() *acols.size());
411 0 : int acol2 = acols[i2];
412 0 : acols.erase(acols.begin() + i2);
413 0 : int i3 = int(rndmPtr->flat() *acols.size());
414 0 : int acol3 = acols[i3];
415 0 : acols.erase(acols.begin() + i3);
416 0 : event.appendJunction(2,acol1,acol2,acol3);
417 : }
418 :
419 : // If we have more than two colour anti colour pairs
420 : // form junction anti junction pair.
421 0 : while (int(acols.size()) > 1) {
422 0 : int i1 = int(rndmPtr->flat() *cols.size());
423 0 : int col1 = cols[i1];
424 0 : cols.erase(cols.begin() + i1);
425 0 : int i2 = int(rndmPtr->flat() *cols.size());
426 0 : int col2 = cols[i2];
427 0 : cols.erase(cols.begin() + i2);
428 0 : int i3 = int(rndmPtr->flat() *acols.size());
429 0 : int acol1 = acols[i3];
430 0 : acols.erase(acols.begin() + i3);
431 0 : int i4 = int(rndmPtr->flat() *acols.size());
432 0 : int acol2 = acols[i4];
433 0 : acols.erase(acols.begin() + i4);
434 0 : int newCol = event.nextColTag();
435 0 : event.appendJunction(1, col1, col2, newCol);
436 0 : event.appendJunction(2, acol1, acol2, newCol);
437 : }
438 :
439 : // If we have one colour and one anti colour, form normal string.
440 0 : if (int(acols.size()) == 1) {
441 : int iCol = -1;
442 0 : for (int iPar = 0; iPar < event.size(); ++iPar)
443 0 : if (event[iPar].isFinal() && event[iPar].col() == cols[0])
444 0 : iCol = iPar;
445 0 : if (iCol == -1) {
446 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunChain:"
447 : "Splitting multiple directly connected junctions failed");
448 0 : return false;
449 : }
450 0 : int iNew = event.copy(iCol, 76);
451 0 : event[iNew].col(acols[0]);
452 0 : }
453 0 : }
454 :
455 : // Delete the junctions from the event record.
456 0 : sort(junRem.begin(),junRem.end());
457 0 : reverse(junRem.begin(),junRem.end());
458 0 : for (int i = 0; i < int(junRem.size()); ++i)
459 0 : event.eraseJunction(junRem[i]);
460 0 : event.saveJunctionSize();
461 :
462 0 : return true;
463 0 : }
464 :
465 : //--------------------------------------------------------------------------
466 :
467 : // Split junction pairs.
468 : // If it has 3 connections just ignore junctions.
469 : // If it has 2 connections colapse into single string.
470 : // If it has 1 connection, depend on the string length.
471 :
472 : bool JunctionSplitting::splitJunPairs(Event& event,
473 : vector<vector< int > >& iPartonJun, vector<vector< int > >& iPartonAntiJun) {
474 :
475 : // Clear old memory.
476 0 : event.saveJunctionSize();
477 0 : vector<int> junRem;
478 :
479 : // Get junction chains.
480 0 : vector<vector<int> > junChains = colTrace.getJunChains(event);
481 :
482 0 : for (int i = 0; i < int(junChains.size()); ++i) {
483 0 : if (junChains[i].size() == 2) {
484 0 : vector<pair<int,int> > matchedLegs;
485 0 : for (int j = 0; j < 3; ++j)
486 0 : for (int k = 0; k < 3; ++k)
487 0 : if (event.colJunction(junChains[i][0],j) ==
488 0 : event.colJunction(junChains[i][1],k))
489 0 : matchedLegs.push_back(make_pair(j,k));
490 :
491 : // For three connected legs, just remove the junctions,
492 : // since the pair contains no energy/momentum.
493 0 : if (matchedLegs.size() == 3) {
494 0 : junRem.push_back(junChains[i][0]);
495 0 : junRem.push_back(junChains[i][1]);
496 :
497 : }
498 :
499 : // Split into string if two legs are combined.
500 0 : else if (matchedLegs.size() == 2) {
501 :
502 : // Find first leg.
503 : int i1 = 0;
504 0 : if (matchedLegs[0].first != 1 && matchedLegs[1].first != 1) i1 = 1;
505 0 : if (matchedLegs[0].first != 2 && matchedLegs[1].first != 2) i1 = 2;
506 :
507 : // Find second leg.
508 : int j1 = 0;
509 0 : if (matchedLegs[0].second != 1 && matchedLegs[1].second != 1) j1 = 1;
510 0 : if (matchedLegs[0].second != 2 && matchedLegs[1].second != 2) j1 = 2;
511 :
512 : // Find corresponding colours.
513 0 : int col = event.colJunction(junChains[i][0],i1);
514 0 : int acol = event.colJunction(junChains[i][1],j1);
515 0 : if (event.kindJunction(junChains[i][1]) % 2 == 1)
516 0 : swap(col,acol);
517 :
518 : // Find index of anti particle.
519 : int iAcol = -1;
520 0 : for (int j = 0;j < event.size();++j)
521 0 : if (event[j].isFinal() && event[j].acol() == acol) {
522 : iAcol = j;
523 0 : break;
524 : }
525 0 : if (iAcol == -1) {
526 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunChain:"
527 : "Anti colour not found when combining two junctions to a string");
528 0 : return false;
529 : }
530 :
531 : // Update anti colour of anti particle.
532 0 : int iNew = event.copy(iAcol,66);
533 0 : event[iNew].acol(col);
534 :
535 : // Remove the junctions from the event record.
536 0 : junRem.push_back(junChains[i][0]);
537 0 : junRem.push_back(junChains[i][1]);
538 0 : }
539 :
540 : // Split into string if two legs are combined.
541 0 : else if (matchedLegs.size() == 1) {
542 :
543 : // store junction numbers.
544 0 : int iJun = junChains[i][0];
545 0 : int iAnti = junChains[i][1];
546 0 : int iLeg = matchedLegs[0].first;
547 0 : int iAntiLeg = matchedLegs[0].second;
548 0 : if (event.kindJunction(iAnti) % 2 == 1) {
549 0 : swap(iJun, iAnti);
550 0 : swap(iLeg, iAntiLeg);
551 0 : }
552 :
553 : // Find the junctions in the parton list.
554 : int iJunList = -1, iAntiList = -1;
555 0 : for (int l = 0;l < int(iPartonJun.size()); ++l)
556 0 : if (- iPartonJun[l][0]/10 - 1 == iJun) {
557 : iJunList = l;
558 0 : break;
559 : }
560 :
561 0 : for (int l = 0;l < int(iPartonAntiJun.size()); ++l)
562 0 : if (- iPartonAntiJun[l][0]/10 - 1 == iAnti) {
563 : iAntiList = l;
564 0 : break;
565 : }
566 :
567 : // Fill in vector of the legs content.
568 0 : vector<vector <int> > iJunLegs;
569 0 : iJunLegs.resize(3);
570 : int leg = -1;
571 :
572 0 : for (int l = 0; l < int(iPartonJun[iJunList].size()); ++l) {
573 0 : if ( iPartonJun[iJunList][l]/10 == iPartonJun[iJunList][0]/10) ++leg;
574 0 : iJunLegs[leg].push_back(iPartonJun[iJunList][l]);
575 : }
576 :
577 : // Fill in vector of the legs content.
578 0 : vector<vector <int> > iAntiLegs;
579 0 : iAntiLegs.resize(3);
580 : leg = -1;
581 0 : for (int l = 0; l < int(iPartonAntiJun[iAntiList].size()); ++l) {
582 0 : if ( iPartonAntiJun[iAntiList][l]/10
583 0 : == iPartonAntiJun[iAntiList][0]/10) ++leg;
584 0 : iAntiLegs[leg].push_back(iPartonAntiJun[iAntiList][l]);
585 : }
586 :
587 : // Identify the two external legs of either junction.
588 0 : vector<int>& iJunLeg0 = (iLeg == 0) ? iJunLegs[1] : iJunLegs[0];
589 0 : vector<int>& iJunLeg1 = (iLeg == 2) ? iJunLegs[1] : iJunLegs[2];
590 0 : vector<int>& iAntiLeg0 = (iAntiLeg == 0) ? iAntiLegs[1] : iAntiLegs[0];
591 0 : vector<int>& iAntiLeg1 = (iAntiLeg == 2) ? iAntiLegs[1] : iAntiLegs[2];
592 :
593 : // Simplified procedure: mainly study first parton on each leg.
594 0 : Vec4 pJunLeg0 = event[ iJunLeg0[1] ].p();
595 0 : Vec4 pJunLeg1 = event[ iJunLeg1[1] ].p();
596 0 : Vec4 pAntiLeg0 = event[ iAntiLeg0[1] ].p();
597 0 : Vec4 pAntiLeg1 = event[ iAntiLeg1[1] ].p();
598 :
599 0 : if ( REtaPhi(pJunLeg0,pJunLeg1) < MINDELTAR ||
600 0 : REtaPhi(pAntiLeg0,pAntiLeg1) < MINDELTAR ||
601 0 : REtaPhi(pJunLeg0,pAntiLeg0) < MINDELTAR ||
602 0 : REtaPhi(pJunLeg0,pAntiLeg1) < MINDELTAR ||
603 0 : REtaPhi(pJunLeg1,pAntiLeg0) < MINDELTAR ||
604 0 : REtaPhi(pJunLeg1,pAntiLeg1) < MINDELTAR) {
605 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunPairs: "
606 : "Parallel junction state not allowed.");
607 0 : return false;
608 : }
609 :
610 : // Starting frame hopefully intermediate to two junction directions.
611 0 : Vec4 pStart = pJunLeg0 / pJunLeg0.e() + pJunLeg1 / pJunLeg1.e()
612 0 : + pAntiLeg0 / pAntiLeg0.e() + pAntiLeg1 / pAntiLeg1.e();
613 :
614 : // Loop over iteration to junction/antijunction rest frames (JRF/ARF).
615 0 : RotBstMatrix MtoJRF, MtoARF;
616 0 : Vec4 pInJRF[3], pInARF[3];
617 0 : for (int iJunLocal = 0; iJunLocal < 2; ++iJunLocal) {
618 0 : int offset = (iJunLocal == 0) ? 0 : 2;
619 :
620 : // Iterate from system rest frame towards the junction rest frame.
621 0 : RotBstMatrix MtoRF, Mstep;
622 0 : MtoRF.bstback(pStart);
623 0 : Vec4 pInRF[4];
624 : int iter = 0;
625 0 : do {
626 0 : ++iter;
627 :
628 : // Find rest-frame momenta on the three sides of the junction.
629 : // Only consider first parton on each leg, for simplicity.
630 0 : pInRF[0 + offset] = pJunLeg0;
631 0 : pInRF[1 + offset] = pJunLeg1;
632 0 : pInRF[2 - offset] = pAntiLeg0;
633 0 : pInRF[3 - offset] = pAntiLeg1;
634 0 : for (int l = 0; l < 4; ++l) pInRF[l].rotbst(MtoRF);
635 :
636 : // For third side add both legs beyond other junction, weighted.
637 : double wt2 = -1e5;
638 0 : if(-pInRF[2].e()/eNormJunction < 10) wt2 = 1. - exp( -pInRF[2].e() / eNormJunction);
639 : double wt3 = -1e5;
640 0 : if(-pInRF[3].e()/eNormJunction < 10) wt3 = 1. - exp( -pInRF[3].e() / eNormJunction);
641 0 : pInRF[2] = wt2 * pInRF[2] + wt3 * pInRF[3];
642 :
643 : // Find new junction rest frame from the set of momenta.
644 0 : Mstep = stringFrag.junctionRestFrame( pInRF[0], pInRF[1], pInRF[2]);
645 0 : MtoRF.rotbst( Mstep );
646 0 : } while (iter < 3 || (Mstep.deviation() > CONVJNREST
647 0 : && iter < NTRYJNREST) );
648 :
649 : // Store final boost and rest-frame (weighted) momenta.
650 0 : if (iJunLocal == 0) {
651 0 : MtoJRF = MtoRF;
652 0 : for (int l = 0; l < 3; ++l) pInJRF[l] = pInRF[l];
653 0 : } else {
654 0 : MtoARF = MtoRF;
655 0 : for (int l = 0; l < 3; ++l) pInARF[l] = pInRF[l];
656 : }
657 0 : }
658 :
659 : // Opposite operations: boost from JRF/ARF to original system.
660 0 : RotBstMatrix MfromJRF = MtoJRF;
661 0 : MfromJRF.invert();
662 0 : RotBstMatrix MfromARF = MtoARF;
663 0 : MfromARF.invert();
664 :
665 : // Velocity vectors of junctions and momentum of legs in lab frame.
666 0 : Vec4 vJun(0., 0., 0., 1.);
667 0 : vJun.rotbst(MfromJRF);
668 0 : Vec4 vAnti(0., 0., 0., 1.);
669 0 : vAnti.rotbst(MfromARF);
670 0 : Vec4 pLabJ[3], pLabA[3];
671 0 : for (int l = 0; l < 3; ++l) {
672 0 : pLabJ[l] = pInJRF[l];
673 0 : pLabJ[l].rotbst(MfromJRF);
674 0 : pLabA[l] = pInARF[l];
675 0 : pLabA[l].rotbst(MfromARF);
676 : }
677 :
678 : // Calculate Lambda-measure length of three possible topologies.
679 0 : double vJvA = vJun * vAnti;
680 : double vJvAe2y = vJvA;
681 0 : if ((vJvA*vJvA - 1.) > 0.) vJvAe2y += sqrt(vJvA*vJvA - 1.);
682 0 : double lambdaJA = stringLength.getJuncLength(pInJRF[0], pInJRF[1],
683 0 : pInARF[0], pInARF[1]);
684 :
685 0 : double lambda00 = stringLength.getStringLength(pLabJ[0], pLabA[0]) +
686 0 : stringLength.getStringLength(pLabJ[1], pLabA[1]);
687 :
688 0 : double lambda01 = stringLength.getStringLength(pLabJ[0], pLabA[1]) +
689 0 : stringLength.getStringLength(pLabJ[1], pLabA[0]);
690 :
691 : // Case when either topology without junctions is the shorter one.
692 0 : if (lambdaJA > min( lambda00, lambda01) && allowDoubleJunRem) {
693 :
694 : // Find indices of particles.
695 0 : int iCol1 = iJunLeg0[1];
696 0 : int iCol2 = iJunLeg1[1];
697 0 : int iAcol1 = iAntiLeg0[1];
698 0 : int iAcol2 = iAntiLeg1[1];
699 0 : if (lambda00 > lambda01)
700 0 : swap(iAcol1, iAcol2);
701 :
702 : // Change the colour index and mark junctions to be removed.
703 0 : int iNew1 = event.copy(iAcol1, 66);
704 0 : event[iNew1].acol(event[iCol1].col());
705 0 : int iNew2 = event.copy(iAcol2, 66);
706 0 : event[iNew2].acol(event[iCol2].col());
707 0 : junRem.push_back(junChains[i][0]);
708 0 : junRem.push_back(junChains[i][1]);
709 : continue;
710 0 : }
711 :
712 : // Case where junction and antijunction to be separated.
713 : // Shuffle (p+/p-) momentum of order <mThad> between systems,
714 : // times 2/3 for 120 degree in JRF, times 1/2 for two legs,
715 : // but not more than half of what nearest parton carries.
716 :
717 : // Only allow to take momentum from non-heavy resonances
718 : // (i.e. id 1-5 and 21 and diquarks). If none is available return false.
719 0 : int idJ0Abs = event[ iJunLeg0[1] ].idAbs();
720 0 : bool acceptJ0 = idJ0Abs < 6 || idJ0Abs == 21 ||
721 0 : (idJ0Abs > 1000 && idJ0Abs < 6000 && (idJ0Abs / 10) % 10 == 0);
722 0 : int idJ1Abs = event[ iJunLeg1[1] ].idAbs();
723 0 : bool acceptJ1 = idJ1Abs < 6 || idJ1Abs == 21 ||
724 0 : (idJ1Abs > 1000 && idJ1Abs < 6000 && (idJ1Abs / 10) % 10 == 0);
725 0 : int idA0Abs = event[ iAntiLeg0[1] ].idAbs();
726 0 : bool acceptA0 = idA0Abs < 6 || idA0Abs == 21 ||
727 0 : (idA0Abs > 1000 && idA0Abs < 6000 && (idA0Abs / 10) % 10 == 0);
728 0 : int idA1Abs = event[ iAntiLeg1[1] ].idAbs();
729 0 : bool acceptA1 = idA1Abs < 6 || idA1Abs == 21 ||
730 0 : (idA1Abs > 1000 && idA1Abs < 6000 && (idA1Abs / 10) % 10 == 0);
731 :
732 0 : if ( !(acceptJ0 || acceptJ1)) {
733 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunPairs: "
734 : "No light quarks available in junction split.");
735 0 : return false;
736 : }
737 :
738 0 : if ( !(acceptA0 || acceptA1)) {
739 0 : infoPtr->errorMsg("Warning in JunctionSplitting::SplitJunPairs: "
740 : "No light quarks available in junction split.");
741 0 : return false;
742 : }
743 :
744 0 : double eShift = MTHAD / (3. * sqrt(vJvAe2y));
745 : double fracJ0 = 0, fracJ1 = 0, fracA0 = 0, fracA1 = 0;
746 0 : if (acceptJ0)
747 0 : fracJ0 = min(0.5, eShift / pInJRF[0].e());
748 0 : if (acceptJ1)
749 0 : fracJ1 = min(0.5, eShift / pInJRF[1].e());
750 0 : Vec4 pFromJun = fracJ0 * pJunLeg0 + fracJ1 * pJunLeg1;
751 0 : if (acceptA0)
752 0 : fracA0 = min(0.5, eShift / pInARF[0].e());
753 0 : if (acceptA1)
754 0 : fracA1 = min(0.5, eShift / pInARF[1].e());
755 0 : Vec4 pFromAnti = fracA0 * pAntiLeg0 + fracA1 * pAntiLeg1;
756 :
757 : // Pick a new quark at random; for simplicity no diquarks.
758 0 : int idQ = flavSel.pickLightQ();
759 :
760 0 : int junMother1 = min(iJunLeg0[1], iJunLeg1[1]);
761 0 : int junMother2 = max(iJunLeg0[1], iJunLeg1[1]);
762 0 : int antiMother1 = min(iAntiLeg0[1], iAntiLeg1[1]);
763 0 : int antiMother2 = max(iAntiLeg0[1], iAntiLeg1[1]);
764 :
765 : // Copy junction partons with scaled-down momenta and update legs.
766 0 : int iJunNew1 = event.copy(iJunLeg0[1], 76);
767 0 : event[iJunNew1].rescale5(1. - fracJ0);
768 0 : iJunLeg0[1] = iJunNew1;
769 0 : event[iJunNew1].mothers(junMother2, junMother1);
770 :
771 0 : int iJunNew2 = event.copy(iJunLeg1[1], 76);
772 0 : event[iJunNew2].rescale5(1. - fracJ1);
773 0 : iJunLeg1[1] = iJunNew2;
774 0 : event[iJunNew2].mothers(junMother2, junMother1);
775 :
776 : // Update antijunction anticolour and store antiquark with junction
777 : // momentum.
778 0 : int acolQ = event.nextColTag();
779 0 : event.endColJunction(iAnti, iAntiLeg, acolQ);
780 0 : event.colJunction(iAnti, iAntiLeg, acolQ);
781 0 : int iNewA = event.append( -idQ, 76, junMother2, junMother1, 0, 0,
782 0 : 0, acolQ, pFromJun, pFromJun.mCalc() );
783 :
784 : // Copy anti junction partons with scaled-down momenta and update legs.
785 0 : int iAntiNew1 = event.copy(iAntiLeg0[1], 76);
786 0 : event[iAntiNew1].rescale5(1. - fracA0);
787 0 : iAntiLeg0[1] = iAntiNew1;
788 0 : event[iAntiNew1].mothers(antiMother2, antiMother1);
789 :
790 0 : int iAntiNew2 = event.copy(iAntiLeg1[1], 76);
791 0 : event[iAntiNew2].rescale5(1. - fracA1);
792 0 : iAntiLeg1[1] = iAntiNew2;
793 0 : event[iAntiNew2].mothers(antiMother2, antiMother1);
794 :
795 : // Update junction colour and store quark with antijunction momentum.
796 0 : int colQ = event.nextColTag();
797 0 : event.endColJunction(iJun, iLeg, colQ);
798 0 : event.colJunction(iJun, iLeg, colQ);
799 0 : int iNewJ = event.append( idQ, 76, antiMother2, antiMother1, 0, 0,
800 0 : colQ, 0, pFromAnti, pFromAnti.mCalc() );
801 :
802 : // Set daughters.
803 0 : event[event[iJunNew1].mother1()].daughters( iJunNew1, iNewA);
804 0 : event[event[iJunNew1].mother2()].daughters( iJunNew1, iNewA);
805 0 : event[event[iAntiNew1].mother1()].daughters( iAntiNew1, iNewJ);
806 0 : event[event[iAntiNew1].mother2()].daughters( iAntiNew1, iNewJ);
807 :
808 : // Done with splitting junction from antijunction.
809 0 : }
810 0 : }
811 : }
812 :
813 : // Delete the junctions from the event record.
814 0 : sort(junRem.begin(),junRem.end());
815 0 : reverse(junRem.begin(),junRem.end());
816 0 : for (int i = 0; i < int(junRem.size()); ++i)
817 0 : event.eraseJunction(junRem[i]);
818 0 : event.saveJunctionSize();
819 :
820 : // Done.
821 0 : return true;
822 0 : }
823 :
824 : //--------------------------------------------------------------------------
825 :
826 : // Get the list of partons connected to the junctions.
827 :
828 : bool JunctionSplitting::getPartonLists(Event& event,
829 : vector<vector< int > > & iPartonJun, vector<vector<int > >& iPartonAntiJun) {
830 :
831 : // Need to try and split junction structures.
832 0 : colTrace.setupColList(event);
833 0 : vector<int> iParton;
834 0 : iPartonJun.clear();
835 0 : iPartonAntiJun.clear();
836 :
837 : // Loop over junctions and collect all junctions.
838 : // Then afterwards collect all anti junctions.
839 : // This ensures that all gluons are collected on the junctions.
840 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
841 0 : if (event.remainsJunction(iJun)) {
842 :
843 0 : int kindJun = event.kindJunction(iJun);
844 0 : iParton.resize(0);
845 :
846 0 : if (kindJun % 2 != 1) continue;
847 :
848 : // Loop over junction legs.
849 0 : for (int iCol = 0; iCol < 3; ++iCol) {
850 0 : int indxCol = event.colJunction(iJun, iCol);
851 0 : iParton.push_back( -(10 + 10 * iJun + iCol) );
852 : // Junctions: find color ends.
853 0 : if (kindJun % 2 == 1 && !colTrace.traceFromAcol(indxCol,
854 : event, iJun, iCol, iParton))
855 0 : return false;
856 0 : }
857 :
858 : // Store the anti junction and junction list.
859 : int nNeg = 0;
860 0 : for (int i = 0; i < int(iParton.size()); ++i) if (iParton[i] < 0)
861 0 : ++nNeg;
862 0 : if (nNeg > 3 )
863 0 : iPartonJun.push_back(iParton);
864 0 : }
865 :
866 : // Loop over all anti junctions.
867 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
868 0 : if (event.remainsJunction(iJun)) {
869 :
870 0 : int kindJun = event.kindJunction(iJun);
871 0 : iParton.resize(0);
872 :
873 0 : if (kindJun % 2 != 0)
874 0 : continue;
875 :
876 : // Loop over junction legs.
877 0 : for (int iCol = 0; iCol < 3; ++iCol) {
878 0 : int indxCol = event.colJunction(iJun, iCol);
879 0 : iParton.push_back( -(10 + 10 * iJun + iCol) );
880 : // Antijunctions: find anticolor ends.
881 0 : if (kindJun % 2 == 0 && !colTrace.traceFromCol(indxCol,
882 : event, iJun, iCol, iParton))
883 0 : return false;
884 0 : }
885 :
886 : // Store the anti junction and junction list.
887 : int nNeg = 0;
888 0 : for (int i = 0; i < int(iParton.size()); ++i) if (iParton[i] < 0)
889 0 : ++nNeg;
890 0 : if (nNeg > 3 )
891 0 : iPartonAntiJun.push_back(iParton);
892 0 : }
893 :
894 : // Done.
895 0 : return true;
896 :
897 0 : }
898 :
899 : //--------------------------------------------------------------------------
900 :
901 : // Change the anticolour of the particle that has acol to be col.
902 :
903 : bool JunctionSplitting::setAcol(Event& event, int col, int acol) {
904 :
905 : // Update anticolour if it belongs to a particle.
906 0 : for (int j = 0;j < event.size(); ++j)
907 0 : if (event[j].isFinal() && event[j].acol() == acol) {
908 0 : int iNew = event.copy(j,66);
909 0 : event[iNew].acol(col);
910 : return true;
911 : }
912 : // Check if anti junction is connected to a junction.
913 0 : for (int j = 0;j < event.sizeJunction(); ++j)
914 0 : for (int jLeg = 0;jLeg < 3; ++jLeg)
915 0 : if (event.colJunction(j, jLeg) == acol) {
916 0 : event.colJunction(j, jLeg, col);
917 0 : return true;
918 : }
919 :
920 : // If no acol was found something went wrong.
921 0 : infoPtr->errorMsg("Warning in JunctionSplitting::setAcol:"
922 : "Anti colour not found when combing two junctions to a string");
923 0 : return false;
924 0 : }
925 :
926 : //==========================================================================
927 :
928 : } // end namespace Pythia8
|
