Line data Source code
1 : // SpaceShower.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 : // SpaceShower class.
8 :
9 : #include "Pythia8/SpaceShower.h"
10 :
11 : namespace 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.
24 : const int SpaceShower::MAXLOOPTINYPDF = 10;
25 :
26 : // Minimal allowed c and b quark masses, for flavour thresholds.
27 : const double SpaceShower::MCMIN = 1.2;
28 : const double SpaceShower::MBMIN = 4.0;
29 :
30 : // Switch to alternative (but equivalent) backwards evolution for
31 : // g -> Q Qbar (Q = c or b) when below QTHRESHOLD * mQ2.
32 : const double SpaceShower::CTHRESHOLD = 2.0;
33 : const double SpaceShower::BTHRESHOLD = 2.0;
34 :
35 : // Renew evaluation of PDF's when the pT2 step is bigger than this
36 : // (in addition to initial scale and c and b thresholds.)
37 : const double SpaceShower::EVALPDFSTEP = 0.1;
38 :
39 : // Lower limit on PDF value in order to avoid division by zero.
40 : const double SpaceShower::TINYPDF = 1e-10;
41 :
42 : // Lower limit on estimated evolution rate, below which stop.
43 : const double SpaceShower::TINYKERNELPDF = 1e-6;
44 :
45 : // Lower limit on pT2, below which branching is rejected.
46 : const double SpaceShower::TINYPT2 = 0.25e-6;
47 :
48 : // No attempt to do backwards evolution of a heavy (c or b) quark
49 : // if evolution starts at a scale pT2 < HEAVYPT2EVOL * mQ2.
50 : const double SpaceShower::HEAVYPT2EVOL = 1.1;
51 :
52 : // No attempt to do backwards evolution of a heavy (c or b) quark
53 : // if evolution starts at a x > HEAVYXEVOL * x_max, where
54 : // x_max is the largest possible x value for a g -> Q Qbar branching.
55 : const double SpaceShower::HEAVYXEVOL = 0.9;
56 :
57 : // When backwards evolution Q -> g + Q creates a heavy quark Q,
58 : // an earlier branching g -> Q + Qbar will restrict kinematics
59 : // to M_{Q Qbar}^2 > EXTRASPACEQ * 4 m_Q^2. (Smarter to be found??)
60 : const double SpaceShower::EXTRASPACEQ = 2.0;
61 :
62 : // Never pick pT so low that alphaS is evaluated too close to Lambda_3.
63 : const double SpaceShower::LAMBDA3MARGIN = 1.1;
64 :
65 : // Do not warn for large PDF ratios at small pT2 scales.
66 : const double SpaceShower::PT2MINWARN = 1.;
67 :
68 : // Cutoff for f_e^e at x < 1 - 10^{-10} to be used in z selection.
69 : // Note: the x_min quantity come from 1 - x_max.
70 : const double SpaceShower::LEPTONXMIN = 1e-10;
71 : const double SpaceShower::LEPTONXMAX = 1. - 1e-10;
72 :
73 : // Stop l -> l gamma evolution slightly above m2l.
74 : const double SpaceShower::LEPTONPT2MIN = 1.2;
75 :
76 : // Enhancement of l -> l gamma trial rate to compensate imperfect modelling.
77 : const double SpaceShower::LEPTONFUDGE = 10.;
78 :
79 : // Overestimation extra factor for t-channel weak ME corrections.
80 : const double SpaceShower::WEAKPSWEIGHT = 5.;
81 :
82 : // Overestimation extra factors by branching type
83 : const double SpaceShower::HEADROOMQ2G = 1.35;
84 : const double SpaceShower::HEADROOMQ2Q = 1.15;
85 : const double SpaceShower::HEADROOMG2Q = 1.35;
86 : const double SpaceShower::HEADROOMG2G = 1.35;
87 :
88 : //--------------------------------------------------------------------------
89 :
90 : // Initialize alphaStrong, alphaEM and related pTmin parameters.
91 :
92 : void SpaceShower::init( BeamParticle* beamAPtrIn,
93 : BeamParticle* beamBPtrIn) {
94 :
95 : // Store input pointers for future use.
96 0 : beamAPtr = beamAPtrIn;
97 0 : beamBPtr = beamBPtrIn;
98 :
99 : // Main flags to switch on and off branchings.
100 0 : doQCDshower = settingsPtr->flag("SpaceShower:QCDshower");
101 0 : doQEDshowerByQ = settingsPtr->flag("SpaceShower:QEDshowerByQ");
102 0 : doQEDshowerByL = settingsPtr->flag("SpaceShower:QEDshowerByL");
103 0 : doWeakShower = settingsPtr->flag("SpaceShower:WeakShower");
104 :
105 : // Matching in pT of hard interaction to shower evolution.
106 0 : pTmaxMatch = settingsPtr->mode("SpaceShower:pTmaxMatch");
107 0 : pTdampMatch = settingsPtr->mode("SpaceShower:pTdampMatch");
108 0 : pTmaxFudge = settingsPtr->parm("SpaceShower:pTmaxFudge");
109 0 : pTmaxFudgeMPI = settingsPtr->parm("SpaceShower:pTmaxFudgeMPI");
110 0 : pTdampFudge = settingsPtr->parm("SpaceShower:pTdampFudge");
111 :
112 : // Optionally force emissions to be ordered in rapidity/angle.
113 0 : doRapidityOrder = settingsPtr->flag("SpaceShower:rapidityOrder");
114 :
115 : // Charm, bottom and lepton mass thresholds.
116 0 : mc = max( MCMIN, particleDataPtr->m0(4));
117 0 : mb = max( MBMIN, particleDataPtr->m0(5));
118 0 : m2c = pow2(mc);
119 0 : m2b = pow2(mb);
120 :
121 : // Parameters of scale choices.
122 0 : renormMultFac = settingsPtr->parm("SpaceShower:renormMultFac");
123 0 : factorMultFac = settingsPtr->parm("SpaceShower:factorMultFac");
124 0 : useFixedFacScale = settingsPtr->flag("SpaceShower:useFixedFacScale");
125 0 : fixedFacScale2 = pow2(settingsPtr->parm("SpaceShower:fixedFacScale"));
126 :
127 : // Parameters of alphaStrong generation.
128 0 : alphaSvalue = settingsPtr->parm("SpaceShower:alphaSvalue");
129 0 : alphaSorder = settingsPtr->mode("SpaceShower:alphaSorder");
130 0 : alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax");
131 0 : alphaSuseCMW = settingsPtr->flag("SpaceShower:alphaSuseCMW");
132 0 : alphaS2pi = 0.5 * alphaSvalue / M_PI;
133 :
134 : // Initialize alpha_strong generation.
135 0 : alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW);
136 :
137 : // Lambda for 5, 4 and 3 flavours.
138 0 : Lambda5flav = alphaS.Lambda5();
139 0 : Lambda4flav = alphaS.Lambda4();
140 0 : Lambda3flav = alphaS.Lambda3();
141 0 : Lambda5flav2 = pow2(Lambda5flav);
142 0 : Lambda4flav2 = pow2(Lambda4flav);
143 0 : Lambda3flav2 = pow2(Lambda3flav);
144 :
145 : // Regularization of QCD evolution for pT -> 0. Can be taken
146 : // same as for multiparton interactions, or be set separately.
147 0 : useSamePTasMPI = settingsPtr->flag("SpaceShower:samePTasMPI");
148 0 : if (useSamePTasMPI) {
149 0 : pT0Ref = settingsPtr->parm("MultipartonInteractions:pT0Ref");
150 0 : ecmRef = settingsPtr->parm("MultipartonInteractions:ecmRef");
151 0 : ecmPow = settingsPtr->parm("MultipartonInteractions:ecmPow");
152 0 : pTmin = settingsPtr->parm("MultipartonInteractions:pTmin");
153 0 : } else {
154 0 : pT0Ref = settingsPtr->parm("SpaceShower:pT0Ref");
155 0 : ecmRef = settingsPtr->parm("SpaceShower:ecmRef");
156 0 : ecmPow = settingsPtr->parm("SpaceShower:ecmPow");
157 0 : pTmin = settingsPtr->parm("SpaceShower:pTmin");
158 : }
159 :
160 : // Calculate nominal invariant mass of events. Set current pT0 scale.
161 0 : sCM = m2( beamAPtr->p(), beamBPtr->p());
162 0 : eCM = sqrt(sCM);
163 0 : pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
164 :
165 : // Restrict pTmin to ensure that alpha_s(pTmin^2 + pT_0^2) does not blow up.
166 0 : double pTminAbs = sqrtpos(pow2(LAMBDA3MARGIN) * Lambda3flav2 / renormMultFac
167 0 : - pT0*pT0);
168 0 : if (pTmin < pTminAbs) {
169 0 : pTmin = pTminAbs;
170 0 : ostringstream newPTmin;
171 0 : newPTmin << fixed << setprecision(3) << pTmin;
172 0 : infoPtr->errorMsg("Warning in SpaceShower::init: pTmin too low",
173 0 : ", raised to " + newPTmin.str() );
174 0 : infoPtr->setTooLowPTmin(true);
175 0 : }
176 :
177 : // Parameters of alphaEM generation.
178 0 : alphaEMorder = settingsPtr->mode("SpaceShower:alphaEMorder");
179 :
180 : // Initialize alphaEM generation.
181 0 : alphaEM.init( alphaEMorder, settingsPtr);
182 :
183 : // Parameters of QED evolution.
184 0 : pTminChgQ = settingsPtr->parm("SpaceShower:pTminchgQ");
185 0 : pTminChgL = settingsPtr->parm("SpaceShower:pTminchgL");
186 :
187 : // Derived parameters of QCD evolution.
188 0 : pT20 = pow2(pT0);
189 0 : pT2min = pow2(pTmin);
190 0 : pT2minChgQ = pow2(pTminChgQ);
191 0 : pT2minChgL = pow2(pTminChgL);
192 :
193 : // Parameters of weak evolution.
194 0 : weakMode = settingsPtr->mode("SpaceShower:weakShowerMode");
195 0 : pTweakCut = settingsPtr->parm("SpaceShower:pTminWeak");
196 0 : pT2weakCut = pow2(pTweakCut);
197 0 : weakEnhancement = settingsPtr->parm("WeakShower:enhancement");
198 0 : singleWeakEmission = settingsPtr->flag("WeakShower:singleEmission");
199 0 : vetoWeakJets = settingsPtr->flag("WeakShower:vetoWeakJets");
200 0 : vetoWeakDeltaR2 = pow2(settingsPtr->parm("weakShower:vetoWeakDeltaR"));
201 :
202 : // Various other parameters.
203 0 : doMEcorrections = settingsPtr->flag("SpaceShower:MEcorrections");
204 0 : doMEafterFirst = settingsPtr->flag("SpaceShower:MEafterFirst");
205 0 : doPhiPolAsym = settingsPtr->flag("SpaceShower:phiPolAsym");
206 0 : doPhiIntAsym = settingsPtr->flag("SpaceShower:phiIntAsym");
207 0 : strengthIntAsym = settingsPtr->parm("SpaceShower:strengthIntAsym");
208 0 : nQuarkIn = settingsPtr->mode("SpaceShower:nQuarkIn");
209 :
210 : // Z0 and W+- properties needed for weak showers.
211 0 : mZ = particleDataPtr->m0(23);
212 0 : gammaZ = particleDataPtr->mWidth(23);
213 0 : thetaWRat = 1. / (16. * coupSMPtr->sin2thetaW()
214 0 : * coupSMPtr->cos2thetaW());
215 0 : mW = particleDataPtr->m0(24);
216 0 : gammaW = particleDataPtr->mWidth(24);
217 :
218 : // Possibility of two predetermined hard emissions in event.
219 0 : doSecondHard = settingsPtr->flag("SecondHard:generate");
220 :
221 : // Optional dampening at small pT's when large multiplicities.
222 0 : enhanceScreening
223 0 : = settingsPtr->mode("MultipartonInteractions:enhanceScreening");
224 0 : if (!useSamePTasMPI) enhanceScreening = 0;
225 :
226 : // Possibility to allow user veto of emission step.
227 0 : canVetoEmission = (userHooksPtr != 0)
228 0 : ? userHooksPtr->canVetoISREmission() : false;
229 :
230 : // Default values for the weak shower.
231 0 : hasWeaklyRadiated = false;
232 0 : weakMaxWt = 1.;
233 :
234 0 : }
235 :
236 : //--------------------------------------------------------------------------
237 :
238 : // Find whether to limit maximum scale of emissions.
239 : // Also allow for dampening at factorization or renormalization scale.
240 :
241 : bool SpaceShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) {
242 :
243 : // Find whether to limit pT. Begin by user-set cases.
244 : bool dopTlimit = false;
245 0 : dopTlimit1 = dopTlimit2 = false;
246 : int nHeavyCol = 0;
247 0 : if (pTmaxMatch == 1) dopTlimit = dopTlimit1 = dopTlimit2 = true;
248 0 : else if (pTmaxMatch == 2) dopTlimit = dopTlimit1 = dopTlimit2 = false;
249 :
250 : // Always restrict SoftQCD processes.
251 0 : else if (infoPtr->isNonDiffractive() || infoPtr->isDiffractiveA()
252 0 : || infoPtr->isDiffractiveB() || infoPtr->isDiffractiveC() )
253 0 : dopTlimit = dopTlimit1 = dopTlimit2 = true;
254 :
255 : // Look if any quark (u, d, s, c, b), gluon or photon in final state.
256 : // Also count number of heavy coloured particles, like top.
257 : else {
258 : int n21 = 0;
259 : int iBegin = 5;
260 0 : if (infoPtr->isHardDiffractive()) iBegin = 9;
261 0 : for (int i = iBegin; i < event.size(); ++i) {
262 0 : if (event[i].status() == -21) ++n21;
263 0 : else if (n21 == 0) {
264 0 : int idAbs = event[i].idAbs();
265 0 : if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit1 = true;
266 0 : if ( (event[i].col() != 0 || event[i].acol() != 0)
267 0 : && idAbs > 5 && idAbs != 21 ) ++nHeavyCol;
268 0 : } else if (n21 == 2) {
269 0 : int idAbs = event[i].idAbs();
270 0 : if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit2 = true;
271 0 : }
272 : }
273 0 : dopTlimit = (doSecondHard) ? (dopTlimit1 && dopTlimit2) : dopTlimit1;
274 : }
275 :
276 : // Dampening at factorization or renormalization scale; only for hardest.
277 0 : dopTdamp = false;
278 0 : pT2damp = 0.;
279 0 : if (!dopTlimit1 && (pTdampMatch == 1 || pTdampMatch == 2)) {
280 0 : dopTdamp = true;
281 0 : pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren);
282 0 : }
283 0 : if (!dopTlimit1 && nHeavyCol > 1 && (pTdampMatch == 3 || pTdampMatch == 4)) {
284 0 : dopTdamp = true;
285 0 : pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 3) ? Q2Fac : Q2Ren);
286 0 : }
287 :
288 : // Done.
289 0 : return dopTlimit;
290 :
291 : }
292 :
293 : //--------------------------------------------------------------------------
294 :
295 : // Prepare system for evolution; identify ME.
296 : // Routine may be called after multiparton interactions, for a new subystem.
297 :
298 : void SpaceShower::prepare( int iSys, Event& event, bool limitPTmaxIn) {
299 :
300 : // Reset W/Z radiation flag and counters at first call for new event.
301 0 : if (iSys == 0) {
302 0 : nRadA.clear();
303 0 : nRadB.clear();
304 0 : hasWeaklyRadiated = false;
305 0 : }
306 :
307 : // Find positions of incoming colliding partons.
308 0 : int in1 = partonSystemsPtr->getInA(iSys);
309 0 : int in2 = partonSystemsPtr->getInB(iSys);
310 :
311 : // Rescattered partons cannot radiate.
312 0 : bool canRadiate1 = !(event[in1].isRescatteredIncoming());
313 0 : bool canRadiate2 = !(event[in2].isRescatteredIncoming());
314 :
315 : // Reset dipole-ends list for first interaction. Also resonances.
316 0 : if (iSys == 0) dipEnd.resize(0);
317 0 : if (iSys == 0) idResFirst = 0;
318 0 : if (iSys == 1) idResSecond = 0;
319 :
320 : // Find matrix element corrections for system.
321 0 : int MEtype = findMEtype( iSys, event, false);
322 :
323 : // In case of DPS overwrite limitPTmaxIn by saved value.
324 0 : if (doSecondHard && iSys == 0) limitPTmaxIn = dopTlimit1;
325 0 : if (doSecondHard && iSys == 1) limitPTmaxIn = dopTlimit2;
326 :
327 : // Maximum pT scale for dipole ends.
328 0 : double pTmax1 = (limitPTmaxIn) ? event[in1].scale() : eCM;
329 0 : double pTmax2 = (limitPTmaxIn) ? event[in2].scale() : eCM;
330 0 : if ( limitPTmaxIn && (iSys == 0 || (iSys == 1 && doSecondHard)) ) {
331 0 : pTmax1 *= pTmaxFudge;
332 0 : pTmax2 *= pTmaxFudge;
333 0 : } else if (limitPTmaxIn && iSys > 0) {
334 0 : pTmax1 *= pTmaxFudgeMPI;
335 0 : pTmax2 *= pTmaxFudgeMPI;
336 0 : }
337 :
338 : // Find dipole ends for QCD radiation.
339 : // Note: colour type can change during evolution, so book also if zero.
340 0 : if (doQCDshower) {
341 0 : int colType1 = event[in1].colType();
342 0 : if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys, 1,
343 0 : in1, in2, pTmax1, colType1, 0, 0, MEtype, canRadiate2) );
344 0 : int colType2 = event[in2].colType();
345 0 : if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys, 2,
346 : in2, in1, pTmax2, colType2, 0, 0, MEtype, canRadiate1) );
347 0 : }
348 :
349 : // Find dipole ends for QED radiation.
350 : // Note: charge type can change during evolution, so book also if zero.
351 0 : if (doQEDshowerByQ || doQEDshowerByL) {
352 0 : int chgType1 = ( (event[in1].isQuark() && doQEDshowerByQ)
353 0 : || (event[in1].isLepton() && doQEDshowerByL) )
354 0 : ? event[in1].chargeType() : 0;
355 : // Special: photons have charge zero, but can evolve (only off Q for now)
356 0 : if (event[in1].id() == 22 && doQEDshowerByQ) chgType1 = 22 ;
357 0 : if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys, -1,
358 0 : in1, in2, pTmax1, 0, chgType1, 0, MEtype, canRadiate2) );
359 0 : int chgType2 = ( (event[in2].isQuark() && doQEDshowerByQ)
360 0 : || (event[in2].isLepton() && doQEDshowerByL) )
361 0 : ? event[in2].chargeType() : 0;
362 : // Special: photons have charge zero, but can evolve (only off Q for now)
363 0 : if (event[in2].id() == 22 && doQEDshowerByQ) chgType2 = 22 ;
364 0 : if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys, -2,
365 : in2, in1, pTmax2, 0, chgType2, 0, MEtype, canRadiate1) );
366 0 : }
367 :
368 : // Find dipole ends for weak radiation. No right-handed W emission.
369 : // Currently leptons are not allow to emit W bosons and only
370 : // emissions from the hard process are included.
371 0 : if (doWeakShower && iSys == 0) {
372 :
373 : // Determine what type of 2 -> 2 process it is.
374 0 : int MEtypeWeak = findMEtype( iSys, event, true);
375 0 : if (MEtypeWeak == 201 || MEtypeWeak == 202 || MEtypeWeak == 203 ||
376 0 : MEtypeWeak == 206 || MEtypeWeak == 207 || MEtypeWeak == 208) {
377 :
378 : // Nonidentical incoming flavours.
379 0 : if (event[in1].id() != event[in2].id()) {
380 0 : if (event[in1].id() == event[in1 + 2].id()) tChannel = true;
381 0 : else if (event[in2].id() == event[in1 + 2].id()) tChannel = false;
382 : // No quark matches the outgoing, choose randomly.
383 0 : else tChannel = (rndmPtr->flat() < 0.5);
384 : // In case of same quark flavours, choose randomly.
385 0 : } else tChannel = (rndmPtr->flat() < 0.5);
386 : }
387 :
388 : // Set up weak dipole ends for first incoming parton.
389 0 : int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
390 0 : if (event[in1].idAbs() < 20) event[in1].pol(weakPol);
391 0 : if (canRadiate1) {
392 0 : if ( (weakMode == 0 || weakMode == 1) && weakPol == -1
393 0 : && event[in1].isQuark() )
394 0 : dipEnd.push_back( SpaceDipoleEnd( iSys, 1, in1, in2, pTmax1,
395 0 : 0, 0, 1, MEtypeWeak, canRadiate2, weakPol) );
396 0 : if ( (weakMode == 0 || weakMode == 2)
397 0 : && (event[in1].isQuark() || event[in1].isLepton()) )
398 0 : dipEnd.push_back( SpaceDipoleEnd( iSys, 1, in1, in2, pTmax1,
399 0 : 0, 0, 2, MEtypeWeak + 5, canRadiate2, weakPol) );
400 : }
401 :
402 : // Set up weak dipole ends for second incoming parton.
403 0 : if (event[in1].id() != - event[in2].id())
404 0 : weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
405 0 : if (event[in2].idAbs() < 20) event[in2].pol(weakPol);
406 0 : if (canRadiate2) {
407 0 : if ( (weakMode == 0 || weakMode == 1) && weakPol == -1
408 0 : && event[in2].isQuark())
409 0 : dipEnd.push_back( SpaceDipoleEnd( iSys, 2, in2, in1, pTmax2,
410 : 0, 0, 1, MEtypeWeak, canRadiate1, weakPol) );
411 0 : if ( (weakMode == 0 || weakMode == 2) &&
412 0 : (event[in2].isQuark() || event[in2].isLepton()) )
413 0 : dipEnd.push_back( SpaceDipoleEnd( iSys, 2, in2, in1, pTmax2,
414 0 : 0, 0, 2, MEtypeWeak + 5, canRadiate1, weakPol) );
415 : }
416 0 : }
417 :
418 : // Store the z and pT2 values of the last previous splitting
419 : // when an event history has already been constructed.
420 0 : if (iSys == 0 && infoPtr->hasHistory()) {
421 0 : double zNow = infoPtr->zNowISR();
422 0 : double pT2Now = infoPtr->pT2NowISR();
423 0 : for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
424 0 : dipEnd[iDipEnd].zOld = zNow;
425 0 : dipEnd[iDipEnd].pT2Old = pT2Now;
426 0 : ++dipEnd[iDipEnd].nBranch;
427 : }
428 0 : }
429 :
430 0 : }
431 :
432 : //--------------------------------------------------------------------------
433 :
434 : // Select next pT in downwards evolution of the existing dipoles.
435 :
436 : double SpaceShower::pTnext( Event& event, double pTbegAll, double pTendAll,
437 : int nRadIn) {
438 :
439 : // Current cm energy, in case it varies between events.
440 0 : sCM = m2( beamAPtr->p(), beamBPtr->p());
441 0 : eCM = sqrt(sCM);
442 0 : pTbegRef = pTbegAll;
443 :
444 : // Starting values: no radiating dipole found.
445 0 : nRad = nRadIn;
446 0 : double pT2sel = pow2(pTendAll);
447 0 : iDipSel = 0;
448 0 : iSysSel = 0;
449 0 : dipEndSel = 0;
450 :
451 : // Loop over all possible dipole ends.
452 0 : for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
453 0 : iDipNow = iDipEnd;
454 0 : dipEndNow = &dipEnd[iDipEnd];
455 0 : iSysNow = dipEndNow->system;
456 0 : dipEndNow->pT2 = 0.;
457 0 : double pTbegDip = min( pTbegAll, dipEndNow->pTmax );
458 :
459 : // Check whether dipole end should be allowed to shower.
460 0 : double pT2begDip = pow2(pTbegDip);
461 0 : if (pT2begDip > pT2sel && ( dipEndNow->colType != 0
462 0 : || dipEndNow->chgType != 0 || dipEndNow->weakType != 0) ) {
463 : double pT2endDip = 0.;
464 :
465 : // Determine lower cut for evolution, for QCD or QED (q or l).
466 0 : if (dipEndNow->colType != 0)
467 0 : pT2endDip = max( pT2sel, pT2min );
468 0 : else if (abs(dipEndNow->weakType) != 0)
469 0 : pT2endDip = max( pT2sel, pT2weakCut);
470 0 : else if (abs(dipEndNow->chgType) != 3 && dipEndNow->chgType != 0)
471 0 : pT2endDip = max( pT2sel, pT2minChgQ );
472 : else
473 0 : pT2endDip = max( pT2sel, pT2minChgL );
474 :
475 : // Find properties of dipole and radiating dipole end.
476 0 : sideA = ( abs(dipEndNow->side) == 1 );
477 0 : BeamParticle& beamNow = (sideA) ? *beamAPtr : *beamBPtr;
478 0 : BeamParticle& beamRec = (sideA) ? *beamBPtr : *beamAPtr;
479 0 : iNow = beamNow[iSysNow].iPos();
480 0 : iRec = beamRec[iSysNow].iPos();
481 0 : idDaughter = beamNow[iSysNow].id();
482 0 : xDaughter = beamNow[iSysNow].x();
483 0 : x1Now = (sideA) ? xDaughter : beamRec[iSysNow].x();
484 0 : x2Now = (sideA) ? beamRec[iSysNow].x() : xDaughter;
485 :
486 : // Note dipole mass correction when recoiler is a rescatter.
487 0 : m2Rec = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
488 0 : m2Dip = x1Now * x2Now * sCM + m2Rec;
489 :
490 : // Now do evolution in pT2, for QCD, QED or weak.
491 0 : if (pT2begDip > pT2endDip) {
492 0 : if (dipEndNow->colType != 0) pT2nextQCD( pT2begDip, pT2endDip);
493 0 : else if (dipEndNow->chgType != 0 || idDaughter == 22)
494 0 : pT2nextQED( pT2begDip, pT2endDip);
495 0 : else if (dipEndNow->weakType != 0) pT2nextWeak( pT2begDip, pT2endDip);
496 :
497 : // Update if found larger pT than current maximum.
498 0 : if (dipEndNow->pT2 > pT2sel) {
499 0 : pT2sel = dipEndNow->pT2;
500 0 : iDipSel = iDipNow;
501 0 : iSysSel = iSysNow;
502 0 : dipEndSel = dipEndNow;
503 0 : }
504 : }
505 0 : }
506 : // End loop over dipole ends.
507 0 : }
508 :
509 : // Return nonvanishing value if found pT is bigger than already found.
510 0 : return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
511 0 : }
512 :
513 : //--------------------------------------------------------------------------
514 :
515 : // Evolve a QCD dipole end.
516 :
517 : void SpaceShower::pT2nextQCD( double pT2begDip, double pT2endDip) {
518 :
519 : // Some properties and kinematical starting values.
520 0 : BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
521 0 : bool isGluon = (idDaughter == 21);
522 0 : bool isValence = beam[iSysNow].isValence();
523 0 : int MEtype = dipEndNow->MEtype;
524 0 : double pT2 = pT2begDip;
525 0 : double xMaxAbs = beam.xMax(iSysNow);
526 0 : double zMinAbs = xDaughter / xMaxAbs;
527 0 : if (xMaxAbs < 0.) {
528 0 : infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
529 : "xMaxAbs negative");
530 0 : return;
531 : }
532 :
533 : // Starting values for handling of massive quarks (c/b), if any.
534 : double idMassive = 0;
535 0 : if ( abs(idDaughter) == 4 ) idMassive = 4;
536 0 : if ( abs(idDaughter) == 5 ) idMassive = 5;
537 0 : bool isMassive = (idMassive > 0);
538 : double m2Massive = 0.;
539 : double mRatio = 0.;
540 0 : double zMaxMassive = 1.;
541 0 : double m2Threshold = pT2;
542 :
543 : // Evolution below scale of massive quark or at large x is impossible.
544 0 : if (isMassive) {
545 0 : m2Massive = (idMassive == 4) ? m2c : m2b;
546 0 : if (pT2 < HEAVYPT2EVOL * m2Massive) return;
547 0 : mRatio = sqrt( m2Massive / m2Dip );
548 0 : zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
549 0 : if (xDaughter > HEAVYXEVOL * zMaxMassive * xMaxAbs) return;
550 :
551 : // Find threshold scale below which only g -> Q + Qbar will be allowed.
552 0 : m2Threshold = (idMassive == 4) ? min( pT2, CTHRESHOLD * m2c)
553 0 : : min( pT2, BTHRESHOLD * m2b);
554 0 : }
555 :
556 : // Variables used inside evolution loop. (Mainly dummy starting values.)
557 : int nFlavour = 3;
558 : double b0 = 4.5;
559 0 : double Lambda2 = Lambda3flav2;
560 0 : double pT2minNow = pT2endDip;
561 : int idMother = 0;
562 : int idSister = 0;
563 0 : double z = 0.;
564 0 : double zMaxAbs = 0.;
565 : double zRootMax = 0.;
566 : double zRootMin = 0.;
567 : double g2gInt = 0.;
568 : double q2gInt = 0.;
569 : double q2qInt = 0.;
570 : double g2qInt = 0.;
571 : double g2Qenhance = 0.;
572 : double xPDFdaughter = 0.;
573 0 : double xPDFmother[21] = {0.};
574 : double xPDFgMother = 0.;
575 : double xPDFmotherSum = 0.;
576 : double kernelPDF = 0.;
577 : double xMother = 0.;
578 : double wt = 0.;
579 : double Q2 = 0.;
580 0 : double mSister = 0.;
581 : double m2Sister = 0.;
582 : double pT2corr = 0.;
583 0 : double pT2PDF = pT2;
584 : bool needNewPDF = true;
585 :
586 : // Begin evolution loop towards smaller pT values.
587 : int loopTinyPDFdau = 0;
588 : bool hasTinyPDFdau = false;
589 0 : do {
590 : wt = 0.;
591 :
592 : // Bad sign if repeated looping with small daughter PDF, so fail.
593 : // (Example: if all PDF's = 0 below Q_0, except for c/b companion.)
594 0 : if (hasTinyPDFdau) ++loopTinyPDFdau;
595 0 : if (loopTinyPDFdau > MAXLOOPTINYPDF) {
596 0 : infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
597 : "small daughter PDF");
598 0 : return;
599 : }
600 :
601 : // Initialize integrals of splitting kernels and evaluate parton
602 : // densities at the beginning. Reinitialize after long evolution
603 : // in pT2 or when crossing c and b flavour thresholds.
604 0 : if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
605 0 : pT2PDF = pT2;
606 : hasTinyPDFdau = false;
607 :
608 : // Determine overestimated z range; switch at c and b masses.
609 0 : if (pT2 > m2b) {
610 : nFlavour = 5;
611 0 : pT2minNow = max( m2b, pT2endDip);
612 : b0 = 23./6.;
613 0 : Lambda2 = Lambda5flav2;
614 0 : } else if (pT2 > m2c) {
615 : nFlavour = 4;
616 0 : pT2minNow = max( m2c, pT2endDip);
617 : b0 = 25./6.;
618 0 : Lambda2 = Lambda4flav2;
619 0 : } else {
620 : nFlavour = 3;
621 0 : pT2minNow = pT2endDip;
622 : b0 = 27./6.;
623 0 : Lambda2 = Lambda3flav2;
624 : }
625 : // A change of renormalization scale expressed by a change of Lambda.
626 0 : Lambda2 /= renormMultFac;
627 0 : zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
628 0 : ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
629 0 : if (isMassive) zMaxAbs = min( zMaxAbs, zMaxMassive);
630 :
631 : // Go to another z range with lower mass scale if current is closed.
632 0 : if (zMinAbs > zMaxAbs) {
633 0 : if (nFlavour == 3 || (idMassive == 4 && nFlavour == 4)
634 0 : || idMassive == 5) return;
635 0 : pT2 = (nFlavour == 4) ? m2c : m2b;
636 0 : continue;
637 : }
638 :
639 : // Parton density of daughter at current scale.
640 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
641 0 : xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2);
642 0 : if (xPDFdaughter < TINYPDF) {
643 : xPDFdaughter = TINYPDF;
644 : hasTinyPDFdau = true;
645 0 : }
646 :
647 : // Integrals of splitting kernels for gluons: g -> g, q -> g.
648 0 : if (isGluon) {
649 : g2gInt = HEADROOMG2G * 6.
650 0 : * log(zMaxAbs * (1.-zMinAbs) / (zMinAbs * (1.-zMaxAbs)));
651 0 : if (doMEcorrections) g2gInt *= calcMEmax(MEtype, 21, 21);
652 : q2gInt = HEADROOMQ2G * (16./3.)
653 0 : * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
654 0 : if (doMEcorrections) q2gInt *= calcMEmax(MEtype, 1, 21);
655 :
656 : // Parton density of potential quark mothers to a g.
657 : xPDFmotherSum = 0.;
658 0 : for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
659 0 : if (i == 0) {
660 0 : xPDFmother[10] = 0.;
661 0 : } else {
662 0 : xPDFmother[i+10] = beam.xfISR(iSysNow, i, xDaughter, pdfScale2);
663 0 : xPDFmotherSum += xPDFmother[i+10];
664 : }
665 : }
666 :
667 : // Total QCD evolution coefficient for a gluon.
668 0 : kernelPDF = g2gInt + q2gInt * xPDFmotherSum / xPDFdaughter;
669 :
670 : // For valence quark only need consider q -> q g branchings.
671 : // Introduce an extra factor sqrt(z) to smooth bumps.
672 0 : } else if (isValence) {
673 0 : zRootMin = (1. + sqrt(zMinAbs)) / (1. - sqrt(zMinAbs));
674 0 : zRootMax = (1. + sqrt(zMaxAbs)) / (1. - sqrt(zMaxAbs));
675 0 : q2qInt = (8./3.) * log( zRootMax / zRootMin );
676 0 : if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
677 : kernelPDF = q2qInt;
678 :
679 : // Integrals of splitting kernels for quarks: q -> q, g -> q.
680 0 : } else {
681 : q2qInt = HEADROOMQ2Q * (8./3.)
682 0 : * log( (1. - zMinAbs) / (1. - zMaxAbs) );
683 0 : if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
684 0 : g2qInt = HEADROOMG2Q * 0.5 * (zMaxAbs - zMinAbs);
685 0 : if (doMEcorrections) g2qInt *= calcMEmax(MEtype, 21, 1);
686 :
687 : // Increase estimated upper weight for g -> Q + Qbar.
688 0 : if (isMassive) {
689 0 : if (alphaSorder == 0) g2Qenhance = log(pT2/m2Massive)
690 0 : / log(m2Threshold/m2Massive);
691 : else {
692 0 : double m2log = log( m2Massive / Lambda2);
693 0 : g2Qenhance = log( log(pT2/Lambda2) / m2log )
694 0 : / log( log(m2Threshold/Lambda2) / m2log );
695 : }
696 0 : g2qInt *= g2Qenhance;
697 0 : }
698 :
699 : // Parton density of a potential gluon mother to a q.
700 0 : xPDFgMother = beam.xfISR(iSysNow, 21, xDaughter, pdfScale2);
701 :
702 : // Total QCD evolution coefficient for a quark.
703 0 : kernelPDF = q2qInt + g2qInt * xPDFgMother / xPDFdaughter;
704 : }
705 :
706 : // End evaluation of splitting kernels and parton densities.
707 : needNewPDF = false;
708 0 : }
709 0 : if (kernelPDF < TINYKERNELPDF) return;
710 :
711 : // Pick pT2 (in overestimated z range), for one of three different cases.
712 : // Assume form alphas(pT0^2 + pT^2) * dpT^2/(pT0^2 + pT^2).
713 : double Q2alphaS;
714 :
715 : // Fixed alpha_strong.
716 0 : if (alphaSorder == 0) {
717 0 : pT2 = (pT2 + pT20) * pow( rndmPtr->flat(),
718 0 : 1. / (alphaS2pi * kernelPDF)) - pT20;
719 :
720 : // First-order alpha_strong.
721 0 : } else if (alphaSorder == 1) {
722 0 : pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
723 0 : pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
724 :
725 : // For second order reject by second term in alpha_strong expression.
726 0 : } else {
727 : do {
728 0 : pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
729 0 : pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
730 0 : Q2alphaS = renormMultFac * max( pT2 + pT20,
731 0 : pow2(LAMBDA3MARGIN) * Lambda3flav2);
732 0 : } while (alphaS.alphaS2OrdCorr(Q2alphaS) < rndmPtr->flat()
733 0 : && pT2 > pT2minNow);
734 : }
735 :
736 : // Check for pT2 values that prompt special action.
737 :
738 : // If fallen into b threshold region, force g -> b + bbar.
739 0 : if (idMassive == 5 && pT2 < m2Threshold) {
740 0 : pT2nearQCDthreshold( beam, m2Massive, m2Threshold, xMaxAbs,
741 0 : zMinAbs, zMaxMassive );
742 0 : return;
743 :
744 : // If crossed b threshold, continue evolution from this threshold.
745 0 : } else if (nFlavour == 5 && pT2 < m2b) {
746 : needNewPDF = true;
747 0 : pT2 = m2b;
748 0 : continue;
749 :
750 : // If fallen into c threshold region, force g -> c + cbar.
751 0 : } else if (idMassive == 4 && pT2 < m2Threshold) {
752 0 : pT2nearQCDthreshold( beam, m2Massive, m2Threshold, xMaxAbs,
753 0 : zMinAbs, zMaxMassive );
754 0 : return;
755 :
756 : // If crossed c threshold, continue evolution from this threshold.
757 0 : } else if (nFlavour == 4 && pT2 < m2c) {
758 : needNewPDF = true;
759 0 : pT2 = m2c;
760 0 : continue;
761 :
762 : // Abort evolution if below cutoff scale, or below another branching.
763 0 : } else if (pT2 < pT2endDip) return;
764 :
765 : // Select z value of branching to g, and corrective weight.
766 0 : if (isGluon) {
767 : // g -> g (+ g).
768 0 : if (rndmPtr->flat() * kernelPDF < g2gInt) {
769 : idMother = 21;
770 : idSister = 21;
771 0 : z = 1. / ( 1. + ((1. - zMinAbs) / zMinAbs) * pow( (zMinAbs *
772 0 : (1. - zMaxAbs)) / (zMaxAbs * (1. - zMinAbs)), rndmPtr->flat() ) );
773 0 : wt = pow2( 1. - z * (1. - z));
774 : // Account for headroom factor used to enhance trial probability
775 0 : wt /= HEADROOMG2G;
776 0 : } else {
777 : // q -> g (+ q): also select flavour.
778 0 : double temp = xPDFmotherSum * rndmPtr->flat();
779 0 : idMother = -nQuarkIn - 1;
780 0 : do { temp -= xPDFmother[(++idMother) + 10]; }
781 0 : while (temp > 0. && idMother < nQuarkIn);
782 : idSister = idMother;
783 0 : z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
784 0 : * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
785 0 : wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z)
786 0 : * xPDFdaughter / xPDFmother[idMother + 10];
787 : // Account for headroom factor used to enhance trial probability
788 0 : wt /= HEADROOMQ2G;
789 : }
790 :
791 : // Select z value of branching to q, and corrective weight.
792 : // Include massive kernel corrections for c and b quarks.
793 : } else {
794 : // q -> q (+ g).
795 0 : if (isValence || rndmPtr->flat() * kernelPDF < q2qInt) {
796 0 : idMother = idDaughter;
797 : idSister = 21;
798 : // Valence more peaked at large z.
799 0 : if (isValence) {
800 0 : double zTmp = zRootMin * pow(zRootMax / zRootMin, rndmPtr->flat() );
801 0 : z = pow2( (1. - zTmp) / (1. + zTmp) );
802 0 : } else {
803 0 : z = 1. - (1. - zMinAbs) * pow( (1. - zMaxAbs) / (1. - zMinAbs),
804 0 : rndmPtr->flat() );
805 : }
806 0 : if (!isMassive) {
807 : wt = 0.5 * (1. + pow2(z));
808 0 : } else {
809 : //?? Bug?? should be 2 more for massive part??
810 : // wt = 0.5 * (1. + pow2(z) - z * pow2(1.-z) * m2Massive / pT2);
811 0 : wt = 0.5 * (1. + pow2(z)) - z * pow2(1.-z) * m2Massive / pT2;
812 : }
813 0 : if (isValence) wt *= sqrt(z);
814 : // Account for headroom factor for sea quarks
815 0 : else wt /= HEADROOMQ2Q;
816 : // g -> q (+ qbar).
817 : } else {
818 : idMother = 21;
819 0 : idSister = - idDaughter;
820 0 : z = zMinAbs + rndmPtr->flat() * (zMaxAbs - zMinAbs);
821 0 : if (!isMassive) {
822 0 : wt = (pow2(z) + pow2(1.-z)) * xPDFdaughter / xPDFgMother;
823 0 : } else {
824 0 : wt = (pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2)
825 0 : * xPDFdaughter / (xPDFgMother * g2Qenhance) ;
826 : }
827 : // Account for headroom factor for gluons
828 0 : wt /= HEADROOMG2Q;
829 : }
830 : }
831 :
832 : // Derive Q2 and x of mother from pT2 and z.
833 0 : Q2 = pT2 / (1. - z);
834 0 : xMother = xDaughter / z;
835 : // Correction to x for massive recoiler from rescattering.
836 0 : if (!dipEndNow->normalRecoil) {
837 0 : if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
838 0 : else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
839 : }
840 0 : if (xMother > xMaxAbs) { wt = 0.; continue; }
841 :
842 : // Forbidden emission if outside allowed z range for given pT2.
843 0 : mSister = particleDataPtr->m0(idSister);
844 0 : m2Sister = pow2(mSister);
845 0 : pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
846 0 : if (pT2corr < TINYPT2) { wt = 0.; continue; }
847 :
848 : // Optionally veto emissions not ordered in rapidity (= angle).
849 0 : if ( doRapidityOrder && dipEndNow->nBranch > 0
850 0 : && pT2 > pow2( (1. - z) / (z * (1. - dipEndNow->zOld)) )
851 0 : * dipEndNow->pT2Old ) { wt = 0.; continue; }
852 :
853 : // If creating heavy quark by Q -> g + Q then next need g -> Q + Qbar.
854 : // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
855 0 : if ( isGluon && ( abs(idMother) == 4 || abs(idMother) == 5 )) {
856 0 : double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
857 0 : double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
858 0 : if (pT2QQcorr < TINYPT2) { wt = 0.; continue; }
859 0 : }
860 :
861 : // Evaluation of ME correction.
862 0 : if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
863 0 : m2Dip, z, Q2, m2Sister) / calcMEmax(MEtype, idMother, idDaughter);
864 :
865 : // Optional dampening of large pT values in first radiation.
866 0 : if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
867 0 : wt *= pT2damp / (pT2 + pT2damp);
868 :
869 : // Idea suggested by Gosta Gustafson: increased screening in events
870 : // with large activity can be simulated by pT0_eff = sqrt(n) * pT0.
871 0 : if (enhanceScreening == 2) {
872 0 : int nSysNow = infoPtr->nMPI() + infoPtr->nISR() + 1;
873 0 : double WTscreen = pow2( (pT2 + pT20) / (pT2 + nSysNow * pT20) );
874 0 : wt *= WTscreen;
875 0 : }
876 :
877 : // Evaluation of new daughter and mother PDF's.
878 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
879 0 : double xPDFdaughterNew = max ( TINYPDF,
880 0 : beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
881 : double xPDFmotherNew =
882 0 : beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
883 0 : wt *= xPDFmotherNew / xPDFdaughterNew;
884 :
885 : // Check that valence step does not cause problem.
886 0 : if (wt > 1. && pT2 > PT2MINWARN) infoPtr->errorMsg("Warning in "
887 : "SpaceShower::pT2nextQCD: weight above unity");
888 :
889 : // Iterate until acceptable pT (or have fallen below pTmin).
890 0 : } while (wt < rndmPtr->flat()) ;
891 :
892 : // Save values for (so far) acceptable branching.
893 0 : dipEndNow->store( idDaughter,idMother, idSister, x1Now, x2Now, m2Dip,
894 0 : pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
895 :
896 0 : }
897 :
898 : //--------------------------------------------------------------------------
899 :
900 : // Evolve a QCD dipole end near threshold, with g -> Q + Qbar enforced.
901 : // Note: No explicit Sudakov factor formalism here. Instead use that
902 : // df_Q(x, pT2) = (alpha_s/2pi) * (dT2/pT2) * ((gluon) * (splitting)).
903 : // This implies that effects of Q -> Q + g are neglected in this range.
904 :
905 : void SpaceShower::pT2nearQCDthreshold( BeamParticle& beam,
906 : double m2Massive, double m2Threshold, double xMaxAbs,
907 : double zMinAbs, double zMaxMassive) {
908 :
909 : // Initial values, to be used in kinematics and weighting.
910 0 : double Lambda2 = (abs(idDaughter) == 4) ? Lambda4flav2 : Lambda5flav2;
911 0 : Lambda2 /= renormMultFac;
912 0 : double logM2Lambda2 = (alphaSorder > 0) ? log( m2Massive / Lambda2 ) : 1.;
913 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2
914 0 : : factorMultFac * m2Threshold;
915 0 : double xPDFmotherOld = beam.xfISR(iSysNow, 21, xDaughter, pdfScale2);
916 :
917 : // Variables used inside evolution loop. (Mainly dummy start values.)
918 : int loop = 0;
919 : double wt = 0.;
920 : double pT2 = 0.;
921 0 : double z = 0.;
922 : double Q2 = 0.;
923 : double pT2corr = 0.;
924 : double xMother = 0.;
925 :
926 : // Begin loop over tries to find acceptable g -> Q + Qbar branching.
927 0 : do {
928 : wt = 0.;
929 :
930 : // Check that not caught in infinite loop with impossible kinematics.
931 0 : if (++loop > 100) {
932 0 : infoPtr->errorMsg("Error in SpaceShower::pT2nearQCDthreshold: "
933 : "stuck in loop");
934 0 : return;
935 : }
936 :
937 : // Pick dpT2/pT2 in range [m2Massive,thresholdRatio * m2Massive].
938 0 : pT2 = m2Massive * pow( m2Threshold / m2Massive, rndmPtr->flat() );
939 :
940 : // Pick z flat in allowed range.
941 0 : z = zMinAbs + rndmPtr->flat() * (zMaxMassive - zMinAbs);
942 :
943 : // Check that kinematically possible choice.
944 0 : Q2 = pT2 / (1.-z) - m2Massive;
945 0 : pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Massive) / m2Dip;
946 0 : if (pT2corr < TINYPT2) continue;
947 :
948 : // Correction factor for running alpha_s. (Only first order for now.)
949 0 : wt = (alphaSorder > 0) ? logM2Lambda2 / log( pT2 / Lambda2 ) : 1.;
950 :
951 : // Correction factor for splitting kernel.
952 0 : wt *= pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2;
953 :
954 : // x, including correction for massive recoiler from rescattering.
955 0 : xMother = xDaughter / z;
956 0 : if (!dipEndNow->normalRecoil) {
957 0 : if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
958 0 : else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
959 : }
960 0 : if (xMother > xMaxAbs) { wt = 0.; continue; }
961 :
962 : // Correction factor for gluon density.
963 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
964 0 : double xPDFmotherNew = beam.xfISR(iSysNow, 21, xMother, pdfScale2);
965 0 : wt *= xPDFmotherNew / xPDFmotherOld;
966 :
967 : // Iterate until acceptable pT and z.
968 0 : } while (wt < rndmPtr->flat()) ;
969 :
970 : // Save values for (so far) acceptable branching.
971 0 : double mSister = (abs(idDaughter) == 4) ? mc : mb;
972 0 : dipEndNow->store( idDaughter, 21, -idDaughter, x1Now, x2Now, m2Dip,
973 0 : pT2, z, xMother, Q2, mSister, pow2(mSister), pT2corr);
974 :
975 0 : }
976 :
977 : //--------------------------------------------------------------------------
978 :
979 : // Evolve a QED dipole end.
980 :
981 : void SpaceShower::pT2nextQED( double pT2begDip, double pT2endDip) {
982 :
983 : // Type of dipole and starting values.
984 0 : BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
985 0 : bool isLeptonBeam = beam.isLepton();
986 0 : int MEtype = dipEndNow->MEtype;
987 0 : bool isPhoton = (idDaughter == 22);
988 : double pT2 = pT2begDip;
989 0 : double m2Lepton = (isLeptonBeam) ? pow2(beam.m()) : 0.;
990 0 : if (isLeptonBeam && pT2begDip < m2Lepton) return;
991 :
992 : // Currently no f -> gamma branching implemented for lepton beams
993 0 : if (isPhoton && isLeptonBeam) return;
994 :
995 : // alpha_em at maximum scale provides upper estimate.
996 0 : double alphaEMmax = alphaEM.alphaEM(renormMultFac * pT2begDip);
997 0 : double alphaEM2pi = alphaEMmax / (2. * M_PI);
998 :
999 : // Maximum x of mother implies minimum z = xDaughter / xMother.
1000 0 : double xMaxAbs = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
1001 0 : double zMinAbs = xDaughter / xMaxAbs;
1002 0 : if (xMaxAbs < 0.) {
1003 0 : infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
1004 : "xMaxAbs negative");
1005 0 : return;
1006 : }
1007 :
1008 : // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
1009 0 : double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
1010 0 : ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
1011 0 : if (isLeptonBeam) {
1012 0 : double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
1013 0 : if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
1014 0 : }
1015 0 : if (zMaxAbs < zMinAbs) return;
1016 :
1017 : // Variables used inside evolution loop. (Mainly dummy start values.)
1018 : int idMother = 0;
1019 : int idSister = 22;
1020 0 : double z = 0.;
1021 : double xMother = 0.;
1022 : double wt = 0.;
1023 : double Q2 = 0.;
1024 0 : double mSister = 0.;
1025 : double m2Sister = 0.;
1026 : double pT2corr = 0.;
1027 :
1028 : // QED evolution of fermions
1029 0 : if (!isPhoton) {
1030 :
1031 : // Integrals of splitting kernels for fermions: f -> f. Use 1 + z^2 < 2.
1032 : // Ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
1033 : double f2fInt = 0.;
1034 0 : double f2fIntA = 2. * log( (1. - zMinAbs) / (1. - zMaxAbs) );
1035 : double f2fIntB = 0.;
1036 0 : if (isLeptonBeam) {
1037 0 : f2fIntB = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
1038 0 : f2fInt = f2fIntA + f2fIntB;
1039 0 : } else f2fInt = pow2(dipEndNow->chgType / 3.) * f2fIntA;
1040 :
1041 : // Upper estimate for evolution equation, including fudge factor.
1042 0 : if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
1043 0 : double kernelPDF = alphaEM2pi * f2fInt;
1044 0 : double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
1045 0 : kernelPDF *= fudge;
1046 0 : if (kernelPDF < TINYKERNELPDF) return;
1047 :
1048 : // Begin evolution loop towards smaller pT values.
1049 : do {
1050 :
1051 : // Pick pT2 (in overestimated z range).
1052 : // For l -> l gamma include extrafactor 1 / ln(pT2 / m2l) in evolution.
1053 0 : double shift = pow(rndmPtr->flat(), 1. / kernelPDF);
1054 0 : if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
1055 0 : else pT2 = pT2 * shift;
1056 :
1057 : // Abort evolution if below cutoff scale, or below another branching.
1058 0 : if (pT2 < pT2endDip) return;
1059 0 : if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
1060 :
1061 : // Select z value of branching f -> f + gamma, and corrective weight.
1062 0 : idMother = idDaughter;
1063 : wt = 1.;
1064 0 : if (isLeptonBeam) {
1065 0 : if (f2fIntA > rndmPtr->flat() * (f2fIntA + f2fIntB)) {
1066 0 : z = 1. - (1. - zMinAbs)
1067 0 : * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
1068 0 : } else {
1069 0 : z = xDaughter + (zMinAbs - xDaughter)
1070 0 : * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter),
1071 0 : rndmPtr->flat() );
1072 : }
1073 0 : wt *= (z - xDaughter) / (1. - xDaughter);
1074 0 : } else {
1075 0 : z = 1. - (1. - zMinAbs)
1076 0 : * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
1077 : }
1078 0 : wt *= 0.5 * (1. + pow2(z));
1079 :
1080 : // Derive Q2 and x of mother from pT2 and z.
1081 0 : Q2 = pT2 / (1. - z);
1082 0 : xMother = xDaughter / z;
1083 : // Correction to x for massive recoiler from rescattering.
1084 0 : if (!dipEndNow->normalRecoil) {
1085 0 : if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1086 0 : else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1087 : }
1088 0 : if (xMother > xMaxAbs) { wt = 0.; continue; }
1089 :
1090 : // Forbidden emission if outside allowed z range for given pT2.
1091 0 : mSister = 0.;
1092 : m2Sister = 0.;
1093 0 : pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1094 0 : if (pT2corr < TINYPT2) { wt = 0.; continue; }
1095 :
1096 : // Correct by ln(pT2 / m2l) and fudge factor.
1097 0 : if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
1098 :
1099 : // Evaluation of ME correction.
1100 0 : if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
1101 0 : m2Dip, z, Q2, m2Sister) / calcMEmax(MEtype, idMother, idDaughter);
1102 :
1103 : // Extra QED correction for f fbar -> W+- gamma. Debug??
1104 0 : if (doMEcorrections && MEtype == 1 && idDaughter == idMother
1105 0 : && ( (iSysNow == 0 && idResFirst == 24)
1106 0 : || (iSysNow == 1 && idResSecond == 24) ) ) {
1107 0 : double tHe = -Q2;
1108 0 : double uHe = Q2 - m2Dip * (1. - z) / z;
1109 0 : double chg1 = abs(dipEndNow->chgType / 3.);
1110 0 : double chg2 = 1. - chg1;
1111 0 : wt *= pow2(chg1 * uHe - chg2 * tHe)
1112 0 : / ( (tHe + uHe) * (pow2(chg1) * uHe + pow2(chg2) * tHe) );
1113 0 : }
1114 :
1115 : // Optional dampening of large pT values in first radiation.
1116 0 : if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
1117 0 : wt *= pT2damp / (pT2 + pT2damp);
1118 :
1119 : // Correct to current value of alpha_EM.
1120 0 : double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
1121 0 : wt *= (alphaEMnow / alphaEMmax);
1122 :
1123 : // Evaluation of new daughter and mother PDF's.
1124 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1125 0 : double xPDFdaughterNew = max ( TINYPDF,
1126 0 : beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
1127 : double xPDFmotherNew =
1128 0 : beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
1129 0 : wt *= xPDFmotherNew / xPDFdaughterNew;
1130 :
1131 : // Iterate until acceptable pT (or have fallen below pTmin).
1132 0 : } while (wt < rndmPtr->flat()) ;
1133 0 : }
1134 :
1135 : // QED evolution of photons (so far only for hadron beams).
1136 : else {
1137 :
1138 : // Initial values
1139 : int nFlavour = 3;
1140 : double kernelPDF = 0.0;
1141 : double xPDFdaughter = 0.;
1142 0 : double xPDFmother[21] = {0.};
1143 : double xPDFmotherSum = 0.0;
1144 : double pT2PDF = pT2;
1145 : double pT2minNow = pT2endDip;
1146 : bool needNewPDF = true;
1147 :
1148 : // Begin evolution loop towards smaller pT values.
1149 : int loopTinyPDFdau = 0;
1150 : bool hasTinyPDFdau = false;
1151 0 : do {
1152 : wt = 0.;
1153 :
1154 : // Bad sign if repeated looping with small daughter PDF, so fail.
1155 0 : if (hasTinyPDFdau) ++loopTinyPDFdau;
1156 0 : if (loopTinyPDFdau > MAXLOOPTINYPDF) {
1157 0 : infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
1158 : "small daughter PDF");
1159 0 : return;
1160 : }
1161 :
1162 : // Initialize integrals of splitting kernels and evaluate parton
1163 : // densities at the beginning. Reinitialize after long evolution
1164 : // in pT2 or when crossing c and b flavour thresholds.
1165 0 : if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
1166 :
1167 : pT2PDF = pT2;
1168 : hasTinyPDFdau = false;
1169 :
1170 : // Determine overestimated z range; switch at c and b masses.
1171 0 : if (pT2 > m2b && nQuarkIn >= 5) {
1172 : nFlavour = 5;
1173 0 : pT2minNow = max( m2b, pT2endDip);
1174 0 : } else if (pT2 > m2c && nQuarkIn >= 4) {
1175 : nFlavour = 4;
1176 0 : pT2minNow = max( m2c, pT2endDip);
1177 0 : } else {
1178 : nFlavour = 3;
1179 0 : pT2minNow = pT2endDip;
1180 : }
1181 :
1182 : // Compute upper z limit
1183 0 : zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
1184 0 : ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
1185 :
1186 : // Parton density of daughter at current scale.
1187 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1188 0 : xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2);
1189 0 : if (xPDFdaughter < TINYPDF) {
1190 : xPDFdaughter = TINYPDF;
1191 : hasTinyPDFdau = true;
1192 0 : }
1193 :
1194 : // Integral over f -> gamma f splitting kernel.
1195 : // Normalized so: 4/3 aS/2pi P(z) -> eq^2 * aEM/2pi P(z).
1196 : // (Charge-weighting happens below.)
1197 0 : double q2gInt = 4. * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
1198 :
1199 : // Charge-weighted Parton density of potential quark mothers.
1200 : xPDFmotherSum = 0.;
1201 0 : for (int i = -nFlavour; i <= nFlavour; ++i) {
1202 0 : if (i == 0) {
1203 0 : xPDFmother[10] = 0.;
1204 0 : } else {
1205 0 : xPDFmother[i+10] = pow2((abs(i+1) % 2 + 1)/3.0)
1206 0 : * beam.xfISR(iSysNow, i, xDaughter, pdfScale2);
1207 0 : xPDFmotherSum += xPDFmother[i+10];
1208 : }
1209 : }
1210 :
1211 : // Total QED evolution coefficient for a photon.
1212 0 : kernelPDF = q2gInt * xPDFmotherSum / xPDFdaughter;
1213 :
1214 : // End evaluation of splitting kernels and parton densities.
1215 : needNewPDF = false;
1216 0 : }
1217 0 : if (kernelPDF < TINYKERNELPDF) return;
1218 :
1219 : // Select pT2 for next trial branching
1220 0 : pT2 *= pow( rndmPtr->flat(), 1. / (alphaEM2pi * kernelPDF));
1221 :
1222 : // If crossed b threshold, continue evolution from this threshold.
1223 0 : if (nFlavour == 5 && pT2 < m2b) {
1224 : needNewPDF = true;
1225 : pT2 = m2b;
1226 0 : continue;
1227 : }
1228 :
1229 : // If crossed c threshold, continue evolution from this threshold.
1230 0 : else if (nFlavour == 4 && pT2 < m2c) {
1231 : needNewPDF = true;
1232 : pT2 = m2c;
1233 0 : continue;
1234 : }
1235 :
1236 : // Abort evolution if below cutoff scale, or below another branching.
1237 0 : else if (pT2 < pT2endDip) return;
1238 :
1239 : // Select flavour for trial branching
1240 0 : double temp = xPDFmotherSum * rndmPtr->flat();
1241 0 : idMother = -nQuarkIn - 1;
1242 0 : do {
1243 0 : temp -= xPDFmother[(++idMother) + 10];
1244 0 : } while (temp > 0. && idMother < nQuarkIn);
1245 :
1246 : // Sister is same as mother, but can have m2 > 0
1247 : idSister = idMother;
1248 0 : mSister = particleDataPtr->m0(idSister);
1249 0 : m2Sister = pow2(mSister);
1250 :
1251 : // Select z for trial branching
1252 0 : z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
1253 0 : * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
1254 :
1255 : // Trial weight: splitting kernel
1256 0 : wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z);
1257 :
1258 : // Trial weight: running alpha_EM
1259 0 : double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
1260 0 : wt *= (alphaEMnow / alphaEMmax);
1261 :
1262 : // Derive Q2 and x of mother from pT2 and z
1263 0 : Q2 = pT2 / (1. - z);
1264 0 : xMother = xDaughter / z;
1265 : // Correction to x for massive recoiler from rescattering.
1266 0 : if (!dipEndNow->normalRecoil) {
1267 0 : if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1268 0 : else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1269 : }
1270 :
1271 : // Compute pT2corr
1272 0 : pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1273 0 : if (pT2corr < TINYPT2) { wt = 0.; continue; }
1274 :
1275 : // If creating heavy quark by Q -> gamma + Q then next g -> Q + Qbar.
1276 : // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
1277 0 : if ( abs(idMother) == 4 || abs(idMother) == 5 ) {
1278 0 : double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
1279 0 : double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
1280 0 : if (pT2QQcorr < TINYPT2) { wt = 0.; continue; }
1281 0 : }
1282 :
1283 : // Optional dampening of large pT values in first radiation.
1284 0 : if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
1285 0 : wt *= pT2damp / (pT2 + pT2damp);
1286 :
1287 : // Evaluation of new daughter PDF
1288 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1289 0 : double xPDFdaughterNew = beam.xfISR(iSysNow, idDaughter, xDaughter,
1290 : pdfScale2);
1291 0 : if (xPDFdaughterNew < TINYPDF) {
1292 : xPDFdaughterNew = TINYPDF;
1293 0 : }
1294 :
1295 : // Evaluation of new charge-weighted mother PDF
1296 0 : double xPDFmotherNew = pow2( (abs(idMother+1) % 2 + 1)/3.0 )
1297 0 : * beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
1298 :
1299 : // Trial weight: divide out old pdf ratio
1300 0 : wt *= xPDFdaughter / xPDFmother[idMother + 10];
1301 :
1302 : // Trial weight: new pdf ratio
1303 0 : wt *= xPDFmotherNew / xPDFdaughterNew;
1304 :
1305 : // Iterate until acceptable pT (or have fallen below pTmin).
1306 0 : } while (wt < rndmPtr->flat()) ;
1307 0 : }
1308 :
1309 : // Save values for (so far) acceptable branching.
1310 0 : dipEndNow->store( idDaughter, idMother, idSister, x1Now, x2Now, m2Dip,
1311 0 : pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
1312 :
1313 0 : }
1314 :
1315 : //--------------------------------------------------------------------------
1316 :
1317 : void SpaceShower::pT2nextWeak( double pT2begDip, double pT2endDip) {
1318 :
1319 : // Type of dipole and starting values.
1320 0 : BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
1321 0 : bool isLeptonBeam = beam.isLepton();
1322 0 : bool isValence = beam[iSysNow].isValence();
1323 0 : int MEtype = dipEndNow->MEtype;
1324 : double pT2 = pT2begDip;
1325 0 : double m2Lepton = (isLeptonBeam) ? pow2(beam.m()) : 0.;
1326 0 : if (isLeptonBeam && pT2begDip < m2Lepton) return;
1327 :
1328 : // alpha_em at maximum scale provides upper estimate.
1329 0 : double alphaEMmax = alphaEM.alphaEM(pT2begDip);
1330 0 : double alphaEM2pi = alphaEMmax / (2. * M_PI);
1331 :
1332 : // Maximum x of mother implies minimum z = xDaughter / xMother.
1333 0 : double xMaxAbs = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
1334 0 : double zMinAbs = xDaughter / xMaxAbs;
1335 0 : if (xMaxAbs < 0.) {
1336 0 : infoPtr->errorMsg("Warning in SpaceShower::pT2nextWeak: "
1337 : "xMaxAbs negative");
1338 0 : return;
1339 : }
1340 :
1341 : // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
1342 0 : double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
1343 0 : ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
1344 0 : if (isLeptonBeam) {
1345 0 : double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
1346 0 : if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
1347 0 : }
1348 0 : if (zMaxAbs < zMinAbs) return;
1349 :
1350 : // Variables used inside evolution loop. (Mainly dummy start values.)
1351 : int idMother = 0;
1352 : int idSister = 0;
1353 : // Check whether emission of W+, W- or Z0.
1354 0 : if (dipEndNow->weakType == 1) {
1355 0 : idSister = (idDaughter > 0) ? -24 : 24;
1356 0 : if (abs(idDaughter) % 2 == 1) idSister = -idSister;
1357 0 : } else if (dipEndNow->weakType == 2) idSister = 23;
1358 : double z = 0.;
1359 : double xMother = 0.;
1360 : double wt = 0.;
1361 : double Q2 = 0.;
1362 0 : double mSister = particleDataPtr->mSel(idSister);
1363 0 : double m2Sister = pow2(mSister);
1364 : double pT2corr = 0.;
1365 :
1366 : // Find maximum z due to massive sister.
1367 : // Still need to prove that this actually is an overestimate.
1368 0 : double mRatio = mSister/ sqrt(m2Dip);
1369 0 : double m2R1 = 1. + pow2(mRatio);
1370 0 : double zMaxMassive = 1. / (m2R1 + pT2endDip/m2Dip);
1371 0 : if (zMaxMassive < zMaxAbs) zMaxAbs = zMaxMassive;
1372 0 : if (zMaxAbs < zMinAbs) return;
1373 :
1374 : // Weak evolution of fermions.
1375 : // Integrals of splitting kernels for fermions: f -> f.
1376 : // Old ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
1377 : // New Ansatz f(z) = (1 + (1+r^2)^2) / (1 - z * (1 + r^2)).
1378 : // This should always be an overestimate for massive emissions.
1379 : // Not yet implemented correctly for lepton beam.
1380 : double f2fInt = 0.;
1381 0 : double f2fIntA = (1. + pow2(zMaxAbs * m2R1)) / m2R1
1382 0 : * log( (1. - m2R1 * zMinAbs) / (1. - m2R1 * zMaxAbs) );
1383 : double f2fIntB = 0.;
1384 0 : double zRootMin = (1. + sqrt(m2R1 * zMinAbs)) / (1. - sqrt(m2R1 * zMinAbs));
1385 0 : double zRootMax = (1. + sqrt(m2R1 * zMaxAbs)) / (1. - sqrt(m2R1 * zMaxAbs));
1386 0 : if (isLeptonBeam) {
1387 0 : f2fIntB = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
1388 0 : f2fInt = f2fIntA + f2fIntB;
1389 0 : } else if (isValence)
1390 0 : f2fInt = (1. + pow2(zMaxAbs) * pow2(m2R1))/ sqrt(m2R1)
1391 0 : * log( zRootMax / zRootMin );
1392 : else
1393 : f2fInt = f2fIntA;
1394 :
1395 : // Calculate the weak coupling: separate for left- and right-handed fermions.
1396 : double weakCoupling = 0;
1397 0 : if (dipEndNow->weakType == 1)
1398 0 : weakCoupling = 2. * alphaEM2pi / (4. * coupSMPtr->sin2thetaW());
1399 0 : else if (dipEndNow->weakType == 2 && dipEndNow->weakPol == -1)
1400 0 : weakCoupling = alphaEM2pi * thetaWRat
1401 0 : * pow2(2. * coupSMPtr->lf( abs(idDaughter) ));
1402 : else
1403 0 : weakCoupling = alphaEM2pi * thetaWRat
1404 0 : * pow2(2. * coupSMPtr->rf(abs(idDaughter) ));
1405 :
1406 : // Find the possible mothers for a W emission. So far quarks only.
1407 0 : vector<int> possibleMothers;
1408 0 : if (abs(idDaughter) % 2 == 0) {
1409 0 : possibleMothers.push_back(1);
1410 0 : possibleMothers.push_back(3);
1411 0 : possibleMothers.push_back(5);
1412 : } else {
1413 0 : possibleMothers.push_back(2);
1414 0 : possibleMothers.push_back(4);
1415 : }
1416 0 : if (idDaughter < 0)
1417 0 : for (unsigned int i = 0;i < possibleMothers.size();i++)
1418 0 : possibleMothers[i] = - possibleMothers[i];
1419 :
1420 : // Check if daughter estimate is 0, return in that case.
1421 : // Only write warning if u, d or g is the daughter.
1422 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2begDip;
1423 0 : double xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2);
1424 0 : if (xPDFdaughter < TINYPDF) {
1425 0 : if (abs(idDaughter) == 1 || abs(idDaughter) == 2 || abs(idDaughter) == 21)
1426 0 : infoPtr->errorMsg("Warning in SpaceShower::pT2nextWeak: "
1427 : "very small PDF");
1428 0 : return;
1429 : }
1430 :
1431 : // PDF and CKM upper estimate needed for W emission.
1432 : double overEstimatePDFCKM = 0.;
1433 0 : if (dipEndNow->weakType == 1) {
1434 0 : for (unsigned int i = 0; i < possibleMothers.size(); i++)
1435 0 : overEstimatePDFCKM += coupSMPtr->V2CKMid(idDaughter, possibleMothers[i])
1436 0 : * beam.xfISR(iSysNow, possibleMothers[i], xDaughter, pdfScale2)
1437 0 : / xPDFdaughter;
1438 0 : }
1439 0 : if (dipEndNow->weakType == 2) overEstimatePDFCKM = 1.;
1440 :
1441 : // Upper estimate for evolution equation, including fudge factor.
1442 0 : if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
1443 0 : double kernelPDF = weakCoupling * f2fInt * weakEnhancement;
1444 :
1445 : // PDF and CKM overestimate.
1446 0 : kernelPDF *= overEstimatePDFCKM;
1447 0 : double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
1448 0 : kernelPDF *= fudge;
1449 0 : if (kernelPDF < TINYKERNELPDF) return;
1450 :
1451 : // Begin evolution loop towards smaller pT values.
1452 : do {
1453 :
1454 : // Pick pT2 (in overestimated z range).
1455 : // For l -> l gamma include extrafactor 1 / ln(pT2 / m2l) in evolution.
1456 0 : double shift = pow(rndmPtr->flat(), 1. / kernelPDF);
1457 0 : if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
1458 0 : else pT2 = pT2 * shift;
1459 :
1460 : // Abort evolution if below cutoff scale, or below another branching.
1461 0 : if (pT2 < pT2endDip) return;
1462 0 : if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
1463 :
1464 : // Abort evolution if below mass treshold.
1465 0 : if (pT2 < HEAVYPT2EVOL * pow2(particleDataPtr->m0(idDaughter))) return;
1466 :
1467 : // Set the id for the mother particle to be equal to the daughter
1468 : // particle. This is correct for Z, and it will later be changed for W.
1469 0 : idMother = idDaughter;
1470 :
1471 : // Select z value of branching f -> f + Z/W, and corrective weight.
1472 : wt = 1.0;
1473 0 : if (isLeptonBeam) {
1474 0 : if (f2fIntA > rndmPtr->flat() * (f2fIntA + f2fIntB)) {
1475 0 : z = 1. - (1. - zMinAbs)
1476 0 : * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
1477 0 : } else {
1478 0 : z = xDaughter + (zMinAbs - xDaughter)
1479 0 : * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter),
1480 0 : rndmPtr->flat() );
1481 : }
1482 0 : wt *= (z - xDaughter) / (1. - xDaughter);
1483 0 : } else if (isValence) {
1484 : // Valence more peaked at large z.
1485 0 : double zTmp = zRootMin * pow(zRootMax / zRootMin, rndmPtr->flat() );
1486 0 : z = pow2( (1. - zTmp) / (1. + zTmp) ) / m2R1;
1487 0 : wt *= sqrt(z);
1488 0 : } else {
1489 0 : z = (1. - (1. - zMinAbs * m2R1) * pow( (1. - zMaxAbs * m2R1)
1490 0 : / (1. - zMinAbs * m2R1), rndmPtr->flat() ) ) / m2R1;
1491 : }
1492 0 : wt *= (1. + pow2(z * m2R1)) / (1. + pow2(zMaxAbs * m2R1));
1493 :
1494 : // Derive Q2 and x of mother from pT2 and z.
1495 0 : Q2 = pT2 / (1. - z);
1496 0 : xMother = xDaughter / z;
1497 :
1498 : // Correction to x for massive recoiler from rescattering.
1499 0 : if (!dipEndNow->normalRecoil) {
1500 0 : if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1501 0 : else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1502 : }
1503 0 : if (xMother > xMaxAbs) { wt = 0.; continue; }
1504 :
1505 : // Forbidden emission if outside allowed z range for given pT2.
1506 0 : pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1507 0 : if (pT2corr < TINYPT2) { wt = 0.; continue; }
1508 :
1509 : // Correct by ln(pT2 / m2l) and fudge factor.
1510 0 : if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
1511 :
1512 : // Evaluation of ME correction.
1513 0 : if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
1514 0 : m2Dip, z, Q2, m2Sister) / calcMEmax(MEtype, idMother, idDaughter);
1515 :
1516 : // Optional dampening of large pT values in first radiation.
1517 : // Allow damping also for corrected matrix elements
1518 0 : if (dopTdamp && iSysNow == 0 && nRad == 0)
1519 0 : wt *= pT2damp / (pT2 + pT2damp);
1520 :
1521 : // Correct to current value of alpha_EM.
1522 0 : double alphaEMnow = alphaEM.alphaEM(pT2);
1523 0 : wt *= (alphaEMnow / alphaEMmax);
1524 :
1525 : // Evaluation of new daughter and mother PDF's for Z.
1526 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1527 0 : double xPDFdaughterNew = max ( TINYPDF,
1528 0 : beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
1529 0 : if (dipEndNow->weakType == 2) {
1530 : double xPDFmotherNew
1531 0 : = beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
1532 0 : wt *= xPDFmotherNew / xPDFdaughterNew;
1533 0 : }
1534 :
1535 : // Evaluation of daughter and mother PDF's for W.
1536 0 : if (dipEndNow->weakType == 1) {
1537 0 : double valPDFCKM[3] = {0.};
1538 : double sumPDFCKM = 0.;
1539 0 : for (unsigned int i = 0; i < possibleMothers.size(); i++) {
1540 0 : valPDFCKM[i] = coupSMPtr->V2CKMid(idDaughter,possibleMothers[i])
1541 0 : * beam.xfISR(iSysNow, possibleMothers[i], xMother, pdfScale2)
1542 0 : / xPDFdaughterNew;
1543 0 : sumPDFCKM += valPDFCKM[i];
1544 : }
1545 0 : wt *= sumPDFCKM / overEstimatePDFCKM;
1546 :
1547 : // Choose id for mother in case of W.
1548 0 : double pickId = sumPDFCKM * rndmPtr->flat();
1549 : int iId = -1;
1550 0 : int pmSize = possibleMothers.size();
1551 0 : do pickId -= valPDFCKM[++iId];
1552 0 : while (pickId > 0. && iId < pmSize);
1553 0 : idMother = possibleMothers[iId];
1554 0 : }
1555 :
1556 : // Warn if too big weight.
1557 0 : if (wt > 1.) infoPtr->errorMsg("Warning in SpaceShower::pT2nextWeak: "
1558 : "weight is above unity.");
1559 :
1560 : // Iterate until acceptable pT (or have fallen below pTmin).
1561 0 : } while (wt < rndmPtr->flat()) ;
1562 :
1563 : // Save values for (so far) acceptable branching.
1564 0 : dipEndNow->store( idDaughter, idMother, idSister, x1Now, x2Now, m2Dip,
1565 0 : pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
1566 :
1567 0 : }
1568 :
1569 : //--------------------------------------------------------------------------
1570 :
1571 : // Kinematics of branching.
1572 : // Construct mother -> daughter + sister, with recoiler on other side.
1573 :
1574 : bool SpaceShower::branch( Event& event) {
1575 :
1576 : // Side on which branching occured.
1577 0 : int side = abs(dipEndSel->side);
1578 0 : double sideSign = (side == 1) ? 1. : -1.;
1579 :
1580 : // Read in flavour and colour variables.
1581 0 : int iDaughter = partonSystemsPtr->getInA(iSysSel);
1582 0 : int iRecoiler = partonSystemsPtr->getInB(iSysSel);
1583 0 : if (side == 2) swap(iDaughter, iRecoiler);
1584 0 : int idDaughterNow = dipEndSel->idDaughter;
1585 0 : int idMother = dipEndSel->idMother;
1586 0 : int idSister = dipEndSel->idSister;
1587 0 : int idRecoiler = event[iRecoiler].id();
1588 0 : int colDaughter = event[iDaughter].col();
1589 0 : int acolDaughter = event[iDaughter].acol();
1590 :
1591 : // Recoil parton may be rescatterer, requiring special processing.
1592 0 : bool normalRecoil = dipEndSel->normalRecoil;
1593 0 : int iRecoilMother = event[iRecoiler].mother1();
1594 :
1595 : // Read in kinematical variables.
1596 0 : double x1 = dipEndSel->x1;
1597 0 : double x2 = dipEndSel->x2;
1598 0 : double xMo = dipEndSel->xMo;
1599 0 : double m2 = dipEndSel->m2Dip;
1600 0 : double m = sqrt(m2);
1601 0 : double pT2 = dipEndSel->pT2;
1602 0 : double z = dipEndSel->z;
1603 0 : double Q2 = dipEndSel->Q2;
1604 0 : double mSister = dipEndSel->mSister;
1605 0 : double m2Sister = dipEndSel->m2Sister;
1606 0 : double pT2corr = dipEndSel->pT2corr;
1607 0 : double x1New = (side == 1) ? xMo : x1;
1608 0 : double x2New = (side == 2) ? xMo : x2;
1609 :
1610 : // Read in MEtype:
1611 0 : int MEtype = dipEndSel->MEtype;
1612 :
1613 : // Rescatter: kinematics may fail; use the rescatterFail flag to tell
1614 : // parton level to try again.
1615 0 : rescatterFail = false;
1616 :
1617 : // Construct kinematics of mother, sister and recoiler in old rest frame.
1618 : // Normally both mother and recoiler are taken massless.
1619 0 : double eNewRec, pzNewRec, pTbranch, pzMother;
1620 0 : if (normalRecoil) {
1621 0 : eNewRec = 0.5 * (m2 + Q2) / m;
1622 0 : pzNewRec = -sideSign * eNewRec;
1623 0 : pTbranch = sqrt(pT2corr) * m2 / ( z * (m2 + Q2) );
1624 0 : pzMother = sideSign * 0.5 * m * ( (m2 - Q2) / ( z * (m2 + Q2) )
1625 0 : + (Q2 + m2Sister) / m2 );
1626 : // More complicated kinematics when recoiler not massless. May fail.
1627 0 : } else {
1628 0 : m2Rec = event[iRecoiler].m2();
1629 0 : double s1Tmp = m2 + Q2 - m2Rec;
1630 0 : double s3Tmp = m2 / z - m2Rec;
1631 0 : double r1Tmp = sqrt(s1Tmp * s1Tmp + 4. * Q2 * m2Rec);
1632 0 : eNewRec = 0.5 * (m2 + m2Rec + Q2) / m;
1633 0 : pzNewRec = -sideSign * 0.5 * r1Tmp / m;
1634 0 : double pT2br = Q2 * s3Tmp * (m2 / z - m2 - Q2)
1635 0 : - m2Sister * s1Tmp * s3Tmp - m2Rec * pow2(Q2 + m2Sister);
1636 0 : if (pT2br <= 0.) return false;
1637 0 : pTbranch = sqrt(pT2br) / r1Tmp;
1638 0 : pzMother = sideSign * (m * s3Tmp
1639 0 : - eNewRec * (m2 / z - Q2 - m2Rec - m2Sister)) / r1Tmp;
1640 0 : }
1641 : // Common final kinematics steps for both normal and rescattering.
1642 0 : double eMother = sqrt( pow2(pTbranch) + pow2(pzMother) );
1643 0 : double pzSister = pzMother + pzNewRec;
1644 0 : double eSister = sqrt( pow2(pTbranch) + pow2(pzSister) + m2Sister );
1645 0 : Vec4 pMother( pTbranch, 0., pzMother, eMother );
1646 0 : Vec4 pSister( pTbranch, 0., pzSister, eSister );
1647 0 : Vec4 pNewRec( 0., 0., pzNewRec, eNewRec );
1648 :
1649 : // Current event and subsystem size.
1650 0 : int eventSizeOld = event.size();
1651 0 : int systemSizeOld = partonSystemsPtr->sizeAll(iSysSel);
1652 :
1653 : // Save properties to be restored in case of user-hook veto of emission.
1654 0 : int beamOff1 = 1 + beamOffset;
1655 0 : int beamOff2 = 2 + beamOffset;
1656 0 : int ev1Dau1V = event[beamOff1].daughter1();
1657 0 : int ev2Dau1V = event[beamOff2].daughter1();
1658 0 : vector<int> statusV, mother1V, mother2V, daughter1V, daughter2V;
1659 :
1660 : // Check if the first emission shoild be checked for removal
1661 0 : bool canMergeFirst = (mergingHooksPtr != 0)
1662 0 : ? mergingHooksPtr->canVetoEmission() : false;
1663 0 : if (canVetoEmission || canMergeFirst || doWeakShower) {
1664 0 : for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1665 0 : int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
1666 0 : statusV.push_back( event[iOldCopy].status());
1667 0 : mother1V.push_back( event[iOldCopy].mother1());
1668 0 : mother2V.push_back( event[iOldCopy].mother2());
1669 0 : daughter1V.push_back( event[iOldCopy].daughter1());
1670 0 : daughter2V.push_back( event[iOldCopy].daughter2());
1671 : }
1672 0 : }
1673 :
1674 : // Take copy of existing system, to be given modified kinematics.
1675 : // Incoming negative status. Rescattered also negative, but after copy.
1676 0 : for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1677 0 : int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
1678 0 : int statusOld = event[iOldCopy].status();
1679 0 : int statusNew = (iOldCopy == iDaughter
1680 0 : || iOldCopy == iRecoiler) ? statusOld : 44;
1681 0 : int iNewCopy = event.copy(iOldCopy, statusNew);
1682 0 : if (statusOld < 0) event[iNewCopy].statusNeg();
1683 : }
1684 :
1685 : // Define colour flow in branching.
1686 : // Default corresponds to f -> f + gamma.
1687 : int colMother = colDaughter;
1688 : int acolMother = acolDaughter;
1689 : int colSister = 0;
1690 : int acolSister = 0;
1691 0 : if (idSister == 22 || idSister == 23 || abs(idSister) == 24) ;
1692 : // q -> q + g and 50% of g -> g + g; need new colour.
1693 0 : else if (idSister == 21 && ( (idMother > 0 && idMother < 9)
1694 0 : || (idMother == 21 && rndmPtr->flat() < 0.5) ) ) {
1695 0 : colMother = event.nextColTag();
1696 : colSister = colMother;
1697 : acolSister = colDaughter;
1698 : // qbar -> qbar + g and other 50% of g -> g + g; need new colour.
1699 0 : } else if (idSister == 21) {
1700 0 : acolMother = event.nextColTag();
1701 : acolSister = acolMother;
1702 : colSister = acolDaughter;
1703 : // q -> g + q.
1704 0 : } else if (idDaughterNow == 21 && idMother > 0) {
1705 : colMother = colDaughter;
1706 : acolMother = 0;
1707 : colSister = acolDaughter;
1708 : // qbar -> g + qbar
1709 0 : } else if (idDaughterNow == 21) {
1710 : acolMother = acolDaughter;
1711 : colMother = 0;
1712 : acolSister = colDaughter;
1713 : // g -> q + qbar.
1714 0 : } else if (idDaughterNow > 0 && idDaughterNow < 9) {
1715 0 : acolMother = event.nextColTag();
1716 : acolSister = acolMother;
1717 : // g -> qbar + q.
1718 0 : } else if (idDaughterNow < 0 && idDaughterNow > -9) {
1719 0 : colMother = event.nextColTag();
1720 : colSister = colMother;
1721 : // q -> gamma + q.
1722 0 : } else if (idDaughterNow == 22 && idMother > 0) {
1723 0 : colMother = event.nextColTag();
1724 : colSister = colMother;
1725 : // qbar -> gamma + qbar.
1726 0 : } else if (idDaughterNow == 22) {
1727 0 : acolMother = event.nextColTag();
1728 : acolSister = acolMother;
1729 0 : }
1730 :
1731 : // Indices of partons involved. Add new sister.
1732 0 : int iMother = eventSizeOld + side - 1;
1733 0 : int iNewRecoiler = eventSizeOld + 2 - side;
1734 0 : int iSister = event.append( idSister, 43, iMother, 0, 0, 0,
1735 0 : colSister, acolSister, pSister, mSister, sqrt(pT2) );
1736 :
1737 : // References to the partons involved.
1738 0 : Particle& daughter = event[iDaughter];
1739 0 : Particle& mother = event[iMother];
1740 0 : Particle& newRecoiler = event[iNewRecoiler];
1741 0 : Particle& sister = event.back();
1742 :
1743 : // Replace old by new mother; update new recoiler.
1744 0 : mother.id( idMother );
1745 0 : mother.status( -41);
1746 0 : mother.cols( colMother, acolMother);
1747 0 : mother.p( pMother);
1748 0 : if (mother.idAbs() == 21 || mother.idAbs() == 22) mother.pol(9);
1749 0 : newRecoiler.status( (normalRecoil) ? -42 : -46 );
1750 0 : newRecoiler.p( pNewRec);
1751 0 : if (!normalRecoil) newRecoiler.m( event[iRecoiler].m() );
1752 :
1753 : // Update mother and daughter pointers; also for beams.
1754 0 : daughter.mothers( iMother, 0);
1755 0 : mother.daughters( iSister, iDaughter);
1756 0 : if (iSysSel == 0) {
1757 0 : event[beamOff1].daughter1( (side == 1) ? iMother : iNewRecoiler );
1758 0 : event[beamOff2].daughter1( (side == 2) ? iMother : iNewRecoiler );
1759 0 : }
1760 :
1761 : // Special checks to set weak particles status equal to 47.
1762 0 : if (sister.idAbs() == 23 || sister.idAbs() == 24) sister.status(47);
1763 :
1764 : // Find boost to old rest frame.
1765 0 : RotBstMatrix Mtot;
1766 0 : if (normalRecoil) Mtot.bst(0., 0., (x2 - x1) / (x1 + x2) );
1767 0 : else if (side == 1)
1768 0 : Mtot.toCMframe( event[iDaughter].p(), event[iRecoiler].p() );
1769 0 : else Mtot.toCMframe( event[iRecoiler].p(), event[iDaughter].p() );
1770 :
1771 : // Initially select phi angle of branching at random.
1772 0 : double phi = 2. * M_PI * rndmPtr->flat();
1773 :
1774 : // Evaluate coefficient of azimuthal asymmetry from gluon polarization.
1775 0 : findAsymPol( event, dipEndSel);
1776 0 : int iFinPol = dipEndSel->iFinPol;
1777 0 : double cPol = dipEndSel->asymPol;
1778 0 : double phiPol = (iFinPol == 0) ? 0. : event[iFinPol].phi();
1779 :
1780 : // If interference: try to match sister (anti)colour to final state.
1781 : int iFinInt = 0;
1782 0 : double cInt = 0.;
1783 : double phiInt = 0.;
1784 0 : if (doPhiIntAsym) {
1785 0 : for (int i = 0; i < partonSystemsPtr->sizeOut(iSysSel); ++ i) {
1786 0 : int iOut = partonSystemsPtr->getOut(iSysSel, i);
1787 0 : if ( (acolSister != 0 && event[iOut].col() == acolSister)
1788 0 : || (colSister != 0 && event[iOut].acol() == colSister) )
1789 0 : iFinInt = iOut;
1790 : }
1791 0 : if (iFinInt != 0) {
1792 : // Boost final-state parton to current frame of new kinematics.
1793 0 : Vec4 pFinTmp = event[iFinInt].p();
1794 0 : pFinTmp.rotbst(Mtot);
1795 0 : double theFin = pFinTmp.theta();
1796 0 : if (side == 2) theFin = M_PI - theFin;
1797 0 : double theSis = pSister.theta();
1798 0 : if (side == 2) theSis = M_PI - theSis;
1799 0 : cInt = strengthIntAsym * 2. * theSis * theFin
1800 0 : / (pow2(theSis) + pow2(theFin));
1801 0 : phiInt = event[iFinInt].phi();
1802 0 : }
1803 : }
1804 :
1805 : // Bias phi distribution for polarization and interference.
1806 0 : if (iFinPol != 0 || iFinInt != 0) {
1807 0 : double cPhiPol, cPhiInt, weight;
1808 0 : do {
1809 0 : phi = 2. * M_PI * rndmPtr->flat();
1810 : weight = 1.;
1811 0 : if (iFinPol !=0 ) {
1812 0 : cPhiPol = cos(phi - phiPol);
1813 0 : weight *= ( 1. + cPol * (2. * pow2(cPhiPol) - 1.) )
1814 0 : / ( 1. + abs(cPol) );
1815 0 : }
1816 0 : if (iFinInt !=0 ) {
1817 0 : cPhiInt = cos(phi - phiInt);
1818 0 : weight *= (1. - cInt) * (1. - cInt * cPhiInt)
1819 0 : / (1. + pow2(cInt) - 2. * cInt * cPhiInt);
1820 0 : }
1821 0 : } while (weight < rndmPtr->flat());
1822 0 : }
1823 :
1824 : // Include rotation -phi on boost to old rest frame.
1825 0 : Mtot.rot(0., -phi);
1826 :
1827 : // Find boost from old rest frame to event cm frame.
1828 0 : RotBstMatrix MfromRest;
1829 : // The boost to the new rest frame.
1830 0 : Vec4 sumNew = pMother + pNewRec;
1831 0 : double betaX = sumNew.px() / sumNew.e();
1832 0 : double betaZ = sumNew.pz() / sumNew.e();
1833 0 : MfromRest.bst( -betaX, 0., -betaZ);
1834 : // Alignment of radiator + recoiler to +- z axis, and rotation +phi.
1835 : // Note: with spacelike (E < 0) recoiler p'_x_mother < 0 can happen!
1836 0 : pMother.rotbst(MfromRest);
1837 0 : double theta = pMother.theta();
1838 0 : if (pMother.px() < 0.) theta = -theta;
1839 0 : if (side == 2) theta += M_PI;
1840 0 : MfromRest.rot(-theta, phi);
1841 : // Boost to radiator + recoiler in event cm frame.
1842 0 : if (normalRecoil) {
1843 0 : MfromRest.bst( 0., 0., (x1New - x2New) / (x1New + x2New) );
1844 0 : } else if (side == 1) {
1845 0 : Vec4 pMotherWanted( 0., 0., 0.5 * eCM * x1New, 0.5 * eCM * x1New);
1846 0 : MfromRest.fromCMframe( pMotherWanted, event[iRecoiler].p() );
1847 :
1848 0 : } else {
1849 0 : Vec4 pMotherWanted( 0., 0., -0.5 * eCM * x2New, 0.5 * eCM * x2New);
1850 0 : MfromRest.fromCMframe( event[iRecoiler].p(), pMotherWanted );
1851 0 : }
1852 0 : Mtot.rotbst(MfromRest);
1853 :
1854 : // ME correction for weak emissions in the t-channel.
1855 : double wt;
1856 0 : if (MEtype == 201 || MEtype == 202 || MEtype == 203 ||
1857 0 : MEtype == 206 || MEtype == 207 || MEtype == 208) {
1858 :
1859 : // Start by finding the correct outgoing particles
1860 : // to use in the ME correction.
1861 0 : Vec4 pA0 = mother.p();
1862 0 : Vec4 pB = newRecoiler.p();
1863 0 : bool sideRad = (abs(side) == 1);
1864 0 : Vec4 p1 = event[5].p();
1865 0 : Vec4 p2 = event[6].p();
1866 0 : int id1 = event[5].id();
1867 0 : int id2 = event[6].id();
1868 0 : if (!tChannel) {swap(p1,p2); swap(id1,id2);}
1869 0 : if (!sideRad) {swap(p1,p2); swap(id1,id2);}
1870 :
1871 : // Rotate with -phi to keep correct for the later +phi rotation.
1872 0 : p1.rot(0., -phi);
1873 0 : p2.rot(0., -phi);
1874 :
1875 : // Calculate the actual weight.
1876 0 : wt = calcMEcorrWeak(MEtype, m2, z, pT2, pA0, pB,
1877 0 : event[iDaughter].p(), event[iRecoiler].p(), p1, p2, sister.p());
1878 0 : if (wt > weakMaxWt) {
1879 0 : weakMaxWt = wt;
1880 0 : infoPtr->errorMsg("Warning in SpaceShower::Branch: "
1881 : "weight is above unity for weak emission.");
1882 0 : }
1883 :
1884 : // If weighting fails then restore event record to state before emission.
1885 0 : if (wt < rndmPtr->flat()) {
1886 0 : event.popBack( event.size() - eventSizeOld);
1887 0 : event[beamOff1].daughter1( ev1Dau1V);
1888 0 : event[beamOff2].daughter1( ev2Dau1V);
1889 0 : for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1890 0 : int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
1891 0 : event[iOldCopy].status( statusV[iCopy]);
1892 0 : event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
1893 0 : event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
1894 : }
1895 0 : return false;
1896 : }
1897 0 : }
1898 :
1899 : // Perform cumulative rotation/boost operation.
1900 : // Mother, recoiler and sister from old rest frame to event cm frame.
1901 0 : mother.rotbst(MfromRest);
1902 0 : newRecoiler.rotbst(MfromRest);
1903 0 : sister.rotbst(MfromRest);
1904 : // The rest from (and to) event cm frame.
1905 0 : for ( int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i)
1906 0 : event[i].rotbst(Mtot);
1907 :
1908 : // Remove double counting. Only implemented for QCD hard processes
1909 : // and for the first emission.
1910 0 : if (infoPtr->nISR() + infoPtr->nFSRinProc() == 0
1911 0 : && infoPtr->code() > 110 && infoPtr->code() < 130
1912 0 : && MEtype >= 200 && MEtype < 210 && vetoWeakJets) {
1913 :
1914 : // Find outgoing particles.
1915 0 : int iP1 = event[5].daughter1();
1916 0 : int iP2 = event[6].daughter1();
1917 0 : Vec4 pP1 = event[iP1].p();
1918 0 : Vec4 pP2 = event[iP2].p();
1919 :
1920 : // Set start pT2 as pT2 of emitted particle and therefore no cut.
1921 0 : double d = sister.pT2();
1922 : bool cut = false;
1923 :
1924 0 : if (pP1.pT2() < d) {d = pP1.pT2(); cut = true;}
1925 0 : if (pP2.pT2() < d) {d = pP2.pT2(); cut = true;}
1926 :
1927 : // Check for angle between weak boson and quarks
1928 : // (require final state particle to be a fermion).
1929 0 : if (event[iP1].idAbs() < 20) {
1930 0 : double dij = min(pP1.pT2(),sister.pT2())
1931 0 : * pow2(RRapPhi(pP1,sister.p()))/vetoWeakDeltaR2;
1932 0 : if (dij < d) {
1933 : d = dij;
1934 : cut = false;
1935 0 : }
1936 0 : }
1937 :
1938 0 : if (event[iP2].idAbs() < 20) {
1939 0 : double dij = min(pP2.pT2(),sister.pT2())
1940 0 : * pow2(RRapPhi(pP2,sister.p()))/vetoWeakDeltaR2;
1941 0 : if (dij < d) {
1942 : d = dij;
1943 : cut = false;
1944 0 : }
1945 0 : }
1946 :
1947 : // Check for angle between recoiler and radiator, if quark anti-quark pair,
1948 : // or if the recoiler is a gluon.
1949 0 : if (event[iP1].idAbs() == 21 || event[iP2].idAbs() == 21 ||
1950 0 : event[iP1].id() == - event[iP2].id()) {
1951 0 : double dij = min(pP1.pT2(),pP2.pT2())
1952 0 : * pow2(RRapPhi(pP1,pP2))/vetoWeakDeltaR2;
1953 0 : if (dij < d) {
1954 : d = dij;
1955 : cut = true;
1956 0 : }
1957 0 : }
1958 :
1959 : // Clean up event if the emission should be removed.
1960 0 : if (cut) {
1961 0 : event.popBack( event.size() - eventSizeOld);
1962 0 : event[beamOff1].daughter1( ev1Dau1V);
1963 0 : event[beamOff2].daughter1( ev2Dau1V);
1964 0 : for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1965 0 : int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
1966 0 : event[iOldCopy].status( statusV[iCopy]);
1967 0 : event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
1968 0 : event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
1969 : }
1970 0 : return false;
1971 : }
1972 0 : }
1973 :
1974 : // Allow veto of branching. If so restore event record to before emission.
1975 0 : if ( (canVetoEmission
1976 0 : && userHooksPtr->doVetoISREmission(eventSizeOld, event, iSysSel))
1977 0 : || (canMergeFirst
1978 0 : && mergingHooksPtr->doVetoEmission( event )) ) {
1979 0 : event.popBack( event.size() - eventSizeOld);
1980 0 : event[beamOff1].daughter1( ev1Dau1V);
1981 0 : event[beamOff2].daughter1( ev2Dau1V);
1982 0 : for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1983 0 : int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
1984 0 : event[iOldCopy].status( statusV[iCopy]);
1985 0 : event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
1986 0 : event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
1987 : }
1988 0 : return false;
1989 : }
1990 :
1991 : // Update list of partons in system; adding newly produced one.
1992 0 : partonSystemsPtr->setInA(iSysSel, eventSizeOld);
1993 0 : partonSystemsPtr->setInB(iSysSel, eventSizeOld + 1);
1994 0 : for (int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
1995 0 : partonSystemsPtr->setOut(iSysSel, iCopy - 2, eventSizeOld + iCopy);
1996 0 : partonSystemsPtr->addOut(iSysSel, eventSizeOld + systemSizeOld);
1997 0 : partonSystemsPtr->setSHat(iSysSel, m2 / z);
1998 :
1999 : // Add dipoles for q -> g q, where the daughter is the gluon.
2000 0 : if (idDaughter == 21 && idMother != 21) {
2001 0 : if (doQEDshowerByQ) {
2002 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2003 0 : iNewRecoiler, pT2, 0, mother.chargeType(), 0, 0, normalRecoil) );
2004 : }
2005 0 : if (doWeakShower && iSysSel == 0) {
2006 : int MEtypeNew = 203;
2007 0 : if (idRecoiler == 21) MEtypeNew = 201;
2008 0 : if (idRecoiler == idMother) MEtypeNew = 202;
2009 : // If original was a Drell-Yan, keep as Drell-Yan.
2010 0 : if( event[3].id() == - event[4].id()) MEtypeNew = 200;
2011 0 : int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
2012 0 : event[iMother].pol(weakPol);
2013 0 : if ((weakMode == 0 || weakMode == 1) && weakPol == -1)
2014 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2015 : iNewRecoiler, pT2, 0, 0, 1, MEtypeNew, normalRecoil, weakPol) );
2016 0 : if (weakMode == 0 || weakMode == 2)
2017 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2018 0 : iNewRecoiler, pT2, 0, 0, 2, MEtypeNew + 5, normalRecoil, weakPol) );
2019 0 : }
2020 : }
2021 :
2022 : // Add dipoles for q -> q gamma, where the daughter is the gamma.
2023 0 : if (idDaughter == 22 && idMother != 22) {
2024 0 : if (doQCDshower && mother.colType() != 0) {
2025 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2026 0 : iNewRecoiler, pT2, mother.colType(), 0, 0, 0, normalRecoil) );
2027 : }
2028 0 : if (doQEDshowerByQ && mother.chargeType() != 3) {
2029 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2030 0 : iNewRecoiler, pT2, 0, mother.chargeType(), 0, 0, normalRecoil) );
2031 : }
2032 0 : if (doQEDshowerByL && mother.chargeType() == 3) {
2033 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2034 0 : iNewRecoiler, pT2, 0, mother.chargeType(), 0, 0, normalRecoil) );
2035 : }
2036 0 : if (doWeakShower && iSysSel == 0) {
2037 : int MEtypeNew = 203;
2038 0 : if (idRecoiler == 21) MEtypeNew = 201;
2039 0 : if (idRecoiler == idMother) MEtypeNew = 202;
2040 : // If original was a Drell-Yan, keep as Drell-Yan.
2041 0 : if( event[3].id() == - event[4].id()) MEtypeNew = 200;
2042 0 : int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
2043 0 : event[iMother].pol(weakPol);
2044 0 : if ((weakMode == 0 || weakMode == 1) && weakPol == -1)
2045 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2046 : iNewRecoiler, pT2, 0, 0, 1, MEtypeNew, normalRecoil, weakPol) );
2047 0 : if (weakMode == 0 || weakMode == 2)
2048 0 : dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2049 0 : iNewRecoiler, pT2, 0, 0, 2, MEtypeNew + 5, normalRecoil, weakPol) );
2050 0 : }
2051 : }
2052 :
2053 : // dipEnd array may have expanded and been moved, so regenerate dipEndSel.
2054 0 : dipEndSel = &dipEnd[iDipSel];
2055 :
2056 : // Set flag to tell that a weak emission has happened.
2057 0 : if (dipEndSel->weakType != 0) hasWeaklyRadiated = true;
2058 :
2059 : // Update list of QCD emissions in side A and B in given iSysSel
2060 : // This is used to veto jets in W/z events.
2061 0 : while (iSysSel >= int(nRadA.size()) || iSysSel >= int(nRadB.size())) {
2062 0 : nRadA.push_back(0);
2063 0 : nRadB.push_back(0);
2064 : }
2065 0 : if (dipEndSel->colType != 0 && side == 1) ++nRadA[iSysSel];
2066 0 : else if (dipEndSel->colType != 0) ++nRadB[iSysSel];
2067 :
2068 : // Update info on radiating dipole ends (QCD, QED or weak).
2069 0 : for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
2070 0 : if ( dipEnd[iDip].system == iSysSel) {
2071 0 : if (abs(dipEnd[iDip].side) == side) {
2072 0 : dipEnd[iDip].iRadiator = iMother;
2073 0 : dipEnd[iDip].iRecoiler = iNewRecoiler;
2074 0 : if (dipEnd[iDip].colType != 0)
2075 0 : dipEnd[iDip].colType = mother.colType();
2076 0 : else if (dipEnd[iDip].chgType != 0) {
2077 0 : dipEnd[iDip].chgType = 0;
2078 0 : if ( (mother.isQuark() && doQEDshowerByQ)
2079 0 : || (mother.isLepton() && doQEDshowerByL) )
2080 0 : dipEnd[iDip].chgType = mother.chargeType();
2081 : }
2082 0 : else if (dipEnd[iDip].weakType != 0) {
2083 : // Kill weak dipole if mother becomes gluon / photon.
2084 0 : if (!(mother.isLepton() || mother.isQuark()))
2085 0 : dipEnd[iDip].weakType = 0;
2086 0 : if (singleWeakEmission && hasWeaklyRadiated)
2087 0 : dipEnd[iDip].weakType = 0;
2088 : }
2089 :
2090 : // Kill ME corrections after first emission for everything
2091 : // but weak showers.
2092 0 : if (dipEnd[iDip].weakType == 0) dipEnd[iDip].MEtype = 0;
2093 :
2094 : // Update info on recoiling dipole ends (QCD or QED).
2095 : } else {
2096 0 : dipEnd[iDip].iRadiator = iNewRecoiler;
2097 0 : dipEnd[iDip].iRecoiler = iMother;
2098 : // Optionally also kill recoiler ME corrections after first emission.
2099 0 : if (!doMEafterFirst && dipEnd[iDip].weakType == 0)
2100 0 : dipEnd[iDip].MEtype = 0;
2101 : // Remove weak dipoles if we only want a single emission.
2102 0 : if (dipEnd[iDip].weakType != 0 && singleWeakEmission
2103 0 : && hasWeaklyRadiated) dipEnd[iDip].weakType = 0;
2104 : }
2105 : }
2106 :
2107 : // Set polarisation of mother for weak emissions.
2108 0 : if (dipEndSel->weakType != 0) mother.pol(dipEndSel->weakPol);
2109 :
2110 : // Update info on beam remnants.
2111 0 : BeamParticle& beamNow = (side == 1) ? *beamAPtr : *beamBPtr;
2112 0 : double xNew = (side == 1) ? x1New : x2New;
2113 0 : beamNow[iSysSel].update( iMother, idMother, xNew);
2114 : // Redo choice of companion kind whenever new flavour.
2115 0 : if (idMother != idDaughterNow) {
2116 0 : pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
2117 0 : beamNow.xfISR( iSysSel, idMother, xNew, pdfScale2);
2118 0 : beamNow.pickValSeaComp();
2119 : }
2120 0 : BeamParticle& beamRec = (side == 1) ? *beamBPtr : *beamAPtr;
2121 0 : beamRec[iSysSel].iPos( iNewRecoiler);
2122 :
2123 : // Store branching values of current dipole. (For rapidity ordering.)
2124 0 : ++dipEndSel->nBranch;
2125 0 : dipEndSel->pT2Old = pT2;
2126 0 : dipEndSel->zOld = z;
2127 :
2128 : // Update history if recoiler rescatters.
2129 0 : if (!normalRecoil)
2130 0 : event[iRecoilMother].daughters( iNewRecoiler, iNewRecoiler);
2131 :
2132 : // Start list of rescatterers that force changed kinematics.
2133 0 : vector<int> iRescatterer;
2134 0 : for ( int i = 0; i < systemSizeOld - 2; ++i) {
2135 0 : int iOutNew = partonSystemsPtr->getOut( iSysSel, i);
2136 0 : if (!event[iOutNew].isFinal()) iRescatterer.push_back(iOutNew);
2137 0 : }
2138 :
2139 : // Start iterate over list of such rescatterers.
2140 : int iRescNow = -1;
2141 0 : while (++iRescNow < int(iRescatterer.size())) {
2142 :
2143 : // Identify partons that induce or are affected by rescatter shift.
2144 : // In following Old is before change of kinematics, New after,
2145 : // Out scatterer in outstate and In in instate of another system.
2146 : // Daughter sequence is (iOutOld ->) iOutNew -> iInNew -> iInOld.
2147 0 : int iOutNew = iRescatterer[iRescNow];
2148 0 : int iInOld = event[iOutNew].daughter1();
2149 0 : int iSysResc = partonSystemsPtr->getSystemOf(iInOld, true);
2150 :
2151 : // Copy incoming partons of rescattered system and hook them up.
2152 0 : int iOldA = partonSystemsPtr->getInA(iSysResc);
2153 0 : int iOldB = partonSystemsPtr->getInB(iSysResc);
2154 0 : bool rescSideA = event[iOldA].isRescatteredIncoming();
2155 0 : int statusNewA = (rescSideA) ? -45 : -42;
2156 0 : int statusNewB = (rescSideA) ? -42 : -45;
2157 0 : int iNewA = event.copy(iOldA, statusNewA);
2158 0 : int iNewB = event.copy(iOldB, statusNewB);
2159 :
2160 : // Copy outgoing partons of rescattered system and hook them up.
2161 0 : int eventSize = event.size();
2162 0 : int sizeOutAB = partonSystemsPtr->sizeOut(iSysResc);
2163 0 : int iOldAB, statusOldAB, iNewAB;
2164 0 : for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB) {
2165 0 : iOldAB = partonSystemsPtr->getOut(iSysResc, iOutAB);
2166 0 : statusOldAB = event[iOldAB].status();
2167 0 : iNewAB = event.copy(iOldAB, 44);
2168 : // Status could be negative for parton that rescatters in its turn.
2169 0 : if (statusOldAB < 0) {
2170 0 : event[iNewAB].statusNeg();
2171 0 : iRescatterer.push_back(iNewAB);
2172 : }
2173 : }
2174 :
2175 : // Hook up new outgoing with new incoming parton.
2176 0 : int iInNew = (rescSideA) ? iNewA : iNewB;
2177 0 : event[iOutNew].daughters( iInNew, iInNew);
2178 0 : event[iInNew].mothers( iOutNew, iOutNew);
2179 :
2180 : // Rescale recoiling incoming parton for correct invariant mass.
2181 0 : event[iInNew].p( event[iOutNew].p() );
2182 0 : double momFac = (rescSideA)
2183 0 : ? event[iInOld].pPos() / event[iInNew].pPos()
2184 0 : : event[iInOld].pNeg() / event[iInNew].pNeg();
2185 0 : int iInRec = (rescSideA) ? iNewB : iNewA;
2186 :
2187 : // Rescatter: A previous boost may cause the light cone momentum of a
2188 : // rescattered parton to change sign. If this happens, tell
2189 : // parton level to try again.
2190 0 : if (momFac < 0.0) {
2191 0 : infoPtr->errorMsg("Warning in SpaceShower::branch: "
2192 : "change in lightcone momentum sign; retrying parton level");
2193 0 : rescatterFail = true;
2194 0 : return false;
2195 : }
2196 0 : event[iInRec].rescale4( momFac);
2197 :
2198 : // Boost outgoing partons to new frame of incoming.
2199 0 : RotBstMatrix MmodResc;
2200 0 : MmodResc.toCMframe( event[iOldA].p(), event[iOldB].p());
2201 0 : MmodResc.fromCMframe(event[iNewA].p(), event[iNewB].p());
2202 0 : for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB)
2203 0 : event[eventSize + iOutAB].rotbst(MmodResc);
2204 :
2205 : // Update list of partons in system.
2206 0 : partonSystemsPtr->setInA(iSysResc, iNewA);
2207 0 : partonSystemsPtr->setInB(iSysResc, iNewB);
2208 0 : for (int iCopy = 0; iCopy < sizeOutAB; ++iCopy)
2209 0 : partonSystemsPtr->setOut(iSysResc, iCopy, eventSize + iCopy);
2210 :
2211 : // Update info on radiating dipole ends (QCD or QED).
2212 0 : for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
2213 0 : if ( dipEnd[iDip].system == iSysResc) {
2214 0 : bool sideAnow = (abs(dipEnd[iDip].side) == 1);
2215 0 : dipEnd[iDip].iRadiator = (sideAnow) ? iNewA : iNewB;
2216 0 : dipEnd[iDip].iRecoiler = (sideAnow) ? iNewB : iNewA;
2217 0 : }
2218 :
2219 : // Update info on beam remnants.
2220 0 : BeamParticle& beamResc = (rescSideA) ? *beamAPtr : *beamBPtr;
2221 0 : beamResc[iSysResc].iPos( iInNew);
2222 0 : beamResc[iSysResc].p( event[iInNew].p() );
2223 0 : beamResc[iSysResc].scaleX( 1. / momFac );
2224 0 : BeamParticle& beamReco = (rescSideA) ? *beamBPtr : *beamAPtr;
2225 0 : beamReco[iSysResc].iPos( iInRec);
2226 0 : beamReco[iSysResc].scaleX( momFac);
2227 :
2228 : // End iterate over list of rescatterers.
2229 0 : }
2230 :
2231 : // Check that beam momentum not used up by rescattered-system boosts.
2232 0 : if (beamAPtr->xMax(-1) < 0.0 || beamBPtr->xMax(-1) < 0.0) {
2233 0 : infoPtr->errorMsg("Warning in SpaceShower::branch: "
2234 : "used up beam momentum; retrying parton level");
2235 0 : rescatterFail = true;
2236 0 : return false;
2237 : }
2238 :
2239 : // Done without any errors.
2240 0 : return true;
2241 :
2242 0 : }
2243 :
2244 : //--------------------------------------------------------------------------
2245 :
2246 : // Find class of ME correction.
2247 :
2248 : int SpaceShower::findMEtype( int iSys, Event& event, bool weakRadiation) {
2249 :
2250 : // Default values and no action.
2251 : int MEtype = 0;
2252 0 : if (!doMEcorrections) return MEtype;
2253 :
2254 : // Identify systems producing a single resonance.
2255 0 : if (partonSystemsPtr->sizeOut( iSys) == 1 && !weakRadiation) {
2256 0 : int idIn1 = event[partonSystemsPtr->getInA(iSys)].id();
2257 0 : int idIn2 = event[partonSystemsPtr->getInA(iSys)].id();
2258 0 : int idRes = event[partonSystemsPtr->getOut(iSys, 0)].id();
2259 0 : if (iSys == 0) idResFirst = abs(idRes);
2260 0 : if (iSys == 1) idResSecond = abs(idRes);
2261 :
2262 : // f + fbar -> vector boson.
2263 0 : if ( (idRes == 23 || abs(idRes) == 24 || idRes == 32
2264 0 : || idRes == 33 || abs(idRes) == 34 || abs(idRes) == 41)
2265 0 : && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 1;
2266 :
2267 : // g + g, gamma + gamma -> Higgs boson.
2268 0 : if ( (idRes == 25 || idRes == 35 || idRes == 36)
2269 0 : && ( ( idIn1 == 21 && idIn2 == 21 )
2270 0 : || ( idIn1 == 22 && idIn2 == 22 ) ) ) MEtype = 2;
2271 :
2272 : // f + fbar -> Higgs boson.
2273 0 : if ( (idRes == 25 || idRes == 35 || idRes == 36)
2274 0 : && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 3;
2275 0 : }
2276 :
2277 : // Weak ME corrections.
2278 0 : if (weakRadiation) {
2279 0 : if (event[3].id() == -event[4].id()
2280 0 : || event[event[3].daughter1()].idAbs() == 24 || infoPtr->nFinal() != 2)
2281 0 : MEtype = 200;
2282 0 : else if (event[3].idAbs() == 21 || event[4].idAbs() == 21) MEtype = 201;
2283 0 : else if (event[3].id() == event[4].id()) MEtype = 202;
2284 : else MEtype = 203;
2285 : }
2286 :
2287 : // Done.
2288 0 : return MEtype;
2289 :
2290 0 : }
2291 :
2292 : //--------------------------------------------------------------------------
2293 :
2294 : // Provide maximum of expected ME weight; for preweighting of evolution.
2295 :
2296 : double SpaceShower::calcMEmax( int MEtype, int idMother, int idDaughterIn) {
2297 :
2298 : // Main non-unity case: g(gamma) f -> V f'.
2299 0 : if (MEtype == 1 && idMother > 20 && idDaughterIn < 20) return 3.;
2300 :
2301 : // Added a case for t-channel W/Z exchange, since the PS is not an
2302 : // overestimate. This does not help fully, but it should only be small
2303 : // pT quarks / gluons that break the overscattering.
2304 0 : if ( MEtype == 201 || MEtype == 202 || MEtype == 203
2305 0 : || MEtype == 206 || MEtype == 207 || MEtype == 208) return WEAKPSWEIGHT;
2306 :
2307 : // Default.
2308 0 : return 1.;
2309 :
2310 0 : }
2311 :
2312 : //--------------------------------------------------------------------------
2313 :
2314 : // Provide actual ME weight for current branching.
2315 : // Note: currently ME corrections are only allowed for first branching
2316 : // on each side, so idDaughter is essentially known and checks overkill.
2317 :
2318 : double SpaceShower::calcMEcorr(int MEtype, int idMother, int idDaughterIn,
2319 : double M2, double z, double Q2, double m2Sister) {
2320 :
2321 : // Convert to Mandelstam variables. Sometimes may need to swap later.
2322 0 : double sH = M2 / z;
2323 0 : double tH = -Q2;
2324 0 : double uH = Q2 - M2 * (1. - z) / z;
2325 0 : int idMabs = abs(idMother);
2326 0 : int idDabs = abs(idDaughterIn);
2327 :
2328 : // Corrections for f + fbar -> s-channel vector boson.
2329 0 : if (MEtype == 1) {
2330 0 : if (idMabs < 20 && idDabs < 20) {
2331 0 : return (tH*tH + uH*uH + 2. * M2 * sH) / (sH*sH + M2*M2);
2332 0 : } else if (idDabs < 20) {
2333 : // g(gamma) f -> V f': -Q2 = (p_g - p_f')^2 in PS while
2334 : // tHat = (p_f - p_f')^2 in ME so need to swap tHat <-> uHat.
2335 0 : swap( tH, uH);
2336 0 : return (sH*sH + uH*uH + 2. * M2 * tH) / (pow2(sH - M2) + M2*M2);
2337 : }
2338 :
2339 : // Corrections for g + g -> Higgs boson.
2340 0 : } else if (MEtype == 2) {
2341 0 : if (idMabs < 20 && idDabs > 20) {
2342 0 : return (sH*sH + uH*uH) / (sH*sH + pow2(sH - M2));
2343 0 : } else if (idDabs > 20) {
2344 0 : return 0.5 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(M2))
2345 0 : / pow2(sH*sH - M2 * (sH - M2));
2346 : }
2347 :
2348 : // Corrections for f + fbar -> Higgs boson (f = b mainly).
2349 0 : } else if (MEtype == 3) {
2350 0 : if (idMabs < 20 && idDabs < 20) {
2351 : // The PS and ME answers agree for f fbar -> H g/gamma.
2352 0 : return 1.;
2353 0 : } else if (idDabs < 20) {
2354 : // Need to swap tHat <-> uHat, cf. vector-boson production above.
2355 0 : swap( tH, uH);
2356 0 : return (sH*sH + uH*uH + 2. * (M2 - uH) * (M2 - sH))
2357 0 : / (pow2(sH - M2) + M2*M2);
2358 : }
2359 :
2360 : // Corrections for f -> f' + W/Z (s-channel).
2361 0 : } else if (MEtype == 200 || MEtype == 205) {
2362 : // Need to redo calculations of uH since we now emit a massive particle.
2363 0 : uH += m2Sister;
2364 0 : double wtME = (uH*uH + tH*tH + 2 * sH * (m2Sister + M2)) / (uH*tH)
2365 0 : - M2 * m2Sister * (1/(tH*tH) + 1/(uH*uH));
2366 0 : double wtPS = (sH*sH + pow2(M2 + m2Sister)) / (tH*uH);
2367 0 : return wtME / wtPS;
2368 0 : } else if (MEtype == 201 || MEtype == 202 || MEtype == 203 ||
2369 0 : MEtype == 206 || MEtype == 207 || MEtype == 208)
2370 0 : return calcMEmax(MEtype, 0, 0);
2371 :
2372 : // Default.
2373 0 : return 1.;
2374 :
2375 0 : }
2376 :
2377 : //--------------------------------------------------------------------------
2378 :
2379 : // Provide actual ME weight for current branching for weak t-channel emissions.
2380 :
2381 : double SpaceShower::calcMEcorrWeak(int MEtype, double m2, double z,
2382 : double pT2, Vec4 pMother, Vec4 pB, Vec4 pDaughter,
2383 : Vec4 pB0, Vec4 p1, Vec4 p2, Vec4 pSister) {
2384 :
2385 : // Find daughter four-momentum in current frame.
2386 0 : Vec4 pA = pMother - pSister;
2387 :
2388 : // Scale outgoing vectors to conserve energy / momentum.
2389 0 : double scaleFactor2 = (pA + pB).m2Calc() / (p1 + p2).m2Calc();
2390 0 : double scaleFactor = sqrt(scaleFactor2);
2391 0 : RotBstMatrix rot2to2frame;
2392 0 : rot2to2frame.bstback(p1 + p2);
2393 0 : p1.rotbst(rot2to2frame);
2394 0 : p2.rotbst(rot2to2frame);
2395 0 : p1 *= scaleFactor;
2396 0 : p2 *= scaleFactor;
2397 :
2398 : // Find 2 to 2 rest frame for incoming particles.
2399 : // This is done before one of the two are made virtual (Q^2 mass).
2400 0 : RotBstMatrix rot2to2frameInc;
2401 0 : rot2to2frameInc.bstback(pDaughter + pB0);
2402 0 : pDaughter.rotbst(rot2to2frameInc);
2403 0 : pB0.rotbst(rot2to2frameInc);
2404 0 : double sHat = (p1 + p2).m2Calc();
2405 0 : double tHat = (p1 - pDaughter).m2Calc();
2406 0 : double uHat = (p1 - pB0).m2Calc();
2407 :
2408 : // Calculate the weak t-channel correction.
2409 0 : double m2R1 = 1. + pSister.m2Calc() / m2;
2410 0 : double wt = 4. * sHat / (pMother + pB).m2Calc() * pT2 * ( 1. - z * m2R1)
2411 0 : / (1. + pow2(z * m2R1)) / (1.-z);
2412 0 : if (MEtype == 201 || MEtype == 206)
2413 0 : wt *= weakShowerMEs.getTchanneluGuGZME(pMother, pB, p2, pSister, p1)
2414 0 : / weakShowerMEs.getTchanneluGuGME(sHat, tHat, uHat);
2415 0 : else if (MEtype == 202 || MEtype == 207)
2416 0 : wt *= weakShowerMEs.getTchannelududZME(pMother, pB, pSister, p2, p1)
2417 0 : / weakShowerMEs.getTchanneluuuuME(sHat, tHat, uHat);
2418 0 : else if (MEtype == 203 || MEtype == 208)
2419 0 : wt *= weakShowerMEs.getTchannelududZME(pMother, pB, pSister, p2, p1)
2420 0 : / weakShowerMEs.getTchannelududME(sHat, tHat, uHat);
2421 :
2422 : // Split of ME into an ISR part and FSR part.
2423 0 : wt *= (pSister + p1).m2Calc() / ( (pSister + p1).m2Calc()
2424 0 : + abs((-pMother + pSister).m2Calc()) );
2425 :
2426 : // Remove the addition weight that was used to get an overestimate.
2427 0 : wt /= calcMEmax(MEtype, 0, 0);
2428 :
2429 0 : return wt;
2430 0 : }
2431 :
2432 : //--------------------------------------------------------------------------
2433 :
2434 : // Find coefficient of azimuthal asymmetry from gluon polarization.
2435 :
2436 : void SpaceShower::findAsymPol( Event& event, SpaceDipoleEnd* dip) {
2437 :
2438 : // Default is no asymmetry. Only gluons are studied.
2439 0 : dip->iFinPol = 0;
2440 0 : dip->asymPol = 0.;
2441 0 : int iRad = dip->iRadiator;
2442 0 : if (!doPhiPolAsym || dip->idDaughter != 21) return;
2443 :
2444 : // At least two particles in final state, whereof at least one coloured.
2445 0 : int systemSizeOut = partonSystemsPtr->sizeOut( iSysSel);
2446 0 : if (systemSizeOut < 2) return;
2447 : bool foundColOut = false;
2448 0 : for (int ii = 0; ii < systemSizeOut; ++ii) {
2449 0 : int i = partonSystemsPtr->getOut( iSysSel, ii);
2450 0 : if (event[i].col() != 0 || event[i].acol() != 0) foundColOut = true;
2451 : }
2452 0 : if (!foundColOut) return;
2453 :
2454 : // Check if granddaughter in final state of hard scattering.
2455 : // (May need to trace across carbon copies to find granddaughters.)
2456 : // If so, only accept 2 -> 2 scatterings with gg or qq in final state.
2457 0 : int iGrandD1 = event[iRad].daughter1();
2458 0 : int iGrandD2 = event[iRad].daughter2();
2459 : bool traceCopy = false;
2460 0 : do {
2461 : traceCopy = false;
2462 0 : if (iGrandD1 > 0 && iGrandD2 == iGrandD1) {
2463 0 : iGrandD1 = event[iGrandD2].daughter1();
2464 0 : iGrandD2 = event[iGrandD2].daughter2();
2465 : traceCopy = true;
2466 0 : }
2467 0 : } while (traceCopy);
2468 0 : int statusGrandD1 = event[ iGrandD1 ].statusAbs();
2469 0 : bool isHardProc = (statusGrandD1 == 23 || statusGrandD1 == 33);
2470 0 : if (isHardProc) {
2471 0 : if (iGrandD2 != iGrandD1 + 1) return;
2472 0 : if (event[iGrandD1].isGluon() && event[iGrandD2].isGluon());
2473 0 : else if (event[iGrandD1].isQuark() && event[iGrandD2].isQuark());
2474 0 : else return;
2475 : }
2476 0 : dip->iFinPol = iGrandD1;
2477 :
2478 : // Coefficient from gluon production.
2479 0 : if (dip->idMother == 21) dip->asymPol = pow2( (1. - dip->z)
2480 0 : / (1. - dip->z * (1. - dip->z) ) );
2481 0 : else dip->asymPol = 2. * (1. - dip->z) / (1. + pow2(1. - dip->z) );
2482 :
2483 : // Coefficients from gluon decay. Put z = 1/2 for hard process.
2484 0 : double zDau = (isHardProc) ? 0.5 : dip->zOld;
2485 0 : if (event[iGrandD1].isGluon()) dip->asymPol *= pow2( (1. - zDau)
2486 0 : / (1. - zDau * (1. - zDau) ) );
2487 0 : else dip->asymPol *= -2. * zDau *( 1. - zDau )
2488 0 : / (1. - 2. * zDau * (1. - zDau) );
2489 :
2490 0 : }
2491 :
2492 : //--------------------------------------------------------------------------
2493 :
2494 : // Remove weak dipoles if FSR already emitted a W/Z
2495 : // and only a single weak emission is permited.
2496 :
2497 : void SpaceShower::update(int , Event &, bool hasWeakRad) {
2498 :
2499 0 : if (hasWeakRad && singleWeakEmission)
2500 0 : for (int i = 0; i < int(dipEnd.size()); i++)
2501 0 : if (dipEnd[i].weakType != 0) dipEnd[i].weakType = 0;
2502 0 : if (hasWeakRad) hasWeaklyRadiated = true;
2503 :
2504 0 : }
2505 :
2506 : //-------------------------------------------------------------------------
2507 :
2508 : // Print the list of dipoles.
2509 :
2510 : void SpaceShower::list(ostream& os) const {
2511 :
2512 : // Header.
2513 0 : os << "\n -------- PYTHIA SpaceShower Dipole Listing -------------- \n"
2514 0 : << "\n i syst side rad rec pTmax col chg ME rec \n"
2515 0 : << fixed << setprecision(3);
2516 :
2517 : // Loop over dipole list and print it.
2518 0 : for (int i = 0; i < int(dipEnd.size()); ++i)
2519 0 : os << setw(5) << i << setw(6) << dipEnd[i].system
2520 0 : << setw(6) << dipEnd[i].side << setw(6) << dipEnd[i].iRadiator
2521 0 : << setw(6) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax
2522 0 : << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType
2523 0 : << setw(5) << dipEnd[i].MEtype << setw(4)
2524 0 : << dipEnd[i].normalRecoil << "\n";
2525 :
2526 : // Done.
2527 0 : os << "\n -------- End PYTHIA SpaceShower Dipole Listing ----------"
2528 0 : << endl;
2529 :
2530 :
2531 0 : }
2532 :
2533 : //==========================================================================
2534 :
2535 : } // end namespace Pythia8
|
