1 | // TimeShower.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
5 | |
---|
6 | // Function definitions (not found in the header) for the TimeShower class. |
---|
7 | |
---|
8 | #include "TimeShower.h" |
---|
9 | |
---|
10 | namespace Pythia8 { |
---|
11 | |
---|
12 | //========================================================================== |
---|
13 | |
---|
14 | // The TimeShower class. |
---|
15 | |
---|
16 | //-------------------------------------------------------------------------- |
---|
17 | |
---|
18 | // Constants: could be changed here if desired, but normally should not. |
---|
19 | // These are of technical nature, as described for each. |
---|
20 | |
---|
21 | // For small x approximate 1 - sqrt(1 - x) by x/2. |
---|
22 | const double TimeShower::SIMPLIFYROOT = 1e-8; |
---|
23 | |
---|
24 | // Do not allow x too close to 0 or 1 in matrix element expressions. |
---|
25 | // Warning: cuts into phase space for E_CM > 2 * pTmin * sqrt(1/XMARGIN), |
---|
26 | // i.e. will become problem roughly for E_CM > 10^6 GeV. |
---|
27 | const double TimeShower::XMARGIN = 1e-12; |
---|
28 | const double TimeShower::XMARGINCOMB = 1e-4; |
---|
29 | |
---|
30 | // Lower limit on PDF value in order to avoid division by zero. |
---|
31 | const double TimeShower::TINYPDF = 1e-10; |
---|
32 | |
---|
33 | // Big starting value in search for smallest invariant-mass pair. |
---|
34 | const double TimeShower::LARGEM2 = 1e20; |
---|
35 | |
---|
36 | // In g -> q qbar or gamma -> f fbar require m2_pair > this * m2_q/f. |
---|
37 | const double TimeShower::THRESHM2 = 4.004; |
---|
38 | |
---|
39 | // Never pick pT so low that alphaS is evaluated too close to Lambda_3. |
---|
40 | const double TimeShower::LAMBDA3MARGIN = 1.1; |
---|
41 | |
---|
42 | // Rescatter: rescattering + ISR + FSR + primordial kT can lead to |
---|
43 | // systems not locally conserving momentum. |
---|
44 | // Fix up momentum in intermediate systems with rescattering |
---|
45 | const bool TimeShower::FIXRESCATTER = true; |
---|
46 | // Veto negative energies when using FIXRESCATTER option. |
---|
47 | const bool TimeShower::VETONEGENERGY = false; |
---|
48 | // Do not allow too large time- or spacelike virtualities in fixing-up. |
---|
49 | const double TimeShower::MAXVIRTUALITYFRACTION = 0.5; |
---|
50 | // Do not allow too large negative spacelike energy in system rest frame. |
---|
51 | const double TimeShower::MAXNEGENERGYFRACTION = 0.7; |
---|
52 | |
---|
53 | //-------------------------------------------------------------------------- |
---|
54 | |
---|
55 | // Initialize alphaStrong, alphaEM and related pTmin parameters. |
---|
56 | |
---|
57 | void TimeShower::init( BeamParticle* beamAPtrIn, |
---|
58 | BeamParticle* beamBPtrIn) { |
---|
59 | |
---|
60 | // Store input pointers for future use. |
---|
61 | beamAPtr = beamAPtrIn; |
---|
62 | beamBPtr = beamBPtrIn; |
---|
63 | |
---|
64 | // Main flags. |
---|
65 | doQCDshower = settingsPtr->flag("TimeShower:QCDshower"); |
---|
66 | doQEDshowerByQ = settingsPtr->flag("TimeShower:QEDshowerByQ"); |
---|
67 | doQEDshowerByL = settingsPtr->flag("TimeShower:QEDshowerByL"); |
---|
68 | doQEDshowerByGamma = settingsPtr->flag("TimeShower:QEDshowerByGamma"); |
---|
69 | doMEcorrections = settingsPtr->flag("TimeShower:MEcorrections"); |
---|
70 | doMEafterFirst = settingsPtr->flag("TimeShower:MEafterFirst"); |
---|
71 | doPhiPolAsym = settingsPtr->flag("TimeShower:phiPolAsym"); |
---|
72 | doInterleave = settingsPtr->flag("TimeShower:interleave"); |
---|
73 | allowBeamRecoil = settingsPtr->flag("TimeShower:allowBeamRecoil"); |
---|
74 | dampenBeamRecoil = settingsPtr->flag("TimeShower:dampenBeamRecoil"); |
---|
75 | recoilToColoured = settingsPtr->flag("TimeShower:recoilToColoured"); |
---|
76 | |
---|
77 | // Matching in pT of hard interaction or MPI to shower evolution. |
---|
78 | pTmaxMatch = settingsPtr->mode("TimeShower:pTmaxMatch"); |
---|
79 | pTdampMatch = settingsPtr->mode("TimeShower:pTdampMatch"); |
---|
80 | pTmaxFudge = settingsPtr->parm("TimeShower:pTmaxFudge"); |
---|
81 | pTmaxFudgeMPI = settingsPtr->parm("TimeShower:pTmaxFudgeMPI"); |
---|
82 | pTdampFudge = settingsPtr->parm("TimeShower:pTdampFudge"); |
---|
83 | |
---|
84 | // Charm and bottom mass thresholds. |
---|
85 | mc = particleDataPtr->m0(4); |
---|
86 | mb = particleDataPtr->m0(5); |
---|
87 | m2c = mc * mc; |
---|
88 | m2b = mb * mb; |
---|
89 | |
---|
90 | // Parameters of scale choices. |
---|
91 | renormMultFac = settingsPtr->parm("TimeShower:renormMultFac"); |
---|
92 | factorMultFac = settingsPtr->parm("TimeShower:factorMultFac"); |
---|
93 | |
---|
94 | // Parameters of alphaStrong generation. |
---|
95 | alphaSvalue = settingsPtr->parm("TimeShower:alphaSvalue"); |
---|
96 | alphaSorder = settingsPtr->mode("TimeShower:alphaSorder"); |
---|
97 | alphaS2pi = 0.5 * alphaSvalue / M_PI; |
---|
98 | |
---|
99 | // Initialize alphaStrong generation. |
---|
100 | alphaS.init( alphaSvalue, alphaSorder); |
---|
101 | |
---|
102 | // Lambda for 5, 4 and 3 flavours. |
---|
103 | Lambda3flav = alphaS.Lambda3(); |
---|
104 | Lambda4flav = alphaS.Lambda4(); |
---|
105 | Lambda5flav = alphaS.Lambda5(); |
---|
106 | Lambda5flav2 = pow2(Lambda5flav); |
---|
107 | Lambda4flav2 = pow2(Lambda4flav); |
---|
108 | Lambda3flav2 = pow2(Lambda3flav); |
---|
109 | |
---|
110 | // Parameters of QCD evolution. Warn if pTmin must be raised. |
---|
111 | nGluonToQuark = settingsPtr->mode("TimeShower:nGluonToQuark"); |
---|
112 | pTcolCutMin = settingsPtr->parm("TimeShower:pTmin"); |
---|
113 | if (pTcolCutMin > LAMBDA3MARGIN * Lambda3flav / sqrt(renormMultFac)) |
---|
114 | pTcolCut = pTcolCutMin; |
---|
115 | else { |
---|
116 | pTcolCut = LAMBDA3MARGIN * Lambda3flav / sqrt(renormMultFac); |
---|
117 | ostringstream newPTcolCut; |
---|
118 | newPTcolCut << fixed << setprecision(3) << pTcolCut; |
---|
119 | infoPtr->errorMsg("Warning in TimeShower::init: pTmin too low", |
---|
120 | ", raised to " + newPTcolCut.str() ); |
---|
121 | infoPtr->setTooLowPTmin(true); |
---|
122 | } |
---|
123 | pT2colCut = pow2(pTcolCut); |
---|
124 | |
---|
125 | // Parameters of alphaEM generation . |
---|
126 | alphaEMorder = settingsPtr->mode("TimeShower:alphaEMorder"); |
---|
127 | |
---|
128 | // Initialize alphaEM generation. |
---|
129 | alphaEM.init( alphaEMorder, settingsPtr); |
---|
130 | |
---|
131 | // Parameters of QED evolution. |
---|
132 | nGammaToQuark = settingsPtr->mode("TimeShower:nGammaToQuark"); |
---|
133 | nGammaToLepton = settingsPtr->mode("TimeShower:nGammaToLepton"); |
---|
134 | pTchgQCut = settingsPtr->parm("TimeShower:pTminChgQ"); |
---|
135 | pT2chgQCut = pow2(pTchgQCut); |
---|
136 | pTchgLCut = settingsPtr->parm("TimeShower:pTminChgL"); |
---|
137 | pT2chgLCut = pow2(pTchgLCut); |
---|
138 | mMaxGamma = settingsPtr->parm("TimeShower:mMaxGamma"); |
---|
139 | m2MaxGamma = pow2(mMaxGamma); |
---|
140 | |
---|
141 | // Consisteny check for gamma -> f fbar variables. |
---|
142 | if (nGammaToQuark <= 0 && nGammaToLepton <= 0) doQEDshowerByGamma = false; |
---|
143 | |
---|
144 | // Possibility of a global recoil stategy, e.g. for MC@NLO. |
---|
145 | globalRecoil = settingsPtr->flag("TimeShower:globalRecoil"); |
---|
146 | nMaxGlobalRecoil = settingsPtr->mode("TimeShower:nMaxGlobalRecoil"); |
---|
147 | |
---|
148 | // Fraction and colour factor of gluon emission off onium octat state. |
---|
149 | octetOniumFraction = settingsPtr->parm("TimeShower:octetOniumFraction"); |
---|
150 | octetOniumColFac = settingsPtr->parm("TimeShower:octetOniumColFac"); |
---|
151 | |
---|
152 | // Z0 properties needed for gamma/Z0 mixing. |
---|
153 | mZ = particleDataPtr->m0(23); |
---|
154 | gammaZ = particleDataPtr->mWidth(23); |
---|
155 | thetaWRat = 1. / (16. * coupSMPtr->sin2thetaW() |
---|
156 | * coupSMPtr->cos2thetaW()); |
---|
157 | |
---|
158 | // May have to fix up recoils related to rescattering. |
---|
159 | allowRescatter = settingsPtr->flag("PartonLevel:MPI") |
---|
160 | && settingsPtr->flag("MultipartonInteractions:allowRescatter"); |
---|
161 | |
---|
162 | // Hidden Valley scenario with further shower activity. |
---|
163 | doHVshower = settingsPtr->flag("HiddenValley:FSR"); |
---|
164 | nCHV = settingsPtr->mode("HiddenValley:Ngauge"); |
---|
165 | alphaHVfix = settingsPtr->parm("HiddenValley:alphaFSR"); |
---|
166 | pThvCut = settingsPtr->parm("HiddenValley:pTminFSR"); |
---|
167 | pT2hvCut = pThvCut * pThvCut; |
---|
168 | CFHV = (nCHV == 1) ? 1. : (nCHV * nCHV - 1.)/(2. * nCHV); |
---|
169 | idHV = (nCHV == 1) ? 4900022 : 4900021; |
---|
170 | mHV = particleDataPtr->m0(idHV); |
---|
171 | brokenHVsym = (nCHV == 1 && mHV > 0.); |
---|
172 | |
---|
173 | // Possibility to allow user veto of emission step. |
---|
174 | canVetoEmission = (userHooksPtr != 0) |
---|
175 | ? userHooksPtr->canVetoFSREmission() : false; |
---|
176 | |
---|
177 | } |
---|
178 | |
---|
179 | //-------------------------------------------------------------------------- |
---|
180 | |
---|
181 | // Find whether to limit maximum scale of emissions. |
---|
182 | // Also allow for dampening at factorization or renormalization scale. |
---|
183 | |
---|
184 | bool TimeShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) { |
---|
185 | |
---|
186 | // Find whether to limit pT. Begin by user-set cases. |
---|
187 | bool dopTlimit = false; |
---|
188 | if (pTmaxMatch == 1) dopTlimit = true; |
---|
189 | else if (pTmaxMatch == 2) dopTlimit = false; |
---|
190 | |
---|
191 | // Look if any quark (u, d, s, c, b), gluon or photon in final state. |
---|
192 | else { |
---|
193 | for (int i = 5; i < event.size(); ++i) |
---|
194 | if (event[i].status() != -21) { |
---|
195 | int idAbs = event[i].idAbs(); |
---|
196 | if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit = true; |
---|
197 | } |
---|
198 | } |
---|
199 | |
---|
200 | // Dampening at factorization or renormalization scale. |
---|
201 | dopTdamp = false; |
---|
202 | pT2damp = 0.; |
---|
203 | if ( !dopTlimit && (pTdampMatch == 1 || pTdampMatch == 2) ) { |
---|
204 | dopTdamp = true; |
---|
205 | pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren); |
---|
206 | } |
---|
207 | |
---|
208 | // Done. |
---|
209 | return dopTlimit; |
---|
210 | |
---|
211 | } |
---|
212 | |
---|
213 | //-------------------------------------------------------------------------- |
---|
214 | |
---|
215 | // Top-level routine to do a full time-like shower in resonance decay. |
---|
216 | |
---|
217 | int TimeShower::shower( int iBeg, int iEnd, Event& event, double pTmax, |
---|
218 | int nBranchMax) { |
---|
219 | |
---|
220 | // Add new system, automatically with two empty beam slots. |
---|
221 | int iSys = partonSystemsPtr->addSys(); |
---|
222 | |
---|
223 | // Loop over allowed range to find all final-state particles. |
---|
224 | Vec4 pSum; |
---|
225 | for (int i = iBeg; i <= iEnd; ++i) if (event[i].isFinal()) { |
---|
226 | partonSystemsPtr->addOut( iSys, i); |
---|
227 | pSum += event[i].p(); |
---|
228 | } |
---|
229 | partonSystemsPtr->setSHat( iSys, pSum.m2Calc() ); |
---|
230 | |
---|
231 | // Let prepare routine do the setup. |
---|
232 | prepare( iSys, event, true); |
---|
233 | |
---|
234 | // Begin evolution down in pT from hard pT scale. |
---|
235 | int nBranch = 0; |
---|
236 | pTLastBranch = 0.; |
---|
237 | do { |
---|
238 | double pTtimes = pTnext( event, pTmax, 0.); |
---|
239 | |
---|
240 | // Do a final-state emission (if allowed). |
---|
241 | if (pTtimes > 0.) { |
---|
242 | if (branch( event)) { |
---|
243 | ++nBranch; |
---|
244 | pTLastBranch = pTtimes; |
---|
245 | } |
---|
246 | pTmax = pTtimes; |
---|
247 | } |
---|
248 | |
---|
249 | // Keep on evolving until nothing is left to be done. |
---|
250 | else pTmax = 0.; |
---|
251 | } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax)); |
---|
252 | |
---|
253 | // Return number of emissions that were performed. |
---|
254 | return nBranch; |
---|
255 | |
---|
256 | } |
---|
257 | |
---|
258 | //-------------------------------------------------------------------------- |
---|
259 | |
---|
260 | // Prepare system for evolution; identify ME. |
---|
261 | |
---|
262 | void TimeShower::prepare( int iSys, Event& event, bool limitPTmaxIn) { |
---|
263 | |
---|
264 | // Reset dipole-ends list for first interaction and for resonance decays. |
---|
265 | int iInA = partonSystemsPtr->getInA(iSys); |
---|
266 | int iInB = partonSystemsPtr->getInB(iSys); |
---|
267 | if (iSys == 0 || iInA == 0) dipEnd.resize(0); |
---|
268 | int dipEndSizeBeg = dipEnd.size(); |
---|
269 | |
---|
270 | // No dipoles for 2 -> 1 processes. |
---|
271 | if (partonSystemsPtr->sizeOut(iSys) < 2) return; |
---|
272 | |
---|
273 | // Loop through final state of system to find possible dipole ends. |
---|
274 | for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { |
---|
275 | int iRad = partonSystemsPtr->getOut( iSys, i); |
---|
276 | if (event[iRad].isFinal() && event[iRad].scale() > 0.) { |
---|
277 | |
---|
278 | // Identify colour octet onium state. Check whether QCD shower allowed. |
---|
279 | int idRad = event[iRad].id(); |
---|
280 | int idRadAbs = abs(idRad); |
---|
281 | bool isOctetOnium |
---|
282 | = ( idRad == 9900441 || idRad == 9900443 || idRad == 9910441 |
---|
283 | || idRad == 9900551 || idRad == 9900553 || idRad == 9910551 ); |
---|
284 | bool doQCD = doQCDshower; |
---|
285 | if (doQCD && isOctetOnium) |
---|
286 | doQCD = (rndmPtr->flat() < octetOniumFraction); |
---|
287 | |
---|
288 | // Find dipole end formed by colour index. |
---|
289 | int colTag = event[iRad].col(); |
---|
290 | if (doQCD && colTag > 0) setupQCDdip( iSys, i, colTag, 1, event, |
---|
291 | isOctetOnium, limitPTmaxIn); |
---|
292 | |
---|
293 | // Find dipole end formed by anticolour index. |
---|
294 | int acolTag = event[iRad].acol(); |
---|
295 | if (doQCD && acolTag > 0) setupQCDdip( iSys, i, acolTag, -1, event, |
---|
296 | isOctetOnium, limitPTmaxIn); |
---|
297 | |
---|
298 | // Find "charge-dipole" and "photon-dipole" ends. |
---|
299 | int chgType = event[iRad].chargeType(); |
---|
300 | bool doChgDip = (chgType != 0) |
---|
301 | && ( ( doQEDshowerByQ && event[iRad].isQuark() ) |
---|
302 | || ( doQEDshowerByL && event[iRad].isLepton() ) ); |
---|
303 | int gamType = (idRad == 22) ? 1 : 0; |
---|
304 | bool doGamDip = (gamType == 1) && doQEDshowerByGamma; |
---|
305 | if (doChgDip || doGamDip) setupQEDdip( iSys, i, chgType, gamType, |
---|
306 | event, limitPTmaxIn); |
---|
307 | |
---|
308 | // Find Hidden Valley dipole ends. |
---|
309 | bool isHVrad = (idRadAbs > 4900000 && idRadAbs < 4900007) |
---|
310 | || (idRadAbs > 4900010 && idRadAbs < 4900017) |
---|
311 | || idRadAbs == 4900101; |
---|
312 | if (doHVshower && isHVrad) setupHVdip( iSys, i, event, limitPTmaxIn); |
---|
313 | |
---|
314 | // End loop over system final state. Have now found the dipole ends. |
---|
315 | } |
---|
316 | } |
---|
317 | |
---|
318 | // Loop through dipole ends to find matrix element corrections. |
---|
319 | for (int iDip = dipEndSizeBeg; iDip < int(dipEnd.size()); ++iDip) |
---|
320 | findMEtype( event, dipEnd[iDip]); |
---|
321 | |
---|
322 | // Update dipole list after a multiparton interactions rescattering. |
---|
323 | if (iSys > 0 && ( (iInA > 0 && event[iInA].status() == -34) |
---|
324 | || (iInB > 0 && event[iInB].status() == -34) ) ) |
---|
325 | rescatterUpdate( iSys, event); |
---|
326 | |
---|
327 | } |
---|
328 | |
---|
329 | //-------------------------------------------------------------------------- |
---|
330 | |
---|
331 | // Update dipole list after a multiparton interactions rescattering. |
---|
332 | void TimeShower::rescatterUpdate( int iSys, Event& event) { |
---|
333 | |
---|
334 | // Loop over two incoming partons in system; find their rescattering mother. |
---|
335 | // (iOut is outgoing from old system = incoming iIn of rescattering system.) |
---|
336 | for (int iResc = 0; iResc < 2; ++iResc) { |
---|
337 | int iIn = (iResc == 0) ? partonSystemsPtr->getInA(iSys) |
---|
338 | : partonSystemsPtr->getInB(iSys); |
---|
339 | if (iIn == 0 || event[iIn].status() != -34) continue; |
---|
340 | int iOut = event[iIn].mother1(); |
---|
341 | |
---|
342 | // Loop over all dipoles. |
---|
343 | int dipEndSize = dipEnd.size(); |
---|
344 | for (int iDip = 0; iDip < dipEndSize; ++iDip) { |
---|
345 | TimeDipoleEnd& dipNow = dipEnd[iDip]; |
---|
346 | |
---|
347 | // Kill dipoles where rescattered parton is radiator. |
---|
348 | if (dipNow.iRadiator == iOut) { |
---|
349 | dipNow.colType = 0; |
---|
350 | dipNow.chgType = 0; |
---|
351 | dipNow.gamType = 0; |
---|
352 | continue; |
---|
353 | } |
---|
354 | // No matrix element for dipoles between scatterings. |
---|
355 | if (dipNow.iMEpartner == iOut) { |
---|
356 | dipNow.MEtype = 0; |
---|
357 | dipNow.iMEpartner = -1; |
---|
358 | } |
---|
359 | |
---|
360 | // Update dipoles where outgoing rescattered parton is recoiler. |
---|
361 | if (dipNow.iRecoiler == iOut) { |
---|
362 | int iRad = dipNow.iRadiator; |
---|
363 | |
---|
364 | // Colour dipole: recoil in final state, initial state or new. |
---|
365 | if (dipNow.colType > 0) { |
---|
366 | int colTag = event[iRad].col(); |
---|
367 | bool done = false; |
---|
368 | for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { |
---|
369 | int iRecNow = partonSystemsPtr->getOut( iSys, i); |
---|
370 | if (event[iRecNow].acol() == colTag) { |
---|
371 | dipNow.iRecoiler = iRecNow; |
---|
372 | dipNow.systemRec = iSys; |
---|
373 | dipNow.MEtype = 0; |
---|
374 | done = true; |
---|
375 | break; |
---|
376 | } |
---|
377 | } |
---|
378 | if (!done) { |
---|
379 | int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys) |
---|
380 | : partonSystemsPtr->getInA(iSys); |
---|
381 | if (event[iIn2].col() == colTag) { |
---|
382 | dipNow.iRecoiler = iIn2; |
---|
383 | dipNow.systemRec = iSys; |
---|
384 | dipNow.MEtype = 0; |
---|
385 | int isrType = event[iIn2].mother1(); |
---|
386 | // This line in case mother is a rescattered parton. |
---|
387 | while (isrType > 2 + beamOffset) |
---|
388 | isrType = event[isrType].mother1(); |
---|
389 | if (isrType > 2) isrType -= beamOffset; |
---|
390 | dipNow.isrType = isrType; |
---|
391 | done = true; |
---|
392 | } |
---|
393 | } |
---|
394 | // If above options failed, then create new dipole. |
---|
395 | if (!done) { |
---|
396 | int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad); |
---|
397 | if (iRadNow != -1) |
---|
398 | setupQCDdip(dipNow.system, iRadNow, event[iRad].col(), 1, |
---|
399 | event, dipNow.isOctetOnium, true); |
---|
400 | else |
---|
401 | infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " |
---|
402 | "failed to locate radiator in system"); |
---|
403 | |
---|
404 | dipNow.colType = 0; |
---|
405 | dipNow.chgType = 0; |
---|
406 | dipNow.gamType = 0; |
---|
407 | |
---|
408 | infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " |
---|
409 | "failed to locate new recoiling colour partner"); |
---|
410 | } |
---|
411 | |
---|
412 | // Anticolour dipole: recoil in final state, initial state or new. |
---|
413 | } else if (dipNow.colType < 0) { |
---|
414 | int acolTag = event[iRad].acol(); |
---|
415 | bool done = false; |
---|
416 | for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { |
---|
417 | int iRecNow = partonSystemsPtr->getOut( iSys, i); |
---|
418 | if (event[iRecNow].col() == acolTag) { |
---|
419 | dipNow.iRecoiler = iRecNow; |
---|
420 | dipNow.systemRec = iSys; |
---|
421 | dipNow.MEtype = 0; |
---|
422 | done = true; |
---|
423 | break; |
---|
424 | } |
---|
425 | } |
---|
426 | if (!done) { |
---|
427 | int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys) |
---|
428 | : partonSystemsPtr->getInA(iSys); |
---|
429 | if (event[iIn2].acol() == acolTag) { |
---|
430 | dipNow.iRecoiler = iIn2; |
---|
431 | dipNow.systemRec = iSys; |
---|
432 | dipNow.MEtype = 0; |
---|
433 | int isrType = event[iIn2].mother1(); |
---|
434 | // This line in case mother is a rescattered parton. |
---|
435 | while (isrType > 2 + beamOffset) |
---|
436 | isrType = event[isrType].mother1(); |
---|
437 | if (isrType > 2) isrType -= beamOffset; |
---|
438 | dipNow.isrType = isrType; |
---|
439 | done = true; |
---|
440 | } |
---|
441 | } |
---|
442 | // If above options failed, then create new dipole. |
---|
443 | if (!done) { |
---|
444 | int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad); |
---|
445 | if (iRadNow != -1) |
---|
446 | setupQCDdip(dipNow.system, iRadNow, event[iRad].acol(), -1, |
---|
447 | event, dipNow.isOctetOnium, true); |
---|
448 | else |
---|
449 | infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " |
---|
450 | "failed to locate radiator in system"); |
---|
451 | |
---|
452 | dipNow.colType = 0; |
---|
453 | dipNow.chgType = 0; |
---|
454 | dipNow.gamType = 0; |
---|
455 | |
---|
456 | infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " |
---|
457 | "failed to locate new recoiling colour partner"); |
---|
458 | } |
---|
459 | |
---|
460 | // Charge or photon dipoles: same flavour in final or initial state. |
---|
461 | } else if (dipNow.chgType != 0 || dipNow.gamType != 0) { |
---|
462 | int idTag = event[dipNow.iRecoiler].id(); |
---|
463 | bool done = false; |
---|
464 | for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { |
---|
465 | int iRecNow = partonSystemsPtr->getOut( iSys, i); |
---|
466 | if (event[iRecNow].id() == idTag) { |
---|
467 | dipNow.iRecoiler = iRecNow; |
---|
468 | dipNow.systemRec = iSys; |
---|
469 | dipNow.MEtype = 0; |
---|
470 | done = true; |
---|
471 | break; |
---|
472 | } |
---|
473 | } |
---|
474 | if (!done) { |
---|
475 | int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys) |
---|
476 | : partonSystemsPtr->getInA(iSys); |
---|
477 | if (event[iIn2].id() == -idTag) { |
---|
478 | dipNow.iRecoiler = iIn2; |
---|
479 | dipNow.systemRec = iSys; |
---|
480 | dipNow.MEtype = 0; |
---|
481 | int isrType = event[iIn2].mother1(); |
---|
482 | // This line in case mother is a rescattered parton. |
---|
483 | while (isrType > 2 + beamOffset) |
---|
484 | isrType = event[isrType].mother1(); |
---|
485 | if (isrType > 2) isrType -= beamOffset; |
---|
486 | dipNow.isrType = isrType; |
---|
487 | done = true; |
---|
488 | } |
---|
489 | } |
---|
490 | // If above options failed, then create new dipole |
---|
491 | if (!done) { |
---|
492 | int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad); |
---|
493 | if (iRadNow != -1) |
---|
494 | setupQEDdip(dipNow.system, iRadNow, dipNow.chgType, |
---|
495 | dipNow.gamType, event, true); |
---|
496 | else |
---|
497 | infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " |
---|
498 | "failed to locate radiator in system"); |
---|
499 | |
---|
500 | dipNow.colType = 0; |
---|
501 | dipNow.chgType = 0; |
---|
502 | dipNow.gamType = 0; |
---|
503 | |
---|
504 | infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " |
---|
505 | "failed to locate new recoiling charge partner"); |
---|
506 | } |
---|
507 | } |
---|
508 | } |
---|
509 | |
---|
510 | // End of loop over dipoles and two incoming sides. |
---|
511 | } |
---|
512 | } |
---|
513 | |
---|
514 | } |
---|
515 | |
---|
516 | //-------------------------------------------------------------------------- |
---|
517 | |
---|
518 | // Update dipole list after each ISR emission (so not used for resonances). |
---|
519 | |
---|
520 | void TimeShower::update( int iSys, Event& event) { |
---|
521 | |
---|
522 | // Start list of rescatterers that gave further changed systems in ISR. |
---|
523 | vector<int> iRescatterer; |
---|
524 | |
---|
525 | // Find new and old positions of incoming partons in the system. |
---|
526 | vector<int> iNew, iOld; |
---|
527 | iNew.push_back( partonSystemsPtr->getInA(iSys) ); |
---|
528 | iOld.push_back( event[iNew[0]].daughter2() ); |
---|
529 | iNew.push_back( partonSystemsPtr->getInB(iSys) ); |
---|
530 | iOld.push_back( event[iNew[1]].daughter2() ); |
---|
531 | |
---|
532 | // Ditto for outgoing partons, except the newly created one. |
---|
533 | int sizeOut = partonSystemsPtr->sizeOut(iSys) - 1; |
---|
534 | for (int i = 0; i < sizeOut; ++i) { |
---|
535 | int iNow = partonSystemsPtr->getOut(iSys, i); |
---|
536 | iNew.push_back( iNow ); |
---|
537 | iOld.push_back( event[iNow].mother1() ); |
---|
538 | // Add non-final to list of rescatterers. |
---|
539 | if (!event[iNow].isFinal()) iRescatterer.push_back( iNow ); |
---|
540 | } |
---|
541 | int iNewNew = partonSystemsPtr->getOut(iSys, sizeOut); |
---|
542 | |
---|
543 | // Swap beams to let 0 be side on which branching occured. |
---|
544 | if (event[iNew[0]].status() != -41) { |
---|
545 | swap( iNew[0], iNew[1]); |
---|
546 | swap( iOld[0], iOld[1]); |
---|
547 | } |
---|
548 | |
---|
549 | // Loop over all dipole ends belonging to the system |
---|
550 | // or to the recoil system, if different. |
---|
551 | for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) |
---|
552 | if (dipEnd[iDip].system == iSys || dipEnd[iDip].systemRec == iSys) { |
---|
553 | TimeDipoleEnd& dipNow = dipEnd[iDip]; |
---|
554 | |
---|
555 | // Replace radiator (always in final state so simple). |
---|
556 | for (int i = 2; i < 2 + sizeOut; ++i) |
---|
557 | if (dipNow.iRadiator == iOld[i]) { |
---|
558 | dipNow.iRadiator = iNew[i]; |
---|
559 | break; |
---|
560 | } |
---|
561 | |
---|
562 | // Replace ME partner (always in final state, if exists, so simple). |
---|
563 | for (int i = 2; i < 2 + sizeOut; ++i) |
---|
564 | if (dipNow.iMEpartner == iOld[i]) { |
---|
565 | dipNow.iMEpartner = iNew[i]; |
---|
566 | break; |
---|
567 | } |
---|
568 | |
---|
569 | // Recoiler: by default pick old one, only moved. Note excluded beam. |
---|
570 | int iRec = 0; |
---|
571 | if (dipNow.systemRec == iSys) { |
---|
572 | for (int i = 1; i < 2 + sizeOut; ++i) |
---|
573 | if (dipNow.iRecoiler == iOld[i]) { |
---|
574 | iRec = iNew[i]; |
---|
575 | break; |
---|
576 | } |
---|
577 | |
---|
578 | // QCD recoiler: check if colour hooks up with new final parton. |
---|
579 | if ( dipNow.colType > 0 |
---|
580 | && event[dipNow.iRadiator].col() == event[iNewNew].acol() ) { |
---|
581 | iRec = iNewNew; |
---|
582 | dipNow.isrType = 0; |
---|
583 | } |
---|
584 | if ( dipNow.colType < 0 |
---|
585 | && event[dipNow.iRadiator].acol() == event[iNewNew].col() ) { |
---|
586 | iRec = iNewNew; |
---|
587 | dipNow.isrType = 0; |
---|
588 | } |
---|
589 | |
---|
590 | // QCD recoiler: check if colour hooks up with new beam parton. |
---|
591 | if ( iRec == 0 && dipNow.colType > 0 |
---|
592 | && event[dipNow.iRadiator].col() == event[iNew[0]].col() ) |
---|
593 | iRec = iNew[0]; |
---|
594 | if ( iRec == 0 && dipNow.colType < 0 |
---|
595 | && event[dipNow.iRadiator].acol() == event[iNew[0]].acol() ) |
---|
596 | iRec = iNew[0]; |
---|
597 | |
---|
598 | // QED/photon recoiler: either to new particle or remains to beam. |
---|
599 | if ( iRec == 0 && (dipNow.chgType != 0 || dipNow.gamType != 0) ) { |
---|
600 | if ( event[iNew[0]].chargeType() == 0 ) { |
---|
601 | iRec = iNewNew; |
---|
602 | dipNow.isrType = 0; |
---|
603 | } else { |
---|
604 | iRec = iNew[0]; |
---|
605 | } |
---|
606 | } |
---|
607 | |
---|
608 | // Recoiler in another system: keep it as is. |
---|
609 | } else iRec = dipNow.iRecoiler; |
---|
610 | |
---|
611 | // Done. Kill dipole if failed to find new recoiler. |
---|
612 | dipNow.iRecoiler = iRec; |
---|
613 | if ( iRec == 0 && (dipNow.colType != 0 || dipNow.chgType != 0 |
---|
614 | || dipNow.gamType != 0) ) { |
---|
615 | dipNow.colType = 0; |
---|
616 | dipNow.chgType = 0; |
---|
617 | dipNow.gamType = 0; |
---|
618 | infoPtr->errorMsg("Error in TimeShower::update: " |
---|
619 | "failed to locate new recoiling partner"); |
---|
620 | } |
---|
621 | } |
---|
622 | |
---|
623 | // Find new dipole end formed by colour index. |
---|
624 | int colTag = event[iNewNew].col(); |
---|
625 | if (doQCDshower && colTag > 0) |
---|
626 | setupQCDdip( iSys, sizeOut, colTag, 1, event, false, true); |
---|
627 | |
---|
628 | // Find new dipole end formed by anticolour index. |
---|
629 | int acolTag = event[iNewNew].acol(); |
---|
630 | if (doQCDshower && acolTag > 0) |
---|
631 | setupQCDdip( iSys, sizeOut, acolTag, -1, event, false, true); |
---|
632 | |
---|
633 | // Find new "charge-dipole" and "photon-dipole" ends. |
---|
634 | int chgType = event[iNewNew].chargeType(); |
---|
635 | bool doChgDip = (chgType != 0) |
---|
636 | && ( ( doQEDshowerByQ && event[iNewNew].isQuark() ) |
---|
637 | || ( doQEDshowerByL && event[iNewNew].isLepton() ) ); |
---|
638 | int gamType = (event[iNewNew].id() == 22) ? 1 : 0; |
---|
639 | bool doGamDip = (gamType == 1) && doQEDshowerByGamma; |
---|
640 | if (doChgDip || doGamDip) |
---|
641 | setupQEDdip( iSys, sizeOut, chgType, gamType, event, true); |
---|
642 | |
---|
643 | // Start iterate over list of rescatterers - may be empty. |
---|
644 | int iRescNow = -1; |
---|
645 | while (++iRescNow < int(iRescatterer.size())) { |
---|
646 | |
---|
647 | // Identify systems that rescatterers belong to. |
---|
648 | int iOutNew = iRescatterer[iRescNow]; |
---|
649 | int iInNew = event[iOutNew].daughter1(); |
---|
650 | int iSysResc = partonSystemsPtr->getSystemOf(iInNew, true); |
---|
651 | |
---|
652 | // Find new and old positions of incoming partons in the system. |
---|
653 | iNew.resize(0); |
---|
654 | iOld.resize(0); |
---|
655 | iNew.push_back( partonSystemsPtr->getInA(iSysResc) ); |
---|
656 | iOld.push_back( event[iNew[0]].daughter1() ); |
---|
657 | iNew.push_back( partonSystemsPtr->getInB(iSysResc) ); |
---|
658 | iOld.push_back( event[iNew[1]].daughter1() ); |
---|
659 | |
---|
660 | // Ditto for outgoing partons. |
---|
661 | sizeOut = partonSystemsPtr->sizeOut(iSysResc); |
---|
662 | for (int i = 0; i < sizeOut; ++i) { |
---|
663 | int iNow = partonSystemsPtr->getOut(iSysResc, i); |
---|
664 | iNew.push_back( iNow ); |
---|
665 | iOld.push_back( event[iNow].mother1() ); |
---|
666 | // Add non-final to list of rescatterers. |
---|
667 | if (!event[iNow].isFinal()) iRescatterer.push_back( iNow ); |
---|
668 | } |
---|
669 | |
---|
670 | // Loop over all dipole ends belonging to the system |
---|
671 | // or to the recoil system, if different. |
---|
672 | for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) |
---|
673 | if (dipEnd[iDip].system == iSysResc |
---|
674 | || dipEnd[iDip].systemRec == iSysResc) { |
---|
675 | TimeDipoleEnd& dipNow = dipEnd[iDip]; |
---|
676 | |
---|
677 | // Replace radiator (always in final state so simple). |
---|
678 | for (int i = 2; i < 2 + sizeOut; ++i) |
---|
679 | if (dipNow.iRadiator == iOld[i]) { |
---|
680 | dipNow.iRadiator = iNew[i]; |
---|
681 | break; |
---|
682 | } |
---|
683 | |
---|
684 | // Replace ME partner (always in final state, if exists, so simple). |
---|
685 | for (int i = 2; i < 2 + sizeOut; ++i) |
---|
686 | if (dipNow.iMEpartner == iOld[i]) { |
---|
687 | dipNow.iMEpartner = iNew[i]; |
---|
688 | break; |
---|
689 | } |
---|
690 | |
---|
691 | // Replace recoiler. |
---|
692 | for (int i = 0; i < 2 + sizeOut; ++i) |
---|
693 | if (dipNow.iRecoiler == iOld[i]) { |
---|
694 | dipNow.iRecoiler = iNew[i]; |
---|
695 | break; |
---|
696 | } |
---|
697 | } |
---|
698 | |
---|
699 | // End iterate over list of rescatterers. |
---|
700 | } |
---|
701 | |
---|
702 | } |
---|
703 | |
---|
704 | //-------------------------------------------------------------------------- |
---|
705 | |
---|
706 | // Setup a dipole end for a QCD colour charge. |
---|
707 | |
---|
708 | void TimeShower::setupQCDdip( int iSys, int i, int colTag, int colSign, |
---|
709 | Event& event, bool isOctetOnium, bool limitPTmaxIn) { |
---|
710 | |
---|
711 | // Initial values. Find if allowed to hook up beams. |
---|
712 | int iRad = partonSystemsPtr->getOut(iSys, i); |
---|
713 | int iRec = 0; |
---|
714 | int sizeAllA = partonSystemsPtr->sizeAll(iSys); |
---|
715 | int sizeOut = partonSystemsPtr->sizeOut(iSys); |
---|
716 | int sizeAll = ( allowBeamRecoil ) ? sizeAllA : sizeOut; |
---|
717 | int sizeIn = sizeAll - sizeOut; |
---|
718 | int sizeInA = sizeAllA - sizeIn - sizeOut; |
---|
719 | int iOffset = i + sizeAllA - sizeOut; |
---|
720 | bool otherSystemRec = false; |
---|
721 | bool allowInitial = (partonSystemsPtr->hasInAB(iSys)) ? true : false; |
---|
722 | // PS dec 2010: possibility to allow for several recoilers and each with |
---|
723 | // flexible normalization |
---|
724 | bool isFlexible = false; |
---|
725 | double flexFactor = 1.0; |
---|
726 | vector<int> iRecVec(0); |
---|
727 | |
---|
728 | // Colour: other end by same index in beam or opposite in final state. |
---|
729 | // Exclude rescattered incoming and not final outgoing. |
---|
730 | if (colSign > 0) |
---|
731 | for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { |
---|
732 | int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); |
---|
733 | if ( ( j < sizeIn && event[iRecNow].col() == colTag |
---|
734 | && !event[iRecNow].isRescatteredIncoming() ) |
---|
735 | || ( j >= sizeIn && event[iRecNow].acol() == colTag |
---|
736 | && event[iRecNow].isFinal() ) ) { |
---|
737 | iRec = iRecNow; |
---|
738 | break; |
---|
739 | } |
---|
740 | } |
---|
741 | |
---|
742 | // Anticolour: other end by same index in beam or opposite in final state. |
---|
743 | // Exclude rescattered incoming and not final outgoing. |
---|
744 | if (colSign < 0) |
---|
745 | for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { |
---|
746 | int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); |
---|
747 | if ( ( j < sizeIn && event[iRecNow].acol() == colTag |
---|
748 | && !event[iRecNow].isRescatteredIncoming() ) |
---|
749 | || ( j >= sizeIn && event[iRecNow].col() == colTag |
---|
750 | && event[iRecNow].isFinal() ) ) { |
---|
751 | iRec = iRecNow; |
---|
752 | break; |
---|
753 | } |
---|
754 | } |
---|
755 | |
---|
756 | // Resonance decays (= no instate): |
---|
757 | // other end to nearest recoiler in same system final state, |
---|
758 | // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j). |
---|
759 | // (junction colours more involved, so keep track if junction colour) |
---|
760 | bool hasJunction = false; |
---|
761 | if (iRec == 0 && !allowInitial) { |
---|
762 | for (int iJun = 0; iJun < event.sizeJunction(); ++ iJun) { |
---|
763 | // For types 1&2, all legs in final state |
---|
764 | // For types 3&4, two legs in final state |
---|
765 | // For types 5&6, one leg in final state |
---|
766 | int iBeg = (event.kindJunction(iJun)-1)/2; |
---|
767 | for (int iLeg = iBeg; iLeg < 3; ++iLeg) |
---|
768 | if (event.endColJunction( iJun, iLeg) == colTag) hasJunction = true; |
---|
769 | } |
---|
770 | double ppMin = LARGEM2; |
---|
771 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
772 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
773 | if (!event[iRecNow].isFinal()) continue; |
---|
774 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
775 | - event[iRecNow].m() * event[iRad].m(); |
---|
776 | if (ppNow < ppMin) { |
---|
777 | iRec = iRecNow; |
---|
778 | ppMin = ppNow; |
---|
779 | } |
---|
780 | } |
---|
781 | } |
---|
782 | |
---|
783 | // If no success then look for matching (anti)colour anywhere in final state. |
---|
784 | if ( iRec == 0 || (!doInterleave && !event[iRec].isFinal()) ) { |
---|
785 | iRec = 0; |
---|
786 | for (int j = 0; j < event.size(); ++j) if (event[j].isFinal()) |
---|
787 | if ( (colSign > 0 && event[j].acol() == colTag) |
---|
788 | || (colSign < 0 && event[j].col() == colTag) ) { |
---|
789 | iRec = j; |
---|
790 | otherSystemRec = true; |
---|
791 | break; |
---|
792 | } |
---|
793 | |
---|
794 | // If no success then look for match to non-rescattered in initial state. |
---|
795 | if (iRec == 0 && allowInitial) { |
---|
796 | for (int iSysR = 0; iSysR < partonSystemsPtr->sizeSys(); ++iSysR) |
---|
797 | if (iSysR != iSys) { |
---|
798 | int j = partonSystemsPtr->getInA(iSysR); |
---|
799 | if (j > 0 && event[j].isRescatteredIncoming()) j = 0; |
---|
800 | if (j > 0 && ( (colSign > 0 && event[j].col() == colTag) |
---|
801 | || (colSign < 0 && event[j].acol() == colTag) ) ) { |
---|
802 | iRec = j; |
---|
803 | otherSystemRec = true; |
---|
804 | break; |
---|
805 | } |
---|
806 | j = partonSystemsPtr->getInB(iSysR); |
---|
807 | if (j > 0 && event[j].isRescatteredIncoming()) j = 0; |
---|
808 | if (j > 0 && ( (colSign > 0 && event[j].col() == colTag) |
---|
809 | || (colSign < 0 && event[j].acol() == colTag) ) ) { |
---|
810 | iRec = j; |
---|
811 | otherSystemRec = true; |
---|
812 | break; |
---|
813 | } |
---|
814 | } |
---|
815 | } |
---|
816 | } |
---|
817 | |
---|
818 | // Junctions (PS&ND dec 2010) |
---|
819 | // For types 1&2: all legs in final state |
---|
820 | // half-strength dipoles between all legs |
---|
821 | // For types 3&4, two legs in final state |
---|
822 | // full-strength dipole between final-state legs |
---|
823 | // For types 5&6, one leg in final state |
---|
824 | // no final-state dipole end |
---|
825 | |
---|
826 | if (hasJunction) { |
---|
827 | for (int iJun = 0; iJun < event.sizeJunction(); ++ iJun) { |
---|
828 | int kindJun = event.kindJunction(iJun); |
---|
829 | int iBeg = (kindJun-1)/2; |
---|
830 | for (int iLeg = iBeg; iLeg < 3; ++iLeg) { |
---|
831 | if (event.endColJunction( iJun, iLeg) == colTag) { |
---|
832 | // For types 5&6, no other leg to recoil against. Switch off if |
---|
833 | // no other particles at all, since radiation then handled by ISR. |
---|
834 | // Example: qq -> ~t* : no radiation off ~t* |
---|
835 | // Allow radiation + recoil if unconnected partners available |
---|
836 | // Example: qq -> ~t* -> tbar ~chi0 : allow radiation off tbar, |
---|
837 | // with ~chi0 as recoiler |
---|
838 | if (kindJun >= 5) { |
---|
839 | if (sizeOut == 1) return; |
---|
840 | else break; |
---|
841 | } |
---|
842 | // For junction types 3 & 4, span one full-strength dipole |
---|
843 | // (only look inside same decay system) |
---|
844 | else if (kindJun >= 3) { |
---|
845 | int iLegRec = 3-iLeg; |
---|
846 | int colTagRec = event.endColJunction( iJun, iLegRec); |
---|
847 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
848 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
849 | if (!event[iRecNow].isFinal()) continue; |
---|
850 | if ( (colSign > 0 && event[iRecNow].col() == colTagRec) |
---|
851 | || (colSign < 0 && event[iRecNow].acol() == colTagRec) ) { |
---|
852 | // Only accept if staying inside same system |
---|
853 | iRec = iRecNow; |
---|
854 | break; |
---|
855 | } |
---|
856 | } |
---|
857 | } |
---|
858 | // For junction types 1 & 2, span two half-strength dipoles |
---|
859 | // (only look inside same decay system) |
---|
860 | else { |
---|
861 | // Loop over two half-strength dipole connections |
---|
862 | for (int jLeg = 1; jLeg <= 2; jLeg++) { |
---|
863 | int iLegRec = (iLeg + jLeg) % 3; |
---|
864 | int colTagRec = event.endColJunction( iJun, iLegRec); |
---|
865 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
866 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
867 | if (!event[iRecNow].isFinal()) continue; |
---|
868 | if ( (colSign > 0 && event[iRecNow].col() == colTagRec) |
---|
869 | || (colSign < 0 && event[iRecNow].acol() == colTagRec) ) { |
---|
870 | // Store recoilers in temporary array |
---|
871 | iRecVec.push_back(iRecNow); |
---|
872 | // Set iRec != 0 for checks below |
---|
873 | iRec = iRecNow; |
---|
874 | } |
---|
875 | } |
---|
876 | } |
---|
877 | |
---|
878 | } // End if-then-else of junction kinds |
---|
879 | |
---|
880 | } // End if leg has right color tag |
---|
881 | } // End of loop over junction legs |
---|
882 | } // End loop over junctions |
---|
883 | |
---|
884 | } // End main junction if |
---|
885 | |
---|
886 | // If fail, then other end to nearest recoiler in same system final state, |
---|
887 | // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j). |
---|
888 | if (iRec == 0) { |
---|
889 | double ppMin = LARGEM2; |
---|
890 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
891 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
892 | if (!event[iRecNow].isFinal()) continue; |
---|
893 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
894 | - event[iRecNow].m() * event[iRad].m(); |
---|
895 | if (ppNow < ppMin) { |
---|
896 | iRec = iRecNow; |
---|
897 | ppMin = ppNow; |
---|
898 | } |
---|
899 | } |
---|
900 | } |
---|
901 | |
---|
902 | // If fail, then other end to nearest recoiler in any system final state, |
---|
903 | // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j). |
---|
904 | if (iRec == 0) { |
---|
905 | double ppMin = LARGEM2; |
---|
906 | for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) |
---|
907 | if (iRecNow != iRad && event[iRecNow].isFinal()) { |
---|
908 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
909 | - event[iRecNow].m() * event[iRad].m(); |
---|
910 | if (ppNow < ppMin) { |
---|
911 | iRec = iRecNow; |
---|
912 | otherSystemRec = true; |
---|
913 | ppMin = ppNow; |
---|
914 | } |
---|
915 | } |
---|
916 | } |
---|
917 | |
---|
918 | // PS dec 2010: make sure iRec is stored in iRecVec |
---|
919 | if (iRecVec.size() == 0 && iRec != 0) iRecVec.push_back(iRec); |
---|
920 | |
---|
921 | // Remove any zero recoilers from normalization |
---|
922 | int nRec = iRecVec.size(); |
---|
923 | for (unsigned int mRec = 0; mRec < iRecVec.size(); ++mRec) |
---|
924 | if (iRecVec[mRec] <= 0) nRec--; |
---|
925 | if (nRec >= 2) { |
---|
926 | isFlexible = true; |
---|
927 | flexFactor = 1.0/nRec; |
---|
928 | } |
---|
929 | |
---|
930 | // Check for failure to locate any recoiler |
---|
931 | if ( nRec <= 0 ) { |
---|
932 | infoPtr->errorMsg("Error in TimeShower::setupQCDdip: " |
---|
933 | "failed to locate any recoiling partner"); |
---|
934 | return; |
---|
935 | } |
---|
936 | |
---|
937 | // Store dipole colour end(s). |
---|
938 | for (unsigned int mRec = 0; mRec < iRecVec.size(); ++mRec) { |
---|
939 | iRec = iRecVec[mRec]; |
---|
940 | if (iRec <= 0) continue; |
---|
941 | // Max scale either by parton scale or by half dipole mass. |
---|
942 | double pTmax = event[iRad].scale(); |
---|
943 | if (limitPTmaxIn) { |
---|
944 | if (iSys == 0) pTmax *= pTmaxFudge; |
---|
945 | if (iSys > 0 && sizeIn > 0) pTmax *= pTmaxFudgeMPI; |
---|
946 | } else pTmax = 0.5 * m( event[iRad], event[iRec]); |
---|
947 | int colType = (event[iRad].id() == 21) ? 2 * colSign : colSign; |
---|
948 | int isrType = (event[iRec].isFinal()) ? 0 : event[iRec].mother1(); |
---|
949 | // This line in case mother is a rescattered parton. |
---|
950 | while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); |
---|
951 | if (isrType > 2) isrType -= beamOffset; |
---|
952 | dipEnd.push_back( TimeDipoleEnd( iRad, iRec, pTmax, |
---|
953 | colType, 0, 0, isrType, iSys, -1, -1, isOctetOnium) ); |
---|
954 | |
---|
955 | // If hooked up with other system then find which. |
---|
956 | if (otherSystemRec) { |
---|
957 | int systemRec = partonSystemsPtr->getSystemOf(iRec, true); |
---|
958 | if (systemRec >= 0) dipEnd.back().systemRec = systemRec; |
---|
959 | dipEnd.back().MEtype = 0; |
---|
960 | } |
---|
961 | |
---|
962 | // PS dec 2010 |
---|
963 | // If non-unity (flexible) normalization, set normalization factor |
---|
964 | if (isFlexible) { |
---|
965 | dipEnd.back().isFlexible = true; |
---|
966 | dipEnd.back().flexFactor = flexFactor; |
---|
967 | } |
---|
968 | } |
---|
969 | |
---|
970 | } |
---|
971 | |
---|
972 | //-------------------------------------------------------------------------- |
---|
973 | |
---|
974 | // Setup a dipole end for a QED colour charge or a photon. |
---|
975 | // No failsafe choice of recoiler, so gradually widen search. |
---|
976 | |
---|
977 | void TimeShower::setupQEDdip( int iSys, int i, int chgType, int gamType, |
---|
978 | Event& event, bool limitPTmaxIn) { |
---|
979 | |
---|
980 | // Initial values. Find if allowed to hook up beams. |
---|
981 | int iRad = partonSystemsPtr->getOut(iSys, i); |
---|
982 | int idRad = event[iRad].id(); |
---|
983 | int iRec = 0; |
---|
984 | int sizeAllA = partonSystemsPtr->sizeAll(iSys); |
---|
985 | int sizeOut = partonSystemsPtr->sizeOut(iSys); |
---|
986 | int sizeAll = ( allowBeamRecoil ) ? sizeAllA : sizeOut; |
---|
987 | int sizeIn = sizeAll - sizeOut; |
---|
988 | int sizeInA = sizeAllA - sizeIn - sizeOut; |
---|
989 | int iOffset = i + sizeAllA - sizeOut; |
---|
990 | double ppMin = LARGEM2; |
---|
991 | bool hasRescattered = false; |
---|
992 | bool otherSystemRec = false; |
---|
993 | |
---|
994 | // Find nearest same- (opposide-) flavour recoiler in initial (final) |
---|
995 | // state of same system, excluding rescattered (in or out) partons. |
---|
996 | // Also find if system is involved in rescattering. |
---|
997 | // Note: (p_i + p_j)2 - (m_i + m_j)2 = 2 (p_i p_j - m_i m_j). |
---|
998 | for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { |
---|
999 | int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); |
---|
1000 | if ( (j < sizeIn && !event[iRecNow].isRescatteredIncoming()) |
---|
1001 | || (j >= sizeIn && event[iRecNow].isFinal()) ) { |
---|
1002 | if ( (j < sizeIn && event[iRecNow].id() == idRad) |
---|
1003 | || (j >= sizeIn && event[iRecNow].id() == -idRad) ) { |
---|
1004 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
1005 | - event[iRecNow].m() * event[iRad].m(); |
---|
1006 | if (ppNow < ppMin) { |
---|
1007 | iRec = iRecNow; |
---|
1008 | ppMin = ppNow; |
---|
1009 | } |
---|
1010 | } |
---|
1011 | } else hasRescattered = true; |
---|
1012 | } |
---|
1013 | |
---|
1014 | // If rescattering then find nearest opposite-flavour recoiler |
---|
1015 | // anywhere in final state. |
---|
1016 | if (iRec == 0 && hasRescattered) { |
---|
1017 | for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) |
---|
1018 | if (event[iRecNow].id() == -idRad && event[iRecNow].isFinal()) { |
---|
1019 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
1020 | - event[iRecNow].m() * event[iRad].m(); |
---|
1021 | if (ppNow < ppMin) { |
---|
1022 | iRec = iRecNow; |
---|
1023 | ppMin = ppNow; |
---|
1024 | otherSystemRec = true; |
---|
1025 | } |
---|
1026 | } |
---|
1027 | } |
---|
1028 | |
---|
1029 | // Find nearest recoiler in same system, charge-squared-weighted, |
---|
1030 | // including initial state, but excluding rescatterer. |
---|
1031 | if (iRec == 0) |
---|
1032 | for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { |
---|
1033 | int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); |
---|
1034 | int chgTypeRecNow = event[iRecNow].chargeType(); |
---|
1035 | if (chgTypeRecNow == 0) continue; |
---|
1036 | if ( (j < sizeIn && !event[iRecNow].isRescatteredIncoming()) |
---|
1037 | || (j >= sizeIn && event[iRecNow].isFinal()) ) { |
---|
1038 | double ppNow = (event[iRecNow].p() * event[iRad].p() |
---|
1039 | - event[iRecNow].m() * event[iRad].m()) |
---|
1040 | / pow2(chgTypeRecNow); |
---|
1041 | if (ppNow < ppMin) { |
---|
1042 | iRec = iRecNow; |
---|
1043 | ppMin = ppNow; |
---|
1044 | } |
---|
1045 | } |
---|
1046 | } |
---|
1047 | |
---|
1048 | // If rescattering then find nearest recoiler in the final state, |
---|
1049 | // charge-squared-weighted. |
---|
1050 | if (iRec == 0 && hasRescattered) { |
---|
1051 | for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) |
---|
1052 | if (iRecNow != iRad && event[iRecNow].isFinal()) { |
---|
1053 | int chgTypeRecNow = event[iRecNow].chargeType(); |
---|
1054 | if (chgTypeRecNow != 0 && event[iRecNow].isFinal()) { |
---|
1055 | double ppNow = (event[iRecNow].p() * event[iRad].p() |
---|
1056 | - event[iRecNow].m() * event[iRad].m()) |
---|
1057 | / pow2(chgTypeRecNow); |
---|
1058 | if (ppNow < ppMin) { |
---|
1059 | iRec = iRecNow; |
---|
1060 | ppMin = ppNow; |
---|
1061 | otherSystemRec = true; |
---|
1062 | } |
---|
1063 | } |
---|
1064 | } |
---|
1065 | } |
---|
1066 | |
---|
1067 | // Find any nearest recoiler in final state of same system. |
---|
1068 | if (iRec == 0) |
---|
1069 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
1070 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
1071 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
1072 | - event[iRecNow].m() * event[iRad].m(); |
---|
1073 | if (ppNow < ppMin) { |
---|
1074 | iRec = iRecNow; |
---|
1075 | ppMin = ppNow; |
---|
1076 | } |
---|
1077 | } |
---|
1078 | |
---|
1079 | // Find any nearest recoiler in final state. |
---|
1080 | if (iRec == 0) |
---|
1081 | for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) |
---|
1082 | if (iRecNow != iRad && event[iRecNow].isFinal()) { |
---|
1083 | double ppNow = event[iRecNow].p() * event[iRad].p() |
---|
1084 | - event[iRecNow].m() * event[iRad].m(); |
---|
1085 | if (ppNow < ppMin) { |
---|
1086 | iRec = iRecNow; |
---|
1087 | ppMin = ppNow; |
---|
1088 | otherSystemRec = true; |
---|
1089 | } |
---|
1090 | } |
---|
1091 | |
---|
1092 | // Fill charge-dipole or photon-dipole end. |
---|
1093 | if (iRec > 0) { |
---|
1094 | // Max scale either by parton scale or by half dipole mass. |
---|
1095 | double pTmax = event[iRad].scale(); |
---|
1096 | if (limitPTmaxIn) { |
---|
1097 | if (iSys == 0) pTmax *= pTmaxFudge; |
---|
1098 | if (iSys > 0 && sizeIn > 0) pTmax *= pTmaxFudgeMPI; |
---|
1099 | } else pTmax = 0.5 * m( event[iRad], event[iRec]); |
---|
1100 | int isrType = (event[iRec].isFinal()) ? 0 : event[iRec].mother1(); |
---|
1101 | // This line in case mother is a rescattered parton. |
---|
1102 | while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); |
---|
1103 | if (isrType > 2) isrType -= beamOffset; |
---|
1104 | dipEnd.push_back( TimeDipoleEnd(iRad, iRec, pTmax, |
---|
1105 | 0, chgType, gamType, isrType, iSys, -1) ); |
---|
1106 | |
---|
1107 | // If hooked up with other system then find which. |
---|
1108 | if (otherSystemRec) { |
---|
1109 | int systemRec = partonSystemsPtr->getSystemOf(iRec); |
---|
1110 | if (systemRec >= 0) dipEnd.back().systemRec = systemRec; |
---|
1111 | dipEnd.back().MEtype = 0; |
---|
1112 | } |
---|
1113 | |
---|
1114 | // Failure to find other end of dipole. |
---|
1115 | } else { |
---|
1116 | infoPtr->errorMsg("Error in TimeShower::setupQEDdip: " |
---|
1117 | "failed to locate any recoiling partner"); |
---|
1118 | } |
---|
1119 | |
---|
1120 | } |
---|
1121 | |
---|
1122 | //-------------------------------------------------------------------------- |
---|
1123 | |
---|
1124 | // Setup a dipole end for a Hidden Valley colour charge. |
---|
1125 | |
---|
1126 | void TimeShower::setupHVdip( int iSys, int i, Event& event, |
---|
1127 | bool limitPTmaxIn) { |
---|
1128 | |
---|
1129 | // Initial values. |
---|
1130 | int iRad = partonSystemsPtr->getOut(iSys, i); |
---|
1131 | int iRec = 0; |
---|
1132 | int idRad = event[iRad].id(); |
---|
1133 | int sizeOut = partonSystemsPtr->sizeOut(iSys); |
---|
1134 | |
---|
1135 | // Hidden Valley colour positive for positive id, and vice versa. |
---|
1136 | // Find opposte HV colour in final state of same system. |
---|
1137 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
1138 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
1139 | int idRec = event[iRecNow].id(); |
---|
1140 | if ( (abs(idRec) > 4900000 && abs(idRec) < 4900017) |
---|
1141 | && idRad * idRec < 0) { |
---|
1142 | iRec = iRecNow; |
---|
1143 | break; |
---|
1144 | } |
---|
1145 | } |
---|
1146 | |
---|
1147 | // Else find heaviest other final-state in same system. |
---|
1148 | // (Intended for decays; should mainly be two-body so unique.) |
---|
1149 | double mMax = -sqrt(LARGEM2); |
---|
1150 | if (iRec == 0) |
---|
1151 | for (int j = 0; j < sizeOut; ++j) if (j != i) { |
---|
1152 | int iRecNow = partonSystemsPtr->getOut(iSys, j); |
---|
1153 | if (event[iRecNow].m() > mMax) { |
---|
1154 | iRec = iRecNow; |
---|
1155 | mMax = event[iRecNow].m(); |
---|
1156 | } |
---|
1157 | } |
---|
1158 | |
---|
1159 | // Set up dipole end, or report failure. |
---|
1160 | if (iRec > 0) { |
---|
1161 | // Max scale either by parton scale or by half dipole mass. |
---|
1162 | double pTmax = event[iRad].scale(); |
---|
1163 | if (limitPTmaxIn) { |
---|
1164 | if (iSys == 0) pTmax *= pTmaxFudge; |
---|
1165 | } else pTmax = 0.5 * m( event[iRad], event[iRec]); |
---|
1166 | int colvType = (event[iRad].id() > 0) ? 1 : -1; |
---|
1167 | dipEnd.push_back( TimeDipoleEnd( iRad, iRec, pTmax, 0, 0, 0, 0, |
---|
1168 | iSys, -1, -1, false, true, colvType) ); |
---|
1169 | } else infoPtr->errorMsg("Error in TimeShower::setupHVdip: " |
---|
1170 | "failed to locate any recoiling partner"); |
---|
1171 | |
---|
1172 | } |
---|
1173 | |
---|
1174 | //-------------------------------------------------------------------------- |
---|
1175 | |
---|
1176 | // Select next pT in downwards evolution of the existing dipoles. |
---|
1177 | |
---|
1178 | double TimeShower::pTnext( Event& event, double pTbegAll, double pTendAll) { |
---|
1179 | |
---|
1180 | // Begin loop over all possible radiating dipole ends. |
---|
1181 | dipSel = 0; |
---|
1182 | iDipSel = -1; |
---|
1183 | double pT2sel = pTendAll * pTendAll; |
---|
1184 | for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) { |
---|
1185 | TimeDipoleEnd& dip = dipEnd[iDip]; |
---|
1186 | |
---|
1187 | // Check if global recoil should be used. |
---|
1188 | useLocalRecoilNow = !(globalRecoil && dip.system == 0 |
---|
1189 | && partonSystemsPtr->sizeOut(0) <= nMaxGlobalRecoil); |
---|
1190 | |
---|
1191 | // Dipole properties; normal local recoil. |
---|
1192 | dip.mRad = event[dip.iRadiator].m(); |
---|
1193 | if (useLocalRecoilNow) { |
---|
1194 | dip.mRec = event[dip.iRecoiler].m(); |
---|
1195 | dip.mDip = m( event[dip.iRadiator], event[dip.iRecoiler] ); |
---|
1196 | |
---|
1197 | // Dipole properties, alternative global recoil. Squares. |
---|
1198 | } else { |
---|
1199 | Vec4 pSumGlobal; |
---|
1200 | for (int i = 0; i < partonSystemsPtr->sizeOut( dip.system); ++i) { |
---|
1201 | int ii = partonSystemsPtr->getOut( dip.system, i); |
---|
1202 | if (ii != dip.iRadiator) pSumGlobal += event[ii].p(); |
---|
1203 | } |
---|
1204 | dip.mRec = pSumGlobal.mCalc(); |
---|
1205 | dip.mDip = m( event[dip.iRadiator].p(), pSumGlobal); |
---|
1206 | } |
---|
1207 | dip.m2Rad = pow2(dip.mRad); |
---|
1208 | dip.m2Rec = pow2(dip.mRec); |
---|
1209 | dip.m2Dip = pow2(dip.mDip); |
---|
1210 | |
---|
1211 | // Find maximum evolution scale for dipole. |
---|
1212 | dip.m2DipCorr = pow2(dip.mDip - dip.mRec) - dip.m2Rad; |
---|
1213 | double pTbegDip = min( pTbegAll, dip.pTmax ); |
---|
1214 | double pT2begDip = min( pow2(pTbegDip), 0.25 * dip.m2DipCorr); |
---|
1215 | |
---|
1216 | // Do QCD, QED or HV evolution if it makes sense. |
---|
1217 | dip.pT2 = 0.; |
---|
1218 | if (pT2begDip > pT2sel) { |
---|
1219 | if (dip.colType != 0) |
---|
1220 | pT2nextQCD(pT2begDip, pT2sel, dip, event); |
---|
1221 | else if (dip.chgType != 0 || dip.gamType != 0) |
---|
1222 | pT2nextQED(pT2begDip, pT2sel, dip, event); |
---|
1223 | else if (dip.colvType != 0) |
---|
1224 | pT2nextHV(pT2begDip, pT2sel, dip, event); |
---|
1225 | |
---|
1226 | // Update if found larger pT than current maximum. |
---|
1227 | if (dip.pT2 > pT2sel) { |
---|
1228 | pT2sel = dip.pT2; |
---|
1229 | dipSel = &dip; |
---|
1230 | iDipSel = iDip; |
---|
1231 | } |
---|
1232 | } |
---|
1233 | } |
---|
1234 | |
---|
1235 | // Return nonvanishing value if found pT bigger than already found. |
---|
1236 | return (dipSel == 0) ? 0. : sqrt(pT2sel); |
---|
1237 | |
---|
1238 | } |
---|
1239 | |
---|
1240 | //-------------------------------------------------------------------------- |
---|
1241 | |
---|
1242 | // Evolve a QCD dipole end. |
---|
1243 | |
---|
1244 | void TimeShower::pT2nextQCD(double pT2begDip, double pT2sel, |
---|
1245 | TimeDipoleEnd& dip, Event& event) { |
---|
1246 | |
---|
1247 | // Lower cut for evolution. Return if no evolution range. |
---|
1248 | double pT2endDip = max( pT2sel, pT2colCut ); |
---|
1249 | if (pT2begDip < pT2endDip) return; |
---|
1250 | |
---|
1251 | // Upper estimate for matrix element weighting and colour factor. |
---|
1252 | // Special cases for triplet recoiling against gluino and octet onia. |
---|
1253 | // Note that g -> g g and g -> q qbar are split on two sides. |
---|
1254 | int colTypeAbs = abs(dip.colType); |
---|
1255 | double wtPSglue = 2.; |
---|
1256 | double colFac = (colTypeAbs == 1) ? 4./3. : 3./2.; |
---|
1257 | if (dip.MEgluinoRec) colFac = 3.; |
---|
1258 | if (dip.isOctetOnium) colFac *= 0.5 * octetOniumColFac; |
---|
1259 | // PS dec 2010. Include possibility for flexible normalization, |
---|
1260 | // e.g., for dipoles stretched to junctions or to switch off radiation. |
---|
1261 | if (dip.isFlexible) colFac *= dip.flexFactor; |
---|
1262 | double wtPSqqbar = (colTypeAbs == 2) ? 0.25 * nGluonToQuark : 0.; |
---|
1263 | |
---|
1264 | // Variables used inside evolution loop. (Mainly dummy start values.) |
---|
1265 | dip.pT2 = pT2begDip; |
---|
1266 | int nFlavour = 3; |
---|
1267 | double zMinAbs = 0.5; |
---|
1268 | double pT2min = pT2endDip; |
---|
1269 | double b0 = 4.5; |
---|
1270 | double Lambda2 = Lambda3flav2; |
---|
1271 | double emitCoefGlue = 0.; |
---|
1272 | double emitCoefQqbar = 0.; |
---|
1273 | double emitCoefTot = 0.; |
---|
1274 | double wt = 0.; |
---|
1275 | bool mustFindRange = true; |
---|
1276 | |
---|
1277 | // Begin evolution loop towards smaller pT values. |
---|
1278 | do { |
---|
1279 | |
---|
1280 | // Initialize evolution coefficients at the beginning and |
---|
1281 | // reinitialize when crossing c and b flavour thresholds. |
---|
1282 | if (mustFindRange) { |
---|
1283 | |
---|
1284 | // Determine overestimated z range; switch at c and b masses. |
---|
1285 | if (dip.pT2 > m2b) { |
---|
1286 | nFlavour = 5; |
---|
1287 | pT2min = max( m2b, pT2endDip); |
---|
1288 | b0 = 23./6.; |
---|
1289 | Lambda2 = Lambda5flav2; |
---|
1290 | } else if (dip.pT2 > m2c) { |
---|
1291 | nFlavour = 4; |
---|
1292 | pT2min = max( m2c, pT2endDip); |
---|
1293 | b0 = 25./6.; |
---|
1294 | Lambda2 = Lambda4flav2; |
---|
1295 | } else { |
---|
1296 | nFlavour = 3; |
---|
1297 | pT2min = pT2endDip; |
---|
1298 | b0 = 27./6.; |
---|
1299 | Lambda2 = Lambda3flav2; |
---|
1300 | } |
---|
1301 | // A change of renormalization scale expressed by a change of Lambda. |
---|
1302 | Lambda2 /= renormMultFac; |
---|
1303 | zMinAbs = 0.5 - sqrtpos( 0.25 - pT2min / dip.m2DipCorr ); |
---|
1304 | if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2min / dip.m2DipCorr; |
---|
1305 | |
---|
1306 | // Find emission coefficients for X -> X g and g -> q qbar. |
---|
1307 | emitCoefGlue = wtPSglue * colFac * log(1. / zMinAbs - 1.); |
---|
1308 | emitCoefTot = emitCoefGlue; |
---|
1309 | if (colTypeAbs == 2 && event[dip.iRadiator].id() == 21) { |
---|
1310 | emitCoefQqbar = wtPSqqbar * (1. - 2. * zMinAbs); |
---|
1311 | emitCoefTot += emitCoefQqbar; |
---|
1312 | } |
---|
1313 | |
---|
1314 | // Initialization done for current range. |
---|
1315 | mustFindRange = false; |
---|
1316 | } |
---|
1317 | |
---|
1318 | // Pick pT2 (in overestimated z range) for fixed alpha_strong. |
---|
1319 | if (alphaSorder == 0) { |
---|
1320 | dip.pT2 = dip.pT2 * pow( rndmPtr->flat(), |
---|
1321 | 1. / (alphaS2pi * emitCoefTot) ); |
---|
1322 | |
---|
1323 | // Ditto for first-order alpha_strong. |
---|
1324 | } else if (alphaSorder == 1) { |
---|
1325 | dip.pT2 = Lambda2 * pow( dip.pT2 / Lambda2, |
---|
1326 | pow( rndmPtr->flat(), b0 / emitCoefTot) ); |
---|
1327 | |
---|
1328 | // For second order reject by second term in alpha_strong expression. |
---|
1329 | } else { |
---|
1330 | do dip.pT2 = Lambda2 * pow( dip.pT2 / Lambda2, |
---|
1331 | pow( rndmPtr->flat(), b0 / emitCoefTot) ); |
---|
1332 | while (alphaS.alphaS2OrdCorr(renormMultFac * dip.pT2) < rndmPtr->flat() |
---|
1333 | && dip.pT2 > pT2min); |
---|
1334 | } |
---|
1335 | wt = 0.; |
---|
1336 | |
---|
1337 | // If crossed c or b thresholds: continue evolution from threshold. |
---|
1338 | if (nFlavour == 5 && dip.pT2 < m2b) { |
---|
1339 | mustFindRange = true; |
---|
1340 | dip.pT2 = m2b; |
---|
1341 | } else if ( nFlavour == 4 && dip.pT2 < m2c) { |
---|
1342 | mustFindRange = true; |
---|
1343 | dip.pT2 = m2c; |
---|
1344 | |
---|
1345 | // Abort evolution if below cutoff scale, or below another branching. |
---|
1346 | } else { |
---|
1347 | if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; } |
---|
1348 | |
---|
1349 | // Pick kind of branching: X -> X g or g -> q qbar. |
---|
1350 | dip.flavour = 21; |
---|
1351 | dip.mFlavour = 0.; |
---|
1352 | if (colTypeAbs == 2 && emitCoefQqbar > rndmPtr->flat() |
---|
1353 | * emitCoefTot) dip.flavour = 0; |
---|
1354 | |
---|
1355 | // Pick z: either dz/(1-z) or flat dz. |
---|
1356 | if (dip.flavour == 21) { |
---|
1357 | dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() ); |
---|
1358 | } else { |
---|
1359 | dip.z = zMinAbs + (1. - 2. * zMinAbs) * rndmPtr->flat(); |
---|
1360 | } |
---|
1361 | |
---|
1362 | // Do not accept branching if outside allowed z range. |
---|
1363 | double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr ); |
---|
1364 | if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr; |
---|
1365 | dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z)); |
---|
1366 | if (dip.z > zMin && dip.z < 1. - zMin |
---|
1367 | && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z) |
---|
1368 | * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) ) { |
---|
1369 | |
---|
1370 | // Flavour choice for g -> q qbar. |
---|
1371 | if (dip.flavour == 0) { |
---|
1372 | dip.flavour = min(5, 1 + int(nGluonToQuark * rndmPtr->flat())); |
---|
1373 | dip.mFlavour = particleDataPtr->m0(dip.flavour); |
---|
1374 | } |
---|
1375 | |
---|
1376 | // No z weight, except threshold, if to do ME corrections later on. |
---|
1377 | if (dip.MEtype > 0) { |
---|
1378 | wt = 1.; |
---|
1379 | if (dip.flavour < 10 && dip.m2 < THRESHM2 * pow2(dip.mFlavour)) |
---|
1380 | wt = 0.; |
---|
1381 | |
---|
1382 | // z weight for X -> X g. |
---|
1383 | } else if (dip.flavour == 21 |
---|
1384 | && (colTypeAbs == 1 || colTypeAbs == 3) ) { |
---|
1385 | wt = (1. + pow2(dip.z)) / wtPSglue; |
---|
1386 | } else if (dip.flavour == 21) { |
---|
1387 | wt = (1. + pow3(dip.z)) / wtPSglue; |
---|
1388 | |
---|
1389 | // z weight for g -> q qbar. |
---|
1390 | } else { |
---|
1391 | double beta = sqrtpos( 1. - 4. * pow2(dip.mFlavour) / dip.m2 ); |
---|
1392 | wt = beta * ( pow2(dip.z) + pow2(1. - dip.z) ); |
---|
1393 | } |
---|
1394 | |
---|
1395 | // Suppression factors for dipole to beam remnant. |
---|
1396 | if (dip.isrType != 0 && useLocalRecoilNow) { |
---|
1397 | BeamParticle& beam = (dip.isrType == 1) ? *beamAPtr : *beamBPtr; |
---|
1398 | int iSysRec = dip.systemRec; |
---|
1399 | double xOld = beam[iSysRec].x(); |
---|
1400 | double xNew = xOld * (1. + (dip.m2 - dip.m2Rad) / |
---|
1401 | (dip.m2Dip - dip.m2Rad)); |
---|
1402 | double xMaxAbs = beam.xMax(iSysRec); |
---|
1403 | if (xMaxAbs < 0.) { |
---|
1404 | infoPtr->errorMsg("Warning in TimeShower::pT2nextQCD: " |
---|
1405 | "xMaxAbs negative"); |
---|
1406 | return; |
---|
1407 | } |
---|
1408 | |
---|
1409 | // Firstly reduce by PDF ratio. |
---|
1410 | if (xNew > xMaxAbs) wt = 0.; |
---|
1411 | else { |
---|
1412 | int idRec = event[dip.iRecoiler].id(); |
---|
1413 | double pdfOld = max ( TINYPDF, |
---|
1414 | beam.xfISR( iSysRec, idRec, xOld, factorMultFac * dip.pT2) ); |
---|
1415 | double pdfNew = |
---|
1416 | beam.xfISR( iSysRec, idRec, xNew, factorMultFac * dip.pT2); |
---|
1417 | wt *= min( 1., pdfNew / pdfOld); |
---|
1418 | } |
---|
1419 | |
---|
1420 | // Secondly optionally reduce by 4 pT2_hard / (4 pT2_hard + m2). |
---|
1421 | if (dampenBeamRecoil) { |
---|
1422 | double pTpT = sqrt(event[dip.iRadiator].pT2() * dip.pT2); |
---|
1423 | wt *= pTpT / (pTpT + dip.m2); |
---|
1424 | } |
---|
1425 | } |
---|
1426 | } |
---|
1427 | } |
---|
1428 | |
---|
1429 | // Iterate until acceptable pT (or have fallen below pTmin). |
---|
1430 | } while (wt < rndmPtr->flat()); |
---|
1431 | |
---|
1432 | } |
---|
1433 | |
---|
1434 | //-------------------------------------------------------------------------- |
---|
1435 | |
---|
1436 | // Evolve a QED dipole end, either charged or photon. |
---|
1437 | |
---|
1438 | void TimeShower::pT2nextQED(double pT2begDip, double pT2sel, |
---|
1439 | TimeDipoleEnd& dip, Event& event) { |
---|
1440 | |
---|
1441 | // Lower cut for evolution. Return if no evolution range. |
---|
1442 | double pT2chgCut = (dip.chgType != 0 && abs(dip.chgType) != 3) |
---|
1443 | ? pT2chgQCut : pT2chgLCut; |
---|
1444 | double pT2endDip = max( pT2sel, pT2chgCut ); |
---|
1445 | if (pT2begDip < pT2endDip) return; |
---|
1446 | |
---|
1447 | // Emission of photon or photon branching. |
---|
1448 | bool hasCharge = (dip.chgType != 0); |
---|
1449 | |
---|
1450 | // Default values. |
---|
1451 | double wtPSgam = 0.; |
---|
1452 | double chg2Sum = 0.; |
---|
1453 | double chg2SumL = 0.; |
---|
1454 | double chg2SumQ = 0.; |
---|
1455 | double zMinAbs = 0.; |
---|
1456 | double emitCoefTot = 0.; |
---|
1457 | |
---|
1458 | // alpha_em at maximum scale provides upper estimate. |
---|
1459 | double alphaEMmax = alphaEM.alphaEM(renormMultFac * pT2begDip); |
---|
1460 | double alphaEM2pi = alphaEMmax / (2. * M_PI); |
---|
1461 | |
---|
1462 | // Emission: upper estimate for matrix element weighting; charge factor. |
---|
1463 | if (hasCharge) { |
---|
1464 | wtPSgam = 2.; |
---|
1465 | double chg2 = pow2(dip.chgType / 3.); |
---|
1466 | |
---|
1467 | // Determine overestimated z range. Find evolution coefficient. |
---|
1468 | zMinAbs = 0.5 - sqrtpos( 0.25 - pT2endDip / dip.m2DipCorr ); |
---|
1469 | if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2endDip / dip.m2DipCorr; |
---|
1470 | emitCoefTot = alphaEM2pi * chg2 * wtPSgam * log(1. / zMinAbs - 1.); |
---|
1471 | |
---|
1472 | // Branching: sum of squared charge factors for lepton and quark daughters. |
---|
1473 | } else { |
---|
1474 | chg2SumL = max(0, min(3, nGammaToLepton)); |
---|
1475 | if (nGammaToQuark > 4) chg2SumQ = 11. / 9.; |
---|
1476 | else if (nGammaToQuark > 3) chg2SumQ = 10. / 9.; |
---|
1477 | else if (nGammaToQuark > 2) chg2SumQ = 6. / 9.; |
---|
1478 | else if (nGammaToQuark > 1) chg2SumQ = 5. / 9.; |
---|
1479 | else if (nGammaToQuark > 0) chg2SumQ = 1. / 9.; |
---|
1480 | |
---|
1481 | // Total sum of squared charge factors. Find evolution coefficient. |
---|
1482 | chg2Sum = chg2SumL + 3. * chg2SumQ; |
---|
1483 | emitCoefTot = alphaEM2pi * chg2Sum; |
---|
1484 | } |
---|
1485 | |
---|
1486 | // Variables used inside evolution loop. |
---|
1487 | dip.pT2 = pT2begDip; |
---|
1488 | double wt; |
---|
1489 | |
---|
1490 | // Begin evolution loop towards smaller pT values. |
---|
1491 | do { |
---|
1492 | |
---|
1493 | // Pick pT2 (in overestimated z range). |
---|
1494 | dip.pT2 = dip.pT2 * pow(rndmPtr->flat(), 1. / emitCoefTot); |
---|
1495 | wt = 0.; |
---|
1496 | |
---|
1497 | // Abort evolution if below cutoff scale, or below another branching. |
---|
1498 | if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; } |
---|
1499 | |
---|
1500 | // Pick z according to dz/(1-z) or flat. |
---|
1501 | if (hasCharge) dip.z = 1. - zMinAbs |
---|
1502 | * pow( 1. / zMinAbs - 1., rndmPtr->flat() ); |
---|
1503 | else dip.z = rndmPtr->flat(); |
---|
1504 | |
---|
1505 | // Do not accept branching if outside allowed z range. |
---|
1506 | double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr ); |
---|
1507 | if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr; |
---|
1508 | dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z)); |
---|
1509 | if (dip.z > zMin && dip.z < 1. - zMin |
---|
1510 | && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z) |
---|
1511 | * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) |
---|
1512 | // For gamma -> f fbar also impose maximum mass. |
---|
1513 | && (hasCharge || dip.m2 < m2MaxGamma) ) { |
---|
1514 | |
---|
1515 | // Photon emission: unique flavour choice. |
---|
1516 | if (hasCharge) { |
---|
1517 | dip.flavour = 22; |
---|
1518 | dip.mFlavour = 0.; |
---|
1519 | |
---|
1520 | // Photon branching: either lepton or quark flavour choice. |
---|
1521 | } else { |
---|
1522 | if (rndmPtr->flat() * chg2Sum < chg2SumL) |
---|
1523 | dip.flavour = 9 + 2 * min(3, 1 + int(chg2SumL * rndmPtr->flat())); |
---|
1524 | else { |
---|
1525 | double rndmQ = 9. * chg2SumQ * rndmPtr->flat(); |
---|
1526 | if (rndmQ < 1.) dip.flavour = 1; |
---|
1527 | else if (rndmQ < 5.) dip.flavour = 2; |
---|
1528 | else if (rndmQ < 6.) dip.flavour = 3; |
---|
1529 | else if (rndmQ < 10.) dip.flavour = 4; |
---|
1530 | else dip.flavour = 5; |
---|
1531 | } |
---|
1532 | dip.mFlavour = particleDataPtr->m0(dip.flavour); |
---|
1533 | } |
---|
1534 | |
---|
1535 | // No z weight, except threshold, if to do ME corrections later on. |
---|
1536 | if (dip.MEtype > 0) { |
---|
1537 | wt = 1.; |
---|
1538 | if (dip.flavour < 20 && dip.m2 < THRESHM2 * pow2(dip.mFlavour)) |
---|
1539 | wt = 0.; |
---|
1540 | |
---|
1541 | // z weight for X -> X gamma. |
---|
1542 | } else if (hasCharge) { |
---|
1543 | wt = (1. + pow2(dip.z)) / wtPSgam; |
---|
1544 | |
---|
1545 | // z weight for gamma -> f fbar. |
---|
1546 | } else { |
---|
1547 | double beta = sqrtpos( 1. - 4. * pow2(dip.mFlavour) / dip.m2 ); |
---|
1548 | wt = beta * ( pow2(dip.z) + pow2(1. - dip.z) ); |
---|
1549 | } |
---|
1550 | |
---|
1551 | // Correct to current value of alpha_EM. |
---|
1552 | double alphaEMnow = alphaEM.alphaEM(renormMultFac * dip.pT2); |
---|
1553 | wt *= (alphaEMnow / alphaEMmax); |
---|
1554 | |
---|
1555 | // Suppression factors for dipole to beam remnant. |
---|
1556 | if (dip.isrType != 0 && useLocalRecoilNow) { |
---|
1557 | BeamParticle& beam = (dip.isrType == 1) ? *beamAPtr : *beamBPtr; |
---|
1558 | int iSys = dip.system; |
---|
1559 | double xOld = beam[iSys].x(); |
---|
1560 | double xNew = xOld * (1. + (dip.m2 - dip.m2Rad) / |
---|
1561 | (dip.m2Dip - dip.m2Rad)); |
---|
1562 | double xMaxAbs = beam.xMax(iSys); |
---|
1563 | if (xMaxAbs < 0.) { |
---|
1564 | infoPtr->errorMsg("Warning in TimeShower::pT2nextQED: " |
---|
1565 | "xMaxAbs negative"); |
---|
1566 | return; |
---|
1567 | } |
---|
1568 | |
---|
1569 | // Firstly reduce by PDF ratio. |
---|
1570 | if (xNew > xMaxAbs) wt = 0.; |
---|
1571 | else { |
---|
1572 | int idRec = event[dip.iRecoiler].id(); |
---|
1573 | double pdfOld = max ( TINYPDF, |
---|
1574 | beam.xfISR( iSys, idRec, xOld, factorMultFac * dip.pT2) ); |
---|
1575 | double pdfNew = |
---|
1576 | beam.xfISR( iSys, idRec, xNew, factorMultFac * dip.pT2); |
---|
1577 | wt *= min( 1., pdfNew / pdfOld); |
---|
1578 | } |
---|
1579 | |
---|
1580 | // Secondly optionally reduce by 4 pT2_hard / (4 pT2_hard + m2). |
---|
1581 | if (dampenBeamRecoil) { |
---|
1582 | double pT24 = 4. * event[dip.iRadiator].pT2(); |
---|
1583 | wt *= pT24 / (pT24 + dip.m2); |
---|
1584 | } |
---|
1585 | } |
---|
1586 | } |
---|
1587 | |
---|
1588 | // Iterate until acceptable pT (or have fallen below pTmin). |
---|
1589 | } while (wt < rndmPtr->flat()); |
---|
1590 | |
---|
1591 | } |
---|
1592 | |
---|
1593 | //-------------------------------------------------------------------------- |
---|
1594 | |
---|
1595 | // Evolve a Hidden Valley dipole end. |
---|
1596 | |
---|
1597 | void TimeShower::pT2nextHV(double pT2begDip, double pT2sel, |
---|
1598 | TimeDipoleEnd& dip, Event& ) { |
---|
1599 | |
---|
1600 | // Lower cut for evolution. Return if no evolution range. |
---|
1601 | double pT2endDip = max( pT2sel, pT2hvCut ); |
---|
1602 | if (pT2begDip < pT2endDip) return; |
---|
1603 | |
---|
1604 | // C_F * alpha_HV/2 pi. |
---|
1605 | int colvTypeAbs = abs(dip.colvType); |
---|
1606 | double colvFac = (colvTypeAbs == 1) ? CFHV : 0.5 * nCHV; |
---|
1607 | double alphaHV2pi = colvFac * (alphaHVfix / (2. * M_PI)); |
---|
1608 | |
---|
1609 | // Determine overestimated z range. Find evolution coefficient. |
---|
1610 | double zMinAbs = 0.5 - sqrtpos( 0.25 - pT2endDip / dip.m2DipCorr ); |
---|
1611 | if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2endDip / dip.m2DipCorr; |
---|
1612 | double emitCoefTot = alphaHV2pi * 2. * log(1. / zMinAbs - 1.); |
---|
1613 | |
---|
1614 | // Variables used inside evolution loop. |
---|
1615 | dip.pT2 = pT2begDip; |
---|
1616 | double wt; |
---|
1617 | |
---|
1618 | // Begin evolution loop towards smaller pT values. |
---|
1619 | do { |
---|
1620 | |
---|
1621 | // Pick pT2 (in overestimated z range). |
---|
1622 | dip.pT2 = dip.pT2 * pow(rndmPtr->flat(), 1. / emitCoefTot); |
---|
1623 | wt = 0.; |
---|
1624 | |
---|
1625 | // Abort evolution if below cutoff scale, or below another branching. |
---|
1626 | if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; } |
---|
1627 | |
---|
1628 | // Pick z according to dz/(1-z). |
---|
1629 | dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() ); |
---|
1630 | |
---|
1631 | // Do not accept branching if outside allowed z range. |
---|
1632 | double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr ); |
---|
1633 | if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr; |
---|
1634 | dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z)); |
---|
1635 | if (dip.z > zMin && dip.z < 1. - zMin |
---|
1636 | && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z) |
---|
1637 | * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) ) { |
---|
1638 | |
---|
1639 | // HV gamma or gluon emission: unique flavour choice. |
---|
1640 | dip.flavour = idHV; |
---|
1641 | dip.mFlavour = mHV; |
---|
1642 | |
---|
1643 | // No z weight, except threshold, if to do ME corrections later on. |
---|
1644 | if (dip.MEtype > 0) wt = 1.; |
---|
1645 | |
---|
1646 | // z weight for X -> X g_HV. |
---|
1647 | else if (colvTypeAbs == 1) wt = (1. + pow2(dip.z)) / 2.; |
---|
1648 | else wt = (1. + pow3(dip.z)) / 2.; |
---|
1649 | } |
---|
1650 | |
---|
1651 | // Iterate until acceptable pT (or have fallen below pTmin). |
---|
1652 | } while (wt < rndmPtr->flat()); |
---|
1653 | |
---|
1654 | } |
---|
1655 | |
---|
1656 | //-------------------------------------------------------------------------- |
---|
1657 | |
---|
1658 | // ME corrections and kinematics that may give failure. |
---|
1659 | // Notation: radBef, recBef = radiator, recoiler before emission, |
---|
1660 | // rad, rec, emt = radiator, recoiler, emitted efter emission. |
---|
1661 | // (rad, emt distinguished by colour flow for g -> q qbar.) |
---|
1662 | |
---|
1663 | bool TimeShower::branch( Event& event, bool isInterleaved) { |
---|
1664 | |
---|
1665 | // Check if global recoil should be used. |
---|
1666 | useLocalRecoilNow = !(globalRecoil && dipSel->system == 0 |
---|
1667 | && partonSystemsPtr->sizeOut(0) <= nMaxGlobalRecoil); |
---|
1668 | |
---|
1669 | // Find initial radiator and recoiler particles in dipole branching. |
---|
1670 | int iRadBef = dipSel->iRadiator; |
---|
1671 | int iRecBef = dipSel->iRecoiler; |
---|
1672 | Particle& radBef = event[iRadBef]; |
---|
1673 | Particle& recBef = event[iRecBef]; |
---|
1674 | |
---|
1675 | // Find their momenta, with special sum for global recoil. |
---|
1676 | Vec4 pRadBef = event[iRadBef].p(); |
---|
1677 | Vec4 pRecBef; |
---|
1678 | vector<int> iGRecBef, iGRec; |
---|
1679 | if (useLocalRecoilNow) pRecBef = event[iRecBef].p(); |
---|
1680 | else { |
---|
1681 | for (int i = 0; i < partonSystemsPtr->sizeOut( dipSel->system); ++i) { |
---|
1682 | int iG = partonSystemsPtr->getOut( dipSel->system, i); |
---|
1683 | if (iG != dipSel->iRadiator) { |
---|
1684 | iGRecBef.push_back(iG); |
---|
1685 | pRecBef += event[iG].p(); |
---|
1686 | } |
---|
1687 | } |
---|
1688 | } |
---|
1689 | |
---|
1690 | // Default flavours and colour tags for new particles in dipole branching. |
---|
1691 | int idRad = radBef.id(); |
---|
1692 | int idEmt = dipSel->flavour; |
---|
1693 | int colRad = radBef.col(); |
---|
1694 | int acolRad = radBef.acol(); |
---|
1695 | int colEmt = 0; |
---|
1696 | int acolEmt = 0; |
---|
1697 | iSysSel = dipSel->system; |
---|
1698 | int iSysSelRec = dipSel->systemRec; |
---|
1699 | |
---|
1700 | // Default OK for photon, photon_HV or gluon_HV emission. |
---|
1701 | if (dipSel->flavour == 22 || dipSel->flavour == idHV) { |
---|
1702 | // New colour tag required for gluon emission. |
---|
1703 | } else if (dipSel->flavour == 21 && dipSel->colType > 0) { |
---|
1704 | colEmt = colRad; |
---|
1705 | colRad = event.nextColTag(); |
---|
1706 | acolEmt = colRad; |
---|
1707 | } else if (dipSel->flavour == 21) { |
---|
1708 | acolEmt = acolRad; |
---|
1709 | acolRad = event.nextColTag(); |
---|
1710 | colEmt = acolRad; |
---|
1711 | // New flavours for g -> q qbar; split colours. |
---|
1712 | } else if (dipSel->colType > 0) { |
---|
1713 | idEmt = dipSel->flavour ; |
---|
1714 | idRad = -idEmt; |
---|
1715 | colEmt = colRad; |
---|
1716 | colRad = 0; |
---|
1717 | } else if (dipSel->colType < 0) { |
---|
1718 | idEmt = -dipSel->flavour ; |
---|
1719 | idRad = -idEmt; |
---|
1720 | acolEmt = acolRad; |
---|
1721 | acolRad = 0; |
---|
1722 | // New flavours for gamma -> f fbar, and maybe also colours. |
---|
1723 | } else if (dipSel->gamType == 1 && rndmPtr->flat() > 0.5) { |
---|
1724 | idEmt = -dipSel->flavour ; |
---|
1725 | idRad = -idEmt; |
---|
1726 | if (idRad < 10) colRad = event.nextColTag(); |
---|
1727 | acolEmt = colRad; |
---|
1728 | } else if (dipSel->gamType == 1) { |
---|
1729 | idEmt = dipSel->flavour ; |
---|
1730 | idRad = -idEmt; |
---|
1731 | if (idEmt < 10) colEmt = event.nextColTag(); |
---|
1732 | acolRad = colEmt; |
---|
1733 | } |
---|
1734 | |
---|
1735 | // Construct kinematics in dipole rest frame: |
---|
1736 | // begin simple (like g -> g g). |
---|
1737 | double pTorig = sqrt( dipSel->pT2); |
---|
1738 | double eRadPlusEmt = 0.5 * (dipSel->m2Dip + dipSel->m2 - dipSel->m2Rec) |
---|
1739 | / dipSel->mDip; |
---|
1740 | double e2RadPlusEmt = pow2(eRadPlusEmt); |
---|
1741 | double pzRadPlusEmt = 0.5 * sqrtpos( pow2(dipSel->m2Dip - dipSel->m2 |
---|
1742 | - dipSel->m2Rec) - 4. * dipSel->m2 * dipSel->m2Rec ) / dipSel->mDip; |
---|
1743 | double pT2corr = dipSel->m2 * (e2RadPlusEmt * dipSel->z * (1. - dipSel->z) |
---|
1744 | - 0.25 * dipSel->m2) / pow2(pzRadPlusEmt); |
---|
1745 | double pTcorr = sqrtpos( pT2corr ); |
---|
1746 | double pzRad = (e2RadPlusEmt * dipSel->z - 0.5 * dipSel->m2) |
---|
1747 | / pzRadPlusEmt; |
---|
1748 | double pzEmt = (e2RadPlusEmt * (1. - dipSel->z) - 0.5 * dipSel->m2) |
---|
1749 | / pzRadPlusEmt; |
---|
1750 | double mRad = dipSel->mRad; |
---|
1751 | double mEmt = 0.; |
---|
1752 | |
---|
1753 | // Kinematics reduction for q -> q gamma_v when m_q > 0 and m_gamma_v > 0. |
---|
1754 | if ( abs(dipSel->colvType) == 1 && dipSel->mFlavour > 0.) { |
---|
1755 | mEmt = dipSel->mFlavour; |
---|
1756 | if (pow2(mRad + mEmt) > dipSel->m2) return false; |
---|
1757 | double m2Emt = pow2(mEmt); |
---|
1758 | double lambda = sqrtpos( pow2(dipSel->m2 - dipSel->m2Rad - m2Emt) |
---|
1759 | - 4. * dipSel->m2Rad * m2Emt ); |
---|
1760 | kRad = 0.5 * (dipSel->m2 - lambda + m2Emt - dipSel->m2Rad) |
---|
1761 | / dipSel->m2; |
---|
1762 | kEmt = 0.5 * (dipSel->m2 - lambda + dipSel->m2Rad - m2Emt) |
---|
1763 | / dipSel->m2; |
---|
1764 | pTorig *= 1. - kRad - kEmt; |
---|
1765 | pTcorr *= 1. - kRad - kEmt; |
---|
1766 | double pzMove = kRad * pzRad - kEmt * pzEmt; |
---|
1767 | pzRad -= pzMove; |
---|
1768 | pzEmt += pzMove; |
---|
1769 | |
---|
1770 | // Kinematics reduction for q -> q g/gamma/g_HV when m_q > 0. |
---|
1771 | } else if (abs(dipSel->colType) == 1 || dipSel->chgType != 0 |
---|
1772 | || abs(dipSel->colvType) == 1) { |
---|
1773 | pTorig *= 1. - dipSel->m2Rad / dipSel->m2; |
---|
1774 | pTcorr *= 1. - dipSel->m2Rad / dipSel->m2; |
---|
1775 | pzRad += pzEmt * dipSel->m2Rad / dipSel->m2; |
---|
1776 | pzEmt *= 1. - dipSel->m2Rad / dipSel->m2; |
---|
1777 | |
---|
1778 | // Kinematics reduction for g -> q qbar or gamma -> f fbar when m_f > 0; |
---|
1779 | } else if (abs(dipSel->flavour) < 20) { |
---|
1780 | mEmt = dipSel->mFlavour; |
---|
1781 | mRad = mEmt; |
---|
1782 | double beta = sqrtpos( 1. - 4. * pow2(mEmt) / dipSel->m2 ); |
---|
1783 | pTorig *= beta; |
---|
1784 | pTcorr *= beta; |
---|
1785 | pzRad = 0.5 * ( (1. + beta) * pzRad + (1. - beta) * pzEmt ); |
---|
1786 | pzEmt = pzRadPlusEmt - pzRad; |
---|
1787 | } |
---|
1788 | |
---|
1789 | // Reject g emission where mass effects have reduced pT below cutoff. |
---|
1790 | if (idEmt == 21 && pTorig < pTcolCut) return false; |
---|
1791 | |
---|
1792 | // Find rest frame and angles of original dipole. |
---|
1793 | RotBstMatrix M; |
---|
1794 | M.fromCMframe(pRadBef, pRecBef); |
---|
1795 | |
---|
1796 | // Evaluate coefficient of azimuthal asymmetry from gluon polarization. |
---|
1797 | findAsymPol( event, dipSel); |
---|
1798 | |
---|
1799 | // Begin construction of new dipole kinematics: pick azimuthal angle. |
---|
1800 | Vec4 pRad, pEmt, pRec; |
---|
1801 | double wtPhi = 1.; |
---|
1802 | do { |
---|
1803 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
1804 | |
---|
1805 | // Define kinematics of branching in dipole rest frame. |
---|
1806 | pRad = Vec4( pTcorr * cos(phi), pTcorr * sin(phi), pzRad, |
---|
1807 | sqrt( pow2(pTcorr) + pow2(pzRad) + pow2(mRad) ) ); |
---|
1808 | pEmt = Vec4( -pRad.px(), -pRad.py(), pzEmt, |
---|
1809 | sqrt( pow2(pTcorr) + pow2(pzEmt) + pow2(mEmt) ) ); |
---|
1810 | pRec = Vec4( 0., 0., -pzRadPlusEmt, sqrt( pow2(pzRadPlusEmt) |
---|
1811 | + dipSel->m2Rec ) ); |
---|
1812 | |
---|
1813 | // Rotate and boost dipole products to the event frame. |
---|
1814 | pRad.rotbst(M); |
---|
1815 | pEmt.rotbst(M); |
---|
1816 | pRec.rotbst(M); |
---|
1817 | |
---|
1818 | // Azimuthal phi weighting: loop to new phi value if required. |
---|
1819 | if (dipSel->asymPol != 0.) { |
---|
1820 | Vec4 pAunt = event[dipSel->iAunt].p(); |
---|
1821 | double cosPhi = cosphi( pRad, pAunt, pRadBef ); |
---|
1822 | wtPhi = ( 1. + dipSel->asymPol * (2. * pow2(cosPhi) - 1.) ) |
---|
1823 | / ( 1. + abs(dipSel->asymPol) ); |
---|
1824 | } |
---|
1825 | } while (wtPhi < rndmPtr->flat()) ; |
---|
1826 | |
---|
1827 | // Kinematics when recoiler is initial-state parton. |
---|
1828 | int isrTypeNow = dipSel->isrType; |
---|
1829 | int isrTypeSave = isrTypeNow; |
---|
1830 | if (!useLocalRecoilNow) isrTypeNow = 0; |
---|
1831 | if (isrTypeNow != 0) pRec = 2. * recBef.p() - pRec; |
---|
1832 | |
---|
1833 | // PS dec 2010: check if radiator has flexible normalization |
---|
1834 | bool isFlexible = dipSel->isFlexible; |
---|
1835 | |
---|
1836 | // Define new particles from dipole branching. |
---|
1837 | double pTsel = sqrt(dipSel->pT2); |
---|
1838 | Particle rad = Particle(idRad, 51, iRadBef, 0, 0, 0, |
---|
1839 | colRad, acolRad, pRad, mRad, pTsel); |
---|
1840 | Particle emt = Particle(idEmt, 51, iRadBef, 0, 0, 0, |
---|
1841 | colEmt, acolEmt, pEmt, mEmt, pTsel); |
---|
1842 | |
---|
1843 | // Recoiler either in final or in initial state |
---|
1844 | Particle rec = (isrTypeNow == 0) |
---|
1845 | ? Particle(recBef.id(), 52, iRecBef, iRecBef, 0, 0, |
---|
1846 | recBef.col(), recBef.acol(), pRec, dipSel->mRec, pTsel) |
---|
1847 | : Particle(recBef.id(), -53, 0, 0, iRecBef, iRecBef, |
---|
1848 | recBef.col(), recBef.acol(), pRec, 0., 0.); |
---|
1849 | |
---|
1850 | // ME corrections can lead to branching being rejected. |
---|
1851 | if (dipSel->MEtype > 0) { |
---|
1852 | Particle& partner = (dipSel->iMEpartner == iRecBef) |
---|
1853 | ? rec : event[dipSel->iMEpartner]; |
---|
1854 | if ( findMEcorr( dipSel, rad, partner, emt) < rndmPtr->flat() ) |
---|
1855 | return false; |
---|
1856 | } |
---|
1857 | |
---|
1858 | // Rescatter: if the recoiling partner is not in the same system |
---|
1859 | // as the radiator, fix up intermediate systems (can lead |
---|
1860 | // to emissions being vetoed) |
---|
1861 | if (allowRescatter && FIXRESCATTER && isInterleaved |
---|
1862 | && iSysSel != iSysSelRec) { |
---|
1863 | Vec4 pNew = rad.p() + emt.p(); |
---|
1864 | if (!rescatterPropagateRecoil(event, pNew)) return false; |
---|
1865 | } |
---|
1866 | |
---|
1867 | // Save properties to be restored in case of user-hook veto of emission. |
---|
1868 | int eventSizeOld = event.size(); |
---|
1869 | int iRadStatusV = event[iRadBef].status(); |
---|
1870 | int iRadDau1V = event[iRadBef].daughter1(); |
---|
1871 | int iRadDau2V = event[iRadBef].daughter2(); |
---|
1872 | int iRecStatusV = event[iRecBef].status(); |
---|
1873 | int iRecMot1V = event[iRecBef].mother1(); |
---|
1874 | int iRecMot2V = event[iRecBef].mother2(); |
---|
1875 | int iRecDau1V = event[iRecBef].daughter1(); |
---|
1876 | int iRecDau2V = event[iRecBef].daughter2(); |
---|
1877 | int beamOff1 = 1 + beamOffset; |
---|
1878 | int beamOff2 = 2 + beamOffset; |
---|
1879 | int ev1Dau1V = event[beamOff1].daughter1(); |
---|
1880 | int ev2Dau1V = event[beamOff2].daughter1(); |
---|
1881 | |
---|
1882 | // Shower may occur at a displaced vertex. |
---|
1883 | if (radBef.hasVertex()) { |
---|
1884 | rad.vProd( radBef.vProd() ); |
---|
1885 | emt.vProd( radBef.vProd() ); |
---|
1886 | } |
---|
1887 | if (recBef.hasVertex()) rec.vProd( recBef.vProd() ); |
---|
1888 | |
---|
1889 | // Put new particles into the event record. |
---|
1890 | // Mark original dipole partons as branched and set daughters/mothers. |
---|
1891 | int iRad = event.append(rad); |
---|
1892 | int iEmt = event.append(emt); |
---|
1893 | event[iRadBef].statusNeg(); |
---|
1894 | event[iRadBef].daughters( iRad, iEmt); |
---|
1895 | int iRec = 0; |
---|
1896 | if (useLocalRecoilNow) { |
---|
1897 | iRec = event.append(rec); |
---|
1898 | if (isrTypeNow == 0) { |
---|
1899 | event[iRecBef].statusNeg(); |
---|
1900 | event[iRecBef].daughters( iRec, iRec); |
---|
1901 | } else { |
---|
1902 | event[iRecBef].mothers( iRec, iRec); |
---|
1903 | event[iRec].mothers( iRecMot1V, iRecMot2V); |
---|
1904 | if (iRecMot1V == beamOff1) event[beamOff1].daughter1( iRec); |
---|
1905 | if (iRecMot1V == beamOff2) event[beamOff2].daughter1( iRec); |
---|
1906 | } |
---|
1907 | |
---|
1908 | // Global recoil: need to find relevant rotation+boost for recoilers: |
---|
1909 | // boost+rotate to rest frame, boost along z axis, rotate+boost back. |
---|
1910 | } else { |
---|
1911 | RotBstMatrix MG = M; |
---|
1912 | MG.invert(); |
---|
1913 | double pzRecBef = -0.5 * sqrtpos( pow2(dipSel->m2Dip - dipSel->m2Rad |
---|
1914 | - dipSel->m2Rec) - 4. * dipSel->m2Rad * dipSel->m2Rec ) / dipSel->mDip; |
---|
1915 | double eRecBef = sqrt( pow2(pzRecBef) + dipSel->m2Rec); |
---|
1916 | double pzRecAft = -pzRadPlusEmt; |
---|
1917 | double eRecAft = sqrt( pow2(pzRecAft) + dipSel->m2Rec); |
---|
1918 | MG.bst( Vec4(0., 0., pzRecBef, eRecBef), Vec4(0., 0., pzRecAft, eRecAft) ); |
---|
1919 | MG.rotbst( M); |
---|
1920 | |
---|
1921 | // Global recoil: copy particles, and rotate+boost momenta (not vertices). |
---|
1922 | for (int iG = 0; iG < int(iGRecBef.size()); ++iG) { |
---|
1923 | iRec = event.copy( iGRecBef[iG], 52); |
---|
1924 | iGRec.push_back( iRec); |
---|
1925 | Vec4 pGRec = event[iRec].p(); |
---|
1926 | pGRec.rotbst( MG); |
---|
1927 | event[iRec].p( pGRec); |
---|
1928 | } |
---|
1929 | } |
---|
1930 | |
---|
1931 | // Allow veto of branching. If so restore event record to before emission. |
---|
1932 | bool inResonance = (partonSystemsPtr->getInA(iSysSel) == 0) ? true : false; |
---|
1933 | if ( canVetoEmission && userHooksPtr->doVetoFSREmission( eventSizeOld, |
---|
1934 | event, iSysSel, inResonance) ) { |
---|
1935 | event.popBack( event.size() - eventSizeOld); |
---|
1936 | event[iRadBef].status( iRadStatusV); |
---|
1937 | event[iRadBef].daughters( iRadDau1V, iRadDau2V); |
---|
1938 | if (useLocalRecoilNow && isrTypeNow == 0) { |
---|
1939 | event[iRecBef].status( iRecStatusV); |
---|
1940 | event[iRecBef].daughters( iRecDau1V, iRecDau2V); |
---|
1941 | } else if (useLocalRecoilNow) { |
---|
1942 | event[iRecBef].mothers( iRecMot1V, iRecMot2V); |
---|
1943 | if (iRecMot1V == beamOff1) event[beamOff1].daughter1( ev1Dau1V); |
---|
1944 | if (iRecMot1V == beamOff2) event[beamOff2].daughter1( ev2Dau1V); |
---|
1945 | } else { |
---|
1946 | for (int iG = 0; iG < int(iGRecBef.size()); ++iG) { |
---|
1947 | event[iGRecBef[iG]].statusPos(); |
---|
1948 | event[iGRecBef[iG]].daughters( 0, 0); |
---|
1949 | } |
---|
1950 | } |
---|
1951 | return false; |
---|
1952 | } |
---|
1953 | |
---|
1954 | // For global recoil restore the one nominal recoiler, for bookkeeping. |
---|
1955 | if (!useLocalRecoilNow) { |
---|
1956 | iRec = iRecBef; |
---|
1957 | for (int iG = 0; iG < int(iGRecBef.size()); ++iG) |
---|
1958 | if (iGRecBef[iG] == iRecBef) iRec = iGRec[iG]; |
---|
1959 | } |
---|
1960 | |
---|
1961 | // For initial-state recoiler also update beam and sHat info. |
---|
1962 | if (isrTypeNow != 0) { |
---|
1963 | BeamParticle& beamRec = (isrTypeNow == 1) ? *beamAPtr : *beamBPtr; |
---|
1964 | double xOld = beamRec[iSysSelRec].x(); |
---|
1965 | double xRec = 2. * pRec.e() / (beamAPtr->e() + beamBPtr->e()); |
---|
1966 | beamRec[iSysSelRec].iPos( iRec); |
---|
1967 | beamRec[iSysSelRec].x( xRec); |
---|
1968 | partonSystemsPtr->setSHat( iSysSelRec, |
---|
1969 | partonSystemsPtr->getSHat(iSysSelRec) * xRec / xOld); |
---|
1970 | } |
---|
1971 | |
---|
1972 | // Photon emission: update to new dipole ends; add new photon "dipole". |
---|
1973 | if (dipSel->flavour == 22) { |
---|
1974 | dipSel->iRadiator = iRad; |
---|
1975 | dipSel->iRecoiler = iRec; |
---|
1976 | // When recoiler was uncharged particle, in resonance decays, |
---|
1977 | // assign recoil to emitted photons. |
---|
1978 | if (recoilToColoured && inResonance && event[iRec].chargeType() == 0) |
---|
1979 | dipSel->iRecoiler = iEmt; |
---|
1980 | dipSel->pTmax = pTsel; |
---|
1981 | if (doQEDshowerByGamma) dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, |
---|
1982 | pTsel, 0, 0, 1, 0, iSysSel, 0)); |
---|
1983 | |
---|
1984 | // Gluon emission: update both dipole ends and add two new ones. |
---|
1985 | } else if (dipSel->flavour == 21) { |
---|
1986 | dipSel->iRadiator = iRad; |
---|
1987 | dipSel->iRecoiler = iEmt; |
---|
1988 | dipSel->systemRec = iSysSel; |
---|
1989 | dipSel->isrType = 0; |
---|
1990 | dipSel->pTmax = pTsel; |
---|
1991 | // Optionally also kill ME corrections after first emission. |
---|
1992 | if (!doMEafterFirst) dipSel->MEtype = 0; |
---|
1993 | // PS dec 2010: check normalization of radiating dipole |
---|
1994 | // Dipole corresponding to the newly created color tag has normal strength |
---|
1995 | double flexFactor = (isFlexible) ? dipSel->flexFactor : 1.0; |
---|
1996 | dipSel->isFlexible = false; |
---|
1997 | for (int i = 0; i < int(dipEnd.size()); ++i) { |
---|
1998 | if (dipEnd[i].iRadiator == iRecBef && dipEnd[i].iRecoiler == iRadBef |
---|
1999 | && dipEnd[i].colType != 0) { |
---|
2000 | dipEnd[i].iRadiator = iRec; |
---|
2001 | dipEnd[i].iRecoiler = iEmt; |
---|
2002 | // Optionally also kill ME corrections after first emission. |
---|
2003 | if (!doMEafterFirst) dipEnd[i].MEtype = 0; |
---|
2004 | // Strive to match colour to anticolour inside closed system. |
---|
2005 | if ( !isFlexible && dipEnd[i].colType * dipSel->colType > 0) |
---|
2006 | dipEnd[i].iRecoiler = iRad; |
---|
2007 | dipEnd[i].pTmax = pTsel; |
---|
2008 | // PS dec 2010: if the (iRadBef,iRecBef) dipole was flexible, the |
---|
2009 | // same should be true for this (opposite) end. If so, this end keeps |
---|
2010 | // the modified normalization, so we shouldn't need to do anything. |
---|
2011 | } |
---|
2012 | } |
---|
2013 | int colType = (dipSel->colType > 0) ? 2 : -2 ; |
---|
2014 | // When recoiler was uncoloured particle, in resonance decays, |
---|
2015 | // assign recoil to coloured particle. |
---|
2016 | int iRecMod = iRec; |
---|
2017 | if (recoilToColoured && inResonance && event[iRec].col() == 0 |
---|
2018 | && event[iRec].acol() == 0) iRecMod = iRad; |
---|
2019 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRecMod, pTsel, |
---|
2020 | colType, 0, 0, isrTypeSave, iSysSel, 0)); |
---|
2021 | dipEnd.back().systemRec = iSysSelRec; |
---|
2022 | // PS dec 2010: the (iEmt,iRec) dipole "inherits" flexible normalization |
---|
2023 | if (isFlexible) { |
---|
2024 | dipEnd.back().isFlexible = true; |
---|
2025 | dipEnd.back().flexFactor = flexFactor; |
---|
2026 | } |
---|
2027 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, |
---|
2028 | -colType, 0, 0, 0, iSysSel, 0)); |
---|
2029 | |
---|
2030 | // Gluon branching to q qbar: update current dipole and other of gluon. |
---|
2031 | } else if (dipSel->colType != 0) { |
---|
2032 | for (int i = 0; i < int(dipEnd.size()); ++i) { |
---|
2033 | // Strive to match colour to anticolour inside closed system. |
---|
2034 | if ( !isFlexible && dipEnd[i].iRecoiler == iRadBef |
---|
2035 | && dipEnd[i].colType * dipSel->colType < 0 ) |
---|
2036 | dipEnd[i].iRecoiler = iEmt; |
---|
2037 | if (dipEnd[i].iRadiator == iRadBef && abs(dipEnd[i].colType) == 2) { |
---|
2038 | dipEnd[i].colType /= 2; |
---|
2039 | if (dipEnd[i].system != dipEnd[i].systemRec) continue; |
---|
2040 | |
---|
2041 | // Note: gluino -> quark + squark gives a deeper radiation dip than |
---|
2042 | // the more obvious alternative photon decay, so is more realistic. |
---|
2043 | dipEnd[i].MEtype = 66; |
---|
2044 | if (&dipEnd[i] == dipSel) dipEnd[i].iMEpartner = iRad; |
---|
2045 | else dipEnd[i].iMEpartner = iEmt; |
---|
2046 | } |
---|
2047 | } |
---|
2048 | dipSel->iRadiator = iEmt; |
---|
2049 | dipSel->iRecoiler = iRec; |
---|
2050 | dipSel->pTmax = pTsel; |
---|
2051 | |
---|
2052 | // Gluon branching to q qbar: also add two charge dipole ends. |
---|
2053 | // Note: gluino -> quark + squark gives a deeper radiation dip than |
---|
2054 | // the more obvious alternative photon decay, so is more realistic. |
---|
2055 | if (doQEDshowerByQ) { |
---|
2056 | int chgType = event[iRad].chargeType(); |
---|
2057 | dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel, |
---|
2058 | 0, chgType, 0, 0, iSysSel, 66, iEmt)); |
---|
2059 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, |
---|
2060 | 0, -chgType, 0, 0, iSysSel, 66, iRad)); |
---|
2061 | } |
---|
2062 | |
---|
2063 | // Photon branching to f fbar: inactivate photon "dipole"; |
---|
2064 | // optionally add new charge and colour dipole ends. |
---|
2065 | } else if (dipSel->gamType != 0) { |
---|
2066 | dipSel->gamType = 0; |
---|
2067 | int chgType = event[iRad].chargeType(); |
---|
2068 | int colType = event[iRad].colType(); |
---|
2069 | // MEtype = 102 for charge in vector decay. |
---|
2070 | if ( chgType != 0 && ( ( doQEDshowerByQ && colType != 0 ) |
---|
2071 | || ( doQEDshowerByL && colType == 0 ) ) ) { |
---|
2072 | dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel, |
---|
2073 | 0, chgType, 0, 0, iSysSel, 102, iEmt)); |
---|
2074 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, |
---|
2075 | 0, -chgType, 0, 0, iSysSel, 102, iRad)); |
---|
2076 | } |
---|
2077 | // MEtype = 11 for colour in vector decay. |
---|
2078 | if (colType != 0 && doQCDshower) { |
---|
2079 | dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel, |
---|
2080 | colType, 0, 0, 0, iSysSel, 11, iEmt)); |
---|
2081 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, |
---|
2082 | -colType, 0, 0, 0, iSysSel, 11, iRad)); |
---|
2083 | } |
---|
2084 | |
---|
2085 | // Photon_HV emission: update to new dipole ends. |
---|
2086 | } else if (dipSel->flavour == 4900022) { |
---|
2087 | dipSel->iRadiator = iRad; |
---|
2088 | dipSel->iRecoiler = iRec; |
---|
2089 | dipSel->pTmax = pTsel; |
---|
2090 | |
---|
2091 | // Gluon_HV emission: update to new dipole ends. |
---|
2092 | } else if (dipSel->flavour == 4900021) { |
---|
2093 | dipSel->iRadiator = iRad; |
---|
2094 | dipSel->iRecoiler = iEmt; |
---|
2095 | dipSel->pTmax = pTsel; |
---|
2096 | for (int i = 0; i < int(dipEnd.size()); ++i) |
---|
2097 | if (dipEnd[i].iRadiator == iRecBef && dipEnd[i].iRecoiler == iRadBef |
---|
2098 | && dipEnd[i].isHiddenValley) { |
---|
2099 | dipEnd[i].iRadiator = iRec; |
---|
2100 | dipEnd[i].iRecoiler = iEmt; |
---|
2101 | dipEnd[i].pTmax = pTsel; |
---|
2102 | } |
---|
2103 | int colvType = (dipSel->colvType > 0) ? 2 : -2 ; |
---|
2104 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRec, pTsel, |
---|
2105 | 0, 0, 0, isrTypeSave, iSysSel, 0, -1, false, true, colvType) ); |
---|
2106 | dipEnd.back().systemRec = iSysSelRec; |
---|
2107 | dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, |
---|
2108 | 0, 0, 0, 0, iSysSel, 0, -1, false, true, -colvType) ); |
---|
2109 | } |
---|
2110 | |
---|
2111 | // Copy or set lifetime for new final state. |
---|
2112 | if (event[iRad].id() == event[iRadBef].id()) { |
---|
2113 | event[iRad].tau( event[iRadBef].tau() ); |
---|
2114 | } else { |
---|
2115 | event[iRad].tau( event[iRad].tau0() * rndmPtr->exp() ); |
---|
2116 | event[iEmt].tau( event[iEmt].tau0() * rndmPtr->exp() ); |
---|
2117 | } |
---|
2118 | event[iRec].tau( event[iRecBef].tau() ); |
---|
2119 | |
---|
2120 | // Now update other dipoles that also involved the radiator or recoiler. |
---|
2121 | for (int i = 0; i < int(dipEnd.size()); ++i) { |
---|
2122 | // PS dec 2010: if radiator was flexible and now is normal, there may |
---|
2123 | // be other flexible dipoles that need updating. |
---|
2124 | if (isFlexible && !dipSel->isFlexible && dipEnd[i].isFlexible) { |
---|
2125 | if (dipEnd[i].iRecoiler == iRadBef) dipEnd[i].iRecoiler = iEmt; |
---|
2126 | if (dipEnd[i].iRadiator == iRadBef) { |
---|
2127 | dipEnd[i].iRadiator = iEmt; |
---|
2128 | if (dipEnd[i].colType == 1 && dipSel->flavour == 21) |
---|
2129 | dipEnd[i].colType = 2; |
---|
2130 | if (dipEnd[i].colType ==-1 && dipSel->flavour == 21) |
---|
2131 | dipEnd[i].colType =-2; |
---|
2132 | } |
---|
2133 | } |
---|
2134 | if (dipEnd[i].iRadiator == iRadBef) dipEnd[i].iRadiator = iRad; |
---|
2135 | if (dipEnd[i].iRecoiler == iRadBef) dipEnd[i].iRecoiler = iRad; |
---|
2136 | if (dipEnd[i].iMEpartner == iRadBef) dipEnd[i].iMEpartner = iRad; |
---|
2137 | if (useLocalRecoilNow) { |
---|
2138 | if (dipEnd[i].iRadiator == iRecBef) dipEnd[i].iRadiator = iRec; |
---|
2139 | if (dipEnd[i].iRecoiler == iRecBef) dipEnd[i].iRecoiler = iRec; |
---|
2140 | if (dipEnd[i].iMEpartner == iRecBef) dipEnd[i].iMEpartner = iRec; |
---|
2141 | } else { |
---|
2142 | for (int iG = 0; iG < int(iGRecBef.size()); ++iG) { |
---|
2143 | if (dipEnd[i].iRadiator == iGRecBef[iG]) |
---|
2144 | dipEnd[i].iRadiator = iGRec[iG]; |
---|
2145 | if (dipEnd[i].iRecoiler == iGRecBef[iG]) |
---|
2146 | dipEnd[i].iRecoiler = iGRec[iG]; |
---|
2147 | if (dipEnd[i].iMEpartner == iGRecBef[iG]) |
---|
2148 | dipEnd[i].iMEpartner = iGRec[iG]; |
---|
2149 | } |
---|
2150 | } |
---|
2151 | } |
---|
2152 | |
---|
2153 | // PS Apr 2011 |
---|
2154 | // Update any junctions downstream of this branching (if necessary) |
---|
2155 | // (This happens, e.g., via LHEF, when adding showers to intermediate |
---|
2156 | // coloured resonances whose decays involved junctions) |
---|
2157 | for (int iJun = 0; iJun < event.sizeJunction(); iJun++) { |
---|
2158 | // Number of incoming colour lines for this junction. |
---|
2159 | int nIncoming = (event.kindJunction(iJun)-1)/2; |
---|
2160 | // Check radiator colour or anticolour, depending on junction kind |
---|
2161 | // (if junction, incoming = anticolours, and vice versa) |
---|
2162 | int colChk = 0; |
---|
2163 | colChk = ( event.kindJunction(iJun) % 2 == 0 ) |
---|
2164 | ? event[iRadBef].col() : event[iRadBef].acol(); |
---|
2165 | // Loop over incoming junction ends |
---|
2166 | for (int iCol = 0; iCol < nIncoming; iCol++) { |
---|
2167 | int colJun = event.colJunction( iJun, iCol); |
---|
2168 | // If match, update junction end with new upstream (anti)colour |
---|
2169 | if (colJun == colChk) { |
---|
2170 | int colNew = 0; |
---|
2171 | if ( event.kindJunction(iJun) % 2 == 0 ) colNew = colRad; |
---|
2172 | else colNew = acolRad; |
---|
2173 | event.colJunction( iJun, iCol, colNew ); |
---|
2174 | } |
---|
2175 | } |
---|
2176 | } |
---|
2177 | |
---|
2178 | // Finally update the list of all partons in all systems. |
---|
2179 | partonSystemsPtr->replace(iSysSel, iRadBef, iRad); |
---|
2180 | partonSystemsPtr->addOut(iSysSel, iEmt); |
---|
2181 | if (useLocalRecoilNow) |
---|
2182 | partonSystemsPtr->replace(iSysSelRec, iRecBef, iRec); |
---|
2183 | else { |
---|
2184 | for (int iG = 0; iG < int(iGRecBef.size()); ++iG) |
---|
2185 | partonSystemsPtr->replace(iSysSel, iGRecBef[iG], iGRec[iG]); |
---|
2186 | } |
---|
2187 | |
---|
2188 | // Done. |
---|
2189 | return true; |
---|
2190 | |
---|
2191 | } |
---|
2192 | |
---|
2193 | //-------------------------------------------------------------------------- |
---|
2194 | |
---|
2195 | // Rescatter: If a dipole stretches between two different systems, those |
---|
2196 | // systems will no longer locally conserve momentum. These |
---|
2197 | // imbalances become problematic when ISR or primordial kT |
---|
2198 | // is switched on as these steps involve Lorentz boosts. |
---|
2199 | // |
---|
2200 | // 'rescatterPropagateRecoil' tries to fix momentum in all |
---|
2201 | // systems by propogating recoil momentum through all |
---|
2202 | // intermediate systems. As the momentum transfer is already |
---|
2203 | // defined, this can lead to internal lines gaining a |
---|
2204 | // virtuality. |
---|
2205 | |
---|
2206 | // Useful definitions for a pair of integers and a vector of pairs |
---|
2207 | typedef pair < int, int > pairInt; |
---|
2208 | typedef vector < pairInt > vectorPairInt; |
---|
2209 | |
---|
2210 | //-------------------------------------------------------------------------- |
---|
2211 | |
---|
2212 | // findParentSystems |
---|
2213 | // Utility routine to find all parent systems of a given system |
---|
2214 | // Returns a vector of pairs of integers with: |
---|
2215 | // a) The system index, including the starting system (negative |
---|
2216 | // if (b) points to a parent system, positive if (b) points |
---|
2217 | // to a daughter system |
---|
2218 | // b) The event record index that is the path out of the system |
---|
2219 | // (if forwards == false, this is an incoming parton to the |
---|
2220 | // system, and is +ve if side A or -ve if side B, |
---|
2221 | // if forwards == true, this is an outgoing parton from the |
---|
2222 | // system). |
---|
2223 | // Returns as empty vector on failure |
---|
2224 | // Note: this assumes single rescattering only and therefore only |
---|
2225 | // one possible parent system |
---|
2226 | |
---|
2227 | inline vectorPairInt findParentSystems(const int sys, |
---|
2228 | Event& event, PartonSystems* partonSystemsPtr, bool forwards) { |
---|
2229 | |
---|
2230 | vectorPairInt parentSystems; |
---|
2231 | parentSystems.reserve(10); |
---|
2232 | |
---|
2233 | int iSysCur = sys; |
---|
2234 | while (true) { |
---|
2235 | // Get two incoming partons |
---|
2236 | int iInA = partonSystemsPtr->getInA(iSysCur); |
---|
2237 | int iInB = partonSystemsPtr->getInB(iSysCur); |
---|
2238 | |
---|
2239 | // Check if either of these links to another system |
---|
2240 | int iIn = 0; |
---|
2241 | if (event[iInA].isRescatteredIncoming()) iIn = iInA; |
---|
2242 | if (event[iInB].isRescatteredIncoming()) iIn = -iInB; |
---|
2243 | |
---|
2244 | // Save the current system to the vector |
---|
2245 | parentSystems.push_back( pairInt(-iSysCur, iIn) ); |
---|
2246 | if (iIn == 0) break; |
---|
2247 | |
---|
2248 | int iInAbs = abs(iIn); |
---|
2249 | int iMother = event[iInAbs].mother1(); |
---|
2250 | iSysCur = partonSystemsPtr->getSystemOf(iMother); |
---|
2251 | if (iSysCur == -1) { |
---|
2252 | parentSystems.clear(); |
---|
2253 | break; |
---|
2254 | } |
---|
2255 | } // while (true) |
---|
2256 | |
---|
2257 | // If forwards is set, change all event record indices to go to daughter |
---|
2258 | // systems rather than parent systems |
---|
2259 | if (forwards) { |
---|
2260 | vectorPairInt::reverse_iterator rit; |
---|
2261 | for (rit = parentSystems.rbegin(); rit < (parentSystems.rend() - 1); |
---|
2262 | ++rit) { |
---|
2263 | pairInt &cur = *rit; |
---|
2264 | pairInt &next = *(rit + 1); |
---|
2265 | cur.first = -cur.first; |
---|
2266 | cur.second = (next.second < 0) ? -event[abs(next.second)].mother1() : |
---|
2267 | event[abs(next.second)].mother1(); |
---|
2268 | } |
---|
2269 | } |
---|
2270 | |
---|
2271 | return parentSystems; |
---|
2272 | } |
---|
2273 | |
---|
2274 | //-------------------------------------------------------------------------- |
---|
2275 | |
---|
2276 | // rescatterPropagateRecoil |
---|
2277 | // Fix up momentum in all intermediate systems when radiator and recoiler |
---|
2278 | // systems are different. The strategy is to look at all parent systems |
---|
2279 | // from the radiator system and the recoiler system and find where they |
---|
2280 | // intersect. |
---|
2281 | |
---|
2282 | bool TimeShower::rescatterPropagateRecoil( Event& event, Vec4& pNew) { |
---|
2283 | |
---|
2284 | // Some useful variables for later |
---|
2285 | int iRadBef = dipSel->iRadiator; |
---|
2286 | iSysSel = dipSel->system; |
---|
2287 | int iSysSelRec = dipSel->systemRec; |
---|
2288 | Vec4 pImbal = pNew - event[iRadBef].p(); |
---|
2289 | |
---|
2290 | // Store changes locally at first in case we veto the branching |
---|
2291 | // eventMod stores index into the event record and new 4-vector |
---|
2292 | vector < pair < int, Vec4 > > eventMod; |
---|
2293 | eventMod.reserve(10); |
---|
2294 | // systemMod stores system index (iSys) and system-parton index (iMem) |
---|
2295 | // iMem >= 0 - index into outgoing partons (iOut) |
---|
2296 | // iMem == -1 - incoming A |
---|
2297 | // iMem == -2 - incoming B |
---|
2298 | vectorPairInt systemMod; |
---|
2299 | systemMod.reserve(10); |
---|
2300 | |
---|
2301 | // Find all parent systems from radiating and recoiling systems |
---|
2302 | vectorPairInt radParent = findParentSystems(iSysSel, event, |
---|
2303 | partonSystemsPtr, false); |
---|
2304 | vectorPairInt recParent = findParentSystems(iSysSelRec, event, |
---|
2305 | partonSystemsPtr, true); |
---|
2306 | if (radParent.size() == 0 || recParent.size() == 0) { |
---|
2307 | // This should never happen |
---|
2308 | infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " |
---|
2309 | "couldn't find parent system; branching vetoed"); |
---|
2310 | return false; |
---|
2311 | } |
---|
2312 | // Find the system that connects radiating and recoiling system |
---|
2313 | bool foundPath = false; |
---|
2314 | unsigned int iRadP = 0; |
---|
2315 | unsigned int iRecP = 0; |
---|
2316 | for (iRadP = 0; iRadP < radParent.size(); iRadP++) { |
---|
2317 | for (iRecP = 0; iRecP < recParent.size(); iRecP++) |
---|
2318 | if (abs(radParent[iRadP].first) == abs(recParent[iRecP].first)) { |
---|
2319 | foundPath = true; |
---|
2320 | break; |
---|
2321 | } |
---|
2322 | if (foundPath) break; |
---|
2323 | } |
---|
2324 | if (!foundPath) { |
---|
2325 | // Can fail e.g. for QED dipoles where there is no connection |
---|
2326 | // between radiator and recoiler systems |
---|
2327 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2328 | "couldn't find recoil path; branching vetoed"); |
---|
2329 | return false; |
---|
2330 | } |
---|
2331 | |
---|
2332 | // Join together to form complete path from radiating system |
---|
2333 | // to recoiling system |
---|
2334 | vectorPairInt path; |
---|
2335 | if (radParent.size() > 1) |
---|
2336 | path.assign(radParent.begin(), radParent.begin() + iRadP); |
---|
2337 | if (recParent.size() > 1) |
---|
2338 | path.insert(path.end(), recParent.rend() - iRecP - 1, |
---|
2339 | recParent.rend() - 1); |
---|
2340 | |
---|
2341 | // Follow the path fixing up momenta as we go |
---|
2342 | for (unsigned int i = 0; i < path.size(); i++) { |
---|
2343 | // Line out of the current system |
---|
2344 | bool isIncoming = (path[i].first < 0) ? true : false; |
---|
2345 | int iSysCur = abs(path[i].first); |
---|
2346 | bool isIncomingA = (path[i].second > 0) ? true : false; |
---|
2347 | int iLink = abs(path[i].second); |
---|
2348 | |
---|
2349 | int iMemCur; |
---|
2350 | if (isIncoming) iMemCur = (isIncomingA) ? -1 : -2; |
---|
2351 | else { |
---|
2352 | iMemCur = -1; |
---|
2353 | for (int j = 0; j < partonSystemsPtr->sizeOut(iSysCur); j++) |
---|
2354 | if (partonSystemsPtr->getOut(iSysCur, j) == iLink) { |
---|
2355 | iMemCur = j; |
---|
2356 | break; |
---|
2357 | } |
---|
2358 | if (iMemCur == -1) { |
---|
2359 | // This should never happen |
---|
2360 | infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " |
---|
2361 | "couldn't find parton system; branching vetoed"); |
---|
2362 | return false; |
---|
2363 | } |
---|
2364 | } |
---|
2365 | |
---|
2366 | Vec4 pMod = (isIncoming) ? event[iLink].p() + pImbal : |
---|
2367 | event[iLink].p() - pImbal; |
---|
2368 | eventMod.push_back(pair <int, Vec4> (iLink, pMod)); |
---|
2369 | systemMod.push_back(pairInt(iSysCur, iMemCur)); |
---|
2370 | |
---|
2371 | // Calculate sHat of iSysCur |
---|
2372 | int iInCurA = partonSystemsPtr->getInA(iSysCur); |
---|
2373 | int iInCurB = partonSystemsPtr->getInB(iSysCur); |
---|
2374 | Vec4 pTotCur = event[iInCurA].p() + event[iInCurB].p(); |
---|
2375 | |
---|
2376 | // If iMemCur is -1 or -2, then we must have changed the sHat of iSysCur |
---|
2377 | if (iMemCur < 0) pTotCur += (isIncoming) ? pImbal : -pImbal; |
---|
2378 | double sHatCur = pTotCur.m2Calc(); |
---|
2379 | |
---|
2380 | // The fixed-up incoming and outgoing partons should not have |
---|
2381 | // too large a virtuality in relation to the system mass-square. |
---|
2382 | if (abs(pMod.m2Calc()) > MAXVIRTUALITYFRACTION * sHatCur) { |
---|
2383 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2384 | "virtuality much larger than sHat; branching vetoed"); |
---|
2385 | return false; |
---|
2386 | } |
---|
2387 | |
---|
2388 | // Outgoing ones should also not have too large negative energy |
---|
2389 | // in the rest frame of the system. |
---|
2390 | if (!isIncoming && pMod * pTotCur < -MAXNEGENERGYFRACTION * sHatCur) { |
---|
2391 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2392 | "rest frame energy too negative; branching vetoed"); |
---|
2393 | return false; |
---|
2394 | } |
---|
2395 | |
---|
2396 | // Veto negative sHat. |
---|
2397 | if (sHatCur < 0.0) { |
---|
2398 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2399 | "sHat became negative; branching vetoed"); |
---|
2400 | return false; |
---|
2401 | } |
---|
2402 | |
---|
2403 | // Line into the new current system |
---|
2404 | iLink = (isIncoming) ? event[iLink].mother1() : |
---|
2405 | event[iLink].daughter1(); |
---|
2406 | iSysCur = partonSystemsPtr->getSystemOf(iLink, true); |
---|
2407 | |
---|
2408 | if (!isIncoming) iMemCur = (isIncomingA) ? -1 : -2; |
---|
2409 | else { |
---|
2410 | iMemCur = -1; |
---|
2411 | for (int j = 0; j < partonSystemsPtr->sizeOut(iSysCur); j++) |
---|
2412 | if (partonSystemsPtr->getOut(iSysCur, j) == iLink) { |
---|
2413 | iMemCur = j; |
---|
2414 | break; |
---|
2415 | } |
---|
2416 | if (iMemCur == -1) { |
---|
2417 | // This should never happen |
---|
2418 | infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " |
---|
2419 | "couldn't find parton system; branching vetoed"); |
---|
2420 | return false; |
---|
2421 | } |
---|
2422 | } |
---|
2423 | |
---|
2424 | pMod = (isIncoming) ? event[iLink].p() + pImbal : |
---|
2425 | event[iLink].p() - pImbal; |
---|
2426 | eventMod.push_back(pair <int, Vec4> (iLink, pMod)); |
---|
2427 | systemMod.push_back(pairInt(iSysCur, iMemCur)); |
---|
2428 | |
---|
2429 | // Calculate sHat of iSysCur |
---|
2430 | iInCurA = partonSystemsPtr->getInA(iSysCur); |
---|
2431 | iInCurB = partonSystemsPtr->getInB(iSysCur); |
---|
2432 | pTotCur = event[iInCurA].p() + event[iInCurB].p(); |
---|
2433 | |
---|
2434 | // If iMemCur is -1 or -2, then we must have changed the sHat of iSysCur |
---|
2435 | if (iMemCur < 0) pTotCur += (isIncoming) ? pImbal : -pImbal; |
---|
2436 | sHatCur = pTotCur.m2Calc(); |
---|
2437 | |
---|
2438 | // The fixed-up incoming and outgoing partons should not have |
---|
2439 | // too large a virtuality in relation to the system mass-square. |
---|
2440 | if (abs(pMod.m2Calc()) > MAXVIRTUALITYFRACTION * sHatCur) { |
---|
2441 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2442 | "virtuality much larger than sHat; branching vetoed"); |
---|
2443 | return false; |
---|
2444 | } |
---|
2445 | |
---|
2446 | // Outgoing ones should also not have too large negative energy |
---|
2447 | // in the rest frame of the system. |
---|
2448 | if (!isIncoming && pMod * pTotCur < -MAXNEGENERGYFRACTION * sHatCur) { |
---|
2449 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2450 | "rest frame energy too negative; branching vetoed"); |
---|
2451 | return false; |
---|
2452 | } |
---|
2453 | |
---|
2454 | // Veto negative sHat |
---|
2455 | if (sHatCur < 0.0) { |
---|
2456 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2457 | "sHat became negative; branching vetoed"); |
---|
2458 | return false; |
---|
2459 | } |
---|
2460 | |
---|
2461 | // Do negative energy veto |
---|
2462 | if (VETONEGENERGY && pMod.e() < 0.0) { |
---|
2463 | infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " |
---|
2464 | "energy became negative; branching vetoed"); |
---|
2465 | return false; |
---|
2466 | } |
---|
2467 | |
---|
2468 | } // for (unsigned int i = 0; i < path.size(); i++) |
---|
2469 | |
---|
2470 | // If no vetos by this point, apply the changes to the event record |
---|
2471 | // An incoming particle with changed momentum is given status code -54, |
---|
2472 | // an outgoing particle with changed momentum is given status code -55 |
---|
2473 | for (unsigned int i = 0; i < eventMod.size(); i++) { |
---|
2474 | int idx = eventMod[i].first; |
---|
2475 | Vec4 &pMod = eventMod[i].second; |
---|
2476 | int iSys = systemMod[i].first; |
---|
2477 | int iMem = systemMod[i].second; |
---|
2478 | |
---|
2479 | // If incoming to a process then set the copy to be the mother |
---|
2480 | if (event[idx].isRescatteredIncoming()) { |
---|
2481 | int mother1 = event[idx].mother1(); |
---|
2482 | idx = event.copy(idx, -54); |
---|
2483 | event[mother1].daughters(idx, idx); |
---|
2484 | |
---|
2485 | // Update beam information if necessary |
---|
2486 | double eCM = sqrt(m2( beamAPtr->p(), beamBPtr->p())); |
---|
2487 | if (iMem == -1) { |
---|
2488 | partonSystemsPtr->setInA(iSys, idx); |
---|
2489 | (*beamAPtr)[iSys].x((pMod.e() + pMod.pz()) / eCM); |
---|
2490 | (*beamAPtr)[iSys].m(pMod.mCalc()); |
---|
2491 | (*beamAPtr)[iSys].p(pMod); |
---|
2492 | (*beamAPtr)[iSys].iPos(idx); |
---|
2493 | } else if (iMem == -2) { |
---|
2494 | partonSystemsPtr->setInB(iSys, idx); |
---|
2495 | (*beamBPtr)[iSys].x((pMod.e() - pMod.pz()) / eCM); |
---|
2496 | (*beamBPtr)[iSys].m(pMod.mCalc()); |
---|
2497 | (*beamBPtr)[iSys].p(pMod); |
---|
2498 | (*beamBPtr)[iSys].iPos(idx); |
---|
2499 | } else { |
---|
2500 | // This should never happen |
---|
2501 | infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " |
---|
2502 | "internal bookeeping error"); |
---|
2503 | } |
---|
2504 | |
---|
2505 | // Otherwise set the new event record entry to be the daughter |
---|
2506 | } else { |
---|
2507 | int daughter1 = event[idx].daughter1(); |
---|
2508 | idx = event.copy(idx, 55); |
---|
2509 | event[idx].statusNeg(); |
---|
2510 | event[daughter1].mothers(idx, idx); |
---|
2511 | |
---|
2512 | partonSystemsPtr->setOut(iSys, iMem, idx); |
---|
2513 | } |
---|
2514 | |
---|
2515 | event[idx].p( eventMod[i].second ); |
---|
2516 | event[idx].m( event[idx].mCalc() ); |
---|
2517 | } |
---|
2518 | |
---|
2519 | return true; |
---|
2520 | } |
---|
2521 | |
---|
2522 | |
---|
2523 | //-------------------------------------------------------------------------- |
---|
2524 | |
---|
2525 | // Find class of QCD ME correction. |
---|
2526 | // MEtype classification follow codes in Norrbin article, |
---|
2527 | // additionally -1 = try to find type, 0 = no ME corrections. |
---|
2528 | // Warning: not yet tried out to do a correct assignment in |
---|
2529 | // arbitrary multiparton configurations! ?? |
---|
2530 | |
---|
2531 | void TimeShower::findMEtype( Event& event, TimeDipoleEnd& dip) { |
---|
2532 | |
---|
2533 | // Initial value. Mark if no ME corrections to be applied. |
---|
2534 | bool setME = true; |
---|
2535 | if (!doMEcorrections) setME = false; |
---|
2536 | int iMother = event[dip.iRadiator].mother1(); |
---|
2537 | int iMother2 = event[dip.iRadiator].mother2(); |
---|
2538 | |
---|
2539 | // Allow ME corrections for Hidden Valley pair in 2 -> 2. |
---|
2540 | if (dip.isHiddenValley && event[dip.iRecoiler].id() |
---|
2541 | == -event[dip.iRadiator].id()); |
---|
2542 | |
---|
2543 | // Else no ME corrections in 2 -> n processes. |
---|
2544 | else { |
---|
2545 | if (iMother2 != iMother && iMother2 != 0) setME = false; |
---|
2546 | if (event[dip.iRecoiler].mother1() != iMother) setME = false; |
---|
2547 | if (event[dip.iRecoiler].mother2() != iMother2) setME = false; |
---|
2548 | } |
---|
2549 | |
---|
2550 | // No ME corrections for recoiler in initial state. |
---|
2551 | if (event[dip.iRecoiler].status() < 0) setME = false; |
---|
2552 | |
---|
2553 | // No ME corrections for recoiler not in same system |
---|
2554 | if (dip.system != dip.systemRec) setME = false; |
---|
2555 | |
---|
2556 | // Done if no ME to be set. |
---|
2557 | if (!setME) { |
---|
2558 | dip.MEtype = 0; |
---|
2559 | return; |
---|
2560 | } |
---|
2561 | |
---|
2562 | // If no ME partner set, assume it is the recoiler. |
---|
2563 | if (dip.iMEpartner < 0) dip.iMEpartner = dip.iRecoiler; |
---|
2564 | |
---|
2565 | // Now begin processing of colour dipole, including Hidden Valley. |
---|
2566 | if (dip.colType != 0 || dip.colvType != 0) { |
---|
2567 | bool isHiddenColour = (dip.colvType != 0); |
---|
2568 | |
---|
2569 | // Find daughter types (may or may not be used later on). |
---|
2570 | int idDau1 = event[dip.iRadiator].id(); |
---|
2571 | int idDau2 = event[dip.iMEpartner].id(); |
---|
2572 | int dau1Type = findMEparticle(idDau1, isHiddenColour); |
---|
2573 | int dau2Type = findMEparticle(idDau2, isHiddenColour); |
---|
2574 | int minDauType = min(dau1Type, dau2Type); |
---|
2575 | int maxDauType = max(dau1Type, dau2Type); |
---|
2576 | |
---|
2577 | // Reorder dipole ends in kinematics. Split ME expression in two sides. |
---|
2578 | dip.MEorder = (dau2Type >= dau1Type); |
---|
2579 | dip.MEsplit = (maxDauType <= 6); |
---|
2580 | dip.MEgluinoRec = false; |
---|
2581 | |
---|
2582 | // If type already set (or set not to have) then done. |
---|
2583 | if (minDauType == 0 && dip.MEtype < 0) dip.MEtype = 0; |
---|
2584 | if (dip.MEtype >= 0) return; |
---|
2585 | dip.MEtype = 0; |
---|
2586 | |
---|
2587 | // For H -> gg -> ggg we found that DGLAP kernels do better than eikonal. |
---|
2588 | if (dau1Type == 4 && dau2Type == 4) return; |
---|
2589 | |
---|
2590 | // Find mother type. |
---|
2591 | int idMother = 0; |
---|
2592 | if ( event[dip.iRecoiler].mother1() == iMother && iMother >= 0) |
---|
2593 | idMother = event[iMother].id(); |
---|
2594 | int motherType = (idMother != 0) |
---|
2595 | ? findMEparticle(idMother, isHiddenColour) : 0; |
---|
2596 | |
---|
2597 | // When a mother if not known then use colour and spin content to guess. |
---|
2598 | if (motherType == 0) { |
---|
2599 | int col1 = event[dip.iRadiator].col(); |
---|
2600 | int acol1 = event[dip.iRadiator].acol(); |
---|
2601 | int col2 = event[dip.iMEpartner].col(); |
---|
2602 | int acol2 = event[dip.iMEpartner].acol(); |
---|
2603 | // spinT = 0/1 = integer or half-integer. |
---|
2604 | int spinT = ( event[dip.iRadiator].spinType() |
---|
2605 | + event[dip.iMEpartner].spinType() )%2; |
---|
2606 | // Colour singlet mother. |
---|
2607 | if ( col1 == acol2 && acol1 == col2 ) |
---|
2608 | motherType = (spinT == 0) ? 7 : 9; |
---|
2609 | // Colour octet mother. |
---|
2610 | else if ( (col1 == acol2 && acol1 != 0 && col2 != 0) |
---|
2611 | || (acol1 == col2 && col1 != 0 && acol2 != 0) ) |
---|
2612 | motherType = (spinT == 0) ? 4 : 5; |
---|
2613 | // Colour triplet mother. |
---|
2614 | else if ( (col1 == acol2 && acol1 != col2) |
---|
2615 | || (acol1 == col2 && col1 != acol2) ) |
---|
2616 | motherType = (spinT == 0) ? 2 : 1; |
---|
2617 | // If no colours are matched then cannot have common mother, so done. |
---|
2618 | else return; |
---|
2619 | } |
---|
2620 | |
---|
2621 | // Now start from default, which is eikonal ME corrections, |
---|
2622 | // and try to find matching ME cases below. |
---|
2623 | int MEkind = 0; |
---|
2624 | int MEcombi = 4; |
---|
2625 | dip.MEmix = 0.5; |
---|
2626 | |
---|
2627 | // Hidden Valley with massive gamma_v covered by two special cases. |
---|
2628 | if (isHiddenColour && brokenHVsym) { |
---|
2629 | MEkind = (dau2Type == 0 || dau2Type > 6) ? 30 : 31; |
---|
2630 | dip.MEtype = 5 * MEkind + 1; |
---|
2631 | return; |
---|
2632 | } |
---|
2633 | |
---|
2634 | // Triplet recoiling against gluino needs enhanced radiation |
---|
2635 | // to match to matrix elements. |
---|
2636 | dip.MEgluinoRec = (dau1Type >= 1 && dau1Type <= 3 && dau2Type == 5); |
---|
2637 | |
---|
2638 | // Vector/axial vector -> q + qbar. |
---|
2639 | if (minDauType == 1 && maxDauType == 1 && |
---|
2640 | (motherType == 4 || motherType == 7) ) { |
---|
2641 | MEkind = 2; |
---|
2642 | if (idMother == 21 || idMother == 22) MEcombi = 1; |
---|
2643 | else if (idMother == 23 || idDau1 + idDau2 == 0) { |
---|
2644 | MEcombi = 3; |
---|
2645 | dip.MEmix = gammaZmix( event, iMother, dip.iRadiator, dip.iRecoiler ); |
---|
2646 | } |
---|
2647 | else if (idMother == 24) MEcombi = 4; |
---|
2648 | } |
---|
2649 | // For chi -> chi q qbar, use V/A -> q qbar as first approximation. |
---|
2650 | else if (minDauType == 1 && maxDauType == 1 && motherType == 9) |
---|
2651 | MEkind = 2; |
---|
2652 | |
---|
2653 | // q -> q + V. |
---|
2654 | else if (minDauType == 1 && maxDauType == 7 && motherType == 1) |
---|
2655 | MEkind = 3; |
---|
2656 | if (idDau1 == 22 || idDau2 == 22) MEcombi = 1; |
---|
2657 | |
---|
2658 | // Scalar/pseudoscalar -> q + qbar; q -> q + S. |
---|
2659 | else if (minDauType == 1 && maxDauType == 1 && motherType == 8) { |
---|
2660 | MEkind = 4; |
---|
2661 | if (idMother == 25 || idMother == 35 || idMother == 37) MEcombi = 1; |
---|
2662 | else if (idMother == 36) MEcombi = 2; |
---|
2663 | } |
---|
2664 | else if (minDauType == 1 && maxDauType == 8 && motherType == 1) |
---|
2665 | MEkind = 5; |
---|
2666 | |
---|
2667 | // V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S. |
---|
2668 | else if (minDauType == 2 && maxDauType == 2 && (motherType == 4 |
---|
2669 | || motherType == 7) ) MEkind = 6; |
---|
2670 | else if (minDauType == 2 && (maxDauType == 4 || maxDauType == 7) |
---|
2671 | && motherType == 2) MEkind = 7; |
---|
2672 | else if (minDauType == 2 && maxDauType == 2 && motherType == 8) |
---|
2673 | MEkind = 8; |
---|
2674 | else if (minDauType == 2 && maxDauType == 8 && motherType == 2) |
---|
2675 | MEkind = 9; |
---|
2676 | |
---|
2677 | // chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi. |
---|
2678 | else if (minDauType == 1 && maxDauType == 2 && motherType == 9) |
---|
2679 | MEkind = 10; |
---|
2680 | else if (minDauType == 1 && maxDauType == 9 && motherType == 2) |
---|
2681 | MEkind = 11; |
---|
2682 | else if (minDauType == 2 && maxDauType == 9 && motherType == 1) |
---|
2683 | MEkind = 12; |
---|
2684 | |
---|
2685 | // ~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g. |
---|
2686 | else if (minDauType == 1 && maxDauType == 2 && motherType == 5) |
---|
2687 | MEkind = 13; |
---|
2688 | else if (minDauType == 1 && maxDauType == 5 && motherType == 2) |
---|
2689 | MEkind = 14; |
---|
2690 | else if (minDauType == 2 && maxDauType == 5 && motherType == 1) |
---|
2691 | MEkind = 15; |
---|
2692 | |
---|
2693 | // In cases where coloured spin 1 particle involved use spin 0. |
---|
2694 | // V_coloured -> q + l. |
---|
2695 | else if (minDauType == 1 && maxDauType == 9 && motherType == 3) |
---|
2696 | MEkind = 11; |
---|
2697 | // q -> V_coloured + l; |
---|
2698 | else if (minDauType == 3 && maxDauType == 9 && motherType == 1) |
---|
2699 | MEkind = 12; |
---|
2700 | |
---|
2701 | // g (+V, S) -> ~g + ~g (eikonal approximation). |
---|
2702 | else if (minDauType == 5 && maxDauType == 5) MEkind = 16; |
---|
2703 | |
---|
2704 | // Save ME type and gamma_5 admixture. |
---|
2705 | dip.MEtype = 5 * MEkind + MEcombi; |
---|
2706 | |
---|
2707 | // Now begin processing of charge dipole - still primitive. |
---|
2708 | } else if (dip.chgType != 0) { |
---|
2709 | |
---|
2710 | // Set defaults for QED case; then possibly done. |
---|
2711 | dip.MEorder = true; |
---|
2712 | dip.MEsplit = true; |
---|
2713 | if (dip.MEtype >= 0) return; |
---|
2714 | |
---|
2715 | // So far only ME corrections for q qbar or l lbar. |
---|
2716 | int idDau1 = event[dip.iRadiator].id(); |
---|
2717 | int idDau2 = event[dip.iMEpartner].id(); |
---|
2718 | if (abs(idDau1) < 9 && abs(idDau2) < 9 && idDau1 * idDau2 < 0) ; |
---|
2719 | else if (abs(idDau1) > 10 && abs(idDau1) < 19 && abs(idDau2) > 10 |
---|
2720 | && abs(idDau2) < 19 && idDau1 * idDau2 < 0) ; |
---|
2721 | else { dip.MEtype = 0; return; } |
---|
2722 | |
---|
2723 | // Distinguish charge sum != 0 or = 0; in latter assume vector source. |
---|
2724 | dip.MEtype = 101; |
---|
2725 | if (idDau1 + idDau2 == 0) dip.MEtype = 102; |
---|
2726 | dip.MEmix = 1.; |
---|
2727 | } |
---|
2728 | |
---|
2729 | } |
---|
2730 | |
---|
2731 | //-------------------------------------------------------------------------- |
---|
2732 | |
---|
2733 | // Find type of particle for ME type: 0 = unknown, 1 = quark, 2 = squark, |
---|
2734 | // 3 = spare triplet, 4 = gluon, 5 = gluino, 6 = spare octet, |
---|
2735 | // 7 = vector boson, 8 = colourless scalar, 9 = colourless spin 1/2. |
---|
2736 | |
---|
2737 | int TimeShower::findMEparticle( int id, bool isHiddenColour) { |
---|
2738 | |
---|
2739 | // find colour and spin of particle. |
---|
2740 | int type = 0; |
---|
2741 | int colType = abs(particleDataPtr->colType(id)); |
---|
2742 | int spinType = particleDataPtr->spinType(id); |
---|
2743 | |
---|
2744 | // For hidden valley particle treat HV colour as normal one. |
---|
2745 | // Note: no need to assign gv/gammav since not in ME. |
---|
2746 | if (isHiddenColour) { |
---|
2747 | colType = 0; |
---|
2748 | int idAbs = abs(id); |
---|
2749 | if ( (idAbs > 4900000 && idAbs < 4900007) |
---|
2750 | || (idAbs > 4900010 && idAbs < 4900017) |
---|
2751 | || idAbs == 4900101) colType = 1; |
---|
2752 | } |
---|
2753 | |
---|
2754 | // Find particle type from colour and spin. |
---|
2755 | if (colType == 1 && spinType == 2) type = 1; |
---|
2756 | else if (colType == 1 && spinType == 1) type = 2; |
---|
2757 | else if (colType == 1) type = 3; |
---|
2758 | else if (colType == 2 && spinType == 3) type = 4; |
---|
2759 | else if (colType == 2 && spinType == 2) type = 5; |
---|
2760 | else if (colType == 2) type = 6; |
---|
2761 | else if (colType == 0 && spinType == 3) type = 7; |
---|
2762 | else if (colType == 0 && spinType == 1) type = 8; |
---|
2763 | else if (colType == 0 && spinType == 2) type = 9; |
---|
2764 | |
---|
2765 | // Done. |
---|
2766 | return type; |
---|
2767 | |
---|
2768 | } |
---|
2769 | |
---|
2770 | //-------------------------------------------------------------------------- |
---|
2771 | |
---|
2772 | // Find mixture of V and A in gamma/Z: energy- and flavour-dependent. |
---|
2773 | |
---|
2774 | double TimeShower::gammaZmix( Event& event, int iRes, int iDau1, int iDau2) { |
---|
2775 | |
---|
2776 | // Try to identify initial flavours; use e+e- as default. |
---|
2777 | int idIn1 = -11; |
---|
2778 | int idIn2 = 11; |
---|
2779 | int iIn1 = (iRes >= 0) ? event[iRes].mother1() : -1; |
---|
2780 | int iIn2 = (iRes >= 0) ? event[iRes].mother2() : -1; |
---|
2781 | if (iIn1 >=0) idIn1 = event[iIn1].id(); |
---|
2782 | if (iIn2 >=0) idIn2 = event[iIn1].id(); |
---|
2783 | |
---|
2784 | // In processes f + g/gamma -> f + Z only need find one fermion. |
---|
2785 | if (idIn1 == 21 || idIn1 == 22) idIn1 = -idIn2; |
---|
2786 | if (idIn2 == 21 || idIn2 == 22) idIn2 = -idIn1; |
---|
2787 | |
---|
2788 | // Initial flavours and couplings; return if don't make sense. |
---|
2789 | if (idIn1 + idIn2 != 0 ) return 0.5; |
---|
2790 | int idInAbs = abs(idIn1); |
---|
2791 | if (idInAbs == 0 || idInAbs > 18 ) return 0.5; |
---|
2792 | double ei = coupSMPtr->ef(idInAbs); |
---|
2793 | double vi = coupSMPtr->vf(idInAbs); |
---|
2794 | double ai = coupSMPtr->af(idInAbs); |
---|
2795 | |
---|
2796 | // Final flavours and couplings; return if don't make sense. |
---|
2797 | if (event[iDau1].id() + event[iDau2].id() != 0) return 0.5; |
---|
2798 | int idOutAbs = abs(event[iDau1].id()); |
---|
2799 | if (idOutAbs == 0 || idOutAbs >18 ) return 0.5; |
---|
2800 | double ef = coupSMPtr->ef(idOutAbs); |
---|
2801 | double vf = coupSMPtr->vf(idOutAbs); |
---|
2802 | double af = coupSMPtr->af(idOutAbs); |
---|
2803 | |
---|
2804 | // Calculate prefactors for interference and resonance part. |
---|
2805 | Vec4 psum = event[iDau1].p() + event[iDau2].p(); |
---|
2806 | double sH = psum.m2Calc(); |
---|
2807 | double intNorm = 2. * thetaWRat * sH * (sH - mZ*mZ) |
---|
2808 | / ( pow2(sH - mZ*mZ) + pow2(sH * gammaZ / mZ) ); |
---|
2809 | double resNorm = pow2(thetaWRat * sH) |
---|
2810 | / ( pow2(sH - mZ*mZ) + pow2(sH * gammaZ / mZ) ); |
---|
2811 | |
---|
2812 | // Calculate vector and axial expressions and find mix. |
---|
2813 | double vect = ei*ei * ef*ef + ei*vi * intNorm * ef*vf |
---|
2814 | + (vi*vi + ai*ai) * resNorm * vf*vf; |
---|
2815 | double axiv = (vi*vi + ai*ai) * resNorm * af*af; |
---|
2816 | return vect / (vect + axiv); |
---|
2817 | } |
---|
2818 | |
---|
2819 | //-------------------------------------------------------------------------- |
---|
2820 | |
---|
2821 | // Set up to calculate QCD ME correction with calcMEcorr. |
---|
2822 | // Normally for primary particles, but also from g/gamma -> f fbar. |
---|
2823 | |
---|
2824 | double TimeShower::findMEcorr(TimeDipoleEnd* dip, Particle& rad, |
---|
2825 | Particle& partner, Particle& emt) { |
---|
2826 | |
---|
2827 | // Initial values and matrix element kind. |
---|
2828 | double wtME = 1.; |
---|
2829 | double wtPS = 1.; |
---|
2830 | int MEkind = dip->MEtype / 5; |
---|
2831 | int MEcombi = dip->MEtype % 5; |
---|
2832 | |
---|
2833 | // Construct ME variables. |
---|
2834 | Vec4 sum = rad.p() + partner.p() + emt.p(); |
---|
2835 | double eCMME = sum.mCalc(); |
---|
2836 | double x1 = 2. * (sum * rad.p()) / pow2(eCMME); |
---|
2837 | double x2 = 2. * (sum * partner.p()) / pow2(eCMME); |
---|
2838 | double r1 = rad.m() / eCMME; |
---|
2839 | double r2 = partner.m() / eCMME; |
---|
2840 | double r3 = 0.; |
---|
2841 | |
---|
2842 | // Evaluate kinematics for Hidden Valley with massive gamma_v. |
---|
2843 | double gammavCorr = 1.; |
---|
2844 | if (dip->colvType != 0 && brokenHVsym) { |
---|
2845 | r3 = emt.m() / eCMME; |
---|
2846 | double x3Tmp = 2. - x1 - x2; |
---|
2847 | gammavCorr = x3Tmp / (x3Tmp - kRad * (x1 + x3Tmp)); |
---|
2848 | // For Q_v Qbar_v pair correct kinematics to common average mass. |
---|
2849 | if (MEkind == 31) { |
---|
2850 | double m2Pair = (rad.p() + partner.p()).m2Calc(); |
---|
2851 | double m2Avg = 0.5 * (rad.m2() + partner.m2()) |
---|
2852 | - 0.25 * pow2(rad.m2() - partner.m2()) / m2Pair; |
---|
2853 | r1 = sqrt(m2Avg) / eCMME; |
---|
2854 | r2 = r1; |
---|
2855 | double xShift = 0.5 * (x1 + x2) * (partner.m2() - rad.m2()) / m2Pair; |
---|
2856 | x1 += xShift; |
---|
2857 | x2 -= xShift; |
---|
2858 | } |
---|
2859 | } |
---|
2860 | |
---|
2861 | // Derived ME variables, suitably protected. |
---|
2862 | double x1minus = max(XMARGIN, 1. + r1*r1 - r2*r2 - x1); |
---|
2863 | double x2minus = max(XMARGIN, 1. + r2*r2 - r1*r1 - x2) ; |
---|
2864 | double x3 = max(XMARGIN, 2. - x1 - x2); |
---|
2865 | |
---|
2866 | // Begin processing of QCD dipoles. |
---|
2867 | if (dip->colType !=0 || dip->colvType != 0) { |
---|
2868 | |
---|
2869 | // Evaluate normal ME, for proper order of particles. |
---|
2870 | if (dip->MEorder) |
---|
2871 | wtME = calcMEcorr(MEkind, MEcombi, dip->MEmix, x1, x2, r1, r2, r3); |
---|
2872 | else wtME = calcMEcorr(MEkind, MEcombi, dip->MEmix, x2, x1, r2, r1, r3); |
---|
2873 | |
---|
2874 | // Split up total ME when two radiating particles. |
---|
2875 | if (dip->MEsplit) wtME = wtME * x1minus / x3; |
---|
2876 | |
---|
2877 | // Evaluate shower rate to be compared with. |
---|
2878 | wtPS = 2. / ( x3 * x2minus ); |
---|
2879 | if (dip->MEgluinoRec) wtPS *= 9./4.; |
---|
2880 | if (dip->colvType != 0 && brokenHVsym) wtPS *= gammavCorr; |
---|
2881 | |
---|
2882 | // For generic charge combination currently only massless expression. |
---|
2883 | // (Masses included only to respect phase space boundaries.) |
---|
2884 | } else if (dip->chgType !=0 && dip->MEtype == 101) { |
---|
2885 | double chg1 = particleDataPtr->charge(rad.id()); |
---|
2886 | double chg2 = particleDataPtr->charge(partner.id()); |
---|
2887 | wtME = (x1*x1 + x2*x2) * pow2( chg1 * x1minus / x3 |
---|
2888 | - chg2 * x2minus / x3 ); |
---|
2889 | wtPS = 2. * ( chg1*chg1 * x1minus / x3 + chg2*chg2 * x2minus / x3 ); |
---|
2890 | |
---|
2891 | // For flavour neutral system assume vector source and include masses. |
---|
2892 | } else if (dip->chgType !=0 && dip->MEtype == 102) { |
---|
2893 | wtME = calcMEcorr(2, 1, dip->MEmix, x1, x2, r1, r2) * x1minus / x3; |
---|
2894 | wtPS = 2. / ( x3 * x2minus ); |
---|
2895 | } |
---|
2896 | if (wtME > wtPS) infoPtr->errorMsg("Warning in TimeShower::findMEcorr: " |
---|
2897 | "ME weight above PS one"); |
---|
2898 | |
---|
2899 | // Return ratio of actual ME to assumed PS rate of emission. |
---|
2900 | return wtME / wtPS; |
---|
2901 | } |
---|
2902 | |
---|
2903 | //-------------------------------------------------------------------------- |
---|
2904 | |
---|
2905 | // Matrix elements for gluon (or photon) emission from |
---|
2906 | // a two-body state; to be used by the parton shower routine. |
---|
2907 | // Here x_i = 2 E_i/E_cm, r_i = m_i/E_cm and |
---|
2908 | // 1/sigma_0 d(sigma)/d(x_1)d(x_2) = (alpha-strong/2 pi) * C_F * (this), |
---|
2909 | // i.e. normalization is such that one recovers the familiar |
---|
2910 | // (x_1^2 + x_2^2)/((1-x_1)*(1-x_2)) for the massless case. |
---|
2911 | // Coupling structure: |
---|
2912 | // kind = 1 : eikonal soft-gluon expression (spin-independent) |
---|
2913 | // = 2 : V -> q qbar (V = vector/axial vector colour singlet) |
---|
2914 | // = 3 : q -> q V |
---|
2915 | // = 4 : S -> q qbar (S = scalar/pseudoscalar colour singlet) |
---|
2916 | // = 5 : q -> q S |
---|
2917 | // = 6 : V -> ~q ~qbar (~q = squark) |
---|
2918 | // = 7 : ~q -> ~q V |
---|
2919 | // = 8 : S -> ~q ~qbar |
---|
2920 | // = 9 : ~q -> ~q S |
---|
2921 | // = 10 : chi -> q ~qbar (chi = neutralino/chargino) |
---|
2922 | // = 11 : ~q -> q chi |
---|
2923 | // = 12 : q -> ~q chi |
---|
2924 | // = 13 : ~g -> q ~qbar |
---|
2925 | // = 14 : ~q -> q ~g |
---|
2926 | // = 15 : q -> ~q ~g |
---|
2927 | // = 16 : (9/4)*(eikonal) for gg -> ~g ~g |
---|
2928 | // = 30 : Dv -> d qv (Dv= hidden valley fermion, qv= valley scalar) |
---|
2929 | // = 31 : S -> Dv Dvbar (S=scalar color singlet) |
---|
2930 | // Note that the order of the decay products is important. |
---|
2931 | // combi = 1 : pure non-gamma5, i.e. vector/scalar/... |
---|
2932 | // = 2 : pure gamma5, i.e. axial vector/pseudoscalar/.... |
---|
2933 | // = 3 : mixture mix*(combi=1) + (1-mix)*(combi=2) |
---|
2934 | // = 4 : mixture (combi=1) +- (combi=2) |
---|
2935 | |
---|
2936 | double TimeShower::calcMEcorr( int kind, int combiIn, double mixIn, |
---|
2937 | double x1, double x2, double r1, double r2, double r3) { |
---|
2938 | |
---|
2939 | // Frequent variable combinations. |
---|
2940 | double x3 = 2. - x1 - x2; |
---|
2941 | double x1s = x1 * x1; |
---|
2942 | double x2s = x2 * x2; |
---|
2943 | double x3s = x3 * x3; |
---|
2944 | double x1c = x1 * x1s; |
---|
2945 | double x2c = x2 * x2s; |
---|
2946 | double x3c = x3 * x3s; |
---|
2947 | double r1s = r1 * r1; |
---|
2948 | double r2s = r2 * r2; |
---|
2949 | double r1c = r1 * r1s; |
---|
2950 | double r2c = r2 * r2s; |
---|
2951 | double r1q = r1s * r1s; |
---|
2952 | double r2q = r2s * r2s; |
---|
2953 | double prop1 = 1. + r1s - r2s - x1; |
---|
2954 | double prop2 = 1. + r2s - r1s - x2; |
---|
2955 | double prop1s = prop1 * prop1; |
---|
2956 | double prop2s = prop2 * prop2; |
---|
2957 | double prop12 = prop1 * prop2; |
---|
2958 | double prop13 = prop1 * x3; |
---|
2959 | double prop23 = prop2 * x3; |
---|
2960 | |
---|
2961 | // Special case: Hidden-Valley massive photon. |
---|
2962 | double r3s = r3 * r3; |
---|
2963 | double prop3 = r3s - x3; |
---|
2964 | double prop3s = prop3 * prop3; |
---|
2965 | if (kind == 30) prop13 = prop1 * prop3; |
---|
2966 | |
---|
2967 | // Check input values. Return zero outside allowed phase space. |
---|
2968 | if (x1 - 2.*r1 < XMARGIN || prop1 < XMARGIN) return 0.; |
---|
2969 | if (x2 - 2.*r2 < XMARGIN || prop2 < XMARGIN) return 0.; |
---|
2970 | // Limits not worked out for r3 > 0. |
---|
2971 | if (kind != 30 && kind != 31) { |
---|
2972 | if (x1 + x2 - 1. - pow2(r1+r2) < XMARGIN) return 0.; |
---|
2973 | // Note: equivalent rewritten form 4. * ( (1. - x1) * (1. - x2) |
---|
2974 | // * (1. - r1s - r2s - x3) + r1s * (1. - x2s - x3) + r2s |
---|
2975 | // * (1. - x1s - x3) - pow2(r1s - r2s) ) gives about same result. |
---|
2976 | if ( (x1s - 4.*r1s) * (x2s - 4.*r2s) |
---|
2977 | - pow2( 2. * (1. - x1 - x2 + r1s + r2s) + x1*x2 ) |
---|
2978 | < XMARGIN * (XMARGINCOMB + r1 + r2) ) return 0.; |
---|
2979 | } |
---|
2980 | |
---|
2981 | // Initial values; phase space. |
---|
2982 | int combi = max(1, min(4, combiIn) ); |
---|
2983 | double mix = max(0., min(1., mixIn) ); |
---|
2984 | bool isSet1 = false; |
---|
2985 | bool isSet2 = false; |
---|
2986 | bool isSet4 = false; |
---|
2987 | double ps = sqrtpos( pow2(1. - r1*r1 - r2*r2) - pow2(2. * r1 * r2) ); |
---|
2988 | double rLO = 0., rFO = 0., rLO1 = 0., rFO1 = 0., rLO2 = 0., |
---|
2989 | rFO2 = 0., rLO4 = 0., rFO4 = 0.; |
---|
2990 | double offset = 0; |
---|
2991 | |
---|
2992 | // Select which kind of ME to use. |
---|
2993 | switch (kind) { |
---|
2994 | |
---|
2995 | // case 1 is equal to default, i.e. eikonal expression. |
---|
2996 | |
---|
2997 | // V -> q qbar (V = gamma*/Z0/W+-/...). |
---|
2998 | case 2: |
---|
2999 | if (combi == 1 || combi == 3) { |
---|
3000 | rLO1 = ps*(2.-r1s-r1q+6.*r1*r2-r2s+2.*r1s*r2s-r2q)/2.; |
---|
3001 | rFO1 = -(3.+6.*r1s+r1q-6.*r1*r2+6.*r1c*r2-2.*r2s-6.*r1s*r2s |
---|
3002 | +6.*r1*r2c+r2q-3.*x1+6.*r1*r2*x1+2.*r2s*x1+x1s-2.*r1s*x1s |
---|
3003 | +3.*r1s*x3+6.*r1*r2*x3-r2s*x3-2.*x1*x3-5.*r1s*x1*x3 |
---|
3004 | +r2s*x1*x3+x1s*x3-3.*x3s-3.*r1s*x3s+r2s*x3s |
---|
3005 | +2.*x1*x3s+x3c-x2) |
---|
3006 | /prop2s |
---|
3007 | -2.*(-3.+r1s-6.*r1*r2+6.*r1c*r2+3.*r2s-4.*r1s*r2s |
---|
3008 | +6.*r1*r2c+2.*x1+3.*r1s*x1+r2s*x1-x1s-r1s*x1s |
---|
3009 | -r2s*x1s+4.*x3+2.*r1s*x3+3.*r1*r2*x3-r2s*x3-3.*x1*x3 |
---|
3010 | -2.*r1s*x1*x3+x1s*x3-x3s-r1s*x3s+r1*r2*x3s+x1*x3s) |
---|
3011 | /prop12 |
---|
3012 | -(-1.+2.*r1s+r1q+6.*r1*r2+6.*r1c*r2-2.*r2s-6.*r1s*r2s |
---|
3013 | +6.*r1*r2c+r2q-x1-2.*r1s*x1-6.*r1*r2*x1+8.*r2s*x1+x1s |
---|
3014 | -2.*r2s*x1s-r1s*x3+r2s*x3-r1s*x1*x3+r2s*x1*x3+x1s*x3+x2) |
---|
3015 | /prop1s; |
---|
3016 | rFO1 = rFO1/2.; |
---|
3017 | isSet1 = true; |
---|
3018 | } |
---|
3019 | if (combi == 2 || combi == 3) { |
---|
3020 | rLO2 = ps*(2.-r1s-r1q-6.*r1*r2-r2s+2.*r1s*r2s-r2q)/2.; |
---|
3021 | rFO2 = -(3.+6.*r1s+r1q+6.*r1*r2-6.*r1c*r2-2.*r2s-6.*r1s*r2s |
---|
3022 | -6.*r1*r2c+r2q-3.*x1-6.*r1*r2*x1+2.*r2s*x1+x1s-2.*r1s*x1s |
---|
3023 | +3.*r1s*x3-6.*r1*r2*x3-r2s*x3-2.*x1*x3-5.*r1s*x1*x3 |
---|
3024 | +r2s*x1*x3+x1s*x3-3.*x3s-3.*r1s*x3s+r2s*x3s+2.*x1*x3s+x3c-x2) |
---|
3025 | /prop2s |
---|
3026 | -2.*(-3+r1s+6.*r1*r2-6.*r1c*r2+3.*r2s-4.*r1s*r2s-6.*r1*r2c |
---|
3027 | +2.*x1+3.*r1s*x1+r2s*x1-x1s-r1s*x1s-r2s*x1s+4.*x3+2.*r1s*x3 |
---|
3028 | -3.*r1*r2*x3-r2s*x3-3.*x1*x3-2.*r1s*x1*x3+x1s*x3-x3s-r1s*x3s |
---|
3029 | -r1*r2*x3s+x1*x3s) |
---|
3030 | /prop12 |
---|
3031 | -(-1.+2.*r1s+r1q-6.*r1*r2-6.*r1c*r2-2.*r2s-6.*r1s*r2s |
---|
3032 | -6.*r1*r2c+r2q-x1-2.*r1s*x1+6.*r1*r2*x1+8.*r2s*x1+x1s |
---|
3033 | -2.*r2s*x1s-r1s*x3+r2s*x3-r1s*x1*x3+r2s*x1*x3+x1s*x3+x2) |
---|
3034 | /prop1s; |
---|
3035 | rFO2 = rFO2/2.; |
---|
3036 | isSet2 = true; |
---|
3037 | } |
---|
3038 | if (combi == 4) { |
---|
3039 | rLO4 = ps*(2.-r1s-r1q-r2s+2.*r1s*r2s-r2q)/2.; |
---|
3040 | rFO4 = (1.-r1q+6.*r1s*r2s-r2q+x1+3.*r1s*x1-9.*r2s*x1-3.*x1s |
---|
3041 | -r1s*x1s+3.*r2s*x1s+x1c-x2-r1s*x2+r2s*x2-r1s*x1*x2+r2s*x1*x2 |
---|
3042 | +x1s*x2) |
---|
3043 | /prop1s |
---|
3044 | -2.*(1.+r1s+r2s-4.*r1s*r2s+r1s*x1+2.*r2s*x1-x1s-r2s*x1s |
---|
3045 | +2.*r1s*x2+r2s*x2-3.*x1*x2+x1s*x2-x2s-r1s*x2s+x1*x2s) |
---|
3046 | /prop12 |
---|
3047 | +(1.-r1q+6.*r1s*r2s-r2q-x1+r1s*x1-r2s*x1+x2-9.*r1s*x2 |
---|
3048 | +3.*r2s*x2+r1s*x1*x2-r2s*x1*x2-3.*x2s+3.*r1s*x2s-r2s*x2s |
---|
3049 | +x1*x2s+x2c) |
---|
3050 | /prop2s; |
---|
3051 | rFO4 = rFO4/2.; |
---|
3052 | isSet4 = true; |
---|
3053 | } |
---|
3054 | break; |
---|
3055 | |
---|
3056 | // q -> q V. |
---|
3057 | case 3: |
---|
3058 | if (combi == 1 || combi == 3) { |
---|
3059 | rLO1 = ps*(1.-2.*r1s+r1q+r2s-6.*r1*r2s+r1s*r2s-2.*r2q); |
---|
3060 | rFO1 = -2.*(-1.+r1-2.*r1s+2.*r1c-r1q+pow5(r1)-r2s+r1*r2s |
---|
3061 | -5.*r1s*r2s+r1c*r2s-2.*r1*r2q+2.*x1-2.*r1*x1+2.*r1s*x1 |
---|
3062 | -2.*r1c*x1+2.*r2s*x1+5.*r1*r2s*x1+r1s*r2s*x1+2.*r2q*x1 |
---|
3063 | -x1s+r1*x1s-r2s*x1s+3.*x2+4.*r1s*x2+r1q*x2+2.*r2s*x2 |
---|
3064 | +2.*r1s*r2s*x2-4.*x1*x2-2.*r1s*x1*x2-r2s*x1*x2+x1s*x2 |
---|
3065 | -2.*x2s-2.*r1s*x2s+x1*x2s) |
---|
3066 | /prop23 |
---|
3067 | +(2.*r2s+6.*r1*r2s-6.*r1s*r2s+6.*r1c*r2s+2.*r2q+6.*r1*r2q |
---|
3068 | -r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2-3.*r2s*x2-6.*r1*r2s*x2 |
---|
3069 | +9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s |
---|
3070 | +2.*r2s*x2s+x1*x2s) |
---|
3071 | /prop2s |
---|
3072 | +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1 |
---|
3073 | +r1q*x1-3.*r2s*x1+6.*r1*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s- |
---|
3074 | 2.*r1s*x1s+x1c+7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2+6.*r1*r2s*x2 |
---|
3075 | +r1s*r2s*x2-2.*r2q*x2-9.*x1*x2-3.*r1s*x1*x2+2.*r2s*x1*x2 |
---|
3076 | +2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s+x1*x2s) |
---|
3077 | /x3s; |
---|
3078 | isSet1 = true; |
---|
3079 | } |
---|
3080 | if (combi == 2 || combi == 3) { |
---|
3081 | rLO2 = ps*(1.-2.*r1s+r1q+r2s+6.*r1*r2s+r1s*r2s-2.*r2q); |
---|
3082 | rFO2 = 2*(1.+r1+2.*r1s+2.*r1c+r1q+pow5(r1)+r2s+r1*r2s |
---|
3083 | +5.*r1s*r2s+r1c*r2s-2.*r1*r2q-2.*x1-2.*r1*x1-2.*r1s*x1 |
---|
3084 | -2.*r1c*x1-2.*r2s*x1+5.*r1*r2s*x1-r1s*r2s*x1-2.*r2q*x1+x1s |
---|
3085 | +r1*x1s+r2s*x1s-3.*x2-4.*r1s*x2-r1q*x2-2.*r2s*x2 |
---|
3086 | -2.*r1s*r2s*x2+4.*x1*x2+2.*r1s*x1*x2+r2s*x1*x2-x1s*x2 |
---|
3087 | +2.*x2s+2.*r1s*x2s-x1*x2s) |
---|
3088 | /prop23 |
---|
3089 | +(2.*r2s-6.*r1*r2s-6.*r1s*r2s-6.*r1c*r2s+2.*r2q-6.*r1*r2q |
---|
3090 | -r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2-3.*r2s*x2+6.*r1*r2s*x2 |
---|
3091 | +9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s |
---|
3092 | +2.*r2s*x2s+x1*x2s) |
---|
3093 | /prop2s |
---|
3094 | +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1 |
---|
3095 | +r1q*x1-3.*r2s*x1-6.*r1*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s |
---|
3096 | -2.*r1s*x1s+x1c+7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2-6.*r1*r2s*x2 |
---|
3097 | +r1s*r2s*x2-2.*r2q*x2-9.*x1*x2-3.*r1s*x1*x2+2.*r2s*x1*x2 |
---|
3098 | +2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s+x1*x2s) |
---|
3099 | /x3s; |
---|
3100 | isSet2 = true; |
---|
3101 | } |
---|
3102 | if (combi == 4) { |
---|
3103 | rLO4 = ps*(1.-2.*r1s+r1q+r2s+r1s*r2s-2.*r2q); |
---|
3104 | rFO4 = 2*(1.+2.*r1s+r1q+r2s+5.*r1s*r2s-2.*x1-2.*r1s*x1 |
---|
3105 | -2.*r2s*x1-r1s*r2s*x1-2.*r2q*x1+x1s+r2s*x1s-3.*x2-4.*r1s*x2 |
---|
3106 | -r1q*x2-2.*r2s*x2-2.*r1s*r2s*x2+4.*x1*x2+2.*r1s*x1*x2+r2s*x1*x2 |
---|
3107 | -x1s*x2+2.*x2s+2.*r1s*x2s-x1*x2s) |
---|
3108 | /prop23 |
---|
3109 | +(2.*r2s-6.*r1s*r2s+2.*r2q-r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2 |
---|
3110 | -3.*r2s*x2+9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s |
---|
3111 | +2.*r2s*x2s+x1*x2s) |
---|
3112 | /prop2s |
---|
3113 | +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1 |
---|
3114 | +r1q*x1-3.*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s-2.*r1s*x1s+x1c |
---|
3115 | +7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2+r1s*r2s*x2-2.*r2q*x2-9.*x1*x2 |
---|
3116 | -3.*r1s*x1*x2+2.*r2s*x1*x2+2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s |
---|
3117 | +x1*x2s) |
---|
3118 | /x3s; |
---|
3119 | isSet4 = true; |
---|
3120 | } |
---|
3121 | break; |
---|
3122 | |
---|
3123 | // S -> q qbar (S = h0/H0/A0/H+-/...). |
---|
3124 | case 4: |
---|
3125 | if (combi == 1 || combi == 3) { |
---|
3126 | rLO1 = ps*(1.-r1s-r2s-2.*r1*r2); |
---|
3127 | rFO1 = -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1 |
---|
3128 | -r1s*x1+2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) |
---|
3129 | /prop1s |
---|
3130 | -2.*(r1s+r1q-2.*r1c*r2+r2s-6.*r1s*r2s-2.*r1*r2c+r2q-r1s*x1 |
---|
3131 | +r1*r2*x1+2.*r2s*x1+2.*r1s*x2+r1*r2*x2-r2s*x2-x1*x2) |
---|
3132 | /prop12 |
---|
3133 | -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1-r1s*x1 |
---|
3134 | +r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2) |
---|
3135 | /prop2s; |
---|
3136 | isSet1 = true; |
---|
3137 | } |
---|
3138 | if (combi == 2 || combi == 3) { |
---|
3139 | rLO2 = ps*(1.-r1s-r2s+2.*r1*r2); |
---|
3140 | rFO2 = -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 |
---|
3141 | -r1s*x1-2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) |
---|
3142 | /prop1s |
---|
3143 | -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 |
---|
3144 | -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2) |
---|
3145 | /prop2s |
---|
3146 | +2.*(-r1s-r1q-2.*r1c*r2-r2s+6.*r1s*r2s-2.*r1*r2c-r2q+r1s*x1 |
---|
3147 | +r1*r2*x1-2.*r2s*x1-2.*r1s*x2+r1*r2*x2+r2s*x2+x1*x2) |
---|
3148 | /prop12; |
---|
3149 | isSet2 = true; |
---|
3150 | } |
---|
3151 | if (combi == 4) { |
---|
3152 | rLO4 = ps*(1.-r1s-r2s); |
---|
3153 | rFO4 = -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+3.*r2s*x1+x2 |
---|
3154 | +r1s*x2-r2s*x2-x1*x2) |
---|
3155 | /prop1s |
---|
3156 | -2.*(r1s+r1q+r2s-6.*r1s*r2s+r2q-r1s*x1 |
---|
3157 | +2.*r2s*x1+2.*r1s*x2-r2s*x2-x1*x2) |
---|
3158 | /prop12 |
---|
3159 | -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1 |
---|
3160 | +x2+3.*r1s*x2-r2s*x2-x1*x2) |
---|
3161 | /prop2s; |
---|
3162 | isSet4 = true; |
---|
3163 | } |
---|
3164 | break; |
---|
3165 | |
---|
3166 | // q -> q S. |
---|
3167 | case 5: |
---|
3168 | if (combi == 1 || combi == 3) { |
---|
3169 | rLO1 = ps*(1.+r1s-r2s+2.*r1); |
---|
3170 | rFO1 = (4.-4.*r1s+4.*r2s-3.*x1-2.*r1*x1+r1s*x1-r2s*x1-5.*x2 |
---|
3171 | -2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3172 | /x3s |
---|
3173 | -2.*(3.-r1-5.*r1s-r1c+3.*r2s+r1*r2s-2.*x1-r1*x1 |
---|
3174 | +r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3175 | /prop23 |
---|
3176 | +(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1 |
---|
3177 | -r2s*x1-3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3178 | /prop2s; |
---|
3179 | isSet1 = true; |
---|
3180 | } |
---|
3181 | if (combi == 2 || combi == 3) { |
---|
3182 | rLO2 = ps*(1.+r1s-r2s-2.*r1); |
---|
3183 | rFO2 = (4.-4.*r1s+4.*r2s-3.*x1+2.*r1*x1+r1s*x1-r2s*x1-5.*x2 |
---|
3184 | +2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3185 | /x3s |
---|
3186 | -2.*(3.+r1-5.*r1s+r1c+3.*r2s-r1*r2s-2.*x1+r1*x1 |
---|
3187 | +r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3188 | /prop23 |
---|
3189 | +(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1 |
---|
3190 | -r2s*x1-3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3191 | /prop2s; |
---|
3192 | isSet2 = true; |
---|
3193 | } |
---|
3194 | if (combi == 4) { |
---|
3195 | rLO4 = ps*(1.+r1s-r2s); |
---|
3196 | rFO4 = (4.-4.*r1s+4.*r2s-3.*x1+r1s*x1-r2s*x1-5.*x2+r1s*x2 |
---|
3197 | -r2s*x2+x1*x2+x2s) |
---|
3198 | /x3s |
---|
3199 | -2.*(3.-5.*r1s+3.*r2s-2.*x1+r1s*x1-4.*x2+2.*r1s*x2 |
---|
3200 | -r2s*x2+x1*x2+x2s) |
---|
3201 | /prop23 |
---|
3202 | +(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2 |
---|
3203 | -r2s*x2+x1*x2+x2s) |
---|
3204 | /prop2s; |
---|
3205 | isSet4 = true; |
---|
3206 | } |
---|
3207 | break; |
---|
3208 | |
---|
3209 | // V -> ~q ~qbar (~q = squark). |
---|
3210 | case 6: |
---|
3211 | rLO1 = ps*(1.-2.*r1s+r1q-2.*r2s-2.*r1s*r2s+r2q); |
---|
3212 | rFO1 = 2.*3.+(1.+r1s+r2s-x1)*(4.*r1s-x1s) |
---|
3213 | /prop1s |
---|
3214 | +2.*(-1.-3.*r1s-r2s+x1+x1s*0.5+x2-x1*x2*0.5) |
---|
3215 | /prop1 |
---|
3216 | +(1.+r1s+r2s-x2)*(4.*r2s-x2s) |
---|
3217 | /prop2s |
---|
3218 | +2.*(-1.-r1s-3.*r2s+x1+x2-x1*x2*0.5+x2s*0.5) |
---|
3219 | /prop2 |
---|
3220 | -(-4.*r1s-4.*r1q-4.*r2s-8.*r1s*r2s-4.*r2q+2.*x1+6.*r1s*x1 |
---|
3221 | +6.*r2s*x1-2.*x1s+2.*x2+6.*r1s*x2+6.*r2s*x2-4.*x1*x2 |
---|
3222 | -2.*r1s*x1*x2-2.*r2s*x1*x2+x1s*x2-2.*x2s+x1*x2s) |
---|
3223 | /prop12; |
---|
3224 | isSet1 = true; |
---|
3225 | break; |
---|
3226 | |
---|
3227 | // ~q -> ~q V. |
---|
3228 | case 7: |
---|
3229 | rLO1 = ps*(1.-2.*r1s+r1q-2.*r2s-2.*r1s*r2s+r2q); |
---|
3230 | rFO1 = 16.*r2s-8.*(4.*r2s+2.*r2s*x1+x2+r1s*x2+r2s*x2-x1*x2 |
---|
3231 | -2.*x2s) |
---|
3232 | /(3.*prop2) |
---|
3233 | +8.*(1.+r1s+r2s-x2)*(4.*r2s-x2s) |
---|
3234 | /(3.*prop2s) |
---|
3235 | +8.*(x1+x2)*(-1.-2.*r1s-r1q-2.*r2s+2.*r1s*r2s-r2q+2.*x1 |
---|
3236 | +2.*r1s*x1+2.*r2s*x1-x1s+2.*x2+2.*r1s*x2+2.*r2s*x2-2.*x1*x2-x2s) |
---|
3237 | /(3.*x3s) |
---|
3238 | +8.*(-1.-r1s+r2s-x1)*(2.*r2s*x1+x2+r1s*x2+r2s*x2-x1*x2-x2s) |
---|
3239 | /(3.*prop2*x3) |
---|
3240 | -8.*(1.+2.*r1s+r1q+2.*r2s-2.*r1s*r2s+r2q-2.*x1-2.*r1s*x1 |
---|
3241 | -4.*r2s*x1+x1s-3.*x2-3.*r1s*x2-3.*r2s*x2+3.*x1*x2+2.*x2s) |
---|
3242 | /(3.*x3); |
---|
3243 | rFO1 = 3.*rFO1/8.; |
---|
3244 | isSet1 = true; |
---|
3245 | break; |
---|
3246 | |
---|
3247 | // S -> ~q ~qbar. |
---|
3248 | case 8: |
---|
3249 | rLO1 = ps; |
---|
3250 | rFO1 = (-1.-2.*r1s-r1q-2.*r2s+2.*r1s*r2s-r2q+2.*x1+2.*r1s*x1 |
---|
3251 | +2.*r2s*x1-x1s-r2s*x1s+2.*x2+2.*r1s*x2+2.*r2s*x2-3.*x1*x2 |
---|
3252 | -r1s*x1*x2-r2s*x1*x2+x1s*x2-x2s-r1s*x2s+x1*x2s) |
---|
3253 | /(prop1s*prop2s); |
---|
3254 | rFO1 = 2.*rFO1; |
---|
3255 | isSet1 = true; |
---|
3256 | break; |
---|
3257 | |
---|
3258 | // ~q -> ~q S. |
---|
3259 | case 9: |
---|
3260 | rLO1 = ps; |
---|
3261 | rFO1 = (-1.-r1s-r2s+x2) |
---|
3262 | /prop2s |
---|
3263 | +(1.+r1s-r2s+x1) |
---|
3264 | /prop23 |
---|
3265 | -(x1+x2) |
---|
3266 | /x3s; |
---|
3267 | isSet1 = true; |
---|
3268 | break; |
---|
3269 | |
---|
3270 | // chi -> q ~qbar (chi = neutralino/chargino). |
---|
3271 | case 10: |
---|
3272 | if (combi == 1 || combi == 3) { |
---|
3273 | rLO1 = ps*(1.+r1s-r2s+2.*r1); |
---|
3274 | rFO1 = (2.*r1+x1)*(-1.-r1s-r2s+x1) |
---|
3275 | /prop1s |
---|
3276 | +2.*(-1.-r1s-2.*r1c-r2s-2.*r1*r2s+3.*x1*0.5+r1*x1 |
---|
3277 | -r1s*x1*0.5-r2s*x1*0.5+x2+r1*x2+r1s*x2-x1*x2*0.5) |
---|
3278 | /prop12 |
---|
3279 | +(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1 |
---|
3280 | -r2s*x1-3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3281 | /prop2s; |
---|
3282 | isSet1 = true; |
---|
3283 | } |
---|
3284 | if (combi == 2 || combi == 3) { |
---|
3285 | rLO2 = ps*(1.-2.*r1+r1s-r2s); |
---|
3286 | rFO2 = (2.*r1-x1)*(1.+r1s+r2s-x1) |
---|
3287 | /prop1s |
---|
3288 | +2.*(-1.-r1s+2.*r1c-r2s+2.*r1*r2s+3.*x1*0.5-r1*x1 |
---|
3289 | -r1s*x1*0.5-r2s*x1*0.5+x2-r1*x2+r1s*x2-x1*x2*0.5) |
---|
3290 | /prop12 |
---|
3291 | +(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1 |
---|
3292 | -r2s*x1-3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)/ |
---|
3293 | prop2s; |
---|
3294 | isSet2 = true; |
---|
3295 | } |
---|
3296 | if (combi == 4) { |
---|
3297 | rLO4 = ps*(1.+r1s-r2s); |
---|
3298 | rFO4 = x1*(-1.-r1s-r2s+x1) |
---|
3299 | /prop1s |
---|
3300 | +2.*(-1.-r1s-r2s+3.*x1*0.5-r1s*x1*0.5-r2s*x1*0.5 |
---|
3301 | +x2+r1s*x2-x1*x2*0.5) |
---|
3302 | /prop12 |
---|
3303 | +(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2 |
---|
3304 | -r2s*x2+x1*x2+x2s) |
---|
3305 | /prop2s; |
---|
3306 | isSet4 = true; |
---|
3307 | } |
---|
3308 | break; |
---|
3309 | |
---|
3310 | // ~q -> q chi. |
---|
3311 | case 11: |
---|
3312 | if (combi == 1 || combi == 3) { |
---|
3313 | rLO1 = ps*(1.-pow2(r1+r2)); |
---|
3314 | rFO1 = (1.+r1s+2.*r1*r2+r2s-x1-x2)*(x1+x2) |
---|
3315 | /x3s |
---|
3316 | -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1 |
---|
3317 | -r1s*x1+r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2) |
---|
3318 | /prop2s |
---|
3319 | +(-1.-2.*r1s-r1q-2.*r1*r2-2.*r1c*r2+2.*r1*r2c+r2q+x1+r1s*x1 |
---|
3320 | -2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) |
---|
3321 | /prop23; |
---|
3322 | isSet1 = true; |
---|
3323 | } |
---|
3324 | if (combi == 2 || combi == 3) { |
---|
3325 | rLO2 = ps*(1.-pow2(r1-r2)); |
---|
3326 | rFO2 = (1.+r1s-2.*r1*r2+r2s-x1-x2)*(x1+x2) |
---|
3327 | /x3s |
---|
3328 | -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 |
---|
3329 | -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2) |
---|
3330 | /prop2s |
---|
3331 | +(-1.-2.*r1s-r1q+2.*r1*r2+2.*r1c*r2-2.*r1*r2c+r2q+x1+r1s*x1 |
---|
3332 | +2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) |
---|
3333 | /prop23; |
---|
3334 | isSet2 = true; |
---|
3335 | } |
---|
3336 | if (combi == 4) { |
---|
3337 | rLO4 = ps*(1.-r1s-r2s); |
---|
3338 | rFO4 = (1.+r1s+r2s-x1-x2)*(x1+x2) |
---|
3339 | /x3s |
---|
3340 | -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1+x2 |
---|
3341 | +3.*r1s*x2-r2s*x2-x1*x2) |
---|
3342 | /prop2s |
---|
3343 | +(-1.-2.*r1s-r1q+r2q+x1+r1s*x1-3.*r2s*x1 |
---|
3344 | +2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) |
---|
3345 | /prop23; |
---|
3346 | isSet4 = true; |
---|
3347 | } |
---|
3348 | break; |
---|
3349 | |
---|
3350 | // q -> ~q chi. |
---|
3351 | case 12: |
---|
3352 | if (combi == 1 || combi == 3) { |
---|
3353 | rLO1 = ps*(1.-r1s+r2s+2.*r2); |
---|
3354 | rFO1 = (2.*r2+x2)*(-1.-r1s-r2s+x2) |
---|
3355 | /prop2s |
---|
3356 | +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1-2.*r2*x1+r2s*x1+x1s |
---|
3357 | -3.*x2-r1s*x2-2.*r2*x2+r2s*x2+x1*x2) |
---|
3358 | /x3s |
---|
3359 | +2.*(-1.-r1s+r2+r1s*r2-r2s-r2c+x1+r2*x1+r2s*x1+2.*x2 |
---|
3360 | +r1s*x2-x1*x2*0.5-x2s*0.5) |
---|
3361 | /prop23; |
---|
3362 | isSet1 = true; |
---|
3363 | } |
---|
3364 | if (combi == 2 || combi == 3) { |
---|
3365 | rLO2 = ps*(1.-r1s+r2s-2.*r2); |
---|
3366 | rFO2 = (2.*r2-x2)*(1.+r1s+r2s-x2) |
---|
3367 | /prop2s |
---|
3368 | +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+2.*r2*x1+r2s*x1+x1s |
---|
3369 | -3.*x2-r1s*x2+2.*r2*x2+r2s*x2+x1*x2) |
---|
3370 | /x3s |
---|
3371 | +2.*(-1.-r1s-r2-r1s*r2-r2s+r2c+x1-r2*x1+r2s*x1+2.*x2 |
---|
3372 | +r1s*x2-x1*x2*0.5-x2s*0.5) |
---|
3373 | /prop23; |
---|
3374 | isSet2 = true; |
---|
3375 | } |
---|
3376 | if (combi == 4) { |
---|
3377 | rLO4 = ps*(1.-r1s+r2s); |
---|
3378 | rFO4 = x2*(-1.-r1s-r2s+x2) |
---|
3379 | /prop2s |
---|
3380 | +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+r2s*x1+x1s |
---|
3381 | -3.*x2-r1s*x2+r2s*x2+x1*x2) |
---|
3382 | /x3s |
---|
3383 | +2.*(-1.-r1s-r2s+x1+r2s*x1+2.*x2 |
---|
3384 | +r1s*x2-x1*x2*0.5-x2s*0.5) |
---|
3385 | /prop23; |
---|
3386 | isSet4 = true; |
---|
3387 | } |
---|
3388 | break; |
---|
3389 | |
---|
3390 | // ~g -> q ~qbar. |
---|
3391 | case 13: |
---|
3392 | if (combi == 1 || combi == 3) { |
---|
3393 | rLO1 = ps*(1.+r1s-r2s+2.*r1); |
---|
3394 | rFO1 = 4.*(2.*r1+x1)*(-1.-r1s-r2s+x1) |
---|
3395 | /(3.*prop1s) |
---|
3396 | -(-1.-r1s-2.*r1c-r2s-2.*r1*r2s+3.*x1*0.5+r1*x1-r1s*x1*0.5 |
---|
3397 | -r2s*x1*0.5+x2+r1*x2+r1s*x2-x1*x2*0.5) |
---|
3398 | /(3.*prop12) |
---|
3399 | +3.*(-1.+r1-r1s-r1c-r2s+r1*r2s+2.*x1+r2s*x1-x1s*0.5+x2+r1*x2 |
---|
3400 | +r1s*x2-x1*x2*0.5) |
---|
3401 | /prop13 |
---|
3402 | +3.*(4.-4.*r1s+4.*r2s-3.*x1-2.*r1*x1+r1s*x1-r2s*x1-5.*x2 |
---|
3403 | -2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3404 | /x3s |
---|
3405 | -3.*(3.-r1-5.*r1s-r1c+3.*r2s+r1*r2s-2.*x1-r1*x1+r1s*x1 |
---|
3406 | -4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3407 | /prop23 |
---|
3408 | +4.*(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1-r2s*x1 |
---|
3409 | -3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3410 | /(3.*prop2s); |
---|
3411 | rFO1 = 3.*rFO1/4.; |
---|
3412 | isSet1 = true; |
---|
3413 | } |
---|
3414 | if (combi == 2 || combi == 3) { |
---|
3415 | rLO2 = ps*(1.+r1s-r2s-2.*r1); |
---|
3416 | rFO2 = 4.*(2.*r1-x1)*(1.+r1s+r2s-x1) |
---|
3417 | /(3.*prop1s) |
---|
3418 | +3.*(-1.-r1-r1s+r1c-r2s-r1*r2s+2.*x1+r2s*x1-x1s*0.5 |
---|
3419 | +x2-r1*x2+r1s*x2-x1*x2*0.5) |
---|
3420 | /prop13 |
---|
3421 | +(2.+2.*r1s-4.*r1c+2.*r2s-4.*r1*r2s-3.*x1+2.*r1*x1 |
---|
3422 | +r1s*x1+r2s*x1-2.*x2+2.*r1*x2-2.*r1s*x2+x1*x2) |
---|
3423 | /(6.*prop12) |
---|
3424 | +3.*(4.-4.*r1s+4.*r2s-3.*x1+2.*r1*x1+r1s*x1-r2s*x1-5.*x2 |
---|
3425 | +2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3426 | /x3s |
---|
3427 | -3.*(3.+r1-5.*r1s+r1c+3.*r2s-r1*r2s-2.*x1+r1*x1+r1s*x1-4.*x2 |
---|
3428 | +2.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3429 | /prop23 |
---|
3430 | +4.*(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1-r2s*x1 |
---|
3431 | -3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3432 | /(3.*prop2s); |
---|
3433 | rFO2 = 3.*rFO2/4.; |
---|
3434 | isSet2 = true; |
---|
3435 | } |
---|
3436 | if (combi == 4) { |
---|
3437 | rLO4 = ps*(1.+r1s-r2s); |
---|
3438 | rFO4 = 8.*x1*(-1.-r1s-r2s+x1) |
---|
3439 | /(3.*prop1s) |
---|
3440 | +6.*(-1-r1s-r2s+2.*x1+r2s*x1-x1s*0.5+x2+r1s*x2-x1*x2*0.5) |
---|
3441 | /prop13 |
---|
3442 | +(2.+2.*r1s+2.*r2s-3.*x1+r1s*x1+r2s*x1-2.*x2-2.*r1s*x2+x1*x2) |
---|
3443 | /(3.*prop12) |
---|
3444 | +6.*(4.-4.*r1s+4.*r2s-3.*x1+r1s*x1-r2s*x1-5.*x2+r1s*x2-r2s*x2 |
---|
3445 | +x1*x2+x2s) |
---|
3446 | /x3s |
---|
3447 | -6.*(3.-5.*r1s+3.*r2s-2.*x1+r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) |
---|
3448 | /prop23 |
---|
3449 | +8.*(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2-r2s*x2 |
---|
3450 | +x1*x2+x2s) |
---|
3451 | /(3.*prop2s); |
---|
3452 | rFO4 = 3.*rFO4/8.; |
---|
3453 | isSet4 = true; |
---|
3454 | } |
---|
3455 | break; |
---|
3456 | |
---|
3457 | // ~q -> q ~g. |
---|
3458 | case 14: |
---|
3459 | if (combi == 1 || combi == 3) { |
---|
3460 | rLO1 = ps*(1.-r1s-r2s-2.*r1*r2); |
---|
3461 | rFO1 = 64.*(1.+r1s+2.*r1*r2+r2s-x1-x2)*(x1+x2) |
---|
3462 | /(9.*x3s) |
---|
3463 | -16.*(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q |
---|
3464 | +x1-r1s*x1+2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) |
---|
3465 | /prop1s |
---|
3466 | -16.*(r1s+r1q-2.*r1c*r2+r2s-6.*r1s*r2s-2.*r1*r2c+r2q-r1s*x1 |
---|
3467 | +r1*r2*x1+2.*r2s*x1+2.*r1s*x2+r1*r2*x2-r2s*x2-x1*x2) |
---|
3468 | /prop12 |
---|
3469 | -64.*(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1 |
---|
3470 | -r1s*x1+r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2) |
---|
3471 | /(9.*prop2s) |
---|
3472 | +8.*(-1.+r1q-2.*r1*r2+2.*r1c*r2-2.*r2s-2.*r1*r2c-r2q-2.*r1s*x1 |
---|
3473 | +2.*r2s*x1+x1s+x2-3.*r1s*x2-2.*r1*r2*x2+r2s*x2+x1*x2) |
---|
3474 | /prop13 |
---|
3475 | -8.*(-1.-2.*r1s-r1q-2.*r1*r2-2.*r1c*r2+2.*r1*r2c+r2q+x1+r1s*x1 |
---|
3476 | -2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) |
---|
3477 | /(9.*prop23); |
---|
3478 | rFO1 = 9.*rFO1/64.; |
---|
3479 | isSet1 = true; |
---|
3480 | } |
---|
3481 | if (combi == 2 || combi == 3) { |
---|
3482 | rLO2 = ps*(1.-r1s-r2s+2.*r1*r2); |
---|
3483 | rFO2 = 64.*(1.+r1s-2.*r1*r2+r2s-x1-x2)*(x1+x2) |
---|
3484 | /(9.*x3s) |
---|
3485 | -16.*(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 |
---|
3486 | -r1s*x1-2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) |
---|
3487 | /prop1s |
---|
3488 | -64.*(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 |
---|
3489 | -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2) |
---|
3490 | /(9.*prop2s) |
---|
3491 | +16.*(-r1s-r1q-2.*r1c*r2-r2s+6.*r1s*r2s-2.*r1*r2c-r2q+r1s*x1 |
---|
3492 | +r1*r2*x1-2.*r2s*x1-2.*r1s*x2+r1*r2*x2+r2s*x2+x1*x2) |
---|
3493 | /prop12 |
---|
3494 | +8.*(-1.+r1q+2.*r1*r2-2.*r1c*r2-2.*r2s+2.*r1*r2c-r2q-2.*r1s*x1 |
---|
3495 | +2.*r2s*x1+x1s+x2-3.*r1s*x2+2.*r1*r2*x2+r2s*x2+x1*x2) |
---|
3496 | /prop13 |
---|
3497 | -8.*(-1.-2.*r1s-r1q+2.*r1*r2+2.*r1c*r2-2.*r1*r2c+r2q+x1+r1s*x1+ |
---|
3498 | 2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) |
---|
3499 | /(9.*prop23); |
---|
3500 | rFO2 = 9.*rFO2/64.; |
---|
3501 | isSet2 = true; |
---|
3502 | } |
---|
3503 | if (combi == 4) { |
---|
3504 | rLO4 = ps*(1.-r1s-r2s); |
---|
3505 | rFO4 = 128.*(1.+r1s+r2s-x1-x2)*(x1+x2) |
---|
3506 | /(9.*x3s) |
---|
3507 | -32*(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+3.*r2s*x1+x2 |
---|
3508 | +r1s*x2-r2s*x2-x1*x2) |
---|
3509 | /prop1s |
---|
3510 | -32.*(r1s+r1q+r2s-6.*r1s*r2s+r2q-r1s*x1+2.*r2s*x1+2.*r1s*x2 |
---|
3511 | -r2s*x2-x1*x2) |
---|
3512 | /prop12 |
---|
3513 | -128.*(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1+x2+3.*r1s*x2 |
---|
3514 | -r2s*x2-x1*x2) |
---|
3515 | /(9.*prop2s) |
---|
3516 | +16.*(-1.+r1q-2.*r2s-r2q-2.*r1s*x1+2.*r2s*x1+x1s |
---|
3517 | +x2-3.*r1s*x2+r2s*x2+x1*x2) |
---|
3518 | /prop13 |
---|
3519 | -16.*(-1.-2.*r1s-r1q+r2q+x1+r1s*x1-3.*r2s*x1 |
---|
3520 | +2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) |
---|
3521 | /(9.*prop23); |
---|
3522 | rFO4 = 9.*rFO4/128.; |
---|
3523 | isSet4 = true; |
---|
3524 | } |
---|
3525 | break; |
---|
3526 | |
---|
3527 | // q -> ~q ~g. |
---|
3528 | case 15: |
---|
3529 | if (combi == 1 || combi == 3) { |
---|
3530 | rLO1 = ps*(1.-r1s+r2s+2.*r2); |
---|
3531 | rFO1 = 32*(2.*r2+x2)*(-1.-r1s-r2s+x2) |
---|
3532 | /(9.*prop2s) |
---|
3533 | +8.*(-1.-r1s-2.*r1s*r2-r2s-2.*r2c+x1+r2*x1+r2s*x1 |
---|
3534 | +3.*x2*0.5-r1s*x2*0.5+r2*x2-r2s*x2*0.5-x1*x2*0.5) |
---|
3535 | /prop12 |
---|
3536 | +8.*(2.+2.*r1s-2.*r2-2.*r1s*r2-6.*r2s-2.*r2c-3.*x1-r1s*x1 |
---|
3537 | +2.*r2*x1+3.*r2s*x1+x1s-x2-r1s*x2+r2s*x2+x1*x2) |
---|
3538 | /prop1s |
---|
3539 | +32.*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1-2.*r2*x1+r2s*x1+x1s |
---|
3540 | -3.*x2-r1s*x2-2.*r2*x2+r2s*x2+x1*x2) |
---|
3541 | /(9.*x3s) |
---|
3542 | -8.*(3.+3.*r1s-r2+r1s*r2-5.*r2s-r2c-4.*x1-r1s*x1 |
---|
3543 | +2.*r2s*x1+x1s-2.*x2-r2*x2+r2s*x2+x1*x2) |
---|
3544 | /prop13 |
---|
3545 | -8.*(-1.-r1s+r2+r1s*r2-r2s-r2c+x1+r2*x1+r2s*x1+2.*x2+r1s*x2 |
---|
3546 | -x1*x2*0.5-x2s*0.5) |
---|
3547 | /(9.*prop23); |
---|
3548 | rFO1 = 9.*rFO1/32.; |
---|
3549 | isSet1 = true; |
---|
3550 | } |
---|
3551 | if (combi == 2 || combi == 3) { |
---|
3552 | rLO2 = ps*(1.-r1s+r2s-2.*r2); |
---|
3553 | rFO2 = 32*(2.*r2-x2)*(1.+r1s+r2s-x2) |
---|
3554 | /(9.*prop2s) |
---|
3555 | +8.*(-1.-r1s+2.*r1s*r2-r2s+2.*r2c+x1-r2*x1+r2s*x1 |
---|
3556 | +3.*x2*0.5-r1s*x2*0.5-r2*x2-r2s*x2*0.5-x1*x2*0.5) |
---|
3557 | /prop12 |
---|
3558 | +8.*(2.+2.*r1s+2.*r2+2.*r1s*r2-6.*r2s+2.*r2c-3.*x1-r1s*x1 |
---|
3559 | -2.*r2*x1+3.*r2s*x1+x1s-x2-r1s*x2+r2s*x2+x1*x2) |
---|
3560 | /prop1s |
---|
3561 | -8.*(3.+3.*r1s+r2-r1s*r2-5.*r2s+r2c-4.*x1-r1s*x1+2.*r2s*x1+x1s |
---|
3562 | -2.*x2+r2*x2+r2s*x2+x1*x2) |
---|
3563 | /prop13 |
---|
3564 | +32*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+2.*r2*x1+r2s*x1 |
---|
3565 | +x1s-3.*x2-r1s*x2+2.*r2*x2+r2s*x2+x1*x2) |
---|
3566 | /(9.*x3s) |
---|
3567 | -8.*(-1.-r1s-r2-r1s*r2-r2s+r2c+x1-r2*x1+r2s*x1+2.*x2+r1s*x2 |
---|
3568 | -x1*x2*0.5-x2s*0.5) |
---|
3569 | /(9.*prop23); |
---|
3570 | rFO2 = 9.*rFO2/32.; |
---|
3571 | isSet2 = true; |
---|
3572 | } |
---|
3573 | if (combi == 4) { |
---|
3574 | rLO4 = ps*(1.-r1s+r2s); |
---|
3575 | rFO4 = 64.*x2*(-1.-r1s-r2s+x2) |
---|
3576 | /(9.*prop2s) |
---|
3577 | +16.*(-1.-r1s-r2s+x1+r2s*x1+3.*x2*0.5-r1s*x2*0.5 |
---|
3578 | -r2s*x2*0.5-x1*x2*0.5) |
---|
3579 | /prop12 |
---|
3580 | -16.*(3.+3.*r1s-5.*r2s-4.*x1-r1s*x1+2.*r2s*x1+x1s-2.*x2+r2s*x2 |
---|
3581 | +x1*x2) |
---|
3582 | /prop13 |
---|
3583 | +64.*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+r2s*x1+x1s-3.*x2 |
---|
3584 | -r1s*x2+r2s*x2+x1*x2) |
---|
3585 | /(9.*x3s) |
---|
3586 | +16.*(2.+2.*r1s-6.*r2s-3.*x1-r1s*x1+3.*r2s*x1+x1s |
---|
3587 | -x2-r1s*x2+r2s*x2+x1*x2) |
---|
3588 | /prop1s |
---|
3589 | -16.*(-1.-r1s-r2s+x1+r2s*x1+2.*x2+r1s*x2-x1*x2*0.5-x2s*0.5) |
---|
3590 | /(9.*prop23); |
---|
3591 | rFO4 = 9.*rFO4/64.; |
---|
3592 | isSet4 = true; |
---|
3593 | } |
---|
3594 | break; |
---|
3595 | |
---|
3596 | // g -> ~g ~g. Use (9/4)*eikonal. May be changed in the future. |
---|
3597 | case 16: |
---|
3598 | rLO = ps; |
---|
3599 | if (combi == 2) offset = x3s; |
---|
3600 | else if (combi == 3) offset = mix * x3s; |
---|
3601 | else if (combi == 4) offset = 0.5 * x3s; |
---|
3602 | rFO = ps * 4.5 * ( (x1+x2-1.+offset-r1s-r2s)/prop12 |
---|
3603 | - r1s/prop2s - r2s/prop1s ); |
---|
3604 | break; |
---|
3605 | |
---|
3606 | // Dv -> qv d. |
---|
3607 | case 30: |
---|
3608 | rLO = ps*(1.-r1s+r2s+2.*r2); |
---|
3609 | rFO = ( 0.5*r3s + 2.*r1q + 0.5*r2s*r3s + r2*r3s - 2.*r1s |
---|
3610 | - 0.5*r1s*r3s - 2.*r1s*r2s - 4.*r1s*r2 ) / prop2s |
---|
3611 | + ( -2. + 2.*r2q + 2.*r1q + 2.*r2s*r3s - 4.*r2 + 2.*r2*r3s |
---|
3612 | + 4.*r2*r2s - 4.*r1s*r2s - 4.*r1s*r2 ) /prop23 |
---|
3613 | + ( -2. - 0.5*r3s - 2.*r2s - 4.*r2 + 2.*r1s ) / prop2 |
---|
3614 | + ( -2. - r3s - 2.*r2s - r2s*r3s - 4.*r2 - 2.*r2*r3s |
---|
3615 | + 2.*r1s + r1s*r3s ) / prop3s |
---|
3616 | + ( -1. - r3s - r2s - 4.*r2 + r1s - x2 ) / prop3 |
---|
3617 | + 1.; |
---|
3618 | break; |
---|
3619 | |
---|
3620 | // S -> Dv Dvbar |
---|
3621 | case 31: |
---|
3622 | rLO = ps*(1.-4.*r1s); |
---|
3623 | rFO = (r3s + 2.*r1s) * (-1. + 4.*r1s) * (1./prop1s + 1./prop2s) |
---|
3624 | + (-1. + 8.*r1s - x2) / prop1 |
---|
3625 | + (-1. + 8.*r1s - x1) / prop2 |
---|
3626 | + 2. * (1. - 6.*r1s + 8.*r1q + 4.*r3s*r1s) / prop12 |
---|
3627 | + 2.; |
---|
3628 | break; |
---|
3629 | |
---|
3630 | // Eikonal expression for kind == 1; also acts as default. |
---|
3631 | default: |
---|
3632 | rLO = ps; |
---|
3633 | if (combi == 2) offset = x3s; |
---|
3634 | else if (combi == 3) offset = mix * x3s; |
---|
3635 | else if (combi == 4) offset = 0.5 * x3s; |
---|
3636 | rFO = ps * 2. * ( (x1+x2-1.+offset-r1s-r2s)/prop12 |
---|
3637 | - r1s/prop2s - r2s/prop1s ); |
---|
3638 | break; |
---|
3639 | |
---|
3640 | // End of ME cases. |
---|
3641 | } |
---|
3642 | |
---|
3643 | // Find relevant leading and first order expressions. |
---|
3644 | if (combi == 1 && isSet1) { |
---|
3645 | rLO = rLO1; |
---|
3646 | rFO = rFO1; } |
---|
3647 | else if (combi == 2 && isSet2) { |
---|
3648 | rLO = rLO2; |
---|
3649 | rFO = rFO2; } |
---|
3650 | else if (combi == 3 && isSet1 && isSet2) { |
---|
3651 | rLO = mix * rLO1 + (1.-mix) * rLO2; |
---|
3652 | rFO = mix * rFO1 + (1.-mix) * rFO2; } |
---|
3653 | else if (isSet4) { |
---|
3654 | rLO = rLO4; |
---|
3655 | rFO = rFO4; } |
---|
3656 | else if (combi == 4 && isSet1 && isSet2) { |
---|
3657 | rLO = 0.5 * (rLO1 + rLO2); |
---|
3658 | rFO = 0.5 * (rFO1 + rFO2); } |
---|
3659 | else if (isSet1) { |
---|
3660 | rLO = rLO1; |
---|
3661 | rFO = rFO1; } |
---|
3662 | |
---|
3663 | // Return ratio of first to leading order cross section. |
---|
3664 | return rFO / rLO; |
---|
3665 | } |
---|
3666 | |
---|
3667 | //-------------------------------------------------------------------------- |
---|
3668 | |
---|
3669 | // Find coefficient of azimuthal asymmetry from gluon polarization. |
---|
3670 | |
---|
3671 | void TimeShower::findAsymPol( Event& event, TimeDipoleEnd* dip) { |
---|
3672 | |
---|
3673 | // Default is no asymmetry. Only gluons are studied. |
---|
3674 | dip->asymPol = 0.; |
---|
3675 | dip->iAunt = 0; |
---|
3676 | int iRad = dip->iRadiator; |
---|
3677 | if (!doPhiPolAsym || event[iRad].id() != 21) return; |
---|
3678 | |
---|
3679 | // Trace grandmother via possibly intermediate recoil copies. |
---|
3680 | int iMother = event.iTopCopy(iRad); |
---|
3681 | int iGrandM = event[iMother].mother1(); |
---|
3682 | |
---|
3683 | // If grandmother in initial state of hard scattering, |
---|
3684 | // then only keep gg and qq initial states. |
---|
3685 | int statusGrandM = event[iGrandM].status(); |
---|
3686 | bool isHardProc = (statusGrandM == -21 || statusGrandM == -31); |
---|
3687 | if (isHardProc) { |
---|
3688 | if (event[iGrandM + 1].status() != statusGrandM) return; |
---|
3689 | if (event[iGrandM].isGluon() && event[iGrandM + 1].isGluon()); |
---|
3690 | else if (event[iGrandM].isQuark() && event[iGrandM + 1].isQuark()); |
---|
3691 | else return; |
---|
3692 | } |
---|
3693 | |
---|
3694 | // Set aunt by history or, for hard scattering, by colour flow. |
---|
3695 | if (isHardProc) dip->iAunt = dip->iRecoiler; |
---|
3696 | else dip->iAunt = (event[iGrandM].daughter1() == iMother) |
---|
3697 | ? event[iGrandM].daughter2() : event[iGrandM].daughter1(); |
---|
3698 | |
---|
3699 | // Coefficient from gluon production (approximate z by energy). |
---|
3700 | // For hard process arbitrarily put z = 1/2. |
---|
3701 | double zProd = (isHardProc) ? 0.5 : event[iRad].e() |
---|
3702 | / (event[iRad].e() + event[dip->iAunt].e()); |
---|
3703 | if (event[iGrandM].isGluon()) dip->asymPol = pow2( (1. - zProd) |
---|
3704 | / (1. - zProd * (1. - zProd) ) ); |
---|
3705 | else dip->asymPol = 2. * (1. - zProd) / (1. + pow2(1. - zProd) ); |
---|
3706 | |
---|
3707 | // Coefficients from gluon decay. |
---|
3708 | if (dip->flavour == 21) dip->asymPol *= pow2( (1. - dip->z) |
---|
3709 | / (1. - dip->z * (1. - dip->z) ) ); |
---|
3710 | else dip->asymPol *= -2. * dip->z *( 1. - dip->z ) |
---|
3711 | / (1. - 2. * dip->z * (1. - dip->z) ); |
---|
3712 | |
---|
3713 | } |
---|
3714 | |
---|
3715 | //-------------------------------------------------------------------------- |
---|
3716 | |
---|
3717 | // Print the list of dipoles. |
---|
3718 | |
---|
3719 | void TimeShower::list(ostream& os) const { |
---|
3720 | |
---|
3721 | // Header. |
---|
3722 | os << "\n -------- PYTHIA TimeShower Dipole Listing ----------------" |
---|
3723 | << "--------------------------------------------- \n \n i rad" |
---|
3724 | << " rec pTmax col chg gam oni hv isr sys sysR typ" |
---|
3725 | << "e MErec mix ord spl ~gR \n" << fixed << setprecision(3); |
---|
3726 | |
---|
3727 | // Loop over dipole list and print it. |
---|
3728 | for (int i = 0; i < int(dipEnd.size()); ++i) |
---|
3729 | os << setw(5) << i << setw(7) << dipEnd[i].iRadiator |
---|
3730 | << setw(7) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax |
---|
3731 | << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType |
---|
3732 | << setw(5) << dipEnd[i].gamType << setw(5) << dipEnd[i].isOctetOnium |
---|
3733 | << setw(5) << dipEnd[i].isHiddenValley << setw(5) << dipEnd[i].isrType |
---|
3734 | << setw(5) << dipEnd[i].system << setw(5) << dipEnd[i].systemRec |
---|
3735 | << setw(5) << dipEnd[i].MEtype << setw(7) << dipEnd[i].iMEpartner |
---|
3736 | << setw(8) << dipEnd[i].MEmix << setw(5) << dipEnd[i].MEorder |
---|
3737 | << setw(5) << dipEnd[i].MEsplit << setw(5) << dipEnd[i].MEgluinoRec |
---|
3738 | << "\n"; |
---|
3739 | |
---|
3740 | // Done. |
---|
3741 | os << "\n -------- End PYTHIA TimeShower Dipole Listing ------------" |
---|
3742 | << "---------------------------------------------" << endl; |
---|
3743 | |
---|
3744 | } |
---|
3745 | |
---|
3746 | //========================================================================== |
---|
3747 | |
---|
3748 | } // end namespace Pythia8 |
---|
3749 | |
---|