]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8130/src/SpaceShower.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SpaceShower.cxx
CommitLineData
5ad4eb21 1// SpaceShower.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 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// SpaceShower class.
8
9#include "SpaceShower.h"
10
11namespace Pythia8 {
12
13//**************************************************************************
14
15// The SpaceShower 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// Leftover companion can give PDF > 0 at small Q2 where other PDF's = 0,
23// and then one can end in infinite loop of impossible kinematics.
24const int SpaceShower::MAXLOOPTINYPDF = 10;
25
26// Switch to alternative (but equivalent) backwards evolution for
27// g -> Q Qbar (Q = c or b) when below QTHRESHOLD * mQ2.
28const double SpaceShower::CTHRESHOLD = 2.0;
29const double SpaceShower::BTHRESHOLD = 2.0;
30
31// Renew evaluation of PDF's when the pT2 step is bigger than this
32// (in addition to initial scale and c and b thresholds.)
33const double SpaceShower::EVALPDFSTEP = 0.1;
34
35// Lower limit on PDF value in order to avoid division by zero.
36const double SpaceShower::TINYPDF = 1e-10;
37
38// Lower limit on estimated evolution rate, below which stop.
39const double SpaceShower::TINYKERNELPDF = 1e-6;
40
41// Lower limit on pT2, below which branching is rejected.
42const double SpaceShower::TINYPT2 = 0.25e-6;
43
44// No attempt to do backwards evolution of a heavy (c or b) quark
45// if evolution starts at a scale pT2 < HEAVYPT2EVOL * mQ2.
46const double SpaceShower::HEAVYPT2EVOL = 1.1;
47
48// No attempt to do backwards evolution of a heavy (c or b) quark
49// if evolution starts at a x > HEAVYXEVOL * x_max, where
50// x_max is the largest possible x value for a g -> Q Qbar branching.
51const double SpaceShower::HEAVYXEVOL = 0.9;
52
53// When backwards evolution Q -> g + Q creates a heavy quark Q,
54// an earlier branching g -> Q + Qbar will restrict kinematics
55// to M_{Q Qbar}^2 > EXTRASPACEQ * 4 m_Q^2. (Smarter to be found??)
56const double SpaceShower::EXTRASPACEQ = 2.0;
57
58// Never pick pT so low that alphaS is evaluated too close to Lambda_3.
59const double SpaceShower::LAMBDA3MARGIN = 1.1;
60
61// Cutoff for f_e^e at x < 1 - 10^{-10} to be used in z selection.
62// Note: the x_min quantity come from 1 - x_max.
63const double SpaceShower::LEPTONXMIN = 1e-10;
64const double SpaceShower::LEPTONXMAX = 1. - 1e-10;
65
66// Stop l -> l gamma evolution slightly above m2l.
67const double SpaceShower::LEPTONPT2MIN = 1.2;
68
69// Enhancement of l -> l gamma trial rate to compensate imperfect modelling.
70const double SpaceShower::LEPTONFUDGE = 10.;
71
72//*********
73
74// Initialize alphaStrong, alphaEM and related pTmin parameters.
75
76void SpaceShower::init( BeamParticle* beamAPtrIn,
77 BeamParticle* beamBPtrIn) {
78
79 // Store input pointers for future use.
80 beamAPtr = beamAPtrIn;
81 beamBPtr = beamBPtrIn;
82
83 // Main flags to switch on and off branchings.
84 doQCDshower = Settings::flag("SpaceShower:QCDshower");
85 doQEDshowerByQ = Settings::flag("SpaceShower:QEDshowerByQ");
86 doQEDshowerByL = Settings::flag("SpaceShower:QEDshowerByL");
87
88 // Matching in pT of hard interaction to shower evolution.
89 pTmaxMatch = Settings::mode("SpaceShower:pTmaxMatch");
90 pTdampMatch = Settings::mode("SpaceShower:pTdampMatch");
91 pTmaxFudge = Settings::parm("SpaceShower:pTmaxFudge");
92 pTdampFudge = Settings::parm("SpaceShower:pTdampFudge");
93
94 // Optionally force emissions to be ordered in rapidity/angle.
95 doRapidityOrder = Settings::flag("SpaceShower:rapidityOrder");
96
97 // Charm, bottom and lepton mass thresholds.
98 mc = ParticleDataTable::m0(4);
99 mb = ParticleDataTable::m0(5);
100 m2c = pow2(mc);
101 m2b = pow2(mb);
102
103 // Parameters of alphaStrong generation.
104 alphaSvalue = Settings::parm("SpaceShower:alphaSvalue");
105 alphaSorder = Settings::mode("SpaceShower:alphaSorder");
106 alphaS2pi = 0.5 * alphaSvalue / M_PI;
107
108 // Initialize alpha_strong generation.
109 alphaS.init( alphaSvalue, alphaSorder);
110
111 // Lambda for 5, 4 and 3 flavours.
112 Lambda5flav = alphaS.Lambda5();
113 Lambda4flav = alphaS.Lambda4();
114 Lambda3flav = alphaS.Lambda3();
115 Lambda5flav2 = pow2(Lambda5flav);
116 Lambda4flav2 = pow2(Lambda4flav);
117 Lambda3flav2 = pow2(Lambda3flav);
118
119 // Regularization of QCD evolution for pT -> 0. Can be taken
120 // same as for multiple interactions, or be set separately.
121 useSamePTasMI = Settings::flag("SpaceShower:samePTasMI");
122 if (useSamePTasMI) {
123 pT0Ref = Settings::parm("MultipleInteractions:pT0Ref");
124 ecmRef = Settings::parm("MultipleInteractions:ecmRef");
125 ecmPow = Settings::parm("MultipleInteractions:ecmPow");
126 pTmin = Settings::parm("MultipleInteractions:pTmin");
127 } else {
128 pT0Ref = Settings::parm("SpaceShower:pT0Ref");
129 ecmRef = Settings::parm("SpaceShower:ecmRef");
130 ecmPow = Settings::parm("SpaceShower:ecmPow");
131 pTmin = Settings::parm("SpaceShower:pTmin");
132 }
133
134 // Calculate nominal invariant mass of events. Set current pT0 scale.
135 sCM = m2( beamAPtr->p(), beamBPtr->p());
136 eCM = sqrt(sCM);
137 pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
138
139 // Restrict pTmin to ensure that alpha_s(pTmin^2 + pT_0^2) does not blow up.
140 pTmin = max( pTmin, sqrtpos(pow2(LAMBDA3MARGIN) * Lambda3flav2 - pT0*pT0) );
141
142 // Parameters of alphaEM generation.
143 alphaEMorder = Settings::mode("SpaceShower:alphaEMorder");
144
145 // Initialize alphaEM generation.
146 alphaEM.init( alphaEMorder);
147
148 // Parameters of QED evolution.
149 pTminChgQ = Settings::parm("SpaceShower:pTminchgQ");
150 pTminChgL = Settings::parm("SpaceShower:pTminchgL");
151
152 // Derived parameters of QCD evolution.
153 pT20 = pow2(pT0);
154 pT2min = pow2(pTmin);
155 pT2minChgQ = pow2(pTminChgQ);
156 pT2minChgL = pow2(pTminChgL);
157
158 // Various other parameters.
159 doMEcorrections = Settings::flag("SpaceShower:MEcorrections");
160 doPhiPolAsym = Settings::flag("SpaceShower:phiPolAsym");
161 nQuarkIn = Settings::mode("SpaceShower:nQuarkIn");
162
163}
164
165//*********
166
167// Find whether to limit maximum scale of emissions.
168
169bool SpaceShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) {
170
171 // Find whether to limit pT. Begin by user-set cases.
172 bool dopTlimit = false;
173 if (pTmaxMatch == 1) dopTlimit = true;
174 else if (pTmaxMatch == 2) dopTlimit = false;
175
176 // Look if any quark (u, d, s, c, b), gluon or photon in final state.
177 else {
178 for (int i = 5; i < event.size(); ++i)
179 if (event[i].status() != -21) {
180 int idAbs = event[i].idAbs();
181 if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit = true;
182 }
183 }
184
185 // Dampening at factorization or renormalization scale.
186 dopTdamp = false;
187 pT2damp = 0.;
188 if ( !dopTlimit && (pTdampMatch == 1 || pTdampMatch == 2) ) {
189 dopTdamp = true;
190 pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren);
191 }
192
193 // Done.
194 return dopTlimit;
195
196}
197
198//*********
199
200// Prepare system for evolution; identify ME.
201// Routine may be called after multiple interactions, for a new subystem.
202
203void SpaceShower::prepare( int iSys, Event& event, bool limitPTmaxIn) {
204
205 // Find positions of incoming colliding partons.
206 int in1 = event.getInSystem( iSys, 0);
207 int in2 = event.getInSystem( iSys, 1);
208
209 // Reset dipole-ends list for first interaction.
210 if (iSys == 0) dipEnd.resize(0);
211
212 // Find matrix element corrections for system.
213 int MEtype = findMEtype( iSys, event);
214
215 // Maximum pT scale for dipole ends.
216 double pTmax1 = (limitPTmaxIn) ? event[in1].scale() : eCM;
217 double pTmax2 = (limitPTmaxIn) ? event[in2].scale() : eCM;
218 if (iSys == 0 && limitPTmaxIn) {
219 pTmax1 *= pTmaxFudge;
220 pTmax2 *= pTmaxFudge;
221 }
222
223 // Find dipole ends for QCD radiation.
224 if (doQCDshower) {
225 int colType1 = event[in1].colType();
226 dipEnd.push_back( SpaceDipoleEnd( iSys, 1, in1, in2, pTmax1,
227 colType1, 0, MEtype) );
228 int colType2 = event[in2].colType();
229 dipEnd.push_back( SpaceDipoleEnd( iSys, 2, in2, in1, pTmax2,
230 colType2, 0, MEtype) );
231 }
232
233 // Find dipole ends for QED radiation.
234 if (doQEDshowerByQ || doQEDshowerByL) {
235 int chgType1 = ( (event[in1].isQuark() && doQEDshowerByQ)
236 || (event[in1].isLepton() && doQEDshowerByL) )
237 ? event[in1].chargeType() : 0;
238 dipEnd.push_back( SpaceDipoleEnd( iSys, -1, in1, in2, pTmax1,
239 0, chgType1, MEtype) );
240 int chgType2 = ( (event[in2].isQuark() && doQEDshowerByQ)
241 || (event[in2].isLepton() && doQEDshowerByL) )
242 ? event[in2].chargeType() : 0;
243 dipEnd.push_back( SpaceDipoleEnd( iSys, -2, in2, in1, pTmax2,
244 0, chgType2, MEtype) );
245 }
246
247}
248
249//*********
250
251// Select next pT in downwards evolution of the existing dipoles.
252
253double SpaceShower::pTnext( Event& , double pTbegAll, double pTendAll,
254 int nRadIn) {
255
256 // Current cm energy, in case it varies between events.
257 sCM = m2( beamAPtr->p(), beamBPtr->p());
258
259 // Starting values: no radiating dipole found.
260 nRad = nRadIn;
261 double pT2sel = pow2(pTendAll);
262 iDipSel = 0;
263 iSysSel = 0;
264 dipEndSel = 0;
265
266 // Loop over all possible dipole ends.
267 for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
268 iDipNow = iDipEnd;
269 dipEndNow = &dipEnd[iDipEnd];
270 iSysNow = dipEndNow->system;
271 dipEndNow->pT2 = 0.;
272
273 // Check whether dipole end should be allowed to shower.
274 double pT2begDip = pow2( min( pTbegAll, dipEndNow->pTmax ));
275 if (pT2begDip > pT2sel
276 && ( dipEndNow->colType != 0 || dipEndNow->chgType != 0 ) ) {
277 double pT2endDip = 0.;
278
279 // Determine lower cut for evolution, for QCD or QED (q or l).
280 if (dipEndNow->colType != 0) pT2endDip = max( pT2sel, pT2min );
281 else if (abs(dipEndNow->chgType) != 3) pT2endDip
282 = max( pT2sel, pT2minChgQ );
283 else pT2endDip = max( pT2sel, pT2minChgL );
284
285 // Find properties of dipole and radiating dipole end.
286 bool ordered = ( abs(dipEndNow->side) == 1 );
287 BeamParticle& beamNow = (ordered) ? *beamAPtr : *beamBPtr;
288 BeamParticle& beamRec = (ordered) ? *beamBPtr : *beamAPtr;
289 iNow = beamNow[iSysNow].iPos();
290 iRec = beamRec[iSysNow].iPos();
291 idDaughter = beamNow[iSysNow].id();
292 xDaughter = beamNow[iSysNow].x();
293 x1Now = (ordered) ? xDaughter : beamRec[iSysNow].x();
294 x2Now = (ordered) ? beamRec[iSysNow].x() : xDaughter;
295 m2Dip = x1Now * x2Now * sCM;
296
297 // Now do evolution in pT2, for QCD or QED
298 if (pT2begDip > pT2endDip) {
299 if (dipEndNow->colType != 0) pT2nextQCD( pT2begDip, pT2endDip);
300 else pT2nextQED( pT2begDip, pT2endDip);
301 }
302
303 // Update if found larger pT than current maximum.
304 if (dipEndNow->pT2 > pT2sel) {
305 pT2sel = dipEndNow->pT2;
306 iDipSel = iDipNow;
307 iSysSel = iSysNow;
308 dipEndSel = dipEndNow;
309 }
310
311 // End loop over dipole ends.
312 }
313 }
314
315 // Return nonvanishing value if found pT is bigger than already found.
316 return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
317}
318
319//*********
320
321// Evolve a QCD dipole end.
322
323void SpaceShower::pT2nextQCD( double pT2begDip, double pT2endDip) {
324
325 // Some properties and kinematical starting values.
326 BeamParticle& beam = ( abs(dipEndNow->side) == 1 )
327 ? *beamAPtr : *beamBPtr;
328 bool isGluon = (idDaughter == 21);
329 bool isValence = beam[iSysNow].isValence();
330 int MEtype = dipEndNow->MEtype;
331 double pT2 = pT2begDip;
332 double xMaxAbs = beam.xMax(iSysNow);
333 double zMinAbs = xDaughter / xMaxAbs;
334
335 // Starting values for handling of massive quarks (c/b), if any.
336 double idMassive = 0;
337 if ( abs(idDaughter) == 4 ) idMassive = 4;
338 if ( abs(idDaughter) == 5 ) idMassive = 5;
339 bool isMassive = (idMassive > 0);
340 double m2Massive = 0.;
341 double mRatio = 0.;
342 double zMaxMassive = 1.;
343 double m2Threshold = pT2;
344
345 // Evolution below scale of massive quark or at large x is impossible.
346 if (isMassive) {
347 m2Massive = (idMassive == 4) ? m2c : m2b;
348 if (pT2 < HEAVYPT2EVOL * m2Massive) return;
349 mRatio = sqrt( m2Massive / m2Dip );
350 zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
351 if (xDaughter > HEAVYXEVOL * zMaxMassive * xMaxAbs) return;
352
353 // Find threshold scale below which only g -> Q + Qbar will be allowed.
354 m2Threshold = (idMassive == 4) ? min( pT2, CTHRESHOLD * m2c)
355 : min( pT2, BTHRESHOLD * m2b);
356 }
357
358 // Variables used inside evolution loop. (Mainly dummy starting values.)
359 int nFlavour = 3;
360 double b0 = 4.5;
361 double Lambda2 = Lambda3flav2;
362 double pT2minNow = pT2endDip;
363 int idMother = 0;
364 int idSister = 0;
365 double z = 0.;
366 double zMaxAbs = 0.;
367 double zRootMax = 0.;
368 double zRootMin = 0.;
369 double g2gInt = 0.;
370 double q2gInt = 0.;
371 double q2qInt = 0.;
372 double g2qInt = 0.;
373 double g2Qenhance = 0.;
374 double xPDFdaughter = 0.;
375 double xPDFmother[21] = {0.};
376 double xPDFgMother = 0.;
377 double xPDFmotherSum = 0.;
378 double kernelPDF = 0.;
379 double xMother = 0.;
380 double wt = 0.;
381 double Q2 = 0.;
382 double mSister = 0.;
383 double m2Sister = 0.;
384 double pT2corr = 0.;
385 double phi = 0.;
386 double pT2PDF = pT2;
387 bool needNewPDF = true;
388
389 // Begin evolution loop towards smaller pT values.
390 int loopTinyPDFdau = 0;
391 bool hasTinyPDFdau = false;
392 do {
393 wt = 0.;
394
395 // Bad sign if repeated looping with small daughter PDF, so fail.
396 // (Example: if all PDF's = 0 below Q_0, except for c/b companion.)
397 if (hasTinyPDFdau) ++loopTinyPDFdau;
398 if (loopTinyPDFdau > MAXLOOPTINYPDF) {
399 infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
400 "small daughter PDF");
401 return;
402 }
403
404 // Initialize integrals of splitting kernels and evaluate parton
405 // densities at the beginning. Reinitialize after long evolution
406 // in pT2 or when crossing c and b flavour thresholds.
407 if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
408 pT2PDF = pT2;
409 hasTinyPDFdau = false;
410
411 // Determine overestimated z range; switch at c and b masses.
412 if (pT2 > m2b) {
413 nFlavour = 5;
414 pT2minNow = m2b;
415 b0 = 23./6.;
416 Lambda2 = Lambda5flav2;
417 } else if (pT2 > m2c) {
418 nFlavour = 4;
419 pT2minNow = m2c;
420 b0 = 25./6.;
421 Lambda2 = Lambda4flav2;
422 } else {
423 nFlavour = 3;
424 pT2minNow = pT2endDip;
425 b0 = 27./6.;
426 Lambda2 = Lambda3flav2;
427 }
428 zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
429 ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
430 if (isMassive) zMaxAbs = min( zMaxAbs, zMaxMassive);
431
432 // Go to another z range with lower mass scale if current is closed.
433 if (zMinAbs > zMaxAbs) {
434 if (nFlavour == 3 || (idMassive == 4 && nFlavour == 4)
435 || idMassive == 5) return;
436 pT2 = (nFlavour == 4) ? m2c : m2b;
437 continue;
438 }
439
440 // Parton density of daughter at current scale.
441 xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pT2);
442 if (xPDFdaughter < TINYPDF) {
443 xPDFdaughter = TINYPDF;
444 hasTinyPDFdau = true;
445 }
446
447 // Integrals of splitting kernels for gluons: g -> g, q -> g.
448 if (isGluon) {
449 g2gInt = 6. * log(zMaxAbs * (1.-zMinAbs)
450 / (zMinAbs * (1.-zMaxAbs)));
451 if (doMEcorrections) g2gInt *= calcMEmax(MEtype, 21, 21);
452 q2gInt = (16./3.) * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
453 if (doMEcorrections) q2gInt *= calcMEmax(MEtype, 1, 21);
454
455 // Parton density of potential quark mothers to a g.
456 xPDFmotherSum = 0.;
457 for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
458 if (i == 0) {
459 xPDFmother[10] = 0.;
460 } else {
461 xPDFmother[i+10] = beam.xfISR(iSysNow, i, xDaughter, pT2);
462 xPDFmotherSum += xPDFmother[i+10];
463 }
464 }
465
466 // Total QCD evolution coefficient for a gluon.
467 kernelPDF = g2gInt + q2gInt * xPDFmotherSum / xPDFdaughter;
468
469 // For valence quark only need consider q -> q g branchings.
470 // Introduce an extra factor sqrt(z) to smooth bumps.
471 } else if (isValence) {
472 zRootMin = (1. + sqrt(zMinAbs)) / (1. - sqrt(zMinAbs));
473 zRootMax = (1. + sqrt(zMaxAbs)) / (1. - sqrt(zMaxAbs));
474 q2qInt = (8./3.) * log( zRootMax / zRootMin );
475 if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
476 g2qInt = 0.;
477 kernelPDF = q2qInt;
478
479 // Integrals of splitting kernels for quarks: q -> q, g -> q.
480 } else {
481 q2qInt = (8./3.) * log( (1. - zMinAbs) / (1. - zMaxAbs) );
482 if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
483 g2qInt = 0.5 * (zMaxAbs - zMinAbs);
484 if (doMEcorrections) g2qInt *= calcMEmax(MEtype, 21, 1);
485
486 // Increase estimated upper weight for g -> Q + Qbar.
487 if (isMassive) {
488 double m2log = log( m2Massive / Lambda2);
489 g2Qenhance = log( log(pT2/Lambda2) / m2log )
490 / log( log(m2Threshold/Lambda2) / m2log );
491 g2qInt *= g2Qenhance;
492 }
493
494 // Parton density of a potential gluon mother to a q.
495 xPDFgMother = beam.xfISR(iSysNow, 21, xDaughter, pT2);
496
497 // Total QCD evolution coefficient for a quark.
498 kernelPDF = q2qInt + g2qInt * xPDFgMother / xPDFdaughter;
499 }
500
501 // End evaluation of splitting kernels and parton densities.
502 needNewPDF = false;
503 }
504 if (kernelPDF < TINYKERNELPDF) return;
505
506 // Pick pT2 (in overestimated z range), for one of three different cases.
507 // Assume form alphas(pT0^2 + pT^2) * dpT^2/(pT0^2 + pT^2).
508 double Q2alphaS;
509
510 // Fixed alpha_strong.
511 if (alphaSorder == 0) {
512 pT2 = (pT2 + pT20) * pow( Rndm::flat(),
513 1. / (alphaS2pi * kernelPDF)) - pT20;
514
515 // First-order alpha_strong.
516 } else if (alphaSorder == 1) {
517 pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
518 pow(Rndm::flat(), b0 / kernelPDF) ) - pT20;
519
520 // For second order reject by second term in alpha_strong expression.
521 } else {
522 do {
523 pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
524 pow(Rndm::flat(), b0 / kernelPDF) ) - pT20;
525 Q2alphaS = max(pT2 + pT20, pow2(LAMBDA3MARGIN) * Lambda3flav2);
526 } while (alphaS.alphaS2OrdCorr(Q2alphaS) < Rndm::flat()
527 && pT2 > pT2minNow);
528 }
529
530 // Check for pT2 values that prompt special action.
531
532 // If fallen into b threshold region, force g -> b + bbar.
533 if (idMassive == 5 && pT2 < m2Threshold) {
534 pT2nearQCDthreshold( beam, m2Massive, m2Threshold, zMinAbs,
535 zMaxMassive );
536 return;
537
538 // If crossed b threshold, continue evolution from this threshold.
539 } else if (nFlavour == 5 && pT2 < m2b) {
540 needNewPDF = true;
541 pT2 = m2b;
542 continue;
543
544 // If fallen into c threshold region, force g -> c + cbar.
545 } else if (idMassive == 4 && pT2 < m2Threshold) {
546 pT2nearQCDthreshold( beam, m2Massive, m2Threshold, zMinAbs,
547 zMaxMassive );
548 return;
549
550 // If crossed c threshold, continue evolution from this threshold.
551 } else if (nFlavour == 4 && pT2 < m2c) {
552 needNewPDF = true;
553 pT2 = m2c;
554 continue;
555
556 // Abort evolution if below cutoff scale, or below another branching.
557 } else if (pT2 < pT2endDip) return;
558
559 // Select z value of branching to g, and corrective weight.
560 if (isGluon) {
561 // g -> g (+ g).
562 if (Rndm::flat() * kernelPDF < g2gInt) {
563 idMother = 21;
564 idSister = 21;
565 z = 1. / ( 1. + ((1. - zMinAbs) / zMinAbs) * pow( (zMinAbs *
566 (1. - zMaxAbs)) / (zMaxAbs * (1. - zMinAbs)), Rndm::flat() ) );
567 wt = pow2( 1. - z * (1. - z));
568 } else {
569 // q -> g (+ q): also select flavour.
570 double temp = xPDFmotherSum * Rndm::flat();
571 idMother = -nQuarkIn - 1;
572 do { temp -= xPDFmother[(++idMother) + 10]; }
573 while (temp > 0. && idMother < nQuarkIn);
574 idSister = idMother;
575 z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + Rndm::flat()
576 * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
577 wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z)
578 * xPDFdaughter / xPDFmother[idMother + 10];
579 }
580
581 // Select z value of branching to q, and corrective weight.
582 // Include massive kernel corrections for c and b quarks.
583 } else {
584 // q -> q (+ g).
585 if (isValence || Rndm::flat() * kernelPDF < q2qInt) {
586 idMother = idDaughter;
587 idSister = 21;
588 // Valence more peaked at large z.
589 if (isValence) {
590 double zTmp = zRootMin * pow(zRootMax / zRootMin, Rndm::flat() );
591 z = pow2( (1. - zTmp) / (1. + zTmp) );
592 } else {
593 z = 1. - (1. - zMinAbs) * pow( (1. - zMaxAbs) / (1. - zMinAbs),
594 Rndm::flat() );
595 }
596 if (!isMassive) {
597 wt = 0.5 * (1. + pow2(z));
598 } else {
599 wt = 0.5 * (1. + pow2(z) - z * pow2(1.-z) * m2Massive / pT2);
600 }
601 if (isValence) wt *= sqrt(z);
602 // g -> q (+ qbar).
603 } else {
604 idMother = 21;
605 idSister = - idDaughter;
606 z = zMinAbs + Rndm::flat() * (zMaxAbs - zMinAbs);
607 if (!isMassive) {
608 wt = (pow2(z) + pow2(1.-z)) * xPDFdaughter / xPDFgMother ;
609 } else {
610 wt = (pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2)
611 * xPDFdaughter / (xPDFgMother * g2Qenhance) ;
612 }
613 }
614 }
615
616 // Derive Q2 and x of mother from pT2 and z.
617 Q2 = pT2 / (1.- z);
618 xMother = xDaughter / z;
619
620 // Forbidden emission if outside allowed z range for given pT2.
621 mSister = ParticleDataTable::m0(idSister);
622 m2Sister = pow2(mSister);
623 pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
624 if(pT2corr < TINYPT2) { wt = 0.; continue; }
625
626 // Optionally veto emissions not ordered in rapidity (= angle).
627 if ( doRapidityOrder && dipEndNow->nBranch > 0
628 && pT2 > pow2( (1. - z) / (z * (1. - dipEndNow->zOld)) )
629 * dipEndNow->pT2Old ) { wt = 0.; continue; }
630
631 // If creating heavy quark by Q -> g + Q then next need g -> Q + Qbar.
632 // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
633 if ( isGluon && ( abs(idMother) == 4 || abs(idMother) == 5 )) {
634 double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
635 double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
636 if(pT2QQcorr < TINYPT2) { wt = 0.; continue; }
637 }
638
639 // Select phi angle of branching at random.
640 phi = 2. * M_PI * Rndm::flat();
641
642 // Evaluation of ME correction.
643 if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
644 m2Dip, z, Q2) / calcMEmax(MEtype, idMother, idDaughter);
645
646 // Optional dampening of large pT values in first radiation.
647 if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
648 wt *= pT2damp / (pT2 + pT2damp);
649
650 // Evaluation of new daughter and mother PDF's.
651 double xPDFdaughterNew = max ( TINYPDF,
652 beam.xfISR(iSysNow, idDaughter, xDaughter, pT2) );
653 double xPDFmotherNew = beam.xfISR(iSysNow, idMother, xMother, pT2);
654 wt *= xPDFmotherNew / xPDFdaughterNew;
655
656 // Check that valence step does not cause problem.
657 if (wt > 1.) infoPtr->errorMsg("Warning in SpaceShower::"
658 "pT2nextQCD: weight above unity");
659
660 // Iterate until acceptable pT (or have fallen below pTmin).
661 } while (wt < Rndm::flat()) ;
662
663 // Save values for (so far) acceptable branching.
664 dipEndNow->store( idDaughter,idMother, idSister, x1Now, x2Now, m2Dip,
665 pT2, z, Q2, mSister, m2Sister, pT2corr, phi);
666
667}
668
669//*********
670
671// Evolve a QCD dipole end near threshold, with g -> Q + Qbar enforced.
672// Note: No explicit Sudakov factor formalism here. Instead use that
673// df_Q(x, pT2) = (alpha_s/2pi) * (dT2/pT2) * ((gluon) * (splitting)).
674// This implies that effects of Q -> Q + g are neglected in this range.
675
676void SpaceShower::pT2nearQCDthreshold( BeamParticle& beam,
677 double m2Massive, double m2Threshold, double zMinAbs,
678 double zMaxMassive) {
679
680 // Initial values, to be used in kinematics and weighting.
681 double Lambda2 = (abs(idDaughter) == 4) ? Lambda4flav2 : Lambda5flav2;
682 double logM2Lambda2 = log( m2Massive / Lambda2 );
683 double xPDFmotherOld = beam.xfISR(iSysNow, 21, xDaughter, m2Threshold);
684
685 // Variables used inside evolution loop. (Mainly dummy start values.)
686 int loop = 0;
687 double wt = 0.;
688 double pT2 = 0.;
689 double z = 0.;
690 double Q2 = 0.;
691 double pT2corr = 0.;
692
693 // Begin loop over tries to find acceptable g -> Q + Qbar branching.
694 do {
695 wt = 0.;
696
697 // Check that not caught in infinite loop with impossible kinematics.
698 if (++loop > 100) {
699 infoPtr->errorMsg("Error in SpaceShower::pT2nearQCDthreshold: "
700 "stuck in loop");
701 return;
702 }
703
704 // Pick dpT2/pT2 in range [m2Massive,thresholdRatio * m2Massive].
705 pT2 = m2Massive * pow( m2Threshold / m2Massive, Rndm::flat() );
706
707 // Pick z flat in allowed range.
708 z = zMinAbs + Rndm::flat() * (zMaxMassive - zMinAbs);
709
710 // Check that kinematically possible choice.
711 Q2 = pT2 / (1.-z) - m2Massive;
712 pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Massive) / m2Dip;
713 if(pT2corr < TINYPT2) continue;
714
715 // Correction factor for running alpha_s. ??
716 wt = logM2Lambda2 / log( pT2 / Lambda2 );
717
718 // Correction factor for splitting kernel.
719 wt *= pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2;
720
721 // Correction factor for gluon density.
722 double xPDFmotherNew = beam.xfISR(iSysNow, 21, xDaughter/z, pT2);
723 wt *= xPDFmotherNew / xPDFmotherOld;
724
725 // Iterate until acceptable pT and z.
726 } while (wt < Rndm::flat()) ;
727
728 // Select phi angle of branching at random.
729 double phi = 2. * M_PI * Rndm::flat();
730
731 // Save values for (so far) acceptable branching.
732 double mSister = (abs(idDaughter) == 4) ? mc : mb;
733 dipEndNow->store( idDaughter, 21, -idDaughter, x1Now, x2Now, m2Dip,
734 pT2, z, Q2, mSister, pow2(mSister), pT2corr, phi);
735
736}
737
738//*********
739
740// Evolve a QED dipole end.
741
742void SpaceShower::pT2nextQED( double pT2begDip, double pT2endDip) {
743
744 // Type of dipole and starting values.
745 BeamParticle& beam = ( abs(dipEndNow->side) == 1 )
746 ? *beamAPtr : *beamBPtr;
747 bool isLeptonBeam = beam.isLepton();
748 int MEtype = dipEndNow->MEtype;
749 bool isPhoton = (idDaughter == 22);
750 double pT2 = pT2begDip;
751 double m2Lepton = (isLeptonBeam) ? pow2(beam.m()) : 0.;
752 if (isLeptonBeam && pT2begDip < m2Lepton) return;
753
754 // Currently no f -> gamma branching implemented. ??
755 if (isPhoton) return;
756
757 // alpha_em at maximum scale provides upper estimate.
758 double alphaEMmax = alphaEM.alphaEM(pT2begDip);
759 double alphaEM2pi = alphaEMmax / (2. * M_PI);
760
761 // Maximum x of mother implies minimum z = xDaughter / xMother.
762 double xMaxAbs = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
763 double zMinAbs = xDaughter / xMaxAbs;
764
765 // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
766 double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
767 ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
768 if (isLeptonBeam) {
769 double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
770 if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
771 }
772 if (zMaxAbs < zMinAbs) return;
773
774 // Integrals of splitting kernels for fermions: f -> f. Use 1 + z^2 < 2.
775 // Ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
776 double f2fInt = 0.;
777 double f2fIntA = 2. * log( (1. - zMinAbs) / (1. - zMaxAbs) );
778 double f2fIntB = 0.;
779 if (isLeptonBeam) {
780 f2fIntB = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
781 f2fInt = f2fIntA + f2fIntB;
782 } else f2fInt = pow2(dipEndNow->chgType / 3.) * f2fIntA;
783
784 // Upper estimate for evolution equation, including fudge factor.
785 if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
786 double kernelPDF = alphaEM2pi * f2fInt;
787 double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
788 kernelPDF *= fudge;
789 if (kernelPDF < TINYKERNELPDF) return;
790
791 // Variables used inside evolution loop. (Mainly dummy start values.)
792 int idMother = 0;
793 double z = 0.;
794 double xMother = 0.;
795 double wt = 0.;
796 double Q2 = 0.;
797 double mSister = 0.;
798 double m2Sister = 0.;
799 double pT2corr = 0.;
800 double phi = 0.;
801
802 // Begin evolution loop towards smaller pT values.
803 do {
804 wt = 0.;
805
806 // Pick pT2 (in overestimated z range).
807 // For l -> l gamma include extrafactor 1 / ln(pT2 / m2l) in evolution.
808 double shift = pow(Rndm::flat(), 1. / kernelPDF);
809 if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
810 else pT2 = pT2 * shift;
811
812 // Abort evolution if below cutoff scale, or below another branching.
813 if (pT2 < pT2endDip) return;
814 if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
815
816 // Select z value of branching f -> f + gamma, and corrective weight.
817 idMother = idDaughter;
818 wt = 0.5 * (1. + pow2(z));
819 if (isLeptonBeam) {
820 if (f2fIntA > Rndm::flat() * (f2fIntA + f2fIntB))
821 z = 1. - (1. - zMinAbs)
822 * pow( (1. - zMaxAbs) / (1. - zMinAbs), Rndm::flat() );
823 else z = xDaughter + (zMinAbs - xDaughter)
824 * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter), Rndm::flat() );
825 wt *= (z - xDaughter) / (1. - xDaughter);
826 } else {
827 z = 1. - (1. - zMinAbs)
828 * pow( (1. - zMaxAbs) / (1. - zMinAbs), Rndm::flat() );
829 }
830
831 // Derive Q2 and x of mother from pT2 and z.
832 Q2 = pT2 / (1. - z);
833 xMother = xDaughter / z;
834
835 // Forbidden emission if outside allowed z range for given pT2.
836 mSister = 0.;
837 m2Sister = 0.;
838 pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
839 if(pT2corr < TINYPT2) { wt = 0.; continue; }
840
841 // Select phi angle of branching at random.
842 phi = 2. * M_PI * Rndm::flat();
843
844 // Correct by ln(pT2 / m2l) and fudge factor.
845 if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
846
847 // Evaluation of ME correction.
848 if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
849 m2Dip, z, Q2) / calcMEmax(MEtype, idMother, idDaughter);
850
851 // Optional dampening of large pT values in first radiation.
852 if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
853 wt *= pT2damp / (pT2 + pT2damp);
854
855 // Correct to current value of alpha_EM.
856 double alphaEMnow = alphaEM.alphaEM(pT2);
857 wt *= (alphaEMnow / alphaEMmax);
858
859 // Evaluation of new daughter and mother PDF's.
860 double xPDFdaughterNew = max ( TINYPDF,
861 beam.xfISR(iSysNow, idDaughter, xDaughter, pT2) );
862 double xPDFmotherNew = beam.xfISR(iSysNow, idMother, xMother, pT2);
863 wt *= xPDFmotherNew / xPDFdaughterNew;
864
865 // Iterate until acceptable pT (or have fallen below pTmin).
866 } while (wt < Rndm::flat()) ;
867
868 // Save values for (so far) acceptable branching.
869 dipEndNow->store( idDaughter,idMother, 22, x1Now, x2Now, m2Dip,
870 pT2, z, Q2, mSister, m2Sister, pT2corr, phi);
871
872}
873
874//*********
875
876// Kinematics of branching.
877// Construct mother -> daughter + sister, with recoiler on other side.
878
879bool SpaceShower::branch( Event& event) {
880
881 // Side on which branching occured.
882 int side = abs(dipEndSel->side);
883 double sideSign = (side == 1) ? 1. : -1.;
884
885 // Read in flavour and colour variables.
886 int iDaughter = event.getInSystem( iSysSel, side - 1);
887 int iRecoiler = event.getInSystem( iSysSel, 2 - side);
888 int idDaughterNow = dipEndSel->idDaughter;
889 int idMother = dipEndSel->idMother;
890 int idSister = dipEndSel->idSister;
891 int colDaughter = event[iDaughter].col();
892 int acolDaughter = event[iDaughter].acol();
893
894 // Read in kinematical variables.
895 double x1 = dipEndSel->x1;
896 double x2 = dipEndSel->x2;
897 double m2 = dipEndSel->m2Dip;
898 double m = sqrt(m2);
899 double pT2 = dipEndSel->pT2;
900 double z = dipEndSel->z;
901 double Q2 = dipEndSel->Q2;
902 double mSister = dipEndSel->mSister;
903 double m2Sister = dipEndSel->m2Sister;
904 double pT2corr = dipEndSel->pT2corr;
905 double phi = dipEndSel->phi;
906
907 // Take copy of existing system, to be given modified kinematics.
908 int eventSizeOld = event.size();
909 int systemSizeOld = event.sizeSystem(iSysSel);
910 for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
911 int iOldCopy = event.getInSystem( iSysSel, iCopy);
912 int statusNew = (iOldCopy == iDaughter
913 || iOldCopy == iRecoiler) ? event[iOldCopy].status() : 44;
914 event.copy(iOldCopy, statusNew);
915 }
916
917 // Define colour flow in branching.
918 // Default corresponds to f -> f + gamma.
919 int colMother = colDaughter;
920 int acolMother = acolDaughter;
921 int colSister = 0;
922 int acolSister = 0;
923 if (idSister == 22) ;
924 // q -> q + g and 50% of g -> g + g; need new colour.
925 else if (idSister == 21 && ( (idMother > 0 && idMother < 9)
926 || (idMother == 21 && Rndm::flat() < 0.5) ) ) {
927 colMother = event.nextColTag();
928 colSister = colMother;
929 acolSister = colDaughter;
930 // qbar -> qbar + g and other 50% of g -> g + g; need new colour.
931 } else if (idSister == 21) {
932 acolMother = event.nextColTag();
933 acolSister = acolMother;
934 colSister = acolDaughter;
935 // q -> g + q.
936 } else if (idDaughterNow == 21 && idMother > 0) {
937 colMother = colDaughter;
938 acolMother = 0;
939 colSister = acolDaughter;
940 // qbar -> g + qbar
941 } else if (idDaughterNow == 21) {
942 acolMother = acolDaughter;
943 colMother = 0;
944 acolSister = colDaughter;
945 // g -> q + qbar.
946 } else if (idDaughterNow > 0 && idDaughterNow < 9) {
947 acolMother = event.nextColTag();
948 acolSister = acolMother;
949 // g -> qbar + q.
950 } else if (idDaughterNow < 0 && idDaughterNow > -9) {
951 colMother = event.nextColTag();
952 colSister = colMother;
953 // q -> gamma + q.
954 } else if (idDaughterNow == 22 && idMother > 0) {
955 colMother = event.nextColTag();
956 colSister = colMother;
957 // qbar -> gamma + qbar.
958 } else if (idDaughterNow == 22) {
959 acolMother = event.nextColTag();
960 acolSister = acolMother;
961 }
962
963 // Construct kinematics of mother, sister and recoiler in old rest frame.
964 double pTbranch = sqrt(pT2corr) * m2 / ( z * (m2 + Q2) );
965 double pzMother = sideSign * 0.5 * m * ( (m2 - Q2) / ( z * (m2 + Q2) )
966 + (Q2 + m2Sister) / m2 );
967 double eMother = sqrt( pow2(pTbranch) + pow2(pzMother) );
968 double pzSister = pzMother - sideSign * 0.5 * (m2 + Q2) / m;
969 double eSister = sqrt( pow2(pTbranch) + pow2(pzSister) + m2Sister );
970 double eNewRecoiler = 0.5 * (m2 + Q2) / m;
971 Vec4 pMother( pTbranch, 0., pzMother, eMother );
972 Vec4 pSister( pTbranch, 0., pzSister, eSister );
973 Vec4 pNewRecoiler( 0., 0., -sideSign * eNewRecoiler, eNewRecoiler);
974
975 // Indices of partons involved. Add new sister.
976 int iMother = eventSizeOld + side - 1;
977 int iNewRecoiler = eventSizeOld + 2 - side;
978 int iSister = event.append( idSister, 43, iMother, 0, 0, 0,
979 colSister, acolSister, pSister, mSister, sqrt(pT2) );
980
981 // References to the partons involved.
982 Particle& daughter = event[iDaughter];
983 Particle& mother = event[iMother];
984 Particle& newRecoiler = event[iNewRecoiler];
985 Particle& sister = event.back();
986
987 // Replace old by new mother; update new recoiler.
988 mother.id( idMother );
989 mother.status( -41);
990 mother.cols( colMother, acolMother);
991 mother.p( pMother);
992 newRecoiler.status( -42);
993 newRecoiler.p( pNewRecoiler);
994
995 // Update mother and daughter pointers; also for beams.
996 daughter.mothers( iMother, 0);
997 mother.daughters( iSister, iDaughter);
998 if (iSysSel == 0) {
999 event[1].daughter1( (side == 1) ? iMother : iNewRecoiler );
1000 event[2].daughter1( (side == 2) ? iMother : iNewRecoiler );
1001 }
1002
1003 // Find boost to old rest frame, and rotation -phi.
1004 RotBstMatrix Mtot;
1005 Mtot.bst(0., 0., (x2 - x1) / (x1 + x2) );
1006 Mtot.rot(0., -phi);
1007
1008 // Find boost from old rest frame to event cm frame.
1009 RotBstMatrix MfromRest;
1010 // The boost to the new rest frame.
1011 Vec4 sumNew = pMother + pNewRecoiler;
1012 double betaX = sumNew.px() / sumNew.e();
1013 double betaZ = sumNew.pz() / sumNew.e();
1014 MfromRest.bst( -betaX, 0., -betaZ);
1015 // Alignment of radiator + recoiler to +- z axis, and rotation +phi.
1016 pMother.rotbst(MfromRest);
1017 double theta = pMother.theta();
1018 if (side == 2) theta += M_PI;
1019 MfromRest.rot(-theta, phi);
1020 // Longitudinal boost to radiator + recoiler in event cm frame.
1021 double x1New = (side == 1) ? x1 / z : x1;
1022 double x2New = (side == 2) ? x2 / z : x2;
1023 MfromRest.bst(0., 0., (x1New - x2New) / (x1New + x2New) );
1024 Mtot.rotbst(MfromRest);
1025
1026 // Perform cumulative rotation/boost operation.
1027 // Mother, recoiler and sister from old rest frame to event cm frame.
1028 mother.rotbst(MfromRest);
1029 newRecoiler.rotbst(MfromRest);
1030 sister.rotbst(MfromRest);
1031 // The rest from (and to) event cm frame.
1032 for ( int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i)
1033 event[i].rotbst(Mtot);
1034
1035 // Update list of partons in system; adding newly produced one.
1036 for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy)
1037 event.setInSystem( iSysSel, iCopy, eventSizeOld + iCopy);
1038 event.addToSystem( iSysSel, eventSizeOld + systemSizeOld);
1039
1040 // Update info on dipole ends (QCD or QED).
1041 for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
1042 if ( dipEnd[iDip].system == iSysSel && abs(dipEnd[iDip].side) == side ) {
1043 dipEnd[iDip].iRadiator = iMother;
1044 dipEnd[iDip].iRecoiler = iNewRecoiler;
1045 if (dipEnd[iDip].side > 0)
1046 dipEnd[iDip].colType = mother.colType();
1047 else {
1048 dipEnd[iDip].chgType = 0;
1049 if ( (mother.isQuark() && doQEDshowerByQ)
1050 || (mother.isLepton() && doQEDshowerByL) )
1051 dipEnd[iDip].chgType = mother.chargeType();
1052 }
1053 // Kill ME corrections after first emission.
1054 dipEnd[iDip].MEtype = 0;
1055 }
1056
1057 // Update info on beam remnants.
1058 BeamParticle& beamNow = (side == 1) ? *beamAPtr : *beamBPtr;
1059 double xNew = (side == 1) ? x1New : x2New;
1060 beamNow[iSysSel].update( iMother, idMother, xNew);
1061 // Redo choice of companion kind whenever new flavour.
1062 if (idMother != idDaughterNow) {
1063 beamNow.xfISR( iSysSel, idMother, xNew, pT2);
1064 beamNow.pickValSeaComp();
1065 }
1066 BeamParticle& beamRec = (side == 1) ? *beamBPtr : *beamAPtr;
1067 beamRec[iSysSel].iPos( iNewRecoiler);
1068
1069 // Store branching values of current dipole. (For rapidity ordering.)
1070 ++dipEndSel->nBranch;
1071 dipEndSel->pT2Old = pT2;
1072 dipEndSel->zOld = z;
1073
1074
1075 // Done without any errors.
1076 return true;
1077
1078}
1079
1080//*********
1081
1082// Find class of ME correction.
1083
1084int SpaceShower::findMEtype( int iSys, Event& event) {
1085
1086 // Default values and no action.
1087 int MEtype = 0;
1088 if (!doMEcorrections) ;
1089
1090 // Identify systems producing a single resonance.
1091 else if (event.sizeSystem( iSys) == 3) {
1092 int idIn1 = event[event.getInSystem( iSys, 0)].id();
1093 int idIn2 = event[event.getInSystem( iSys, 1)].id();
1094 int idRes = event[event.getInSystem( iSys, 2)].id();
1095
1096 // f + fbar -> vector boson.
1097 if ( (idRes == 23 || abs(idRes) == 24 || idRes == 32
1098 || idRes == 33 || abs(idRes) == 34 || abs(idRes) == 41)
1099 && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 1;
1100
1101 // g + g, gamma + gamma -> Higgs boson.
1102 if ( (idRes == 25 || idRes == 35 || idRes == 36)
1103 && ( ( idIn1 == 21 && idIn2 == 21 )
1104 || ( idIn1 == 22 && idIn2 == 22 ) ) ) MEtype = 2;
1105
1106 // f + fbar -> Higgs boson.
1107 if ( (idRes == 25 || idRes == 35 || idRes == 36)
1108 && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 3;
1109 }
1110
1111 // Done.
1112 return MEtype;
1113
1114}
1115
1116//*********
1117
1118// Provide maximum of expected ME weight; for preweighting of evolution.
1119
1120double SpaceShower::calcMEmax( int MEtype, int idMother, int idDaughterIn) {
1121
1122 // Currently only one non-unity case: g(gamma) f -> V f'.
1123 if (MEtype == 1 && idMother > 20 && idDaughterIn < 20) return 3.;
1124 return 1.;
1125
1126}
1127
1128//*********
1129
1130// Provide actual ME weight for current branching.
1131// Note: currently ME corrections are only allowed for first branching
1132// on each side, so idDaughter is essentially known and checks overkill.
1133
1134double SpaceShower::calcMEcorr(int MEtype, int idMother, int idDaughterIn,
1135 double M2, double z, double Q2) {
1136
1137 // Convert to Mandelstam variables. Sometimes may need to swap later.
1138 double sH = M2 / z;
1139 double tH = -Q2;
1140 double uH = Q2 - M2 * (1. - z) / z;
1141 int idMabs = abs(idMother);
1142 int idDabs = abs(idDaughterIn);
1143
1144 // Corrections for f + fbar -> s-channel vector boson.
1145 if (MEtype == 1) {
1146 if (idMabs < 20 && idDabs < 20) {
1147 return (tH*tH + uH*uH + 2. * M2 * sH) / (sH*sH + M2*M2);
1148 } else if (idDabs < 20) {
1149 // g(gamma) f -> V f': -Q2 = (p_g - p_f')^2 in PS while
1150 // tHat = (p_f - p_f')^2 in ME so need to swap tHat <-> uHat.
1151 swap( tH, uH);
1152 return (sH*sH + uH*uH + 2. * M2 * tH) / (pow2(sH - M2) + M2*M2);
1153 }
1154
1155 // Corrections for g + g -> Higgs boson.
1156 } else if (MEtype == 2) {
1157 if (idMabs < 20 && idDabs > 20) {
1158 return (sH*sH + uH*uH) / (sH*sH + pow2(sH - M2));
1159 } else if (idDabs > 20) {
1160 return 0.5 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(M2))
1161 / pow2(sH*sH - M2 * (sH - M2));
1162 }
1163
1164 // Corrections for f + fbar -> Higgs boson (f = b mainly).
1165 } else if (MEtype == 3) {
1166 if (idMabs < 20 && idDabs < 20) {
1167 // The PS and ME answers agree for f fbar -> H g/gamma.
1168 return 1.;
1169 } else if (idDabs < 20) {
1170 // Need to swap tHat <-> uHat, cf. vector-boson production above.
1171 swap( tH, uH);
1172 return (sH*sH + uH*uH + 2. * (M2 - uH) * (M2 - sH))
1173 / (pow2(sH - M2) + M2*M2);
1174 }
1175 }
1176
1177 return 1.;
1178
1179}
1180
1181//*********
1182
1183// Print the list of dipoles.
1184
1185void SpaceShower::list(ostream& os) {
1186
1187 // Header.
1188 os << "\n -------- PYTHIA SpaceShower Dipole Listing ---------- \n"
1189 << "\n i syst side rad rec pTmax col chg type \n"
1190 << fixed << setprecision(3);
1191
1192 // Loop over dipole list and print it.
1193 for (int i = 0; i < int(dipEnd.size()); ++i)
1194 os << setw(5) << i << setw(6) << dipEnd[i].system
1195 << setw(6) << dipEnd[i].side << setw(6) << dipEnd[i].iRadiator
1196 << setw(6) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax
1197 << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType
1198 << setw(5) << dipEnd[i].MEtype << "\n";
1199
1200 // Done.
1201 os << "\n -------- End PYTHIA SpaceShower Dipole Listing ------"
1202 << endl;
1203
1204
1205}
1206
1207//**************************************************************************
1208
1209} // end namespace Pythia8