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