1
// HadronLevel.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.
6
// Function definitions (not found in the header) for the HadronLevel class.
8
#include "HadronLevel.h"
12
//==========================================================================
14
// The HadronLevel class.
16
//--------------------------------------------------------------------------
18
// Constants: could be changed here if desired, but normally should not.
19
// These are of technical nature, as described for each.
21
// For breaking J-J string, pick a Gamma by taking a step with fictitious mass.
22
const double HadronLevel::JJSTRINGM2MAX = 25.;
23
const double HadronLevel::JJSTRINGM2FRAC = 0.1;
25
// Iterate junction rest frame boost until convergence or too many tries.
26
const double HadronLevel::CONVJNREST = 1e-5;
27
const int HadronLevel::NTRYJNREST = 20;
29
// Typical average transvere primary hadron mass <mThad>.
30
const double HadronLevel::MTHAD = 0.9;
32
//--------------------------------------------------------------------------
34
// Find settings. Initialize HadronLevel classes as required.
36
bool HadronLevel::init(Info* infoPtrIn, Settings& settings,
37
ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
38
Couplings* couplingsPtrIn, TimeShower* timesDecPtr,
39
RHadrons* rHadronsPtrIn, DecayHandler* decayHandlePtr,
40
vector<int> handledParticles) {
44
particleDataPtr = particleDataPtrIn;
46
couplingsPtr = couplingsPtrIn;
47
rHadronsPtr = rHadronsPtrIn;
50
doHadronize = settings.flag("HadronLevel:Hadronize");
51
doDecay = settings.flag("HadronLevel:Decay");
52
doBoseEinstein = settings.flag("HadronLevel:BoseEinstein");
54
// Boundary mass between string and ministring handling.
55
mStringMin = settings.parm("HadronLevel:mStringMin");
57
// For junction processing.
58
eNormJunction = settings.parm("StringFragmentation:eNormJunction");
60
// Allow R-hadron formation.
61
allowRH = settings.flag("RHadrons:allow");
63
// Particles that should decay or not before Bose-Einstein stage.
64
widthSepBE = settings.parm("BoseEinstein:widthSep");
66
// Hadron scattering --rjc
67
doHadronScatter = settings.flag("HadronScatter:scatter");
68
hsAfterDecay = settings.flag("HadronScatter:afterDecay");
70
// Initialize auxiliary fragmentation classes.
71
flavSel.init(settings, rndmPtr);
72
pTSel.init(settings, *particleDataPtr, rndmPtr);
73
zSel.init(settings, *particleDataPtr, rndmPtr);
75
// Initialize auxiliary administrative class.
76
colConfig.init(infoPtr, settings, &flavSel);
78
// Initialize string and ministring fragmentation.
79
stringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr,
80
&flavSel, &pTSel, &zSel);
81
ministringFrag.init(infoPtr, settings, particleDataPtr, rndmPtr,
82
&flavSel, &pTSel, &zSel);
84
// Initialize particle decays.
85
decays.init(infoPtr, settings, particleDataPtr, rndmPtr, couplingsPtr,
86
timesDecPtr, &flavSel, decayHandlePtr, handledParticles);
88
// Initialize BoseEinstein.
89
boseEinstein.init(infoPtr, settings, *particleDataPtr);
91
// Initialize HadronScatter --rjc
93
hadronScatter.init(infoPtr, settings, rndmPtr, particleDataPtr);
95
// Initialize Hidden-Valley fragmentation, if necessary.
96
useHiddenValley = hiddenvalleyFrag.init(infoPtr, settings,
97
particleDataPtr, rndmPtr);
99
// Send flavour and z selection pointers to R-hadron machinery.
100
rHadronsPtr->fragPtrs( &flavSel, &zSel);
107
//--------------------------------------------------------------------------
109
// Hadronize and decay the next parton-level.
111
bool HadronLevel::next( Event& event) {
113
// Do Hidden-Valley fragmentation, if necessary.
114
if (useHiddenValley) hiddenvalleyFrag.fragment(event);
116
// Colour-octet onia states must be decayed to singlet + gluon.
117
if (!decayOctetOnia(event)) return false;
119
// Possibility of hadronization inside decay, but then no BE second time.
120
// Hadron scattering, first pass only --rjc
121
bool moreToDo, firstPass = true;
122
bool doBoseEinsteinNow = doBoseEinstein;
126
// First part: string fragmentation.
129
// Find the complete colour singlet configuration of the event.
130
if (!findSinglets( event)) return false;
132
// Fragment off R-hadrons, if necessary.
133
if (allowRH && !rHadronsPtr->produce( colConfig, event))
136
// Process all colour singlet (sub)system
137
for (int iSub = 0; iSub < colConfig.size(); ++iSub) {
139
// Collect sequentially all partons in a colour singlet subsystem.
140
colConfig.collect(iSub, event);
142
// String fragmentation of each colour singlet (sub)system.
143
if ( colConfig[iSub].massExcess > mStringMin ) {
144
if (!stringFrag.fragment( iSub, colConfig, event)) return false;
146
// Low-mass string treated separately. Tell if diffractive system.
148
bool isDiff = infoPtr->isDiffractiveA()
149
|| infoPtr->isDiffractiveB();
150
if (!ministringFrag.fragment( iSub, colConfig, event, isDiff))
156
// Hadron scattering --rjc
157
if (doHadronScatter && !hsAfterDecay && firstPass)
158
hadronScatter.scatter(event);
160
// Second part: sequential decays of short-lived particles (incl. K0).
163
// Loop through all entries to find those that should decay.
166
Particle& decayer = event[iDec];
167
if ( decayer.isFinal() && decayer.canDecay() && decayer.mayDecay()
168
&& (decayer.mWidth() > widthSepBE || decayer.idAbs() == 311) ) {
169
decays.decay( iDec, event);
170
if (decays.moreToDo()) moreToDo = true;
172
} while (++iDec < event.size());
175
// Hadron scattering --rjc
176
if (doHadronScatter && hsAfterDecay && firstPass)
177
hadronScatter.scatter(event);
179
// Third part: include Bose-Einstein effects among current particles.
180
if (doBoseEinsteinNow) {
181
if (!boseEinstein.shiftEvent(event)) return false;
182
doBoseEinsteinNow = false;
185
// Fourth part: sequential decays also of long-lived particles.
188
// Loop through all entries to find those that should decay.
191
Particle& decayer = event[iDec];
192
if ( decayer.isFinal() && decayer.canDecay() && decayer.mayDecay() ) {
193
decays.decay( iDec, event);
194
if (decays.moreToDo()) moreToDo = true;
196
} while (++iDec < event.size());
199
// Normally done first time around, but sometimes not (e.g. Upsilon).
207
//--------------------------------------------------------------------------
209
// Allow more decays if on/off switches changed.
210
// Note: does not do sequential hadronization, e.g. for Upsilon.
212
bool HadronLevel::moreDecays( Event& event) {
214
// Colour-octet onia states must be decayed to singlet + gluon.
215
if (!decayOctetOnia(event)) return false;
217
// Loop through all entries to find those that should decay.
220
if ( event[iDec].isFinal() && event[iDec].canDecay()
221
&& event[iDec].mayDecay() ) decays.decay( iDec, event);
222
} while (++iDec < event.size());
229
//--------------------------------------------------------------------------
231
// Decay colour-octet onium states.
233
bool HadronLevel::decayOctetOnia(Event& event) {
235
// Onium states to be decayed.
236
int idOnium[6] = { 9900443, 9900441, 9910441,
237
9900553, 9900551, 9910551 };
239
// Loop over particles and identify onia.
240
for (int iDec = 0; iDec < event.size(); ++iDec)
241
if (event[iDec].isFinal()) {
242
int id = event[iDec].id();
243
bool isOnium = false;
244
for (int j = 0; j < 6; ++j) if (id == idOnium[j]) isOnium = true;
246
// Decay any onia encountered.
248
if (!decays.decay( iDec, event)) return false;
250
// Set colour flow by hand: gluon inherits octet-onium state.
251
int iGlu = event.size() - 1;
252
event[iGlu].cols( event[iDec].col(), event[iDec].acol() );
261
//--------------------------------------------------------------------------
263
// Trace colour flow in the event to form colour singlet subsystems.
265
bool HadronLevel::findSinglets(Event& event) {
267
// Find a list of final partons and of all colour ends and gluons.
270
iColAndAcol.resize(0);
271
for (int i = 0; i < event.size(); ++i) if (event[i].isFinal()) {
272
if (event[i].col() > 0 && event[i].acol() > 0) iColAndAcol.push_back(i);
273
else if (event[i].col() > 0) iColEnd.push_back(i);
274
else if (event[i].acol() > 0) iAcolEnd.push_back(i);
277
// Begin arrange the partons into separate colour singlets.
279
iPartonJun.resize(0);
280
iPartonAntiJun.resize(0);
282
// Junctions: loop over them, and identify kind.
283
for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
284
if (event.remainsJunction(iJun)) {
285
event.remainsJunction(iJun, false);
286
int kindJun = event.kindJunction(iJun);
289
// Loop over junction legs.
290
for (int iCol = 0; iCol < 3; ++iCol) {
291
int indxCol = event.colJunction(iJun, iCol);
292
iParton.push_back( -(10 + 10 * iJun + iCol) );
293
// Junctions: find color ends.
294
if (kindJun % 2 == 1 && !traceFromAcol(indxCol, event, iJun, iCol))
296
// Antijunctions: find anticolor ends.
297
if (kindJun % 2 == 0 && !traceFromCol(indxCol, event, iJun, iCol))
301
// Reject triple- and higher-junction systems (physics not implemented).
303
for (unsigned int i=0; i<iParton.size(); ++i) if (iParton[i]<0) ++nJun;
305
infoPtr->errorMsg("Error in HadronLevel::findSinglets: "
306
"too many junction-junction connections");
310
// Keep in memory a junction hooked up with an antijunction,
311
// else store found single-junction system.
313
for (int i = 0; i < int(iParton.size()); ++i) if (iParton[i] < 0)
315
if (nNeg > 3 && kindJun % 2 == 1) {
316
for (int i = 0; i < int(iParton.size()); ++i)
317
iPartonJun.push_back(iParton[i]);
318
} else if (nNeg > 3 && kindJun % 2 == 0) {
319
for (int i = 0; i < int(iParton.size()); ++i)
320
iPartonAntiJun.push_back(iParton[i]);
322
// A junction may be eliminated by insert if two quarks are nearby.
323
int nJunOld = event.sizeJunction();
324
if (!colConfig.insert(iParton, event)) return false;
325
if (event.sizeJunction() < nJunOld) --iJun;
329
// Split junction-antijunction system into two, and store those.
330
// (Only one system in extreme cases, and then second empty.)
331
if (iPartonJun.size() > 0 && iPartonAntiJun.size() > 0) {
332
if (!splitJunctionPair(event)) return false;
333
if (!colConfig.insert(iPartonJun, event)) return false;
334
if (iPartonAntiJun.size() > 0)
335
if (!colConfig.insert(iPartonAntiJun, event)) return false;
336
// Error if only one of junction and antijuction left here.
337
} else if (iPartonJun.size() > 0 || iPartonAntiJun.size() > 0) {
338
infoPtr->errorMsg("Error in HadronLevel::findSinglets: "
339
"unmatched (anti)junction");
343
// Open strings: pick up each colour end and trace to its anticolor end.
344
for (int iEnd = 0; iEnd < int(iColEnd.size()); ++iEnd) {
346
iParton.push_back( iColEnd[iEnd] );
347
int indxCol = event[ iColEnd[iEnd] ].col();
348
if (!traceFromCol(indxCol, event)) return false;
350
// Store found open string system. Analyze its properties.
351
if (!colConfig.insert(iParton, event)) return false;
354
// Closed strings : begin at any gluon and trace until back at it.
355
while (iColAndAcol.size() > 0) {
357
iParton.push_back( iColAndAcol[0] );
358
int indxCol = event[ iColAndAcol[0] ].col();
359
int indxAcol = event[ iColAndAcol[0] ].acol();
360
iColAndAcol[0] = iColAndAcol.back();
361
iColAndAcol.pop_back();
362
if (!traceInLoop(indxCol, indxAcol, event)) return false;
364
// Store found closed string system. Analyze its properties.
365
if (!colConfig.insert(iParton, event)) return false;
373
//--------------------------------------------------------------------------
375
// Trace a colour line, from a colour to an anticolour.
377
bool HadronLevel::traceFromCol(int indxCol, Event& event, int iJun,
380
// Junction kind, if any.
381
int kindJun = (iJun >= 0) ? event.kindJunction(iJun) : 0;
383
// Begin to look for a matching anticolour.
385
int loopMax = iColAndAcol.size() + 2;
386
bool hasFound = false;
391
// First check list of matching anticolour ends.
392
for (int i = 0; i < int(iAcolEnd.size()); ++i)
393
if (event[ iAcolEnd[i] ].acol() == indxCol) {
394
iParton.push_back( iAcolEnd[i] );
396
iAcolEnd[i] = iAcolEnd.back();
402
// Then check list of intermediate gluons.
404
for (int i = 0; i < int(iColAndAcol.size()); ++i)
405
if (event[ iColAndAcol[i] ].acol() == indxCol) {
406
iParton.push_back( iColAndAcol[i] );
408
// Update to new colour. Remove gluon.
409
indxCol = event[ iColAndAcol[i] ].col();
410
if (kindJun > 0) event.endColJunction(iJun, iCol, indxCol);
411
iColAndAcol[i] = iColAndAcol.back();
412
iColAndAcol.pop_back();
417
// In a pinch, check list of opposite-sign junction end colours.
418
// Store in iParton list as -(10 + 10 * iAntiJun + iAntiLeg).
419
if (!hasFound && kindJun % 2 == 0 && event.sizeJunction() > 1)
420
for (int iAntiJun = 0; iAntiJun < event.sizeJunction(); ++iAntiJun)
421
if (iAntiJun != iJun && event.kindJunction(iAntiJun) %2 == 1)
422
for (int iColAnti = 0; iColAnti < 3; ++iColAnti)
423
if (event.endColJunction(iAntiJun, iColAnti) == indxCol) {
424
iParton.push_back( -(10 + 10 * iAntiJun + iColAnti) );
430
// Keep on tracing via gluons until reached end of leg.
431
} while (hasFound && indxCol > 0 && loop < loopMax);
433
// Something went wrong in colour tracing.
434
if (!hasFound || loop == loopMax) {
435
infoPtr->errorMsg("Error in HadronLevel::traceFromCol: "
436
"colour tracing failed");
445
//--------------------------------------------------------------------------
447
// Trace a colour line, from an anticolour to a colour.
449
bool HadronLevel::traceFromAcol(int indxCol, Event& event, int iJun,
452
// Junction kind, if any.
453
int kindJun = (iJun >= 0) ? event.kindJunction(iJun) : 0;
455
// Begin to look for a matching colour.
457
int loopMax = iColAndAcol.size() + 2;
458
bool hasFound = false;
463
// First check list of matching colour ends.
464
for (int i = 0; i < int(iColEnd.size()); ++i)
465
if (event[ iColEnd[i] ].col() == indxCol) {
466
iParton.push_back( iColEnd[i] );
468
iColEnd[i] = iColEnd.back();
474
// Then check list of intermediate gluons.
476
for (int i = 0; i < int(iColAndAcol.size()); ++i)
477
if (event[ iColAndAcol[i] ].col() == indxCol) {
478
iParton.push_back( iColAndAcol[i] );
479
// Update to new colour. Remove gluon.
480
indxCol = event[ iColAndAcol[i] ].acol();
481
if (kindJun > 0) event.endColJunction(iJun, iCol, indxCol);
482
iColAndAcol[i] = iColAndAcol.back();
483
iColAndAcol.pop_back();
488
// In a pinch, check list of opposite-sign junction end colours.
489
// Store in iParton list as -(10 + 10 * iAntiJun + iLeg).
490
if (!hasFound && kindJun % 2 == 1 && event.sizeJunction() > 1)
491
for (int iAntiJun = 0; iAntiJun < event.sizeJunction(); ++iAntiJun)
492
if (iAntiJun != iJun && event.kindJunction(iAntiJun) % 2 == 0)
493
for (int iColAnti = 0; iColAnti < 3; ++iColAnti)
494
if (event.endColJunction(iAntiJun, iColAnti) == indxCol) {
495
iParton.push_back( -(10 + 10 * iAntiJun + iColAnti) );
501
// Keep on tracing via gluons until reached end of leg.
502
} while (hasFound && indxCol > 0 && loop < loopMax);
504
// Something went wrong in colour tracing.
505
if (!hasFound || loop == loopMax) {
506
infoPtr->errorMsg("Error in HadronLevel::traceFromAcol: "
507
"colour tracing failed");
516
//--------------------------------------------------------------------------
518
// Trace a colour loop, from a colour back to the anticolour of the same.
520
bool HadronLevel::traceInLoop(int indxCol, int indxAcol, Event& event) {
522
// Move around until back where begun.
524
int loopMax = iColAndAcol.size() + 2;
525
bool hasFound = false;
530
// Check list of gluons.
531
for (int i = 0; i < int(iColAndAcol.size()); ++i)
532
if (event[ iColAndAcol[i] ].acol() == indxCol) {
533
iParton.push_back( iColAndAcol[i] );
534
indxCol = event[ iColAndAcol[i] ].col();
535
iColAndAcol[i] = iColAndAcol.back();
536
iColAndAcol.pop_back();
540
} while (hasFound && indxCol != indxAcol && loop < loopMax);
542
// Something went wrong in colour tracing.
543
if (!hasFound || loop == loopMax) {
544
infoPtr->errorMsg("Error in HadronLevel::traceInLoop: "
545
"colour tracing failed");
554
//--------------------------------------------------------------------------
556
// Split junction-antijunction system into two, or simplify other way.
558
bool HadronLevel::splitJunctionPair(Event& event) {
560
// Construct separate index arrays for the three junction legs.
561
int identJun = (-iPartonJun[0])/10;
566
for (int i = 0; i < int(iPartonJun.size()); ++ i) {
567
if ( (-iPartonJun[i])/10 == identJun) ++leg;
568
if (leg == 0) iJunLegA.push_back( iPartonJun[i] );
569
else if (leg == 1) iJunLegB.push_back( iPartonJun[i] );
570
else iJunLegC.push_back( iPartonJun[i] );
573
// Construct separate index arrays for the three antijunction legs.
574
int identAnti = (-iPartonAntiJun[0])/10;
579
for (int i = 0; i < int(iPartonAntiJun.size()); ++ i) {
580
if ( (-iPartonAntiJun[i])/10 == identAnti) ++leg;
581
if (leg == 0) iAntiLegA.push_back( iPartonAntiJun[i] );
582
else if (leg == 1) iAntiLegB.push_back( iPartonAntiJun[i] );
583
else iAntiLegC.push_back( iPartonAntiJun[i] );
586
// Find interjunction legs, i.e. between junction and antijunction.
588
int legJun[3], legAnti[3], nGluLeg[3];
589
if (iJunLegA.back() < 0) { legJun[nMatch] = 0;
590
legAnti[nMatch] = (-iJunLegA.back())%10; ++nMatch;}
591
if (iJunLegB.back() < 0) { legJun[nMatch] = 1;
592
legAnti[nMatch] = (-iJunLegB.back())%10; ++nMatch;}
593
if (iJunLegC.back() < 0) { legJun[nMatch] = 2;
594
legAnti[nMatch] = (-iJunLegC.back())%10; ++nMatch;}
596
// Loop over interjunction legs.
597
for (int iMatch = 0; iMatch < nMatch; ++iMatch) {
598
vector<int>& iJunLeg = (legJun[iMatch] == 0) ? iJunLegA
599
: ( (legJun[iMatch] == 1) ? iJunLegB : iJunLegC );
600
vector<int>& iAntiLeg = (legAnti[iMatch] == 0) ? iAntiLegA
601
: ( (legAnti[iMatch] == 1) ? iAntiLegB : iAntiLegC );
603
// Find number of gluons on each. Do nothing for now if none.
604
nGluLeg[iMatch] = iJunLeg.size() + iAntiLeg.size() - 4;
605
if (nGluLeg[iMatch] == 0) continue;
607
// Else pick up the gluons on the interjunction leg in order.
609
for (int i = 1; i < int(iJunLeg.size()) - 1; ++i)
610
iGluLeg.push_back( iJunLeg[i] );
611
for (int i = int(iAntiLeg.size()) - 2; i > 0; --i)
612
iGluLeg.push_back( iAntiLeg[i] );
614
// Remove those gluons from the junction/antijunction leg lists.
618
// Pick a new quark at random; for simplicity no diquarks.
619
int idQ = flavSel.pickLightQ();
622
// If one gluon on leg, split it into a collinear q-qbar pair.
623
if (iGluLeg.size() == 1) {
625
// Store the new q qbar pair, sharing gluon colour and momentum.
626
colQ = event[ iGluLeg[0] ].col();
627
acolQ = event[ iGluLeg[0] ].acol();
628
Vec4 pQ = 0.5 * event[ iGluLeg[0] ].p();
629
double mQ = 0.5 * event[ iGluLeg[0] ].m();
630
int iQ = event.append( idQ, 75, iGluLeg[0], 0, 0, 0, colQ, 0, pQ, mQ );
631
int iQbar = event.append( -idQ, 75, iGluLeg[0], 0, 0, 0, 0, acolQ,
634
// Mark split gluon and update junction and antijunction legs.
635
event[ iGluLeg[0] ].statusNeg();
636
event[ iGluLeg[0] ].daughters( iQ, iQbar);
637
iJunLeg.push_back(iQ);
638
iAntiLeg.push_back(iQbar);
640
// If several gluons on the string, decide which g-g region to split up.
643
// Evaluate mass-squared for all adjacent gluon pairs.
646
for (int i = 0; i < int(iGluLeg.size()) - 1; ++i) {
647
double m2Now = 0.5 * event[ iGluLeg[i] ].p()
648
* event[ iGluLeg[i + 1] ].p();
649
m2Pair.push_back(m2Now);
653
// Pick breakup region with probability proportional to mass-squared.
654
double m2Reg = m2Sum * rndmPtr->flat();
656
do m2Reg -= m2Pair[++iReg];
657
while (m2Reg > 0. && iReg < int(iGluLeg.size()) - 1);
658
m2Reg = m2Pair[iReg];
660
// Pick breaking point of string in chosen region (symmetrically).
661
double m2Temp = min( JJSTRINGM2MAX, JJSTRINGM2FRAC * m2Reg);
665
double zTemp = zSel.zFrag( idQ, 0, m2Temp);
667
xNeg = m2Temp / (zTemp * m2Reg);
669
if (rndmPtr->flat() > 0.5) swap(xPos, xNeg);
671
// Pick up two "mother" gluons of breakup. Mark them decayed.
672
Particle& gJun = event[ iGluLeg[iReg] ];
673
Particle& gAnti = event[ iGluLeg[iReg + 1] ];
676
int dau1 = event.size();
677
gJun.daughters(dau1, dau1 + 3);
678
gAnti.daughters(dau1, dau1 + 3);
679
int mother1 = min( iGluLeg[iReg], iGluLeg[iReg + 1]);
680
int mother2 = max( iGluLeg[iReg], iGluLeg[iReg + 1]);
682
// Can keep one of old colours but need one new so unambiguous.
684
acolQ = event.nextColTag();
686
// Store copied gluons with reduced momenta.
687
int iGjun = event.append( 21, 75, mother1, mother2, 0, 0,
688
gJun.col(), gJun.acol(), (1. - 0.5 * xPos) * gJun.p(),
689
(1. - 0.5 * xPos) * gJun.m());
690
int iGanti = event.append( 21, 75, mother1, mother2, 0, 0,
691
acolQ, gAnti.acol(), (1. - 0.5 * xNeg) * gAnti.p(),
692
(1. - 0.5 * xNeg) * gAnti.m());
694
// Store the new q qbar pair with remaining momenta.
695
int iQ = event.append( idQ, 75, mother1, mother2, 0, 0,
696
colQ, 0, 0.5 * xNeg * gAnti.p(), 0.5 * xNeg * gAnti.m() );
697
int iQbar = event.append( -idQ, 75, mother1, mother2, 0, 0,
698
0, acolQ, 0.5 * xPos * gJun.p(), 0.5 * xPos * gJun.m() );
700
// Update junction and antijunction legs with gluons and quarks.
701
for (int i = 0; i < iReg; ++i)
702
iJunLeg.push_back( iGluLeg[i] );
703
iJunLeg.push_back(iGjun);
704
iJunLeg.push_back(iQ);
705
for (int i = int(iGluLeg.size()) - 1; i > iReg + 1; --i)
706
iAntiLeg.push_back( iGluLeg[i] );
707
iAntiLeg.push_back(iGanti);
708
iAntiLeg.push_back(iQbar);
711
// Update end colours for both g -> q qbar and g g -> g g q qbar.
712
event.endColJunction(identJun - 1, legJun[iMatch], colQ);
713
event.endColJunction(identAnti - 1, legAnti[iMatch], acolQ);
716
// Update list of interjunction legs after splittings above.
718
while (iMatchUp < nMatch) {
719
if (nGluLeg[iMatchUp] > 0) {
720
for (int i = iMatchUp; i < nMatch - 1; ++i) {
721
legJun[i] = legJun[i + 1];
722
legAnti[i] = legAnti[i + 1];
723
nGluLeg[i] = nGluLeg[i + 1];
728
// Should not ever have three empty interjunction legs.
730
infoPtr->errorMsg("Error in HadronLevel::splitJunctionPair: "
731
"three empty junction-junction legs");
735
// If two legs are empty, then collapse system to a single string.
737
int legJunLeft = 3 - legJun[0] - legJun[1];
738
int legAntiLeft = 3 - legAnti[0] - legAnti[1];
739
vector<int>& iJunLeg = (legJunLeft == 0) ? iJunLegA
740
: ( (legJunLeft == 1) ? iJunLegB : iJunLegC );
741
vector<int>& iAntiLeg = (legAntiLeft == 0) ? iAntiLegA
742
: ( (legAntiLeft == 1) ? iAntiLegB : iAntiLegC );
743
iPartonJun.resize(0);
744
for (int i = int(iJunLeg.size()) - 1; i > 0; --i)
745
iPartonJun.push_back( iJunLeg[i] );
746
for (int i = 1; i < int(iAntiLeg.size()); ++i)
747
iPartonJun.push_back( iAntiLeg[i] );
749
// Match up the colours where the strings are joined.
750
int iColJoin = iJunLeg[1];
751
int iAcolJoin = iAntiLeg[1];
752
event[iAcolJoin].acol( event[iColJoin].col() );
754
// Other string system empty. Remove junctions from their list. Done.
755
iPartonAntiJun.resize(0);
756
event.eraseJunction( max(identJun, identAnti) - 1);
757
event.eraseJunction( min(identJun, identAnti) - 1);
761
// If one leg is empty then, depending on string length, either
762
// (a) annihilate junction and antijunction into two simple strings, or
763
// (b) split the empty leg by borrowing energy from nearby legs.
766
// Identify the two external legs of either junction.
767
vector<int>& iJunLeg0 = (legJun[0] == 0) ? iJunLegB : iJunLegA;
768
vector<int>& iJunLeg1 = (legJun[0] == 2) ? iJunLegB : iJunLegC;
769
vector<int>& iAntiLeg0 = (legAnti[0] == 0) ? iAntiLegB : iAntiLegA;
770
vector<int>& iAntiLeg1 = (legAnti[0] == 2) ? iAntiLegB : iAntiLegC;
772
// Simplified procedure: mainly study first parton on each leg.
773
Vec4 pJunLeg0 = event[ iJunLeg0[1] ].p();
774
Vec4 pJunLeg1 = event[ iJunLeg1[1] ].p();
775
Vec4 pAntiLeg0 = event[ iAntiLeg0[1] ].p();
776
Vec4 pAntiLeg1 = event[ iAntiLeg1[1] ].p();
778
// Starting frame hopefully intermediate to two junction directions.
779
Vec4 pStart = pJunLeg0 / pJunLeg0.e() + pJunLeg1 / pJunLeg1.e()
780
+ pAntiLeg0 / pAntiLeg0.e() + pAntiLeg1 / pAntiLeg1.e();
782
// Loop over iteration to junction/antijunction rest frames (JRF/ARF).
783
RotBstMatrix MtoJRF, MtoARF;
784
Vec4 pInJRF[3], pInARF[3];
785
for (int iJun = 0; iJun < 2; ++iJun) {
786
int offset = (iJun == 0) ? 0 : 2;
788
// Iterate from system rest frame towards the junction rest frame.
789
RotBstMatrix MtoRF, Mstep;
790
MtoRF.bstback(pStart);
796
// Find rest-frame momenta on the three sides of the junction.
797
// Only consider first parton on each leg, for simplicity.
798
pInRF[0 + offset] = pJunLeg0;
799
pInRF[1 + offset] = pJunLeg1;
800
pInRF[2 - offset] = pAntiLeg0;
801
pInRF[3 - offset] = pAntiLeg1;
802
for (int i = 0; i < 4; ++i) pInRF[i].rotbst(MtoRF);
804
// For third side add both legs beyond other junction, weighted.
805
double wt2 = 1. - exp( -pInRF[2].e() / eNormJunction);
806
double wt3 = 1. - exp( -pInRF[3].e() / eNormJunction);
807
pInRF[2] = wt2 * pInRF[2] + wt3 * pInRF[3];
809
// Find new junction rest frame from the set of momenta.
810
Mstep = stringFrag.junctionRestFrame( pInRF[0], pInRF[1], pInRF[2]);
811
MtoRF.rotbst( Mstep );
812
} while (iter < 3 || (Mstep.deviation() > CONVJNREST
813
&& iter < NTRYJNREST) );
815
// Store final boost and rest-frame (weighted) momenta.
818
for (int i = 0; i < 3; ++i) pInJRF[i] = pInRF[i];
821
for (int i = 0; i < 3; ++i) pInARF[i] = pInRF[i];
825
// Opposite operations: boost from JRF/ARF to original system.
826
RotBstMatrix MfromJRF = MtoJRF;
828
RotBstMatrix MfromARF = MtoARF;
831
// Velocity vectors of junctions and momentum of legs in lab frame.
832
Vec4 vJun(0., 0., 0., 1.);
833
vJun.rotbst(MfromJRF);
834
Vec4 vAnti(0., 0., 0., 1.);
835
vAnti.rotbst(MfromARF);
836
Vec4 pLabJ[3], pLabA[3];
837
for (int i = 0; i < 3; ++i) {
838
pLabJ[i] = pInJRF[i];
839
pLabJ[i].rotbst(MfromJRF);
840
pLabA[i] = pInARF[i];
841
pLabA[i].rotbst(MfromARF);
844
// Calculate Lambda-measure length of three possible topologies.
845
double vJvA = vJun * vAnti;
846
double vJvAe2y = vJvA + sqrt(vJvA*vJvA - 1.);
847
double LambdaJA = (2. * pInJRF[0].e()) * (2. * pInJRF[1].e())
848
* (2. * pInARF[0].e()) * (2. * pInARF[1].e()) * vJvAe2y;
849
double Lambda00 = (2. * pLabJ[0] * pLabA[0])
850
* (2. * pLabJ[1] * pLabA[1]);
851
double Lambda01 = (2. * pLabJ[0] * pLabA[1])
852
* (2. * pLabJ[1] * pLabA[0]);
854
// Case when either topology without junctions is the shorter one.
855
if (LambdaJA > min( Lambda00, Lambda01)) {
856
vector<int>& iAntiMatch0 = (Lambda00 < Lambda01)
857
? iAntiLeg0 : iAntiLeg1;
858
vector<int>& iAntiMatch1 = (Lambda00 < Lambda01)
859
? iAntiLeg1 : iAntiLeg0;
861
// Define two quark-antiquark strings.
862
iPartonJun.resize(0);
863
for (int i = int(iJunLeg0.size()) - 1; i > 0; --i)
864
iPartonJun.push_back( iJunLeg0[i] );
865
for (int i = 1; i < int(iAntiMatch0.size()); ++i)
866
iPartonJun.push_back( iAntiMatch0[i] );
867
iPartonAntiJun.resize(0);
868
for (int i = int(iJunLeg1.size()) - 1; i > 0; --i)
869
iPartonAntiJun.push_back( iJunLeg1[i] );
870
for (int i = 1; i < int(iAntiMatch1.size()); ++i)
871
iPartonAntiJun.push_back( iAntiMatch1[i] );
873
// Match up the colours where the strings are joined.
874
int iColJoin = iJunLeg0[1];
875
int iAcolJoin = iAntiMatch0[1];
876
event[iAcolJoin].acol( event[iColJoin].col() );
877
iColJoin = iJunLeg1[1];
878
iAcolJoin = iAntiMatch1[1];
879
event[iAcolJoin].acol( event[iColJoin].col() );
881
// Remove junctions from their list. Done.
882
event.eraseJunction( max(identJun, identAnti) - 1);
883
event.eraseJunction( min(identJun, identAnti) - 1);
887
// Case where junction and antijunction to be separated.
888
// Shuffle (p+/p-) momentum of order <mThad> between systems,
889
// times 2/3 for 120 degree in JRF, times 1/2 for two legs,
890
// but not more than half of what nearest parton carries.
891
double eShift = MTHAD / (3. * sqrt(vJvAe2y));
892
double fracJ0 = min(0.5, eShift / pInJRF[0].e());
893
double fracJ1 = min(0.5, eShift / pInJRF[0].e());
894
Vec4 pFromJun = fracJ0 * pJunLeg0 + fracJ1 * pJunLeg1;
895
double fracA0 = min(0.5, eShift / pInARF[0].e());
896
double fracA1 = min(0.5, eShift / pInARF[0].e());
897
Vec4 pFromAnti = fracA0 * pAntiLeg0 + fracA1 * pAntiLeg1;
899
// Copy partons with scaled-down momenta and update legs.
900
int iNew = event.copy(iJunLeg0[1], 76);
901
event[iNew].rescale5(1. - fracJ0);
903
iNew = event.copy(iJunLeg1[1], 76);
904
event[iNew].rescale5(1. - fracJ1);
906
iNew = event.copy(iAntiLeg0[1], 76);
907
event[iNew].rescale5(1. - fracA0);
909
iNew = event.copy(iAntiLeg1[1], 76);
910
event[iNew].rescale5(1. - fracA1);
913
// Pick a new quark at random; for simplicity no diquarks.
914
int idQ = flavSel.pickLightQ();
916
// Update junction colours for new quark and antiquark.
917
int colQ = event.nextColTag();
918
int acolQ = event.nextColTag();
919
event.endColJunction(identJun - 1, legJun[0], colQ);
920
event.endColJunction(identAnti - 1, legAnti[0], acolQ);
922
// Store new quark and antiquark with momentum from other junction.
923
int mother1 = min(iJunLeg0[1], iJunLeg1[1]);
924
int mother2 = max(iJunLeg0[1], iJunLeg1[1]);
925
int iNewJ = event.append( idQ, 76, mother1, mother2, 0, 0,
926
colQ, 0, pFromAnti, pFromAnti.mCalc() );
927
mother1 = min(iAntiLeg0[1], iAntiLeg1[1]);
928
mother2 = max(iAntiLeg0[1], iAntiLeg1[1]);
929
int iNewA = event.append( -idQ, 76, mother1, mother2, 0, 0,
930
0, acolQ, pFromJun, pFromJun.mCalc() );
932
// Bookkeep new quark and antiquark on third legs.
933
if (legJun[0] == 0) iJunLegA[1] = iNewJ;
934
else if (legJun[0] == 1) iJunLegB[1] = iNewJ;
935
else iJunLegC[1] = iNewJ;
936
if (legAnti[0] == 0) iAntiLegA[1] = iNewA;
937
else if (legAnti[0] == 1) iAntiLegB[1] = iNewA;
938
else iAntiLegC[1] = iNewA;
940
// Done with splitting junction from antijunction.
943
// Put together new junction parton list.
944
iPartonJun.resize(0);
945
for (int i = 0; i < int(iJunLegA.size()); ++i)
946
iPartonJun.push_back( iJunLegA[i] );
947
for (int i = 0; i < int(iJunLegB.size()); ++i)
948
iPartonJun.push_back( iJunLegB[i] );
949
for (int i = 0; i < int(iJunLegC.size()); ++i)
950
iPartonJun.push_back( iJunLegC[i] );
952
// Put together new antijunction parton list.
953
iPartonAntiJun.resize(0);
954
for (int i = 0; i < int(iAntiLegA.size()); ++i)
955
iPartonAntiJun.push_back( iAntiLegA[i] );
956
for (int i = 0; i < int(iAntiLegB.size()); ++i)
957
iPartonAntiJun.push_back( iAntiLegB[i] );
958
for (int i = 0; i < int(iAntiLegC.size()); ++i)
959
iPartonAntiJun.push_back( iAntiLegC[i] );
961
// Now the two junction systems are separated and can be stored.
966
//==========================================================================
968
} // end namespace Pythia8
added
removed
src/HadronLevel.cc
