StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SpaceShower.cc
1 // SpaceShower.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 const double SpaceShower::HEADROOMHQG = 1.35;
88 
89 // Limit on size of number of rejections for uncertainty variations.
90 const double SpaceShower::REJECTFACTOR = 0.1;
91 
92 // Limit on probability for uncertainty variations.
93 const double SpaceShower::PROBLIMIT = 0.99;
94 
95 //--------------------------------------------------------------------------
96 
97 // Initialize alphaStrong, alphaEM and related pTmin parameters.
98 
99 void SpaceShower::init( BeamParticle* beamAPtrIn,
100  BeamParticle* beamBPtrIn) {
101 
102  // Store input pointers for future use.
103  beamAPtr = beamAPtrIn;
104  beamBPtr = beamBPtrIn;
105 
106  // Main flags to switch on and off branchings.
107  doQCDshower = settingsPtr->flag("SpaceShower:QCDshower");
108  doQEDshowerByQ = settingsPtr->flag("SpaceShower:QEDshowerByQ");
109  doQEDshowerByL = settingsPtr->flag("SpaceShower:QEDshowerByL");
110  doWeakShower = settingsPtr->flag("SpaceShower:WeakShower");
111 
112  // Matching in pT of hard interaction to shower evolution.
113  pTmaxMatch = settingsPtr->mode("SpaceShower:pTmaxMatch");
114  pTdampMatch = settingsPtr->mode("SpaceShower:pTdampMatch");
115  pTmaxFudge = settingsPtr->parm("SpaceShower:pTmaxFudge");
116  pTmaxFudgeMPI = settingsPtr->parm("SpaceShower:pTmaxFudgeMPI");
117  pTdampFudge = settingsPtr->parm("SpaceShower:pTdampFudge");
118 
119  // Optionally force emissions to be ordered in rapidity/angle.
120  doRapidityOrder = settingsPtr->flag("SpaceShower:rapidityOrder");
121  doRapidityOrderMPI = settingsPtr->flag("SpaceShower:rapidityOrderMPI");
122 
123  // Charm, bottom and lepton mass thresholds.
124  mc = max( MCMIN, particleDataPtr->m0(4));
125  mb = max( MBMIN, particleDataPtr->m0(5));
126  m2c = pow2(mc);
127  m2b = pow2(mb);
128 
129  // Parameters of scale choices.
130  renormMultFac = settingsPtr->parm("SpaceShower:renormMultFac");
131  factorMultFac = settingsPtr->parm("SpaceShower:factorMultFac");
132  useFixedFacScale = settingsPtr->flag("SpaceShower:useFixedFacScale");
133  fixedFacScale2 = pow2(settingsPtr->parm("SpaceShower:fixedFacScale"));
134 
135  // Parameters of alphaStrong generation.
136  alphaSvalue = settingsPtr->parm("SpaceShower:alphaSvalue");
137  alphaSorder = settingsPtr->mode("SpaceShower:alphaSorder");
138  alphaSnfmax = settingsPtr->mode("StandardModel:alphaSnfmax");
139  alphaSuseCMW = settingsPtr->flag("SpaceShower:alphaSuseCMW");
140  alphaS2pi = 0.5 * alphaSvalue / M_PI;
141 
142  // Initialize alpha_strong generation.
143  alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW);
144 
145  // Lambda for 5, 4 and 3 flavours.
146  Lambda5flav = alphaS.Lambda5();
147  Lambda4flav = alphaS.Lambda4();
148  Lambda3flav = alphaS.Lambda3();
149  Lambda5flav2 = pow2(Lambda5flav);
150  Lambda4flav2 = pow2(Lambda4flav);
151  Lambda3flav2 = pow2(Lambda3flav);
152 
153  // Regularization of QCD evolution for pT -> 0. Can be taken
154  // same as for multiparton interactions, or be set separately.
155  useSamePTasMPI = settingsPtr->flag("SpaceShower:samePTasMPI");
156  if (useSamePTasMPI) {
157 
158  // Different parametrization for photon-photon collisions.
159  if (beamAPtr->isGamma() && beamBPtr->isGamma()) {
160  pT0paramMode = settingsPtr->mode("PhotonPhoton:pT0parametrization");
161  pT0Ref = settingsPtr->parm("PhotonPhoton:pT0Ref");
162  ecmRef = settingsPtr->parm("PhotonPhoton:ecmRef");
163  ecmPow = settingsPtr->parm("PhotonPhoton:ecmPow");
164  pTmin = settingsPtr->parm("PhotonPhoton:pTmin");
165  } else {
166  pT0paramMode
167  = settingsPtr->mode("MultipartonInteractions:pT0parametrization");
168  pT0Ref = settingsPtr->parm("MultipartonInteractions:pT0Ref");
169  ecmRef = settingsPtr->parm("MultipartonInteractions:ecmRef");
170  ecmPow = settingsPtr->parm("MultipartonInteractions:ecmPow");
171  pTmin = settingsPtr->parm("MultipartonInteractions:pTmin");
172  }
173 
174  } else {
175  pT0paramMode = settingsPtr->mode("SpaceShower:pT0parametrization");
176  pT0Ref = settingsPtr->parm("SpaceShower:pT0Ref");
177  ecmRef = settingsPtr->parm("SpaceShower:ecmRef");
178  ecmPow = settingsPtr->parm("SpaceShower:ecmPow");
179  pTmin = settingsPtr->parm("SpaceShower:pTmin");
180  }
181 
182  // Calculate nominal invariant mass of events.
183  sCM = m2( beamAPtr->p(), beamBPtr->p());
184  eCM = sqrt(sCM);
185 
186  // Set current pT0 scale according to the chosen parametrization.
187  if (pT0paramMode == 0) pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
188  else pT0 = pT0Ref + ecmPow * log (eCM / ecmRef);
189 
190  // Restrict pTmin to ensure that alpha_s(pTmin^2 + pT_0^2) does not blow up.
191  double pTminAbs = sqrtpos(pow2(LAMBDA3MARGIN) * Lambda3flav2 / renormMultFac
192  - pT0*pT0);
193  if (pTmin < pTminAbs) {
194  pTmin = pTminAbs;
195  ostringstream newPTmin;
196  newPTmin << fixed << setprecision(3) << pTmin;
197  infoPtr->errorMsg("Warning in SpaceShower::init: pTmin too low",
198  ", raised to " + newPTmin.str() );
199  infoPtr->setTooLowPTmin(true);
200  }
201 
202  // Parameters of alphaEM generation.
203  alphaEMorder = settingsPtr->mode("SpaceShower:alphaEMorder");
204 
205  // Initialize alphaEM generation.
206  alphaEM.init( alphaEMorder, settingsPtr);
207 
208  // Parameters of QED evolution.
209  pTminChgQ = settingsPtr->parm("SpaceShower:pTminchgQ");
210  pTminChgL = settingsPtr->parm("SpaceShower:pTminchgL");
211 
212  // Derived parameters of QCD evolution.
213  pT20 = pow2(pT0);
214  pT2min = pow2(pTmin);
215  pT2minChgQ = pow2(pTminChgQ);
216  pT2minChgL = pow2(pTminChgL);
217 
218  // Parameters of weak evolution.
219  weakMode = settingsPtr->mode("SpaceShower:weakShowerMode");
220  pTweakCut = settingsPtr->parm("SpaceShower:pTminWeak");
221  pT2weakCut = pow2(pTweakCut);
222  weakEnhancement = settingsPtr->parm("WeakShower:enhancement");
223  singleWeakEmission = settingsPtr->flag("WeakShower:singleEmission");
224  vetoWeakJets = settingsPtr->flag("WeakShower:vetoWeakJets");
225  vetoWeakDeltaR2 = pow2(settingsPtr->parm("weakShower:vetoWeakDeltaR"));
226  weakExternal = settingsPtr->flag("WeakShower:externalSetup");
227 
228  // Various other parameters.
229  doMEcorrections = settingsPtr->flag("SpaceShower:MEcorrections");
230  doMEafterFirst = settingsPtr->flag("SpaceShower:MEafterFirst");
231  doPhiPolAsym = settingsPtr->flag("SpaceShower:phiPolAsym");
232  doPhiPolAsymHard = settingsPtr->flag("SpaceShower:phiPolAsymHard");
233  doPhiIntAsym = settingsPtr->flag("SpaceShower:phiIntAsym");
234  strengthIntAsym = settingsPtr->parm("SpaceShower:strengthIntAsym");
235  nQuarkIn = settingsPtr->mode("SpaceShower:nQuarkIn");
236 
237  // Do not do phiIntAsym if dipoleRecoil is on, to avoid doublecounting.
238  doDipoleRecoil = settingsPtr->flag("SpaceShower:dipoleRecoil");
239  if (doDipoleRecoil) doPhiIntAsym = false;
240 
241  // Z0 and W+- properties needed for weak showers.
242  mZ = particleDataPtr->m0(23);
243  gammaZ = particleDataPtr->mWidth(23);
244  thetaWRat = 1. / (16. * coupSMPtr->sin2thetaW()
245  * coupSMPtr->cos2thetaW());
246  mW = particleDataPtr->m0(24);
247  gammaW = particleDataPtr->mWidth(24);
248 
249  // Possibility of two predetermined hard emissions in event.
250  doSecondHard = settingsPtr->flag("SecondHard:generate");
251 
252  // gamma->qqbar splittings handled differently with and without MPIs.
253  doMPI = settingsPtr->flag("PartonLevel:MPI");
254  gamma2qqbar = false;
255 
256  // Optional dampening at small pT's when large multiplicities.
257  enhanceScreening
258  = settingsPtr->mode("MultipartonInteractions:enhanceScreening");
259  if (!useSamePTasMPI) enhanceScreening = 0;
260 
261  // Possibility to allow user veto of emission step.
262  hasUserHooks = (userHooksPtr != 0);
263  canVetoEmission = hasUserHooks && userHooksPtr->canVetoISREmission();
264 
265  // Default values for the weak shower.
266  hasWeaklyRadiated = false;
267  weakMaxWt = 1.;
268 
269  // Disallow simultaneous splitting and trial emission enhancements.
270  canEnhanceEmission = hasUserHooks && userHooksPtr->canEnhanceEmission();
271  canEnhanceTrial = hasUserHooks && userHooksPtr->canEnhanceTrial();
272  if (canEnhanceEmission && canEnhanceTrial) {
273  infoPtr->errorMsg("Error in SpaceShower::init: Enhance for both actual "
274  "and trial emissions not possible. Both switched off.");
275  canEnhanceEmission = false;
276  canEnhanceTrial = false;
277  }
278 
279  // Properties for enhanced emissions.
280  splittingNameSel = "";
281  splittingNameNow = "";
282  enhanceFactors.clear();
283 
284  // Enable automated uncertainty variations.
285  nVarQCD = 0;
286  doUncertainties = settingsPtr->flag("UncertaintyBands:doVariations")
287  && initUncertainties();
288  doUncertaintiesNow = doUncertainties;
289  uVarNflavQ = settingsPtr->mode("UncertaintyBands:nFlavQ");
290  uVarMPIshowers = settingsPtr->flag("UncertaintyBands:MPIshowers");
291  cNSpTmin = settingsPtr->parm("UncertaintyBands:cNSpTmin");
292  uVarpTmin2 = pow2(pT0Ref);
293  uVarpTmin2 *= settingsPtr->parm("UncertaintyBands:FSRpTmin2Fac");
294  overFactor = settingsPtr->parm("UncertaintyBands:overSampleISR");
295 
296  // Possibility to set parton vertex information.
297  doPartonVertex = settingsPtr->flag("PartonVertex:setVertex")
298  && (partonVertexPtr != 0);
299 
300 }
301 
302 //--------------------------------------------------------------------------
303 
304 // Find whether to limit maximum scale of emissions.
305 // Also allow for dampening at factorization or renormalization scale.
306 
307 bool SpaceShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) {
308 
309  // Find whether to limit pT. Begin by user-set cases.
310  bool dopTlimit = false;
311  dopTlimit1 = dopTlimit2 = false;
312  int nHeavyCol = 0;
313  if (pTmaxMatch == 1) dopTlimit = dopTlimit1 = dopTlimit2 = true;
314  else if (pTmaxMatch == 2) dopTlimit = dopTlimit1 = dopTlimit2 = false;
315 
316  // Always restrict SoftQCD processes.
317  else if (infoPtr->isNonDiffractive() || infoPtr->isDiffractiveA()
318  || infoPtr->isDiffractiveB() || infoPtr->isDiffractiveC() )
319  dopTlimit = dopTlimit1 = dopTlimit2 = true;
320 
321  // Look if any quark (u, d, s, c, b), gluon or photon in final state.
322  // Also count number of heavy coloured particles, like top.
323  else {
324  int n21 = 0;
325  int iBegin = 5 + beamOffset;
326  for (int i = iBegin; i < event.size(); ++i) {
327  if (event[i].status() == -21) ++n21;
328  else if (n21 == 0) {
329  int idAbs = event[i].idAbs();
330  if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit1 = true;
331  if ( (event[i].col() != 0 || event[i].acol() != 0)
332  && idAbs > 5 && idAbs != 21 ) ++nHeavyCol;
333  } else if (n21 == 2) {
334  int idAbs = event[i].idAbs();
335  if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit2 = true;
336  }
337  }
338  dopTlimit = (doSecondHard) ? (dopTlimit1 && dopTlimit2) : dopTlimit1;
339  }
340 
341  // Dampening at factorization or renormalization scale; only for hardest.
342  dopTdamp = false;
343  pT2damp = 0.;
344  if (!dopTlimit1 && (pTdampMatch == 1 || pTdampMatch == 2)) {
345  dopTdamp = true;
346  pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren);
347  }
348  if (!dopTlimit1 && nHeavyCol > 1 && (pTdampMatch == 3 || pTdampMatch == 4)) {
349  dopTdamp = true;
350  pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 3) ? Q2Fac : Q2Ren);
351  }
352 
353  // Done.
354  return dopTlimit;
355 
356 }
357 
358 //--------------------------------------------------------------------------
359 
360 // Prepare system for evolution; identify ME.
361 // Routine may be called after multiparton interactions, for a new subystem.
362 
363 void SpaceShower::prepare( int iSys, Event& event, bool limitPTmaxIn) {
364 
365  // Reset W/Z radiation flag and counters at first call for new event.
366  if (iSys == 0) {
367  nRadA.clear();
368  nRadB.clear();
369  hasWeaklyRadiated = false;
370  }
371 
372  // Find positions of incoming colliding partons.
373  int in1 = partonSystemsPtr->getInA(iSys);
374  int in2 = partonSystemsPtr->getInB(iSys);
375 
376  // Rescattered partons cannot radiate.
377  bool canRadiate1 = !(event[in1].isRescatteredIncoming());
378  bool canRadiate2 = !(event[in2].isRescatteredIncoming());
379 
380  // Reset dipole-ends list for first interaction. Also resonances.
381  if (iSys == 0) dipEnd.resize(0);
382  if (iSys == 0) idResFirst = 0;
383  if (iSys == 1) idResSecond = 0;
384 
385  // Find matrix element corrections for system.
386  int MEtype = findMEtype( iSys, event, false);
387 
388  // In case of DPS overwrite limitPTmaxIn by saved value.
389  if (doSecondHard && iSys == 0) limitPTmaxIn = dopTlimit1;
390  if (doSecondHard && iSys == 1) limitPTmaxIn = dopTlimit2;
391 
392  // Maximum pT scale for dipole ends.
393  double pTmax1 = (limitPTmaxIn) ? event[in1].scale() : eCM;
394  double pTmax2 = (limitPTmaxIn) ? event[in2].scale() : eCM;
395  if ( limitPTmaxIn && (iSys == 0 || (iSys == 1 && doSecondHard)) ) {
396  pTmax1 *= pTmaxFudge;
397  pTmax2 *= pTmaxFudge;
398  } else if (limitPTmaxIn && iSys > 0) {
399  pTmax1 *= pTmaxFudgeMPI;
400  pTmax2 *= pTmaxFudgeMPI;
401  }
402 
403  // Find dipole ends for QCD radiation.
404  // Note: colour type can change during evolution, so book also if zero.
405  if (doQCDshower) {
406  int colType1 = event[in1].colType();
407  if (canRadiate1) {
408  // Look if there is an IF dipole in case of dipole recoil.
409  int iColPartner = (doDipoleRecoil)
410  ? findColPartner(event, in1, in2, iSys) : 0;
411  int idColPartner = (iColPartner != 0) ? event[iColPartner].id() : 0;
412  dipEnd.push_back( SpaceDipoleEnd( iSys, 1, in1, in2, pTmax1,
413  colType1, 0, 0, MEtype, canRadiate2, 0, iColPartner, idColPartner) );
414  }
415  int colType2 = event[in2].colType();
416  if (canRadiate2) {
417  // Look if there is an IF dipole in case of dipole recoil.
418  int iColPartner = (doDipoleRecoil)
419  ? findColPartner(event, in2, in1, iSys) : 0;
420  int idColPartner = (iColPartner != 0) ? event[iColPartner].id() : 0;
421  dipEnd.push_back( SpaceDipoleEnd( iSys, 2, in2, in1, pTmax2,
422  colType2, 0, 0, MEtype, canRadiate1, 0, iColPartner, idColPartner) );
423  }
424  }
425 
426  // Find dipole ends for QED radiation.
427  // Note: charge type can change during evolution, so book also if zero.
428  if (doQEDshowerByQ || doQEDshowerByL) {
429  int chgType1 = ( (event[in1].isQuark() && doQEDshowerByQ)
430  || (event[in1].isLepton() && doQEDshowerByL) )
431  ? event[in1].chargeType() : 0;
432  // Special: photons have charge zero, but can evolve (only off Q for now)
433  if (event[in1].id() == 22 && doQEDshowerByQ) chgType1 = 22 ;
434  if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys, -1,
435  in1, in2, pTmax1, 0, chgType1, 0, MEtype, canRadiate2) );
436  int chgType2 = ( (event[in2].isQuark() && doQEDshowerByQ)
437  || (event[in2].isLepton() && doQEDshowerByL) )
438  ? event[in2].chargeType() : 0;
439  // Special: photons have charge zero, but can evolve (only off Q for now)
440  if (event[in2].id() == 22 && doQEDshowerByQ) chgType2 = 22 ;
441  if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys, -2,
442  in2, in1, pTmax2, 0, chgType2, 0, MEtype, canRadiate1) );
443  }
444 
445  // Find dipole ends for weak radiation. No right-handed W emission.
446  // Currently leptons are not allow to emit W bosons and only
447  // emissions from the hard process are included.
448  if (doWeakShower && iSys == 0) {
449  // Normal internal setup.
450  if (!weakExternal) {
451  // Determine what type of 2 -> 2 process it is.
452  int MEtypeWeak = findMEtype( iSys, event, true);
453  if (MEtypeWeak == 201 || MEtypeWeak == 202 || MEtypeWeak == 203 ||
454  MEtypeWeak == 206 || MEtypeWeak == 207 || MEtypeWeak == 208) {
455 
456  // Nonidentical incoming flavours.
457  if (event[in1].id() != event[in2].id()) {
458  if (event[in1].id() == event[in1 + 2].id()) tChannel = true;
459  else if (event[in2].id() == event[in1 + 2].id()) tChannel = false;
460  // No quark matches the outgoing, choose randomly.
461  else tChannel = (rndmPtr->flat() < 0.5);
462  // In case of same quark flavours, choose randomly.
463  } else tChannel = (rndmPtr->flat() < 0.5);
464  }
465 
466  // Set up weak dipole ends for first incoming parton.
467  int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
468  if (event[in1].idAbs() < 20) event[in1].pol(weakPol);
469  if (canRadiate1) {
470  if ( (weakMode == 0 || weakMode == 1) && weakPol == -1
471  && event[in1].isQuark() )
472  dipEnd.push_back( SpaceDipoleEnd( iSys, 1, in1, in2, pTmax1,
473  0, 0, 1, MEtypeWeak, canRadiate2, weakPol) );
474  if ( (weakMode == 0 || weakMode == 2)
475  && (event[in1].isQuark() || event[in1].isLepton()) )
476  dipEnd.push_back( SpaceDipoleEnd( iSys, 1, in1, in2, pTmax1,
477  0, 0, 2, MEtypeWeak + 5, canRadiate2, weakPol) );
478  }
479 
480  // Set up weak dipole ends for second incoming parton.
481  if (event[in1].id() != - event[in2].id())
482  weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
483  if (event[in2].idAbs() < 20) event[in2].pol(weakPol);
484  if (canRadiate2) {
485  if ( (weakMode == 0 || weakMode == 1) && weakPol == -1
486  && event[in2].isQuark())
487  dipEnd.push_back( SpaceDipoleEnd( iSys, 2, in2, in1, pTmax2,
488  0, 0, 1, MEtypeWeak, canRadiate1, weakPol) );
489  if ( (weakMode == 0 || weakMode == 2) &&
490  (event[in2].isQuark() || event[in2].isLepton()) )
491  dipEnd.push_back( SpaceDipoleEnd( iSys, 2, in2, in1, pTmax2,
492  0, 0, 2, MEtypeWeak + 5, canRadiate1, weakPol) );
493  }
494  // External setup from infoPtr.
495  } else {
496  // Get information.
497  vector<pair<int,int> > weakDipoles = infoPtr->getWeakDipoles();
498  vector<int> weakModes = infoPtr->getWeakModes();
499  weakMomenta = infoPtr->getWeakMomenta();
500  tChannel = true;
501  // Loop over dipoles.
502  for (int i = 0; i < int(weakDipoles.size()); ++i) {
503  // Only consider ISR dipoles.
504 
505  if (event[weakDipoles[i].first].status() < 0) {
506  // Find ME.
507  int iRadLocal = weakDipoles[i].first;
508  int iRecLocal = weakDipoles[i].second;
509  int side = (iRadLocal == 3) ? 1 : 2;
510  double pTmax = (side == 1) ? pTmax1 : pTmax2;
511 
512  // Find MEtype.
513  int MEtypeWeak = 0;
514  if (weakModes[weakDipoles[i].first] == 1) MEtypeWeak = 200;
515  else if (weakModes[weakDipoles[i].first] == 2) MEtypeWeak = 201;
516  else if (weakModes[weakDipoles[i].first] == 3) MEtypeWeak = 202;
517  else MEtypeWeak = 203;
518 
519  // Find correct polarization, if it is already set use it.
520  // Otherwise pick randomly.
521  int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
522  if (event[weakDipoles[i].first].intPol() != 9)
523  weakPol = event[weakDipoles[i].first].intPol();
524  event[weakDipoles[i].first].pol(weakPol);
525 
526  // Add the dipoles.
527  if ( (weakMode == 0 || weakMode == 1) && weakPol == -1)
528  dipEnd.push_back( SpaceDipoleEnd( iSys, side, iRadLocal, iRecLocal,
529  pTmax, 0, 0, 1, MEtypeWeak, true, weakPol) );
530  if (weakMode == 0 || weakMode == 2)
531  dipEnd.push_back( SpaceDipoleEnd( iSys, side, iRadLocal, iRecLocal,
532  pTmax, 0, 0, 2, MEtypeWeak + 5, true, weakPol) );
533  }
534  }
535  }
536  }
537 
538  // Store the z and pT2 values of the last previous splitting
539  // when an event history has already been constructed.
540  if (iSys == 0 && infoPtr->hasHistory()) {
541  double zNow = infoPtr->zNowISR();
542  double pT2Now = infoPtr->pT2NowISR();
543  for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
544  dipEnd[iDipEnd].zOld = zNow;
545  dipEnd[iDipEnd].pT2Old = pT2Now;
546  ++dipEnd[iDipEnd].nBranch;
547  }
548  }
549 
550 }
551 
552 //--------------------------------------------------------------------------
553 
554 // Select next pT in downwards evolution of the existing dipoles.
555 
556 double SpaceShower::pTnext( Event& event, double pTbegAll, double pTendAll,
557  int nRadIn, bool doTrialIn) {
558 
559  // Current cm energy, in case it varies between events.
560  sCM = m2( beamAPtr->p(), beamBPtr->p());
561  eCM = sqrt(sCM);
562  pTbegRef = pTbegAll;
563 
564  // Starting values: no radiating dipole found.
565  nRad = nRadIn;
566  double pT2sel = pow2(pTendAll);
567  iDipSel = 0;
568  iSysSel = 0;
569  dipEndSel = 0;
570 
571  // Check if enhanced emissions should be applied.
572  doTrialNow = doTrialIn;
573  canEnhanceET = (!doTrialNow && canEnhanceEmission)
574  || ( doTrialNow && canEnhanceTrial);
575 
576  // Starting values for enhanced emissions.
577  splittingNameSel = "";
578  splittingNameNow = "";
579  enhanceFactors.clear();
580  if (hasUserHooks) userHooksPtr->setEnhancedTrial(0., 1.);
581 
582  // Loop over all possible dipole ends.
583  for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
584  iDipNow = iDipEnd;
585  dipEndNow = &dipEnd[iDipEnd];
586  iSysNow = dipEndNow->system;
587  dipEndNow->pT2 = 0.;
588  dipEndNow->pAccept = 1.0;
589  double pTbegDip = min( pTbegAll, dipEndNow->pTmax );
590 
591  // Check whether dipole end should be allowed to shower.
592  double pT2begDip = pow2(pTbegDip);
593  if (pT2begDip > pT2sel && ( dipEndNow->colType != 0
594  || dipEndNow->chgType != 0 || dipEndNow->weakType != 0) ) {
595  double pT2endDip = 0.;
596 
597  // Determine lower cut for evolution, for QCD or QED (q or l).
598  if (dipEndNow->colType != 0)
599  pT2endDip = max( pT2sel, pT2min );
600  else if (abs(dipEndNow->weakType) != 0)
601  pT2endDip = max( pT2sel, pT2weakCut);
602  else if (abs(dipEndNow->chgType) != 3 && dipEndNow->chgType != 0)
603  pT2endDip = max( pT2sel, pT2minChgQ );
604  else
605  pT2endDip = max( pT2sel, pT2minChgL );
606 
607  // Find properties of dipole and radiating dipole end.
608  sideA = ( abs(dipEndNow->side) == 1 );
609  BeamParticle& beamNow = (sideA) ? *beamAPtr : *beamBPtr;
610  BeamParticle& beamRec = (sideA) ? *beamBPtr : *beamAPtr;
611  iNow = beamNow[iSysNow].iPos();
612  iRec = beamRec[iSysNow].iPos();
613  idDaughter = beamNow[iSysNow].id();
614  xDaughter = beamNow[iSysNow].x();
615  x1Now = (sideA) ? xDaughter : beamRec[iSysNow].x();
616  x2Now = (sideA) ? beamRec[iSysNow].x() : xDaughter;
617 
618  // If reconstructed back to the beam photon, no further ISR emissions.
619  if ( beamNow.isGamma() && !(beamNow.resolvedGamma()) ) continue;
620 
621  // Note dipole mass correction when recoiler is a rescatter.
622  m2Rec = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
623  m2Dip = x1Now * x2Now * sCM + m2Rec;
624 
625  // Prepare kinematics for final-state dipole recoil.
626  m2ColPair = (dipEndNow->iColPartner == 0) ? 0.
627  : m2( event[iNow].p(), event[dipEndNow->iColPartner].p() );
628  mColPartner = (dipEndNow->iColPartner == 0) ? 0.
629  : event[dipEndNow->iColPartner].m();
630  m2ColPartner = pow2(mColPartner);
631 
632  // Stop if m2ColPair is negative.
633  if (m2ColPair < 0.) return 0.;
634 
635  // Now do evolution in pT2, for QCD, QED or weak.
636  if (pT2begDip > pT2endDip) {
637  if (dipEndNow->colType != 0) pT2nextQCD( pT2begDip, pT2endDip);
638  else if (dipEndNow->chgType != 0 || idDaughter == 22)
639  pT2nextQED( pT2begDip, pT2endDip);
640  else if (dipEndNow->weakType != 0) pT2nextWeak( pT2begDip, pT2endDip);
641 
642  // Update if found larger pT than current maximum.
643  if (dipEndNow->pT2 > pT2sel) {
644  pT2sel = dipEndNow->pT2;
645  iDipSel = iDipNow;
646  iSysSel = iSysNow;
647  dipEndSel = dipEndNow;
648  splittingNameSel = splittingNameNow;
649  }
650  }
651  }
652  // End loop over dipole ends.
653  }
654 
655  // Return nonvanishing value if found pT is bigger than already found.
656  return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
657 }
658 
659 //--------------------------------------------------------------------------
660 
661 // Evolve a QCD dipole end.
662 
663 void SpaceShower::pT2nextQCD( double pT2begDip, double pT2endDip) {
664 
665  // Some properties and kinematical starting values.
666  BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
667  bool isGluon = (idDaughter == 21);
668  bool isValence = beam[iSysNow].isValence();
669  int MEtype = dipEndNow->MEtype;
670  int iColPartner = dipEndNow->iColPartner;
671  int idColPartner = dipEndNow->idColPartner;
672  double pT2 = pT2begDip;
673  double xMaxAbs = beam.xMax(iSysNow);
674  double zMinAbs = xDaughter / xMaxAbs;
675  if (xMaxAbs < 0.) {
676  infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
677  "xMaxAbs negative");
678  return;
679  }
680 
681  // Starting values for handling of massive quarks (c/b), if any.
682  double idMassive = 0;
683  if ( abs(idDaughter) == 4 ) idMassive = 4;
684  if ( abs(idDaughter) == 5 ) idMassive = 5;
685  bool isMassive = (idMassive > 0);
686  double m2Massive = 0.;
687  double mRatio = 0.;
688  double zMaxMassive = 1.;
689  double m2Threshold = pT2;
690 
691  // Evolution below scale of massive quark or at large x is impossible.
692  if (isMassive) {
693  m2Massive = (idMassive == 4) ? m2c : m2b;
694  if (pT2 < HEAVYPT2EVOL * m2Massive) return;
695  if (iColPartner == 0) {
696  mRatio = sqrt( m2Massive / m2Dip );
697  zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
698  } else {
699  double m2Red = m2ColPair - m2ColPartner;
700  zMaxMassive = m2Red / (m2Red + m2Massive);
701  }
702  if (xDaughter > HEAVYXEVOL * zMaxMassive * xMaxAbs) return;
703 
704  // Find threshold scale below which only g -> Q + Qbar will be allowed.
705  m2Threshold = (idMassive == 4) ? min( pT2, CTHRESHOLD * m2c)
706  : min( pT2, BTHRESHOLD * m2b);
707  }
708 
709  // Variables used inside evolution loop. (Mainly dummy starting values.)
710  int nFlavour = 3;
711  double b0 = 4.5;
712  double Lambda2 = Lambda3flav2;
713  double pT2minNow = pT2endDip;
714  int idMother = 0;
715  int idSister = 0;
716  double z = 0.;
717  double zMaxAbs = 0.;
718  double zRootMax = 0.;
719  double zRootMin = 0.;
720  double g2gInt = 0.;
721  double q2gInt = 0.;
722  double q2qInt = 0.;
723  double g2qInt = 0.;
724  double g2Qenhance = 0.;
725  double xPDFdaughter = 0.;
726  double xPDFmother[21] = {0.};
727  double xPDFgMother = 0.;
728  double xPDFmotherSum = 0.;
729  double kernelPDF = 0.;
730  double xMother = 0.;
731  double wt = 0.;
732  double Q2 = 0.;
733  double mSister = 0.;
734  double m2Sister = 0.;
735  double pT2corr = 0.;
736  double pT2PDF = pT2;
737  bool needNewPDF = true;
738 
739  // Add more headroom if doing uncertainty variations
740  // (to ensure at least a minimal number of failed branchings).
741  doUncertaintiesNow = doUncertainties;
742  if (!uVarMPIshowers && iSysNow != 0) doUncertaintiesNow = false;
743  double overFac = doUncertaintiesNow ? overFactor : 1.0;
744 
745  // For dipole recoil: other-end colour factor correction in q-g dipole.
746  double coefColRec = (iColPartner != 0 && idColPartner == 21) ? 9./8. : 1.;
747 
748  // Set default values for enhanced emissions.
749  bool isEnhancedQ2QG, isEnhancedG2QQ, isEnhancedQ2GQ, isEnhancedG2GG;
750  isEnhancedQ2QG = isEnhancedG2QQ = isEnhancedQ2GQ = isEnhancedG2GG = false;
751  double enhanceNow = 1.;
752  string nameNow = "";
753 
754  // Begin evolution loop towards smaller pT values.
755  int loopTinyPDFdau = 0;
756  bool hasTinyPDFdau = false;
757  do {
758 
759  // Default values for current tentative emission.
760  wt = 0.;
761  isEnhancedQ2QG = isEnhancedG2QQ = isEnhancedQ2GQ = isEnhancedG2GG = false;
762  enhanceNow = 1.;
763  nameNow = "";
764 
765  // Bad sign if repeated looping with small daughter PDF, so fail.
766  // (Example: if all PDF's = 0 below Q_0, except for c/b companion.)
767  if (hasTinyPDFdau) ++loopTinyPDFdau;
768  if (loopTinyPDFdau > MAXLOOPTINYPDF) {
769  infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
770  "small daughter PDF");
771  return;
772  }
773 
774  // Initialize integrals of splitting kernels and evaluate parton
775  // densities at the beginning. Reinitialize after long evolution
776  // in pT2 or when crossing c and b flavour thresholds.
777  if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
778  pT2PDF = pT2;
779  hasTinyPDFdau = false;
780 
781  // Determine overestimated z range; switch at c and b masses.
782  if (pT2 > m2b) {
783  nFlavour = 5;
784  pT2minNow = max( m2b, pT2endDip);
785  b0 = 23./6.;
786  Lambda2 = Lambda5flav2;
787  } else if (pT2 > m2c) {
788  nFlavour = 4;
789  pT2minNow = max( m2c, pT2endDip);
790  b0 = 25./6.;
791  Lambda2 = Lambda4flav2;
792  } else {
793  nFlavour = 3;
794  pT2minNow = pT2endDip;
795  b0 = 27./6.;
796  Lambda2 = Lambda3flav2;
797  }
798 
799  // A change of renormalization scale expressed by a change of Lambda.
800  Lambda2 /= renormMultFac;
801 
802  // Upper limit on z range: global or local recoil.
803  if (iColPartner == 0) zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip)
804  * ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
805  else {
806  double m2Red = m2ColPair - m2ColPartner;
807  double pT2Ext = sqrt(pT2minNow * (pT2minNow + 4. * m2ColPartner));
808  zMaxAbs = (m2Red + 0.5 * (pT2minNow - pT2Ext))
809  / (m2Red + pT2minNow * (1. - m2ColPartner / m2Red));
810  }
811  if (isMassive) zMaxAbs = min( zMaxAbs, zMaxMassive);
812 
813  // Go to another z range with lower mass scale if current is closed.
814  if (zMinAbs > zMaxAbs) {
815  if (nFlavour == 3 || (idMassive == 4 && nFlavour == 4)
816  || idMassive == 5) return;
817  pT2 = (nFlavour == 4) ? m2c : m2b;
818  continue;
819  }
820 
821  // Parton density of daughter at current scale.
822  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
823  xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2);
824  if (xPDFdaughter < TINYPDF) {
825  xPDFdaughter = TINYPDF;
826  hasTinyPDFdau = true;
827  }
828 
829  // Integrals of splitting kernels for gluons: g -> g, q -> g.
830  if (isGluon) {
831  g2gInt = overFac * HEADROOMG2G * 6.
832  * log(zMaxAbs * (1.-zMinAbs) / (zMinAbs * (1.-zMaxAbs)));
833  if (doMEcorrections) g2gInt *= calcMEmax(MEtype, 21, 21);
834  // Optionally enhanced branching rate.
835  if (canEnhanceET) g2gInt *= userHooksPtr->enhanceFactor("isr:G2GG");
836  q2gInt = overFac * HEADROOMQ2G * (16./3.)
837  * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
838  if (doMEcorrections) q2gInt *= calcMEmax(MEtype, 1, 21);
839  // Optionally enhanced branching rate.
840  if (canEnhanceET) q2gInt *= userHooksPtr->enhanceFactor("isr:Q2GQ");
841 
842  // Parton density of potential quark mothers to a g.
843  xPDFmotherSum = 0.;
844  for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
845  if (i == 0) {
846  xPDFmother[10] = 0.;
847  } else {
848  xPDFmother[i+10] = beam.xfISR(iSysNow, i, xDaughter, pdfScale2);
849  xPDFmotherSum += xPDFmother[i+10];
850  }
851  }
852 
853  // Total QCD evolution coefficient for a gluon.
854  kernelPDF = g2gInt + q2gInt * xPDFmotherSum / xPDFdaughter;
855 
856  // For valence quark only need consider q -> q g branchings.
857  // Introduce an extra factor sqrt(z) to smooth bumps.
858  } else if (isValence) {
859  zRootMin = (1. + sqrt(zMinAbs)) / (1. - sqrt(zMinAbs));
860  zRootMax = (1. + sqrt(zMaxAbs)) / (1. - sqrt(zMaxAbs));
861  q2qInt = coefColRec * overFac * (8./3.) * log( zRootMax / zRootMin );
862  if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
863  // Optionally enhanced branching rate.
864  if (canEnhanceET) q2qInt *= userHooksPtr->enhanceFactor("isr:Q2QG");
865  kernelPDF = q2qInt;
866 
867  // Integrals of splitting kernels for quarks: q -> q, g -> q.
868  } else {
869  q2qInt = coefColRec * overFac * HEADROOMQ2Q * (8./3.)
870  * log( (1. - zMinAbs) / (1. - zMaxAbs) );
871  if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
872  // Optionally enhanced branching rate.
873  if (canEnhanceET) q2qInt *= userHooksPtr->enhanceFactor("isr:Q2QG");
874  g2qInt = overFac * HEADROOMG2Q * 0.5 * (zMaxAbs - zMinAbs);
875  if (doMEcorrections) g2qInt *= calcMEmax(MEtype, 21, 1);
876  // Optionally enhanced branching rate.
877  if (canEnhanceET) g2qInt *= userHooksPtr->enhanceFactor("isr:G2QQ");
878 
879  // Increase the upper weight for heavy quarks in photon beam
880  // due to different behavior of the PDFs.
881  if (beam.isGamma() && isMassive) q2qInt *= HEADROOMHQG;
882 
883  // Increase estimated upper weight for g -> Q + Qbar.
884  if (isMassive) {
885  if (alphaSorder == 0) g2Qenhance = log(pT2/m2Massive)
886  / log(m2Threshold/m2Massive);
887  else {
888  double m2log = log( m2Massive / Lambda2);
889  g2Qenhance = log( log(pT2/Lambda2) / m2log )
890  / log( log(m2Threshold/Lambda2) / m2log );
891  }
892  g2qInt *= g2Qenhance;
893  }
894 
895  // Parton density of a potential gluon mother to a q.
896  xPDFgMother = beam.xfISR(iSysNow, 21, xDaughter, pdfScale2);
897 
898  // Total QCD evolution coefficient for a quark.
899  kernelPDF = q2qInt + g2qInt * xPDFgMother / xPDFdaughter;
900  }
901 
902  // End evaluation of splitting kernels and parton densities.
903  needNewPDF = false;
904  }
905  if (kernelPDF < TINYKERNELPDF) return;
906 
907  // Pick pT2 (in overestimated z range), for one of three different cases.
908  // Assume form alphas(pT0^2 + pT^2) * dpT^2/(pT0^2 + pT^2).
909  double Q2alphaS;
910 
911  // Fixed alpha_strong.
912  if (alphaSorder == 0) {
913  pT2 = (pT2 + pT20) * pow( rndmPtr->flat(),
914  1. / (alphaS2pi * kernelPDF)) - pT20;
915 
916  // First-order alpha_strong.
917  } else if (alphaSorder == 1) {
918  pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
919  pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
920 
921  // For second order reject by second term in alpha_strong expression.
922  } else {
923  do {
924  pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
925  pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
926  Q2alphaS = renormMultFac * max( pT2 + pT20,
927  pow2(LAMBDA3MARGIN) * Lambda3flav2);
928  } while (alphaS.alphaS2OrdCorr(Q2alphaS) < rndmPtr->flat()
929  && pT2 > pT2minNow);
930  }
931 
932  // Check for pT2 values that prompt special action.
933 
934  // If fallen into b threshold region, force g -> b + bbar.
935  if (idMassive == 5 && pT2 < m2Threshold) {
936  pT2nearThreshold( beam, m2Massive, m2Threshold, xMaxAbs,
937  zMinAbs, zMaxMassive, iColPartner );
938  return;
939 
940  // If crossed b threshold, continue evolution from this threshold.
941  } else if (nFlavour == 5 && pT2 < m2b) {
942  needNewPDF = true;
943  pT2 = m2b;
944  continue;
945 
946  // If fallen into c threshold region, force g -> c + cbar.
947  } else if (idMassive == 4 && pT2 < m2Threshold) {
948  pT2nearThreshold( beam, m2Massive, m2Threshold, xMaxAbs,
949  zMinAbs, zMaxMassive, iColPartner );
950  return;
951 
952  // If crossed c threshold, continue evolution from this threshold.
953  } else if (nFlavour == 4 && pT2 < m2c) {
954  needNewPDF = true;
955  pT2 = m2c;
956  continue;
957 
958  // Abort evolution if below cutoff scale, or below another branching.
959  } else if (pT2 < pT2endDip) return;
960 
961  // Select z value of branching to g, and corrective weight.
962  if (isGluon) {
963 
964  // g -> g (+ g).
965  if (rndmPtr->flat() * kernelPDF < g2gInt) {
966  idMother = 21;
967  idSister = 21;
968  z = 1. / ( 1. + ((1. - zMinAbs) / zMinAbs) * pow( (zMinAbs *
969  (1. - zMaxAbs)) / (zMaxAbs * (1. - zMinAbs)), rndmPtr->flat() ) );
970  wt = pow2( 1. - z * (1. - z));
971  // For dipole recoil: other-end colour factor correction in g-q dipole.
972  if (iColPartner != 0 && idColPartner != 21)
973  wt *= (m2ColPair * pow2(1. - z) + z * pT2 * 8./9.)
974  / (m2ColPair * pow2(1. - z) + z * pT2);
975  // Account for headroom factor used to enhance trial probability
976  wt /= HEADROOMG2G;
977  // Optionally enhanced branching rate.
978  nameNow = "isr:G2GG";
979  if (canEnhanceET) {
980  double enhance = userHooksPtr->enhanceFactor(nameNow);
981  if (enhance != 1.) {
982  enhanceNow = enhance;
983  isEnhancedG2GG = true;
984  }
985  }
986 
987  // q -> g (+ q): also select flavour.
988  } else {
989  double temp = xPDFmotherSum * rndmPtr->flat();
990  idMother = -nQuarkIn - 1;
991  do { temp -= xPDFmother[(++idMother) + 10]; }
992  while (temp > 0. && idMother < nQuarkIn);
993  idSister = idMother;
994  z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
995  * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
996  wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z)
997  * xPDFdaughter / xPDFmother[idMother + 10];
998  // Account for headroom factor used to enhance trial probability
999  wt /= HEADROOMQ2G;
1000  // Optionally enhanced branching rate.
1001  nameNow = "isr:Q2GQ";
1002  if (canEnhanceET) {
1003  double enhance = userHooksPtr->enhanceFactor(nameNow);
1004  if (enhance != 1.) {
1005  enhanceNow = enhance;
1006  isEnhancedQ2GQ = true;
1007  }
1008  }
1009  }
1010 
1011  // Select z value of branching to q, and corrective weight.
1012  // Include massive kernel corrections for c and b quarks.
1013  } else {
1014 
1015  // q -> q (+ g).
1016  if (isValence || rndmPtr->flat() * kernelPDF < q2qInt) {
1017  idMother = idDaughter;
1018  idSister = 21;
1019  // Valence more peaked at large z.
1020  if (isValence) {
1021  double zTmp = zRootMin * pow(zRootMax / zRootMin, rndmPtr->flat() );
1022  z = pow2( (1. - zTmp) / (1. + zTmp) );
1023  } else {
1024  z = 1. - (1. - zMinAbs) * pow( (1. - zMaxAbs) / (1. - zMinAbs),
1025  rndmPtr->flat() );
1026  }
1027  if (!isMassive) {
1028  wt = 0.5 * (1. + pow2(z));
1029  } else {
1030  wt = 0.5 * (1. + pow2(z)) - z * pow2(1.-z) * m2Massive / pT2;
1031  }
1032  if (isValence) wt *= sqrt(z);
1033  // Account for headroom factor for sea quarks
1034  else wt /= HEADROOMQ2Q;
1035  // Account for headroom factor for heavy quarks in photon beam.
1036  if (beam.isGamma() && isMassive) wt /= HEADROOMHQG;
1037  // For dipole recoil: other-end colour factor correction in q-g dipole.
1038  if (iColPartner != 0 && idColPartner == 21)
1039  wt *= (m2ColPair * pow2(1. - z) + z * pT2 * 9./8.)
1040  / ((m2ColPair * pow2(1. - z) + z * pT2) * coefColRec);
1041  // Optionally enhanced branching rate.
1042  nameNow = "isr:Q2QG";
1043  if (canEnhanceET) {
1044  double enhance = userHooksPtr->enhanceFactor(nameNow);
1045  if (enhance != 1.) {
1046  enhanceNow = enhance;
1047  isEnhancedQ2QG = true;
1048  }
1049  }
1050 
1051  // g -> q (+ qbar).
1052  } else {
1053  idMother = 21;
1054  idSister = - idDaughter;
1055  z = zMinAbs + rndmPtr->flat() * (zMaxAbs - zMinAbs);
1056  if (!isMassive) {
1057  wt = (pow2(z) + pow2(1.-z)) * xPDFdaughter / xPDFgMother;
1058  } else {
1059  wt = (pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2)
1060  * xPDFdaughter / (xPDFgMother * g2Qenhance) ;
1061  }
1062  // Account for headroom factor for gluons
1063  wt /= HEADROOMG2Q;
1064  // Optionally enhanced branching rate.
1065  if (abs(idSister) < 4) nameNow = "isr:G2QQ";
1066  else if (abs(idSister) == 4) nameNow = "isr:G2QQ:cc";
1067  else nameNow = "isr:G2QQ:bb";
1068  if (canEnhanceET) {
1069  double enhance = userHooksPtr->enhanceFactor(nameNow);
1070  if (enhance != 1.) {
1071  enhanceNow = enhance;
1072  isEnhancedG2QQ = true;
1073  }
1074  }
1075  }
1076  }
1077 
1078  // Cancel out uncertainty-band extra headroom factors.
1079  wt /= overFac;
1080 
1081  // Derive Q2 and x of mother from pT2 and z.
1082  Q2 = pT2 / (1. - z);
1083  xMother = xDaughter / z;
1084 
1085  // Correction to x for massive recoiler from rescattering.
1086  if (!dipEndNow->normalRecoil) {
1087  if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1088  else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1089  }
1090  if (xMother > xMaxAbs) { wt = 0.; continue; }
1091 
1092  // Forbidden emission if outside allowed z range for given pT2.
1093  mSister = particleDataPtr->m0(idSister);
1094  m2Sister = pow2(mSister);
1095  if (iColPartner == 0) {
1096  pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1097  } else {
1098  double tmpRat = z * (Q2 + m2Sister) / (m2ColPair - m2ColPartner);
1099  pT2corr = ((1. - z) * Q2 - z * m2Sister) * (1. - tmpRat)
1100  - m2ColPartner * pow2(tmpRat);
1101  }
1102  if (pT2corr < TINYPT2) { wt = 0.; continue; }
1103 
1104  // For emissions in the hard scattering system, optionally veto
1105  // emissions not ordered in rapidity (= angle).
1106  if ( iSysNow == 0 && doRapidityOrder && dipEndNow->nBranch > 0
1107  && pT2 > pow2( (1. - z) / (z * (1. - dipEndNow->zOld)) )
1108  * dipEndNow->pT2Old ) { wt = 0.; continue; }
1109 
1110  // For emissions in any secondary scattering system, optionally veto
1111  // emissions not ordered in rapidity (= angle).
1112  if ( iSysNow != 0 && doRapidityOrderMPI && dipEndNow->nBranch > 0
1113  && pT2 > pow2( (1. - z) / (z * (1. - dipEndNow->zOld)) )
1114  * dipEndNow->pT2Old ) { wt = 0.; continue; }
1115 
1116  // If creating heavy quark by Q -> g + Q then next need g -> Q + Qbar.
1117  // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
1118  if ( isGluon && ( abs(idMother) == 4 || abs(idMother) == 5 )) {
1119  double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
1120  double pT2QQcorr = 0.;
1121  if (iColPartner == 0) {
1122  pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
1123  } else {
1124  double tmpRat = z * (Q2 + m2QQsister) / (m2ColPair - m2ColPartner);
1125  pT2QQcorr = ((1. - z) * Q2 - z * m2QQsister) * (1. - tmpRat)
1126  - m2ColPartner * pow2(tmpRat);
1127  }
1128  if (pT2QQcorr < TINYPT2) { wt = 0.; continue; }
1129  }
1130 
1131  // Check that room left for beam remnants with photon beam.
1132  if ( beam.isGamma() ) {
1133  BeamParticle& beamOther = (sideA) ? *beamBPtr : *beamAPtr;
1134  if ( !beamOther.resolvedGamma() ) {
1135  // Check that room left for 1 beam remnant if other fixed
1136  if ( !beam.roomFor1Remnant( idMother, xMother, eCM) ) {
1137  wt = 0.;
1138  continue;
1139  }
1140  // Otherwise check that room left for 2 beam remnants.
1141  } else {
1142  if ( !beamOther.roomFor2Remnants( idMother, xMother, eCM ) ) {
1143  wt = 0.;
1144  continue;
1145  }
1146  }
1147  }
1148 
1149  // Evaluation of ME correction.
1150  if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
1151  m2Dip, z, Q2, m2Sister) / calcMEmax(MEtype, idMother, idDaughter);
1152 
1153  // Optional dampening of large pT values in first radiation.
1154  if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
1155  wt *= pT2damp / (pT2 + pT2damp);
1156 
1157  // Idea suggested by Gosta Gustafson: increased screening in events
1158  // with large activity can be simulated by pT0_eff = sqrt(n) * pT0.
1159  if (enhanceScreening == 2) {
1160  int nSysNow = infoPtr->nMPI() + infoPtr->nISR() + 1;
1161  double WTscreen = pow2( (pT2 + pT20) / (pT2 + nSysNow * pT20) );
1162  wt *= WTscreen;
1163  }
1164 
1165  // Evaluation of new daughter and mother PDF's.
1166  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1167  double xPDFdaughterNew = max ( TINYPDF,
1168  beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
1169  double xPDFmotherNew =
1170  beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
1171  wt *= xPDFmotherNew / xPDFdaughterNew;
1172 
1173  // If doing uncertainty variations, postpone accept/reject to branch()
1174  if (wt > 0. && pT2 > pT2min && doUncertaintiesNow ) {
1175  dipEndNow->pAccept = wt;
1176  wt = 1.0;
1177  }
1178 
1179  // Check that valence step does not cause problem.
1180  if (wt > 1. && pT2 > PT2MINWARN) infoPtr->errorMsg("Warning in "
1181  "SpaceShower::pT2nextQCD: weight above unity");
1182 
1183  // Iterate until acceptable pT (or have fallen below pTmin).
1184  } while (wt < rndmPtr->flat()) ;
1185 
1186  // Store outcome of enhanced branching rate analysis.
1187  splittingNameNow = nameNow;
1188  if (canEnhanceET) {
1189  if (isEnhancedQ2QG) storeEnhanceFactor(pT2,"isr:Q2QG", enhanceNow);
1190  if (isEnhancedG2QQ) storeEnhanceFactor(pT2,"isr:G2QQ", enhanceNow);
1191  if (isEnhancedQ2GQ) storeEnhanceFactor(pT2,"isr:Q2GQ", enhanceNow);
1192  if (isEnhancedG2GG) storeEnhanceFactor(pT2,"isr:G2GG", enhanceNow);
1193  }
1194 
1195  // Save values for (so far) acceptable branching.
1196  dipEndNow->store( idDaughter,idMother, idSister, x1Now, x2Now, m2Dip,
1197  pT2, z, xMother, Q2, mSister, m2Sister, pT2corr, iColPartner, m2ColPair,
1198  mColPartner);
1199 
1200 }
1201 
1202 //--------------------------------------------------------------------------
1203 
1204 // Evolve a QCD dipole end near threshold, with g -> Q + Qbar enforced.
1205 // Note: No explicit Sudakov factor formalism here. Instead use that
1206 // df_Q(x, pT2) = (alpha_s/2pi) * (dT2/pT2) * ((gluon) * (splitting)).
1207 // This implies that effects of Q -> Q + g are neglected in this range.
1208 // Now includes also threshold behaviour with photon beams.
1209 
1210 void SpaceShower::pT2nearThreshold( BeamParticle& beam,
1211  double m2Massive, double m2Threshold, double xMaxAbs,
1212  double zMinAbs, double zMaxMassive, int iColPartner) {
1213 
1214  // Initial values, to be used in kinematics and weighting.
1215  double Lambda2 = (abs(idDaughter) == 4) ? Lambda4flav2 : Lambda5flav2;
1216  Lambda2 /= renormMultFac;
1217  double logM2Lambda2 = (alphaSorder > 0) ? log( m2Massive / Lambda2 ) : 1.;
1218  pdfScale2 = (useFixedFacScale) ? fixedFacScale2
1219  : factorMultFac * m2Threshold;
1220  double xPDFmotherOld = beam.xfISR(iSysNow, 21, xDaughter, pdfScale2);
1221  // Check that xPDF is not vanishing.
1222  if ( xPDFmotherOld < TINYPDF ) {
1223  infoPtr->errorMsg("Error in SpaceShower::pT2nearThreshold: xPDF = 0");
1224  return;
1225  }
1226 
1227  // Variables used inside evolution loop. (Mainly dummy start values.)
1228  int loop = 0;
1229  double wt = 0.;
1230  double pT2 = 0.;
1231  double z = 0.;
1232  double Q2 = 0.;
1233  double pT2corr = 0.;
1234  double xMother = 0.;
1235 
1236  // Check if photon beam is being evolved.
1237  bool isGammaBeam = beam.isGamma();
1238  if( isGammaBeam ){
1239  BeamParticle& beamOther = (sideA) ? *beamBPtr : *beamAPtr;
1240  // Allow splitting only if room for remnants in the other side.
1241  if ( !beamOther.roomFor1Remnant(eCM) ) return;
1242  }
1243 
1244  // Begin loop over tries to find acceptable g -> Q + Qbar branching.
1245  do {
1246  wt = 0.;
1247 
1248  // Check that not caught in infinite loop with impossible kinematics.
1249  if (++loop > 100) {
1250  infoPtr->errorMsg("Error in SpaceShower::pT2nearThreshold: "
1251  "stuck in loop");
1252  return;
1253  }
1254 
1255  // Pick dpT2/pT2 in range [m2Massive,thresholdRatio * m2Massive].
1256  pT2 = m2Massive * pow( m2Threshold / m2Massive, rndmPtr->flat() );
1257 
1258  // For photon beams kinematics are fixed.
1259  if (isGammaBeam) {
1260  xMother = 1.0;
1261  z = xDaughter/xMother;
1262  // Pick z flat in allowed range.
1263  } else {
1264  z = zMinAbs + rndmPtr->flat() * (zMaxMassive - zMinAbs);
1265  }
1266 
1267  // Check that kinematically possible choice.
1268  Q2 = pT2 / (1.-z) - m2Massive;
1269  if (iColPartner == 0) {
1270  pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Massive) / m2Dip;
1271  } else {
1272  double tmpRat = z * (Q2 + m2Massive) / (m2ColPair - m2ColPartner);
1273  pT2corr = ((1. - z) * Q2 - z * m2Massive) * (1. - tmpRat)
1274  - m2ColPartner * pow2(tmpRat);
1275  }
1276  if (pT2corr < TINYPT2) continue;
1277 
1278  // Correction factor for splitting kernel.
1279  wt = pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2;
1280 
1281  // Sample the kinematics for hadron beams.
1282  if (!isGammaBeam) {
1283 
1284  // Correction factor for running alpha_s. (Only first order for now.)
1285  wt *= (alphaSorder > 0) ? logM2Lambda2 / log( pT2 / Lambda2 ) : 1.;
1286 
1287  // x, including correction for massive recoiler from rescattering.
1288  xMother = xDaughter / z;
1289  if (!dipEndNow->normalRecoil) {
1290  if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1291  else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1292  }
1293  if (xMother > xMaxAbs) { wt = 0.; continue; }
1294 
1295  // Correction factor for gluon density.
1296  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1297  double xPDFmotherNew = beam.xfISR(iSysNow, 21, xMother, pdfScale2);
1298  wt *= xPDFmotherNew / xPDFmotherOld;
1299 
1300  }
1301 
1302  // If doing uncertainty variations, postpone accept/reject to branch().
1303  if (wt > 0. && pT2 > pT2min && doUncertaintiesNow ) {
1304  dipEndNow->pAccept = wt;
1305  wt = 1.0;
1306  }
1307 
1308  // Iterate until acceptable pT and z.
1309  } while (wt < rndmPtr->flat()) ;
1310 
1311  // Select correct mother for the splitting.
1312  int idMother = isGammaBeam ? 22 : 21;
1313 
1314  // Save values for (so far) acceptable branching.
1315  double mSister = (abs(idDaughter) == 4) ? mc : mb;
1316 
1317  if ( isGammaBeam ) splittingNameNow = "isr:A2QQ";
1318  else splittingNameNow = "isr:G2QQ";
1319  dipEndNow->store( idDaughter, idMother, -idDaughter, x1Now, x2Now, m2Dip,
1320  pT2, z, xMother, Q2, mSister, pow2(mSister), pT2corr, iColPartner,
1321  m2ColPair, mColPartner);
1322 
1323 }
1324 
1325 //--------------------------------------------------------------------------
1326 
1327 // Evolve a QED dipole end.
1328 
1329 void SpaceShower::pT2nextQED( double pT2begDip, double pT2endDip) {
1330 
1331  // Type of dipole and starting values.
1332  BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
1333  bool isLeptonBeam = beam.isLepton();
1334  bool isGammaBeam = beam.isGamma();
1335  int MEtype = dipEndNow->MEtype;
1336  bool isPhoton = (idDaughter == 22);
1337  double pT2 = pT2begDip;
1338  double m2Lepton = (isLeptonBeam) ? pow2(beam.m()) : 0.;
1339  if (isLeptonBeam && pT2begDip < m2Lepton) return;
1340 
1341  // Currently no f -> gamma branching implemented for lepton or photon beams.
1342  if (isPhoton && (isLeptonBeam || isGammaBeam) ) return;
1343 
1344  // alpha_em at maximum scale provides upper estimate.
1345  double alphaEMmax = alphaEM.alphaEM(renormMultFac * pT2begDip);
1346  double alphaEM2pi = alphaEMmax / (2. * M_PI);
1347 
1348  // Maximum x of mother implies minimum z = xDaughter / xMother.
1349  double xMaxAbs = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
1350  double zMinAbs = xDaughter / xMaxAbs;
1351  if (xMaxAbs < 0.) {
1352  infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
1353  "xMaxAbs negative");
1354  return;
1355  }
1356 
1357  // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
1358  double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
1359  ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
1360  if (isLeptonBeam) {
1361  double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
1362  if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
1363  }
1364  if (zMaxAbs < zMinAbs) return;
1365 
1366  // Similar threshold behaviour as in QCD evolution for photon beams.
1367  double idMassive = 0;
1368  if ( isGammaBeam && abs(idDaughter) == 4 ) idMassive = 4;
1369  if ( isGammaBeam && abs(idDaughter) == 5 ) idMassive = 5;
1370  bool isMassive = (idMassive > 0);
1371  double m2Massive = 0.;
1372  double mRatio = 0.;
1373  double zMaxMassive = 1.;
1374  double m2Threshold = pT2;
1375 
1376  // Evolution below scale of massive quark or at large x is impossible.
1377  if (isMassive) {
1378  m2Massive = (idMassive == 4) ? m2c : m2b;
1379  if (pT2 < HEAVYPT2EVOL * m2Massive) return;
1380  mRatio = sqrt( m2Massive / m2Dip );
1381  zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
1382  if (xDaughter > HEAVYXEVOL * zMaxMassive * xMaxAbs) return;
1383 
1384  // Find threshold scale below which only gamma -> Q + Qbar will be allowed.
1385  m2Threshold = (idMassive == 4) ? min( pT2, CTHRESHOLD * m2c)
1386  : min( pT2, BTHRESHOLD * m2b);
1387  }
1388 
1389  // Variables used inside evolution loop. (Mainly dummy start values.)
1390  int idMother = 0;
1391  int idSister = 22;
1392  double z = 0.;
1393  double xMother = 0.;
1394  double wt = 0.;
1395  double Q2 = 0.;
1396  double mSister = 0.;
1397  double m2Sister = 0.;
1398  double pT2corr = 0.;
1399 
1400  // Set default values for enhanced emissions.
1401  bool isEnhancedQ2QA, isEnhancedQ2AQ, isEnhancedA2QQ;
1402  isEnhancedQ2QA = isEnhancedQ2AQ = isEnhancedA2QQ = false;
1403  double enhanceNow = 1.;
1404  string nameNow = "";
1405 
1406  // QED evolution of fermions
1407  if (!isPhoton) {
1408 
1409  // Integrals of splitting kernels for fermions: f -> f. Use 1 + z^2 < 2.
1410  // Ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
1411  double f2fInt = 0.;
1412  double f2fIntA = 2. * log( (1. - zMinAbs) / (1. - zMaxAbs) );
1413  double f2fIntB = 0.;
1414  if (isLeptonBeam) {
1415  f2fIntB = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
1416  f2fInt = f2fIntA + f2fIntB;
1417  } else f2fInt = pow2(dipEndNow->chgType / 3.) * f2fIntA;
1418 
1419  // gamma -> q qbar splitting where gamma is the original beam particle.
1420  // P(z) = 3e_q^2(z^2 + (1-z)^2)
1421  // Correct for the weight function (currently constant*log(pT2/pT2ref)).
1422  double gamma2f = 0;
1423  double pT2ref = 0;
1424  double xfApprox = 0;
1425  if (isGammaBeam) {
1426 
1427  // For photon beams approximate PDFs to estimate the splitting
1428  // probability.
1429  pT2ref = beam.gammaPDFRefScale(idDaughter);
1430  xfApprox = beam.gammaPDFxDependence(idDaughter, xDaughter);
1431 
1432  // Allow splitting only if room for remnants at the other side.
1433  BeamParticle& beamOther = (sideA) ? *beamBPtr : *beamAPtr;
1434  if ( beamOther.roomFor1Remnant(eCM) ) {
1435  gamma2f = alphaEM2pi * pow2(double(dipEndNow->chgType) / 3.)* 3.
1436  * (pow2(xDaughter) + pow2(1. - xDaughter))/xfApprox;
1437  }
1438  }
1439 
1440  // Upper estimate for evolution equation, including fudge factor.
1441  if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
1442  double kernelPDF = alphaEM2pi * f2fInt;
1443  double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
1444  kernelPDF *= fudge;
1445  // Check the kernelPDF value before possible enhancement but do not
1446  // enhance gamma -> q qbar splittings.
1447  if ( (kernelPDF + gamma2f) < TINYKERNELPDF ) return;
1448 
1449  // Optionally enhanced branching rate.
1450  if (canEnhanceET) kernelPDF *= userHooksPtr->enhanceFactor("isr:Q2QA");
1451 
1452  // Optionally enhanced branching rate.
1453  if (canEnhanceET) gamma2f *= userHooksPtr->enhanceFactor("isr:A2QQ");
1454 
1455  // Add gamma -> q qbar splittings to kernelPDF for photon beam.
1456  kernelPDF += gamma2f;
1457 
1458  // Begin evolution loop towards smaller pT values.
1459  do {
1460 
1461  // Default values for current tentative emission.
1462  isEnhancedQ2QA = isEnhancedA2QQ = false;
1463  enhanceNow = 1.;
1464  nameNow = "";
1465 
1466  // gamma -> f fbar splitting with photon beam.
1467  if( (rndmPtr->flat() * kernelPDF) < gamma2f ){
1468 
1469  // Sample pT2 using pT2ref.
1470  pT2 = pT2ref*pow(pT2/pT2ref, pow(rndmPtr->flat(), 1. / kernelPDF));
1471 
1472  // If fallen into b or c threshold region, force gamma -> q + qbar.
1473  if ( isMassive && (pT2 < m2Threshold ) ) {
1474  pT2nearThreshold( beam, m2Massive, m2Threshold, xMaxAbs,
1475  zMinAbs, zMaxMassive, 0 );
1476  return;
1477  } else if (pT2 < pT2endDip) return;
1478 
1479  // Fix the flavors and masses for the partons.
1480  idMother = 22;
1481  xMother = 1.0;
1482  z = xDaughter/xMother;
1483  idSister = -idDaughter;
1484  mSister = particleDataPtr->m0(idSister);
1485  m2Sister = pow2(mSister);
1486 
1487  // Derive Q2 from pT2 and z.
1488  Q2 = pT2 / (1. - z);
1489 
1490  // Forbidden emission if outside allowed z range for given pT2.
1491  pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1492  if (pT2corr < TINYPT2) { wt = 0.; continue; }
1493 
1494  // Weight with the x*log(pT2/pT2ref)/(x*pdf).
1495  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1496  double daughterPDF = max ( TINYPDF,
1497  beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
1498  wt = xDaughter*log(pT2/pT2ref)/daughterPDF;
1499 
1500  // Correct for the x-dependence of PDFs
1501  // (Currently only constant depending on the quark flavor).
1502  wt *= xfApprox;
1503 
1504  // Correct to current value of alpha_EM.
1505  double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
1506  wt *= (alphaEMnow / alphaEMmax);
1507 
1508  // Optionally enhanced branching rate.
1509  nameNow = "isr:A2QQ";
1510  if (canEnhanceET) {
1511  double enhance = userHooksPtr->enhanceFactor(nameNow);
1512  if (enhance != 1.) {
1513  enhanceNow = enhance;
1514  isEnhancedA2QQ = true;
1515  }
1516  }
1517 
1518  // Check that gamma -> q qbar step does not cause problem.
1519  if (wt > 1. && pT2 > PT2MINWARN){
1520  infoPtr->errorMsg("Warning in "
1521  "SpaceShower::pT2nextQED: weight above unity");
1522  }
1523 
1524  // f -> f gamma splittings
1525  } else {
1526 
1527  // Pick pT2 (in overestimated z range).
1528  // For l -> l gamma include extrafactor 1 / ln(pT2/m2l) in evolution.
1529  double shift = pow(rndmPtr->flat(), 1. / kernelPDF);
1530  if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
1531  else pT2 = pT2 * shift;
1532 
1533  // If fallen into b or c threshold region, force gamma -> Q + Qbar.
1534  if ( isGammaBeam && isMassive && (pT2 < m2Threshold ) ) {
1535  pT2nearThreshold( beam, m2Massive, m2Threshold, xMaxAbs,
1536  zMinAbs, zMaxMassive, 0 );
1537  return;
1538  }
1539 
1540  // Abort evolution if below cutoff scale, or below another branching.
1541  if (pT2 < pT2endDip) return;
1542  if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
1543 
1544  // Set the IDs for f -> f + gamma splittings.
1545  idSister = 22;
1546  idMother = idDaughter;
1547 
1548  // Select z value of branching f -> f + gamma, and corrective weight.
1549  wt = 1.;
1550  if (isLeptonBeam) {
1551  if (f2fIntA > rndmPtr->flat() * (f2fIntA + f2fIntB)) {
1552  z = 1. - (1. - zMinAbs)
1553  * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
1554  } else {
1555  z = xDaughter + (zMinAbs - xDaughter)
1556  * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter),
1557  rndmPtr->flat() );
1558  }
1559  wt *= (z - xDaughter) / (1. - xDaughter);
1560  } else {
1561  z = 1. - (1. - zMinAbs)
1562  * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
1563  }
1564 
1565  // The same mass corrections as for QCD q -> q g
1566  if (!isMassive) {
1567  wt *= 0.5 * (1. + pow2(z));
1568  } else {
1569  wt *= 0.5 * (1. + pow2(z)) - z * pow2(1.-z) * m2Massive / pT2;
1570  }
1571 
1572  // Optionally enhanced branching rate.
1573  nameNow = "isr:Q2QA";
1574  if (canEnhanceET) {
1575  double enhance = userHooksPtr->enhanceFactor(nameNow);
1576  if (enhance != 1.) {
1577  enhanceNow = enhance;
1578  isEnhancedQ2QA = true;
1579  }
1580  }
1581 
1582  // Derive Q2 and x of mother from pT2 and z.
1583  Q2 = pT2 / (1. - z);
1584  xMother = xDaughter / z;
1585  // Correction to x for massive recoiler from rescattering.
1586  if (!dipEndNow->normalRecoil) {
1587  if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1588  else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1589  }
1590  if (xMother > xMaxAbs) { wt = 0.; continue; }
1591 
1592  // Forbidden emission if outside allowed z range for given pT2.
1593  mSister = 0.;
1594  m2Sister = 0.;
1595  pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1596  if (pT2corr < TINYPT2) { wt = 0.; continue; }
1597 
1598  // Check that room left for beam remnants with photon beams.
1599  if ( beam.isGamma() ) {
1600  BeamParticle& beamOther = (sideA) ? *beamBPtr : *beamAPtr;
1601  // Room for one remnant, other already fixed by ISR.
1602  if( !beamOther.resolvedGamma() ){
1603  if ( !beam.roomFor1Remnant( idMother, xMother, eCM) ) {
1604  wt = 0.;
1605  continue;
1606  }
1607  // Room left for two beam remnants.
1608  } else {
1609  if( !beamOther.roomFor2Remnants( idMother, xMother, eCM ) ) {
1610  wt = 0.;
1611  continue;
1612  }
1613  }
1614  }
1615 
1616  // Correct by ln(pT2 / m2l) and fudge factor.
1617  if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
1618 
1619  // Evaluation of ME correction.
1620  if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
1621  m2Dip, z, Q2, m2Sister) / calcMEmax(MEtype, idMother, idDaughter);
1622 
1623  // Extra QED correction for f fbar -> W+- gamma. Debug??
1624  if (doMEcorrections && MEtype == 1 && idDaughter == idMother
1625  && ( (iSysNow == 0 && idResFirst == 24)
1626  || (iSysNow == 1 && idResSecond == 24) ) ) {
1627  double tHe = -Q2;
1628  double uHe = Q2 - m2Dip * (1. - z) / z;
1629  double chg1 = abs(dipEndNow->chgType / 3.);
1630  double chg2 = 1. - chg1;
1631  wt *= pow2(chg1 * uHe - chg2 * tHe)
1632  / ( (tHe + uHe) * (pow2(chg1) * uHe + pow2(chg2) * tHe) );
1633  }
1634 
1635  // Optional dampening of large pT values in first radiation.
1636  if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
1637  wt *= pT2damp / (pT2 + pT2damp);
1638 
1639  // Correct to current value of alpha_EM.
1640  double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
1641  wt *= (alphaEMnow / alphaEMmax);
1642 
1643  // Evaluation of new daughter and mother PDF's.
1644  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1645  double xPDFdaughterNew = max ( TINYPDF,
1646  beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
1647  double xPDFmotherNew =
1648  beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
1649  wt *= xPDFmotherNew / xPDFdaughterNew;
1650  }
1651 
1652  // Iterate until acceptable pT (or have fallen below pTmin).
1653  } while (wt < rndmPtr->flat()) ;
1654 
1655  }
1656 
1657  // QED evolution of photons (so far only for hadron beams).
1658  else {
1659 
1660  // Initial values
1661  int nFlavour = 3;
1662  double kernelPDF = 0.0;
1663  double xPDFdaughter = 0.;
1664  double xPDFmother[21] = {0.};
1665  double xPDFmotherSum = 0.0;
1666  double pT2PDF = pT2;
1667  double pT2minNow = pT2endDip;
1668  bool needNewPDF = true;
1669 
1670  // Begin evolution loop towards smaller pT values.
1671  int loopTinyPDFdau = 0;
1672  bool hasTinyPDFdau = false;
1673  do {
1674 
1675  // Default values for current tentative emission.
1676  wt = 0.;
1677  isEnhancedQ2AQ = false;
1678  enhanceNow = 1.;
1679  nameNow = "";
1680 
1681  // Bad sign if repeated looping with small daughter PDF, so fail.
1682  if (hasTinyPDFdau) ++loopTinyPDFdau;
1683  if (loopTinyPDFdau > MAXLOOPTINYPDF) {
1684  infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
1685  "small daughter PDF");
1686  return;
1687  }
1688 
1689  // Initialize integrals of splitting kernels and evaluate parton
1690  // densities at the beginning. Reinitialize after long evolution
1691  // in pT2 or when crossing c and b flavour thresholds.
1692  if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
1693 
1694  pT2PDF = pT2;
1695  hasTinyPDFdau = false;
1696 
1697  // Determine overestimated z range; switch at c and b masses.
1698  if (pT2 > m2b && nQuarkIn >= 5) {
1699  nFlavour = 5;
1700  pT2minNow = max( m2b, pT2endDip);
1701  } else if (pT2 > m2c && nQuarkIn >= 4) {
1702  nFlavour = 4;
1703  pT2minNow = max( m2c, pT2endDip);
1704  } else {
1705  nFlavour = 3;
1706  pT2minNow = pT2endDip;
1707  }
1708 
1709  // Compute upper z limit
1710  zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
1711  ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
1712 
1713  // Parton density of daughter at current scale.
1714  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1715  xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2);
1716  if (xPDFdaughter < TINYPDF) {
1717  xPDFdaughter = TINYPDF;
1718  hasTinyPDFdau = true;
1719  }
1720 
1721  // Integral over f -> gamma f splitting kernel.
1722  // Normalized so: 4/3 aS/2pi P(z) -> eq^2 * aEM/2pi P(z).
1723  // (Charge-weighting happens below.)
1724  double q2gInt = 4. * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
1725  // Optionally enhanced branching rate.
1726  if (canEnhanceET) q2gInt *= userHooksPtr->enhanceFactor("isr:Q2QA");
1727 
1728 
1729  // Charge-weighted Parton density of potential quark mothers.
1730  xPDFmotherSum = 0.;
1731  for (int i = -nFlavour; i <= nFlavour; ++i) {
1732  if (i == 0) {
1733  xPDFmother[10] = 0.;
1734  } else {
1735  xPDFmother[i+10] = pow2((abs(i+1) % 2 + 1)/3.0)
1736  * beam.xfISR(iSysNow, i, xDaughter, pdfScale2);
1737  xPDFmotherSum += xPDFmother[i+10];
1738  }
1739  }
1740 
1741  // Total QED evolution coefficient for a photon.
1742  kernelPDF = q2gInt * xPDFmotherSum / xPDFdaughter;
1743 
1744  // End evaluation of splitting kernels and parton densities.
1745  needNewPDF = false;
1746  }
1747  if (kernelPDF < TINYKERNELPDF) return;
1748 
1749  // Select pT2 for next trial branching
1750  pT2 *= pow( rndmPtr->flat(), 1. / (alphaEM2pi * kernelPDF));
1751 
1752  // If crossed b threshold, continue evolution from this threshold.
1753  if (nFlavour == 5 && pT2 < m2b) {
1754  needNewPDF = true;
1755  pT2 = m2b;
1756  continue;
1757  }
1758 
1759  // If crossed c threshold, continue evolution from this threshold.
1760  else if (nFlavour == 4 && pT2 < m2c) {
1761  needNewPDF = true;
1762  pT2 = m2c;
1763  continue;
1764  }
1765 
1766  // Abort evolution if below cutoff scale, or below another branching.
1767  else if (pT2 < pT2endDip) return;
1768 
1769  // Select flavour for trial branching
1770  double temp = xPDFmotherSum * rndmPtr->flat();
1771  idMother = -nQuarkIn - 1;
1772  do {
1773  temp -= xPDFmother[(++idMother) + 10];
1774  } while (temp > 0. && idMother < nQuarkIn);
1775 
1776  // Sister is same as mother, but can have m2 > 0
1777  idSister = idMother;
1778  mSister = particleDataPtr->m0(idSister);
1779  m2Sister = pow2(mSister);
1780 
1781  // Select z for trial branching
1782  z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
1783  * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
1784 
1785  // Trial weight: splitting kernel
1786  wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z);
1787 
1788  // Trial weight: running alpha_EM
1789  double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
1790  wt *= (alphaEMnow / alphaEMmax);
1791 
1792  // Optionally enhanced branching rate.
1793  nameNow = "isr:Q2AQ";
1794  if (canEnhanceET) {
1795  double enhance = userHooksPtr->enhanceFactor(nameNow);
1796  if (enhance != 1.) {
1797  enhanceNow = enhance;
1798  isEnhancedQ2AQ = true;
1799  }
1800  }
1801 
1802  // Derive Q2 and x of mother from pT2 and z
1803  Q2 = pT2 / (1. - z);
1804  xMother = xDaughter / z;
1805  // Correction to x for massive recoiler from rescattering.
1806  if (!dipEndNow->normalRecoil) {
1807  if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1808  else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1809  }
1810 
1811  // Compute pT2corr
1812  pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1813  if (pT2corr < TINYPT2) { wt = 0.; continue; }
1814 
1815  // If creating heavy quark by Q -> gamma + Q then next g -> Q + Qbar.
1816  // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
1817  if ( abs(idMother) == 4 || abs(idMother) == 5 ) {
1818  double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
1819  double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
1820  if (pT2QQcorr < TINYPT2) { wt = 0.; continue; }
1821  }
1822 
1823  // Optional dampening of large pT values in first radiation.
1824  if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
1825  wt *= pT2damp / (pT2 + pT2damp);
1826 
1827  // Evaluation of new daughter PDF
1828  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
1829  double xPDFdaughterNew = beam.xfISR(iSysNow, idDaughter, xDaughter,
1830  pdfScale2);
1831  if (xPDFdaughterNew < TINYPDF) {
1832  xPDFdaughterNew = TINYPDF;
1833  }
1834 
1835  // Evaluation of new charge-weighted mother PDF
1836  double xPDFmotherNew = pow2( (abs(idMother+1) % 2 + 1)/3.0 )
1837  * beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
1838 
1839  // Trial weight: divide out old pdf ratio
1840  wt *= xPDFdaughter / xPDFmother[idMother + 10];
1841 
1842  // Trial weight: new pdf ratio
1843  wt *= xPDFmotherNew / xPDFdaughterNew;
1844 
1845  // Iterate until acceptable pT (or have fallen below pTmin).
1846  } while (wt < rndmPtr->flat()) ;
1847  }
1848 
1849  // Store outcome of enhanced branching rate analysis.
1850  splittingNameNow = nameNow;
1851  if (canEnhanceET) {
1852  if (isEnhancedQ2QA) storeEnhanceFactor(pT2,"isr:Q2QA", enhanceNow);
1853  if (isEnhancedQ2AQ) storeEnhanceFactor(pT2,"isr:Q2AQ", enhanceNow);
1854  if (isEnhancedA2QQ) storeEnhanceFactor(pT2,"isr:A2QQ", enhanceNow);
1855  }
1856 
1857  // Save values for (so far) acceptable branching.
1858  dipEndNow->store( idDaughter, idMother, idSister, x1Now, x2Now, m2Dip,
1859  pT2, z, xMother, Q2, mSister, m2Sister, pT2corr, 0, m2ColPair,
1860  mColPartner);
1861 
1862 }
1863 
1864 //--------------------------------------------------------------------------
1865 
1866 void SpaceShower::pT2nextWeak( double pT2begDip, double pT2endDip) {
1867 
1868  // Type of dipole and starting values.
1869  BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
1870  bool isLeptonBeam = beam.isLepton();
1871  bool isValence = beam[iSysNow].isValence();
1872  int MEtype = dipEndNow->MEtype;
1873  double pT2 = pT2begDip;
1874  double m2Lepton = (isLeptonBeam) ? pow2(beam.m()) : 0.;
1875  if (isLeptonBeam && pT2begDip < m2Lepton) return;
1876 
1877  // alpha_em at maximum scale provides upper estimate.
1878  double alphaEMmax = alphaEM.alphaEM(pT2begDip);
1879  double alphaEM2pi = alphaEMmax / (2. * M_PI);
1880 
1881  // Maximum x of mother implies minimum z = xDaughter / xMother.
1882  double xMaxAbs = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
1883  double zMinAbs = xDaughter / xMaxAbs;
1884  if (xMaxAbs < 0.) {
1885  infoPtr->errorMsg("Warning in SpaceShower::pT2nextWeak: "
1886  "xMaxAbs negative");
1887  return;
1888  }
1889 
1890  // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
1891  double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
1892  ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
1893  if (isLeptonBeam) {
1894  double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
1895  if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
1896  }
1897  if (zMaxAbs < zMinAbs) return;
1898 
1899  // Variables used inside evolution loop. (Mainly dummy start values.)
1900  int idMother = 0;
1901  int idSister = 0;
1902  // Check whether emission of W+, W- or Z0.
1903  if (dipEndNow->weakType == 1) {
1904  idSister = (idDaughter > 0) ? -24 : 24;
1905  if (abs(idDaughter) % 2 == 1) idSister = -idSister;
1906  } else if (dipEndNow->weakType == 2) idSister = 23;
1907  double z = 0.;
1908  double xMother = 0.;
1909  double wt = 0.;
1910  double Q2 = 0.;
1911  double mSister = particleDataPtr->mSel(idSister);
1912  double m2Sister = pow2(mSister);
1913  double pT2corr = 0.;
1914 
1915  // Find maximum z due to massive sister.
1916  // Still need to prove that this actually is an overestimate.
1917  double mRatio = mSister/ sqrt(m2Dip);
1918  double m2R1 = 1. + pow2(mRatio);
1919  double zMaxMassive = 1. / (m2R1 + pT2endDip/m2Dip);
1920  if (zMaxMassive < zMaxAbs) zMaxAbs = zMaxMassive;
1921  if (zMaxAbs < zMinAbs) return;
1922 
1923  // Set default values for enhanced emissions.
1924  bool isEnhancedQ2QW;
1925  isEnhancedQ2QW = false;
1926  double enhanceNow = 1.;
1927  string nameNow = "";
1928 
1929  // Weak evolution of fermions.
1930  // Integrals of splitting kernels for fermions: f -> f.
1931  // Old ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
1932  // New Ansatz f(z) = (1 + (1+r^2)^2) / (1 - z * (1 + r^2)).
1933  // This should always be an overestimate for massive emissions.
1934  // Not yet implemented correctly for lepton beam.
1935  double f2fInt = 0.;
1936  double f2fIntA = (1. + pow2(zMaxAbs * m2R1)) / m2R1
1937  * log( (1. - m2R1 * zMinAbs) / (1. - m2R1 * zMaxAbs) );
1938  double f2fIntB = 0.;
1939  double zRootMin = (1. + sqrt(m2R1 * zMinAbs)) / (1. - sqrt(m2R1 * zMinAbs));
1940  double zRootMax = (1. + sqrt(m2R1 * zMaxAbs)) / (1. - sqrt(m2R1 * zMaxAbs));
1941  if (isLeptonBeam) {
1942  f2fIntB = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
1943  f2fInt = f2fIntA + f2fIntB;
1944  } else if (isValence)
1945  f2fInt = (1. + pow2(zMaxAbs) * pow2(m2R1))/ sqrt(m2R1)
1946  * log( zRootMax / zRootMin );
1947  else
1948  f2fInt = f2fIntA;
1949 
1950  // Calculate the weak coupling: separate for left- and right-handed fermions.
1951  double weakCoupling = 0;
1952  if (dipEndNow->weakType == 1)
1953  weakCoupling = 2. * alphaEM2pi / (4. * coupSMPtr->sin2thetaW());
1954  else if (dipEndNow->weakType == 2 && dipEndNow->weakPol == -1)
1955  weakCoupling = alphaEM2pi * thetaWRat
1956  * pow2(2. * coupSMPtr->lf( abs(idDaughter) ));
1957  else
1958  weakCoupling = alphaEM2pi * thetaWRat
1959  * pow2(2. * coupSMPtr->rf(abs(idDaughter) ));
1960 
1961  // Find the possible mothers for a W emission. So far quarks only.
1962  vector<int> possibleMothers;
1963  if (abs(idDaughter) % 2 == 0) {
1964  possibleMothers.push_back(1);
1965  possibleMothers.push_back(3);
1966  possibleMothers.push_back(5);
1967  } else {
1968  possibleMothers.push_back(2);
1969  possibleMothers.push_back(4);
1970  }
1971  if (idDaughter < 0)
1972  for (unsigned int i = 0;i < possibleMothers.size();i++)
1973  possibleMothers[i] = - possibleMothers[i];
1974 
1975  // Check if daughter estimate is 0, return in that case.
1976  // Only write warning if u, d or g is the daughter.
1977  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2begDip;
1978  double xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2);
1979  if (xPDFdaughter < TINYPDF) {
1980  if (abs(idDaughter) == 1 || abs(idDaughter) == 2 || abs(idDaughter) == 21)
1981  infoPtr->errorMsg("Warning in SpaceShower::pT2nextWeak: "
1982  "very small PDF");
1983  return;
1984  }
1985 
1986  // PDF and CKM upper estimate needed for W emission.
1987  double overEstimatePDFCKM = 0.;
1988  if (dipEndNow->weakType == 1) {
1989  for (unsigned int i = 0; i < possibleMothers.size(); i++)
1990  overEstimatePDFCKM += coupSMPtr->V2CKMid(idDaughter, possibleMothers[i])
1991  * beam.xfISR(iSysNow, possibleMothers[i], xDaughter, pdfScale2)
1992  / xPDFdaughter;
1993  }
1994  if (dipEndNow->weakType == 2) overEstimatePDFCKM = 1.;
1995 
1996  // Upper estimate for evolution equation, including fudge factor.
1997  if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
1998  double kernelPDF = weakCoupling * f2fInt * weakEnhancement;
1999 
2000  // PDF and CKM overestimate.
2001  kernelPDF *= overEstimatePDFCKM;
2002  double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
2003  kernelPDF *= fudge;
2004  if (kernelPDF < TINYKERNELPDF) return;
2005  // Optionally enhanced branching rate.
2006  if (canEnhanceET) kernelPDF *= userHooksPtr->enhanceFactor("isr:Q2QW");
2007 
2008  // Begin evolution loop towards smaller pT values.
2009  do {
2010 
2011  // Default values for current tentative emission.
2012  isEnhancedQ2QW = false;
2013  enhanceNow = 1.;
2014  nameNow = "";
2015 
2016  // Pick pT2 (in overestimated z range).
2017  // For l -> l gamma include extrafactor 1 / ln(pT2 / m2l) in evolution.
2018  double shift = pow(rndmPtr->flat(), 1. / kernelPDF);
2019  if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
2020  else pT2 = pT2 * shift;
2021 
2022  // Abort evolution if below cutoff scale, or below another branching.
2023  if (pT2 < pT2endDip) return;
2024  if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
2025 
2026  // Abort evolution if below mass treshold.
2027  if (pT2 < HEAVYPT2EVOL * pow2(particleDataPtr->m0(idDaughter))) return;
2028 
2029  // Set the id for the mother particle to be equal to the daughter
2030  // particle. This is correct for Z, and it will later be changed for W.
2031  idMother = idDaughter;
2032 
2033  // Select z value of branching f -> f + Z/W, and corrective weight.
2034  wt = 1.0;
2035  if (isLeptonBeam) {
2036  if (f2fIntA > rndmPtr->flat() * (f2fIntA + f2fIntB)) {
2037  z = 1. - (1. - zMinAbs)
2038  * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
2039  } else {
2040  z = xDaughter + (zMinAbs - xDaughter)
2041  * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter),
2042  rndmPtr->flat() );
2043  }
2044  wt *= (z - xDaughter) / (1. - xDaughter);
2045  } else if (isValence) {
2046  // Valence more peaked at large z.
2047  double zTmp = zRootMin * pow(zRootMax / zRootMin, rndmPtr->flat() );
2048  z = pow2( (1. - zTmp) / (1. + zTmp) ) / m2R1;
2049  wt *= sqrt(z);
2050  } else {
2051  z = (1. - (1. - zMinAbs * m2R1) * pow( (1. - zMaxAbs * m2R1)
2052  / (1. - zMinAbs * m2R1), rndmPtr->flat() ) ) / m2R1;
2053  }
2054  wt *= (1. + pow2(z * m2R1)) / (1. + pow2(zMaxAbs * m2R1));
2055 
2056  // Optionally enhanced branching rate.
2057  nameNow = "isr:Q2QW";
2058  if (canEnhanceET) {
2059  double enhance = userHooksPtr->enhanceFactor(nameNow);
2060  if (enhance != 1.) {
2061  enhanceNow = enhance;
2062  isEnhancedQ2QW = true;
2063  }
2064  }
2065 
2066  // Derive Q2 and x of mother from pT2 and z.
2067  Q2 = pT2 / (1. - z);
2068  xMother = xDaughter / z;
2069 
2070  // Correction to x for massive recoiler from rescattering.
2071  if (!dipEndNow->normalRecoil) {
2072  if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
2073  else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
2074  }
2075  if (xMother > xMaxAbs) { wt = 0.; continue; }
2076 
2077  // Forbidden emission if outside allowed z range for given pT2.
2078  pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
2079  if (pT2corr < TINYPT2) { wt = 0.; continue; }
2080 
2081  // Correct by ln(pT2 / m2l) and fudge factor.
2082  if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
2083 
2084  // Evaluation of ME correction.
2085  if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
2086  m2Dip, z, Q2, m2Sister) / calcMEmax(MEtype, idMother, idDaughter);
2087 
2088  // Optional dampening of large pT values in first radiation.
2089  // Allow damping also for corrected matrix elements
2090  if (dopTdamp && iSysNow == 0 && nRad == 0)
2091  wt *= pT2damp / (pT2 + pT2damp);
2092 
2093  // Correct to current value of alpha_EM.
2094  double alphaEMnow = alphaEM.alphaEM(pT2);
2095  wt *= (alphaEMnow / alphaEMmax);
2096 
2097  // Evaluation of new daughter and mother PDF's for Z.
2098  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
2099  double xPDFdaughterNew = max ( TINYPDF,
2100  beam.xfISR(iSysNow, idDaughter, xDaughter, pdfScale2) );
2101  if (dipEndNow->weakType == 2) {
2102  double xPDFmotherNew
2103  = beam.xfISR(iSysNow, idMother, xMother, pdfScale2);
2104  wt *= xPDFmotherNew / xPDFdaughterNew;
2105  }
2106 
2107  // Evaluation of daughter and mother PDF's for W.
2108  if (dipEndNow->weakType == 1) {
2109  double valPDFCKM[3] = {0.};
2110  double sumPDFCKM = 0.;
2111  for (unsigned int i = 0; i < possibleMothers.size(); i++) {
2112  valPDFCKM[i] = coupSMPtr->V2CKMid(idDaughter,possibleMothers[i])
2113  * beam.xfISR(iSysNow, possibleMothers[i], xMother, pdfScale2)
2114  / xPDFdaughterNew;
2115  sumPDFCKM += valPDFCKM[i];
2116  }
2117  wt *= sumPDFCKM / overEstimatePDFCKM;
2118 
2119  // Choose id for mother in case of W.
2120  double pickId = sumPDFCKM * rndmPtr->flat();
2121  int iId = -1;
2122  int pmSize = possibleMothers.size();
2123  do pickId -= valPDFCKM[++iId];
2124  while (pickId > 0. && iId < pmSize);
2125  idMother = possibleMothers[iId];
2126  }
2127 
2128  // Warn if too big weight.
2129  if (wt > 1.) infoPtr->errorMsg("Warning in SpaceShower::pT2nextWeak: "
2130  "weight is above unity.");
2131 
2132  // Iterate until acceptable pT (or have fallen below pTmin).
2133  } while (wt < rndmPtr->flat()) ;
2134 
2135  // Store outcome of enhanced branching rate analysis.
2136  splittingNameNow = nameNow;
2137  if (canEnhanceET && isEnhancedQ2QW)
2138  storeEnhanceFactor(pT2,"isr:Q2QW", enhanceNow);
2139 
2140  // Save values for (so far) acceptable branching.
2141  dipEndNow->store( idDaughter, idMother, idSister, x1Now, x2Now, m2Dip,
2142  pT2, z, xMother, Q2, mSister, m2Sister, pT2corr, 0, m2ColPair,
2143  mColPartner);
2144 
2145 }
2146 
2147 //--------------------------------------------------------------------------
2148 
2149 // Kinematics of branching.
2150 // Construct mother -> daughter + sister, with recoiler on other side.
2151 
2152 bool SpaceShower::branch( Event& event) {
2153 
2154  // Side on which branching occured.
2155  int side = abs(dipEndSel->side);
2156  double sideSign = (side == 1) ? 1. : -1.;
2157 
2158  // Read in flavour and colour variables.
2159  int iDaughter = partonSystemsPtr->getInA(iSysSel);
2160  int iRecoiler = partonSystemsPtr->getInB(iSysSel);
2161  if (side == 2) swap(iDaughter, iRecoiler);
2162  int idDaughterNow = dipEndSel->idDaughter;
2163  int idMother = dipEndSel->idMother;
2164  int idSister = dipEndSel->idSister;
2165  int idRecoiler = event[iRecoiler].id();
2166  int colDaughter = event[iDaughter].col();
2167  int acolDaughter = event[iDaughter].acol();
2168  int iColPartner = dipEndSel->iColPartner;
2169 
2170  // Recoil parton may be rescatterer, requiring special processing.
2171  bool normalRecoil = dipEndSel->normalRecoil;
2172  int iRecoilMother = event[iRecoiler].mother1();
2173 
2174  // Read in kinematical variables.
2175  double x1 = dipEndSel->x1;
2176  double x2 = dipEndSel->x2;
2177  double xMo = dipEndSel->xMo;
2178  double m2II = dipEndSel->m2Dip;
2179  double mII = sqrt(m2II);
2180  double pT2 = dipEndSel->pT2;
2181  double z = dipEndSel->z;
2182  double Q2 = dipEndSel->Q2;
2183  double mSister = dipEndSel->mSister;
2184  double m2Sister = dipEndSel->m2Sister;
2185  double pT2corr = dipEndSel->pT2corr;
2186  double x1New = (side == 1) ? xMo : x1;
2187  double x2New = (side == 2) ? xMo : x2;
2188 
2189  // Flag for gamma -> q qbar splittings.
2190  gamma2qqbar = false;
2191 
2192  // Current beam particle.
2193  BeamParticle& beamNow = (side == 1) ? *beamAPtr : *beamBPtr;
2194 
2195  // If gamma->qqbar splitting and nMPI > 1 then only save pT2 value
2196  // and constrcut the kinematics in beamRemnants.
2197  if ( beamNow.isGamma() && idMother == 22 && infoPtr->nMPI() > 1 ) {
2198  gamma2qqbar = true;
2199  beamNow.resolvedGamma(false);
2200  beamNow.pT2gamma2qqbar(pT2corr);
2201  beamNow.gamVal(iSysSel);
2202  return true;
2203  }
2204 
2205  // Read in MEtype. Four-vectors to reconstruct.
2206  int MEtype = dipEndSel->MEtype;
2207  Vec4 pMother, pSister, pNewRec, pNewColPartner;
2208 
2209  // Rescatter: kinematics may fail; use the rescatterFail flag to tell
2210  // parton level to try again.
2211  rescatterFail = false;
2212 
2213  // Kinematics for II dipole.
2214  // Construct kinematics of mother, sister and recoiler in old rest frame.
2215  // Normally both mother and recoiler are taken massless.
2216  if (iColPartner == 0) {
2217  double eNewRec, pzNewRec, pTbranch, pzMother;
2218  if (normalRecoil) {
2219  eNewRec = 0.5 * (m2II + Q2) / mII;
2220  pzNewRec = -sideSign * eNewRec;
2221  pTbranch = sqrt(pT2corr) * m2II / ( z * (m2II + Q2) );
2222  pzMother = sideSign * 0.5 * mII * ( (m2II - Q2)
2223  / ( z * (m2II + Q2) ) + (Q2 + m2Sister) / m2II );
2224  // More complicated kinematics when recoiler not massless. May fail.
2225  } else {
2226  m2Rec = event[iRecoiler].m2();
2227  double s1Tmp = m2II + Q2 - m2Rec;
2228  double s3Tmp = m2II / z - m2Rec;
2229  double r1Tmp = sqrt(s1Tmp * s1Tmp + 4. * Q2 * m2Rec);
2230  eNewRec = 0.5 * (m2II + m2Rec + Q2) / mII;
2231  pzNewRec = -sideSign * 0.5 * r1Tmp / mII;
2232  double pT2br = Q2 * s3Tmp * (m2II / z - m2II - Q2)
2233  - m2Sister * s1Tmp * s3Tmp - m2Rec * pow2(Q2 + m2Sister);
2234  if (pT2br <= 0.) return false;
2235  pTbranch = sqrt(pT2br) / r1Tmp;
2236  pzMother = sideSign * (mII * s3Tmp
2237  - eNewRec * (m2II / z - Q2 - m2Rec - m2Sister)) / r1Tmp;
2238  }
2239  // Common final kinematics steps for both normal and rescattering.
2240  double eMother = sqrt( pow2(pTbranch) + pow2(pzMother) );
2241  double pzSister = pzMother + pzNewRec;
2242  double eSister = sqrt( pow2(pTbranch) + pow2(pzSister) + m2Sister );
2243  pMother.p( pTbranch, 0., pzMother, eMother );
2244  pSister.p( pTbranch, 0., pzSister, eSister );
2245  pNewRec.p( 0., 0., pzNewRec, eNewRec );
2246 
2247  // Kinematics of IF dipole.
2248  // Construct kinematics in old rest frame of daughter + colour partner.
2249  // Mother and recoiler massless but massive colour partner.
2250  } else {
2251  double m2IF = dipEndSel->m2IF;
2252  double mIF = sqrt(m2IF);
2253  m2ColPartner = pow2( dipEndSel->mColPartner );
2254 
2255  // Construct kinematics.
2256  double m2IFred = m2IF - m2ColPartner;
2257  pMother.p( 0., 0., 0.5 * m2IFred / (z * mIF), 0.5 * m2IFred / (z * mIF) );
2258  Vec4 pShift( 0., 0.,
2259  -0.5 * Q2 / mIF - z * (Q2 + m2Sister) * m2ColPartner / (mIF * m2IFred),
2260  -0.5 * Q2 / mIF + z * (Q2 + m2Sister) / mIF );
2261  pSister = (1. - z) * pMother + pShift;
2262  pNewColPartner = Vec4( 0., 0., -0.5 * m2IFred /mIF,
2263  0.5 * (m2IF + m2ColPartner) / mIF ) - pShift;
2264  // Do not change the momentum of the recoiler (in event frame).
2265  pNewRec.p( event[iRecoiler].p() );
2266 
2267  // Flat azimuthal angle.
2268  double phi = 2. * M_PI * rndmPtr->flat();
2269 
2270  // Evaluate coefficient of azimuthal asymmetry from gluon polarization.
2271  findAsymPol( event, dipEndSel);
2272  int iFinPol = dipEndSel->iFinPol;
2273  double cPol = dipEndSel->asymPol;
2274  double phiPol = (iFinPol == 0) ? 0. : event[iFinPol].phi();
2275 
2276  // Bias phi distribution for polarization.
2277  if (iFinPol != 0) {
2278  double cPhiPol, weight;
2279  // Boost back to the rest frame of daughter + recoiler.
2280  RotBstMatrix M;
2281  M.fromCMframe( event[iDaughter].p(), event[iColPartner].p() );
2282  double pTcorr = sqrt(pT2corr);
2283  do {
2284  phi = 2. * M_PI * rndmPtr->flat();
2285  weight = 1.;
2286  Vec4 pSisTmp = pSister + pTcorr * Vec4( cos(phi), sin(phi), 0., 0.);
2287  pSisTmp.rotbst(M);
2288  cPhiPol = cos(pSisTmp.phi() - phiPol);
2289  weight *= ( 1. + cPol * (2. * pow2(cPhiPol) - 1.) )
2290  / ( 1. + abs(cPol) );
2291  } while (weight < rndmPtr->flat());
2292  }
2293 
2294  // Add azimuthal part to the kinematics.
2295  pSister.px( sqrt(pT2corr) * cos(phi) );
2296  pSister.py( sqrt(pT2corr) * sin(phi) );
2297  pNewColPartner.px( - pSister.px() );
2298  pNewColPartner.py( - pSister.py() );
2299  }
2300 
2301  // Current event and subsystem size.
2302  int eventSizeOld = event.size();
2303  int systemSizeOld = partonSystemsPtr->sizeAll(iSysSel);
2304 
2305  // Save properties to be restored in case of user-hook veto of emission.
2306  int beamOff1 = 1 + beamOffset;
2307  int beamOff2 = 2 + beamOffset;
2308  int ev1Dau1V = event[beamOff1].daughter1();
2309  int ev2Dau1V = event[beamOff2].daughter1();
2310  vector<int> statusV, mother1V, mother2V, daughter1V, daughter2V;
2311 
2312  // Check if the first emission should be checked for removal.
2313  bool canMergeFirst = (mergingHooksPtr != 0)
2314  ? mergingHooksPtr->canVetoEmission() : false;
2315 
2316  // Check if doing uncertainty variations
2317  doUncertaintiesNow = doUncertainties;
2318  if (!uVarMPIshowers && iSysSel != 0) doUncertaintiesNow = false;
2319  if (pT2 < uVarpTmin2) doUncertaintiesNow = false;
2320 
2321  // Save further properties to be restored.
2322  if (canVetoEmission || canMergeFirst || canEnhanceET || doWeakShower
2323  || doUncertainties) {
2324  for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
2325  int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
2326  statusV.push_back( event[iOldCopy].status());
2327  mother1V.push_back( event[iOldCopy].mother1());
2328  mother2V.push_back( event[iOldCopy].mother2());
2329  daughter1V.push_back( event[iOldCopy].daughter1());
2330  daughter2V.push_back( event[iOldCopy].daughter2());
2331  }
2332  }
2333 
2334  // Take copy of existing system, to be given modified kinematics.
2335  // Incoming negative status. Rescattered also negative, but after copy.
2336  int iNewColPartner = 0;
2337  for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
2338  int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
2339  int statusOld = event[iOldCopy].status();
2340  int statusNew = (iOldCopy == iDaughter
2341  || iOldCopy == iRecoiler) ? statusOld : 44;
2342  int iNewCopy = event.copy(iOldCopy, statusNew);
2343  if (statusOld < 0) event[iNewCopy].statusNeg();
2344  if (iOldCopy == iColPartner) iNewColPartner = iNewCopy;
2345  }
2346 
2347  // Define colour flow in branching.
2348  // Default corresponds to f -> f + gamma.
2349  int colMother = colDaughter;
2350  int acolMother = acolDaughter;
2351  int colSister = 0;
2352  int acolSister = 0;
2353  if (idSister == 22 || idSister == 23 || abs(idSister) == 24) ;
2354  // q -> q + g and 50% of g -> g + g; need new colour.
2355  else if (idSister == 21 && ( (idMother > 0 && idMother < 9)
2356  || (idMother == 21 && rndmPtr->flat() < 0.5) ) ) {
2357  colMother = event.nextColTag();
2358  colSister = colMother;
2359  acolSister = colDaughter;
2360  // qbar -> qbar + g and other 50% of g -> g + g; need new colour.
2361  } else if (idSister == 21) {
2362  acolMother = event.nextColTag();
2363  acolSister = acolMother;
2364  colSister = acolDaughter;
2365  // q -> g + q.
2366  } else if (idDaughterNow == 21 && idMother > 0) {
2367  colMother = colDaughter;
2368  acolMother = 0;
2369  colSister = acolDaughter;
2370  // qbar -> g + qbar
2371  } else if (idDaughterNow == 21) {
2372  acolMother = acolDaughter;
2373  colMother = 0;
2374  acolSister = colDaughter;
2375  // g -> q + qbar but not gamma -> q + qbar.
2376  } else if ( (idDaughterNow > 0 && idDaughterNow < 9) && idMother != 22) {
2377  acolMother = event.nextColTag();
2378  acolSister = acolMother;
2379  // g -> qbar + q but not gamma -> qbar + q.
2380  } else if ( (idDaughterNow < 0 && idDaughterNow > -9) && idMother != 22) {
2381  colMother = event.nextColTag();
2382  colSister = colMother;
2383  // q -> gamma + q.
2384  } else if (idDaughterNow == 22 && idMother > 0) {
2385  colMother = event.nextColTag();
2386  colSister = colMother;
2387  // qbar -> gamma + qbar.
2388  } else if (idDaughterNow == 22) {
2389  acolMother = event.nextColTag();
2390  acolSister = acolMother;
2391  // gamma -> q + qbar.
2392  } else if ( (idDaughterNow > 0 && idDaughterNow < 9) && idMother == 22) {
2393  colMother = 0;
2394  acolMother = 0;
2395  acolSister = colDaughter;
2396  gamma2qqbar = true;
2397  // gamma -> qbar + q.
2398  } else if ( (idDaughterNow < 0 && idDaughterNow > -9) && idMother == 22) {
2399  colMother = 0;
2400  acolMother = 0;
2401  colSister = acolDaughter;
2402  gamma2qqbar = true;
2403  }
2404 
2405  // Indices of partons involved. Add new sister.
2406  int iMother = eventSizeOld + side - 1;
2407  int iNewRec = eventSizeOld + 2 - side;
2408  int iSister = event.append( idSister, 43, iMother, 0, 0, 0,
2409  colSister, acolSister, pSister, mSister, sqrt(pT2) );
2410 
2411  // References to the partons involved.
2412  Particle& daughter = event[iDaughter];
2413  Particle& mother = event[iMother];
2414  Particle& newRecoiler = event[iNewRec];
2415  Particle& sister = event.back();
2416 
2417  // Allow setting of vertex for daughter parton, recoiler and sister.
2418  if (doPartonVertex) {
2419  if (!daughter.hasVertex()) partonVertexPtr->vertexISR( iDaughter, event);
2420  if (!newRecoiler.hasVertex()) partonVertexPtr->vertexISR( iNewRec, event);
2421  if (!sister.hasVertex()) partonVertexPtr->vertexISR( iSister, event);
2422  }
2423 
2424  // Replace old by new mother; update new recoiler.
2425  mother.id( idMother );
2426  mother.status( -41);
2427  mother.cols( colMother, acolMother);
2428  mother.p( pMother);
2429  if (mother.idAbs() == 21 || mother.idAbs() == 22) mother.pol(9);
2430  newRecoiler.status( (normalRecoil) ? -42 : -46 );
2431  newRecoiler.p( pNewRec);
2432  if (!normalRecoil) newRecoiler.m( event[iRecoiler].m() );
2433  // Update the colour partner in case of dipole recoil.
2434  if (iNewColPartner != 0) event[iNewColPartner].p( pNewColPartner );
2435 
2436  // Update mother and daughter pointers; also for beams.
2437  daughter.mothers( iMother, 0);
2438  mother.daughters( iSister, iDaughter);
2439  if (iSysSel == 0) {
2440  event[beamOff1].daughter1( (side == 1) ? iMother : iNewRec );
2441  event[beamOff2].daughter1( (side == 2) ? iMother : iNewRec );
2442  }
2443 
2444  // Special checks to set weak particles status equal to 47.
2445  if (sister.idAbs() == 23 || sister.idAbs() == 24) sister.status(47);
2446 
2447  // Normal azimuth and boost procedure for II dipole.
2448  if (iNewColPartner == 0) {
2449 
2450  // Find boost to old rest frame.
2451  RotBstMatrix Mtot;
2452  if (normalRecoil) Mtot.bst(0., 0., (x2 - x1) / (x1 + x2) );
2453  else if (side == 1)
2454  Mtot.toCMframe( event[iDaughter].p(), event[iRecoiler].p() );
2455  else Mtot.toCMframe( event[iRecoiler].p(), event[iDaughter].p() );
2456 
2457  // Initially select phi angle of branching at random.
2458  double phi = 2. * M_PI * rndmPtr->flat();
2459 
2460  // Evaluate coefficient of azimuthal asymmetry from gluon polarization.
2461  findAsymPol( event, dipEndSel);
2462  int iFinPol = dipEndSel->iFinPol;
2463  double cPol = dipEndSel->asymPol;
2464  double phiPol = (iFinPol == 0) ? 0. : event[iFinPol].phi();
2465 
2466  // If interference: try to match sister (anti)colour to final state.
2467  int iFinInt = 0;
2468  double cInt = 0.;
2469  double phiInt = 0.;
2470  if (doPhiIntAsym) {
2471  for (int i = 0; i < partonSystemsPtr->sizeOut(iSysSel); ++ i) {
2472  int iOut = partonSystemsPtr->getOut(iSysSel, i);
2473  if ( (acolSister != 0 && event[iOut].col() == acolSister)
2474  || (colSister != 0 && event[iOut].acol() == colSister) )
2475  iFinInt = iOut;
2476  }
2477  if (iFinInt != 0) {
2478  // Boost final-state parton to current frame of new kinematics.
2479  Vec4 pFinTmp = event[iFinInt].p();
2480  pFinTmp.rotbst(Mtot);
2481  double theFin = pFinTmp.theta();
2482  if (side == 2) theFin = M_PI - theFin;
2483  double theSis = pSister.theta();
2484  if (side == 2) theSis = M_PI - theSis;
2485  cInt = strengthIntAsym * 2. * theSis * theFin
2486  / (pow2(theSis) + pow2(theFin));
2487  phiInt = event[iFinInt].phi();
2488  }
2489  }
2490 
2491  // Bias phi distribution for polarization and interference.
2492  if (iFinPol != 0 || iFinInt != 0) {
2493  double cPhiPol, cPhiInt, weight;
2494  do {
2495  phi = 2. * M_PI * rndmPtr->flat();
2496  weight = 1.;
2497  if (iFinPol !=0 ) {
2498  cPhiPol = cos(phi - phiPol);
2499  weight *= ( 1. + cPol * (2. * pow2(cPhiPol) - 1.) )
2500  / ( 1. + abs(cPol) );
2501  }
2502  if (iFinInt !=0 ) {
2503  cPhiInt = cos(phi - phiInt);
2504  weight *= (1. - cInt) * (1. - cInt * cPhiInt)
2505  / (1. + pow2(cInt) - 2. * cInt * cPhiInt);
2506  }
2507  } while (weight < rndmPtr->flat());
2508  }
2509 
2510  // Include rotation -phi on boost to old rest frame.
2511  Mtot.rot(0., -phi);
2512 
2513  // Find boost from old rest frame to event cm frame.
2514  RotBstMatrix MfromRest;
2515  // The boost to the new rest frame.
2516  Vec4 sumNew = pMother + pNewRec;
2517  double betaX = sumNew.px() / sumNew.e();
2518  double betaZ = sumNew.pz() / sumNew.e();
2519  MfromRest.bst( -betaX, 0., -betaZ);
2520  // Alignment of radiator + recoiler to +- z axis, and rotation +phi.
2521  // Note: with spacelike (E < 0) recoiler p'_x_mother < 0 can happen!
2522  pMother.rotbst(MfromRest);
2523  double theta = pMother.theta();
2524  if (pMother.px() < 0.) theta = -theta;
2525  if (side == 2) theta += M_PI;
2526  MfromRest.rot(-theta, phi);
2527  // Boost to radiator + recoiler in event cm frame.
2528  if (normalRecoil) {
2529  MfromRest.bst( 0., 0., (x1New - x2New) / (x1New + x2New) );
2530  } else if (side == 1) {
2531  Vec4 pMotherWanted( 0., 0., 0.5 * eCM * x1New, 0.5 * eCM * x1New);
2532  MfromRest.fromCMframe( pMotherWanted, event[iRecoiler].p() );
2533  } else {
2534  Vec4 pMotherWanted( 0., 0., -0.5 * eCM * x2New, 0.5 * eCM * x2New);
2535  MfromRest.fromCMframe( event[iRecoiler].p(), pMotherWanted );
2536  }
2537  Mtot.rotbst(MfromRest);
2538 
2539  // ME correction for weak emissions in the t-channel.
2540  double wt;
2541  if (MEtype == 201 || MEtype == 202 || MEtype == 203 ||
2542  MEtype == 206 || MEtype == 207 || MEtype == 208) {
2543 
2544  // Start by finding the correct outgoing particles for the ME correction.
2545  Vec4 pA0 = mother.p();
2546  Vec4 pB = newRecoiler.p();
2547  bool sideRad = (abs(side) == 1);
2548  Vec4 p1,p2;
2549  if (weakExternal) {
2550  p1 = weakMomenta[2];
2551  p2 = weakMomenta[3];
2552  } else {
2553  p1 = event[5].p();
2554  p2 = event[6].p();
2555  }
2556  if (!tChannel) swap(p1,p2);
2557  if (!sideRad) swap(p1,p2);
2558 
2559  // Rotate with -phi to keep correct for the later +phi rotation.
2560  p1.rot(0., -phi);
2561  p2.rot(0., -phi);
2562 
2563  // Calculate the actual weight.
2564  wt = calcMEcorrWeak(MEtype, m2II, z, pT2, pA0, pB,
2565  event[iDaughter].p(), event[iRecoiler].p(), p1, p2, sister.p());
2566  if (wt > weakMaxWt) {
2567  weakMaxWt = wt;
2568  infoPtr->errorMsg("Warning in SpaceShower::Branch: "
2569  "weight is above unity for weak emission.");
2570  }
2571 
2572  // If weighting fails then restore event record to state before emission.
2573  if (wt < rndmPtr->flat()) {
2574  event.popBack( event.size() - eventSizeOld);
2575  event[beamOff1].daughter1( ev1Dau1V);
2576  event[beamOff2].daughter1( ev2Dau1V);
2577  for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
2578  int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
2579  event[iOldCopy].status( statusV[iCopy]);
2580  event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
2581  event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
2582  }
2583  return false;
2584  }
2585  }
2586 
2587  // Perform cumulative rotation/boost operation.
2588  // Mother, recoiler and sister from old rest frame to event cm frame.
2589  mother.rotbst(MfromRest);
2590  newRecoiler.rotbst(MfromRest);
2591  sister.rotbst(MfromRest);
2592  // The rest from (and to) event cm frame.
2593  for ( int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i)
2594  event[i].rotbst(Mtot);
2595 
2596  // Simpler special boost procedure for IF dipole.
2597  } else {
2598 
2599  // Find boost from daughter + colour partner rest frame to the event frame.
2600  RotBstMatrix MfromRest;
2601  MfromRest.fromCMframe( event[iDaughter].p(), event[iColPartner].p() );
2602  // Boost mother, sister and new colour partner (recoiler already fine).
2603  mother.rotbst(MfromRest);
2604  sister.rotbst(MfromRest);
2605  event[iNewColPartner].rotbst(MfromRest);
2606  }
2607 
2608  // Remove double counting. Only implemented for QCD hard processes
2609  // and for the first emission.
2610  if (infoPtr->nISR() + infoPtr->nFSRinProc() == 0
2611  && infoPtr->code() > 110 && infoPtr->code() < 130
2612  && MEtype >= 200 && MEtype < 210 && vetoWeakJets) {
2613 
2614  // Find outgoing particles.
2615  int iP1 = event[5].daughter1();
2616  int iP2 = event[6].daughter1();
2617  Vec4 pP1 = event[iP1].p();
2618  Vec4 pP2 = event[iP2].p();
2619 
2620  // Set start pT2 as pT2 of emitted particle and therefore no cut.
2621  double d = sister.pT2();
2622  bool cut = false;
2623 
2624  if (pP1.pT2() < d) {d = pP1.pT2(); cut = true;}
2625  if (pP2.pT2() < d) {d = pP2.pT2(); cut = true;}
2626 
2627  // Check for angle between weak boson and quarks
2628  // (require final state particle to be a fermion).
2629  if (event[iP1].idAbs() < 20) {
2630  double dij = min(pP1.pT2(),sister.pT2())
2631  * pow2(RRapPhi(pP1,sister.p()))/vetoWeakDeltaR2;
2632  if (dij < d) {
2633  d = dij;
2634  cut = false;
2635  }
2636  }
2637 
2638  if (event[iP2].idAbs() < 20) {
2639  double dij = min(pP2.pT2(),sister.pT2())
2640  * pow2(RRapPhi(pP2,sister.p()))/vetoWeakDeltaR2;
2641  if (dij < d) {
2642  d = dij;
2643  cut = false;
2644  }
2645  }
2646 
2647  // Check for angle between recoiler and radiator, if quark anti-quark pair,
2648  // or if the recoiler is a gluon.
2649  if (event[iP1].idAbs() == 21 || event[iP2].idAbs() == 21 ||
2650  event[iP1].id() == - event[iP2].id()) {
2651  double dij = min(pP1.pT2(),pP2.pT2())
2652  * pow2(RRapPhi(pP1,pP2))/vetoWeakDeltaR2;
2653  if (dij < d) {
2654  d = dij;
2655  cut = true;
2656  }
2657  }
2658 
2659  // Clean up event if the emission should be removed.
2660  if (cut) {
2661  event.popBack( event.size() - eventSizeOld);
2662  event[beamOff1].daughter1( ev1Dau1V);
2663  event[beamOff2].daughter1( ev2Dau1V);
2664  for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
2665  int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
2666  event[iOldCopy].status( statusV[iCopy]);
2667  event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
2668  event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
2669  }
2670  return false;
2671  }
2672  }
2673 
2674  // Allow veto of branching. If so restore event record to before emission.
2675  if ( (canVetoEmission
2676  && userHooksPtr->doVetoISREmission(eventSizeOld, event, iSysSel))
2677  || (canMergeFirst
2678  && mergingHooksPtr->doVetoEmission( event )) ) {
2679  event.popBack( event.size() - eventSizeOld);
2680  event[beamOff1].daughter1( ev1Dau1V);
2681  event[beamOff2].daughter1( ev2Dau1V);
2682  for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
2683  int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
2684  event[iOldCopy].status( statusV[iCopy]);
2685  event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
2686  event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
2687  }
2688  return false;
2689  }
2690 
2691  // Recover delayed shower-accept probability for uncertainty variations.
2692  // This should occur after ISR emission veto, because that is a
2693  // phase space restriction.
2694  double pAccept = dipEndSel->pAccept;
2695 
2696  // Decide if we are going to accept or reject this branching.
2697  // (Without wasting time generating random numbers if pAccept = 1.)
2698  bool acceptEvent = true;
2699  if (pAccept < 1.0) acceptEvent = (rndmPtr->flat() < pAccept);
2700 
2701  // Default values for uncertainty calculations
2702  double weight = 1.;
2703  double vp = 0.;
2704  bool vetoedEnhancedEmission = false;
2705 
2706  // Calculate event weight for enhanced emission rate.
2707  if (canEnhanceET) {
2708  // Check if emission weight was enhanced. Get enhance weight factor.
2709  bool foundEnhance = false;
2710  // Move backwards as last elements have highest pT, thus are chosen
2711  // splittings.
2712  for ( map<double,pair<string,double> >::reverse_iterator
2713  it = enhanceFactors.rbegin();
2714  it != enhanceFactors.rend(); ++it ){
2715  if (splittingNameSel.find(it->second.first) != string::npos
2716  && abs(it->second.second-1.0) > 1e-9) {
2717  foundEnhance = true;
2718  weight = it->second.second;
2719  vp = userHooksPtr->vetoProbability(splittingNameSel);
2720  break;
2721  }
2722  }
2723 
2724  // Check emission veto.
2725  if (foundEnhance && rndmPtr->flat() < vp ) vetoedEnhancedEmission = true;
2726  // Calculate new event weight.
2727  double rwgt = 1.;
2728  if (foundEnhance && vetoedEnhancedEmission) rwgt *= (1.-1./weight)/vp;
2729  else if (foundEnhance) rwgt *= 1./((1.-vp)*weight);
2730 
2731  // Reset enhance factors after usage.
2732  enhanceFactors.clear();
2733 
2734  // Set events weights, so that these could be used externally.
2735  double wtOld = userHooksPtr->getEnhancedEventWeight();
2736  if (!doTrialNow && canEnhanceEmission && !doUncertaintiesNow)
2737  userHooksPtr->setEnhancedEventWeight(wtOld*rwgt);
2738  if ( doTrialNow && canEnhanceTrial)
2739  userHooksPtr->setEnhancedTrial(sqrt(pT2), weight);
2740  // Increment counter of rejected splittings.
2741  if (vetoedEnhancedEmission && canEnhanceEmission) infoPtr->addCounter(40);
2742  }
2743 
2744  if (vetoedEnhancedEmission) acceptEvent = false;
2745 
2746  // If doing uncertainty variations, calculate accept/reject reweightings.
2747  if (doUncertaintiesNow) calcUncertainties( acceptEvent, pAccept, pT20,
2748  weight, vp, dipEndSel, &mother, &sister);
2749 
2750  // Veto if necessary.
2751  // Return false if we decided to reject this branching.
2752  if ( (doUncertainties && !acceptEvent)
2753  || (vetoedEnhancedEmission && canEnhanceEmission) ) {
2754  // Restore kinematics before returning.
2755  event.popBack( event.size() - eventSizeOld);
2756  event[beamOff1].daughter1( ev1Dau1V);
2757  event[beamOff2].daughter1( ev2Dau1V);
2758  for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
2759  int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
2760  event[iOldCopy].status( statusV[iCopy]);
2761  event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
2762  event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
2763  }
2764  return false;
2765  }
2766 
2767 
2768 
2769  // Update list of partons in system; adding newly produced one.
2770  partonSystemsPtr->setInA(iSysSel, eventSizeOld);
2771  partonSystemsPtr->setInB(iSysSel, eventSizeOld + 1);
2772  for (int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
2773  partonSystemsPtr->setOut(iSysSel, iCopy - 2, eventSizeOld + iCopy);
2774  partonSystemsPtr->addOut(iSysSel, eventSizeOld + systemSizeOld);
2775  partonSystemsPtr->setSHat(iSysSel, m2II / z);
2776 
2777  // Add dipoles for q -> g q, where the daughter is the gluon.
2778  if (idDaughter == 21 && idMother != 21) {
2779  if (doQEDshowerByQ) {
2780  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2781  iNewRec, pT2, 0, mother.chargeType(), 0, 0, normalRecoil) );
2782  }
2783  if (doWeakShower && iSysSel == 0) {
2784  int MEtypeNew = 203;
2785  if (idRecoiler == 21) MEtypeNew = 201;
2786  if (idRecoiler == idMother) MEtypeNew = 202;
2787  // If original was a Drell-Yan, keep as Drell-Yan.
2788  if( event[3].id() == - event[4].id()) MEtypeNew = 200;
2789  int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
2790  event[iMother].pol(weakPol);
2791  if ((weakMode == 0 || weakMode == 1) && weakPol == -1)
2792  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2793  iNewRec, pT2, 0, 0, 1, MEtypeNew, normalRecoil, weakPol) );
2794  if (weakMode == 0 || weakMode == 2)
2795  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2796  iNewRec, pT2, 0, 0, 2, MEtypeNew + 5, normalRecoil, weakPol) );
2797  }
2798  }
2799 
2800  // Add dipoles for q -> q gamma, where the daughter is the gamma.
2801  if (idDaughter == 22 && idMother != 22) {
2802  if (doQCDshower && mother.colType() != 0) {
2803  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2804  iNewRec, pT2, mother.colType(), 0, 0, 0, normalRecoil) );
2805  }
2806  if (doQEDshowerByQ && mother.chargeType() != 3) {
2807  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2808  iNewRec, pT2, 0, mother.chargeType(), 0, 0, normalRecoil) );
2809  }
2810  if (doQEDshowerByL && mother.chargeType() == 3) {
2811  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2812  iNewRec, pT2, 0, mother.chargeType(), 0, 0, normalRecoil) );
2813  }
2814  if (doWeakShower && iSysSel == 0) {
2815  int MEtypeNew = 203;
2816  if (idRecoiler == 21) MEtypeNew = 201;
2817  if (idRecoiler == idMother) MEtypeNew = 202;
2818  // If original was a Drell-Yan, keep as Drell-Yan.
2819  if( event[3].id() == - event[4].id()) MEtypeNew = 200;
2820  int weakPol = (rndmPtr->flat() > 0.5) ? -1 : 1;
2821  event[iMother].pol(weakPol);
2822  if ((weakMode == 0 || weakMode == 1) && weakPol == -1)
2823  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2824  iNewRec, pT2, 0, 0, 1, MEtypeNew, normalRecoil, weakPol) );
2825  if (weakMode == 0 || weakMode == 2)
2826  dipEnd.push_back( SpaceDipoleEnd( iSysSel, side, iMother,
2827  iNewRec, pT2, 0, 0, 2, MEtypeNew + 5, normalRecoil, weakPol) );
2828  }
2829  }
2830 
2831  // dipEnd array may have expanded and been moved, so regenerate dipEndSel.
2832  dipEndSel = &dipEnd[iDipSel];
2833 
2834  // Set flag to tell that a weak emission has happened.
2835  if (dipEndSel->weakType != 0) hasWeaklyRadiated = true;
2836 
2837  // Update list of QCD emissions in side A and B in given iSysSel
2838  // This is used to veto jets in W/z events.
2839  while (iSysSel >= int(nRadA.size()) || iSysSel >= int(nRadB.size())) {
2840  nRadA.push_back(0);
2841  nRadB.push_back(0);
2842  }
2843  if (dipEndSel->colType != 0 && side == 1) ++nRadA[iSysSel];
2844  else if (dipEndSel->colType != 0) ++nRadB[iSysSel];
2845 
2846  // Update info on radiating dipole ends (QCD, QED or weak).
2847  for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
2848  if ( dipEnd[iDip].system == iSysSel) {
2849  if (abs(dipEnd[iDip].side) == side) {
2850  dipEnd[iDip].iRadiator = iMother;
2851  dipEnd[iDip].iRecoiler = iNewRec;
2852  // Look if there is a IF dipole in case of dipole recoil.
2853  dipEnd[iDip].iColPartner = (doDipoleRecoil) ? findColPartner( event,
2854  dipEnd[iDip].iRadiator, dipEnd[iDip].iRecoiler, iSysSel) : 0;
2855  dipEnd[iDip].idColPartner = (dipEnd[iDip].iColPartner != 0)
2856  ? event[dipEnd[iDip].iColPartner].id() : 0;
2857  if (dipEnd[iDip].colType != 0)
2858  dipEnd[iDip].colType = mother.colType();
2859  else if (dipEnd[iDip].chgType != 0) {
2860  dipEnd[iDip].chgType = 0;
2861  if ( (mother.isQuark() && doQEDshowerByQ)
2862  || (mother.isLepton() && doQEDshowerByL) )
2863  dipEnd[iDip].chgType = mother.chargeType();
2864  }
2865  else if (dipEnd[iDip].weakType != 0) {
2866  // Kill weak dipole if mother becomes gluon / photon.
2867  if (!(mother.isLepton() || mother.isQuark()))
2868  dipEnd[iDip].weakType = 0;
2869  if (singleWeakEmission && hasWeaklyRadiated)
2870  dipEnd[iDip].weakType = 0;
2871  }
2872 
2873  // Kill ME corrections after first emission for everything
2874  // but weak showers.
2875  if (dipEnd[iDip].weakType == 0) dipEnd[iDip].MEtype = 0;
2876 
2877  // Update info on recoiling dipole ends (QCD or QED).
2878  } else {
2879  dipEnd[iDip].iRadiator = iNewRec;
2880  dipEnd[iDip].iRecoiler = iMother;
2881  // Look if there is an IF dipole in case of dipole recoil.
2882  dipEnd[iDip].iColPartner = (doDipoleRecoil) ? findColPartner( event,
2883  dipEnd[iDip].iRadiator, dipEnd[iDip].iRecoiler, iSysSel) : 0;
2884  dipEnd[iDip].idColPartner = (dipEnd[iDip].iColPartner != 0)
2885  ? event[dipEnd[iDip].iColPartner].id() : 0;
2886  // Optionally also kill recoiler ME corrections after first emission.
2887  if (!doMEafterFirst && dipEnd[iDip].weakType == 0)
2888  dipEnd[iDip].MEtype = 0;
2889  // Remove weak dipoles if we only want a single emission.
2890  if (dipEnd[iDip].weakType != 0 && singleWeakEmission
2891  && hasWeaklyRadiated) dipEnd[iDip].weakType = 0;
2892  }
2893  }
2894 
2895  // Set polarisation of mother for weak emissions.
2896  if (dipEndSel->weakType != 0) mother.pol(dipEndSel->weakPol);
2897 
2898  // Update info on beam remnants.
2899  double xNew = (side == 1) ? x1New : x2New;
2900  beamNow[iSysSel].update( iMother, idMother, xNew);
2901  // Redo choice of companion kind whenever new flavour.
2902  if (idMother != idDaughterNow) {
2903  pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
2904  beamNow.xfISR( iSysSel, idMother, xNew, pdfScale2);
2905  beamNow.pickValSeaComp();
2906  }
2907  BeamParticle& beamRec = (side == 1) ? *beamBPtr : *beamAPtr;
2908  beamRec[iSysSel].iPos( iNewRec);
2909 
2910  // Store branching values of current dipole. (For rapidity ordering.)
2911  ++dipEndSel->nBranch;
2912  dipEndSel->pT2Old = pT2;
2913  dipEndSel->zOld = z;
2914 
2915  // Update history if recoiler rescatters.
2916  if (!normalRecoil)
2917  event[iRecoilMother].daughters( iNewRec, iNewRec);
2918 
2919  // Start list of rescatterers that force changed kinematics.
2920  vector<int> iRescatterer;
2921  for ( int i = 0; i < systemSizeOld - 2; ++i) {
2922  int iOutNew = partonSystemsPtr->getOut( iSysSel, i);
2923  if (!event[iOutNew].isFinal()) iRescatterer.push_back(iOutNew);
2924  }
2925 
2926  // Start iterate over list of such rescatterers.
2927  int iRescNow = -1;
2928  while (++iRescNow < int(iRescatterer.size())) {
2929 
2930  // Identify partons that induce or are affected by rescatter shift.
2931  // In following Old is before change of kinematics, New after,
2932  // Out scatterer in outstate and In in instate of another system.
2933  // Daughter sequence is (iOutOld ->) iOutNew -> iInNew -> iInOld.
2934  int iOutNew = iRescatterer[iRescNow];
2935  int iInOld = event[iOutNew].daughter1();
2936  int iSysResc = partonSystemsPtr->getSystemOf(iInOld, true);
2937 
2938  // Copy incoming partons of rescattered system and hook them up.
2939  int iOldA = partonSystemsPtr->getInA(iSysResc);
2940  int iOldB = partonSystemsPtr->getInB(iSysResc);
2941  bool rescSideA = event[iOldA].isRescatteredIncoming();
2942  int statusNewA = (rescSideA) ? -45 : -42;
2943  int statusNewB = (rescSideA) ? -42 : -45;
2944  int iNewA = event.copy(iOldA, statusNewA);
2945  int iNewB = event.copy(iOldB, statusNewB);
2946 
2947  // Copy outgoing partons of rescattered system and hook them up.
2948  int eventSize = event.size();
2949  int sizeOutAB = partonSystemsPtr->sizeOut(iSysResc);
2950  int iOldAB, statusOldAB, iNewAB;
2951  for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB) {
2952  iOldAB = partonSystemsPtr->getOut(iSysResc, iOutAB);
2953  statusOldAB = event[iOldAB].status();
2954  iNewAB = event.copy(iOldAB, 44);
2955  // Status could be negative for parton that rescatters in its turn.
2956  if (statusOldAB < 0) {
2957  event[iNewAB].statusNeg();
2958  iRescatterer.push_back(iNewAB);
2959  }
2960  }
2961 
2962  // Hook up new outgoing with new incoming parton.
2963  int iInNew = (rescSideA) ? iNewA : iNewB;
2964  event[iOutNew].daughters( iInNew, iInNew);
2965  event[iInNew].mothers( iOutNew, iOutNew);
2966 
2967  // Rescale recoiling incoming parton for correct invariant mass.
2968  event[iInNew].p( event[iOutNew].p() );
2969  double momFac = (rescSideA)
2970  ? event[iInOld].pPos() / event[iInNew].pPos()
2971  : event[iInOld].pNeg() / event[iInNew].pNeg();
2972  int iInRec = (rescSideA) ? iNewB : iNewA;
2973 
2974  // Rescatter: A previous boost may cause the light cone momentum of a
2975  // rescattered parton to change sign. If this happens, tell
2976  // parton level to try again.
2977  if (momFac < 0.0) {
2978  infoPtr->errorMsg("Warning in SpaceShower::branch: "
2979  "change in lightcone momentum sign; retrying parton level");
2980  rescatterFail = true;
2981  return false;
2982  }
2983  event[iInRec].rescale4( momFac);
2984 
2985  // Boost outgoing partons to new frame of incoming.
2986  RotBstMatrix MmodResc;
2987  MmodResc.toCMframe( event[iOldA].p(), event[iOldB].p());
2988  MmodResc.fromCMframe(event[iNewA].p(), event[iNewB].p());
2989  for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB)
2990  event[eventSize + iOutAB].rotbst(MmodResc);
2991 
2992  // Update list of partons in system.
2993  partonSystemsPtr->setInA(iSysResc, iNewA);
2994  partonSystemsPtr->setInB(iSysResc, iNewB);
2995  for (int iCopy = 0; iCopy < sizeOutAB; ++iCopy)
2996  partonSystemsPtr->setOut(iSysResc, iCopy, eventSize + iCopy);
2997 
2998  // Update info on radiating dipole ends (QCD or QED).
2999  for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
3000  if ( dipEnd[iDip].system == iSysResc) {
3001  bool sideAnow = (abs(dipEnd[iDip].side) == 1);
3002  dipEnd[iDip].iRadiator = (sideAnow) ? iNewA : iNewB;
3003  dipEnd[iDip].iRecoiler = (sideAnow) ? iNewB : iNewA;
3004  }
3005 
3006  // Update info on beam remnants.
3007  BeamParticle& beamResc = (rescSideA) ? *beamAPtr : *beamBPtr;
3008  beamResc[iSysResc].iPos( iInNew);
3009  beamResc[iSysResc].p( event[iInNew].p() );
3010  beamResc[iSysResc].scaleX( 1. / momFac );
3011  BeamParticle& beamReco = (rescSideA) ? *beamBPtr : *beamAPtr;
3012  beamReco[iSysResc].iPos( iInRec);
3013  beamReco[iSysResc].scaleX( momFac);
3014 
3015  // End iterate over list of rescatterers.
3016  }
3017 
3018  // Check that beam momentum not used up by rescattered-system boosts.
3019  if ( ( beamAPtr->xMax(-1) < 0.0 && !(beamAPtr->isUnresolved()) )
3020  || (beamBPtr->xMax(-1) < 0.0 && !(beamBPtr->isUnresolved()) ) ) {
3021  infoPtr->errorMsg("Warning in SpaceShower::branch: "
3022  "used up beam momentum; retrying parton level");
3023  rescatterFail = true;
3024  return false;
3025  }
3026 
3027  // If gamma -> q qbar valid with photon beam no need for remnants.
3028  if ( beamNow.isGamma() && beamNow.resolvedGamma() && gamma2qqbar)
3029  beamNow.resolvedGamma(false);
3030 
3031  // Done without any errors.
3032  return true;
3033 
3034 }
3035 
3036 //--------------------------------------------------------------------------
3037 
3038 // Initialize the choices of uncertainty variations of the shower.
3039 
3040 bool SpaceShower::initUncertainties() {
3041 
3042  // Populate lists of uncertainty variations for SpaceShower, by keyword.
3043  uVarMuSoftCorr = settingsPtr->flag("UncertaintyBands:muSoftCorr");
3044  dASmax = settingsPtr->parm("UncertaintyBands:deltaAlphaSmax");
3045 
3046  // Reset uncertainty variation maps.
3047  varG2GGmuRfac.clear(); varG2GGcNS.clear();
3048  varQ2QGmuRfac.clear(); varQ2QGcNS.clear();
3049  varQ2GQmuRfac.clear(); varQ2GQcNS.clear();
3050  varX2XGmuRfac.clear(); varX2XGcNS.clear();
3051  varG2QQmuRfac.clear(); varG2QQcNS.clear();
3052  // Maps that must be known by TimeShower
3053  varPDFplus = &infoPtr->varPDFplus;
3054  varPDFminus = &infoPtr->varPDFminus;
3055  varPDFmember = &infoPtr->varPDFmember;
3056  varPDFplus->clear(); varPDFminus->clear();
3057  varPDFmember->clear();
3058 
3059  vector<string> keys;
3060  // List of keywords recognised by SpaceShower.
3061  keys.push_back("isr:murfac");
3062  keys.push_back("isr:g2gg:murfac");
3063  keys.push_back("isr:q2qg:murfac");
3064  keys.push_back("isr:q2gq:murfac");
3065  keys.push_back("isr:x2xg:murfac");
3066  keys.push_back("isr:g2qq:murfac");
3067  keys.push_back("isr:cns");
3068  keys.push_back("isr:g2gg:cns");
3069  keys.push_back("isr:q2qg:cns");
3070  keys.push_back("isr:q2gq:cns");
3071  keys.push_back("isr:x2xg:cns");
3072  keys.push_back("isr:g2qq:cns");
3073  keys.push_back("isr:pdf:plus");
3074  keys.push_back("isr:pdf:minus");
3075  keys.push_back("isr:pdf:member");
3076 
3077  // Store number of QCD variations (as separator to QED ones).
3078  int nKeysQCD=keys.size();
3079 
3080  // Get uncertainty variations from Settings (as list of strings to parse).
3081  vector<string> uVars = settingsPtr->wvec("UncertaintyBands:List");
3082  size_t varSize = uVars.size();
3083  nUncertaintyVariations = int(uVars.size());
3084  if (nUncertaintyVariations == 0) return false;
3085  vector<string> uniqueVars;
3086 
3087  // Expand uVars if PDFmembers has been chosen
3088  string tmpKey("isr:pdf:family");
3089  // Parse each string in uVars to look for recognized keywords.
3090  for (size_t iWeight = 0; iWeight < varSize; ++iWeight) {
3091  // Convert to lowercase (to be case-insensitive). Also remove "=" signs
3092  // and extra spaces, so "key=value", "key = value" mapped to "key value"
3093  string uVarString = toLower(uVars[iWeight]);
3094  while (uVarString.find(" ") == 0) uVarString.erase( 0, 1);
3095  int iEnd = uVarString.find(" ", 0);
3096  uVarString.erase(0,iEnd+1);
3097  while (uVarString.find("=") != string::npos) {
3098  int firstEqual = uVarString.find_first_of("=");
3099  string testString = uVarString.substr(0, firstEqual);
3100  iEnd = uVarString.find_first_of(" ", 0);
3101  if( iEnd<0 ) iEnd = uVarString.length();
3102  string insertString = uVarString.substr(0,iEnd);
3103  // does the key match an fsr one?
3104  if( find(keys.begin(), keys.end(), testString) != keys.end() ) {
3105  if( uniqueVars.size() == 0 ) {
3106  uniqueVars.push_back(insertString);
3107  } else if ( find(uniqueVars.begin(), uniqueVars.end(), insertString)
3108  == uniqueVars.end() ) {
3109  uniqueVars.push_back(insertString);
3110  }
3111  } else if ( testString == tmpKey ) {
3112  int nMembers(0);
3113  BeamParticle& beam = *beamAPtr;
3114  nMembers = beam.nMembers();
3115  for(int iMem=1; iMem<nMembers; ++iMem) {
3116  ostringstream iss;
3117  iss << iMem;
3118  string tmpString("isr:pdf:member="+iss.str());
3119  if (find(uniqueVars.begin(), uniqueVars.end(), tmpString)
3120  == uniqueVars.end() ) {
3121  uniqueVars.push_back(tmpString);
3122  }
3123  }
3124  }
3125  uVarString.erase(0,iEnd+1);
3126  }
3127  }
3128 
3129  nUncertaintyVariations = int(uniqueVars.size());
3130 
3131  if ( nUncertaintyVariations > 0 ) {
3132  int nWeights = infoPtr->nWeights();
3133  infoPtr->setNWeights( nWeights + nUncertaintyVariations );
3134  int newSize = infoPtr->nWeights();
3135  for(int iWeight = nWeights; iWeight < newSize; ++iWeight) {
3136  string uVarString = uniqueVars[iWeight - nWeights];
3137  infoPtr->setWeightLabel(iWeight, uVarString);
3138  // Parse each string in uVars to look for recognised keywords.
3139  // Convert to lowercase (to be case-insensitive). Also remove "=" signs
3140  // and extra spaces, so "key=value", "key = value" mapped to "key value"
3141 
3142  while (uVarString.find("=") != string::npos) {
3143  int firstEqual = uVarString.find_first_of("=");
3144  uVarString.replace(firstEqual, 1, " ");
3145  }
3146  while (uVarString.find(" ") != string::npos)
3147  uVarString.erase( uVarString.find(" "), 1);
3148  if (uVarString == "" || uVarString == " ") continue;
3149 
3150  // Loop over all keywords.
3151  int nRecognizedQCD = 0;
3152  for (int iWord = 0; iWord < int(keys.size()); ++iWord) {
3153  // Transform string to lowercase to avoid case-dependence.
3154  string key = toLower(keys[iWord]);
3155  // Extract variation value/factor.
3156  int iKey = uVarString.find(key);
3157  int iBeg = uVarString.find(" ", iKey) + 1;
3158  int iEnd = uVarString.find(" ", iBeg);
3159  string valueString = uVarString.substr(iBeg, iEnd - iBeg);
3160  stringstream ss(valueString);
3161  double value;
3162  ss >> value;
3163  if (!ss) continue;
3164 
3165  // Store (iWeight,value) pairs
3166  // RECALL: use lowercase for all keys here (since converted above).
3167  if (key == "isr:murfac" || key == "isr:g2gg:murfac")
3168  varG2GGmuRfac[iWeight] = value;
3169  if (key == "isr:murfac" || key == "isr:q2qg:murfac")
3170  varQ2QGmuRfac[iWeight] = value;
3171  if (key == "isr:murfac" || key == "isr:q2gq:murfac")
3172  varQ2GQmuRfac[iWeight] = value;
3173  if (key == "isr:murfac" || key == "isr:x2xg:murfac")
3174  varX2XGmuRfac[iWeight] = value;
3175  if (key == "isr:murfac" || key == "isr:g2qq:murfac")
3176  varG2QQmuRfac[iWeight] = value;
3177  if (key == "isr:cns" || key == "isr:g2gg:cns")
3178  varG2GGcNS[iWeight] = value;
3179  if (key == "isr:cns" || key == "isr:q2qg:cns")
3180  varQ2QGcNS[iWeight] = value;
3181  if (key == "isr:cns" || key == "isr:q2gq:cns")
3182  varQ2GQcNS[iWeight] = value;
3183  if (key == "isr:cns" || key == "isr:x2xg:cns")
3184  varX2XGcNS[iWeight] = value;
3185  if (key == "isr:cns" || key == "isr:g2qq:cns")
3186  varG2QQcNS[iWeight] = value;
3187  if (key == "isr:pdf:plus") varPDFplus->operator[](iWeight) = 1;
3188  if (key == "isr:pdf:minus") varPDFminus->operator[](iWeight) = 1;
3189  if (key == "isr:pdf:member")
3190  varPDFmember->operator[](iWeight) = int(value);
3191  // Tell that we found at least one recognized and parseable keyword.
3192  if (iWord < nKeysQCD) nRecognizedQCD++;
3193  } // End loop over QCD keywords
3194 
3195  // Tell whether this uncertainty variation contained >= 1 QCD variation.
3196  if (nRecognizedQCD > 0) ++nVarQCD;
3197  } // End loop over UVars.
3198  }
3199 
3200  infoPtr->initUncertainties(&uVars,true);
3201  // Let the calling function know if we found anything.
3202  return (nVarQCD > 0);
3203 }
3204 
3205 //--------------------------------------------------------------------------
3206 
3207 // Calculate uncertainties for the current event.
3208 
3209 void SpaceShower::calcUncertainties(bool accept, double pAccept, double pT20in,
3210  double enhance, double vp, SpaceDipoleEnd* dip, Particle* motPtr,
3211  Particle* sisPtr) {
3212 
3213  // Sanity check.
3214  if (!doUncertainties || !doUncertaintiesNow || nUncertaintyVariations <= 0)
3215  return;
3216 
3217  // Define pointer and iterator to loop over the contents of each
3218  // (iWeight,value) map.
3219  map<int,double>* varPtr=0;
3220  map<int,double>::iterator itVar;
3221  // Make sure we have a dummy to point to if no map to be used.
3222  map<int,double> dummy; dummy.clear();
3223 
3224  int numWeights = infoPtr->nWeights();
3225  // Store uncertainty variation factors, initialised to unity.
3226  // Make vector sizes + 1 since 0 = default and variations start at 1.
3227  vector<double> uVarFac(numWeights, 1.0);
3228  vector<bool> doVar(numWeights, false);
3229  // When performing biasing, the nominal weight need not be unity.
3230  doVar[0] = true;
3231  uVarFac[0] = 1.0;
3232 
3233  // Extract IDs, with standard ISR nomenclature: mot -> dau(Q2) + sis
3234  int idSis = sisPtr->id();
3235  int idMot = motPtr->id();
3236 
3237  // PDF variations
3238  if ( !varPDFplus->empty() || !varPDFminus->empty() || !varPDFmember->empty()
3239  ) {
3240  // Evaluation of new daughter and mother PDF's.
3241  double scale2 = (useFixedFacScale) ? fixedFacScale2
3242  : factorMultFac * dip->pT2;
3243  double xMother = dip->xMo;
3244  double xDau = dip->z * xMother;
3245  BeamParticle& beam = (abs(dip->side) == 1) ? *beamAPtr : *beamBPtr;
3246  int valSea = (beam[iSysSel].isValence()) ? 1 : 0;
3247  if( beam[iSysSel].isUnmatched() ) valSea = 2;
3248  beam.calcPDFEnvelope( make_pair(dip->idMother,dip->idDaughter),
3249  make_pair(xMother,xDau), scale2, valSea);
3250  PDF::PDFEnvelope ratioPDFEnv = beam.getPDFEnvelope( );
3251  //
3252  varPtr = varPDFplus;
3253  for (itVar = varPtr->begin(); itVar != varPtr->end(); ++itVar) {
3254  int iWeight = itVar->first;
3255  uVarFac[iWeight] *= 1.0 + min(ratioPDFEnv.errplusPDF
3256  / ratioPDFEnv.centralPDF, 0.5);
3257  doVar[iWeight] = true;
3258  }
3259  //
3260  varPtr = varPDFminus;
3261  for (itVar = varPtr->begin(); itVar != varPtr->end(); ++itVar) {
3262  int iWeight = itVar->first;
3263  uVarFac[iWeight] *= max(.01,1.0 - min(ratioPDFEnv.errminusPDF
3264  / ratioPDFEnv.centralPDF, 0.5));
3265  doVar[iWeight] = true;
3266  }
3267  varPtr = varPDFmember;
3268  for (itVar = varPtr->begin(); itVar != varPtr->end(); ++itVar) {
3269  int iWeight = itVar->first;
3270  int member = int( itVar->second );
3271  uVarFac[iWeight] *= max(.01,ratioPDFEnv.pdfMemberVars[member]
3272  / ratioPDFEnv.centralPDF);
3273  doVar[iWeight] = true;
3274  }
3275  }
3276 
3277  // QCD variations.
3278  if (dip->colType != 0) {
3279 
3280  // QCD renormalization-scale variations.
3281  if (alphaSorder == 0) varPtr = &dummy;
3282  else if (idMot == 21 && idSis == 21) varPtr = &varG2GGmuRfac;
3283  else if (idMot == 21 && abs(idSis) <= nQuarkIn) varPtr = &varG2QQmuRfac;
3284  else if (abs(idMot) <= nQuarkIn) {
3285  if (abs(idMot) <= uVarNflavQ) varPtr = &varQ2QGmuRfac;
3286  else varPtr = &varX2XGmuRfac;
3287  }
3288  else varPtr = &dummy;
3289  double Q2 = dip->pT2;
3290  double muR2 = renormMultFac * (Q2 + pT20in);
3291  double alphaSbaseline = alphaS.alphaS(muR2);
3292  for (itVar = varPtr->begin(); itVar != varPtr->end(); ++itVar) {
3293  int iWeight = itVar->first;
3294  double valFac = itVar->second;
3295  // Correction-factor alphaS.
3296  double muR2var = max(1.1 * Lambda3flav2, pow2(valFac) * muR2);
3297  double alphaSratio = alphaS.alphaS(muR2var) / alphaSbaseline;
3298  // Apply soft correction factor only for (on-shell) gluon emission
3299  double facCorr = 1.;
3300  if (idSis == 21 && uVarMuSoftCorr) {
3301  // Use smallest alphaS and b0, to make the compensation conservative.
3302  int nf = 5;
3303  if (dip->pT2 < pow2(mc)) nf = 3;
3304  else if (dip->pT2 < pow2(mb)) nf = 4;
3305  double alphaScorr = alphaS.alphaS(dip->m2Dip);
3306  double facSoft = alphaScorr * (33. - 2. * nf) / (6. * M_PI);
3307  // Zeta is energy fraction of emitted (on-shell) gluon = 1 - z.
3308  double zeta = 1. - dip->z;
3309  facCorr = 1. + (1. - zeta) * facSoft * log(valFac);
3310  }
3311  // Apply correction factor here for emission processes.
3312  double alphaSfac = alphaSratio * facCorr;
3313  // Limit absolute variation to +/- deltaAlphaSmax
3314  if (alphaSfac > 1.)
3315  alphaSfac = min(alphaSfac, (alphaSbaseline + dASmax) / alphaSbaseline);
3316  else if (alphaSbaseline > dASmax)
3317  alphaSfac = max(alphaSfac, (alphaSbaseline - dASmax) / alphaSbaseline);
3318  uVarFac[iWeight] *= alphaSfac;
3319  doVar[iWeight] = true;
3320  }
3321 
3322  // QCD finite-term variations (only when no MECs and above pT threshold).
3323  if (dip->MEtype != 0 || dip->pT2 < pow2(cNSpTmin) ) varPtr = &dummy;
3324  else if (idMot == 21 && idSis == 21) varPtr = &varG2GGcNS;
3325  else if (idMot == 21 && abs(idSis) <= nQuarkIn) varPtr = &varG2QQcNS;
3326  else if (abs(idMot) <= nQuarkIn) {
3327  if (abs(idMot) <= uVarNflavQ) varPtr = &varQ2QGcNS;
3328  else varPtr = &varX2XGcNS;
3329  }
3330  else varPtr = &dummy;
3331  double z = dip->z;
3332  for (itVar = varPtr->begin(); itVar != varPtr->end(); ++itVar) {
3333  int iWeight = itVar->first;
3334  double valFac = itVar->second;
3335  // Correction-factor alphaS.
3336  // Virtuality for off-shell massive quarks.
3337  if (idMot == 21 && abs(idSis) >= 4 && idSis != 21)
3338  Q2 = max(1., Q2+pow2(sisPtr->m0()));
3339  else if (idSis == 21 && abs(idMot) >= 4 && idMot != 21)
3340  Q2 = max(1., Q2+pow2(motPtr->m0()));
3341  double yQ = Q2 / dip->m2Dip;
3342  double num = yQ * valFac;
3343  double denom = 1.;
3344  // G->GG.
3345  if (idSis == 21 && idMot == 21)
3346  denom = pow2(1. - z * (1.-z)) / (z*(1.-z));
3347  // Q->QG.
3348  else if (idSis == 21)
3349  denom = (1. + pow2(z)) / (1. - z);
3350  // Q->GQ.
3351  else if (idMot == idSis)
3352  denom = (1. + pow2(1. - z)) / z;
3353  // G->QQ.
3354  else
3355  denom = pow2(z) + pow2(1. - z);
3356  // Compute reweight ratio.
3357  double minReWeight = max( 1. + num / denom, REJECTFACTOR );
3358  uVarFac[iWeight] *= minReWeight;
3359  doVar[iWeight] = true;
3360  }
3361  }
3362 
3363  // Ensure 0 < PacceptPrime < 1 (with small margins).
3364  // Skip the central weight, so as to avoid confusion
3365  for (int iWeight = 1; iWeight<numWeights; ++iWeight) {
3366  if (!doVar[iWeight]) continue;
3367  double pAcceptPrime = pAccept * uVarFac[iWeight];
3368  if (pAcceptPrime > PROBLIMIT && dip->colType != 0) {
3369  uVarFac[iWeight] *= PROBLIMIT / pAcceptPrime;
3370  }
3371  }
3372 
3373  // Apply reject or accept reweighting factors according to input decision.
3374  for (int iWeight = 0; iWeight < numWeights; ++iWeight) {
3375  if (!doVar[iWeight]) continue;
3376  // If trial accepted: apply ratio of accept probabilities.
3377  if (accept) infoPtr->reWeight( iWeight,
3378  uVarFac[iWeight] / ((1.0 - vp) * enhance) );
3379  // If trial rejected : apply Sudakov reweightings.
3380  else {
3381  // Check for near-singular denominators (indicates too few failures,
3382  // and hence would need to increase headroom).
3383  double denom = 1. - pAccept * (1.0 - vp);
3384  if (denom < REJECTFACTOR) {
3385  stringstream message;
3386  message << iWeight;
3387  infoPtr->errorMsg("Warning in SpaceShower: reject denom for"
3388  " iWeight = ", message.str());
3389  }
3390  // Force reweighting factor > 0.
3391  double reWtFail = max(0.01, (1. - uVarFac[iWeight] * pAccept / enhance )
3392  / denom);
3393  infoPtr->reWeight(iWeight, reWtFail);
3394  }
3395  }
3396 }
3397 
3398 //--------------------------------------------------------------------------
3399 
3400 // Find class of ME correction.
3401 
3402  int SpaceShower::findMEtype( int iSys, Event& event, bool weakRadiation) {
3403 
3404  // Default values and no action.
3405  int MEtype = 0;
3406  if (!doMEcorrections) return MEtype;
3407 
3408  // Identify systems producing a single resonance.
3409  if (partonSystemsPtr->sizeOut( iSys) == 1 && !weakRadiation) {
3410  int idIn1 = event[partonSystemsPtr->getInA(iSys)].id();
3411  int idIn2 = event[partonSystemsPtr->getInA(iSys)].id();
3412  int idRes = event[partonSystemsPtr->getOut(iSys, 0)].id();
3413  if (iSys == 0) idResFirst = abs(idRes);
3414  if (iSys == 1) idResSecond = abs(idRes);
3415 
3416  // f + fbar -> vector boson.
3417  if ( (idRes == 23 || abs(idRes) == 24 || idRes == 32
3418  || idRes == 33 || abs(idRes) == 34 || abs(idRes) == 41)
3419  && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 1;
3420 
3421  // g + g, gamma + gamma -> Higgs boson.
3422  if ( (idRes == 25 || idRes == 35 || idRes == 36)
3423  && ( ( idIn1 == 21 && idIn2 == 21 )
3424  || ( idIn1 == 22 && idIn2 == 22 ) ) ) MEtype = 2;
3425 
3426  // f + fbar -> Higgs boson.
3427  if ( (idRes == 25 || idRes == 35 || idRes == 36)
3428  && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 3;
3429  }
3430 
3431  // Weak ME corrections.
3432  if (weakRadiation) {
3433  if (event[3].id() == -event[4].id()
3434  || event[event[3].daughter1()].idAbs() == 24 || infoPtr->nFinal() != 2)
3435  MEtype = 200;
3436  else if (event[3].idAbs() == 21 || event[4].idAbs() == 21) MEtype = 201;
3437  else if (event[3].id() == event[4].id()) MEtype = 202;
3438  else MEtype = 203;
3439  }
3440 
3441  // Done.
3442  return MEtype;
3443 
3444 }
3445 
3446 //--------------------------------------------------------------------------
3447 
3448 // Provide maximum of expected ME weight; for preweighting of evolution.
3449 
3450 double SpaceShower::calcMEmax( int MEtype, int idMother, int idDaughterIn) {
3451 
3452  // Main non-unity case: g(gamma) f -> V f'.
3453  if (MEtype == 1 && idMother > 20 && idDaughterIn < 20) return 3.;
3454 
3455  // Added a case for t-channel W/Z exchange, since the PS is not an
3456  // overestimate. This does not help fully, but it should only be small
3457  // pT quarks / gluons that break the overscattering.
3458  if ( MEtype == 201 || MEtype == 202 || MEtype == 203
3459  || MEtype == 206 || MEtype == 207 || MEtype == 208) return WEAKPSWEIGHT;
3460 
3461  // Default.
3462  return 1.;
3463 
3464 }
3465 
3466 //--------------------------------------------------------------------------
3467 
3468 // Provide actual ME weight for current branching.
3469 // Note: currently ME corrections are only allowed for first branching
3470 // on each side, so idDaughter is essentially known and checks overkill.
3471 
3472 double SpaceShower::calcMEcorr(int MEtype, int idMother, int idDaughterIn,
3473  double M2, double z, double Q2, double m2Sister) {
3474 
3475  // Convert to Mandelstam variables. Sometimes may need to swap later.
3476  double sH = M2 / z;
3477  double tH = -Q2;
3478  double uH = Q2 - M2 * (1. - z) / z;
3479  int idMabs = abs(idMother);
3480  int idDabs = abs(idDaughterIn);
3481 
3482  // Corrections for f + fbar -> s-channel vector boson.
3483  if (MEtype == 1) {
3484  if (idMabs < 20 && idDabs < 20) {
3485  return (tH*tH + uH*uH + 2. * M2 * sH) / (sH*sH + M2*M2);
3486  } else if (idDabs < 20) {
3487  // g(gamma) f -> V f': -Q2 = (p_g - p_f')^2 in PS while
3488  // tHat = (p_f - p_f')^2 in ME so need to swap tHat <-> uHat.
3489  swap( tH, uH);
3490  return (sH*sH + uH*uH + 2. * M2 * tH) / (pow2(sH - M2) + M2*M2);
3491  }
3492 
3493  // Corrections for g + g -> Higgs boson.
3494  } else if (MEtype == 2) {
3495  if (idMabs < 20 && idDabs > 20) {
3496  return (sH*sH + uH*uH) / (sH*sH + pow2(sH - M2));
3497  } else if (idDabs > 20) {
3498  return 0.5 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(M2))
3499  / pow2(sH*sH - M2 * (sH - M2));
3500  }
3501 
3502  // Corrections for f + fbar -> Higgs boson (f = b mainly).
3503  } else if (MEtype == 3) {
3504  if (idMabs < 20 && idDabs < 20) {
3505  // The PS and ME answers agree for f fbar -> H g/gamma.
3506  return 1.;
3507  } else if (idDabs < 20) {
3508  // Need to swap tHat <-> uHat, cf. vector-boson production above.
3509  swap( tH, uH);
3510  return (sH*sH + uH*uH + 2. * (M2 - uH) * (M2 - sH))
3511  / (pow2(sH - M2) + M2*M2);
3512  }
3513 
3514  // Corrections for f -> f' + W/Z (s-channel).
3515  } else if (MEtype == 200 || MEtype == 205) {
3516  // Need to redo calculations of uH since we now emit a massive particle.
3517  uH += m2Sister;
3518  double wtME = (uH*uH + tH*tH + 2 * sH * (m2Sister + M2)) / (uH*tH)
3519  - M2 * m2Sister * (1/(tH*tH) + 1/(uH*uH));
3520  double wtPS = (sH*sH + pow2(M2 + m2Sister)) / (tH*uH);
3521  return wtME / wtPS;
3522  } else if (MEtype == 201 || MEtype == 202 || MEtype == 203 ||
3523  MEtype == 206 || MEtype == 207 || MEtype == 208)
3524  return calcMEmax(MEtype, 0, 0);
3525 
3526  // Default.
3527  return 1.;
3528 
3529 }
3530 
3531 //--------------------------------------------------------------------------
3532 
3533 // Provide actual ME weight for current branching for weak t-channel emissions.
3534 
3535 double SpaceShower::calcMEcorrWeak(int MEtype, double m2, double z,
3536  double pT2, Vec4 pMother, Vec4 pB, Vec4 pDaughter,
3537  Vec4 pB0, Vec4 p1, Vec4 p2, Vec4 pSister) {
3538 
3539  // Find daughter four-momentum in current frame.
3540  Vec4 pA = pMother - pSister;
3541 
3542  // Scale outgoing vectors to conserve energy / momentum.
3543  double scaleFactor2 = (pA + pB).m2Calc() / (p1 + p2).m2Calc();
3544  double scaleFactor = sqrt(scaleFactor2);
3545  RotBstMatrix rot2to2frame;
3546  rot2to2frame.bstback(p1 + p2);
3547  p1.rotbst(rot2to2frame);
3548  p2.rotbst(rot2to2frame);
3549  p1 *= scaleFactor;
3550  p2 *= scaleFactor;
3551 
3552  // Find 2 to 2 rest frame for incoming particles.
3553  // This is done before one of the two are made virtual (Q^2 mass).
3554  RotBstMatrix rot2to2frameInc;
3555  rot2to2frameInc.bstback(pDaughter + pB0);
3556  pDaughter.rotbst(rot2to2frameInc);
3557  pB0.rotbst(rot2to2frameInc);
3558  double sHat = (p1 + p2).m2Calc();
3559  double tHat = (p1 - pDaughter).m2Calc();
3560  double uHat = (p1 - pB0).m2Calc();
3561 
3562  // Calculate the weak t-channel correction.
3563  double m2R1 = 1. + pSister.m2Calc() / m2;
3564  double wt = 4. * sHat / (pMother + pB).m2Calc() * pT2 * ( 1. - z * m2R1)
3565  / (1. + pow2(z * m2R1)) / (1.-z);
3566  if (MEtype == 201 || MEtype == 206)
3567  wt *= weakShowerMEs.getMEqg2qgZ(pMother, pB, p2, pSister, p1)
3568  / weakShowerMEs.getMEqg2qg(sHat, tHat, uHat);
3569  else if (MEtype == 202 || MEtype == 207)
3570  wt *= weakShowerMEs.getMEqq2qqZ(pMother, pB, pSister, p2, p1)
3571  / weakShowerMEs.getMEqq2qq(sHat, tHat, uHat, true);
3572  else if (MEtype == 203 || MEtype == 208)
3573  wt *= weakShowerMEs.getMEqq2qqZ(pMother, pB, pSister, p2, p1)
3574  / weakShowerMEs.getMEqq2qq(sHat, tHat, uHat, false);
3575 
3576  // Split of ME into an ISR part and FSR part.
3577  wt *= (pSister + p1).m2Calc() / ( (pSister + p1).m2Calc()
3578  + abs((-pMother + pSister).m2Calc()) );
3579 
3580  // Remove the addition weight that was used to get an overestimate.
3581  wt /= calcMEmax(MEtype, 0, 0);
3582 
3583  return wt;
3584 }
3585 
3586 //--------------------------------------------------------------------------
3587 
3588 // Find coefficient of azimuthal asymmetry from gluon polarization.
3589 
3590 void SpaceShower::findAsymPol( Event& event, SpaceDipoleEnd* dip) {
3591 
3592  // Default is no asymmetry. Only gluons are studied.
3593  dip->iFinPol = 0;
3594  dip->asymPol = 0.;
3595  int iRad = dip->iRadiator;
3596  if (!doPhiPolAsym || dip->idDaughter != 21) return;
3597 
3598  // At least two particles in final state, whereof at least one coloured.
3599  int systemSizeOut = partonSystemsPtr->sizeOut( iSysSel);
3600  if (systemSizeOut < 2) return;
3601  bool foundColOut = false;
3602  for (int ii = 0; ii < systemSizeOut; ++ii) {
3603  int i = partonSystemsPtr->getOut( iSysSel, ii);
3604  if (event[i].col() != 0 || event[i].acol() != 0) foundColOut = true;
3605  }
3606  if (!foundColOut) return;
3607 
3608  // Check if granddaughter in final state of hard scattering.
3609  // (May need to trace across carbon copies to find granddaughters.)
3610  // If so, at most accept 2 -> 2 scatterings with gg or qq in final state.
3611  int iGrandD1 = event[iRad].daughter1();
3612  int iGrandD2 = event[iRad].daughter2();
3613  bool traceCopy = false;
3614  do {
3615  traceCopy = false;
3616  if (iGrandD1 > 0 && iGrandD2 == iGrandD1) {
3617  iGrandD1 = event[iGrandD2].daughter1();
3618  iGrandD2 = event[iGrandD2].daughter2();
3619  traceCopy = true;
3620  }
3621  } while (traceCopy);
3622  int statusGrandD1 = event[ iGrandD1 ].statusAbs();
3623  bool isHardProc = (statusGrandD1 == 23 || statusGrandD1 == 33);
3624  if (isHardProc) {
3625  if (!doPhiPolAsymHard) return;
3626  if (iGrandD2 != iGrandD1 + 1) return;
3627  if (event[iGrandD1].isGluon() && event[iGrandD2].isGluon());
3628  else if (event[iGrandD1].isQuark() && event[iGrandD2].isQuark());
3629  else return;
3630  }
3631  dip->iFinPol = iGrandD1;
3632 
3633  // Coefficient from gluon production.
3634  if (dip->idMother == 21) dip->asymPol = pow2( (1. - dip->z)
3635  / (1. - dip->z * (1. - dip->z) ) );
3636  else dip->asymPol = 2. * (1. - dip->z) / (1. + pow2(1. - dip->z) );
3637 
3638  // Coefficients from gluon decay. Put z = 1/2 for hard process.
3639  double zDau = (isHardProc) ? 0.5 : dip->zOld;
3640  if (event[iGrandD1].isGluon()) dip->asymPol *= pow2( zDau * (1. - zDau)
3641  / (1. - zDau * (1. - zDau) ) );
3642  else dip->asymPol *= -2. * zDau *( 1. - zDau )
3643  / (1. - 2. * zDau * (1. - zDau) );
3644 
3645 }
3646 
3647 //--------------------------------------------------------------------------
3648 
3649 // Find a possible colour partner in the case of dipole recoil.
3650 
3651 int SpaceShower::findColPartner(Event& event, int iSideA, int iSideB,
3652  int iSystem) {
3653 
3654  int iColPartner = 0;
3655  int colSideA = event[iSideA].col();
3656  int acolSideA = event[iSideA].acol();
3657 
3658  // Check if the parton on the other side is a colour partner.
3659  if ( (colSideA != 0 && event[iSideB].acol() == colSideA)
3660  || (acolSideA != 0 && event[iSideB].col() == acolSideA) ) {
3661 
3662  // Another possible colour partner among the outgoing partons
3663  // in the case of a gluon.
3664  if (event[iSideA].id() == 21)
3665  for (int i = 0; i < partonSystemsPtr->sizeOut(iSystem); ++i) {
3666  int iOut = partonSystemsPtr->getOut(iSystem, i);
3667  if ( event[iOut].col() == colSideA
3668  || event[iOut].acol() == acolSideA )
3669  // 50% for II and 50% for IF.
3670  if (rndmPtr->flat() < 0.5) iColPartner = iOut;
3671  }
3672 
3673  // Otherwise, check within the set of outgoing partons.
3674  } else if (colSideA != 0 || acolSideA != 0) {
3675  for (int i = 0; i < partonSystemsPtr->sizeOut(iSystem); ++ i) {
3676  int iOut = partonSystemsPtr->getOut(iSystem, i);
3677  if ( (colSideA != 0 && event[iOut].col() == colSideA)
3678  || (acolSideA != 0 && event[iOut].acol() == acolSideA) ) {
3679  if (iColPartner == 0) iColPartner = iOut;
3680  // 50% for each IF in the case of a gluon.
3681  else if (rndmPtr->flat() < 0.5) iColPartner = iOut;
3682  }
3683  }
3684  }
3685  return iColPartner;
3686 
3687 }
3688 
3689 //--------------------------------------------------------------------------
3690 
3691 // Remove weak dipoles if FSR already emitted a W/Z
3692 // and only a single weak emission is permited.
3693 // Update colour partner in case of dipole recoil.
3694 
3695 void SpaceShower::update(int iSys, Event& event, bool hasWeakRad) {
3696 
3697  if (hasWeakRad && singleWeakEmission)
3698  for (int i = 0; i < int(dipEnd.size()); i++)
3699  if (dipEnd[i].weakType != 0) dipEnd[i].weakType = 0;
3700  if (hasWeakRad) hasWeaklyRadiated = true;
3701 
3702  // Update colour partner in case of dipole recoil.
3703  if (doDipoleRecoil)
3704  for (int i = 0; i < int(dipEnd.size()); i++)
3705  if (dipEnd[i].system == iSys) {
3706  dipEnd[i].iColPartner = findColPartner(event, dipEnd[i].iRadiator,
3707  dipEnd[i].iRecoiler, iSys);
3708  dipEnd[i].idColPartner = (dipEnd[i].iColPartner != 0)
3709  ? event[dipEnd[i].iColPartner].id() : 0;
3710  }
3711 
3712 }
3713 
3714 //-------------------------------------------------------------------------
3715 
3716 // Print the list of dipoles.
3717 
3718 void SpaceShower::list() const {
3719 
3720  // Header.
3721  cout << "\n -------- PYTHIA SpaceShower Dipole Listing -------------- \n"
3722  << "\n i syst side rad rec pTmax col chg ME rec \n"
3723  << fixed << setprecision(3);
3724 
3725  // Loop over dipole list and print it.
3726  for (int i = 0; i < int(dipEnd.size()); ++i)
3727  cout << setw(5) << i << setw(6) << dipEnd[i].system
3728  << setw(6) << dipEnd[i].side << setw(6) << dipEnd[i].iRadiator
3729  << setw(6) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax
3730  << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType
3731  << setw(5) << dipEnd[i].MEtype << setw(4)
3732  << dipEnd[i].normalRecoil << "\n";
3733 
3734  // Done.
3735  cout << "\n -------- End PYTHIA SpaceShower Dipole Listing ----------"
3736  << endl;
3737 
3738 }
3739 
3740 //==========================================================================
3741 
3742 } // end namespace Pythia8
Definition: beam.h:43
Definition: AgUStep.h:26