StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireSplittingsQED.cc
1 // DireSplittingsQED.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Stefan Prestel, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) for the
7 // DireSplittingQED and derived classes.
8 
9 #include "Pythia8/DireSplittingsQED.h"
10 #include "Pythia8/DireSpace.h"
11 #include "Pythia8/DireTimes.h"
12 
13 namespace Pythia8 {
14 
15 //==========================================================================
16 
17 double chgprefac = 0.0;
18 
19 // The SplittingQED class.
20 
21 //--------------------------------------------------------------------------
22 
23 void DireSplittingQED::init() {
24 
25  int nGammaToQuark = settingsPtr->mode("TimeShower:nGammaToQuark");
26  int nGammaToLepton = settingsPtr->mode("TimeShower:nGammaToLepton");
27 
28  sumCharge2L = max(0, min(3, nGammaToLepton));
29  sumCharge2Q = 0.;
30  if (nGammaToQuark > 4) sumCharge2Q = 11. / 9.;
31  else if (nGammaToQuark > 3) sumCharge2Q = 10. / 9.;
32  else if (nGammaToQuark > 2) sumCharge2Q = 6. / 9.;
33  else if (nGammaToQuark > 1) sumCharge2Q = 5. / 9.;
34  else if (nGammaToQuark > 0) sumCharge2Q = 1. / 9.;
35  sumCharge2Tot = sumCharge2L + 3. * sumCharge2Q;
36 
37  // Parameters of alphaEM.
38  int alphaEMorder = settingsPtr->mode("SpaceShower:alphaEMorder");
39  // Initialize alphaEM.
40  alphaEM.init( alphaEMorder, settingsPtr);
41 
42  aem0 = settingsPtr->parm("StandardModel:alphaEM0");
43 
44  enhance = settingsPtr->parm("Enhance:"+ id);
45 
46  doQEDshowerByQ = (is_fsr) ? settingsPtr->flag("TimeShower:QEDshowerByQ")
47  : settingsPtr->flag("SpaceShower:QEDshowerByQ");
48  doQEDshowerByL = (is_fsr) ? settingsPtr->flag("TimeShower:QEDshowerByL")
49  : settingsPtr->flag("SpaceShower:QEDshowerByL");
50  doForcePos = settingsPtr->flag("Dire:QED:doForcePosChgCorrelators");
51  pT2minForcePos = pow2(settingsPtr->parm("Dire:QED:pTminForcePos"));
52 
53  pT2min = pow2(settingsPtr->parm("TimeShower:pTmin"));
54  pT2minL = pow2(settingsPtr->parm("TimeShower:pTminChgL"));
55  pT2minQ = pow2(settingsPtr->parm("TimeShower:pTminChgQ"));
56 
57 }
58 
59 //--------------------------------------------------------------------------
60 
61 // Function to calculate the correct alphaem/2*Pi value, including
62 // renormalisation scale variations + threshold matching.
63 
64 double DireSplittingQED::aem2Pi( double pT2, int ) {
65 
66  double scale = pT2*renormMultFac;
67 
68  // Get alphaEM(k*pT^2) and subtractions.
69  double aemPT2pi = alphaEM.alphaEM(scale) / (2.*M_PI);
70 
71  // Done.
72  return aemPT2pi;
73 
74 }
75 
76 //--------------------------------------------------------------------------
77 
78 bool DireSplittingQED::aboveCutoff( double t, const Particle& radBef,
79  const Particle&, int iSys, PartonSystems* partonSystemsPtr) {
80 
81  if (particleDataPtr->isLepton(radBef.id()) && t < pT2minL )
82  return false;
83  if (particleDataPtr->isQuark(radBef.id()) && t < pT2minQ )
84  return false;
85  if ((iSys == 0 || partonSystemsPtr->hasInAB(iSys)) && t < pT2min)
86  return false;
87 
88  return true;
89 }
90 
91 //==========================================================================
92 
93 // Class inheriting from SplittingQED class.
94 
95 // SplittingQED function Q->QG (FSR)
96 
97 // Return true if this kernel should partake in the evolution.
98 bool Dire_fsr_qed_Q2QA::canRadiate ( const Event& state, pair<int,int> ints,
99  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
100  return ( state[ints.first].isFinal()
101  && state[ints.first].isQuark() && state[ints.second].isCharged()
102  && bools["doQEDshowerByQ"]);
103 }
104 
105 bool Dire_fsr_qed_Q2QA::canRadiate ( const Event& state, int iRadBef,
106  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
107  return ( state[iRadBef].isFinal()
108  && state[iRadBef].isQuark() && state[iRecBef].isCharged()
109  && doQEDshowerByQ);
110 }
111 
112 int Dire_fsr_qed_Q2QA::kinMap() { return 1;}
113 int Dire_fsr_qed_Q2QA::motherID(int idDaughter) { return idDaughter;}
114 int Dire_fsr_qed_Q2QA::sisterID(int) { return 22;}
115 
116 double Dire_fsr_qed_Q2QA::gaugeFactor ( int idRadBef, int idRecBef) {
117  double chgRad = particleDataPtr->charge(idRadBef);
118  double chgRec = particleDataPtr->charge(idRecBef);
119  double charge = -1.*chgRad*chgRec;
120  if (!splitInfo.radBef()->isFinal) charge *= -1.;
121  if (!splitInfo.recBef()->isFinal) charge *= -1.;
122  if (idRadBef != 0 && idRecBef != 0) return charge;
123  // Set probability to zero.
124  return 0.;
125 }
126 
127 double Dire_fsr_qed_Q2QA::symmetryFactor ( int, int ) { return 1.;}
128 
129 int Dire_fsr_qed_Q2QA::radBefID(int idRad, int idEmt) {
130  if (particleDataPtr->isQuark(idRad) && idEmt == 22 ) return idRad;
131  return 0;
132 }
133 
134 pair<int,int> Dire_fsr_qed_Q2QA::radBefCols(
135  int colRadAfter, int acolRadAfter,
136  int, int) { return make_pair(colRadAfter,acolRadAfter); }
137 
138 vector<int>Dire_fsr_qed_Q2QA::recPositions(const Event& state, int iRad,
139  int iEmt) {
140 
141  vector<int> recs;
142  if ( !state[iRad].isFinal()
143  || !state[iRad].isQuark()
144  || state[iEmt].id() != 22) return recs;
145 
146  // Particles to exclude as recoilers.
147  vector<int> iExc(createvector<int>(iRad)(iEmt));
148  // Find charged particles.
149  for (int i=0; i < state.size(); ++i) {
150  if ( find(iExc.begin(), iExc.end(), i) != iExc.end() ) continue;
151  if ( state[i].isCharged() ) {
152  if (state[i].isFinal())
153  recs.push_back(i);
154  if (state[i].mother1() == 1 && state[i].mother2() == 0)
155  recs.push_back(i);
156  if (state[i].mother1() == 2 && state[i].mother2() == 0)
157  recs.push_back(i);
158  }
159  }
160  // Done.
161  return recs;
162 
163 }
164 
165 // Pick z for new splitting.
166 double Dire_fsr_qed_Q2QA::zSplit(double zMinAbs, double, double m2dip) {
167  double Rz = rndmPtr->flat();
168  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/m2dip;
169  double p = pow( 1. + pow2(1-zMinAbs)/kappa2, Rz );
170  double res = 1. - sqrt( p - 1. )*sqrt(kappa2);
171  return res;
172 }
173 
174 // New overestimates, z-integrated versions.
175 double Dire_fsr_qed_Q2QA::overestimateInt(double zMinAbs, double,
176  double, double m2dip, int) {
177  double wt = 0.;
178  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
179  double preFac = symmetryFactor() * abs(charge);
180  // Q -> QG, soft part (currently also used for collinear part).
181  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/m2dip;
182  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa2);
183  return wt;
184 }
185 
186 // Return overestimate for new splitting.
187 double Dire_fsr_qed_Q2QA::overestimateDiff(double z, double m2dip, int) {
188  double wt = 0.;
189  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
190  double preFac = symmetryFactor() * abs(charge);
191  double kappaOld2 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/m2dip;
192  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappaOld2);
193  return wt;
194 }
195 
196 // Return kernel for new splitting.
197 bool Dire_fsr_qed_Q2QA::calc(const Event& state, int orderNow) {
198 
199  // Dummy statement to avoid compiler warnings.
200  if (false) cout << state[0].e() << orderNow << endl;
201 
202  // Read all splitting variables.
203  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
204  m2dip(splitInfo.kinematics()->m2Dip),
205  m2RadBef(splitInfo.kinematics()->m2RadBef),
206  m2Rad(splitInfo.kinematics()->m2RadAft),
207  m2Rec(splitInfo.kinematics()->m2Rec),
208  m2Emt(splitInfo.kinematics()->m2EmtAft);
209  int splitType(splitInfo.type);
210 
211  // Calculate kernel.
212  // Note: We are calculating the z <--> 1-z symmetrised kernel here,
213  // and later multiply with z to project out Q->QQ,
214  // i.e. the gluon is soft and the quark is identified.
215  double wt = 0.;
216 
217  double chargeFac = gaugeFactor(splitInfo.radBef()->id,
218  splitInfo.recBef()->id);
219  vector <int> in, out;
220  for (int i=0; i < state.size(); ++i) {
221  if (state[i].isFinal()) out.push_back(state[i].id());
222  if (state[i].mother1() == 1 && state[i].mother2() == 0)
223  in.push_back(state[i].id());
224  if (state[i].mother1() == 2 && state[i].mother2() == 0)
225  in.push_back(state[i].id());
226  }
227  out.push_back(22);
228  bool hasME = pT2 > pow2(settingsPtr->parm("Dire:pTminMECs"))
229  && doMECs && fsr->weights->hasME(in,out);
230  if (hasME && chargeFac < 0.0) chargeFac = abs(chargeFac);
231 
232  if ( doForcePos
233  && (chargeFac < 0. || splitInfo.radBef()->id != splitInfo.recBef()->id)
234  && (hasME || pT2 > pT2minForcePos)) chargeFac = chgprefac*abs(chargeFac);
235 
236  double preFac = symmetryFactor() * chargeFac;
237  double kappa2 = pT2/m2dip;
238  wt = preFac * ( 2. * z * (1.-z) / ( pow2(1.-z) + kappa2) );
239 
240  // Correction for massive splittings.
241  bool doMassive = (abs(splitType) == 2);
242 
243  // Add collinear term for massless splittings.
244  if (!doMassive && orderNow >= 0) wt += preFac * ( 1.-z );
245 
246  // Add collinear term for massive splittings.
247  if (doMassive && orderNow >= 0) {
248 
249  double pipj = 0., vijkt = 1., vijk = 1.;
250 
251  // splitType == 2 -> Massive FF
252  if (splitType == 2) {
253 
254  // Calculate CS variables.
255  double yCS = kappa2 / (1.-z);
256  double nu2RadBef = m2RadBef/m2dip;
257  double nu2Rad = m2Rad/m2dip;
258  double nu2Emt = m2Emt/m2dip;
259  double nu2Rec = m2Rec/m2dip;
260  vijk = pow2(1.-yCS) - 4.*(yCS+nu2Rad+nu2Emt)*nu2Rec;
261  double Q2mass = m2dip + m2Rad + m2Rec + m2Emt;
262  vijkt = pow2(Q2mass/m2dip - nu2RadBef - nu2Rec)
263  - 4.*nu2RadBef*nu2Rec;
264  vijk = sqrt(vijk) / (1-yCS);
265  vijkt = sqrt(vijkt)/ (Q2mass/m2dip - nu2RadBef - nu2Rec);
266  pipj = m2dip * yCS/2.;
267  // splitType ==-2 -> Massive FI
268  } else if (splitType ==-2) {
269 
270  // Calculate CS variables.
271  double xCS = 1 - kappa2/(1.-z);
272  vijk = 1.;
273  vijkt = 1.;
274  pipj = m2dip/2. * (1-xCS)/xCS;
275  }
276 
277  // Add B1 for massive splittings.
278  double massCorr = vijkt/vijk*( 1. - z - m2RadBef/pipj);
279  wt += preFac*massCorr;
280 
281  }
282 
283  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
284 
285  // Now multiply with z to project out Q->QG,
286  // i.e. the gluon is soft and the quark is identified.
287  wt *= z;
288 
289  // Trivial map of values, since kernel does not depend on coupling.
290  unordered_map<string,double> wts;
291  wts.insert( make_pair("base", wt ));
292  if (doVariations) {
293  // Create muR-variations.
294  if (settingsPtr->parm("Variations:muRfsrDown") != 1.)
295  wts.insert( make_pair("Variations:muRfsrDown", wt ));
296  if (settingsPtr->parm("Variations:muRfsrUp") != 1.)
297  wts.insert( make_pair("Variations:muRfsrUp", wt ));
298  }
299 
300  // Store kernel values.
301  clearKernels();
302  for ( unordered_map<string,double>::iterator it = wts.begin();
303  it != wts.end(); ++it )
304  kernelVals.insert(make_pair( it->first, it->second ));
305 
306  return true;
307 
308 }
309 
310 //==========================================================================
311 
312 // Class inheriting from SplittingQED class.
313 
314 // SplittingQED function Q->GQ (FSR)
315 // At leading order, this can be combined with Q->QG because of symmetry. Since
316 // this is no longer possible at NLO, we keep the kernels separately.
317 
318 // Return true if this kernel should partake in the evolution.
319 bool Dire_fsr_qed_Q2AQ::canRadiate ( const Event& state, pair<int,int> ints,
320  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
321  return ( state[ints.first].isFinal()
322  && state[ints.first].isQuark() && state[ints.second].isCharged()
323  && bools["doQEDshowerByQ"]);
324 }
325 
326 bool Dire_fsr_qed_Q2AQ::canRadiate ( const Event& state, int iRadBef,
327  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
328  return ( state[iRadBef].isFinal()
329  && state[iRadBef].isQuark() && state[iRecBef].isCharged()
330  && doQEDshowerByQ);
331 }
332 
333 int Dire_fsr_qed_Q2AQ::kinMap() { return 1;}
334 int Dire_fsr_qed_Q2AQ::motherID(int idDaughter) { return idDaughter;}
335 int Dire_fsr_qed_Q2AQ::sisterID(int) { return 22;}
336 
337 double Dire_fsr_qed_Q2AQ::gaugeFactor ( int idRadBef, int idRecBef) {
338  double chgRad = particleDataPtr->charge(idRadBef);
339  double chgRec = particleDataPtr->charge(idRecBef);
340  double charge = -1.*chgRad*chgRec;
341  if (!splitInfo.radBef()->isFinal) charge *= -1.;
342  if (!splitInfo.recBef()->isFinal) charge *= -1.;
343  if (idRadBef != 0 && idRecBef != 0) return charge;
344  // Set probability to zero.
345  return 0.;
346 }
347 
348 double Dire_fsr_qed_Q2AQ::symmetryFactor ( int, int ) { return 1.;}
349 
350 int Dire_fsr_qed_Q2AQ::radBefID(int idRad, int idEmt) {
351  if (idRad == 22 && particleDataPtr->isQuark(idEmt)) return idEmt;
352  if (idEmt == 22 && particleDataPtr->isQuark(idRad)) return idRad;
353  return 0;
354 }
355 
356 pair<int,int> Dire_fsr_qed_Q2AQ::radBefCols(
357  int colRadAfter, int acolRadAfter,
358  int, int) { return make_pair(colRadAfter,acolRadAfter); }
359 
360 vector<int>Dire_fsr_qed_Q2AQ::recPositions(const Event& state, int iRad,
361  int iEmt) {
362 
363  vector<int> recs;
364  if ( !state[iRad].isFinal()
365  || !state[iRad].isQuark()
366  || state[iEmt].id() != 22) return recs;
367 
368  // Particles to exclude as recoilers.
369  vector<int> iExc(createvector<int>(iRad)(iEmt));
370  // Find charged particles.
371  for (int i=0; i < state.size(); ++i) {
372  if ( find(iExc.begin(), iExc.end(), i) != iExc.end() ) continue;
373  if ( state[i].isCharged() ) {
374  if (state[i].isFinal())
375  recs.push_back(i);
376  if (state[i].mother1() == 1 && state[i].mother2() == 0)
377  recs.push_back(i);
378  if (state[i].mother1() == 2 && state[i].mother2() == 0)
379  recs.push_back(i);
380  }
381  }
382  // Done.
383  return recs;
384 
385 }
386 
387 // Pick z for new splitting.
388 double Dire_fsr_qed_Q2AQ::zSplit(double zMinAbs, double, double m2dip) {
389  double Rz = rndmPtr->flat();
390  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/m2dip;
391  double p = pow( 1. + pow2(1-zMinAbs)/kappa2, Rz );
392  double res = 1. - sqrt( p - 1. )*sqrt(kappa2);
393  return res;
394 }
395 
396 // New overestimates, z-integrated versions.
397 double Dire_fsr_qed_Q2AQ::overestimateInt(double zMinAbs, double,
398  double, double m2dip, int) {
399  double wt = 0.;
400  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
401  double preFac = symmetryFactor() * abs(charge);
402  // Q -> QG, soft part (currently also used for collinear part).
403  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/m2dip;
404  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa2);
405  return wt;
406 }
407 
408 // Return overestimate for new splitting.
409 double Dire_fsr_qed_Q2AQ::overestimateDiff(double z, double m2dip, int) {
410  double wt = 0.;
411  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
412  double preFac = symmetryFactor() * abs(charge);
413  double kappaOld2 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/m2dip;
414  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappaOld2);
415  return wt;
416 }
417 
418 // Return kernel for new splitting.
419 bool Dire_fsr_qed_Q2AQ::calc(const Event& state, int orderNow) {
420 
421  // Dummy statement to avoid compiler warnings.
422  if (false) cout << state[0].e() << orderNow << endl;
423 
424  // Read all splitting variables.
425  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
426  m2dip(splitInfo.kinematics()->m2Dip),
427  m2RadBef(splitInfo.kinematics()->m2RadBef),
428  m2Rad(splitInfo.kinematics()->m2RadAft),
429  m2Rec(splitInfo.kinematics()->m2Rec),
430  m2Emt(splitInfo.kinematics()->m2EmtAft);
431  int splitType(splitInfo.type);
432 
433  // Calculate kernel.
434  // Note: We are calculating the z <--> 1-z symmetrised kernel here,
435  // and later multiply with 1-z to project out Q->GQ,
436  // i.e. the quark is soft and the gluon is identified.
437  double wt = 0.;
438 
439  double chargeFac = gaugeFactor(splitInfo.radBef()->id,
440  splitInfo.recBef()->id);
441  vector <int> in, out;
442  for (int i=0; i < state.size(); ++i) {
443  if (state[i].isFinal()) out.push_back(state[i].id());
444  if (state[i].mother1() == 1 && state[i].mother2() == 0)
445  in.push_back(state[i].id());
446  if (state[i].mother1() == 2 && state[i].mother2() == 0)
447  in.push_back(state[i].id());
448  }
449  out.push_back(22);
450  bool hasME = pT2 > pow2(settingsPtr->parm("Dire:pTminMECs"))
451  && doMECs && fsr->weights->hasME(in,out);
452  if (hasME && chargeFac < 0.0) chargeFac = abs(chargeFac);
453 
454  if ( doForcePos
455  && (chargeFac < 0. || splitInfo.radBef()->id != splitInfo.recBef()->id)
456  && (hasME || pT2 > pT2minForcePos)) chargeFac = chgprefac*abs(chargeFac);
457 
458  double preFac = symmetryFactor() * chargeFac;
459  double kappa2 = pT2/m2dip;
460  wt = preFac * ( 2. * z * (1.-z) / ( pow2(1.-z) + kappa2) );
461 
462  // Correction for massive splittings.
463  bool doMassive = (abs(splitType) == 2);
464 
465  // Add collinear term for massless splittings.
466  if (!doMassive && orderNow >= 0) wt += preFac * ( 1.-z );
467 
468  // Add collinear term for massive splittings.
469  if (doMassive && orderNow >= 0) {
470 
471  double pipj = 0., vijkt = 1., vijk = 1.;
472 
473  // splitType == 2 -> Massive FF
474  if (splitType == 2) {
475 
476  // Calculate CS variables.
477  double yCS = kappa2 / (1.-z);
478  double nu2RadBef = m2RadBef/m2dip;
479  double nu2Rad = m2Rad/m2dip;
480  double nu2Emt = m2Emt/m2dip;
481  double nu2Rec = m2Rec/m2dip;
482  vijk = pow2(1.-yCS) - 4.*(yCS+nu2Rad+nu2Emt)*nu2Rec;
483  double Q2mass = m2dip + m2Rad + m2Rec + m2Emt;
484  vijkt = pow2(Q2mass/m2dip - nu2RadBef - nu2Rec)
485  - 4.*nu2RadBef*nu2Rec;
486  vijk = sqrt(vijk) / (1-yCS);
487  vijkt = sqrt(vijkt)/ (Q2mass/m2dip - nu2RadBef - nu2Rec);
488  pipj = m2dip * yCS/2.;
489  // splitType ==-2 -> Massive FI
490  } else if (splitType ==-2) {
491 
492  // Calculate CS variables.
493  double xCS = 1 - kappa2/(1.-z);
494  vijk = 1.;
495  vijkt = 1.;
496  pipj = m2dip/2. * (1-xCS)/xCS;
497  }
498 
499  // Add B1 for massive splittings.
500  double massCorr = vijkt/vijk*( 1. - z - m2RadBef/pipj);
501  wt += preFac*massCorr;
502 
503  }
504 
505  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
506 
507  // Now multiply with (1-z) to project out Q->GQ,
508  // i.e. the quark is soft and the gluon is identified.
509  wt *= ( 1. - z );
510 
511  // Trivial map of values, since kernel does not depend on coupling.
512  unordered_map<string,double> wts;
513  wts.insert( make_pair("base", wt ));
514  if (doVariations) {
515  // Create muR-variations.
516  if (settingsPtr->parm("Variations:muRfsrDown") != 1.)
517  wts.insert( make_pair("Variations:muRfsrDown", wt ));
518  if (settingsPtr->parm("Variations:muRfsrUp") != 1.)
519  wts.insert( make_pair("Variations:muRfsrUp", wt ));
520  }
521 
522  // Store kernel values.
523  clearKernels();
524  for ( unordered_map<string,double>::iterator it = wts.begin();
525  it != wts.end(); ++it )
526  kernelVals.insert(make_pair( it->first, it->second ));
527 
528  return true;
529 
530 }
531 
532 //==========================================================================
533 
534 // Class inheriting from SplittingQED class.
535 
536 // SplittingQED function Q->QG (FSR)
537 
538 // Return true if this kernel should partake in the evolution.
539 bool Dire_fsr_qed_L2LA::canRadiate ( const Event& state, pair<int,int> ints,
540  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
541  return ( state[ints.first].isFinal()
542  && state[ints.first].isLepton() && state[ints.first].isCharged()
543  && state[ints.second].isCharged()
544  && bools["doQEDshowerByL"]);
545 }
546 
547 bool Dire_fsr_qed_L2LA::canRadiate ( const Event& state, int iRadBef,
548  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
549  return ( state[iRadBef].isFinal()
550  && state[iRadBef].isLepton() && state[iRadBef].isCharged()
551  && state[iRecBef].isCharged()
552  && doQEDshowerByL);
553 }
554 
555 int Dire_fsr_qed_L2LA::kinMap() { return 1;}
556 int Dire_fsr_qed_L2LA::motherID(int idDaughter) { return idDaughter;}
557 int Dire_fsr_qed_L2LA::sisterID(int) { return 22;}
558 
559 double Dire_fsr_qed_L2LA::gaugeFactor ( int idRadBef, int idRecBef) {
560  double chgRad = particleDataPtr->charge(idRadBef);
561  double chgRec = particleDataPtr->charge(idRecBef);
562  double charge = -1.*chgRad*chgRec;
563  if (!splitInfo.radBef()->isFinal) charge *= -1.;
564  if (!splitInfo.recBef()->isFinal) charge *= -1.;
565  if (idRadBef != 0 && idRecBef != 0) return charge;
566  // Set probability to zero.
567  return 0.;
568 }
569 
570 double Dire_fsr_qed_L2LA::symmetryFactor ( int, int ) { return 1.;}
571 
572 int Dire_fsr_qed_L2LA::radBefID(int idRad, int idEmt) {
573  if (idEmt == 22 && particleDataPtr->isLepton(idRad)
574  && particleDataPtr->charge(idRad) != 0) return idRad;
575  return 0;
576 }
577 
578 pair<int,int> Dire_fsr_qed_L2LA::radBefCols(
579  int colRadAfter, int acolRadAfter,
580  int, int) { return make_pair(colRadAfter,acolRadAfter); }
581 
582 vector<int>Dire_fsr_qed_L2LA::recPositions(const Event& state, int iRad,
583  int iEmt) {
584 
585  vector<int> recs;
586  if ( !state[iRad].isFinal()
587  || !(state[iRad].isLepton() && state[iRad].isCharged())
588  || state[iEmt].id() != 22) return recs;
589 
590  // Particles to exclude as recoilers.
591  vector<int> iExc(createvector<int>(iRad)(iEmt));
592  // Find charged particles.
593  for (int i=0; i < state.size(); ++i) {
594  if ( find(iExc.begin(), iExc.end(), i) != iExc.end() ) continue;
595  if ( state[i].isCharged() ) {
596  if (state[i].isFinal())
597  recs.push_back(i);
598  if (state[i].mother1() == 1 && state[i].mother2() == 0)
599  recs.push_back(i);
600  if (state[i].mother1() == 2 && state[i].mother2() == 0)
601  recs.push_back(i);
602  }
603  }
604  // Done.
605  return recs;
606 
607 }
608 
609 // Pick z for new splitting.
610 double Dire_fsr_qed_L2LA::zSplit(double zMinAbs, double, double m2dip) {
611  double Rz = rndmPtr->flat();
612  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/m2dip;
613  double p = pow( 1. + pow2(1-zMinAbs)/kappa2, Rz );
614  double res = 1. - sqrt( p - 1. )*sqrt(kappa2);
615  return res;
616 }
617 
618 // New overestimates, z-integrated versions.
619 double Dire_fsr_qed_L2LA::overestimateInt(double zMinAbs, double,
620  double, double m2dip, int) {
621  double wt = 0.;
622  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
623  double preFac = symmetryFactor() * abs(charge);
624  // Q -> QG, soft part (currently also used for collinear part).
625  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/m2dip;
626  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa2);
627  return wt;
628 }
629 
630 // Return overestimate for new splitting.
631 double Dire_fsr_qed_L2LA::overestimateDiff(double z, double m2dip, int) {
632  double wt = 0.;
633  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
634  double preFac = symmetryFactor() * abs(charge);
635  double kappaOld2 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/m2dip;
636  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappaOld2);
637  return wt;
638 }
639 
640 // Return kernel for new splitting.
641 bool Dire_fsr_qed_L2LA::calc(const Event& state, int orderNow) {
642 
643  // Dummy statement to avoid compiler warnings.
644  if (false) cout << state[0].e() << orderNow << endl;
645 
646  // Read all splitting variables.
647  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
648  m2dip(splitInfo.kinematics()->m2Dip),
649  m2RadBef(splitInfo.kinematics()->m2RadBef),
650  m2Rad(splitInfo.kinematics()->m2RadAft),
651  m2Rec(splitInfo.kinematics()->m2Rec),
652  m2Emt(splitInfo.kinematics()->m2EmtAft);
653  int splitType(splitInfo.type);
654 
655  // Calculate kernel.
656  // Note: We are calculating the z <--> 1-z symmetrised kernel here,
657  // and later multiply with z to project out Q->QQ,
658  // i.e. the gluon is soft and the quark is identified.
659  double wt = 0.;
660 
661  double chargeFac = gaugeFactor(splitInfo.radBef()->id,
662  splitInfo.recBef()->id);
663  vector <int> in, out;
664  for (int i=0; i < state.size(); ++i) {
665  if (state[i].isFinal()) out.push_back(state[i].id());
666  if (state[i].mother1() == 1 && state[i].mother2() == 0)
667  in.push_back(state[i].id());
668  if (state[i].mother1() == 2 && state[i].mother2() == 0)
669  in.push_back(state[i].id());
670  }
671  out.push_back(22);
672  bool hasME = pT2 > pow2(settingsPtr->parm("Dire:pTminMECs"))
673  && doMECs && fsr->weights->hasME(in,out);
674  if (hasME && chargeFac < 0.0) chargeFac = abs(chargeFac);
675 
676  if ( doForcePos
677  && (chargeFac < 0. || splitInfo.radBef()->id != splitInfo.recBef()->id)
678  && (hasME || pT2 > pT2minForcePos)) chargeFac = chgprefac*abs(chargeFac);
679 
680  double preFac = symmetryFactor() * chargeFac;
681  double kappa2 = pT2/m2dip;
682  wt = preFac * ( 2. * z * (1.-z) / ( pow2(1.-z) + kappa2) );
683 
684  // Correction for massive splittings.
685  bool doMassive = (abs(splitType) == 2);
686 
687  // Add collinear term for massless splittings.
688  if (!doMassive && orderNow >= 0) wt += preFac * ( 1.-z );
689 
690  // Add collinear term for massive splittings.
691  if (doMassive && orderNow >= 0) {
692 
693  double pipj = 0., vijkt = 1., vijk = 1.;
694 
695  // splitType == 2 -> Massive FF
696  if (splitType == 2) {
697 
698  // Calculate CS variables.
699  double yCS = kappa2 / (1.-z);
700  double nu2RadBef = m2RadBef/m2dip;
701  double nu2Rad = m2Rad/m2dip;
702  double nu2Emt = m2Emt/m2dip;
703  double nu2Rec = m2Rec/m2dip;
704  vijk = pow2(1.-yCS) - 4.*(yCS+nu2Rad+nu2Emt)*nu2Rec;
705  double Q2mass = m2dip + m2Rad + m2Rec + m2Emt;
706  vijkt = pow2(Q2mass/m2dip - nu2RadBef - nu2Rec)
707  - 4.*nu2RadBef*nu2Rec;
708  vijk = sqrt(vijk) / (1-yCS);
709  vijkt = sqrt(vijkt)/ (Q2mass/m2dip - nu2RadBef - nu2Rec);
710  pipj = m2dip * yCS/2.;
711  // splitType ==-2 -> Massive FI
712  } else if (splitType ==-2) {
713 
714  // Calculate CS variables.
715  double xCS = 1 - kappa2/(1.-z);
716  vijk = 1.;
717  vijkt = 1.;
718  pipj = m2dip/2. * (1-xCS)/xCS;
719  }
720 
721  // Add B1 for massive splittings.
722  double massCorr = vijkt/vijk*( 1. - z - m2RadBef/pipj);
723  wt += preFac*massCorr;
724 
725  }
726 
727  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
728 
729  // Now multiply with z to project out Q->QG,
730  // i.e. the gluon is soft and the quark is identified.
731  wt *= z;
732 
733  // Trivial map of values, since kernel does not depend on coupling.
734  unordered_map<string,double> wts;
735  wts.insert( make_pair("base", wt ));
736  if (doVariations) {
737  // Create muR-variations.
738  if (settingsPtr->parm("Variations:muRfsrDown") != 1.)
739  wts.insert( make_pair("Variations:muRfsrDown", wt ));
740  if (settingsPtr->parm("Variations:muRfsrUp") != 1.)
741  wts.insert( make_pair("Variations:muRfsrUp", wt ));
742  }
743 
744  // Store kernel values.
745  clearKernels();
746  for ( unordered_map<string,double>::iterator it = wts.begin();
747  it != wts.end(); ++it )
748  kernelVals.insert(make_pair( it->first, it->second ));
749 
750  return true;
751 
752 }
753 
754 //==========================================================================
755 
756 // Class inheriting from SplittingQED class.
757 
758 // SplittingQED function Q->GQ (FSR)
759 // At leading order, this can be combined with Q->QG because of symmetry. Since
760 // this is no longer possible at NLO, we keep the kernels separately.
761 
762 // Return true if this kernel should partake in the evolution.
763 bool Dire_fsr_qed_L2AL::canRadiate ( const Event& state, pair<int,int> ints,
764  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
765  return ( state[ints.first].isFinal()
766  && state[ints.first].isLepton() && state[ints.first].isCharged()
767 // && state[ints.first].isLepton() && state[ints.second].isLepton()
768  && state[ints.second].isCharged()
769  && bools["doQEDshowerByL"]);
770 }
771 
772 bool Dire_fsr_qed_L2AL::canRadiate ( const Event& state, int iRadBef,
773  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
774  return ( state[iRadBef].isFinal()
775  && state[iRadBef].isLepton() && state[iRadBef].isCharged()
776  && state[iRecBef].isCharged()
777  && doQEDshowerByL);
778 }
779 
780 int Dire_fsr_qed_L2AL::kinMap() { return 1;}
781 int Dire_fsr_qed_L2AL::motherID(int idDaughter) { return idDaughter;}
782 int Dire_fsr_qed_L2AL::sisterID(int) { return 22;}
783 
784 double Dire_fsr_qed_L2AL::gaugeFactor ( int idRadBef, int idRecBef) {
785  double chgRad = particleDataPtr->charge(idRadBef);
786  double chgRec = particleDataPtr->charge(idRecBef);
787  double charge = -1.*chgRad*chgRec;
788  if (!splitInfo.radBef()->isFinal) charge *= -1.;
789  if (!splitInfo.recBef()->isFinal) charge *= -1.;
790  if (idRadBef != 0 && idRecBef != 0) return charge;
791  // Set probability to zero.
792  return 0.;
793 }
794 
795 double Dire_fsr_qed_L2AL::symmetryFactor ( int, int ) { return 1.;}
796 
797 int Dire_fsr_qed_L2AL::radBefID(int idRad, int idEmt) {
798  if (idRad == 22 && particleDataPtr->isLepton(idEmt)
799  && particleDataPtr->charge(idEmt) != 0) return idEmt;
800  if (idEmt == 22 && particleDataPtr->isLepton(idRad)
801  && particleDataPtr->charge(idRad) != 0) return idRad;
802  return 0;
803 }
804 
805 pair<int,int> Dire_fsr_qed_L2AL::radBefCols(
806  int colRadAfter, int acolRadAfter,
807  int, int) { return make_pair(colRadAfter,acolRadAfter); }
808 
809 vector<int>Dire_fsr_qed_L2AL::recPositions(const Event& state, int iRad,
810  int iEmt) {
811 
812  vector<int> recs;
813  if ( !state[iRad].isFinal()
814  || !(state[iRad].isLepton() && state[iRad].isCharged())
815  || state[iEmt].id() != 22) return recs;
816 
817  // Particles to exclude as recoilers.
818  vector<int> iExc(createvector<int>(iRad)(iEmt));
819  // Find charged particles.
820  for (int i=0; i < state.size(); ++i) {
821  if ( find(iExc.begin(), iExc.end(), i) != iExc.end() ) continue;
822  if ( state[i].isCharged() ) {
823  if (state[i].isFinal())
824  recs.push_back(i);
825  if (state[i].mother1() == 1 && state[i].mother2() == 0)
826  recs.push_back(i);
827  if (state[i].mother1() == 2 && state[i].mother2() == 0)
828  recs.push_back(i);
829  }
830  }
831  // Done.
832  return recs;
833 
834 }
835 
836 // Pick z for new splitting.
837 double Dire_fsr_qed_L2AL::zSplit(double zMinAbs, double, double m2dip) {
838  double Rz = rndmPtr->flat();
839  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/m2dip;
840  double p = pow( 1. + pow2(1-zMinAbs)/kappa2, Rz );
841  double res = 1. - sqrt( p - 1. )*sqrt(kappa2);
842  return res;
843 }
844 
845 // New overestimates, z-integrated versions.
846 double Dire_fsr_qed_L2AL::overestimateInt(double zMinAbs, double,
847  double, double m2dip, int) {
848  double wt = 0.;
849  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
850  double preFac = symmetryFactor() * abs(charge);
851  // Q -> QG, soft part (currently also used for collinear part).
852  double kappa2 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/m2dip;
853  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa2);
854  return wt;
855 }
856 
857 // Return overestimate for new splitting.
858 double Dire_fsr_qed_L2AL::overestimateDiff(double z, double m2dip, int) {
859  double wt = 0.;
860  double charge = gaugeFactor(splitInfo.radBef()->id, splitInfo.recBef()->id);
861  double preFac = symmetryFactor() * abs(charge);
862  double kappaOld2 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/m2dip;
863  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappaOld2);
864  return wt;
865 }
866 
867 // Return kernel for new splitting.
868 bool Dire_fsr_qed_L2AL::calc(const Event& state, int orderNow) {
869 
870  // Dummy statement to avoid compiler warnings.
871  if (false) cout << state[0].e() << orderNow << endl;
872 
873  // Read all splitting variables.
874  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
875  m2dip(splitInfo.kinematics()->m2Dip),
876  m2RadBef(splitInfo.kinematics()->m2RadBef),
877  m2Rad(splitInfo.kinematics()->m2RadAft),
878  m2Rec(splitInfo.kinematics()->m2Rec),
879  m2Emt(splitInfo.kinematics()->m2EmtAft);
880  int splitType(splitInfo.type);
881 
882  // Calculate kernel.
883  // Note: We are calculating the z <--> 1-z symmetrised kernel here,
884  // and later multiply with 1-z to project out Q->GQ,
885  // i.e. the quark is soft and the gluon is identified.
886  double wt = 0.;
887 
888  double chargeFac = gaugeFactor(splitInfo.radBef()->id,
889  splitInfo.recBef()->id);
890  vector <int> in, out;
891  for (int i=0; i < state.size(); ++i) {
892  if (state[i].isFinal()) out.push_back(state[i].id());
893  if (state[i].mother1() == 1 && state[i].mother2() == 0)
894  in.push_back(state[i].id());
895  if (state[i].mother1() == 2 && state[i].mother2() == 0)
896  in.push_back(state[i].id());
897  }
898  out.push_back(22);
899  bool hasME = pT2 > pow2(settingsPtr->parm("Dire:pTminMECs"))
900  && doMECs && fsr->weights->hasME(in,out);
901  if (hasME && chargeFac < 0.0) chargeFac = abs(chargeFac);
902 
903  if ( doForcePos
904  && (chargeFac < 0. || splitInfo.radBef()->id != splitInfo.recBef()->id)
905  && (hasME || pT2 > pT2minForcePos)) chargeFac = chgprefac*abs(chargeFac);
906 
907  double preFac = symmetryFactor() * chargeFac;
908  double kappa2 = pT2/m2dip;
909  wt = preFac * ( 2. * z * (1.-z) / ( pow2(1.-z) + kappa2) );
910 
911  // Correction for massive splittings.
912  bool doMassive = (abs(splitType) == 2);
913 
914  // Add collinear term for massless splittings.
915  if (!doMassive && orderNow >= 0) wt += preFac * ( 1. - z );
916 
917  // Add collinear term for massive splittings.
918  if (doMassive && orderNow >= 0) {
919 
920  double pipj = 0., vijkt = 1., vijk = 1.;
921 
922  // splitType == 2 -> Massive FF
923  if (splitType == 2) {
924 
925  // Calculate CS variables.
926  double yCS = kappa2 / (1.-z);
927  double nu2RadBef = m2RadBef/m2dip;
928  double nu2Rad = m2Rad/m2dip;
929  double nu2Emt = m2Emt/m2dip;
930  double nu2Rec = m2Rec/m2dip;
931  vijk = pow2(1.-yCS) - 4.*(yCS+nu2Rad+nu2Emt)*nu2Rec;
932  double Q2mass = m2dip + m2Rad + m2Rec + m2Emt;
933  vijkt = pow2(Q2mass/m2dip - nu2RadBef - nu2Rec)
934  - 4.*nu2RadBef*nu2Rec;
935  vijk = sqrt(vijk) / (1-yCS);
936  vijkt = sqrt(vijkt)/ (Q2mass/m2dip - nu2RadBef - nu2Rec);
937  pipj = m2dip * yCS/2.;
938  // splitType ==-2 -> Massive FI
939  } else if (splitType ==-2) {
940 
941  // Calculate CS variables.
942  double xCS = 1 - kappa2/(1.-z);
943  vijk = 1.;
944  vijkt = 1.;
945  pipj = m2dip/2. * (1-xCS)/xCS;
946  }
947 
948  // Add B1 for massive splittings.
949  double massCorr = vijkt/vijk*( 1. - z - m2RadBef/pipj);
950  wt += preFac*massCorr;
951 
952  }
953 
954  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
955 
956  // Now multiply with (1-z) to project out Q->GQ,
957  // i.e. the quark is soft and the gluon is identified.
958  wt *= ( 1. - z );
959 
960  // Trivial map of values, since kernel does not depend on coupling.
961  unordered_map<string,double> wts;
962  wts.insert( make_pair("base", wt ));
963  if (doVariations) {
964  // Create muR-variations.
965  if (settingsPtr->parm("Variations:muRfsrDown") != 1.)
966  wts.insert( make_pair("Variations:muRfsrDown", wt ));
967  if (settingsPtr->parm("Variations:muRfsrUp") != 1.)
968  wts.insert( make_pair("Variations:muRfsrUp", wt ));
969  }
970 
971  // Store kernel values.
972  clearKernels();
973  for ( unordered_map<string,double>::iterator it = wts.begin();
974  it != wts.end(); ++it )
975  kernelVals.insert(make_pair( it->first, it->second ));
976 
977  return true;
978 
979 }
980 
981 //==========================================================================
982 
983 // Class inheriting from SplittingQED class.
984 
985 // SplittingQED function Q->QG (ISR)
986 
987 // Return true if this kernel should partake in the evolution.
988 bool Dire_isr_qed_Q2QA::canRadiate ( const Event& state, pair<int,int> ints,
989  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
990  return (!state[ints.first].isFinal()
991  && state[ints.first].isQuark()
992  && state[ints.second].isCharged()
993  && bools["doQEDshowerByQ"] );
994 }
995 
996 bool Dire_isr_qed_Q2QA::canRadiate ( const Event& state, int iRadBef,
997  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
998  return (!state[iRadBef].isFinal()
999  && state[iRadBef].isQuark()
1000  && state[iRecBef].isCharged()
1001  && doQEDshowerByQ);
1002 }
1003 
1004 int Dire_isr_qed_Q2QA::kinMap() { return 1;}
1005 int Dire_isr_qed_Q2QA::motherID(int idDaughter) { return idDaughter;}
1006 int Dire_isr_qed_Q2QA::sisterID(int) { return 22;}
1007 
1008 double Dire_isr_qed_Q2QA::gaugeFactor ( int idRadBef, int idRecBef) {
1009  double chgRad = particleDataPtr->charge(idRadBef);
1010  double chgRec = particleDataPtr->charge(idRecBef);
1011  double charge = -1.*chgRad*chgRec;
1012  if (!splitInfo.radBef()->isFinal) charge *= -1.;
1013  if (!splitInfo.recBef()->isFinal) charge *= -1.;
1014  if (idRadBef != 0 && idRecBef != 0) return charge;
1015  // Set probability to zero.
1016  return 0.;
1017 }
1018 
1019 double Dire_isr_qed_Q2QA::symmetryFactor ( int, int ) { return 1.;}
1020 
1021 int Dire_isr_qed_Q2QA::radBefID(int idRad, int idEmt) {
1022  if (particleDataPtr->isQuark(idRad) && idEmt == 22 ) return idRad;
1023  return 0;
1024 }
1025 
1026 pair<int,int> Dire_isr_qed_Q2QA::radBefCols( int colRadAfter, int acolRadAfter,
1027  int, int) {
1028  bool isQuark = (colRadAfter > 0);
1029  if (isQuark) return make_pair(colRadAfter,0);
1030  return make_pair(0,acolRadAfter);
1031 }
1032 
1033 vector<int>Dire_isr_qed_Q2QA::recPositions(const Event& state, int iRad,
1034  int iEmt) {
1035 
1036  vector<int> recs;
1037  if ( state[iRad].isFinal() || !state[iRad].isQuark()
1038  || state[iEmt].id() != 22) return recs;
1039 
1040  // Particles to exclude as recoilers.
1041  vector<int> iExc(createvector<int>(iRad)(iEmt));
1042  // Find charged particles.
1043  for (int i=0; i < state.size(); ++i) {
1044  if ( find(iExc.begin(), iExc.end(), i) != iExc.end() ) continue;
1045  if ( state[i].isCharged() ) {
1046  if (state[i].isFinal())
1047  recs.push_back(i);
1048  if (state[i].mother1() == 1 && state[i].mother2() == 0)
1049  recs.push_back(i);
1050  if (state[i].mother1() == 2 && state[i].mother2() == 0)
1051  recs.push_back(i);
1052  }
1053  }
1054  // Done.
1055  return recs;
1056 
1057 }
1058 
1059 // Pick z for new splitting.
1060 double Dire_isr_qed_Q2QA::zSplit(double zMinAbs, double, double m2dip) {
1061  double Rz = rndmPtr->flat();
1062  double kappa2 = pow2(settingsPtr->parm("SpaceShower:pTminChgQ"))/m2dip;
1063  double p = pow( 1. + pow2(1-zMinAbs)/kappa2, Rz );
1064  double res = 1. - sqrt( p - 1. )*sqrt(kappa2);
1065  return res;
1066 }
1067 
1068 // New overestimates, z-integrated versions.
1069 double Dire_isr_qed_Q2QA::overestimateInt(double zMinAbs, double,
1070  double, double m2dip, int) {
1071  double wt = 0.;
1072  double preFac = symmetryFactor() * abs(gaugeFactor(splitInfo.radBef()->id,
1073  splitInfo.recBef()->id));
1074  double kappa2 = pow2(settingsPtr->parm("SpaceShower:pTminChgQ"))/m2dip;
1075  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa2);
1076 
1077  return wt;
1078 }
1079 
1080 // Return overestimate for new splitting.
1081 double Dire_isr_qed_Q2QA::overestimateDiff(double z, double m2dip, int) {
1082  double wt = 0.;
1083  double preFac = symmetryFactor() * abs(gaugeFactor(splitInfo.radBef()->id,
1084  splitInfo.recBef()->id));
1085  double kappaOld2 = pow2(settingsPtr->parm("SpaceShower:pTminChgQ"))/m2dip;
1086  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappaOld2);
1087  return wt;
1088 }
1089 
1090 // Return kernel for new splitting.
1091 bool Dire_isr_qed_Q2QA::calc(const Event& state, int orderNow) {
1092 
1093  // Dummy statement to avoid compiler warnings.
1094  if (false) cout << state[0].e() << orderNow << endl;
1095 
1096  // Read all splitting variables.
1097  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
1098  m2dip(splitInfo.kinematics()->m2Dip);
1099 
1100  double wt = 0.;
1101 
1102  double chargeFac = gaugeFactor(splitInfo.radBef()->id,
1103  splitInfo.recBef()->id);
1104  vector <int> in, out;
1105  for (int i=0; i < state.size(); ++i) {
1106  if (state[i].isFinal()) out.push_back(state[i].id());
1107  if (state[i].mother1() == 1 && state[i].mother2() == 0)
1108  in.push_back(state[i].id());
1109  if (state[i].mother1() == 2 && state[i].mother2() == 0)
1110  in.push_back(state[i].id());
1111  }
1112  out.push_back(22);
1113  bool hasME = pT2 > pow2(settingsPtr->parm("Dire:pTminMECs"))
1114  && doMECs && isr->weights->hasME(in,out);
1115  if (hasME && chargeFac < 0.0) chargeFac = abs(chargeFac);
1116 
1117  if ( doForcePos
1118  && (chargeFac < 0. || splitInfo.radBef()->id != splitInfo.recBef()->id)
1119  && (hasME || pT2 > pT2minForcePos)) chargeFac = chgprefac*abs(chargeFac);
1120 
1121  double preFac = symmetryFactor() * chargeFac;
1122  double kappa2 = pT2/m2dip;
1123  wt = preFac * ( 2. * z * (1.-z) / ( pow2(1.-z) + kappa2) );
1124  if (orderNow >= 0) wt += preFac * (1.-z);
1125  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
1126 
1127  // Trivial map of values, since kernel does not depend on coupling.
1128  unordered_map<string,double> wts;
1129  wts.insert( make_pair("base", wt ));
1130  if (doVariations) {
1131  // Create muR-variations.
1132  if (settingsPtr->parm("Variations:muRisrDown") != 1.)
1133  wts.insert( make_pair("Variations:muRisrDown", wt ));
1134  if (settingsPtr->parm("Variations:muRisrUp") != 1.)
1135  wts.insert( make_pair("Variations:muRisrUp", wt ));
1136  }
1137 
1138  // Store kernel values.
1139  clearKernels();
1140  for ( unordered_map<string,double>::iterator it = wts.begin();
1141  it != wts.end(); ++it )
1142  kernelVals.insert(make_pair( it->first, it->second ));
1143 
1144  return true;
1145 
1146 }
1147 
1148 //==========================================================================
1149 
1150 // Class inheriting from SplittingQED class.
1151 
1152 // SplittingQED function G->QQ (ISR)
1153 
1154 // Return true if this kernel should partake in the evolution.
1155 bool Dire_isr_qed_A2QQ::canRadiate ( const Event& state, pair<int,int> ints,
1156  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1157  return (!state[ints.first].isFinal()
1158  && state[ints.first].isQuark()
1159  && bools["doQEDshowerByQ"] );
1160 }
1161 
1162 bool Dire_isr_qed_A2QQ::canRadiate ( const Event& state, int iRadBef,
1163  int, Settings*, PartonSystems*, BeamParticle*){
1164  return (!state[iRadBef].isFinal()
1165  && state[iRadBef].isQuark()
1166  && doQEDshowerByQ);
1167 }
1168 
1169 int Dire_isr_qed_A2QQ::kinMap() { return 1;}
1170 int Dire_isr_qed_A2QQ::motherID(int) { return 22;}
1171 int Dire_isr_qed_A2QQ::sisterID(int idDaughter) { return -idDaughter;}
1172 double Dire_isr_qed_A2QQ::gaugeFactor ( int, int ) { return 1.;}
1173 double Dire_isr_qed_A2QQ::symmetryFactor ( int, int ) { return 1.;}
1174 
1175 int Dire_isr_qed_A2QQ::radBefID(int, int idEA){ return -idEA;}
1176 pair<int,int> Dire_isr_qed_A2QQ::radBefCols( int, int, int colEmtAfter,
1177  int acolEmtAfter) {
1178  if ( acolEmtAfter > 0 ) return make_pair(acolEmtAfter,0);
1179  return make_pair(0, colEmtAfter);
1180 }
1181 
1182 // Pick z for new splitting.
1183 double Dire_isr_qed_A2QQ::zSplit(double zMinAbs, double zMaxAbs, double) {
1184  // Note: Combined with PDF ratio, flat overestimate performs
1185  // better than using the full splitting kernel as overestimate.
1186  double res = zMinAbs + rndmPtr->flat() * (zMaxAbs - zMinAbs);
1187  return res;
1188 }
1189 
1190 // New overestimates, z-integrated versions.
1191 double Dire_isr_qed_A2QQ::overestimateInt(double zMinAbs, double zMaxAbs,
1192  double, double, int) {
1193  double wt = 0.;
1194  double preFac = symmetryFactor() * gaugeFactor();
1195  // Note: Combined with PDF ratio, flat overestimate performs
1196  // better than using the full splitting kernel as overestimate.
1197  wt = enhance * preFac
1198  * 2. * ( zMaxAbs - zMinAbs);
1199  return wt;
1200 }
1201 
1202 // Return overestimate for new splitting.
1203 double Dire_isr_qed_A2QQ::overestimateDiff(double, double, int) {
1204  double wt = 0.;
1205  double preFac = symmetryFactor() * gaugeFactor();
1206  // Note: Combined with PDF ratio, flat overestimate performs
1207  // better than using the full splitting kernel as overestimate.
1208  wt = enhance * preFac
1209  * 2.;
1210  return wt;
1211 }
1212 
1213 // Return kernel for new splitting.
1214 bool Dire_isr_qed_A2QQ::calc(const Event& state, int orderNow) {
1215 
1216  // Dummy statement to avoid compiler warnings.
1217  if (false) cout << state[0].e() << orderNow << endl;
1218 
1219  // Read all splitting variables.
1220  double z(splitInfo.kinematics()->z);
1221 
1222  double wt = 0.;
1223  double preFac = symmetryFactor() * gaugeFactor();
1224  wt = preFac * (pow(1.-z,2.) + pow(z,2.));
1225 
1226  if (orderNow >= 0) wt = 0.;
1227 
1228  // Trivial map of values, since kernel does not depend on coupling.
1229  unordered_map<string,double> wts;
1230  wts.insert( make_pair("base", wt ));
1231  if (doVariations) {
1232  // Create muR-variations.
1233  if (settingsPtr->parm("Variations:muRisrDown") != 1.)
1234  wts.insert( make_pair("Variations:muRisrDown", wt ));
1235  if (settingsPtr->parm("Variations:muRisrUp") != 1.)
1236  wts.insert( make_pair("Variations:muRisrUp", wt ));
1237  }
1238 
1239  // Store kernel values.
1240  clearKernels();
1241  for ( unordered_map<string,double>::iterator it = wts.begin();
1242  it != wts.end(); ++it )
1243  kernelVals.insert(make_pair( it->first, it->second ));
1244 
1245  return true;
1246 
1247 }
1248 
1249 //==========================================================================
1250 
1251 // Class inheriting from SplittingQED class.
1252 
1253 // SplittingQED function Q->AQ (ISR)
1254 
1255 // Return true if this kernel should partake in the evolution.
1256 bool Dire_isr_qed_Q2AQ::canRadiate ( const Event& state, pair<int,int> ints,
1257  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1258  return (!state[ints.first].isFinal()
1259  && state[ints.first].id() == 22
1260  && bools["doQEDshowerByQ"] );
1261 }
1262 
1263 bool Dire_isr_qed_Q2AQ::canRadiate ( const Event& state, int iRadBef,
1264  int, Settings*, PartonSystems*, BeamParticle*){
1265  return (!state[iRadBef].isFinal()
1266  && state[iRadBef].id() == 22
1267  && doQEDshowerByQ);
1268 }
1269 
1270 int Dire_isr_qed_Q2AQ::kinMap() { return 1;}
1271 int Dire_isr_qed_Q2AQ::motherID(int) { return 1;} // Use 1 as dummy
1272 int Dire_isr_qed_Q2AQ::sisterID(int) { return 1;} // Use 1 as dummy
1273 double Dire_isr_qed_Q2AQ::gaugeFactor ( int, int ) { return 1.;}
1274 double Dire_isr_qed_Q2AQ::symmetryFactor ( int, int ) { return 0.5;}
1275 
1276 int Dire_isr_qed_Q2AQ::radBefID(int, int){ return 22;}
1277 pair<int,int> Dire_isr_qed_Q2AQ::radBefCols( int, int, int, int) {
1278  return make_pair(0,0); }
1279 
1280 // Pick z for new splitting.
1281 double Dire_isr_qed_Q2AQ::zSplit(double zMinAbs, double, double) {
1282  double R = rndmPtr->flat();
1283  double res = pow(zMinAbs,3./4.)
1284  / ( pow(1. + R*(-1. + pow(zMinAbs,-3./8.)),2./3.)
1285  *pow(R - (-1. + R)*pow(zMinAbs,3./8.),2.));
1286  return res;
1287 }
1288 
1289 // New overestimates, z-integrated versions.
1290 double Dire_isr_qed_Q2AQ::overestimateInt(double zMinAbs, double,
1291  double, double, int) {
1292  double wt = 0.;
1293  double preFac = symmetryFactor() * gaugeFactor();
1294  wt = enhance * preFac * 2./3. * (8.*(-1. + pow(zMinAbs,-3./8.)));
1295  return wt;
1296 }
1297 
1298 // Return overestimate for new splitting.
1299 double Dire_isr_qed_Q2AQ::overestimateDiff(double z, double, int) {
1300  double wt = 0.;
1301  double preFac = symmetryFactor() * gaugeFactor();
1302  wt = enhance * preFac * 2. / pow(z,11./8.);
1303  return wt;
1304 }
1305 
1306 // Return kernel for new splitting.
1307 bool Dire_isr_qed_Q2AQ::calc(const Event& state, int orderNow) {
1308 
1309  // Dummy statement to avoid compiler warnings.
1310  if (false) cout << state[0].e() << orderNow << endl;
1311 
1312  // Read all splitting variables.
1313  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
1314  m2dip(splitInfo.kinematics()->m2Dip),
1315  m2Rec(splitInfo.kinematics()->m2Rec);
1316  int splitType(splitInfo.type);
1317 
1318  double wt = 0.;
1319  double preFac = symmetryFactor() * gaugeFactor();
1320  double kappa2 = pT2 / m2dip;
1321  wt = preFac * 2.*z*(1.-z) / (pow2(z)+kappa2);
1322  if (orderNow >= 0) wt += preFac*z;
1323 
1324  // Correction for massive IF splittings.
1325  bool doMassive = ( m2Rec > 0. && splitType == 2);
1326 
1327  if (doMassive && orderNow >= 0) {
1328  // Construct CS variables.
1329  double uCS = kappa2 / (1-z);
1330 
1331  double massCorr = -2. * m2Rec / m2dip * uCS / (1.-uCS);
1332  // Add correction.
1333  wt += preFac * massCorr;
1334 
1335  }
1336 
1337  // Trivial map of values, since kernel does not depend on coupling.
1338  unordered_map<string,double> wts;
1339  wts.insert( make_pair("base", wt ));
1340  if (doVariations) {
1341  // Create muR-variations.
1342  if (settingsPtr->parm("Variations:muRisrDown") != 1.)
1343  wts.insert( make_pair("Variations:muRisrDown", wt ));
1344  if (settingsPtr->parm("Variations:muRisrUp") != 1.)
1345  wts.insert( make_pair("Variations:muRisrUp", wt ));
1346  }
1347 
1348  // Store kernel values.
1349  clearKernels();
1350  for ( unordered_map<string,double>::iterator it = wts.begin();
1351  it != wts.end(); ++it )
1352  kernelVals.insert(make_pair( it->first, it->second ));
1353 
1354  return true;
1355 
1356 }
1357 
1358 //==========================================================================
1359 
1360 // Class inheriting from SplittingQED class.
1361 
1362 // SplittingQED function L->LA (ISR)
1363 
1364 // Return true if this kernel should partake in the evolution.
1365 bool Dire_isr_qed_L2LA::canRadiate ( const Event& state, pair<int,int> ints,
1366  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1367  return (!state[ints.first].isFinal()
1368  && state[ints.first].isLepton() && state[ints.first].isCharged()
1369  && state[ints.second].isCharged()
1370  && bools["doQEDshowerByL"]);
1371 }
1372 
1373 bool Dire_isr_qed_L2LA::canRadiate ( const Event& state, int iRadBef,
1374  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
1375  return (!state[iRadBef].isFinal()
1376  && state[iRadBef].isLepton()
1377  && state[iRadBef].isCharged()
1378  && state[iRecBef].isCharged()
1379  && doQEDshowerByL);
1380 }
1381 
1382 int Dire_isr_qed_L2LA::kinMap() { return 1;}
1383 int Dire_isr_qed_L2LA::motherID(int idDaughter) { return idDaughter;}
1384 int Dire_isr_qed_L2LA::sisterID(int) { return 22;}
1385 
1386 double Dire_isr_qed_L2LA::gaugeFactor ( int idRadBef, int idRecBef) {
1387  double chgRad = particleDataPtr->charge(idRadBef);
1388  double chgRec = particleDataPtr->charge(idRecBef);
1389  double charge = -1.*chgRad*chgRec;
1390  if (!splitInfo.radBef()->isFinal) charge *= -1.;
1391  if (!splitInfo.recBef()->isFinal) charge *= -1.;
1392  if (idRadBef != 0 && idRecBef != 0) return charge;
1393  // Set probability to zero.
1394  return 0.;
1395 }
1396 
1397 double Dire_isr_qed_L2LA::symmetryFactor ( int, int ) { return 1.;}
1398 
1399 int Dire_isr_qed_L2LA::radBefID(int idRad, int idEmt) {
1400  if (particleDataPtr->isLepton(idRad) && particleDataPtr->charge(idRad) != 0
1401  && idEmt == 22 ) return idRad;
1402  return 0;
1403 }
1404 
1405 
1406 pair<int,int> Dire_isr_qed_L2LA::radBefCols( int colRadAfter, int acolRadAfter,
1407  int, int) {
1408  bool isQuark = (colRadAfter > 0);
1409  if (isQuark) return make_pair(colRadAfter,0);
1410  return make_pair(0,acolRadAfter);
1411 }
1412 
1413 vector<int>Dire_isr_qed_L2LA::recPositions(const Event& state, int iRad,
1414  int iEmt) {
1415 
1416  vector<int> recs;
1417  if ( state[iRad].isFinal()
1418  || !(state[iRad].isLepton() && state[iRad].isCharged())
1419  || state[iEmt].id() != 22) return recs;
1420 
1421  // Particles to exclude as recoilers.
1422  vector<int> iExc(createvector<int>(iRad)(iEmt));
1423  // Find charged particles.
1424  for (int i=0; i < state.size(); ++i) {
1425  if ( find(iExc.begin(), iExc.end(), i) != iExc.end() ) continue;
1426  if ( state[i].isCharged() ) {
1427  if (state[i].isFinal())
1428  recs.push_back(i);
1429  if (state[i].mother1() == 1 && state[i].mother2() == 0)
1430  recs.push_back(i);
1431  if (state[i].mother1() == 2 && state[i].mother2() == 0)
1432  recs.push_back(i);
1433  }
1434  }
1435  // Done.
1436  return recs;
1437 
1438 }
1439 
1440 // Pick z for new splitting.
1441 double Dire_isr_qed_L2LA::zSplit(double zMinAbs, double, double m2dip) {
1442  double Rz = rndmPtr->flat();
1443  double kappa2 = pow2(settingsPtr->parm("SpaceShower:pTminChgL"))/m2dip;
1444  double p = pow( 1. + pow2(1-zMinAbs)/kappa2, Rz );
1445  double res = 1. - sqrt( p - 1. )*sqrt(kappa2);
1446  return res;
1447 }
1448 
1449 // New overestimates, z-integrated versions.
1450 double Dire_isr_qed_L2LA::overestimateInt(double zMinAbs, double,
1451  double, double m2dip, int) {
1452  double wt = 0.;
1453  double preFac = symmetryFactor() * abs(gaugeFactor(splitInfo.radBef()->id,
1454  splitInfo.recBef()->id));
1455  double kappa2 = pow2(settingsPtr->parm("SpaceShower:pTminChgL"))/m2dip;
1456  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa2);
1457  return wt;
1458 }
1459 
1460 // Return overestimate for new splitting.
1461 double Dire_isr_qed_L2LA::overestimateDiff(double z, double m2dip, int) {
1462  double wt = 0.;
1463  double preFac = symmetryFactor() * abs(gaugeFactor(splitInfo.radBef()->id,
1464  splitInfo.recBef()->id));
1465  double kappaOld2 = pow2(settingsPtr->parm("SpaceShower:pTminChgL"))/m2dip;
1466  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappaOld2);
1467  return wt;
1468 }
1469 
1470 // Return kernel for new splitting.
1471 bool Dire_isr_qed_L2LA::calc(const Event& state, int orderNow) {
1472 
1473  // Dummy statement to avoid compiler warnings.
1474  if (false) cout << state[0].e() << orderNow << endl;
1475 
1476  // Read all splitting variables.
1477  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
1478  m2dip(splitInfo.kinematics()->m2Dip);
1479 
1480  double wt = 0.;
1481 
1482  double chargeFac = gaugeFactor(splitInfo.radBef()->id,
1483  splitInfo.recBef()->id);
1484  vector <int> in, out;
1485  for (int i=0; i < state.size(); ++i) {
1486  if (state[i].isFinal()) out.push_back(state[i].id());
1487  if (state[i].mother1() == 1 && state[i].mother2() == 0)
1488  in.push_back(state[i].id());
1489  if (state[i].mother1() == 2 && state[i].mother2() == 0)
1490  in.push_back(state[i].id());
1491  }
1492  out.push_back(22);
1493  bool hasME = pT2 > pow2(settingsPtr->parm("Dire:pTminMECs"))
1494  && doMECs && isr->weights->hasME(in,out);
1495  if (hasME && chargeFac < 0.0) chargeFac = abs(chargeFac);
1496 
1497  if ( doForcePos
1498  && (chargeFac < 0. || splitInfo.radBef()->id != splitInfo.recBef()->id)
1499  && (hasME || pT2 > pT2minForcePos)) chargeFac = chgprefac*abs(chargeFac);
1500 
1501  double preFac = symmetryFactor() * chargeFac;
1502  double kappa2 = pT2/m2dip;
1503  wt = preFac * ( 2. * z * (1.-z) / ( pow2(1.-z) + kappa2) );
1504  if (orderNow >= 0) wt += preFac * (1.-z);
1505 
1506  // Trivial map of values, since kernel does not depend on coupling.
1507  unordered_map<string,double> wts;
1508  wts.insert( make_pair("base", wt ));
1509  if (doVariations) {
1510  // Create muR-variations.
1511  if (settingsPtr->parm("Variations:muRisrDown") != 1.)
1512  wts.insert( make_pair("Variations:muRisrDown", wt ));
1513  if (settingsPtr->parm("Variations:muRisrUp") != 1.)
1514  wts.insert( make_pair("Variations:muRisrUp", wt ));
1515  }
1516 
1517  // Store kernel values.
1518  clearKernels();
1519  for ( unordered_map<string,double>::iterator it = wts.begin();
1520  it != wts.end(); ++it )
1521  kernelVals.insert(make_pair( it->first, it->second ));
1522 
1523  return true;
1524 
1525 }
1526 
1527 //==========================================================================
1528 
1529 // Class inheriting from SplittingQED class.
1530 
1531 // SplittingQED function A->LL (ISR)
1532 
1533 // Return true if this kernel should partake in the evolution.
1534 bool Dire_isr_qed_A2LL::canRadiate ( const Event& state, pair<int,int> ints,
1535  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1536  return (!state[ints.first].isFinal()
1537  && state[ints.first].isLepton() && state[ints.first].isCharged()
1538  && bools["doQEDshowerByL"]);
1539 }
1540 
1541 bool Dire_isr_qed_A2LL::canRadiate ( const Event& state, int iRadBef,
1542  int, Settings*, PartonSystems*, BeamParticle*){
1543  return (!state[iRadBef].isFinal()
1544  && state[iRadBef].isLepton()
1545  && state[iRadBef].isCharged()
1546  && doQEDshowerByL);
1547 }
1548 
1549 int Dire_isr_qed_A2LL::kinMap() { return 1;}
1550 int Dire_isr_qed_A2LL::motherID(int) { return 22;}
1551 int Dire_isr_qed_A2LL::sisterID(int idDaughter) { return -idDaughter;}
1552 double Dire_isr_qed_A2LL::gaugeFactor ( int, int ) { return 1.;}
1553 double Dire_isr_qed_A2LL::symmetryFactor ( int, int ) { return 1.;}
1554 
1555 int Dire_isr_qed_A2LL::radBefID(int, int idEA){ return -idEA;}
1556 pair<int,int> Dire_isr_qed_A2LL::radBefCols( int, int, int colEmtAfter,
1557  int acolEmtAfter) {
1558  if ( acolEmtAfter > 0 ) return make_pair(acolEmtAfter,0);
1559  return make_pair(0, colEmtAfter);
1560 }
1561 
1562 // Pick z for new splitting.
1563 double Dire_isr_qed_A2LL::zSplit(double zMinAbs, double zMaxAbs, double) {
1564  // Note: Combined with PDF ratio, flat overestimate performs
1565  // better than using the full splitting kernel as overestimate.
1566  double res = zMinAbs + rndmPtr->flat() * (zMaxAbs - zMinAbs);
1567  return res;
1568 }
1569 
1570 // New overestimates, z-integrated versions.
1571 double Dire_isr_qed_A2LL::overestimateInt(double zMinAbs, double zMaxAbs,
1572  double, double, int) {
1573  double wt = 0.;
1574  double preFac = symmetryFactor() * gaugeFactor();
1575  // Note: Combined with PDF ratio, flat overestimate performs
1576  // better than using the full splitting kernel as overestimate.
1577  wt = enhance * preFac
1578  * 2. * ( zMaxAbs - zMinAbs);
1579 
1580  return wt;
1581 }
1582 
1583 // Return overestimate for new splitting.
1584 double Dire_isr_qed_A2LL::overestimateDiff(double, double, int) {
1585  double wt = 0.;
1586  double preFac = symmetryFactor() * gaugeFactor();
1587  // Note: Combined with PDF ratio, flat overestimate performs
1588  // better than using the full splitting kernel as overestimate.
1589  wt = enhance * preFac
1590  * 2.;
1591  return wt;
1592 }
1593 
1594 // Return kernel for new splitting.
1595 bool Dire_isr_qed_A2LL::calc(const Event& state, int orderNow) {
1596 
1597  // Dummy statement to avoid compiler warnings.
1598  if (false) cout << state[0].e() << orderNow << endl;
1599 
1600  // Read all splitting variables.
1601  double z(splitInfo.kinematics()->z);
1602 
1603  double wt = 0.;
1604  double preFac = symmetryFactor() * gaugeFactor();
1605  wt = preFac * (pow(1.-z,2.) + pow(z,2.));
1606 
1607  if (orderNow == -1) wt = 0.0;
1608 
1609  // Trivial map of values, since kernel does not depend on coupling.
1610  unordered_map<string,double> wts;
1611  wts.insert( make_pair("base", wt ));
1612  if (doVariations) {
1613  // Create muR-variations.
1614  if (settingsPtr->parm("Variations:muRisrDown") != 1.)
1615  wts.insert( make_pair("Variations:muRisrDown", wt ));
1616  if (settingsPtr->parm("Variations:muRisrUp") != 1.)
1617  wts.insert( make_pair("Variations:muRisrUp", wt ));
1618  }
1619 
1620  // Store kernel values.
1621  clearKernels();
1622  for ( unordered_map<string,double>::iterator it = wts.begin();
1623  it != wts.end(); ++it )
1624  kernelVals.insert(make_pair( it->first, it->second ));
1625 
1626  return true;
1627 
1628 }
1629 
1630 //==========================================================================
1631 
1632 // Class inheriting from SplittingQED class.
1633 
1634 // SplittingQED function L->AL (ISR)
1635 
1636 // Return true if this kernel should partake in the evolution.
1637 bool Dire_isr_qed_L2AL::canRadiate ( const Event& state, pair<int,int> ints,
1638  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1639  return (!state[ints.first].isFinal()
1640  && state[ints.first].id() == 22
1641  && bools["doQEDshowerByL"]);
1642 }
1643 
1644 bool Dire_isr_qed_L2AL::canRadiate ( const Event& state, int iRadBef,
1645  int, Settings*, PartonSystems*, BeamParticle*){
1646  return (!state[iRadBef].isFinal()
1647  && state[iRadBef].id() == 22
1648  && doQEDshowerByL);
1649 }
1650 
1651 int Dire_isr_qed_L2AL::kinMap() { return 1;}
1652 int Dire_isr_qed_L2AL::motherID(int) { return 1;} // Use 1 as dummy
1653 int Dire_isr_qed_L2AL::sisterID(int) { return 1;} // Use 1 as dummy
1654 double Dire_isr_qed_L2AL::gaugeFactor ( int, int ) { return 1.;}
1655 double Dire_isr_qed_L2AL::symmetryFactor ( int, int ) { return 0.5;}
1656 
1657 int Dire_isr_qed_L2AL::radBefID(int, int){ return 22;}
1658 pair<int,int> Dire_isr_qed_L2AL::radBefCols( int, int, int, int) {
1659  return make_pair(0,0); }
1660 
1661 // Pick z for new splitting.
1662 double Dire_isr_qed_L2AL::zSplit(double zMinAbs, double, double) {
1663  double R = rndmPtr->flat();
1664  double res = pow(zMinAbs,3./4.)
1665  / ( pow(1. + R*(-1. + pow(zMinAbs,-3./8.)),2./3.)
1666  *pow(R - (-1. + R)*pow(zMinAbs,3./8.),2.));
1667  return res;
1668 }
1669 
1670 // New overestimates, z-integrated versions.
1671 double Dire_isr_qed_L2AL::overestimateInt(double zMinAbs, double,
1672  double, double, int) {
1673  double wt = 0.;
1674  double preFac = symmetryFactor() * gaugeFactor();
1675  wt = enhance * preFac * 2./3. * (8.*(-1. + pow(zMinAbs,-3./8.)));
1676  return wt;
1677 }
1678 
1679 // Return overestimate for new splitting.
1680 double Dire_isr_qed_L2AL::overestimateDiff(double z, double, int) {
1681  double wt = 0.;
1682  double preFac = symmetryFactor() * gaugeFactor();
1683  wt = enhance * preFac * 2. / pow(z,11./8.);
1684  return wt;
1685 }
1686 
1687 // Return kernel for new splitting.
1688 bool Dire_isr_qed_L2AL::calc(const Event& state, int orderNow) {
1689 
1690  // Dummy statement to avoid compiler warnings.
1691  if (false) cout << state[0].e() << orderNow << endl;
1692 
1693  // Read all splitting variables.
1694  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
1695  m2dip(splitInfo.kinematics()->m2Dip),
1696  m2Rec(splitInfo.kinematics()->m2Rec);
1697  int splitType(splitInfo.type);
1698 
1699  double wt = 0.;
1700  double preFac = symmetryFactor() * gaugeFactor();
1701  double kappa2 = pT2 / m2dip;
1702 
1703  wt = preFac * 2.*z*(1.-z) / (pow2(z)+kappa2);
1704  if (orderNow >= 0) wt += preFac*z;
1705 
1706  // Correction for massive IF splittings.
1707  bool doMassive = ( m2Rec > 0. && splitType == 2);
1708 
1709  if (doMassive && orderNow >= 0) {
1710  // Construct CS variables.
1711  double uCS = kappa2 / (1-z);
1712 
1713  double massCorr = -2. * m2Rec / m2dip * uCS / (1.-uCS);
1714  // Add correction.
1715  wt += preFac * massCorr;
1716 
1717  }
1718 
1719  // Trivial map of values, since kernel does not depend on coupling.
1720  unordered_map<string,double> wts;
1721  wts.insert( make_pair("base", wt ));
1722  if (doVariations) {
1723  // Create muR-variations.
1724  if (settingsPtr->parm("Variations:muRisrDown") != 1.)
1725  wts.insert( make_pair("Variations:muRisrDown", wt ));
1726  if (settingsPtr->parm("Variations:muRisrUp") != 1.)
1727  wts.insert( make_pair("Variations:muRisrUp", wt ));
1728  }
1729 
1730  // Store kernel values.
1731  clearKernels();
1732  for ( unordered_map<string,double>::iterator it = wts.begin();
1733  it != wts.end(); ++it )
1734  kernelVals.insert(make_pair( it->first, it->second ));
1735 
1736  return true;
1737 
1738 }
1739 
1740 //==========================================================================
1741 
1742 // Return true if this kernel should partake in the evolution.
1743 bool Dire_fsr_qed_Q2QA_notPartial::canRadiate ( const Event& state,
1744  pair<int,int> ints,
1745  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1746  return ( state[ints.first].isFinal()
1747  && state[ints.first].isQuark() && !state[ints.second].isCharged()
1748  && bools["doQEDshowerByQ"]);
1749 }
1750 
1751 bool Dire_fsr_qed_Q2QA_notPartial::canRadiate (const Event& state, int iRadBef,
1752  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
1753  return ( state[iRadBef].isFinal()
1754  && state[iRadBef].isQuark() && !state[iRecBef].isCharged()
1755  && doQEDshowerByQ);
1756 }
1757 
1758 int Dire_fsr_qed_Q2QA_notPartial::kinMap() {return 1;}
1759 int Dire_fsr_qed_Q2QA_notPartial::motherID(int idDaughter) {return idDaughter;}
1760 int Dire_fsr_qed_Q2QA_notPartial::sisterID(int) {return 22;}
1761 
1762 double Dire_fsr_qed_Q2QA_notPartial::gaugeFactor ( int idRadBef, int) {
1763  if (idRadBef != 0) return pow2(particleDataPtr->charge(idRadBef));
1764  // Set probability to zero.
1765  return 0.;
1766 }
1767 
1768 double Dire_fsr_qed_Q2QA_notPartial::symmetryFactor ( int, int ) { return 1.;}
1769 
1770 int Dire_fsr_qed_Q2QA_notPartial::radBefID(int idRad, int idEmt) {
1771  if (particleDataPtr->isQuark(idRad) && idEmt == 22 ) return idRad;
1772  return 0;
1773 }
1774 
1775 pair<int,int> Dire_fsr_qed_Q2QA_notPartial::radBefCols(
1776  int colRadAfter, int acolRadAfter,
1777  int, int) { return make_pair(colRadAfter,acolRadAfter); }
1778 
1779 vector<int>Dire_fsr_qed_Q2QA_notPartial::recPositions(const Event&, int, int) {
1780  // Not implemented yet.
1781  vector<int> recs;
1782  return recs;
1783 }
1784 
1785 // Pick z for new splitting.
1786 double Dire_fsr_qed_Q2QA_notPartial::zSplit(double zMinAbs, double,
1787  double m2dip) {
1788  double Rz = rndmPtr->flat();
1789  double kappa4 = pow4(settingsPtr->parm("TimeShower:pTminChgQ"))/pow2(m2dip);
1790  double p = pow( 1. + pow2(1-zMinAbs)/kappa4, Rz );
1791  double res = 1. - sqrt( p - 1. )*sqrt(kappa4);
1792  return res;
1793 }
1794 
1795 // New overestimates, z-integrated versions.
1796 double Dire_fsr_qed_Q2QA_notPartial::overestimateInt(double zMinAbs, double,
1797  double, double m2dip, int) {
1798  double wt = 0.;
1799  double charge = gaugeFactor(splitInfo.radBef()->id);
1800  double preFac = symmetryFactor() * abs(charge);
1801  // Q -> QG, soft part (currently also used for collinear part).
1802  double kappa4 = pow4(settingsPtr->parm("TimeShower:pTminChgQ"))/pow2(m2dip);
1803  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa4);
1804  return wt;
1805 }
1806 
1807 // Return overestimate for new splitting.
1808 double Dire_fsr_qed_Q2QA_notPartial::overestimateDiff(double z, double m2dip,
1809  int) {
1810  double wt = 0.;
1811  double charge = gaugeFactor(splitInfo.radBef()->id);
1812  double preFac = symmetryFactor() * abs(charge);
1813  double kappa4 = pow2(settingsPtr->parm("TimeShower:pTminChgQ"))/pow2(m2dip);
1814  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappa4);
1815  return wt;
1816 }
1817 
1818 // Return kernel for new splitting.
1819 bool Dire_fsr_qed_Q2QA_notPartial::calc(const Event& state, int orderNow) {
1820 
1821  // Dummy statement to avoid compiler warnings.
1822  if (false) cout << state[0].e() << orderNow << endl;
1823 
1824  // Read all splitting variables.
1825  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
1826  m2dip(splitInfo.kinematics()->m2Dip),
1827  m2RadBef(splitInfo.kinematics()->m2RadBef),
1828  m2Rad(splitInfo.kinematics()->m2RadAft),
1829  m2Rec(splitInfo.kinematics()->m2Rec),
1830  m2Emt(splitInfo.kinematics()->m2EmtAft);
1831  int splitType(splitInfo.type);
1832 
1833  // Calculate kernel.
1834  // Note: We are calculating the z <--> 1-z symmetrised kernel here,
1835  // and later multiply with z to project out Q->QQ,
1836  // i.e. the gluon is soft and the quark is identified.
1837  double wt = 0.;
1838 
1839  double chargeFac = gaugeFactor(splitInfo.radBef()->id);
1840 
1841  double preFac = symmetryFactor() * chargeFac;
1842  double kappa2 = pT2/m2dip;
1843  wt = preFac * 2. * z / (1.-z);
1844 
1845  // Correction for massive splittings.
1846  bool doMassive = (abs(splitType) == 2);
1847 
1848  // Add collinear term for massless splittings.
1849  if (!doMassive && orderNow >= 0) wt += preFac * ( 1.-z );
1850 
1851  // Add collinear term for massive splittings.
1852  if (doMassive && orderNow >= 0) {
1853 
1854  double pipj = 0., vijkt = 1., vijk = 1.;
1855 
1856  // splitType == 2 -> Massive FF
1857  if (splitType == 2) {
1858 
1859  // Calculate CS variables.
1860  double yCS = kappa2 / (1.-z);
1861  double nu2RadBef = m2RadBef/m2dip;
1862  double nu2Rad = m2Rad/m2dip;
1863  double nu2Emt = m2Emt/m2dip;
1864  double nu2Rec = m2Rec/m2dip;
1865  vijk = pow2(1.-yCS) - 4.*(yCS+nu2Rad+nu2Emt)*nu2Rec;
1866  double Q2mass = m2dip + m2Rad + m2Rec + m2Emt;
1867  vijkt = pow2(Q2mass/m2dip - nu2RadBef - nu2Rec)
1868  - 4.*nu2RadBef*nu2Rec;
1869  vijk = sqrt(vijk) / (1-yCS);
1870  vijkt = sqrt(vijkt)/ (Q2mass/m2dip - nu2RadBef - nu2Rec);
1871  pipj = m2dip * yCS/2.;
1872  // splitType ==-2 -> Massive FI
1873  } else if (splitType ==-2) {
1874 
1875  // Calculate CS variables.
1876  double xCS = 1 - kappa2/(1.-z);
1877  vijk = 1.;
1878  vijkt = 1.;
1879  pipj = m2dip/2. * (1-xCS)/xCS;
1880  }
1881 
1882  // Add B1 for massive splittings.
1883  double massCorr = vijkt/vijk*( 1. - z - m2RadBef/pipj);
1884  wt += preFac*massCorr;
1885  }
1886 
1887  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
1888 
1889  // Trivial map of values, since kernel does not depend on coupling.
1890  unordered_map<string,double> wts;
1891  wts.insert( make_pair("base", wt ));
1892  if (doVariations) {
1893  // Create muR-variations.
1894  if (settingsPtr->parm("Variations:muRfsrDown") != 1.)
1895  wts.insert( make_pair("Variations:muRfsrDown", wt ));
1896  if (settingsPtr->parm("Variations:muRfsrUp") != 1.)
1897  wts.insert( make_pair("Variations:muRfsrUp", wt ));
1898  }
1899 
1900  // Store kernel values.
1901  clearKernels();
1902  for ( unordered_map<string,double>::iterator it = wts.begin();
1903  it != wts.end(); ++it )
1904  kernelVals.insert(make_pair( it->first, it->second ));
1905 
1906  return true;
1907 
1908 }
1909 
1910 //==========================================================================
1911 
1912 // Class inheriting from SplittingQED class.
1913 
1914 // SplittingQED function Q->QG (FSR)
1915 
1916 // Return true if this kernel should partake in the evolution.
1917 bool Dire_fsr_qed_L2LA_notPartial::canRadiate ( const Event& state,
1918  pair<int,int> ints,
1919  unordered_map<string,bool> bools, Settings*, PartonSystems*, BeamParticle*){
1920  return ( state[ints.first].isFinal()
1921  && state[ints.first].isLepton() && state[ints.first].isCharged()
1922  && !state[ints.second].isCharged() && bools["doQEDshowerByL"]);
1923 }
1924 
1925 bool Dire_fsr_qed_L2LA_notPartial::canRadiate (const Event& state, int iRadBef,
1926  int iRecBef, Settings*, PartonSystems*, BeamParticle*){
1927  return ( state[iRadBef].isFinal()
1928  && state[iRadBef].isLepton() && state[iRadBef].isCharged()
1929  && !state[iRecBef].isCharged() && doQEDshowerByL);
1930 }
1931 
1932 int Dire_fsr_qed_L2LA_notPartial::kinMap() {return 1;}
1933 int Dire_fsr_qed_L2LA_notPartial::motherID(int idDaughter) {return idDaughter;}
1934 int Dire_fsr_qed_L2LA_notPartial::sisterID(int) {return 22;}
1935 
1936 double Dire_fsr_qed_L2LA_notPartial::gaugeFactor (int idRadBef, int) {
1937  if (idRadBef != 0) return pow2(particleDataPtr->charge(idRadBef));
1938  return 0.;
1939 }
1940 
1941 double Dire_fsr_qed_L2LA_notPartial::symmetryFactor ( int, int ) { return 1.;}
1942 
1943 int Dire_fsr_qed_L2LA_notPartial::radBefID(int idRad, int idEmt) {
1944  if (idEmt == 22 && particleDataPtr->isLepton(idRad)
1945  && particleDataPtr->charge(idRad) != 0) return idRad;
1946  return 0;
1947 }
1948 
1949 
1950 pair<int,int> Dire_fsr_qed_L2LA_notPartial::radBefCols(
1951  int colRadAfter, int acolRadAfter,
1952  int, int) { return make_pair(colRadAfter,acolRadAfter); }
1953 
1954 vector<int>Dire_fsr_qed_L2LA_notPartial::recPositions(const Event&, int, int) {
1955  // Not implemented yet.
1956  vector<int> recs;
1957  return recs;
1958 }
1959 
1960 // Pick z for new splitting.
1961 double Dire_fsr_qed_L2LA_notPartial::zSplit(double zMinAbs, double,
1962  double m2dip) {
1963  double Rz = rndmPtr->flat();
1964  double kappa4 = pow4(settingsPtr->parm("TimeShower:pTminChgL"))/pow2(m2dip);
1965  double p = pow( 1. + pow2(1-zMinAbs)/kappa4, Rz );
1966  double res = 1. - sqrt( p - 1. )*sqrt(kappa4);
1967  return res;
1968 }
1969 
1970 // New overestimates, z-integrated versions.
1971 double Dire_fsr_qed_L2LA_notPartial::overestimateInt(double zMinAbs, double,
1972  double, double m2dip, int) {
1973  double wt = 0.;
1974  double charge = gaugeFactor(splitInfo.radBef()->id);
1975  double preFac = symmetryFactor() * abs(charge);
1976  // Q -> QG, soft part (currently also used for collinear part).
1977  double kappa4 = pow4(settingsPtr->parm("TimeShower:pTminChgL"))/pow2(m2dip);
1978  wt = enhance * preFac * 2. * 0.5 * log( 1. + pow2(1.-zMinAbs)/kappa4);
1979  return wt;
1980 }
1981 
1982 // Return overestimate for new splitting.
1983 double Dire_fsr_qed_L2LA_notPartial::overestimateDiff(double z, double m2dip,
1984  int) {
1985  double wt = 0.;
1986  double charge = gaugeFactor(splitInfo.radBef()->id);
1987  double preFac = symmetryFactor() * abs(charge);
1988  double kappa4 = pow2(settingsPtr->parm("TimeShower:pTminChgL"))/pow2(m2dip);
1989  wt = enhance * preFac * 2.* (1.-z) / ( pow2(1.-z) + kappa4);
1990  return wt;
1991 }
1992 
1993 // Return kernel for new splitting.
1994 bool Dire_fsr_qed_L2LA_notPartial::calc(const Event& state, int orderNow) {
1995 
1996  // Dummy statement to avoid compiler warnings.
1997  if (false) cout << state[0].e() << orderNow << endl;
1998 
1999  // Read all splitting variables.
2000  double z(splitInfo.kinematics()->z), pT2(splitInfo.kinematics()->pT2),
2001  m2dip(splitInfo.kinematics()->m2Dip),
2002  m2RadBef(splitInfo.kinematics()->m2RadBef),
2003  m2Rad(splitInfo.kinematics()->m2RadAft),
2004  m2Rec(splitInfo.kinematics()->m2Rec),
2005  m2Emt(splitInfo.kinematics()->m2EmtAft);
2006  int splitType(splitInfo.type);
2007 
2008  // Calculate kernel.
2009  // Note: We are calculating the z <--> 1-z symmetrised kernel here,
2010  // and later multiply with z to project out Q->QQ,
2011  // i.e. the gluon is soft and the quark is identified.
2012  double wt = 0.;
2013 
2014  double chargeFac = gaugeFactor(splitInfo.radBef()->id);
2015 
2016  double preFac = symmetryFactor() * chargeFac;
2017  double kappa2 = pT2/m2dip;
2018  wt = preFac * 2. * z / (1.-z);
2019 
2020  // Correction for massive splittings.
2021  bool doMassive = (abs(splitType) == 2);
2022 
2023  // Add collinear term for massless splittings.
2024  if (!doMassive && orderNow >= 0) wt += preFac * ( 1.-z );
2025 
2026  // Add collinear term for massive splittings.
2027  if (doMassive && orderNow >= 0) {
2028 
2029  double pipj = 0., vijkt = 1., vijk = 1.;
2030 
2031  // splitType == 2 -> Massive FF
2032  if (splitType == 2) {
2033 
2034  // Calculate CS variables.
2035  double yCS = kappa2 / (1.-z);
2036  double nu2RadBef = m2RadBef/m2dip;
2037  double nu2Rad = m2Rad/m2dip;
2038  double nu2Emt = m2Emt/m2dip;
2039  double nu2Rec = m2Rec/m2dip;
2040  vijk = pow2(1.-yCS) - 4.*(yCS+nu2Rad+nu2Emt)*nu2Rec;
2041  double Q2mass = m2dip + m2Rad + m2Rec + m2Emt;
2042  vijkt = pow2(Q2mass/m2dip - nu2RadBef - nu2Rec)
2043  - 4.*nu2RadBef*nu2Rec;
2044  vijk = sqrt(vijk) / (1-yCS);
2045  vijkt = sqrt(vijkt)/ (Q2mass/m2dip - nu2RadBef - nu2Rec);
2046  pipj = m2dip * yCS/2.;
2047  // splitType ==-2 -> Massive FI
2048  } else if (splitType ==-2) {
2049 
2050  // Calculate CS variables.
2051  double xCS = 1 - kappa2/(1.-z);
2052  vijk = 1.;
2053  vijkt = 1.;
2054  pipj = m2dip/2. * (1-xCS)/xCS;
2055  }
2056 
2057  // Add B1 for massive splittings.
2058  double massCorr = vijkt/vijk*( 1. - z - m2RadBef/pipj);
2059  wt += preFac*massCorr;
2060  }
2061 
2062  if (orderNow < 0 && chargeFac < 0.) wt = 0.;
2063 
2064  // Trivial map of values, since kernel does not depend on coupling.
2065  unordered_map<string,double> wts;
2066  wts.insert( make_pair("base", wt ));
2067  if (doVariations) {
2068  // Create muR-variations.
2069  if (settingsPtr->parm("Variations:muRfsrDown") != 1.)
2070  wts.insert( make_pair("Variations:muRfsrDown", wt ));
2071  if (settingsPtr->parm("Variations:muRfsrUp") != 1.)
2072  wts.insert( make_pair("Variations:muRfsrUp", wt ));
2073  }
2074 
2075  // Store kernel values.
2076  clearKernels();
2077  for ( unordered_map<string,double>::iterator it = wts.begin();
2078  it != wts.end(); ++it )
2079  kernelVals.insert(make_pair( it->first, it->second ));
2080 
2081  return true;
2082 
2083 }
2084 
2085 //==========================================================================
2086 
2087 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5