StRoot  1
VinciaQED.cc
1 // VinciaQED.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Peter Skands, 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
7 // shower class and related auxiliary methods. Main author is Rob
8 // Verheyen.
9
10 #include "Pythia8/VinciaQED.h"
11
12 namespace Pythia8 {
13
14 //==========================================================================
15
16 // Class for the "Hungarian" pairing algorithm.
17
18 //--------------------------------------------------------------------------
19
20 // A single function wrapper for solving assignment problem.
21
22 double HungarianAlgorithm::solve(std::vector <std::vector<double> >&
23  distMatrix, std::vector<int>& assignment) {
24
25  unsigned int nRows = distMatrix.size();
26  unsigned int nCols = distMatrix[0].size();
27  double *distMatrixIn = new double[nRows * nCols];
28  int *solution = new int[nRows];
29  double cost = 0.0;
30
31  // Fill in the distMatrixIn. Mind the index is "i + nRows * j".
32  // Here the cost matrix of size MxN is defined as a double precision
33  // array of N*M elements. In the solving functions matrices are seen
34  // to be saved MATLAB-internally in row-order. (i.e. the matrix [1
35  // 2; 3 4] will be stored as a vector [1 3 2 4], NOT [1 2 3 4]).
36  for (unsigned int i = 0; i < nRows; i++)
37  for (unsigned int j = 0; j < nCols; j++)
38  distMatrixIn[i + nRows * j] = distMatrix[i][j];
39
40  // Call solving function.
41  optimal(solution, &cost, distMatrixIn, nRows, nCols);
42  assignment.clear();
43  for (unsigned int r = 0; r < nRows; r++)
44  assignment.push_back(solution[r]);
45  delete[] distMatrixIn;
46  delete[] solution;
47  return cost;
48
49 }
50
51 //--------------------------------------------------------------------------
52
53 // Solve optimal solution for assignment.
54
55 void HungarianAlgorithm::optimal(int *assignment, double *cost,
56  double *distMatrixIn, int nOfRows, int nOfColumns) {
57
58  // Initialization.
59  double *distMatrix, *distMatrixTemp, *distMatrixEnd, *columnEnd,
60  value, minValue;
61  bool *coveredColumns, *coveredRows, *starMatrix, *newStarMatrix,
62  *primeMatrix;
63  int nOfElements, minDim, row, col;
64  *cost = 0;
65  for (row = 0; row<nOfRows; row++) assignment[row] = -1;
66
67  // Generate working copy of distance matrix. Check if all matrix
68  // elements are positive.
69  nOfElements = nOfRows * nOfColumns;
70  distMatrix = (double *)malloc(nOfElements * sizeof(double));
71  distMatrixEnd = distMatrix + nOfElements;
72  for (row = 0; row<nOfElements; row++) {
73  value = distMatrixIn[row];
74  if (value < 0)
75  std::cerr << "HungarianAlgorithm::assigmentoptimal(): All"
76  << " matrix elements have to be non-negative." << std::endl;
77  distMatrix[row] = value;
78  }
79
80  // Memory allocation.
81  coveredColumns = (bool *)calloc(nOfColumns, sizeof(bool));
82  coveredRows = (bool *)calloc(nOfRows, sizeof(bool));
83  starMatrix = (bool *)calloc(nOfElements, sizeof(bool));
84  primeMatrix = (bool *)calloc(nOfElements, sizeof(bool));
85  // Used in step4.
86  newStarMatrix = (bool *)calloc(nOfElements, sizeof(bool));
87
88  // Preliminary steps.
89  if (nOfRows <= nOfColumns) {
90  minDim = nOfRows;
91  for (row = 0; row<nOfRows; row++) {
92  // Find the smallest element in the row.
93  distMatrixTemp = distMatrix + row;
94  minValue = *distMatrixTemp;
95  distMatrixTemp += nOfRows;
96  while (distMatrixTemp < distMatrixEnd) {
97  value = *distMatrixTemp;
98  if (value < minValue)
99  minValue = value;
100  distMatrixTemp += nOfRows;
101  }
102
103  // Subtract the smallest element from each element of the row.
104  distMatrixTemp = distMatrix + row;
105  while (distMatrixTemp < distMatrixEnd) {
106  *distMatrixTemp -= minValue;
107  distMatrixTemp += nOfRows;
108  }
109  }
110
111  // Steps 1 and 2a.
112  for (row = 0; row<nOfRows; row++)
113  for (col = 0; col<nOfColumns; col++)
114  if (abs(distMatrix[row + nOfRows*col]) < DBL_EPSILON)
115  if (!coveredColumns[col]) {
116  starMatrix[row + nOfRows*col] = true;
117  coveredColumns[col] = true;
118  break;
119  }
120  } else {
121  minDim = nOfColumns;
122  for (col = 0; col<nOfColumns; col++) {
123  // Find the smallest element in the column.
124  distMatrixTemp = distMatrix + nOfRows*col;
125  columnEnd = distMatrixTemp + nOfRows;
126  minValue = *distMatrixTemp++;
127  while (distMatrixTemp < columnEnd) {
128  value = *distMatrixTemp++;
129  if (value < minValue)
130  minValue = value;
131  }
132
133  // Subtract the smallest element from each element of the column.
134  distMatrixTemp = distMatrix + nOfRows*col;
135  while (distMatrixTemp < columnEnd)
136  *distMatrixTemp++ -= minValue;
137  }
138
139  // Steps 1 and 2a.
140  for (col = 0; col<nOfColumns; col++)
141  for (row = 0; row<nOfRows; row++)
142  if (abs(distMatrix[row + nOfRows*col]) < DBL_EPSILON)
143  if (!coveredRows[row]) {
144  starMatrix[row + nOfRows*col] = true;
145  coveredColumns[col] = true;
146  coveredRows[row] = true;
147  break;
148  }
149  for (row = 0; row<nOfRows; row++)
150  coveredRows[row] = false;
151  }
152
153  // Move to step 2b.
154  step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix,
155  coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
156
157  // Compute cost and remove invalid assignments.
158  calcCost(assignment, cost, distMatrixIn, nOfRows);
159
160  // Free allocated memory.
161  free(distMatrix);
162  free(coveredColumns);
163  free(coveredRows);
164  free(starMatrix);
165  free(primeMatrix);
166  free(newStarMatrix);
167  return;
168
169 }
170
171 //--------------------------------------------------------------------------
172
173 // Build the assignment vector.
174
175 void HungarianAlgorithm::vect(int *assignment,
176  bool *starMatrix, int nOfRows, int nOfColumns) {
177  int row, col;
178  for (row = 0; row<nOfRows; row++)
179  for (col = 0; col<nOfColumns; col++)
180  if (starMatrix[row + nOfRows*col]) {
181  assignment[row] = col;
182  break;
183  }
184 }
185
186 //--------------------------------------------------------------------------
187
188 // Calculate the assignment cost.
189
190 void HungarianAlgorithm::calcCost(int *assignment, double *cost,
191  double *distMatrix, int nOfRows) {
192  int row, col;
193  for (row = 0; row<nOfRows; row++) {
194  col = assignment[row];
195  if (col >= 0) *cost += distMatrix[row + nOfRows*col];
196  }
197 }
198
199 //--------------------------------------------------------------------------
200
201 // Factorized step 2a of the algorithm.
202
203 void HungarianAlgorithm::step2a(int *assignment, double *distMatrix,
204  bool *starMatrix, bool *newStarMatrix, bool *primeMatrix,
205  bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns,
206  int minDim) {
207
208  // Cover every column containing a starred zero.
209  bool *starMatrixTemp, *columnEnd;
210  int col;
211  for (col = 0; col<nOfColumns; col++) {
212  starMatrixTemp = starMatrix + nOfRows*col;
213  columnEnd = starMatrixTemp + nOfRows;
214  while (starMatrixTemp < columnEnd) {
215  if (*starMatrixTemp++) {
216  coveredColumns[col] = true;
217  break;
218  }
219  }
220  }
221
222  // Move to step 2b (note, original comment changed by Skands).
223  step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix,
224  coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
225
226 }
227
228 //--------------------------------------------------------------------------
229
230 // Factorized step 2b of the algorithm.
231
232 void HungarianAlgorithm::step2b(int *assignment, double *distMatrix,
233  bool *starMatrix, bool *newStarMatrix, bool *primeMatrix,
234  bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns,
235  int minDim) {
236
237  // Count covered columns.
238  int col, nOfCoveredColumns;
239  nOfCoveredColumns = 0;
240  for (col = 0; col<nOfColumns; col++)
241  if (coveredColumns[col]) nOfCoveredColumns++;
242
243  // Algorithm finished.
244  if (nOfCoveredColumns == minDim)
245  vect(assignment, starMatrix, nOfRows, nOfColumns);
246  // Move to step 3.
247  else step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix,
248  coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
249
250 }
251
252 //--------------------------------------------------------------------------
253
254 // Factorized step 3 of the algorithm.
255
256 void HungarianAlgorithm::step3(int *assignment, double *distMatrix,
257  bool *starMatrix, bool *newStarMatrix, bool *primeMatrix,
258  bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns,
259  int minDim) {
260
261  bool zerosFound;
262  int row, col, starCol;
263  zerosFound = true;
264  while (zerosFound) {
265  zerosFound = false;
266  for (col = 0; col<nOfColumns; col++)
267  if (!coveredColumns[col])
268  for (row = 0; row<nOfRows; row++)
269  if ((!coveredRows[row]) && (abs(distMatrix[row + nOfRows*col])
270  < DBL_EPSILON)) {
271  // Prime zero.
272  primeMatrix[row + nOfRows*col] = true;
273
274  // Find starred zero in current row.
275  for (starCol = 0; starCol<nOfColumns; starCol++)
276  if (starMatrix[row + nOfRows*starCol])
277  break;
278
279  // No starred zero found, move to step 4.
280  if (starCol == nOfColumns) {
281  step4(assignment, distMatrix, starMatrix, newStarMatrix,
282  primeMatrix, coveredColumns, coveredRows, nOfRows,
283  nOfColumns, minDim, row, col);
284  return;
285  } else {
286  coveredRows[row] = true;
287  coveredColumns[starCol] = false;
288  zerosFound = true;
289  break;
290  }
291  }
292  }
293
294  // Move to step 5.
295  step5(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix,
296  coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
297
298 }
299
300 //--------------------------------------------------------------------------
301
302 // Factorized step 4 of the algorithm.
303
304 void HungarianAlgorithm::step4(int *assignment, double *distMatrix,
305  bool *starMatrix, bool *newStarMatrix, bool *primeMatrix,
306  bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns,
307  int minDim, int row, int col) {
308
309  // Generate temporary copy of starMatrix.
310  int n, starRow, starCol, primeRow, primeCol;
311  int nOfElements = nOfRows*nOfColumns;
312  for (n = 0; n<nOfElements; n++) newStarMatrix[n] = starMatrix[n];
313  // Star current zero.
314  newStarMatrix[row + nOfRows*col] = true;
315  // Find starred zero in current column.
316  starCol = col;
317  for (starRow = 0; starRow<nOfRows; starRow++)
318  if (starMatrix[starRow + nOfRows*starCol])
319  break;
320  while (starRow < nOfRows) {
321  // Unstar the starred zero.
322  newStarMatrix[starRow + nOfRows*starCol] = false;
323  // Find primed zero in current row.
324  primeRow = starRow;
325  for (primeCol = 0; primeCol<nOfColumns; primeCol++)
326  if (primeMatrix[primeRow + nOfRows*primeCol])
327  break;
328  // Star the primed zero.
329  newStarMatrix[primeRow + nOfRows*primeCol] = true;
330  // Find starred zero in current column.
331  starCol = primeCol;
332  for (starRow = 0; starRow<nOfRows; starRow++)
333  if (starMatrix[starRow + nOfRows*starCol])
334  break;
335  }
336
337  // Use temporary copy as new starMatrix, delete all primes, uncover
338  // all rows.
339  for (n = 0; n<nOfElements; n++) {
340  primeMatrix[n] = false;
341  starMatrix[n] = newStarMatrix[n];
342  }
343  for (n = 0; n<nOfRows; n++) coveredRows[n] = false;
344
345  // Move to step 2a.
346  step2a(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix,
347  coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
348
349 }
350
351 //--------------------------------------------------------------------------
352
353 // Factorized step 5 of the algorithm.
354
355 void HungarianAlgorithm::step5(int *assignment, double *distMatrix,
356  bool *starMatrix, bool *newStarMatrix, bool *primeMatrix,
357  bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns,
358  int minDim) {
359
360  // Find smallest uncovered element h.
361  double h, value;
362  int row, col;
363  h = DBL_MAX;
364  for (row = 0; row<nOfRows; row++)
365  if (!coveredRows[row])
366  for (col = 0; col<nOfColumns; col++)
367  if (!coveredColumns[col]) {
368  value = distMatrix[row + nOfRows*col];
369  if (value < h)
370  h = value;
371  }
372
373  // Add h to each covered row.
374  for (row = 0; row<nOfRows; row++)
375  if (coveredRows[row])
376  for (col = 0; col<nOfColumns; col++)
377  distMatrix[row + nOfRows*col] += h;
378
379  // Subtract h from each uncovered column.
380  for (col = 0; col<nOfColumns; col++)
381  if (!coveredColumns[col])
382  for (row = 0; row<nOfRows; row++)
383  distMatrix[row + nOfRows*col] -= h;
384
385  // Move to step 3.
386  step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix,
387  coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
388
389 }
390
391 //==========================================================================
392
393 // Class for QED emissions.
394
395 //--------------------------------------------------------------------------
396
397 // Initialize the pointers.
398
399 void QEDemitElemental::initPtr(Rndm* rndmPtrIn,
400  PartonSystems* partonSystemsPtrIn) {
401  rndmPtr = rndmPtrIn;
402  partonSystemsPtr = partonSystemsPtrIn;
403  isInitPtr = true;
404 }
405
406 //--------------------------------------------------------------------------
407
408 // Initialize.
409
410 void QEDemitElemental::init(Event &event, int xIn, int yIn, double shhIn,
411  double verboseIn) {
412
413  if (!isInitPtr) printOut(__METHOD_NAME__, "initPtr not called");
414  x = xIn;
415  y = yIn;
416  shh = shhIn;
417  hasTrial = false;
418  isII = false;
419  isIF = false;
420  isFF = false;
421  isRF = false;
422  isIA = false;
423  isDip = false;
424
425  // If an II antenna, make sure x is the positive pz state.
426  if (!event[x].isFinal() && !event[y].isFinal() && event[x].pz() < 0)
427  swap(x,y);
428
429  // If an IF/RF antenna, make sure x is the initial state.
430  if (event[x].isFinal() && !event[y].isFinal()) swap(x,y);
431
432  // If a dipole, make sure x is the emitting object.
433  if (event[x].isFinal() && event[y].isFinal())
434  if (!event[x].isCharged() || event[y].isCharged()) swap(x,y);
435
436  idx = event[x].id();
437  idy = event[y].id();
438  mx2 = event[x].m2();
439  my2 = event[y].m2();
440  ex = event[x].e();
441  ey = event[y].e();
442  m2Ant = m2(event[x], event[y]);
443  sAnt = 2*dot4(event[x], event[y]);
444  QQ = - event[x].charge() * event[y].charge();
445
446  // II.
447  if (!event[x].isFinal() && !event[y].isFinal()) isII = true;
448
449  // IF/RF.
450  if (!event[x].isFinal() && event[y].isFinal()) {
451  // QQ is flipped for IF antennae.
452  QQ = -QQ;
453  // Check if initial state is in a beam.
454  int mother1 = event[x].mother1();
455  // Check if initial particle is A or B.
456  if (mother1 <= 2) {
457  isIF = true;
458  if (event[x].pz() > 0) isIA = true;
459  // Otherwise it's a resonance decay.
460  } else isRF = true;
461  }
462
463  // FF.
464  if (event[x].isFinal() && event[y].isFinal()) isFF = true;
465  isInit = true;
466  verbose = verboseIn;
467
468 }
469
470 //--------------------------------------------------------------------------
471
472 // Initialize.
473
474 void QEDemitElemental::init(Event &event, int xIn, vector<int> iRecoilIn,
475  double shhIn, double verboseIn) {
476
477  x = xIn;
478  iRecoil = iRecoilIn;
479  shh = shhIn;
480  hasTrial = false;
481  isII = false;
482  isIF = false;
483  isFF = false;
484  isRF = false;
485  isIA = false;
486  isDip = true;
487  idx = event[x].id();
488  mx2 = event[x].m2();
489
490  // Compute total recoiler momentum.
491  Vec4 pRecoil;
492  for (int i = 0; i < (int)iRecoil.size(); i++)
493  pRecoil += event[iRecoil[i]].p();
494  my2 = pRecoil.m2Calc();
495  m2Ant = (pRecoil + event[xIn].p()).m2Calc();
496  sAnt = 2*pRecoil*event[xIn].p();
497  QQ = 1;
498  isInit = true;
499  verbose = verboseIn;
500
501 }
502
503 //--------------------------------------------------------------------------
504
505 // Generate a trial point.
506
507 double QEDemitElemental::generateTrial(Event &event, double q2Start,
508  double q2Low, double alphaIn, double cIn) {
509
510  if (hasTrial) return q2Sav;
511  q2Sav = q2Low;
512  alpha = alphaIn;
513  c = cIn;
514
515  // FF.
516  if (isFF || isDip) {
518  q2Start = min(q2Start, sAnt/4.);
519  if (q2Start < q2Low) return q2Low;
520
521  // Compute phase space constants.
522  double lambda = m2Ant*m2Ant + mx2*mx2 + my2*my2
523  - 2.*m2Ant*mx2 - 2.*m2Ant*my2 - 2.*mx2*my2;
524  // zMin is identical for all instances.
525  double zMin = (4*q2Low/sAnt < 1E-8)
526  ? q2Low/sAnt
527  : 0.5*(1. - sqrt(1. - 4*q2Low/sAnt));
528
529  // Generate scale for eikonal piece.
530  if (true) {
531  double Iz = (zMin < 1E-8)
532  ? -2*log(zMin) - 2*zMin - pow2(zMin)
533  : 2*log((1-zMin)/zMin);
534  double comFac = 2*M_PI*sqrt(lambda)/alpha/Iz/c/sAnt;
535  double q2New = q2Start*pow(rndmPtr->flat(), comFac);
536  if (q2New > q2Sav) {
537  q2Sav = q2New;
538  zetaSav = 1/(exp(Iz*(0.5 - rndmPtr->flat())) + 1);
539  sxjSav = zetaSav != 1 ? sqrt(sAnt*q2Sav*zetaSav/((1-zetaSav))) : 0;
540  syjSav = zetaSav != 0 ? sqrt(sAnt*q2Sav*(1-zetaSav)/(zetaSav)) :
541  std::numeric_limits<double>::infinity();
542  }
543  }
544  // Generate scale for additional W piece on x.
545  if (isFF && abs(idx) == 24) {
546  double Iz = (zMin < 1E-8) ?
547  -log(zMin) - zMin - pow2(zMin)/2. : log((1-zMin)/zMin);
548  double comFac = 3.*M_PI*sqrt(lambda)/alpha/Iz/c/sAnt/2.;
549  double q2New = q2Start*pow(rndmPtr->flat(), comFac);
550  if (q2New > q2Sav) {
551  q2Sav = q2New;
552  double r = rndmPtr->flat();
553  zetaSav = (zMin < 1E-8)
554  ? 1 - pow(zMin,r)*(1. - (1.-r)*zMin)
555  : 1 - pow(zMin,r)*pow(1.-zMin, 1.-r);
556  sxjSav = q2Sav/zetaSav;
557  syjSav = zetaSav*sAnt;
558  }
559  }
560  // Generate scale for additional W piece on y.
561  if (isFF && abs(idy) == 24) {
562  double Iz = (zMin < 1E-8)
563  ? -log(zMin) - zMin - pow2(zMin)/2.
564  : log((1-zMin)/zMin);
565  double comFac = 3.*M_PI*sqrt(lambda)/alpha/Iz/c/sAnt/2.;
566  double q2New = q2Start*pow(rndmPtr->flat(), comFac);
567  if (q2New > q2Sav) {
568  q2Sav = q2New;
569  double r = rndmPtr->flat();
570  zetaSav = (zMin < 1E-8)
571  ? 1 - pow(zMin,r)*(1. - (1.-r)*zMin)
572  : 1 - pow(zMin,r)*pow(1.-zMin, 1.-r);
573  sxjSav = zetaSav*sAnt;
574  syjSav = q2Sav/zetaSav;
575  }
576  }
577  }
578
579  // IF.
580  if (isIF) {
581  // Compute exmax and sjkMax.
582  double exUsed = 0;
583  int nSys = partonSystemsPtr->sizeSys();
584  for (int i=0; i<nSys; i++) {
585  int iEv;
586  if (isIA) iEv = partonSystemsPtr->getInA(i);
587  else iEv = partonSystemsPtr->getInB(i);
588  exUsed += event[iEv].e();
589  }
590  double exMax = sqrt(shh)/2.0 - (exUsed-ex);
591  double sjkMax = sAnt*(exMax-ex)/ex;
592
594  q2Start = min(q2Start, sAnt*(exMax - ex)/ex);
595  if (q2Start < q2Low) return q2Low;
596  double zMax = sjkMax/(sjkMax + my2);
597  double zMin = q2Low/sjkMax;
598
599  // Check if there is any phase space available.
600  if (zMin < zMax) {
601  // Generate scale for eikonal piece.
602  if (true) {
603  double Iz = log(zMax/zMin);
604  double Rpdf = 1.;
605  double comFac = M_PI/alpha/Iz/c/Rpdf;
606  double q2New = q2Start*pow(rndmPtr->flat(), comFac);
607  if (q2New > q2Sav) {
608  q2Sav = q2New;
609  zetaSav = zMin*pow(zMax/zMin, rndmPtr->flat());
610  sxjSav = sAnt*zetaSav + q2Sav;
611  syjSav = q2Sav/zetaSav;
612  }
613  }
614
615  // Generate scale for additional W piece on y. The veto
616  // probability for this antenna piece includes an additional
617  // factor which is incorporated by a veto locally.
618  if (abs(idy) == 24) {
619  double Iz = log((1-zMin)/(1-zMax));
620  double Rpdf = 1.;
621  double comFac = 3.*M_PI/alpha/Iz/c/Rpdf/2.;
622  double q2New = q2Start;
623  double zetaNew, sxjNew, syjNew;
624  while (true) {
625  q2New *= pow(rndmPtr->flat(), comFac);
626  if (q2New < q2Sav) {break;}
627  zetaNew = 1. - (1-zMin)*pow((1-zMax)/(1-zMin),rndmPtr->flat());
628  sxjNew = sAnt*zetaNew + q2New;
629  syjNew = q2New/zetaNew;
630
631  // Veto probability.
632  double pVeto = sAnt/(sAnt + syjNew);
633  if (rndmPtr->flat() < pVeto) {
634  q2Sav = q2New;
635  zetaSav = zetaNew;
636  sxjSav = sxjNew;
637  syjSav = syjNew;
638  break;
639  }
640  }
641  }
642  }
643  }
644
645  // II.
646  if (isII) {
648  q2Start = min(q2Start, pow2(shh-sAnt)/shh/4.);
649  if (q2Start < q2Low) {return q2Low;}
650
651  // Generate scale for eikonal piece.
652  if (true) {
653  double zMin = 0.5*(shh-sAnt -
654  sqrt((shh-sAnt)*(shh-sAnt) - (4.*shh*q2Low)))/shh;
655  double zMax = 0.5*(shh-sAnt +
656  sqrt((shh-sAnt)*(shh-sAnt) - (4.*shh*q2Low)))/shh;
657  if (4.*shh*q2Low/pow2((shh-sAnt)) < 1e-8)
658  zMin = q2Low/(shh-sAnt);
659  double Iz = log(zMax*(1-zMin)/(1-zMax)/zMin);
660  double Rpdf = 1.;
661  double comFac = M_PI/alpha/Iz/c/Rpdf;
662  double q2New = q2Start*pow(rndmPtr->flat(), comFac);
663  if (q2New > q2Sav) {
664  q2Sav = q2New;
665  double r = rndmPtr->flat();
666  double w = pow(zMax/(1-zMax), r) * pow(zMin/(1-zMin), 1.-r);
667  zetaSav = w/(1.+w);
668  sxjSav = (q2Sav + sAnt*zetaSav)/(1.-zetaSav);
669  syjSav = q2Sav/zetaSav;
670  }
671  }
672  }
673
674  // RF.
675  if (isRF) {
676  // Compute phase space constants.
677  double mr2 = abs((event[x].p() - event[y].p()).m2Calc());
678  double mx = sqrt(mx2);
679  double my = sqrt(my2);
680  double mr = sqrt(mr2);
681  double lambda = mr2*mr2 + mx2*mx2 + my2*my2
682  - 2.*mr2*mx2 - 2.*mr2*my2 - 2.*mx2*my2;
683  double sjkMax = pow2(mx - mr) - my2;
684  double sajMax = mx2 - pow2(my + mr);
686  q2Start = min(q2Start, sajMax*sjkMax/(sAnt + sjkMax));
687
688  // Generate scale for eikonal piece.
689  if (true) {
690  double zMin = q2Low/sjkMax;
691  double zMax = sajMax/sAnt;
692  if(zMin < zMax){
693  double Iz = log(zMax/zMin);
694  double comFac = M_PI*sqrt(lambda)*sAnt/alpha/Iz/c/pow2(sAnt+sjkMax);
695  double q2New = q2Start;
696  double zetaNew, sxjNew, syjNew;
697  while (true) {
698  q2New *= pow(rndmPtr->flat(), comFac);
699  if (q2New < q2Sav) {break;}
700  zetaNew = zMin*pow(zMax/zMin, rndmPtr->flat());
701  sxjNew = sAnt*zetaNew + q2New;
702  syjNew = q2New/zetaNew;
703
704  // Veto probability.
705  double pVeto = pow2(syjNew+sAnt)/pow2(sjkMax+sAnt);
706  if (rndmPtr->flat() < pVeto) {
707  q2Sav = q2New;
708  zetaSav = zetaNew;
709  sxjSav = sxjNew;
710  syjSav = syjNew;
711  break;
712  }
713  }
714  }
715  }
716
717  // Generate scale for W in initial state.
718  if (abs(idx) == 24) {
719  double zMin = q2Low/(sajMax - q2Low);
720  double zMax = sjkMax/sAnt;
721  if(zMin < zMax && zMin > 0){
722  double Iz = pow2(zMax) + (1./3.)*pow3(zMax)
723  - pow2(zMin) - (1./3.)*pow3(zMin);
724  double comFac = 3.*M_PI*sqrt(lambda)/alpha/Iz/c/sAnt/2.;
725  double q2New = q2Start*pow(rndmPtr->flat(), comFac);
726
727  if (q2New > q2Sav) {
728  double a = rndmPtr->flat()*Iz + pow2(zMin) + (1./3.)*pow3(zMin);
729  // Solve for zeta using Newton-Raphson.
730  int n = 0;
731  zetaSav = zMin;
732  while(true) {
733  n++;
734  double f = pow2(zetaSav) + pow3(zetaSav)/3. - a;
735  double fPrime = 2.*zetaSav + pow2(zetaSav);
736  double zetaNew = zetaSav - f/fPrime;
737  if (zetaNew > zMax) {zetaSav = zMax; continue;}
738  if (zetaNew < zMin) {zetaSav = zMin; continue;}
739  if (abs(zetaNew - zetaSav) < 1E-8*zetaNew) {
740  zetaSav = zetaNew;
741  break;
742  }
743  if (n > 500) {
744  printOut(__METHOD_NAME__,
745  "RF(W) failed to find zeta with Newton-Raphson");
746  break;
747  }
748  zetaSav = zetaNew;
749  }
750  q2Sav = q2New;
751  sxjSav = (1.+zetaSav)*q2Sav/zetaSav;
752  syjSav = sAnt*zetaSav;
753  }
754  }
755  }
756
757  // Generate scale for W in final state.
758  if (abs(idy) == 24) {
759  double zMin = q2Low/sjkMax;
760  double zMax = sajMax/sAnt;
761  if (zMin < zMax) {
762  double Iz = log((1-zMin)/(1-zMax));
763  double comFac = 3.*M_PI*sqrt(lambda)/alpha/Iz/c/(sAnt+sjkMax)/2.;
764  double q2New = q2Start;
765  double zetaNew, sxjNew, syjNew;
766  while (true) {
767  q2New *= pow(rndmPtr->flat(), comFac);
768  if (q2New < q2Sav) {break;}
769  zetaNew = 1. - (1-zMin)*pow((1-zMax)/(1-zMin),rndmPtr->flat());
770  sxjNew = sAnt*zetaNew + q2New;
771  syjNew = q2New/zetaNew;
772
773  // Veto probability.
774  double pVeto = (syjNew+sAnt)/(sjkMax+sAnt);
775  if (rndmPtr->flat() < pVeto) {
776  q2Sav = q2New;
777  zetaSav = zetaNew;
778  sxjSav = sxjNew;
779  syjSav = syjNew;
780  break;
781  }
782  }
783  }
784  }
785  }
786  phiSav = 2.*M_PI*rndmPtr->flat();
787  hasTrial = true;
788  return q2Sav;
789
790 }
791
792 //==========================================================================
793
794 // Class for a QED emission system.
795
796 //--------------------------------------------------------------------------
797
798 // Initialize the pointers.
799
800 void QEDemitSystem::initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) {
801  infoPtr = infoPtrIn;
802  particleDataPtr = infoPtr->particleDataPtr;
803  partonSystemsPtr = infoPtr->partonSystemsPtr;
804  rndmPtr = infoPtr->rndmPtr;
805  settingsPtr = infoPtr->settingsPtr;
806  vinComPtr = vinComPtrIn;
807  isInitPtr = true;
808 }
809
810 //--------------------------------------------------------------------------
811
812 // Initialize settings for current run.
813
814 void QEDemitSystem::init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
815  int verboseIn) {
816
817  // Verbose setting.
818  if (!isInitPtr)
819  printOut(__METHOD_NAME__,"QEDemitSystem:initPtr not called");
820  verbose = verboseIn;
821
822  // Set beam pointers.
823  beamAPtr = beamAPtrIn;
824  beamBPtr = beamBPtrIn;
825
826  // Settings.
827  mode = settingsPtr->mode("Vincia:photonEmissionMode");
828  useFullWkernel = settingsPtr->flag("Vincia:fullWkernel");
830
831  // Constants.
832  TINYPDF = pow(10,-10);
833
834  // Initialized.
835  isInit = true;
836
837 }
838
839 //--------------------------------------------------------------------------
840
841 // Prepare a QED system.
842
843 void QEDemitSystem::prepare(int iSysIn, Event &event, double q2CutIn,
844  bool isBelowHadIn, vector<double> evolutionWindowsIn, AlphaEM alIn) {
845
846  if (!isInit) {
847  infoPtr->errorMsg("Error in "+__METHOD_NAME__+": Not initialised.");
848  return;
849  }
850
851  // Verbose output.
852  if (verbose >= louddebug) printOut(__METHOD_NAME__, "begin --------------");
853
854  // Input.
855  iSys = iSysIn;
856  shh = infoPtr->s();
857  q2Cut = q2CutIn;
859  evolutionWindows = evolutionWindowsIn;
860  al = alIn;
861
862  // Build internal system.
863  buildSystem(event);
864  if (verbose >= louddebug) printOut(__METHOD_NAME__, "end --------------");
865
866 }
867
868 //--------------------------------------------------------------------------
869
870 // Trial antenna function.
871
872 double QEDemitSystem::aTrial(QEDemitElemental* ele, double sxj, double syj,
873  double sxy) {
874  int idx = ele->idx;
875  int idy = ele->idy;
876  double ant = 0;
877
878  // FF.
879  if (ele->isFF || ele->isDip) {
880  double s = sxj + syj + sxy;
881  ant += 4*s/sxj/syj;
882  if (ele->isFF && abs(idx) == 24) ant += 8.*s/sxj/(s - syj)/3.;
883  if (ele->isFF && abs(idy) == 24) ant += 8.*s/syj/(s - sxj)/3.;
884  }
885
886  // IF.
887  if (ele->isIF) {
888  double s = sxj + sxy - syj;
889  ant += 4*pow2(s+syj)/(s*sxj*syj);
890  if (abs(idy) == 24) ant += 8.*(s + syj)/syj/(s + syj - sxj)/3.;
891  }
892
893  // II.
894  if (ele->isII) {
895  double s = sxy - sxj - syj;
896  ant += 4*sxy*sxy/s/sxj/syj;
897  }
898
899  // RF.
900  if (ele->isRF) {
901  double s = sxj + sxy - syj;
902  ant += 4*pow2(s+syj)/s/sxj/syj;
903  if (abs(idx) == 24) ant += 8*(2.*syj/s + pow2(syj)/pow2(s))/sxj/3.;
904  if (abs(idy) == 24) ant += 8.*(s + syj)/syj/(s + syj - sxj)/3.;
905  }
906  return ant;
907
908 }
909
910 //--------------------------------------------------------------------------
911
912 // Physical antenna function.
913
914 double QEDemitSystem::aPhys(QEDemitElemental* ele, double sxj, double syj,
915  double sxy) {
916  double mx2 = ele->mx2;
917  double my2 = ele->my2;
918  int idx = ele->idx;
919  int idy = ele->idy;
920  double ant = 0;
921
922  // FF.
923  if (ele->isFF) {
924  double s = sxj + syj + sxy;
925  // Eikonal.
926  ant += 4.*sxy/sxj/syj - 4.*mx2/sxj/sxj - 4.*my2/syj/syj;
927
928  // Check if x is a W or a fermion.
929  if (abs(idx) == 24 && useFullWkernel)
930  ant += (4./3.)*(syj/(s - syj) + syj*(s - syj)/s/s)/sxj;
931  else
932  ant += 2.*syj/sxj/s;
933
934  // Check if y is a W or a fermion.
935  if (abs(idy) == 24 && useFullWkernel)
936  ant += (4./3.)*(sxj/(s - sxj) + sxj*(s - sxj)/s/s)/syj;
937  else
938  ant += 2.*sxj/syj/s;
939  }
940
941  // FF (dipole).
942  if (ele->isDip) {
943  double s = sxj + syj + sxy;
944  ant += 4.*sxy/sxj/(sxj+syj) - 4.*mx2/sxj/sxj + 2.*syj/sxj/s;
945  }
946
947  // IF.
948  if (ele->isIF) {
949  double s = sxj + sxy - syj;
950  // Eikonal + initial state fermion.
951  // The initial state is never a W and has no mass.
952  ant += 4.*sxy/sxj/syj - 4.*my2/syj/syj + 2.*syj/sxj/s;
953
954  if (abs(idy) == 24 && useFullWkernel)
955  ant += (8./3.)*( sxj/(sxy + syj) + sxj/(s + syj)
956  - pow2(sxj)/pow2(s + syj) )/syj;
957  else
958  ant += 2.*sxj/s/syj;
959  }
960
961  // II.
962  if (ele->isII) {
963  double s = sxy - sxj - syj;
964  // Eikonal + fermion.
965  ant = 4.*sxy/sxj/syj + 2.*(sxj/syj + syj/sxj)/s;
966  }
967
968  // RF.
969  if (ele->isRF) {
970  double s = sxj + sxy - syj;
971  // Eikonal.
972  ant = 4.*sxy/sxj/syj - 4.*mx2/sxj/sxj - 4.*my2/syj/syj;
973
974  // Check if x is a W or a fermion
975  if (abs(idx) == 24 && useFullWkernel)
976  ant += (8./3.)*( syj/(s+syj) + syj/s + pow2(syj)/pow2(s) )/sxj;
977  else
978  ant += 2.*syj/sxj/s;
979
980  // Check if y is a W or a fermion.
981  if (abs(idy) == 24 && useFullWkernel)
982  ant += (8./3.)*( sxj/(sxy + syj) + sxj/(s + syj)
983  - pow2(sxj)/pow2(s + syj) )/syj;
984  else
985  ant += 2.*sxj/syj/s;
986  }
987  return ant;
988
989 }
990
991 //--------------------------------------------------------------------------
992
993 // Ratio between PDFs.
994
995 double QEDemitSystem::PDFratio(bool isA, double eOld, double eNew, int id,
996  double Qt2) {
997  double xOld = eOld/(sqrt(shh)/2.0);
998  double xNew = eNew/(sqrt(shh)/2.0);
999  double newPDF, oldPDF;
1000  if (isA) {
1001  newPDF = beamAPtr->xfISR(iSys, id, xNew, Qt2)/xNew;
1002  oldPDF = beamAPtr->xfISR(iSys, id, xOld, Qt2)/xOld;
1003  if (abs(newPDF) < TINYPDF) newPDF = TINYPDF;
1004  if (abs(oldPDF) < TINYPDF) oldPDF = TINYPDF;
1005  } else {
1006  newPDF = beamBPtr->xfISR(iSys, id, xNew, Qt2)/xNew;
1007  oldPDF = beamBPtr->xfISR(iSys, id, xOld, Qt2)/xOld;
1008  if (abs(newPDF) < TINYPDF) newPDF = TINYPDF;
1009  if (abs(oldPDF) < TINYPDF) oldPDF = TINYPDF;
1010  }
1011  return newPDF/oldPDF;
1012 }
1013
1014 //--------------------------------------------------------------------------
1015
1016 // Set up antenna pairing for incoherent mode.
1017
1018 void QEDemitSystem::buildSystem(Event &event) {
1019
1020  // Clear previous antennae.
1021  eleVec.clear();
1022  eleMat.clear();
1023  iCoh.clear();
1024
1025  // Construct hungarian algorithm solver.
1026  HungarianAlgorithm ha;
1029  map<int, vector<int> > posMap, negMap;
1030  vector<Vec4> posMoms, negMoms;
1031
1032  // Find all (final-state) quarks and leptons.
1033  vector<int> iQuarks, iLeptons;
1034  int sysSize = partonSystemsPtr->sizeOut(iSys);
1035  for (int i = 0; i < sysSize; i++) {
1036  int iEv = partonSystemsPtr->getOut(iSys, i);
1037  if (event[iEv].col() != 0 && event[iEv].acol()==0 &&
1038  event[iEv].isFinal()) {
1039  // For now, ignore quarks that are connected to junctions. In
1040  // principle, we could add them, and any antijunction dittos.
1041  bool isJun = false;
1042  for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
1043  for (int iLeg = 0; iLeg < 3; ++iLeg) {
1044  if (event[iEv].col() == event.endColJunction(iJun,iLeg)) {
1045  isJun = true;
1046  break;
1047  }
1048  }
1049  }
1050  if (!isJun) iQuarks.push_back(iEv);
1051  }
1052  if (event[iEv].isLepton() && event[iEv].isCharged())
1053  iLeptons.push_back(iEv);
1054  }
1055
1056  // Currently no showering below hadronisation scale if no leptons.
1057  if (iLeptons.size() == 0) return;
1058
1059  // Sort all leptons into maps.
1060  for (int i = 0; i < (int)iLeptons.size(); i++) {
1061  int iEv = iLeptons[i];
1062  vector<int> iLeptonVec;
1063  iLeptonVec.push_back(iEv);
1064  if (event[iEv].chargeType() == 3) {
1065  posMoms.push_back(event[iEv].p());
1066  posMap[posMoms.size()-1] = iLeptonVec;
1067  }
1068  if (event[iEv].chargeType() == -3) {
1069  negMoms.push_back(event[iEv].p());
1070  negMap[negMoms.size()-1] = iLeptonVec;
1071  }
1072  }
1073  // Find all colour strings.
1074  for (int i = 0; i < (int)iQuarks.size(); i++) {
1075  // Get initial quark and add to pseudo particle.
1076  Vec4 pPseudo;
1077  int iEv = iQuarks[i];
1078  vector<int> iPseudoVec;
1079  iPseudoVec.push_back(iEv);
1080  pPseudo += event[iEv].p();
1081
1082  // Find next colour-connected particle.
1083  do {
1084  int colTag = event[iEv].col();
1085  for (int j = 0; j < sysSize; j++) {
1086  int jEv = partonSystemsPtr->getOut(iSys, j);
1087  if (event[jEv].acol() == colTag && event[jEv].isFinal()) {
1088  iEv = jEv;
1089  break;
1090  }
1091  }
1092  if (iEv == iPseudoVec.back()) {
1093  infoPtr->errorMsg("Error in "+__METHOD_NAME__
1094  +": Colour tracing failed.");
1095  break;
1096  }
1097  iPseudoVec.push_back(iEv);
1098  pPseudo += event[iEv].p();
1099  }
1100  while(!event[iEv].isQuark()&&!event[iEv].isDiquark());
1101
1102  // Get charge of pseudoparticle and sort into maps.
1103  int chargeTypePseudo = event[iPseudoVec.front()].chargeType()
1104  + event[iPseudoVec.back()].chargeType();
1105  // Strings with only quarks are total charge 1 or -1.
1106  if (chargeTypePseudo == 3) {
1107  posMoms.push_back(pPseudo);
1108  posMap[posMoms.size()-1] = iPseudoVec;
1109  }
1110  if (chargeTypePseudo == -3) {
1111  negMoms.push_back(pPseudo);
1112  negMap[negMoms.size()-1] = iPseudoVec;
1113  }
1114  // Strings with a diquark can be charge 2 or -2. Add these
1115  // twice to list of recoilers.
1116  if (chargeTypePseudo == 6) {
1117  posMoms.push_back(pPseudo);
1118  posMap[posMoms.size()-1] = iPseudoVec;
1119  posMoms.push_back(pPseudo);
1120  posMap[posMoms.size()-1] = iPseudoVec;
1121  }
1122  if (chargeTypePseudo == -6) {
1123  negMoms.push_back(pPseudo);
1124  negMap[negMoms.size()-1] = iPseudoVec;
1125  negMoms.push_back(pPseudo);
1126  negMap[negMoms.size()-1] = iPseudoVec;
1127  }
1128  }
1129
1130  // If no leptons and overall hadronic system has charge = 0, do nothing.
1131  if (posMoms.size() == 0) return;
1132
1133  // Solve assignment problem.
1134  vector<vector<double> > weights;
1135  weights.resize(posMoms.size());
1136  for (int i=0; i<(int)posMoms.size(); i++) {
1137  weights[i].resize(negMoms.size());
1138  for (int j=0; j<(int)negMoms.size(); j++) {
1139  double w = posMoms[i]*negMoms[j]
1140  - posMoms[i].mCalc()*negMoms[j].mCalc();
1141  weights[i][j] = w;
1142  }
1143  }
1144  vector<int> assignment;
1145  ha.solve(weights, assignment);
1146
1147  for (int i = 0; i < (int)posMoms.size(); i++) {
1148  int iPos = i;
1149  int iNeg = assignment[i];
1150  // Only keep antennae with at least one lepton.
1151  if (posMap[iPos].size() == 1 || negMap[iNeg].size() == 1) {
1152  eleVec.push_back(QEDemitElemental());
1153  eleVec.back().initPtr(rndmPtr, partonSystemsPtr);
1154  // If two leptons, add regular antenna.
1155  if (posMap[iPos].size() == 1 && negMap[iNeg].size() == 1)
1156  eleVec.back().init(event, posMap[iPos][0], negMap[iNeg][0], shh,
1157  verbose);
1158  // If lepton + pseudoparticle, add dipole.
1159  if (posMap[iPos].size() == 1 && negMap[iNeg].size() != 1)
1160  eleVec.back().init(event, posMap[iPos][0], negMap[iNeg], shh,
1161  verbose);
1162  if (posMap[iPos].size()!=1 && negMap[iNeg].size()==1)
1163  eleVec.back().init(event, negMap[iNeg][0], posMap[iPos], shh,
1164  verbose);
1165  }
1166  }
1167
1170  // Collect relevant particles.
1171  int sysSize = partonSystemsPtr->sizeAll(iSys);
1172  for (int i = 0; i < sysSize; i++) {
1173  int iEv = partonSystemsPtr->getAll(iSys, i);
1174  if (event[iEv].isCharged()) iCoh.push_back(iEv);
1175  }
1176
1177  // Catch cases (like hadron->partons decays) where an explicit
1178  // charged mother may not have been added to the partonSystem as a
1179  // resonance.
1180  if (partonSystemsPtr->getInA(iSys) == 0 &&
1181  partonSystemsPtr->getInB(iSys) == 0 &&
1182  partonSystemsPtr->getInRes(iSys) == 0) {
1183  // Guess that the decaying particle is mother of first parton.
1184  int iRes = event[partonSystemsPtr->getOut(iSys, 0)].mother1();
1185  if (iRes != 0 && event[iRes].isCharged()) {
1186  // Check daughter list consistent with whole system.
1187  int ida1 = event[iRes].daughter1();
1188  int ida2 = event[iRes].daughter2();
1189  if (ida2 > ida1) {
1190  bool isOK = true;
1191  for (int i=0; i<partonSystemsPtr->sizeOut(iSys); ++i)
1192  if (partonSystemsPtr->getOut(iSys,i) < ida1
1193  || partonSystemsPtr->getOut(iSys,i) > ida2) isOK = false;
1194  if (isOK) {iCoh.push_back(iRes);}
1195  }
1196  }
1197  }
1198
1199  // First check charge conservation.
1200  int chargeTypeTot = 0;
1201  for (int i = 0; i < (int)iCoh.size(); i++) {
1202  double cType = event[iCoh[i]].chargeType();
1203  chargeTypeTot += (event[iCoh[i]].isFinal() ? cType : -cType);
1204  }
1205
1206  if (chargeTypeTot != 0) {
1207  infoPtr->errorMsg("Error in "+__METHOD_NAME__
1208  +": Charge not conserved above hadronization scale");
1209  if (verbose >= superdebug) {
1210  printOut(__METHOD_NAME__, "Printing events and systems");
1211  event.list();
1212  partonSystemsPtr->list();
1213  }
1214  }
1215
1216  // Pairing algorithm.
1217  if (mode == 1) {
1218  vector<vector<int> > posChargeTypes;
1219  posChargeTypes.resize(3);
1220  vector<vector<int> > negChargeTypes;
1221  negChargeTypes.resize(3);
1222
1223  for (int i = 0; i < (int)iCoh.size(); i++) {
1224  int iEv = iCoh[i];
1225  // Separate particles into charge types.
1226  double Q = event[iEv].charge();
1227  // Get index in pos/negChargeTypes.
1228  int n = abs(event[iEv].chargeType()) - 1;
1229  // Flip charge contribution of initial state.
1230  if (!event[iEv].isFinal()) {Q = -Q;}
1231  if (Q > 0) posChargeTypes[n].push_back(iEv);
1232  else negChargeTypes[n].push_back(iEv);
1233  }
1234
1235  // Clear list of charged particles.
1236  iCoh.clear();
1237
1238  // Solve assignment problems.
1239  for (int i=0; i<3; i++) {
1240  int posSize = posChargeTypes[i].size();
1241  int negSize = negChargeTypes[i].size();
1242  int maxSize = max(posSize,negSize);
1243  if (maxSize > 0) {
1244  vector<vector<double> > weights;
1245  weights.resize(maxSize);
1246  // Set up matrix of weights.
1247  for (int x = 0; x < maxSize; x++) {
1248  weights[x].resize(maxSize);
1249  for (int y = 0; y < maxSize; y++) {
1250  // If either index is out of range. Add some random
1251  // large weight.
1252  double wIn = (0.9 + 0.2*rndmPtr->flat())*1E300;
1253  if (x < posSize && y < negSize) {
1254  int xEv = posChargeTypes[i][x];
1255  int yEv = negChargeTypes[i][y];
1256  wIn = event[xEv].p()*event[yEv].p()
1257  - event[xEv].m()*event[yEv].m();
1258  }
1259  weights[x][y] = wIn;
1260  }
1261  }
1262
1263  // Find solution.
1264  vector<int> assignment;
1265  ha.solve(weights, assignment);
1266
1267  // Add pairings to list of emitElementals.
1268  // Add unpaired particles to index list for coherent algorithm.
1269  for (int j = 0; j < maxSize; j++) {
1270  int x = j;
1271  int y = assignment[j];
1272  if (x < posSize && y < negSize) {
1273  int xEv = posChargeTypes[i][x];
1274  int yEv = negChargeTypes[i][y];
1275  eleVec.push_back(QEDemitElemental());
1276  eleVec.back().initPtr(rndmPtr, partonSystemsPtr);
1277  eleVec.back().init(event, xEv, yEv, shh, verbose);
1278  } else if (x < posSize) {
1279  int xEv = posChargeTypes[i][x];
1280  iCoh.push_back(xEv);
1281  } else if (y < negSize) {
1282  int yEv = negChargeTypes[i][y];
1283  iCoh.push_back(yEv);
1284  }
1285  }
1286  }
1287  }
1288  }
1289
1290  // Create eleMat.
1291  eleMat.resize(iCoh.size());
1292  for (int i = 0; i < (int)iCoh.size(); i++) {
1293  eleMat[i].resize(i);
1294  for (int j = 0; j < i; j++) {
1295  eleMat[i][j].initPtr(rndmPtr, partonSystemsPtr);
1296  eleMat[i][j].init(event, iCoh[i], iCoh[j], shh, verbose);
1297  }
1298  }
1299
1300  // Compute overestimate constant.
1301  cMat = 0;
1302  for (int i = 0; i < (int)eleMat.size(); i++)
1303  for (int j = 0; j < i; j++) cMat += max(eleMat[i][j].QQ, 0.);
1304  }
1305
1306 }
1307
1308 //--------------------------------------------------------------------------
1309
1310 // Generate a trial scale.
1311
1312 double QEDemitSystem::generateTrialScale(Event &event, double q2Start) {
1313
1314  // Check if qTrial is below the cutoff.
1315  if (q2Start < q2Cut || evolutionWindows.size() == 0) {return 0;}
1316
1317  // Find lower value from evolution window.
1318  int iEvol = evolutionWindows.size() - 1;
1319  while (iEvol >= 1 && q2Start <= evolutionWindows[iEvol]) iEvol--;
1320  double q2Low = evolutionWindows[iEvol];
1321  if (q2Low < 0)
1322  infoPtr->errorMsg("Error in "+__METHOD_NAME__+": Evolution window < 0");
1323  double q2Trial = 0;
1324
1325  // Generate a scale.
1326  double alphaMax = al.alphaEM(q2Start);
1327
1328  // Pull scales from eleVec.
1329  for (int i = 0; i < (int)eleVec.size(); i++) {
1330  double c = eleVec[i].QQ;
1331  double q2New = eleVec[i].generateTrial(event, q2Start, q2Low, alphaMax, c);
1332  if (q2New > q2Low && q2New > q2Trial) {
1333  q2Trial = q2New;
1334  eleTrial = &eleVec[i];
1335  trialIsVec = true;
1336  }
1337  }
1338
1339  // Pull scales from eleMat.
1340  for (int i = 0; i < (int)eleMat.size(); i++) {
1341  for (int j = 0; j < i; j++) {
1342  double q2New = eleMat[i][j].generateTrial(event, q2Start, q2Low,
1343  alphaMax, cMat);
1344  if (q2New > q2Low && q2New > q2Trial) {
1345  q2Trial = q2New;
1346  eleTrial = &eleMat[i][j];
1347  trialIsVec = false;
1348  }
1349  }
1350  }
1351
1352  // Check if evolution window was crossed.
1353  if (q2Trial < q2Low) {
1354  if (iEvol == 0) {return 0;}
1355  // Reset all trials.
1356  for (int i = 0; i < (int)eleVec.size(); i++) eleVec[i].hasTrial = false;
1357  for (int i=0; i<(int)eleMat.size(); i++)
1358  for (int j=0; j<i; j++) eleMat[i][j].hasTrial = false;
1359  return generateTrialScale(event, q2Low);
1360  }
1361
1362  // Otherwise return trial scale.
1363  return q2Trial;
1364
1365 }
1366
1367 //--------------------------------------------------------------------------
1368
1369 // Check the veto.
1370
1371 bool QEDemitSystem::checkVeto(Event &event) {
1372
1373  // Mark trial as used.
1374  if (verbose >= debug) printOut(__METHOD_NAME__, "begin --------------");
1375  eleTrial->hasTrial = false;
1376
1377  // Pre- and post-branching momenta.
1378  vector<Vec4> pOld, pNew;
1379
1380  // Global recoil momenta.
1381  vector<Vec4> pRec;
1382  vector<int> iRec;
1383
1384  // II.
1385  if (eleTrial->isII) {
1386  double saj = eleTrial->sxjSav;
1387  double sbj = eleTrial->syjSav;
1388  double phi = eleTrial->phiSav;
1389  double sAB = eleTrial->sAnt;
1390  double sab = sAB + saj + sbj;
1391
1392  // Pre-branching momenta.
1393  pOld.push_back(event[eleTrial->x].p());
1394  pOld.push_back(event[eleTrial->y].p());
1395
1396  // Collect the recoiling final state particles.
1397  int sysSize = partonSystemsPtr->sizeAll(iSys);
1398  for (int i = 0; i < sysSize; i++) {
1399  int iEv = partonSystemsPtr->getAll(iSys, i);
1400  if (iEv < 0 || !event[iEv].isFinal()) continue;
1401  if (iEv == eleTrial->x || iEv == eleTrial->y) continue;
1402  pRec.push_back(event[iEv].p());
1403  iRec.push_back(iEv);
1404  }
1405
1406  // Kinematics.
1407  if (!vinComPtr->map2to3IImassless(pNew, pRec, pOld, sAB,saj,sbj,sab, phi))
1408  return false;
1409
1410  // Check if new energies don't exceed hadronic maxima.
1411  double eaUsed = 0, ebUsed = 0;
1412  int nSys = partonSystemsPtr->sizeSys();
1413  for (int i = 0; i < nSys; i++) {
1414  eaUsed += event[partonSystemsPtr->getInA(i)].e();
1415  ebUsed += event[partonSystemsPtr->getInB(i)].e();
1416  }
1417  if ((eaUsed - pOld[0].e() + pNew[0].e()) > 0.98*sqrt(shh)/2.) return false;
1418  if ((ebUsed - pOld[1].e() + pNew[2].e()) > 0.98*sqrt(shh)/2.) return false;
1419  }
1420
1421  // IF.
1422  else if (eleTrial->isIF) {
1423  double saj = eleTrial->sxjSav;
1424  double sjk = eleTrial->syjSav;
1425  double phi = eleTrial->phiSav;
1426  double sAK = eleTrial->sAnt;
1427  double sak = sAK + sjk - saj;
1428  double mK2 = eleTrial->my2;
1429
1430  // Check phase space.
1431  if (sak < 0 || saj*sjk*sak - saj*saj*mK2 < 0) {return false;}
1432
1433  // Pre-branching momenta.
1434  pOld.push_back(event[eleTrial->x].p());
1435  pOld.push_back(event[eleTrial->y].p());
1436
1437  // Kinematics. (TODO: check if need for global kinematics map here).
1438  if (!vinComPtr->map2to3IFlocal(pNew, pOld, sAK, saj, sjk, sak, phi,
1439  mK2, 0, mK2)) return false;
1440
1441  // Check if new energy doesn't exceed the hadronic maximum.
1442  double eaUsed = 0;
1443  int nSys = partonSystemsPtr->sizeSys();
1444  for (int i = 0; i < nSys; i++) {
1445  int iEv;
1446  if (eleTrial->isIA) iEv = partonSystemsPtr->getInA(i);
1447  else iEv = partonSystemsPtr->getInB(i);
1448  eaUsed += event[iEv].e();
1449  }
1450  if ((eaUsed - pOld[0].e() + pNew[0].e()) > 0.98*sqrt(shh)/2.) return false;
1451  }
1452
1453  // RF.
1454  else if (eleTrial->isRF) {
1455  double saj = eleTrial->sxjSav;
1456  double sjk = eleTrial->syjSav;
1457  double sAK = eleTrial->sAnt;
1458  double sak = sAK + sjk - saj;
1459  double phi = eleTrial->phiSav;
1460  double mA2 = eleTrial->mx2;
1461  double mK2 = eleTrial->my2;
1462
1463  // Check phase space.
1464  if (sak < 0 || saj*sjk*sak - saj*saj*mK2 - sjk*sjk*mA2 < 0) return false;
1465
1466  // Pre-branching momenta.
1467  pOld.push_back(event[eleTrial->x].p());
1468  pOld.push_back(event[eleTrial->y].p());
1469
1470  // Collect the recoiling final state particles.
1471  int sysSize = partonSystemsPtr->sizeAll(iSys);
1472  for (int i = 0; i < sysSize; i++) {
1473  int iEv = partonSystemsPtr->getAll(iSys, i);
1474  if (iEv < 0 || !event[iEv].isFinal()) {continue;}
1475  if (iEv == eleTrial->x || iEv == eleTrial->y) {continue;}
1476  pRec.push_back(event[iEv].p());
1477  iRec.push_back(iEv);
1478  }
1479
1480  // Do kinematics.
1481  vector<double> masses;
1482  masses.push_back(sqrt(mA2));
1483  masses.push_back(0.);
1484  masses.push_back(sqrt(mK2));
1485  masses.push_back(sqrtpos(mA2+mK2-sAK));
1486  vector<double> invariants;
1487  invariants.push_back(sAK);
1488  invariants.push_back(saj);
1489  invariants.push_back(sjk);
1490  invariants.push_back(sak);
1491  vector<Vec4> pAfter;
1492  vector<Vec4> pBefore = pOld;
1493  pBefore.insert(pBefore.end(), pRec.begin(), pRec.end());
1494  if (!vinComPtr->map2toNRFmassive(pAfter, pBefore, 0, 1, invariants, phi,
1495  masses)) return false;
1496  pNew.push_back(pAfter[0]);
1497  pNew.push_back(pAfter[1]);
1498  pNew.push_back(pAfter[2]);
1499
1500  // Replace momenta with boosted counterpart.
1501  pRec.clear();
1502  for (int i = 3; i < (int)pAfter.size(); i++) pRec.push_back(pAfter[i]);
1503
1504  // Check if nothing got messed up along the way.
1505  if (pRec.size() != iRec.size()) {
1506  infoPtr->errorMsg("Error in "+__METHOD_NAME__
1507  +": inconsistent recoilers in RF kinematics.");
1508  return false;
1509  }
1510  }
1511
1512  // FF>
1513  else if (eleTrial->isFF) {
1514  double sIK = eleTrial->m2Ant - eleTrial->mx2 - eleTrial->my2;
1515  double sij = eleTrial->sxjSav;
1516  double sjk = eleTrial->syjSav;
1517  double sik = sIK - sij - sjk;
1518  double mi = sqrt(eleTrial->mx2);
1519  double mk = sqrt(eleTrial->my2);
1520  double phi = eleTrial->phiSav;
1521
1522  vector<double> invariants;
1523  invariants.push_back(sIK);
1524  invariants.push_back(sij);
1525  invariants.push_back(sjk);
1526
1527  vector<double> masses;
1528  masses.push_back(mi);
1529  masses.push_back(0);
1530  masses.push_back(mk);
1531
1532  // Check phase space.
1533  if (sik < 0) return false;
1534  if (sij*sjk*sik - pow2(sij)*pow2(mk) - pow2(sjk)*pow2(mi) < 0)
1535  return false;
1536
1537  // Pre-branching momenta.
1538  pOld.push_back(event[eleTrial->x].p());
1539  pOld.push_back(event[eleTrial->y].p());
1540
1541  // Kinematics.
1542  if (!vinComPtr->map2to3FF(pNew, pOld, 3, invariants, phi, masses))
1543  return false;
1544  }
1545
1546  // Dipole.
1547  else if (eleTrial->isDip) {
1548  // Construct recoiler momentum.
1549  Vec4 pk;
1550  for (int i = 0; i < (int)eleTrial->iRecoil.size(); i++)
1551  pk += event[eleTrial->iRecoil[i]].p();
1552  double sIK = eleTrial->m2Ant - eleTrial->mx2 - eleTrial->my2;
1553  double sij = eleTrial->sxjSav;
1554  double sjk = eleTrial->syjSav;
1555  double sik = sIK - sij - sjk;
1556  double mi = sqrt(eleTrial->mx2);
1557  double mk = pk.mCalc();
1558  double phi = eleTrial->phiSav;
1559
1560  vector<double> invariants;
1561  invariants.push_back(sIK);
1562  invariants.push_back(sij);
1563  invariants.push_back(sjk);
1564
1565  vector<double> masses;
1566  masses.push_back(mi);
1567  masses.push_back(0);
1568  masses.push_back(mk);
1569
1570  // Check phase space.
1571  if (sik < 0) {return false;}
1572  if (sij*sjk*sik - pow2(sij)*pow2(mk) - pow2(sjk)*pow2(mi) < 0)
1573  return false;
1574
1575  // Pre-branching momenta.
1576  pOld.push_back(event[eleTrial->x].p());
1577  pOld.push_back(pk);
1578
1579  // Kinematics.
1580  if (!vinComPtr->map2to3FF(pNew, pOld, 3, invariants, phi, masses))
1581  return false;
1582  }
1583  Vec4 pPhot = pNew[1];
1584  Vec4 px = pNew[0];
1585  Vec4 py = pNew[2];
1586  int x = eleTrial->x;
1587  int y = eleTrial->y;
1588  double sxj = eleTrial->sxjSav;
1589  double syj = eleTrial->syjSav;
1590  double sxy = px*py*2.;
1591
1592  // Compute veto probability.
1593  double pVeto = 1.;
1594
1596  pVeto *= al.alphaEM(eleTrial->q2Sav) / eleTrial->alpha;
1597  if (pVeto > 1) printOut(__METHOD_NAME__, "Alpha increased");
1598
1599  // Add antenna veto. Simple veto for eleTrial in eleVec.
1600  if (trialIsVec) {
1601  // Note that charge factor is included at generation step.
1602  double aTrialNow = aTrial(eleTrial, sxj, syj, sxy);
1603  double aPhysNow = aPhys(eleTrial, sxj, syj, sxy);
1604
1605  if (aPhysNow/aTrialNow > 1.001) {
1606  stringstream ss1;
1607  ss1 << "Incorrect overestimate (eleVec) " << aPhysNow/aTrialNow;
1608  printOut(__METHOD_NAME__, ss1.str());
1609  if (verbose > louddebug) {
1610  stringstream ss2, ss3;
1611  if (eleTrial->isFF) ss2 << "Antenna is FF";
1612  if (eleTrial->isIF) ss2 << "Antenna is IF";
1613  if (eleTrial->isRF) ss2 << "Antenna is RF";
1614  if (eleTrial->isII) ss2 << "Antenna is II";
1615  printOut(__METHOD_NAME__, ss2.str());
1616  printOut(__METHOD_NAME__, ss3.str());
1617  }
1618  }
1619  pVeto *= aPhysNow/aTrialNow;
1620
1621  // Construct full branching kernel for eleTrial in eleMat. Perform
1622  // sector check too.
1623  } else {
1624  double aSectorNow = aPhys(eleTrial, sxj, syj, sxy);
1625  double aTrialFull = eleTrial->c*aTrial(eleTrial, sxj, syj, sxy);
1626  double aPhysFull = 0;
1627
1628  // Build map of momenta & invariants with new photon.
1629  map<int, double> sj;
1630  map<int, Vec4> p;
1631  // Loop over the first column in eleMat.
1632  for (int i = 0; i < (int)iCoh.size(); i++) {
1633  int iEv = iCoh[i];
1634  // If the particle is in eleTrial, use shower variables.
1635  if (iEv == x) {
1636  p[iEv] = px;
1637  sj[iEv] = sxj;
1638  } else if (iEv == y) {
1639  p[iEv] = py;
1640  sj[iEv] = syj;
1641  // Otherwise get the momenta elsewhere
1642  } else {
1643  // If global recoil, get them from pRec.
1644  if (eleTrial->isII) {
1645  // Find index.
1646  for (int j = 0; j < (int)iRec.size(); j++) {
1647  if (iEv == iRec[j]) {
1648  p[iEv] = pRec[j];
1649  sj[iEv] = 2.*pRec[j]*pPhot;
1650  break;
1651  }
1652  }
1653  // Otherwise use momentum from event.
1654  } else {
1655  p[iEv] = event[iEv].p();
1656  sj[iEv] = 2.*event[iEv].p()*pPhot;
1657  }
1658  }
1659  }
1660
1661  // Then build aPhys.
1662  for (int v=0; v<(int)eleMat.size(); v++) {
1663  for (int w=0; w<v; w++) {
1664  double sxjNow = sj[eleMat[v][w].x];
1665  double syjNow = sj[eleMat[v][w].y];
1666  double sxyNow = 2.*p[eleMat[v][w].x]*p[eleMat[v][w].y];
1667  double aPhysNow = aPhys(&eleMat[v][w], sxjNow, syjNow, sxyNow);
1668
1669  // Sector veto.
1670  if (aPhysNow > aSectorNow) return false;
1671
1672  // Add aPhysNew to aPhys.
1673  aPhysFull += eleMat[v][w].QQ*aPhysNow;
1674  }
1675  }
1676  // Set aPhys to zeto if below zero.
1677  if (aPhysFull < 0) {aPhysFull = 0;}
1678
1679  // Check overestimate.
1680  if (aPhysFull/aTrialFull > 1) {
1681  stringstream ss1;
1682  ss1 << "Incorrect overestimate (eleVec) " << aPhysFull/aTrialFull;
1683  printOut(__METHOD_NAME__, ss1.str());
1684  }
1686  pVeto *= aPhysFull/aTrialFull;
1687  }
1688
1690  if (eleTrial->isIF) {
1691  pVeto *= PDFratio(eleTrial->isIA, pOld[0].e(), pNew[0].e(),
1692  eleTrial->idx, eleTrial->q2Sav);
1693  }
1694  if (eleTrial->isII) {
1695  pVeto *= PDFratio(true, pOld[0].e(), pNew[0].e(),
1696  eleTrial->idx, eleTrial->q2Sav);
1697  pVeto *= PDFratio(false, pOld[1].e(), pNew[2].e(),
1698  eleTrial->idy, eleTrial->q2Sav);
1699  }
1700
1701  // Perform veto.
1702  if (rndmPtr->flat() > pVeto) return false;
1703
1704  // Accepted, now fix the event. Different procedures for dipole and
1705  // antennae.
1706
1707  // If it is a dipole.
1708  if (eleTrial->isDip) {
1709  // Set up new particles.
1710  Particle newPhoton(22, 51, 0, 0, 0, 0, 0, 0, pPhot);
1711  Particle newPartx = event[x];
1712  newPartx.p(px);
1713
1715  int xNew = event.append(newPartx);
1716  int jNew = event.append(newPhoton);
1717
1718  // Update parton system.
1719  partonSystemsPtr->replace(iSys, x, xNew);
1721
1722  // Set old particles to negative.
1723  event[x].statusNeg();
1724
1725  // Update mother-daughter structure.
1726  event[xNew].mothers(x,0);
1727  event[jNew].mothers(x,0);
1728  event[x].daughters(xNew, jNew);
1729  event[xNew].daughters(0,0);
1730  event[jNew].daughters(0,0);
1731  event[xNew].statusCode(51);
1732  event[jNew].statusCode(51);
1733
1734  // Boost momenta and update.
1735  for (int i = 0; i < (int)eleTrial->iRecoil.size(); i++) {
1736  int iDipRec = eleTrial->iRecoil[i];
1737  Vec4 pDipRec = event[iDipRec].p();
1738  pDipRec.bstback(pOld[1]);
1739  pDipRec.bst(pNew[2]);
1740  // Copy the recoiler.
1741  int iDipRecNew = event.copy(iDipRec, 52);
1742  // Change the momentum.
1743  event[iDipRecNew].p(pDipRec);
1744  // Update parton system.
1745  partonSystemsPtr->replace(iSys, iDipRec, iDipRecNew);
1746  }
1747
1748  } else if (eleTrial->isRF) {
1749
1750  // Set up new particles.
1751  Particle newPhoton(22, 51, 0, 0, 0, 0, 0, 0, pPhot);
1752  Particle newParty = event[y];
1753  newParty.p(py);
1754
1755  // Add branched particles to event.
1756  int jNew, yNew;
1757  if (x < y) {
1758  jNew = event.append(newPhoton);
1759  yNew = event.append(newParty);
1760  } else {
1761  yNew = event.append(newParty);
1762  jNew = event.append(newPhoton);
1763  }
1764
1765  // Update parton system.
1767  partonSystemsPtr->replace(iSys, y, yNew);
1768
1769  // Set old particles to negative.
1770  event[y].statusNeg();
1771  event[jNew].mothers(y,0);
1772  event[yNew].mothers(y,0);
1773  event[jNew].daughters(0,0);
1774  event[yNew].daughters(0,0);
1775  event[y].daughters(yNew,jNew);
1776  event[yNew].statusCode(51);
1777
1778  // Update event for global recoil.
1779  vector<pair<int,int> > iRecNew;
1780  iRecNew.clear();
1781  iRecNew.resize(0);
1782  for (int j = 0; j < event.size(); j++)
1783  if (event[j].isFinal())
1784  for (int k = 0; k < (int)iRec.size(); k++)
1785  if (iRec[k] == j) {
1786  // Copy the recoiler.
1787  int inew = event.copy(j, 52);
1788  // Change the momentum.
1789  event[inew].p(pRec[k]);
1790  iRecNew.push_back(make_pair(iRec[k], inew));
1791  }
1792
1793  // Update parton system.
1794  for (int k=0; k<(int)iRecNew.size(); k++)
1795  partonSystemsPtr->replace(iSys, iRecNew[k].first, iRecNew[k].second);
1796
1797  // If it is an antenna
1798  } else {
1799  // Set up new particles.
1800  Particle newPhoton(22, 51, 0, 0, 0, 0, 0, 0, pPhot);
1801  Particle newPartx = event[x];
1802  newPartx.p(px);
1803  Particle newParty = event[y];
1804  newParty.p(py);
1805
1806  // Add branched particles to event.
1807  int xNew, jNew, yNew;
1808  if (x < y) {
1809  xNew = event.append(newPartx);
1810  jNew = event.append(newPhoton);
1811  yNew = event.append(newParty);
1812  } else {
1813  yNew = event.append(newParty);
1814  jNew = event.append(newPhoton);
1815  xNew = event.append(newPartx);
1816  }
1817
1818  // Update parton system.
1819  partonSystemsPtr->replace(iSys, x, xNew);
1821  partonSystemsPtr->replace(iSys, y, yNew);
1822
1823  // Set old particles to negative.
1824  event[x].statusNeg();
1825  event[y].statusNeg();
1826
1827  // Update everything.
1828  if (eleTrial->isII) {
1829  event[xNew].mothers(event[x].mother1(), event[x].mother2());
1830  event[jNew].mothers(xNew, yNew);
1831  event[yNew].mothers(event[y].mother1(), event[y].mother2());
1832  event[x].mothers(xNew, 0);
1833  event[y].mothers(yNew, 0);
1834  event[xNew].daughters(jNew, x);
1835  event[yNew].daughters(jNew, y);
1836  event[jNew].daughters(0,0);
1837  event[xNew].status(-41);
1838  event[yNew].status(-41);
1839  event[jNew].status(43);
1840
1841  // Update beam daughters.
1842  if (iSys == 0) {
1843  bool founda = false;
1844  bool foundb = false;
1845  for (int i = 0; i < (int)event.size(); i++) {
1846  if (!founda)
1847  if (event[i].daughter1() == x) {
1848  event[i].daughters(xNew, 0);
1849  founda = true;
1850  }
1851  if (!foundb)
1852  if (event[i].daughter1() == y) {
1853  event[i].daughters(yNew, 0);
1854  foundb = true;
1855  }
1856  if (founda && foundb) break;
1857  }
1858  }
1859
1860  // Update event for global recoil.
1861  vector<pair<int,int> > iRecNew;
1862  iRecNew.clear();
1863  iRecNew.resize(0);
1864  for (int j=0; j<event.size(); j++)
1865  if (event[j].isFinal())
1866  for (int k=0; k<(int)iRec.size(); k++)
1867  if (iRec[k] == j) {
1868  // Copy the recoiler.
1869  int inew = event.copy(j, 44);
1870  // Change the momentum.
1871  event[inew].p(pRec[k]);
1872  iRecNew.push_back(make_pair(iRec[k], inew));
1873  }
1874
1875
1876  // Update parton system.
1877  for (int k = 0; k < (int)iRecNew.size(); k++)
1878  partonSystemsPtr->replace(iSys, iRecNew[k].first, iRecNew[k].second);
1879  partonSystemsPtr->setInA(iSys, xNew);
1880  partonSystemsPtr->setInB(iSys, yNew);
1881
1882  // Update beams.
1883  BeamParticle& beam1 = *beamAPtr;
1884  BeamParticle& beam2 = *beamBPtr;
1885
1886  // Check that x is always a with pz>0.
1887  if (event[xNew].pz() < 0) {
1888  printOut(__METHOD_NAME__, "Swapped II antenna");
1889  beam1 = *beamBPtr;
1890  beam2 = *beamAPtr;
1891  }
1892  beam1[iSys].update(xNew, event[xNew].id(), event[xNew].e()/beam1.e());
1893  beam2[iSys].update(yNew, event[yNew].id(), event[yNew].e()/beam2.e());
1894  }
1895
1896  if (eleTrial->isIF) {
1897  event[xNew].mothers(event[x].mother1(), event[x].mother2());
1898  event[jNew].mothers(y,xNew);
1899  event[yNew].mothers(y,0);
1900  event[x].mothers(xNew,0);
1901  event[xNew].daughters(jNew,x);
1902  event[jNew].daughters(0,0);
1903  event[yNew].daughters(0,0);
1904  event[y].daughters(jNew, yNew);
1905  event[xNew].status(-41);
1906  event[yNew].status(43);
1907  event[jNew].status(43);
1908
1909  // Update beam daughter.
1910  if (iSys == 0)
1911  for (int i=0; i<(int)event.size(); i++)
1912  if (event[i].daughter1() == x) {
1913  event[i].daughters(xNew, 0);
1914  break;
1915  }
1916
1917  // Update parton system.
1918  if (eleTrial->isIA) partonSystemsPtr->setInA(iSys, xNew);
1919  else partonSystemsPtr->setInB(iSys, xNew);
1920
1921  // Update beams.
1922  BeamParticle& beam = (eleTrial->isIA ? *beamAPtr : *beamBPtr);
1923  beam[iSys].update(xNew, event[xNew].id(), event[xNew].e()/beam.e());
1924  }
1925
1926  if (eleTrial->isFF) {
1927  event[xNew].mothers(x,0);
1928  event[jNew].mothers(x,y);
1929  event[yNew].mothers(y,0);
1930  event[x].daughters(xNew, jNew);
1931  event[y].daughters(yNew, jNew);
1932  event[xNew].daughters(0,0);
1933  event[jNew].daughters(0,0);
1934  event[yNew].daughters(0,0);
1935  event[xNew].statusCode(51);
1936  event[jNew].statusCode(51);
1937  event[yNew].statusCode(51);
1938  }
1939
1940  // Update event pointers.
1941  event.restorePtrs();
1942
1943  // Fix sHat for parton system.
1944  double shat = (event[partonSystemsPtr->getInA(iSys)].p() +
1945  event[partonSystemsPtr->getInB(iSys)].p()).m2Calc();
1946  partonSystemsPtr->setSHat(iSys, shat);
1947  }
1948  if (verbose >= debug) {
1949  if (verbose >= superdebug) {
1950  event.list();
1951  partonSystemsPtr->list();
1952  }
1953  printOut(__METHOD_NAME__, "end --------------");
1954  }
1955  return true;
1956
1957 }
1958
1959 //--------------------------------------------------------------------------
1960
1961 // Print the QED emit internal system.
1962
1963 void QEDemitSystem::print() {
1964  cout << "Printing QEDemit internal system" << endl;
1965  cout << "Pairing elementals" << endl;
1966  for (int i = 0; i < (int)eleVec.size(); i++) {
1967  if (eleVec[i].isDip) {
1968  cout << "Dipole: x = ";
1969  cout << eleVec[i].x << " Recoilers: (";
1970  for (int j = 0; j < (int)eleVec[i].iRecoil.size(); j++) {
1971  cout << eleVec[i].iRecoil[j] << ", ";
1972  if (j == (int)eleVec[i].iRecoil.size()-1) cout << ")";
1973  else cout << ", ";
1974  }
1975  } else cout << "Antennae: x = " << eleVec[i].x << ", y = " << eleVec[i].y;
1976  cout << ", QQ = " << eleVec[i].QQ << ", s = " << eleVec[i].sAnt << endl;
1977  }
1978  cout << "Coherent elementals" << endl;
1979  for (int i = 0; i < (int)eleMat.size(); i++)
1980  for (int j = 0; j < i; j++)
1981  cout << "x = " << eleMat[i][j].x << ", y = " << eleMat[i][j].y
1982  << ", QQ = " << eleMat[i][j].QQ << ", s = " << eleMat[i][j].sAnt
1983  << endl;
1984 }
1985
1986 //==========================================================================
1987
1988 // Class for a QED splitting system.
1989
1990 //--------------------------------------------------------------------------
1991
1992 // Initialize pointers.
1993
1994 void QEDsplitSystem::initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) {
1995  infoPtr = infoPtrIn;
1996  particleDataPtr = infoPtr->particleDataPtr;
1997  partonSystemsPtr = infoPtr->partonSystemsPtr;
1998  rndmPtr = infoPtr->rndmPtr;
1999  settingsPtr = infoPtr->settingsPtr;
2000  vinComPtr = vinComPtrIn;
2001  isInitPtr = true;
2002 }
2003
2004 //--------------------------------------------------------------------------
2005
2006 // Initialize.
2007
2008 void QEDsplitSystem::init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
2009  int verboseIn) {
2010  if (!isInitPtr) printOut(__METHOD_NAME__, "initPtr not called");
2011  verbose = verboseIn;
2012  q2Max = pow2(settingsPtr->parm("Vincia:mMaxGamma"));
2013  nLepton = settingsPtr->mode("Vincia:nGammaToLepton");
2014  nQuark = settingsPtr->mode("Vincia:nGammaToQuark");
2015  beamAPtr = beamAPtrIn;
2016  beamBPtr = beamBPtrIn;
2017  isInit = true;
2018 }
2019
2020 //--------------------------------------------------------------------------
2021
2022 // Prepare list of final-state photons - with recoilers - for splittings.
2023
2024 void QEDsplitSystem::prepare(int iSysIn, Event &event, double q2CutIn,
2025  bool isBelowHadIn, vector<double> evolutionWindowsIn, AlphaEM alIn) {
2026
2027  if (!isInit) {
2028  infoPtr->errorMsg("Error in "+__METHOD_NAME__+": Not initialised.");
2029  return;
2030  }
2031  if (verbose >= louddebug) printOut(__METHOD_NAME__, "begin --------------");
2032
2033  // Input.
2034  iSys = iSysIn;
2035  q2Cut = q2CutIn;
2037  evolutionWindows = evolutionWindowsIn;
2038  al = alIn;
2039
2040  // Set up weights for splitting flavours.
2041  ids.clear();
2042  idWeights.clear();
2043  totIdWeight = 0;
2044  maxIdWeight = 0;
2045
2046  // Splittings for gamma->lepton+lepton-.
2047  for (int i = 0; i < nLepton; i++) {
2048  ids.push_back(11 + 2*i);
2049  idWeights.push_back(1);
2050  }
2051  // Only include gamma->qqbar if above hadronisation scale.
2053  for (int i = 1; i <= nQuark; i++) {
2054  ids.push_back(i);
2055  idWeights.push_back((i%2==0 ? 4./3. : 1./3.));
2056  }
2057  }
2058  // Total weights.
2059  for (int i=0; i<(int)ids.size(); i++) {
2060  totIdWeight += idWeights[i];
2061  if (idWeights[i] > maxIdWeight) maxIdWeight = idWeights[i];
2062  }
2063
2064  // Build internal system.
2065  buildSystem(event);
2066  if (verbose >= louddebug) printOut(__METHOD_NAME__, "end --------------");
2067
2068 }
2069
2070 //--------------------------------------------------------------------------
2071
2072 // Build the splitting system.
2073
2074 void QEDsplitSystem::buildSystem(Event &event) {
2075
2076  // Get rid of saved trial and clear all antennae.
2077  hasTrial = false;
2078  eleVec.clear();
2079
2080  // Build lists of particles.
2081  vector<int> photList, chSpecList, uchSpecList;
2082  int sysSize = partonSystemsPtr->sizeAll(iSys);
2083  for (int i = 0; i < sysSize; i++) {
2084  int iEv = partonSystemsPtr->getAll(iSys, i);
2085  if (iEv > 0) {
2086  // Only involve final state particles.
2087  if (event[iEv].isFinal()) {
2088  // Find photons.
2089  if (event[iEv].id()==22) photList.push_back(iEv);
2090  // Find recoilers.
2091  if (event[iEv].isCharged()) chSpecList.push_back(iEv);
2092  else uchSpecList.push_back(iEv);
2093  }
2094  }
2095  }
2096
2097  // If no charged and no uncharged spectators, return.
2098  if (chSpecList.empty() && uchSpecList.empty()) return;
2099
2100  // Loop over photons.
2101  for (int i = 0; i < (int)photList.size(); i++) {
2102  int iPhot = photList[i];
2103  // If no charged spectators, use uncharged.
2104  if (chSpecList.empty()) {
2105  // Check if there is another spectator than the current photon.
2106  bool otherSpecAvail = false;
2107  for (int j = 0; j < (int)uchSpecList.size(); j++)
2108  if (uchSpecList[j] != iPhot) {otherSpecAvail = true; break;}
2109  // Continue to next photon if no spectator is available.
2110  if (!otherSpecAvail) continue;
2111
2112  // Select one at random that's not the photon itself.
2113  int iSpec;
2114  while (true) {
2115  iSpec = uchSpecList[rndmPtr->flat()*uchSpecList.size()];
2116  if (iSpec != iPhot) break;
2117  }
2118  eleVec.push_back(QEDsplitElemental(event, iPhot, iSpec));
2119  eleVec.back().ariWeight = 1.;
2120
2121  // Else use charged spectators.
2122  } else {
2123  double ariNorm = 0;
2124  vector<QEDsplitElemental> tempEleVec;
2125  for (int j = 0; j < (int)chSpecList.size(); j++) {
2126  int iSpec = chSpecList[j];
2127  tempEleVec.push_back(QEDsplitElemental(event, iPhot, iSpec));
2128  ariNorm += 1./tempEleVec.back().m2Ant;
2129  }
2130  // Set up Ariadne factors.
2131  for (int j = 0; j < (int)tempEleVec.size(); j++)
2132  tempEleVec[j].ariWeight = 1./(tempEleVec[j].m2Ant*ariNorm);
2133  eleVec.insert(eleVec.end(), tempEleVec.begin(), tempEleVec.end());
2134  }
2135  }
2136 }
2137
2138 //--------------------------------------------------------------------------
2139
2140 // Generate a scale for the system.
2141
2142 double QEDsplitSystem::generateTrialScale(Event &event, double q2Start) {
2143
2144  // Return saved trial.
2145  if (hasTrial) return q2Trial;
2146
2147  // Check if there are any photons left.
2148  if (eleVec.size() == 0) return 0;
2149
2150  // Starting scale - account for cut on mGammaMax.
2151  q2Trial = min(q2Max, q2Start);
2152
2153  // Check if qTrial is below the cutoff.
2154  if (q2Trial <= q2Cut) return 0;
2155
2156  // Find lower value from evolution window.
2157  int iEvol = evolutionWindows.size() - 1;
2158  while (q2Start <= evolutionWindows[iEvol]) iEvol--;
2159  double q2Low = evolutionWindows[iEvol];
2160
2161  // Compute weights.
2162  vector<double> weightVec;
2163  double totWeight(0), maxWeight(0);
2164  for (int i = 0; i < (int)eleVec.size(); i++) {
2165  double Iz = q2Low > eleVec[i].m2Ant ? 0 : 1. - q2Low/eleVec[i].m2Ant;
2166  double w = totIdWeight*eleVec[i].ariWeight*Iz*eleVec[i].getKallen();
2167  weightVec.push_back(w);
2168  totWeight += w;
2169  if (w > maxWeight) maxWeight = w;
2170  }
2171
2172  // If no antennae are active, don't generate new scale.
2173  if (totWeight < TINY) q2Trial = 0;
2174
2175  // Generate scale and do alpha veto.
2176  else {
2177  while (true) {
2178  double alphaMax = al.alphaEM(q2Trial);
2179  q2Trial *= pow(rndmPtr->flat(), M_PI/totWeight/alphaMax);
2180  double alphaNew = al.alphaEM(q2Trial);
2181  if (rndmPtr->flat() < alphaNew/alphaMax) break;
2182  }
2183  }
2184
2185  // Check if evolution window was crossed.
2186  if (q2Trial < q2Low) {
2187  if (iEvol == 0) return 0;
2188  return generateTrialScale(event, q2Low);
2189  }
2190
2191  // Select antenna.
2192  while (true) {
2193  int i = rndmPtr->flat()*weightVec.size();
2194  if (rndmPtr->flat() < weightVec[i]/maxWeight) {
2195  eleTrial = &eleVec[i];
2196  break;
2197  }
2198  }
2199
2200  // Select splitting ID.
2201  while (true) {
2202  int idIndex = rndmPtr->flat()*ids.size();
2203  idTrial = ids[idIndex];
2204  if (rndmPtr->flat() < idWeights[idIndex]/maxIdWeight) break;
2205  }
2206
2207  // Generate value of zeta and phi.
2208  zTrial = (1. - q2Low/eleTrial->m2Ant)*rndmPtr->flat();
2209  phiTrial = rndmPtr->flat()*2*M_PI;
2210  hasTrial = true;
2211  return q2Trial;
2212
2213 }
2214
2215 //--------------------------------------------------------------------------
2216
2217 // Check the veto.
2218
2219 bool QEDsplitSystem::checkVeto(Event &event) {
2220
2221  // Mark trial as used
2222  if (verbose >= debug) printOut(__METHOD_NAME__, "begin --------------");
2223  hasTrial = false;
2224
2225  // Set up some shorthands.
2226  int iPhot = eleTrial->iPhot;
2227  int iSpec = eleTrial->iSpec;
2228  double m2Ant = eleTrial->m2Ant;
2229
2230  // New momenta.
2231  vector<Vec4> pOld, pNew;
2232  pOld.push_back(event[iPhot].p());
2233  pOld.push_back(event[iSpec].p());
2234
2235  // ij is the new pair, k is the spectator.
2236  double mFerm = particleDataPtr->m0(idTrial);
2237  double mSpec = sqrt(eleTrial->m2Spec);
2238  double sIK = m2Ant - 2*pow2(mFerm) - pow2(mSpec);
2239  double sij = q2Trial - 2*pow2(mFerm);
2240  double sjk = zTrial*m2Ant;
2241  double sik = m2Ant - sij - sjk - 2*pow2(mFerm) - pow2(mSpec);
2242
2243  // Check phase space.
2244  if (sik < 0) return false;
2245  if (sij*sjk*sik - pow2(sij)*pow2(mSpec)
2246  - (pow2(sjk) + pow2(sik))*pow2(mFerm) < 0) return false;
2247
2248  // Kernel veto.
2249  double pVeto = ( (pow2(sik) + pow2(sjk))/m2Ant + 2.*pow2(mFerm)/q2Trial)/2.;
2250  if (rndmPtr->flat() > pVeto) return false;
2251  vector<double> invariants;
2252  invariants.push_back(sIK);
2253  invariants.push_back(sij);
2254  invariants.push_back(sjk);
2255  vector<double> masses;
2256  masses.push_back(mFerm);
2257  masses.push_back(mFerm);
2258  masses.push_back(mSpec);
2259
2260  // Kinematics.
2261  if (!vinComPtr->map2to3FF(pNew, pOld, 3, invariants, phiTrial, masses))
2262  return false;
2263
2264  // Set up the new fermions. Stochastic colour tag.
2265  int colTag = idTrial < 10 ? 10*(event.nextColTag()/10 + 1) + 1
2266  + rndmPtr->flat()*10 : 0;
2267  Particle partFermNew(idTrial, 51, iPhot, 0, 0, 0, colTag, 0, pNew[0], mFerm);
2268  Particle partAFermNew(-idTrial,51,iPhot, 0, 0, 0, 0, colTag, pNew[1], mFerm);
2269  Particle partSpecNew = event[iSpec];
2270  partSpecNew.mothers(iSpec, iSpec);
2271  partSpecNew.p(pNew[2]);
2272  partSpecNew.statusCode(52);
2273
2274  // Change the event - add new particles.
2275  int iFermNew = event.append(partFermNew);
2276  int iAFermNew = event.append(partAFermNew);
2277  int iSpecNew = event.append(partSpecNew);
2278
2280  event[iPhot].statusNeg();
2281  event[iPhot].daughters(iFermNew, iAFermNew);
2282  event[iSpec].statusNeg();
2283  event[iSpec].daughters(iSpecNew, 0);
2284
2285  // Update parton system.
2286  partonSystemsPtr->replace(iSys, iPhot, iFermNew);
2288  partonSystemsPtr->replace(iSys, iSpec, iSpecNew);
2289  event.restorePtrs();
2290  if (verbose >= debug) printOut(__METHOD_NAME__, "end --------------");
2291  return true;
2292
2293 }
2294
2295 //--------------------------------------------------------------------------
2296
2297 // Print the system.
2298
2299 void QEDsplitSystem::print() {
2300  cout << "Splitting" << endl;
2301  for (int i = 0; i < (int)eleVec.size(); i++)
2302  cout << "(" << eleVec[i].iPhot << " " << eleVec[i].iSpec << ") "
2303  << "s = " << eleVec[i].m2Ant << " ariFac = " << eleVec[i].ariWeight
2304  << endl;
2305 }
2306
2307 //==========================================================================
2308
2309 // Class for a QED conversion system.
2310
2311 //--------------------------------------------------------------------------
2312
2313 // Initialize the pointers.
2314
2315 void QEDconvSystem::initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) {
2316  infoPtr = infoPtrIn;
2317  particleDataPtr = infoPtr->particleDataPtr;
2318  partonSystemsPtr = infoPtr->partonSystemsPtr;
2319  rndmPtr = infoPtr->rndmPtr;
2320  settingsPtr = infoPtr->settingsPtr;
2321  vinComPtr = vinComPtrIn;
2322  isInitPtr = true;
2323 }
2324
2325 //--------------------------------------------------------------------------
2326
2327 // Initialize the system.
2328 void QEDconvSystem::init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
2329  int verboseIn) {
2330
2331  // Verbosity setting.
2332  if (!isInitPtr) printOut(__METHOD_NAME__, "initPtr not called");
2333  verbose = verboseIn;
2334
2335  // Settings, number of incoming flavours to allow conversions to.
2336  // Could be extended to allow top quarks in future; for now up to b.
2337  nQuark = 5;
2338  if (!settingsPtr->flag("Vincia:convertGammaToQuark")) nQuark = 0;
2339
2340  // Set trial pdf ratios.
2341  Rhat[1] = 77;
2342  Rhat[2] = 140;
2343  Rhat[3] = 30;
2344  Rhat[4] = 22;
2345  Rhat[5] = 15;
2346  Rhat[-1] = 63;
2347  Rhat[-2] = 65;
2348  Rhat[-3] = 30;
2349  Rhat[-4] = 30;
2350  Rhat[-5] = 16;
2351
2352  // Constants.
2353  TINYPDF = pow(10, -10);
2354
2355  // Beam pointers.
2356  beamAPtr = beamAPtrIn;
2357  beamBPtr = beamBPtrIn;
2358  isInit = true;
2359
2360 }
2361
2362 //--------------------------------------------------------------------------
2363
2364 // Prepare for backwards-evolution of photons.
2365
2366 void QEDconvSystem::prepare(int iSysIn, Event &event, double q2CutIn,
2367  bool isBelowHadIn, vector<double> evolutionWindowsIn, AlphaEM alIn) {
2368
2369  if (!isInit) {
2370  infoPtr->errorMsg("Error in "+__METHOD_NAME__+": Not initialised.");
2371  return;
2372  }
2373  if (verbose >= louddebug) printOut(__METHOD_NAME__, "begin --------------");
2374
2375  // Input.
2376  iSys = iSysIn;
2377  shh = infoPtr->s();
2379  q2Cut = q2CutIn;
2380  evolutionWindows = evolutionWindowsIn;
2381  al = alIn;
2382
2383  // Set up weights for conversion flavours.
2384  ids.clear();
2385  idWeights.clear();
2386  totIdWeight = 0;
2387  maxIdWeight = 0;
2388
2389  // If switched off, do nothing.
2390  if (nQuark == 0) return;
2391
2392  // Only do conversions to quarks if above hadronisation scale.
2394  for (int i = 1; i <= nQuark; i++) {
2395  ids.push_back(i);
2396  ids.push_back(-i);
2397  idWeights.push_back((i%2==0 ? 4./9. : 1./9.)*Rhat[i]);
2398  idWeights.push_back((i%2==0 ? 4./9. : 1./9.)*Rhat[-i]);
2399  }
2400  // Total weights.
2401  for (int i = 0; i < (int)idWeights.size(); i++) {
2402  totIdWeight += idWeights[i];
2403  if (idWeights[i] > maxIdWeight) maxIdWeight = idWeights[i];
2404  }
2405
2406  // Build system.
2407  buildSystem(event);
2408  if (verbose >= louddebug) printOut(__METHOD_NAME__, "end --------------");
2409
2410 }
2411
2412 //--------------------------------------------------------------------------
2413
2414 // Build the system.
2415
2416 void QEDconvSystem::buildSystem(Event &event) {
2417
2418  // Get rid of saved trial.
2419  hasTrial = false;
2420
2421  // Get initial states.
2422  iA = partonSystemsPtr->getInA(iSys);
2423  iB = partonSystemsPtr->getInB(iSys);
2424  isAPhot = event[iA].id() == 22;
2425  isBPhot = event[iB].id() == 22;
2426  s = (event[iA].p() + event[iB].p()).m2Calc();
2427
2428 }
2429
2430 //--------------------------------------------------------------------------
2431
2432 // Generate a trial scale.
2433
2434 double QEDconvSystem::generateTrialScale(Event &event, double q2Start) {
2435
2436  // Return if saved trial.
2437  if (hasTrial) return q2Trial;
2438
2439  // Check if there are any photons.
2440  if (!isAPhot && !isBPhot) {return 0;}
2441  double totWeight = 1.;
2442
2443  // Select a photon.
2444  if (isAPhot && !isBPhot) {iPhotTrial = iA; iSpecTrial = iB;}
2445  else if (isBPhot && !isAPhot) {iPhotTrial = iB; iSpecTrial = iA;}
2446  else {
2447  if (rndmPtr->flat() < 0.5) {iPhotTrial = iA; iSpecTrial = iB;}
2448  else {iPhotTrial = iB; iSpecTrial = iA;}
2449  // Two photon antennae -> twice the weight.
2450  totWeight *= 2.;
2451  }
2452
2453  // Starting scale.
2454  q2Trial = q2Start;
2455
2456  // Check if qTrial is below the cutoff.
2457  if (q2Trial <= q2Cut) return 0;
2458
2459  // Find lower value from evolution window.
2460  int iEvol = evolutionWindows.size() - 1;
2461  while(q2Start <= evolutionWindows[iEvol]) {iEvol--;}
2462  double q2Low = evolutionWindows[iEvol];
2463
2464  // Iz integral.
2465  double zPlus = shh/s;
2466  double zMin = 1 + q2Low/s;
2467  if (zPlus < zMin) {return 0;}
2468  double Iz = log(zPlus/zMin);
2469  totWeight *= totIdWeight*Iz;
2470
2471  // If no antennae are active, don't generate new scale.
2472  if (totWeight < TINY) return 0;
2473
2474  // Generate scale and do alpha veto.
2475  else
2476  while(true) {
2477  double alphaMax = al.alphaEM(q2Trial);
2478  q2Trial *= pow(rndmPtr->flat(), M_PI/totWeight/alphaMax);
2479  double alphaNew = al.alphaEM(q2Trial);
2480  if (rndmPtr->flat() < alphaNew/alphaMax) break;
2481  }
2482
2483  // Check if evolution window was crossed.
2484  if (q2Trial < q2Low) {
2485  if (iEvol==0) return 0;
2486  return generateTrialScale(event, q2Low);
2487  }
2488
2489  // Select conversion ID.
2490  while( true) {
2491  int idIndex = rndmPtr->flat()*ids.size();
2492  idTrial = ids[idIndex];
2493  if (rndmPtr->flat() < idWeights[idIndex]/maxIdWeight) break;
2494  }
2495  zTrial = zMin*pow(zPlus/zMin, rndmPtr->flat());
2496  phiTrial = rndmPtr->flat()*2*M_PI;
2497  hasTrial = true;
2498  return q2Trial;
2499
2500 }
2501
2502 //--------------------------------------------------------------------------
2503
2504 // Check the veto.
2505
2506 bool QEDconvSystem::checkVeto(Event &event) {
2507
2508  // Mark trial as used
2509  if (verbose >= debug) printOut(__METHOD_NAME__, "begin --------------");
2510  hasTrial = false;
2511
2512  // Conversion mass.
2513  double mf2 = pow2(particleDataPtr->m0(idTrial));
2514
2515  // Spectator ID.
2516  int idSpec = event[iSpecTrial].id();
2517
2518  // New momenta.
2519  vector<Vec4> pOld, pNew;
2520  pOld.push_back(event[iPhotTrial].p());
2521  pOld.push_back(event[iSpecTrial].p());
2522
2523  // Note that we treat the initial state as massless and the final
2524  // state as massive. q2 is defined as saj - 2*mf2, but otherwise we
2525  // adhere to the proper kinematics.
2526  double saj = q2Trial + 2*mf2;
2527  double sbj = (zTrial-1)*s - saj + mf2;
2528  double sab = s + saj + sbj - mf2;
2529
2530  // Check kinematic boundaries.
2531  if (sbj < 0) return false;
2532  if (sab*saj*sbj - mf2*sab*sab < 0) return false;
2533
2534  // Check if photon is in beam a or b.
2535  bool isPhotA = (iPhotTrial == iA) ? true : false;
2536
2537  // Global recoil momenta.
2538  vector<Vec4> pRec;
2539  vector<int> iRec;
2540  int sysSize = partonSystemsPtr->sizeAll(iSys);
2541  for (int i=0; i<sysSize; i++) {
2542  int iEv = partonSystemsPtr->getAll(iSys, i);
2543  if (iEv < 0 || !event[iEv].isFinal()) continue;
2544  pRec.push_back(event[iEv].p());
2545  iRec.push_back(iEv);
2546  }
2547
2548  // Kinematics.
2549  if (!vinComPtr->map2to3II(pNew, pRec, pOld, s, saj, sbj, sab,
2550  phiTrial, mf2)) return false;
2551
2552  // Check if new energies don't exceed hadronic maxima.
2553  double eaUsed = 0, ebUsed = 0;
2554  int nSys = partonSystemsPtr->sizeSys();
2555  for (int i=0; i<nSys; i++) {
2556  eaUsed += event[partonSystemsPtr->getInA(i)].e();
2557  ebUsed += event[partonSystemsPtr->getInB(i)].e();
2558  }
2559  if (isPhotA) {
2560  if ((eaUsed - pOld[0].e() + pNew[0].e()) > 0.98*sqrt(shh)/2.) return false;
2561  if ((ebUsed - pOld[1].e() + pNew[2].e()) > 0.98*sqrt(shh)/2.) return false;
2562  } else {
2563  if ((ebUsed - pOld[0].e() + pNew[0].e()) > 0.98*sqrt(shh)/2.) return false;
2564  if ((eaUsed - pOld[1].e() + pNew[2].e()) > 0.98*sqrt(shh)/2.) return false;
2565  }
2566
2567  // Kernel veto probability.
2568  double pVeto = 0.5*(1. + pow2(sbj)/pow2(sab)
2569  - 2.*mf2*pow2(s)/pow2(sab)/(saj-2*mf2));
2570
2571  // Compute pdf ratios.
2572  double Rpdf = 1.;
2573  double xPhotOld = pOld[0].e()/(sqrt(shh)/2.);
2574  double xPhotNew = pNew[0].e()/(sqrt(shh)/2.);
2575  double xSpecOld = pOld[1].e()/(sqrt(shh)/2.);
2576  double xSpecNew = pNew[2].e()/(sqrt(shh)/2.);
2577
2578  if (isPhotA) {
2579  // Photon pdf.
2580  double newPDFPhot = beamAPtr->xfISR(iSys, idTrial, xPhotNew, q2Trial);
2581  double oldPDFPhot = beamAPtr->xfISR(iSys, 22, xPhotOld, q2Trial);
2582  if (abs(newPDFPhot) < TINYPDF) newPDFPhot = TINYPDF;
2583  if (abs(oldPDFPhot) < TINYPDF) oldPDFPhot = TINYPDF;
2584  Rpdf *= newPDFPhot/oldPDFPhot;
2585
2586  // Spectator pdf.
2587  double newPDFSpec = beamBPtr->xfISR(iSys, idSpec, xSpecNew, q2Trial);
2588  double oldPDFSpec = beamBPtr->xfISR(iSys, idSpec, xSpecOld, q2Trial);
2589  if (abs(newPDFSpec) < TINYPDF) newPDFSpec = TINYPDF;
2590  if (abs(oldPDFSpec) < TINYPDF) oldPDFSpec = TINYPDF;
2591  Rpdf *= newPDFSpec/oldPDFSpec;
2592  } else {
2593  // Photon pdf.
2594  double newPDFPhot = beamBPtr->xfISR(iSys, idTrial, xPhotNew, q2Trial);
2595  double oldPDFPhot = beamBPtr->xfISR(iSys, 22, xPhotOld, q2Trial);
2596  if (abs(newPDFPhot) < TINYPDF) newPDFPhot = TINYPDF;
2597  if (abs(oldPDFPhot) < TINYPDF) oldPDFPhot = TINYPDF;
2598  Rpdf *= newPDFPhot/oldPDFPhot;
2599
2600  // Spectator pdf.
2601  double newPDFSpec = beamAPtr->xfISR(iSys, idSpec, xSpecNew, q2Trial);
2602  double oldPDFSpec = beamAPtr->xfISR(iSys, idSpec, xSpecOld, q2Trial);
2603  if (abs(newPDFSpec) < TINYPDF) newPDFSpec = TINYPDF;
2604  if (abs(oldPDFSpec) < TINYPDF) oldPDFSpec = TINYPDF;
2605  Rpdf *= newPDFSpec/oldPDFSpec;
2606  }
2607
2608  if (Rpdf > Rhat[idTrial]) {
2609  stringstream ss;
2610  ss << "Incorrect pdf overestimate " << "id = " << idTrial << "ratio = "
2611  << Rpdf/Rhat[idTrial];
2612  printOut(__METHOD_NAME__, ss.str());
2613  }
2614
2615  // Pdf ratio veto probability.
2616  pVeto *= (Rpdf/Rhat[idTrial]);
2617
2618  // Do veto.
2619  if (rndmPtr->flat() > pVeto) return false;
2620
2621  // Passed veto, set up new particles.
2622  Particle partSpecNew = event[iSpecTrial];
2623  partSpecNew.p(pNew[2]);
2624  partSpecNew.statusCode(41);
2625
2626  // Stochastic colour tag.
2627  int colTag = 10*(event.nextColTag()/10 + 1) + 1 + rndmPtr->flat()*10;
2628  Particle partBeamNew (idTrial, -41, 0, 0, 0, 0, idTrial > 0 ?
2629  colTag : 0, idTrial > 0 ? 0 : colTag, pNew[0]);
2630  Particle partFinalNew (idTrial, 43, 0, 0, 0, 0, idTrial > 0 ?
2631  colTag : 0, idTrial > 0 ? 0 : colTag, pNew[1], sqrt(mf2));
2632  int iBeamNew = event.append(partBeamNew);
2633  int iFinalNew = event.append(partFinalNew);
2634  int iSpecNew = event.append(partSpecNew);
2635  event[iPhotTrial].statusNeg();
2636  event[iSpecTrial].statusNeg();
2637  event[iBeamNew].mothers(event[iPhotTrial].mother1(),
2638  event[iPhotTrial].mother2());
2639  event[iFinalNew].mothers(iBeamNew, 0);
2640  event[iSpecNew].mothers(event[iSpecTrial].mother1(),
2641  event[iSpecTrial].mother2());
2642  event[iPhotTrial].mothers(iBeamNew, 0);
2643  event[iSpecTrial].mothers(iSpecNew, 0);
2644  event[iBeamNew].daughters(iFinalNew, iPhotTrial);
2645  event[iFinalNew].daughters(0,0);
2646  event[iSpecNew].daughters(iSpecTrial);
2647
2648  // Change daughters of beams for hard process.
2649  if (iSys == 0) {
2650  bool foundPhot = false;
2651  bool foundSpec = false;
2652  for (int i = 0; i < (int)event.size(); i++) {
2653  if (!foundPhot)
2654  if (event[i].daughter1() == iPhotTrial) {
2655  event[i].daughters(iBeamNew, 0);
2656  foundPhot = true;
2657  }
2658  if (!foundSpec)
2659  if (event[i].daughter1() == iSpecTrial) {
2660  event[i].daughters(iSpecNew, 0);
2661  foundSpec = true;
2662  }
2663  if (foundPhot && foundSpec) break;
2664  }
2665  }
2666
2667  // Update event for global recoil.
2668  vector<pair<int,int> > iRecNew;
2669  iRecNew.clear();
2670  iRecNew.resize(0);
2671  for (int j = 0; j < event.size(); j++)
2672  if (event[j].isFinal())
2673  for (int k = 0; k < (int)iRec.size(); k++)
2674  if (iRec[k] == j) {
2675  // Copy the recoiler and change the momentum.
2676  int inew = event.copy(j,44);
2677  event[inew].p(pRec[k]);
2678  iRecNew.push_back(make_pair(iRec[k], inew));
2679  }
2680
2681  // Update parton system.
2682  partonSystemsPtr->replace(iSys, iPhotTrial, iBeamNew);
2684  partonSystemsPtr->replace(iSys, iSpecTrial, iSpecNew);
2685
2686  // Recoilers and initial particles.
2687  for (int k=0; k<(int)iRecNew.size(); k++) {
2688  partonSystemsPtr->replace(iSys, iRecNew[k].first, iRecNew[k].second);
2689  }
2690  if (isPhotA) {
2691  partonSystemsPtr->setInA(iSys, iBeamNew);
2692  partonSystemsPtr->setInB(iSys, iSpecNew);
2693  } else {
2694  partonSystemsPtr->setInA(iSys, iSpecNew);
2695  partonSystemsPtr->setInB(iSys, iBeamNew);
2696  }
2697  double shat = (event[partonSystemsPtr->getInA(iSys)].p() +
2698  event[partonSystemsPtr->getInB(iSys)].p()).m2Calc();
2699  partonSystemsPtr->setSHat(iSys, shat);
2700
2701  // Update beams.
2702  BeamParticle& beam1 = *beamAPtr;
2703  BeamParticle& beam2 = *beamBPtr;
2704  if (isPhotA) {
2705  beam1[iSys].update(iBeamNew, event[iBeamNew].id(),
2706  event[iBeamNew].e()/beam1.e());
2707  beam2[iSys].update(iSpecNew, event[iSpecNew].id(),
2708  event[iSpecNew].e()/beam2.e());
2709  } else {
2710  beam1[iSys].update(iSpecNew, event[iSpecNew].id(),
2711  event[iSpecNew].e()/beam1.e());
2712  beam2[iSys].update(iBeamNew, event[iBeamNew].id(),
2713  event[iBeamNew].e()/beam2.e());
2714  }
2715  if (verbose >= debug) printOut(__METHOD_NAME__, "end --------------");
2716  return true;
2717
2718 }
2719
2720 //--------------------------------------------------------------------------
2721
2722 // Print.
2723
2724 void QEDconvSystem::print() {
2725  cout << "Conversion" << endl;
2726  cout << "s = " << s << endl;
2727 }
2728
2729 //==========================================================================
2730
2731 // Class for performing QED showers.
2732
2733 //--------------------------------------------------------------------------
2734
2735 // Initialize the pointers.
2736
2737 void QEDShower::initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) {
2738  infoPtr = infoPtrIn;
2739  particleDataPtr = infoPtr->particleDataPtr;
2740  partonSystemsPtr = infoPtr->partonSystemsPtr;
2741  rndmPtr = infoPtr->rndmPtr;
2742  settingsPtr = infoPtr->settingsPtr;
2743  vinComPtr = vinComPtrIn;
2744  isInitPtr = true;
2745 }
2746
2747 //--------------------------------------------------------------------------
2748
2749 // Initialize settings for current run.
2750
2751 void QEDShower::init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn) {
2752
2753  // Initialize alphaEM
2754  verbose = settingsPtr->mode("Vincia:verbose");
2755  double alpEM0Vincia = settingsPtr->parm("Vincia:alphaEM0");
2756  double alpEMmzVincia = settingsPtr->parm("Vincia:alphaEMmz");
2757  double alpEM0Pythia = settingsPtr->parm("StandardModel:alphaEM0");
2758  double alpEMmzPythia = settingsPtr->parm("StandardModel:alphaEMmZ");
2759  int alphaEMorder = settingsPtr->mode("Vincia:alphaEMorder");
2760
2761  // Change Pythia settings, initialize, then change them back.
2762  settingsPtr->parm("StandardModel:alphaEM0", alpEM0Vincia);
2763  settingsPtr->parm("StandardModel:alphaEMmZ", alpEMmzVincia);
2764  al.init(alphaEMorder, settingsPtr);
2765  settingsPtr->parm("StandardModel:alphaEM0", alpEM0Pythia);
2766  settingsPtr->parm("StandardModel:alphaEMmz", alpEMmzPythia);
2767
2768  // Get settings.
2769  doQED = settingsPtr->flag("Vincia:doQED");
2770  doEmission = doQED;
2771  nGammaToLepton = settingsPtr->mode("Vincia:nGammaToLepton");
2772  nGammaToQuark = settingsPtr->mode("Vincia:nGammaToQuark") >= 1;
2773  doConvertGamma = settingsPtr->flag("Vincia:convertGammaToQuark");
2774  doConvertQuark = settingsPtr->flag("Vincia:convertQuarkToGamma");
2775
2776  // QED cutoff for coloured particles/hadronisation scale.
2777  q2minColouredSav = pow2(settingsPtr->parm("Vincia:QminChgQ"));
2778  q2minSav = pow2(settingsPtr->parm("Vincia:QminChgL"));
2779
2780  // Set beam pointers.
2781  beamAPtr = beamAPtrIn;
2782  beamBPtr = beamBPtrIn;
2783  isInitSav = true;
2784
2785 }
2786
2787 //--------------------------------------------------------------------------
2788
2789 // Prepare to shower a system.
2790
2791 void QEDShower::prepare(int iSysIn, Event &event, bool isBelowHad) {
2792
2793  // Verbose output
2794  if (!doQED) return;
2795  if (verbose >= debug) {
2796  printOut(__METHOD_NAME__, "begin --------------");
2797  if (verbose >= superdebug) event.list();
2798  }
2799
2800  // Above or below hadronisation scale.
2801  double q2cut = (isBelowHad) ? q2minSav : q2minColouredSav;
2802
2803  // Initialize windows for the hard systen and
2804  // the final after beam remnants system.
2805  if (iSysIn == 0 || iSysIn == -1) {
2806  // The cutoff scale is the lowest window boundary, then step up to
2807  // q2Max successively by factor winStep.
2808  double q2BiggestEver = infoPtr->s();
2809  double q2Window = q2cut;
2810  double winStep = 100.0;
2811  evolutionWindows.clear();
2812  do {
2813  evolutionWindows.push_back(q2Window);
2814  q2Window *= winStep;
2815  } while(q2Window < q2BiggestEver);
2816  }
2817
2818  // Find out if showering a resonance decay.
2819  bool isResDecay = false;
2820  if (iSysIn >= 0)
2821  if (partonSystemsPtr->hasInRes(iSysIn) ||
2822  (partonSystemsPtr->getInA(iSysIn)==0 &&
2823  partonSystemsPtr->getInB(iSysIn)==0) ) isResDecay = true;
2824
2825  // If showering a resonance decay or below hadronization scale clear
2826  // out all previous systems.
2827  if ( iSysIn <= 0 || isResDecay || isBelowHad ) {
2828  iSystems.clear();
2829  emitSystems.clear();
2830  splitSystems.clear();
2831  convSystems.clear();
2832  }
2833
2834  // Special case: iSysIn = -1 implies below hadronisation scale.
2835  // Collect all final-state particles into one new system. Note:
2836  // this system will have sum(charge) != 0 if the sum of the beam
2837  // charges is nonzero.
2838  if (iSysIn == -1) {
2840  // Loop over all particles in event rather than all parton systems
2841  // since beam remnant partons not part of any partonSystem.
2842  for (int i = 1; i < event.size(); ++i) {
2843  if (!event[i].isFinal()) continue;
2845  }
2846  }
2847
2849  iSystems.push_back(iSysIn);
2850  emitSystems.push_back(QEDemitSystem());
2851  splitSystems.push_back(QEDsplitSystem());
2852  convSystems.push_back(QEDconvSystem());
2853
2854  // Initialze pointers.
2855  emitSystems.back().initPtr(infoPtr, vinComPtr);
2856  splitSystems.back().initPtr(infoPtr, vinComPtr);
2857  convSystems.back().initPtr(infoPtr, vinComPtr);
2858
2859  // Initialize. Note, these calls read from settings database, should
2860  // be rewritten to require settings to be read only once.
2861  emitSystems.back().init(beamAPtr, beamBPtr,verbose);
2862  splitSystems.back().init(beamAPtr, beamBPtr,verbose);
2863  convSystems.back().init(beamAPtr, beamBPtr,verbose);
2864
2865  // Prepare and build QED systems.
2867  evolutionWindows, al);
2869  evolutionWindows, al);
2871  evolutionWindows, al);
2872
2873  // Verbose output.
2874  if(verbose >= debug) printOut(__METHOD_NAME__, "end --------------");
2875
2876 }
2877
2878 //--------------------------------------------------------------------------
2879
2880 // Update QED shower system(s) each time something has changed in event.
2881
2882 void QEDShower::update(Event &event, int iSys) {
2883
2884  // Find index of the system.
2885  if (verbose >= debug) printOut(__METHOD_NAME__, "begin --------------");
2886  for (int i = 0; i < (int)iSystems.size(); i++) {
2887  if (iSystems[i] == iSys) {
2888  // Rebuild systems.
2889  emitSystems[i].buildSystem(event);
2890  splitSystems[i].buildSystem(event);
2891  convSystems[i].buildSystem(event);
2892  break;
2893  }
2894  }
2895  if (verbose >= debug) {
2896  if (verbose >= superdebug) event.list();
2897  printOut(__METHOD_NAME__, "end --------------");
2898  }
2899 }
2900
2901 //--------------------------------------------------------------------------
2902
2903 // Generate a trial scale.
2904
2905 double QEDShower::generateTrialScale(Event &event, double q2Start) {
2906
2907  // Get a trial from every system.
2908  q2Trial = 0;
2909  isTrialEmit = false;
2910  isTrialSplit = false;
2911  isTrialConv = false;
2912  if (!doQED) return 0.0;
2913  if (verbose >= louddebug) {
2914  printOut(__METHOD_NAME__, "begin --------------");
2915  if (verbose >= superdebug)
2916  cout << " QEDShower::generateTrialScale() q2Start = " << q2Start
2917  << " doEmit = " << bool2str(doEmission)
2918  << " nSplitGamToLep = " << num2str(nGammaToLepton)
2919  << " nSplitGamToQuark = " << num2str(nGammaToQuark)
2920  << " doConv = " << bool2str(doConvertGamma) << endl;
2921  }
2922
2923  // Emissions.
2924  if (doEmission) {
2925  for (int i = 0; i < (int)emitSystems.size(); i++) {
2926  double q2TrialEmitNew =
2927  emitSystems[i].generateTrialScale(event, q2Start);
2928  if (q2TrialEmitNew > q2Trial) {
2929  q2Trial = q2TrialEmitNew;
2930  iSysTrial = iSystems[i];
2931  iSysIndexTrial = i;
2932  isTrialEmit = true;
2933  isTrialSplit = false;
2934  isTrialConv = false;
2935  }
2936  }
2937  }
2938
2939  // Splittings.
2940  if (nGammaToLepton + nGammaToQuark > 0) {
2941  for (int i = 0; i < (int)splitSystems.size(); i++) {
2942  double q2TrialSplitNew =
2943  splitSystems[i].generateTrialScale(event, q2Start);
2944  if (q2TrialSplitNew > q2Trial) {
2945  q2Trial = q2TrialSplitNew;
2946  iSysTrial = iSystems[i];
2947  iSysIndexTrial = i;
2948  isTrialEmit = false;
2949  isTrialSplit = true;
2950  isTrialConv = false;
2951  }
2952  }
2953  }
2954
2955  // Conversions.
2956  if (doConvertGamma) {
2957  for (int i = 0; i < (int)convSystems.size(); i++) {
2958  double q2TrialConvNew =
2959  convSystems[i].generateTrialScale(event, q2Start);
2960  if (q2TrialConvNew > q2Trial) {
2961  q2Trial = q2TrialConvNew;
2962  iSysTrial = iSystems[i];
2963  iSysIndexTrial = i;
2964  isTrialEmit = false;
2965  isTrialSplit = false;
2966  isTrialConv = true;
2967  }
2968  }
2969  }
2970  if (verbose >= louddebug) printOut(__METHOD_NAME__, "end --------------");
2971  return q2Trial;
2972
2973 }
2974
2975 //--------------------------------------------------------------------------
2976
2977 // Check the veto.
2978
2979 bool QEDShower::checkVeto(Event &event) {
2980  if (verbose >= debug) printOut(__METHOD_NAME__, "begin --------------");
2981  bool doVeto = false;
2982  if (isTrialEmit) doVeto = emitSystems[iSysIndexTrial].checkVeto(event);
2983  if (isTrialSplit) doVeto = splitSystems[iSysIndexTrial].checkVeto(event);
2984  if (isTrialConv) doVeto = convSystems[iSysIndexTrial].checkVeto(event);
2985  if (verbose >= debug) printOut(__METHOD_NAME__, "end --------------");
2986  return doVeto;
2987 }
2988
2989 //==========================================================================
2990
2991 } // end namespace Pythia8
Definition: beam.h:43
Definition: AgUStep.h:26