StRoot  1
VinciaCommon.cc
1 // VinciaCommon.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 // VinciaCommon classes, and related auxiliary methods.
8
9 #include "Pythia8/VinciaCommon.h"
10
11 namespace Pythia8 {
12
13 //==========================================================================
14
15 // Rambo phase space generator.
16
17 //--------------------------------------------------------------------------
18
19 // Massless flat phase space generator. Generate a random (uniformly
20 // distributed) massless PS point with nOut particles and total sqrt(s) = eCM.
21
22 double Rambo::genPoint(double eCM, int nOut, vector<Vec4>& pOut) {
23
24  // Set size of output vector
25  pOut.resize(nOut);
26  // Create momentum-sum four-vector
27  Vec4 R;
28  // Generate nParticles independent massless 4-momenta with isotropic angles
29  for (int i = 0; i < nOut; ++i) {
30  // Cos(theta), sin(theta), and phi
31  double c = 2.0*rndmPtr->flat() - 1.0;
32  double s = sqrt(1.0-pow2(c));
33  double phi = 2.0*M_PI*rndmPtr->flat();
34  // Norm
35  double r12 = 0.0;
36  while (r12 == 0.0) {
37  double r1 = rndmPtr->flat();
38  double r2 = rndmPtr->flat();
39  r12 = r1*r2;
40  }
41  double En = -log(r12);
42  pOut[i].e(En);
43  pOut[i].pz(En*c);
44  pOut[i].py(En*s*cos(phi));
45  pOut[i].px(En*s*sin(phi));
47  R += pOut[i];
48  }
49  // Compute ECM and normalise to unity (with sign flip)
50  double Rmass = R.mCalc();
51  R /= -Rmass;
52  // Transform momenta so add up to (eCM, 0, 0, 0)
53  double a = 1.0/(1.0-R.e());
54  double x = eCM/Rmass;
55  for (int i = 0; i < nOut; ++i) {
56  double bq = dot3(R, pOut[i]);
57  pOut[i].px( x * (pOut[i].px()+R.px()*(pOut[i].e()+a*bq)) );
58  pOut[i].py( x * (pOut[i].py()+R.py()*(pOut[i].e()+a*bq)) );
59  pOut[i].pz( x * (pOut[i].pz()+R.pz()*(pOut[i].e()+a*bq)) );
60  pOut[i].e( x * (-R.e()*pOut[i].e()+bq) );
61  }
62  // The weight is always unity for the massless algorithm.
63  return 1.0;
64
65 }
66
67 //--------------------------------------------------------------------------
68
69 // Massive flat phase space generator, generalised according to the
70 // original paper. The momenta are not distributed flat in phase
71 // space anymore, returns the weight of the phase space configutation.
72
73 double Rambo::genPoint(double eCM, vector<double> mIn, vector<Vec4>& pOut) {
74
75  // Call the massless genPoint, initializing weight.
76  int nOut = mIn.size();
77  if (nOut <= 1 || eCM <= 0.) return 0.;
78  double weight = genPoint(eCM, nOut, pOut);
79  bool massesnonzero = false;
80
81  // Set up the function determining the rescaling parameter xi.
82  vector<double> energies;
83  for (int i = 0; i < nOut; i++) {
84  energies.push_back(pOut[i].e());
85  if (pow2(mIn[i]/eCM) > TINY) massesnonzero = true;
86  }
87
88  // If none of the reduced masses is > TINY, return
89  if (!massesnonzero) return weight;
90
91  // Rescale all the momenta.
92  TXiFunctor rhs = TXiFunctor(mIn, energies);
93  double xi = zbrent(rhs, eCM, 0., 1., 1e-10);
94  for (int iMom = 0; iMom < nOut; iMom++) {
95  pOut[iMom].rescale3(xi);
96  pOut[iMom].e( sqrt(pow2(mIn[iMom]) + pow2(xi)*pow2(pOut[iMom].e())) );
97  }
98
99  // Determine the quantities needed for the calculation of the weight.
100  double sumP(0.), prodPdivE(1.), sumP2divE(0.);
101  for (int iMom = 0; iMom < nOut; iMom++) {
102  double pAbs2 = pOut[iMom].pAbs2();
103  double pAbs = sqrt(pAbs2);
104  sumP += pAbs;
105  prodPdivE *= pAbs/pOut[iMom].e();
106  sumP2divE += pAbs2/pOut[iMom].e();
107  }
108
109  // There's a typo in eq. 4.11 of the Rambo paper by Kleiss, Stirling
110  // and Ellis, the Ecm below is not present there.
111  weight *= pow(sumP/eCM,2*nOut-3)*prodPdivE*eCM/sumP2divE;
112  return weight;
113 }
114
115 //==========================================================================
116
117 // The Colour class.
118
119 //--------------------------------------------------------------------------
120
121 // Initialize.
122
123 bool Colour::init() {
124
125  // Sanity check.
126  if (!isInitPtr) return false;
127
128  // Set verbosity level.
129  verbose = settingsPtr->mode("Vincia:verbose");
130
131  // Set parameters. CR disabled in this version.
132  inheritMode = settingsPtr->mode("Vincia:CRinheritMode");
133
134  // Sett initialization.
135  isInit = true;
136  return isInit;
137
138 }
139
140 //--------------------------------------------------------------------------
141
142 // Map a list of particles with ordinary Les Houches colour tags into
143 // a new set of Les Houches colour tags with the last index running
144 // from 1-9.
145 //
146 // Note on coloured resonances. These should have their colours
147 // defined by whatever system they were produced by, and should not
148 // be recoloured when decaying. This means that only colour lines for
149 // which both the colour AND the corresponding anticolour line are
150 // found inside the same system should be redefined (assuming the
151 // coloured resonance itself does not appear explicitly in its
152 // corresponding decay partonSystem.
153
154 bool Colour::colourise(int iSys, Event& event) {
155
156  // Sanity checks.
157  if (!isInit) {
158  printOut("Colour::colourise","ERROR! Colour not initialised");
159  return false;
160  }
161  else if (partonSystemsPtr->sizeAll(iSys) <= 1) return false;
162
163  // Construct colour map and assign new tags.
164  map<int, int> colourMap;
165  int startTag = event.lastColTag();
166  int nTries = 0;
167  bool accept = false;
168
169  // Do not recolour decaying resonances (as they will already have
170  // been coloured when they were produced).
171  int colRes(0), acolRes(0);
172  if (partonSystemsPtr->hasInRes(iSys)) {
173  int iRes = partonSystemsPtr->getInRes(iSys);
174  colRes = event[iRes].col();
175  acolRes = event[iRes].acol();
176  }
177
178  // Reject assignments that would produce (subleading) singlet gluons.
179  while (!accept && ++nTries < 10) {
180  colourMap.clear();
181
182  // Example: if startTag is 220-229, nextTagBase is 230.
183  int nextTagBase = 10*int((startTag/10)+1);
184  for (int i=0; i<partonSystemsPtr->sizeAll(iSys); ++i) {
185  int i1 = partonSystemsPtr->getAll(iSys,i);
186  if (i1 <= 0) continue;
187  Particle* partonPtr = &event[i1];
188  if (partonPtr->colType() == 0) continue;
189  int col, acol;
190
191  // Cross initial-state colours.
192  if (i < partonSystemsPtr->sizeAll(iSys)
193  - partonSystemsPtr->sizeOut(iSys)) {
194  acol = partonPtr->col();
195  col = partonPtr->acol();
196  if (col == acolRes) col = 0;
197  if (acol == colRes) acol = 0;
198  }
199  else {
200  col = partonPtr->col();
201  acol = partonPtr->acol();
202  if (col == colRes) col = 0;
203  if (acol == acolRes) acol = 0;
204  }
205  if (col == 0 && acol == 0) continue;
206  int colIndx = colourMap[col];
207  int acolIndx = colourMap[acol];
208  if (col != 0) {
209
210  // First time a tag is encountered, mark it negative -> one end found.
211  if (colIndx == 0) {
212  // Ensure gluons have different col and acol indices.
213  while (colIndx == 0 || colIndx == acolIndx) {
214  colIndx = nextTagBase + int(rndmPtr->flat()*9) + 1;
215  }
216  colourMap[col] = -colIndx;
217  }
218  // Second time mark it positive -> both ends found
219  else colourMap[col] = abs(colourMap[col]);
220  }
221
222  if (acol != 0) {
223  // First time a tag is encountered, mark it negative -> one end found
224  if (acolIndx == 0) {
225  // Ensure gluons have different col and acol indices
226  while (acolIndx == 0 || colIndx == acolIndx) {
227  acolIndx = nextTagBase + int(rndmPtr->flat()*9) + 1;
228  }
229  colourMap[acol] = -acolIndx;
230  }
231  // Second time mark it positive -> both ends found
232  else colourMap[acol] = abs(colourMap[acol]);
233  }
234  // Update nextTagBase
235  nextTagBase += 10;
236  }
237
238  // Check if these assignments would produce any singlet gluons
239  accept = true;
240  for (int i=0; i<partonSystemsPtr->sizeAll(iSys); ++i) {
241  int i1 = partonSystemsPtr->getAll(iSys,i);
242  Particle* partonPtr = &event[i1];
243  if (partonPtr->colType() != 2) continue;
244  int colIndexNew = colourMap[partonPtr->col()] % 10;
245  int acolIndexNew = colourMap[partonPtr->acol()] % 10;
246  if (colIndexNew == acolIndexNew) {
247  accept=false;
248  break;
249  }
250  }
251  }
252
253  // Check for failure to find acceptable conf.
254  if (!accept) {
255  if (verbose >= 3) event.list();
256  printOut("Colour::colourise","Warning! failed to avoid singlet gluon(s).");
257  }
258
259  // Update event.
260  for (int i = 0; i < partonSystemsPtr->sizeAll(iSys); ++i) {
261  int ip = partonSystemsPtr->getAll(iSys,i);
262  Particle* partonPtr = &event[ip];
263  if (partonPtr->colType() == 0) continue;
264  if ( colourMap[partonPtr->col()] > 0 )
265  partonPtr->col(colourMap[partonPtr->col()]);
266  if ( colourMap[partonPtr->acol()] > 0 )
267  partonPtr->acol(colourMap[partonPtr->acol()]);
268
269  // Update max used colour tag.
270  int lastTag = event.lastColTag();
271  int colMax = max(abs(partonPtr->col()),abs(partonPtr->acol()));
272  while (colMax > lastTag) lastTag = event.nextColTag();
273  }
274
275  // Return successful.
276  return true;
277 }
278
279 //--------------------------------------------------------------------------
280
281 // Order a list of partons in colour sequence.
282
283 vector<int> Colour::colourSort(vector<Particle*> partons) {
284
285  // Output vector (starts empty).
286  vector<int> iSorted;
287
288  // Shorthand for final-state parton multiplicities.
289  int nPartons=partons.size();
290  if (nPartons <= 1) return iSorted;
291
292  // Find string endpoints and colour types of particles.
293  vector<int> iTrip, iAnti, iOct, iOtherIn, iOtherOut;
294
295  // Definition of colType (classified by multiplet up to total charge
296  // p+q = 4).
297  // +- 1 : triplet (e.g., Red) [1,0] / [0,1] {quark, antidiquark}
298  // 2 : octet (e.g., R-Gbar) [1,1] {
299  // gluon, incoherent qqbar, partially coherent gluon-gluon
300  // eg R-(Bbar-B)-Gbar -> R-Gbar (no junction)
301  // or (R-B)-(Bbar-Gbar) -> Gbar-R (junction-antijunction)}
302  // +- 3 : sextet (e.g., 2 x Red) [2,0] / [0,2] {2 incoherent quarks}
303  // +- 4 : fifteen (e.g., R-R-Gbar) [2,1] / [1,2] {incoherent qg}
304  // +- 5 : decuplet (e.g., 3 x Red) [3,0] / [0,3] {
305  // 3 incoherent quarks / partially coherent gluon-gluon
306  // eg R-Gbar-R-Bbar -> R-R-R}
307  // 6 : vigintiseptet (e.g., 2 x R-Gbar) [2,2] {2 incoherent gluons}
308  // +- 7 : fifteen' (e.g., 4 x Red) [4,0] / [0,4] {4 incoherent quarks}
309  // +- 8 : vigintiquartet (e.g., R-R-R-Gbar) [3,1] / [1,3]
310
311  map<int, int> iOfAcol;
312  for (int i=partons.size()-1; i>=0; --i) {
313  int sign = (partons[i]->isFinal() ? 1 : -1);
314  int colType = particleDataPtr->colType(partons[i]->id());
315
316  // Store indices of anticolour partners.
317  if (sign == 1 && ( colType == -1 || colType == 2))
318  iOfAcol[partons[i]->acol()] = i;
319  else if (sign == -1 && ( colType == 1 || colType == 2 ))
320  iOfAcol[partons[i]->col()] = i;
321
322  // Construct list of triplets (= starting points).
323  if (colType * sign == 1) iTrip.push_back(i);
324
325  // List of antitriplets.
326  else if (colType * sign == -1) iAnti.push_back(i);
327  // List of octets.
328  else if (colType == 2) iOct.push_back(i);
329
330  // Higher representations.
331  else if (abs(colType) >= 3) {
332  cout << "colourSort(): ERROR! handling of coloured particles in "
333  << "representations higher than triplet or octet is not implemented"
334  << endl;
335  }
336
337  // Colourless particles.
338  else if (sign == -1) iOtherIn.push_back(i);
339  else iOtherOut.push_back(i);
340  }
341
342  // Now sort particles.
343  int i1 = -1;
344  bool beginNewChain = true;
345
346  // Keep looping until we have sorted all particles.
347  while (iSorted.size() < partons.size()) {
348
349  // Start new piece (also add colourless particles at front and end).
350  if (beginNewChain) {
351
352  // Insert any incoming colourless particles at front of iSorted.
353  if (iOtherIn.size() > 0) {
354  iSorted.push_back(iOtherIn.back());
355  iOtherIn.pop_back();
356
357  // Triplet starting point (charge += 1).
358  } else if (iTrip.size() > 0) {
359  beginNewChain = false;
360  iSorted.push_back(iTrip.back());
361  iTrip.pop_back();
362
363  // Octet starting point if no triplets/sextets available.
364  } else if (iOct.size() > 0) {
365  beginNewChain = false;
366  iSorted.push_back(iOct.back());
367  iOct.pop_back();
368
369  } else if (iOtherOut.size() > 0) {
370  iSorted.push_back(iOtherOut.back());
371  iOtherOut.pop_back();
372  }
373
374  // Index of current starting parton.
375  i1 = iSorted.back();
376
377  // Step to next parton in this chain.
378  } else {
379  bool isFinal = partons[iSorted.back()]->isFinal();
380  int col = (isFinal) ? partons[iSorted.back()]->col()
381  : partons[iSorted.back()]->acol();
382  int iNext = iOfAcol[col];
383
384  // Sanity check.
385  if (iNext < 0) {
386  cout << "colourSort(): ERROR! cannot step to < 0" << endl;
387  beginNewChain = true;
388
389  // Catch close of gluon ring.
390  } else if (iNext == i1) {
391  beginNewChain = true;
392
393  // Step to next parton; end if not gluon (antiquark or IS quark).
394  } else {
395  // Add to sorted list.
396  iSorted.push_back(iNext);
397  // If endpoint reached, begin new chain.
398  if (particleDataPtr->colType(partons[iNext]->id()) != 2)
399  beginNewChain = true;
400  // Octet: continue chain and erase this octet from list.
401  else {
402  beginNewChain = false;
403  // Erase this endpoint from list.
404  for (int i=0; i<(int)iOct.size(); ++i) {
405  if (iOct[i] == iNext) {
406  iOct.erase(iOct.begin()+i);
407  break;
408  }
409  }
410  }
411  }
412  } // End step to next parton.
413  }
414
415  // Normal return.
416  return iSorted;
417
418 }
419
420 //--------------------------------------------------------------------------
421
422 // Make colour maps and construct list of parton pairs that form QCD dipoles.
423
424 void Colour::makeColourMaps(const int iSysIn, const Event& event,
425  map<int,int>& indexOfAcol, map<int,int>& indexOfCol,
426  vector< pair<int,int> >& antLC, const bool findFF, const bool findIX) {
427
428  // Loop over all parton systems.
429  int iSysBeg = (iSysIn >= 0) ? iSysIn : 0;
430  int iSysEnd = (iSysIn >= 0) ? iSysIn + 1: partonSystemsPtr->sizeSys();
431  for (int iSys = iSysBeg; iSys < iSysEnd; ++iSys) {
432
433  // Loop over a single parton system.
434  int sizeSystem = partonSystemsPtr->sizeAll(iSys);
435  for (int i = 0; i < sizeSystem; ++i) {
436  int i1 = partonSystemsPtr->getAll( iSys, i);
437  if ( i1 <= 0 ) continue;
438
439  // Save to colour maps.
440  int col = event[i1].col();
441  int acol = event[i1].acol();
442
443  // Switch colours for initial partons.
444  if (!event[i1].isFinal()) {
445  col = acol;
446  acol = event[i1].col();
447  }
448
449  // Save colours (taking negative-index sextets into account).
450  if (col > 0) indexOfCol[col] = i1;
451  else if (col < 0) indexOfAcol[-col] = i1;
452  if (acol > 0) indexOfAcol[acol] = i1;
453  else if (acol < 0) indexOfCol[-acol] = i1;
454
455  // Look for partner on colour side.
456  if (col > 0 && indexOfAcol.count(col) == 1) {
457  int i2 = indexOfAcol[col];
458  if ( event[i1].isFinal() && event[i2].isFinal() ) {
459  if (findFF) antLC.push_back( make_pair(i1,i2) );
460  } else if (findIX) antLC.push_back( make_pair(i1,i2) );
461  }
462
463  // Look for partner on anticolour side.
464  if (acol > 0 && indexOfCol.count(acol) == 1) {
465  int i2 = indexOfCol[acol];
466  // Coloured parton first -> i2, i1 instead of i1, i2)
467  if (event[i1].isFinal() && event[i2].isFinal()) {
468  if (findFF) antLC.push_back( make_pair(i2, i1) );
469  } else if (findIX) antLC.push_back( make_pair(i2,i1) );
470  }
471
472  // Allow for sextets: negative acol -> extra positive col.
473  if (acol < 0 && indexOfAcol.count(-acol) == 1) {
474  int i2 = indexOfAcol[-acol];
475  if (event[i1].isFinal() && event[i2].isFinal()) {
476  if (findFF) antLC.push_back( make_pair(i1,i2) );
477  } else if (findIX) antLC.push_back( make_pair(i1,i2) );
478  }
479  if (col < 0 && indexOfCol.count(-col) == 1) {
480  int i2 = indexOfAcol[-acol];
481  if (event[i1].isFinal() && event[i2].isFinal()) {
482  if (findFF) antLC.push_back( make_pair(i1,i2) );
483  } else if (findIX) antLC.push_back( make_pair(i1,i2) );
484  }
485  }
486  }
487  return;
488
489 }
490
491 //--------------------------------------------------------------------------
492
493 // Determine which of two antennae inherits the old colour tag after a
494 // branching. Default is that the largest invariant has the largest
495 // probability to inherit, with a few alternatives also implemented.
496
497 bool Colour::inherit01(double s01, double s12) {
498  // Initialization check.
499  if (!isInit) {
500  printOut("VinciaColour::inherit01",
501  "ERROR! Colour not initialised");
502  if (isInitPtr && rndmPtr->flat() < 0.5) return false;
503  else return true;
504  }
505
506  // Mode 0: Purely random.
507  if (inheritMode == 0) {
508  if (rndmPtr->flat() < 0.5) return true;
509  else return false;
510  }
511
512  // Safety checks: small, or approximately equal s01, s12.
513  double a12 = abs(s01);
514  double a23 = abs(s12);
515
516  // Inverted mode (smallest invariant inherits - should only be used
517  // for extreme variation checks).
518  if (inheritMode < 0) {
519  a12 = abs(s12);
520  a23 = abs(s01);
521  inheritMode = abs(inheritMode);
522  }
523
524  // Winner-takes-all mode.
525  if (inheritMode == 2) {
526  if (a12 > a23) return true;
527  else return false;
528  }
529  double p12 = 0.5;
530  if ( max(a12,a23) > TINY ) {
531  if ( a12 < TINY ) p12 = 0.;
532  else if ( a23 < TINY ) p12 = 1.;
533  else {
534  double r = a23/a12;
535  if (r < TINY) p12 = 1. - r;
536  else if (r > 1./TINY) p12 = 1./r;
537  else p12 = 1./(1. + r);
538  }
539  }
540  if (rndmPtr->flat() < p12) return true;
541  else return false;
542
543 }
544
545 //==========================================================================
546
547 // The Resolution class.
548
549 //--------------------------------------------------------------------------
550
551 // Initialize.
552
553 bool Resolution::init() {
554
555  // Check that pointers initialized.
556  if (!isInitPtr) {
557  printOut("Resolution::init","Cannot initialize, pointers not set.");
558  return false;
559  }
560
561  // Set members.
562  verbose = settingsPtr->mode("Vincia:verbose");
563  nFlavZeroMassSav = settingsPtr->mode("Vincia:nFlavZeroMass");
564  isInit = true;
565  return isInit;
566
567 }
568
569 //--------------------------------------------------------------------------
570
571 // Sector resolution functions.
572
573 double Resolution::q2sector2to3(const Particle* a, const Particle* b,
574  const Particle* j, bool) {
575
576  // Construct basic 4-products.
577  double saj = 2*a->p()*j->p();
578  double sjb = 2*b->p()*j->p();
579  double sab = 2*a->p()*b->p();
580  // Gluon emission.
581  if (j->id() == 21) {
582  // FF emission.
583  if (a->isFinal() && b->isFinal()) {
584  return saj * sjb / (saj + sjb + sab) ;
585
586  // IF emission.
587  } else if (b->isFinal()) {
588  return saj * sjb / (saj + sab) ;
589
590  // IF emission (a <-> b).
591  } else if (a->isFinal()) {
592  return saj * sjb / (sjb + sab) ;
593
594  // II emission.
595  } else {
596  return saj * sjb / sab;
597  }
598
599  // FF Gluon splitting.
600  } else if (a->isFinal() && b->isFinal()) {
601  // Assume b is recoiler.
602  double m2j = pow2(j->m());
603  double m2qq = saj + 2*m2j;
604  // Find colour-connected invariant.
605  double sX = 0;
606  if (j->col() != 0 && j->col() == b->acol())
607  sX = sjb + m2j;
608  else {
609  sX = sab + m2j;
610  }
611  // Normalisation.
612  double sAnt = saj + sjb + sab + 2 * m2j;
613
614  // Return the g->qq sector variable defined in arXiv:1109.3608.
615  return m2qq * sqrt(sX/sAnt);
616  } else {
617  cout << "Sector criterion not implemented for II/IF splittings/conversions"
618  << endl;
619  }
620  return -1;
621
622 }
623
624 //--------------------------------------------------------------------------
625
626 // Sector resolution function for 3->4 branchings (currently only used
627 // for gluon splitting, with m2qq as the measure).
628
629 double Resolution::q2sector3to4(const Particle*, const Particle* ,
630  const Particle* j1, const Particle* j2) {
631  Vec4 pqq = j1->p() + j2->p();
632  double m2qq = pqq.m2Calc();
633  return m2qq;
634 }
635
636 //--------------------------------------------------------------------------
637
638 // Sector resolution function for 2->4 branchings (double emission).
639 // Assume j1 and j2 are colour connected, with a and b hard recoilers.
640
641 double Resolution::q2sector2to4(const Particle* a, const Particle* b,
642  const Particle* j1, const Particle* j2) {
643  return min(q2sector2to3(a,j2,j1),q2sector2to3(j1,b,j2));
644 }
645
646 //--------------------------------------------------------------------------
647
648 // Sector resolution function for 3->5 branchings
649 // (emission + splitting).
650
651 double Resolution::q2sector3to5(Particle* a, Particle* b,
652  Particle* j1, Particle* j2, Particle* j3) {
653
654  // j3 is gluon.
655  Particle* gluPtr;
656  Particle* qPtr;
657  Particle* qBarPtr;
658  if (j1->id() == 21) {
659  gluPtr = j1;
660  qPtr = (j2->id() > 0 ? j2 : j3);
661  qBarPtr = (j2->id() < 0 ? j2 : j3);
662  } else if (j2->id() == 21) {
663  gluPtr = j2;
664  qPtr = (j1->id() > 0 ? j1 : j3);
665  qBarPtr = (j1->id() < 0 ? j1 : j3);
666  } else if (j3->id() == 21) {
667  gluPtr = j3;
668  qPtr = (j2->id() > 0 ? j2 : j1);
669  qBarPtr = (j2->id() < 0 ? j2 : j1);
670  } else {
671  cout << " q2sector3to5: unable to identify branching type" << endl;
672  return 1.e19;
673  }
674  Vec4 pqq = qPtr->p() + qBarPtr->p();
675  double m2qq = pqq.m2Calc();
676  Particle* colPtr = a;
677  if (a->col() != gluPtr->acol()) colPtr = j1;
678  if (j1->col() != gluPtr->acol()) colPtr = j2;
679  if (j2->col() != gluPtr->acol()) colPtr = j3;
680  if (j3->col() != gluPtr->acol()) colPtr = b;
681  Particle* acolPtr = b;
682  if (b->acol() != gluPtr->col()) acolPtr = j3;
683  if (j3->acol() != gluPtr->col()) acolPtr = j2;
684  if (j2->acol() != gluPtr->col()) acolPtr = j1;
685  if (j1->acol() != gluPtr->col()) acolPtr = a;
686  double q2emit = q2sector2to3(colPtr,acolPtr,gluPtr);
687  return min(q2emit,m2qq);
688
689 }
690
691 //--------------------------------------------------------------------------
692
693 // Sector accept function. Optionally prevent g->qq clusterings if
694 // that would reduce the number of fermion lines below some minimum
695 // (cheap way to indicate that Z->qq/ll always has at least one
696 // fermion pair).
697
698 double Resolution::findSector(vector<int>& iSec, vector<Particle> state,
699  int nFmin) {
700
701  map<int, int> colMap, acolMap;
702  map<int, vector<int> > flavMap;
703  vector<int> helX;
704  int nFerm(0), nIn(0);
705  for (int i=0; i<(int)state.size(); ++i) {
706  colMap[state[i].isFinal() ? state[i].col() : state[i].acol()] = i;
707  acolMap[state[i].isFinal() ? state[i].acol() : state[i].col()] = i;
708  flavMap[state[i].isFinal() ? state[i].id() :
709  -state[i].id()].push_back(i);
710  helX.push_back(state[i].isFinal() ? state[i].pol() : -state[i].pol());
711  if (state[i].isQuark() || state[i].isLepton()) ++nFerm;
712  if (!state[i].isFinal()) ++nIn;
713  }
714  nFerm /= 2;
715
716  // Initialise.
717  double q2min = 1e19;
718  iSec.clear();
719
720  // Do all possible 2->3 clusterings. Note, needs modification for
721  // 3->4 and 2->4 showers.
722  for (int ir = 0; ir < (int)state.size(); ++ir) {
723  // There are no sectors for emission into the initial state.
724  if (!state[ir].isFinal()) continue;
725  // If a gluon, compute LC pT from colour partners.
726  if (state[ir].isGluon()) {
727  int ia = colMap[state[ir].acol()];
728  int ib = acolMap[state[ir].col()];
729  double q2this = q2sector2to3(&state[ia],&state[ib],&state[ir]);
730  if (q2this < q2min) {
731  q2min = q2this;
732  iSec.clear();
733  iSec.push_back(ia);
734  iSec.push_back(ib);
735  iSec.push_back(ir);
736  }
737  }
738
739  // If fermion, check if we are allowed to look at this clustering.
740  if (state[ir].isQuark()) {
741  // Skip if this clustering would reduce below number requested.
742  if (nFerm <= nFmin) continue;
743  // Cluster quark: recoiler is the anticolour partner.
744  // Cluster antiquark: recoiler is colour partner.
745  int ib = (state[ir].id() > 0) ? acolMap[state[ir].col()]
746  : colMap[state[ir].acol()];
747
748  // Loop over all same-flavour (anti)quarks (must in principle
749  // also require opposite helicities).
750  vector<int> aList;
751  aList = flavMap[-state[ir].id()];
752  for (int ifa = 0; ifa < (int)aList.size(); ++ifa) {
753  int ia = aList[ifa];
754  if (helX[ir] != helX[ia] || helX[ia] == 9 || helX[ir] == 9) {
755  double q2this = q2sector2to3(&state[ia],&state[ib],&state[ir]);
756  if (q2this < q2min) {
757  q2min = q2this;
758  iSec.clear();
759  iSec.push_back(ia);
760  iSec.push_back(ib);
761  iSec.push_back(ir);
762  }
763  }
764  }
765  }
766  // Other sectors (QED, EW, ...) to be implemented here.
767  } // End loop over ir.
768
769  // Return smallest scale found.
770  return q2min;
771
772 }
773
774 //==========================================================================
775
776 // The VinciaCommon class.
777
778 //--------------------------------------------------------------------------
779
780 // Initialize the class.
781
782 bool VinciaCommon::init() {
783
784  // Check initPtr.
785  if (!isInitPtr) {
786  if (verbose >= 1)
787  printOut("VinciaCommon::init","Error! pointers not initialized");
788  return false;
789  }
790
791  // Verbosity level and checks.
792  verbose = settingsPtr->mode("Vincia:verbose");
793  epTolErr = settingsPtr->parm("Check:epTolErr");
794  epTolWarn = settingsPtr->parm("Check:epTolWarn");
795  mTolErr = settingsPtr->parm("Check:mTolErr");
796  mTolWarn = settingsPtr->parm("Check:mTolWarn");
797
798  // Counters
799  nUnkownPDG = 0;
800  nIncorrectCol = 0;
801  nNAN = 0;
802  nVertex = 0;
803  nChargeCons = 0;
804  nMotDau = 0;
805  nUnmatchedMass.resize(2);
806  nEPcons.resize(2);
807  for (int i=0; i<2; i++) {
808  nUnmatchedMass[i] = 0;
809  nEPcons[i] = 0;
810  }
811
812  // Quark masses
813  mt = particleDataPtr->m0(6);
814  if (mt < TINY) mt = 171.0;
815  mb = min(mt,particleDataPtr->m0(5));
816  if (mb < TINY) mb = min(mt,4.8);
817  mc = min(mb,particleDataPtr->m0(4));
818  if (mc < TINY) mc = min(mb,1.5);
819  ms = min(mc,particleDataPtr->m0(3));
820  if (ms < TINY) ms = min(mc,0.1);
821
822  // Number of flavours to treat as massless in clustering and
823  // kinematics maps.
824  nFlavZeroMass = settingsPtr->mode("Vincia:nFlavZeroMass");
825
826  // Default alphaS, with and without CMW.
827  double alphaSvalue = settingsPtr->parmDefault("Vincia:alphaSvalue");
828  int alphaSorder = settingsPtr->modeDefault("Vincia:alphaSorder");
829  int alphaSnfmax = settingsPtr->modeDefault("Vincia:alphaSnfmax");
830  bool useCMW = settingsPtr->flagDefault("Vincia:useCMW");
831  alphaStrongDef.init( alphaSvalue, alphaSorder, alphaSnfmax, false);
832  alphaStrongDefCMW.init( alphaSvalue, alphaSorder, alphaSnfmax, true);
833
834  // Strong coupling for use in merging.
835  alphaSvalue = settingsPtr->parm("Vincia:alphaSvalue");
836  alphaSorder = settingsPtr->mode("Vincia:alphaSorder");
837  alphaSnfmax = settingsPtr->mode("Vincia:alphaSnfmax");
838  useCMW = settingsPtr->flag("Vincia:useCMW");
839  alphaS.init(alphaSvalue, alphaSorder, alphaSnfmax, useCMW);
840
841  // User alphaS, with and without CMW.
842  alphaSvalue = settingsPtr->parm("Vincia:alphaSvalue");
843  alphaSorder = settingsPtr->mode("Vincia:alphaSorder");
844  alphaSnfmax = settingsPtr->mode("Vincia:alphaSnfmax");
845  useCMW = settingsPtr->flag("Vincia:useCMW");
846  alphaStrong.init( alphaSvalue, alphaSorder, alphaSnfmax, false);
847  alphaStrongCMW.init( alphaSvalue, alphaSorder, alphaSnfmax, true);
848
849  // Freeze and minimum scales.
850  mu2freeze = pow2(settingsPtr->parm("Vincia:alphaSmuFreeze"));
851  alphaSmax = settingsPtr->parm("Vincia:alphaSmax");
852
853  // Find the overall minimum scale. Take into account the freezeout
854  // scale, Lambda pole, and alphaSmax.
855  double muMin = max(sqrt(mu2freeze),1.05*alphaS.Lambda3());
856  double muMinASmax;
857  if (alphaStrong.alphaS(mu2min) < alphaSmax) {
858  muMinASmax = muMin;
859  } else if (settingsPtr->mode("Vincia:alphaSorder") == 0) {
860  muMinASmax = muMin;
861  } else {
862  muMinASmax = muMin;
863  while (true) {
864  if (alphaS.alphaS(pow2(muMinASmax)) < alphaSmax) break;
865  muMinASmax += 0.001;
866  }
867  }
868  mu2min = pow2( max(muMinASmax, muMin) );
869
870  // EM coupling for use in merging. Dummy, as no EW clusterings.
871  alphaEM.init(1, settingsPtr);
872
873  // Return.
874  isInit = true;
875  return true;
876
877 }
878
879 //--------------------------------------------------------------------------
880
881 // Function to find the lowest meson mass for a parton pair treating
882 // gluons as down quarks. Used to determine hadronisation boundaries
883 // consistent with assumption of string length > 1 meson.
884
885 double VinciaCommon::mHadMin(const int id1in, const int id2in) {
886  int id1 = abs(id1in);
887  if (id1 == 21 || id1 <= 2) id1 = 1;
888  int id2 = abs(id2in);
889  if (id2 == 21 || id1 <= 2) id2 = 1;
890  int idMes = max(id1,id2)*100 + min(id1,id2)*10 + 1;
891
892  // Special for ssbar, use eta rather than eta'.
893  if (idMes == 331) idMes = 221;
894  return particleDataPtr->m0(idMes);
895
896 }
897
898 //--------------------------------------------------------------------------
899
900 // Function to check the event after each branching, mostly copied
901 // from Pythia8.
902
903 bool VinciaCommon::showerChecks(Event& event, bool ISR) {
904
905  // Only for verbose >= 4.
906  if (verbose < 4) return true;
907
908  // Count incoming partons (beam daughters) with negative momentum and charge.
909  Vec4 pSum;
910  double chargeSum = 0.0;
912  for (int i = 0; i < event.size(); ++i){
914  if (!event[i].isFinal()) {
915  if ( (event[i].mother1() == 1) || (event[i].mother1() == 2) ) {
916  pSum -= event[i].p();
917  chargeSum -= event[i].charge();
918  }
919  }
920  if(abs(event[i].status())==63){
922  pSum = -(event[1].p() + event[2].p());
923  chargeSum = -(event[1].charge()+event[2].charge());
924  }
925  }
926  }
927  double eLab = abs(pSum.e());
928
929  // Loop over particles in the event.
930  for (int i = 0; i < event.size(); ++i) {
931
932  // Look for any unrecognized particle codes.
933  int id = event[i].id();
934  if (id == 0 || !particleDataPtr->isParticle(id)) {
935  nUnkownPDG++;
936  if (nUnkownPDG == 1){
937  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
938  << ": unknown particle code"
939  << ", i = " << i << ", id = " << id << endl;
940  return false;
941  }
942  }
943
944  // Check that colour assignments are the expected ones.
945  else {
946  int colType = event[i].colType();
947  int col = event[i].col();
948  int acol = event[i].acol();
949  if ( (colType == 0 && (col > 0 || acol > 0))
950  || (colType == 1 && (col <= 0 || acol > 0))
951  || (colType == -1 && (col > 0 || acol <= 0))
952  || (colType == 2 && (col <= 0 || acol <= 0)) ) {
953  nIncorrectCol++;
954  if (nIncorrectCol == 1){
955  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
956  << ": incorrect colours"
957  << ", i = " << i << ", id = " << id << " cols = " << col
958  << " " << acol << endl;
959  return false;
960  }
961  }
962  }
963
964  // Look for particles with mismatched or not-a-number energy/momentum/mass.
965  if (abs(event[i].px()) >= 0.0 && abs(event[i].py()) >= 0.0
966  && abs(event[i].pz()) >= 0.0 && abs(event[i].e()) >= 0.0
967  && abs(event[i].m()) >= 0.0 ) {
968  double errMass = abs(event[i].mCalc() - event[i].m()) /
969  max( 1.0, event[i].e());
970
971  if (errMass > mTolErr) {
972  nUnmatchedMass[0]++;
973  if (nUnmatchedMass[0] == 1){
974  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
975  << ": unmatched particle energy/momentum/mass"
976  << ", i = " << i << ", id = " << id << endl;
977  return false;
978  }
979  } else if (errMass > mTolWarn) {
980  nUnmatchedMass[1]++;
981  if (nUnmatchedMass[1] == 1){
982  cout << "WARNING in Vincia::ShowerChecks"
983  << (ISR ? "(ISR)" : "(FSR)")
984  << ": not quite matched particle energy/momentum/mass"
985  << ", i = " << i << ", id = " << id << endl;
986  }
987  }
988  } else {
989  nNAN++;
990  if (nNAN == 1){
991  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
992  << ": not-a-number energy/momentum/mass"
993  << ", i = " << i << ", id = " << id << endl;
994  return false;
995  }
996  }
997
998  // Look for particles with not-a-number vertex/lifetime.
999  if (abs(event[i].xProd()) >= 0.0 && abs(event[i].yProd()) >= 0.0
1000  && abs(event[i].zProd()) >= 0.0 && abs(event[i].tProd()) >= 0.0
1001  && abs(event[i].tau()) >= 0.0) {
1002  } else {
1003  nVertex++;
1004  if (nVertex == 1){
1005  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
1007  << ", i = " << i << ", id = " << id << endl;
1008  return false;
1009  }
1010  }
1011
1012  // Add final-state four-momentum and charge.
1013  if (event[i].isFinal()) {
1014  pSum += event[i].p();
1015  chargeSum += event[i].charge();
1016  }
1017  } // End of particle loop.
1018
1019  // Check energy-momentum/charge conservation.
1020  double epDev = abs( pSum.e()) + abs(pSum.px()) + abs(pSum.py())
1021  + abs(pSum.pz() );
1022  if (epDev > epTolErr * eLab) {
1023  nEPcons[0]++;
1024  if (nEPcons[0] == 1) {
1025  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
1026  << ": energy-momentum not conserved" << endl;
1027  cout << " epDev = " << epDev << " Ein = " << eLab
1028  << " pSum = " << pSum << endl;
1029  return false;
1030  }
1031  } else if (epDev > epTolWarn * eLab) {
1032  nEPcons[1]++;
1033  if (nEPcons[1] == 1)
1034  cout << "WARNING in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
1035  << ": energy-momentum not quite conserved" << endl;
1036  }
1037  if (abs(chargeSum) > 0.1) {
1038  nChargeCons++;
1039  if (nChargeCons == 1){
1040  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
1041  << ": charge not conserved" << endl;
1042  return false;
1043  }
1044  }
1045
1046  // Check that mother and daughter information match for each particle.
1047  vector<int> noMot, noDau;
1048  vector< pair<int,int> > noMotDau;
1049
1050  // Loop through the event and check that there are beam particles.
1051  bool hasBeams = false;
1052  for (int i = 0; i < event.size(); ++i) {
1053  int status = event[i].status();
1054  if (abs(status) == 12) hasBeams = true;
1055
1056  // Check that mother and daughter lists not empty where not expected to.
1057  vector<int> mList = event[i].motherList();
1058  vector<int> dList = event[i].daughterList();
1059  if (mList.size() == 0 && abs(status) != 11 && abs(status) != 12)
1060  noMot.push_back(i);
1061  if (dList.size() == 0 && status < 0 && status != -11)
1062  noDau.push_back(i);
1063
1064  // Check that the particle appears in the daughters list of each mother.
1065  for (int j = 0; j < int(mList.size()); ++j) {
1066  if (event[mList[j]].daughter1() <= i
1067  && event[mList[j]].daughter2() >= i) continue;
1068  vector<int> dmList = event[mList[j]].daughterList();
1069  bool foundMatch = false;
1070  for (int k = 0; k < int(dmList.size()); ++k)
1071  if (dmList[k] == i) {
1072  foundMatch = true;
1073  break;
1074  }
1075  if (!hasBeams && mList.size() == 1 && mList[0] == 0) foundMatch = true;
1076  if (!foundMatch) {
1077  bool oldPair = false;
1078  for (int k = 0; k < int(noMotDau.size()); ++k)
1079  if (noMotDau[k].first == mList[j] && noMotDau[k].second == i) {
1080  oldPair = true;
1081  break;
1082  }
1083  if (!oldPair) noMotDau.push_back( make_pair( mList[j], i) );
1084  }
1085  }
1086
1087  // Check that the particle appears in the mothers list of each daughter.
1088  for (int j = 0; j < int(dList.size()); ++j) {
1089  if (event[dList[j]].statusAbs() > 80
1090  && event[dList[j]].statusAbs() < 90
1091  && event[dList[j]].mother1() <= i
1092  && event[dList[j]].mother2() >= i ) continue;
1093  vector<int> mdList = event[dList[j]].motherList();
1094  bool foundMatch = false;
1095  for (int k = 0; k < int(mdList.size()); ++k)
1096  if (mdList[k] == i) {
1097  foundMatch = true;
1098  break;
1099  }
1100  if (!foundMatch) {
1101  bool oldPair = false;
1102  for (int k = 0; k < int(noMotDau.size()); ++k)
1103  if (noMotDau[k].first == i && noMotDau[k].second == dList[j]) {
1104  oldPair = true;
1105  break;
1106  }
1107  if (!oldPair) noMotDau.push_back( make_pair( i, dList[j]) );
1108  }
1109  } // End loop through the event.
1110  }
1111
1112  // Warn if any errors were found.
1113  if (noMot.size() > 0 || noDau.size() > 0 || noMotDau.size() > 0) {
1114  nMotDau++;
1115  if (nMotDau == 1) {
1116  cout << "ERROR in Vincia::ShowerChecks" << (ISR ? "(ISR)" : "(FSR)")
1117  << ": mismatch in daughter and mother lists" << endl;
1119  if (noMot.size() > 0) {
1120  cout << " missing mothers for particles ";
1121  for (int i = 0; i < int(noMot.size()); ++i) cout << noMot[i] << " ";
1122  cout << endl;
1123  }
1124  if (noDau.size() > 0) {
1125  cout << " missing daughters for particles ";
1126  for (int i = 0; i < int(noDau.size()); ++i) cout << noDau[i] << " ";
1127  cout << endl;
1128  }
1129  if (noMotDau.size() > 0) {
1130  cout << " inconsistent history for (mother,daughter) pairs ";
1131  for (int i = 0; i < int(noMotDau.size()); ++i)
1132  cout << "(" << noMotDau[i].first << ","
1133  << noMotDau[i].second << ") ";
1134  cout << endl;
1135  }
1136  return false;
1137  }
1138  }
1139
1140  //Made it to here: no major problems.
1141  return true;
1142
1143 }
1144
1145 //--------------------------------------------------------------------------
1146
1147 // Get the shower starting scale.
1148
1149 double VinciaCommon::getShowerStartingScale(int iSys,
1150  PartonSystems* partonSystemsPtr, const Event& event, double sbbSav) {
1151
1152  // Depending on user choice shower starts at q2maxFudge *
1153  // factorization scale of phase space maximum.
1154  int qMaxMatch = settingsPtr->mode("Vincia:QmaxMatch");
1155  double q2maxFudge = pow2(settingsPtr->parm("Vincia:QmaxFudge"));
1156  bool hasFSpartons = false;
1157  int nOut = partonSystemsPtr->sizeOut(iSys);
1158  vector<int> iFS;
1159  for (int iOut = 0; iOut < nOut; ++iOut) {
1160  int i = partonSystemsPtr->getOut(iSys,iOut);
1161  int idAbs = event[i].idAbs();
1162  if ((idAbs >= 1 && idAbs <= 5) || idAbs == 21) hasFSpartons = true;
1163  iFS.push_back(i);
1164  }
1165  if (qMaxMatch == 1 || (qMaxMatch == 0 && hasFSpartons) ) {
1166  double Q2facSav = sbbSav;
1167  double Q2facFix = settingsPtr->parm("SigmaProcess:factorFixScale");
1168  double Q2facMult = settingsPtr->parm("SigmaProcess:factorMultFac");
1169
1171  if (nOut == 1) {
1172  int factorScale1 = settingsPtr->mode("SigmaProcess:factorScale1");
1173  Q2facSav = ( factorScale1 == 1 ? Q2facMult*event[iFS[0]].m2Calc() :
1174  Q2facFix );
1175
1177  } else if (iFS.size() == 2) {
1178  int factorScale2 = settingsPtr->mode("SigmaProcess:factorScale2");
1179  double mT21 = event[iFS[0]].mT2(), mT22 = event[iFS[1]].mT2();
1180  double sHat = m2(event[iFS[0]],event[iFS[1]]);
1181  double tHat = m2(event[3],event[iFS[0]]);
1182  if (factorScale2 == 1) Q2facSav = min(mT21,mT22);
1183  else if (factorScale2 == 2) Q2facSav = sqrt(mT21*mT22);
1184  else if (factorScale2 == 3) Q2facSav = 0.5*(mT21+mT22);
1185  else if (factorScale2 == 4) Q2facSav = sHat;
1186  else if (factorScale2 == 5) Q2facSav = Q2facFix;
1187  else if (factorScale2 == 6) Q2facSav = abs(-tHat);
1188  if (factorScale2 != 5) Q2facSav *= Q2facMult;
1189
1191  } else if (iFS.size() == 3) {
1192  int factorScale3 = settingsPtr->mode("SigmaProcess:factorScale3");
1193  double mT21 = event[iFS[0]].mT2(), mT22 = event[iFS[1]].mT2(),
1194  mT23 = event[iFS[2]].mT2();
1195  double mT2min1 = min(mT21,min(mT22,mT23));
1196  double mT2min2 = max(max(min(mT21,mT22),min(mT21,mT23)),min(mT22,mT23));
1197  double sHat = m2(event[iFS[0]],event[iFS[1]],event[iFS[2]]);
1198  if (factorScale3 == 1) Q2facSav = mT2min1;
1199  else if (factorScale3 == 2) Q2facSav = sqrt(mT2min1*mT2min2);
1200  else if (factorScale3 == 3) Q2facSav = pow(mT21*mT22*mT23,1.0/3.0);
1201  else if (factorScale3 == 4) Q2facSav = (mT21+mT22+mT23)/3.0;
1202  else if (factorScale3 == 5) Q2facSav = sHat;
1203  else if (factorScale3 == 6) Q2facSav = Q2facFix;
1204  if (factorScale3 != 6) Q2facSav *= Q2facMult;
1205
1206  // Unknown, leave as is, all emissions allowed now.
1207  } else {;}
1208  return q2maxFudge*Q2facSav;
1209  }
1210  return sbbSav;
1211
1212 }
1213
1214 //--------------------------------------------------------------------------
1215
1216 // FF clustering map(s) for massless partons. Inputs are as follows:
1217 // kMapType = map number ( 1 : ARIADNE, 2 : PS, 3 : Kosower)
1218 // pIn = Vec4 list (arbitrary length, but greater than 3)
1219 // a,r,b = indices of 3 particles in pIn to be clustered
1220 // Outputs are as follows:
1221 // pClu = pIn but with the a and b momenta replaced by the clustered
1222 // aHat and bHat and with r erased.
1223 // For example:
1224 // pIn(...,a,...,r-1,r,r+1,...,b,...) ->
1225 // pOut(...,aHat,...,r-1,r+1,...,bHat,...)
1226 // with {a,r,b} -> {aHat,bHat} using kinematics map kMapType.
1227
1228 bool VinciaCommon::map3to2FFmassless(vector<Vec4>& pClu, vector<Vec4> pIn,
1229  int kMapType, int a, int r, int b) {
1230
1231  // Intialize and sanity check.
1232  pClu=pIn;
1233  if (max(max(a,r),b) > int(pIn.size()) || min(min(a,r),b) < 0) {
1234  if (verbose >= 3)
1235  printOut("VinciaCommon::map3to2FFmassless",
1236  "Error! Unable to cluster (a,r,b) = "+
1237  num2str(a)+num2str(r)+num2str(b)+" p.size ="
1238  +num2str(int(pIn.size())));
1239  return false;
1240  }
1241
1242  // Verbose output.
1243  if (verbose >= superdebug) {
1244  printOut("VinciaCommon:map3to2FFmassless", "called with ");
1245  cout << "pi = " << pIn[a];
1246  cout << "pj = " << pIn[r];
1247  cout << "pk = " << pIn[b];
1248  }
1249
1250  // Compute total invariant mass squared.
1251  Vec4 pSum = pIn[a] + pIn[r] + pIn[b];
1252  double m2Ant = pSum.m2Calc();
1253  if (m2Ant < 1e-20) {
1254  printOut("VinciaCommon::map3to2FFmassless",
1255  "Massless or spacelike system. Cannot find rest frame");
1256  return false;
1257  }
1258
1259  // ARIADNE and PS maps (recoded from old F77 maps for v.1.2.01)
1260  // (maps -1 and -2 are special: force A or B to take all recoil).
1261  if (kMapType == 1 || kMapType == 2 || kMapType == -1 || kMapType == -2) {
1262
1263  // Make copies of PA, PB, and compute sum of LAB momenta and CM energy.
1265  Vec4 pbDum = pIn[b];
1266  double eCM = sqrt(m2Ant);
1268  pbDum.bstback(pSum);
1269
1270  // Rotate so a goes into upper half of (x,z) plane.
1273  pbDum.rot(0.,-phiA);
1274
1275  // Rotate so a goes onto z axis.
1277  pbDum.rot(-theta,0.);
1278
1279  // Rotate so (r,b) go into (x,z) plane.
1280  double phiB = pbDum.phi();
1281
1282  // Compute psi(a,ahat) from A, B energies and theta(AB).
1283  double thetaAB = pbDum.theta();
1284  double psi = 0.0;
1285
1286  // ARIADNE angle (smoothly varying antenna recoil).
1287  if (kMapType == 1) {
1289  * (M_PI - thetaAB);
1290
1291  // PS angle (direction a fixed if sar > srb, and vice versa). Use
1292  // org momenta, since new ones may not all have been rotated.
1293  } else if (kMapType == 2) {
1294  Vec4 pAR = pIn[a] + pIn[r];
1295  Vec4 pRB = pIn[r] + pIn[b];
1296  double sar = pAR.m2Calc();
1297  double srb = pRB.m2Calc();
1298  if (sar > srb) psi = 0.0;
1299  else psi = (M_PI - thetaAB);
1300
1301  // Force A to be the emitter. B recoils longitudinally.
1302  } else if (kMapType == -1) {
1303  psi = M_PI - thetaAB;
1304
1305  // Force B to be the emitter. A recoils longitudinally.
1306  } else if (kMapType == -2) {
1307  psi = 0.0;
1308  }
1309  // Now we know everything:
1310  // CM -> LAB : -PSI, PHIB, THETA, PHIA, BOOST
1311
1312  // Set up initial massless AHAT, BHAT with AHAT along z axis.
1313  pClu[a] = Vec4(0.,0.,eCM/2.,eCM/2.);
1314  pClu[b] = Vec4(0.,0.,-eCM/2.,eCM/2.);
1315
1316  // 1) Take into account antenna recoil, and rotate back in phiB.
1317  pClu[a].rot(-psi,phiB);
1318  pClu[b].rot(-psi,phiB);
1319
1320  // 2) Rotate back in theta and phiA.
1321  pClu[a].rot(theta,phiA);
1322  pClu[b].rot(theta,phiA);
1323
1324  // 3) Boost back to LAB.
1325  pClu[a].bst(pSum);
1326  pClu[b].bst(pSum);
1327
1328  // kMapType = 3, 4 (and catchall for undefined kMapType
1329  // values). Implementation of the massless Kosower antenna map(s).
1330  } else {
1331
1332  // Compute invariants.
1333  double s01 = 2*pIn[a]*pIn[r];
1334  double s12 = 2*pIn[r]*pIn[b];
1335  double s02 = 2*pIn[a]*pIn[b];
1336
1337  // Check whether the arguments need to be reversed for kMapType == 4.
1338  if (kMapType == 4 && ( ! (s01 < s12) ) ) {
1339  if (verbose >= superdebug) {
1340  printOut("VinciaCommon::map3to2FFmassless",
1341  "choose parton i as the recoiler");
1342  }
1343  // Call the function with reverse arguments, then return.
1344  return map3to2FFmassless(pClu, pIn, kMapType, b, r, a);
1345  }
1346  double sAnt = s01 + s12 + s02;
1347
1348  // Compute coefficients.
1349  // kMapType = 3: GGG choice
1350  // kMapType >= 4: r=1
1351  double rMap = kMapType == 3 ? s12/(s01 + s12) : 1;
1352  double rho = sqrt(1.0+(4*rMap*(1.0-rMap)*s01*s12)/sAnt/s02);
1353  double x = 0.5/(s01+s02)*((1+rho)*(s01+s02)+(1+rho-2*rMap)*s12);
1354  double z = 0.5/(s12+s02)*((1-rho)*sAnt-2*rMap*s01);
1355
1356  // Compute reclustered momenta.
1357  pClu[a] = x*pIn[a] + rMap*pIn[r] + z*pIn[b];
1358  pClu[b] = (1-x)*pIn[a] + (1-rMap)*pIn[r] + (1-z)*pIn[b];
1359  }
1360
1361  // A dimensionless quantitiy to compare with TINY.
1362  if (pClu[a].m2Calc()/m2Ant >= TINY || pClu[b].m2Calc()/m2Ant >= TINY) {
1363  if (verbose >= 3)
1364  printOut("VinciaCommon::map3to2FFmassless",
1365  "on-shell check failed. m2I/sIK ="
1366  + num2str(pClu[a].m2Calc()/m2Ant)+" m2K/m2Ant ="
1367  + num2str(pClu[b].m2Calc()/m2Ant)+" m2Ant = "+num2str(m2Ant));
1368  return false;
1369  }
1370
1371  // Erase the clustered momentum and return.
1372  pClu.erase(pClu.begin()+r);
1373  return true;
1374
1375 }
1376
1377 //--------------------------------------------------------------------------
1378
1379 // Implementations of FF clustering maps for massive partons. See
1380 // VinciaCommon::map3to2FFmassless for details but with the additional
1381 // input of masses mI and mK for the first and second parent
1382 // particles.
1383
1384 bool VinciaCommon::map3to2FFmassive(vector<Vec4>& pClu, vector<Vec4> pIn,
1385  int kMapType, int a, int r, int b, double mI, double mK) {
1386
1387  // If both parent masses are negligible and all the momenta are
1388  // measure off-shellness normalised to average energy of the partons
1389  // to be clustered, avoids small denominator for collinear massless
1390  // p_a, p_r, p_b.
1391  double eAv = 1.0/3.0*( pIn[a].e() + pIn[r].e() + pIn[b].e() );
1392  if (mI/eAv < TINY && mK/eAv < TINY && pIn[a].mCalc()/eAv < TINY
1393  && pIn[r].mCalc()/eAv < TINY && pIn[b].mCalc()/eAv < TINY ) {
1394  return map3to2FFmassless(pClu, pIn, kMapType, a, r, b);
1395  }
1396
1397  // Intialize and sanity check.
1398  pClu = pIn;
1399  if (max(max(a,r),b) > int(pIn.size()) || min(min(a,r),b) < 0) return false;
1400
1401  // Verbose output.
1402  if (verbose >= superdebug) {
1403  printOut("VinciaCommon:map3to2FFmassive", "called with ");
1404  cout << "p0 = " << pIn[a];
1405  cout << "p1 = " << pIn[r];
1406  cout << "p2 = " << pIn[b];
1407  }
1408
1409  // ARIADNE map not defined for massive partons; use Kosower map instead.
1410  if (kMapType == 1) kMapType = 3;
1411
1412  // Longitudinal map; use Kosower map with r = 1.
1413  if (kMapType == 2) kMapType = 4;
1414
1415  // Forced-longitudinal maps.
1416  if (kMapType < 0) {
1417  printOut("VinciaCommon::map3to2FFmassive", "longitudinal clustering maps "
1418  "not yet implemented for massive partons.");
1419  return false;
1420  }
1421
1422  // Compute invariants.
1423  double m0 = pIn[a].mCalc();
1424  double m1 = pIn[r].mCalc();
1425  double m2 = pIn[b].mCalc();
1426  double s01 = 2*pIn[a]*pIn[r];
1427  double s12 = 2*pIn[r]*pIn[b];
1428  double s02 = 2*pIn[a]*pIn[b];
1429
1430  // Check whether the arguments need to be reversed for mapType == 4.
1431  if (kMapType == 4 && (! (s01 < s12) ) ) {
1432  return map3to2FFmassive(pClu, pIn, kMapType, b, r, a, mK, mI);
1433  }
1434  Vec4 pAnt = pIn[a] + pIn[r] + pIn[b];
1435  double m2Ant = pAnt.m2Calc();
1436  double mAnt = sqrt(m2Ant);
1437
1438  // Rewrite the determination in terms of dimensionless variables.
1439  // Note normalisation choice here is mAnt, rather than sAnt.
1440  double mu0 = m0/mAnt;
1441  double mu1 = m1/mAnt;
1442  double mu2 = m2/mAnt;
1443  double y01 = s01/m2Ant;
1444  double y12 = s12/m2Ant;
1445  double y02 = s02/m2Ant;
1446  double y01min = 2*mu0*mu1;
1447  double y12min = 2*mu1*mu2;
1448  double y02min = 2*mu0*mu2;
1449  double muI = mI/mAnt;
1450  double muK = mK/mAnt;
1451  double yIK = 1. - pow2(muI) - pow2(muK);
1452  double yIKmin = 2*muI*muK;
1453  double sig2 = 1 + pow2(muI) - pow2(muK);
1454  double gdetdimless = gramDet(y01,y12,y02,mu0,mu1,mu2);
1455  double gdetdimless01 = (y02*y12-2.*pow2(mu2)*y01)/4.;
1456  double gdetdimless12 = (y02*y01-2.*pow2(mu0)*y12)/4.;
1457  double rMap = 1.;
1458  if ( kMapType == 3) {
1459  rMap =
1460  (
1461  sig2
1462  + sqrt( pow2(yIK) - pow2(yIKmin) )
1463  *( y12 - y12min - ( y01 - y01min ) )
1464  /( y12 - y12min + ( y01 - y01min ) )
1465  )/2.;
1466
1467  // Fallback: map with massless r -> 1.
1468  } else {
1469  double mu01squa = pow2(mu0) + pow2(mu1) + y01;
1470  double lambda = 1 + pow2(mu01squa) + pow4(mu2) - 2*mu01squa - 2*pow2(mu2)
1471  - 2*mu01squa*pow2(mu2);
1472  rMap = (sig2 + sqrt(pow2(yIK) - pow2(yIKmin))
1473  *(1 - pow2(mu0) - pow2(mu1) + pow2(mu2) - y01)/sqrt(lambda))/2.;
1474  }
1475
1476  // Compute reclustered momenta.
1477  double bigsqr = sqrt(
1478  16.*( rMap*(1.-rMap) - (1.-rMap)*pow2(muI) - rMap*pow2(muK) )*gdetdimless
1479  + ( pow2(y02) - pow2(y02min) )*( pow2(yIK) - pow2(yIKmin) ));
1480  double x = (
1481  sig2*( pow2(y02) - pow2(y02min) + 4.*gdetdimless01)
1482  + 8.*rMap*( gdetdimless - gdetdimless01 )
1483  + bigsqr*( 1. - pow2(mu0) - pow2(mu1) + pow2(mu2) - y01)
1484  )/(2.*( 4.*gdetdimless + pow2(y02) - pow2(y02min) ));
1485  double z = (
1486  sig2*( pow2(y02) - pow2(y02min) + 4.*gdetdimless12)
1487  + 8.*rMap*( gdetdimless - gdetdimless12 )
1488  - bigsqr*( 1. + pow2(mu0) - pow2(mu1) - pow2(mu2) - y12)
1489  )/(2.*( 4.*gdetdimless + pow2(y02) - pow2(y02min)));
1490  pClu[a] = x*pIn[a] + rMap*pIn[r] + z*pIn[b];
1491  pClu[b] = (1-x)*pIn[a] + (1-rMap)*pIn[r] + (1-z)*pIn[b];
1492
1493  // Check if on-shell.
1494  double offshellnessI = abs(pClu[a].m2Calc() - pow2(mI))/m2Ant;
1495  double offshellnessK = abs(pClu[b].m2Calc() - pow2(mK))/m2Ant;
1496  if (offshellnessI > TINY || offshellnessK > TINY) {
1497  if (verbose >= 3) {
1498  printOut("VinciaCommon::map3to2FFmassive","on-shell check failed");
1499  }
1500  return false;
1501  }
1502
1503  // Erase the clustered parton and return.
1504  pClu.erase(pClu.begin()+r);
1505  return true;
1506
1507 }
1508
1509 //--------------------------------------------------------------------------
1510
1511 // Implementations of IF clustering maps for massive partons.
1512
1513 bool VinciaCommon::map3to2IFmassive(vector<Vec4>& pClu, vector<Vec4>& pIn,
1514  double saj, double sjk, double sak) {
1515  double sAK = saj + sak - sjk;
1516  Vec4 pA = pIn[0];
1517  pA.rescale4(sAK/(sAK + sjk));
1518  Vec4 pK = pA - pIn[0] + pIn[1] + pIn[2];
1519  pClu.push_back(pA);
1520  pClu.push_back(pK);
1521  return true;
1522
1523 }
1524
1525 //--------------------------------------------------------------------------
1526
1527 // Implementations of II clustering maps for massive partons.
1528
1529 bool VinciaCommon::map3to2IImassive(vector<Vec4>& pClu, vector<Vec4>& pIn,
1530  vector<Vec4>& pRec, double saj, double sjb, double sab, bool doBoost) {
1531
1532  // Scale factors and momenta.
1533  double sAB = sab - saj - sjb;
1534  double rescaleFacA = 1./sqrt(sab/sAB * (sab - saj)/(sab - sjb));
1535  double rescaleFacB = 1./sqrt(sab/sAB * (sab - sjb)/(sab - saj));
1536  Vec4 pA = pIn[0];
1537  pA.rescale4(rescaleFacA);
1538  pClu.push_back(pA);
1539  Vec4 pB = pIn[2];
1540  pB.rescale4(rescaleFacB);
1541  pClu.push_back(pB);
1542  Vec4 pInSum = pIn[0] + pIn[2] - pIn[1];
1543  Vec4 pCluSum = pClu[0] + pClu[1];
1544
1545  // Perform boost - if doBoost, we boost back to the lab frame.
1546  if (doBoost) {
1547  for (int i=0; i<(int)pRec.size(); i++) {
1548  pRec[i].bstback(pInSum);
1549  pRec[i].bst(pCluSum);
1550  }
1551
1552  // Otherwise stay in the current frame. Adjust clustered momenta.
1553  } else {
1554  for (int i=0; i<(int)pClu.size(); i++) {
1555  pClu[i].bstback(pCluSum);
1556  pClu[i].bst(pInSum);
1557  }
1558  }
1559  return true;
1560
1561 }
1562
1563 //--------------------------------------------------------------------------
1564
1565 // 2->3 branching kinematics: output arguments first, then input
1566 // arguments. The output is p[i,j,k] whil the inputs are p[I,K],
1567 // kMapType, inviariants[sIK,sij,sjk], phi, and
1568 // masses[mi,mj,mk]. Note, sab defined as 2*pa.pb.
1569
1570 bool VinciaCommon::map2to3FFmassive(vector<Vec4>& pThree,
1571  const vector<Vec4>& pTwo, int kMapType, const vector<double>& invariants,
1572  double phi, vector<double> masses) {
1573
1574  // Hand off to massless map if all masses << sIK.
1575  if (masses.size() < 3 ||
1576  (masses[0] <= TINY && masses[1] <= TINY && masses[2] <= TINY))
1577  return map2to3FFmassless(pThree,pTwo,kMapType,invariants,phi);
1578
1579  // Antenna invariant mass and sIK = 2*pI.pK.
1580  double m2Ant = m2(pTwo[0],pTwo[1]);
1581  double mAnt = sqrt(m2Ant);
1582  double sAnt = invariants[0];
1583  if (sAnt <= 0.0) return false;
1584
1585  // Masses and invariants
1586  double mass0 = max(0.,masses[0]);
1587  double mass1 = max(0.,masses[1]);
1588  double mass2 = max(0.,masses[2]);
1589  // Check for totally closed phase space. Should normally have
1590  // happened before generation of invariants but put last-resort
1591  // check here since not caught by Gram determinant.
1592  if (mAnt < mass0 + mass1 + mass2) {
1593  cout <<" (VinciaCommon::map2to3FFmassive:) "
1594  <<"ERROR! unphysical phase space\n";
1595  return false;
1596  }
1597  double s01 = max(0.,invariants[1]);
1598  double s12 = max(0.,invariants[2]);
1599  double s02 = m2Ant - s01 - s12 - pow2(mass0) - pow2(mass1) - pow2(mass2);
1600  if (s02 <= 0.) return false;
1601
1602  // Check whether we are inside massive phase space.
1603  double gDet = gramDet(s01, s12, s02, mass0, mass1, mass2);
1604  if (gDet <= 0.) {
1605  if (verbose >= 9)
1606  cout << " map2to3FFmassive: failed massive phase space check" << endl;
1607  return false;
1608  }
1609
1610  // Verbose output.
1611  if (verbose >= 7) {
1612  cout << " (VinciaCommon::map2to3FFmassive:) m = " << num2str(mAnt)
1613  << " sqrtsIK = " << num2str(sqrt(sAnt)) << " sqrts(ij,jk,ik) ="
1614  << num2str(sqrt(s01)) << " " << num2str(sqrt(s12)) << " "
1615  << num2str(sqrt(s02)) << " m(i,j,k) =" << num2str(mass0) << " "
1616  << num2str(mass1) << " " << num2str(mass2) << " D = " << gDet << endl;
1617  RotBstMatrix M;
1618  Vec4 p1cm = pTwo[0];
1619  Vec4 p2cm = pTwo[1];
1620  M.toCMframe(p1cm,p2cm);
1621  p1cm.rotbst(M);
1622  p2cm.rotbst(M);
1623  Vec4 tot = p1cm+p2cm;
1624  cout << " (VinciaCommon::map2to3FFmassive:) starting dipole in CM\n"
1625  << " p1cm = " << p1cm << " p2cm = " << p2cm
1626  << " total = " << tot << endl;
1627  }
1628
1629  // Set up kinematics in rest frame.
1630  double E0 = (pow2(mass0) + s01/2 + s02/2)/mAnt;
1631  double E1 = (pow2(mass1) + s12/2 + s01/2)/mAnt;
1632  double E2 = (pow2(mass2) + s02/2 + s12/2)/mAnt;
1633
1634  // Make sure energies > masses (should normally be ensured by
1635  // combination of open phase space and positive Gram determinant).
1636  if (E0 < mass0 || E1 < mass1 || E2 < mass2) {
1637  cout <<" (VinciaCommon::map2to3FFmassive:) ERROR! "
1638  <<"Unphysical energy value(s)\n";
1639  return false;
1640  }
1641  double ap0 = sqrt( pow2(E0) - pow2(mass0) );
1642  double ap1 = sqrt( pow2(E1) - pow2(mass1) );
1643  double ap2 = sqrt( pow2(E2) - pow2(mass2) );
1644  double cos01 = (E0*E1 - s01/2)/(ap0*ap1);
1645  double cos02 = (E0*E2 - s02/2)/(ap0*ap2);
1646
1647  // Protection: num. precision loss for small (ultracollinear) invariants.
1648  if ( 1-abs(cos01) < 1e-15 ) cos01 = cos01 > 0 ? 1. : -1.;
1649  if ( 1-abs(cos02) < 1e-15 ) cos02 = cos02 > 0 ? 1. : -1.;
1650
1651  // Use positive square root for sine.
1652  double sin01 = (abs(cos01) < 1) ? sqrt(abs(1.0 - pow2(cos01))) : 0.0;
1653  double sin02 = (abs(cos02) < 1) ? sqrt(abs(1.0 - pow2(cos02))) : 0.0;
1654
1655  // Set momenta in CMz frame (frame with 1 oriented along positive z
1656  // axis and event in (x,z) plane).
1657  Vec4 p1(0.0,0.0,ap0,E0);
1658  Vec4 p2(-ap1*sin01,0.0,ap1*cos01,E1);
1659  Vec4 p3(ap2*sin02,0.0,ap2*cos02,E2);
1660
1661  // Verbose output.
1662  if (verbose >= superdebug) {
1663  Vec4 tot = p1+p2+p3;
1664  cout << " (map2to3FFmassive:) configuration in CM* (def: 1 along z)\n";
1665  cout << " k1* = " << p1 << " k2* = " << p2 << " k3* = " << p3
1666  << " total = " << tot << endl;
1667  }
1668
1669  // Choose global rotation around axis perpendicular to event plane.
1670  double psi;
1671
1672  // kMapType = -2(-1): force A(B) to be purely longitudinal recoiler.
1673  if (kMapType == -2) psi = 0.0;
1674  else if (kMapType == -1) psi = M_PI - acos(cos02);
1675  // Else general antenna-like recoils.
1676  else {
1677  double m2I = max(0.0,m2(pTwo[0]));
1678  double m2K = max(0.0,m2(pTwo[1]));
1679  double sig2 = m2Ant + m2I - m2K;
1680  double sAntMin = 2*sqrt(m2I*m2K);
1681  double s01min = max(0.0,2*mass0*mass1);
1682  double s12min = max(0.0,2*mass1*mass2);
1683  double s02min = max(0.0,2*mass0*mass2);
1684
1685  // The r and R parameters in arXiv:1108.6172.
1686  double rAntMap = ( sig2 + sqrt( pow2(sAnt) - pow2(sAntMin) )
1687  * ( s12-s12min - (s01-s01min) )
1688  / ( s01-s01min + s12-s12min ) ) / (2*m2Ant);
1689  double bigRantMap2 = 16*gDet * ( m2Ant*rAntMap * (1.-rAntMap)
1690  - (1.-rAntMap)*m2I - rAntMap*m2K )
1691  + ( pow2(s02) - pow2(s02min) )
1692  * ( pow2(sAnt) - pow2(sAntMin) );
1693  if(bigRantMap2 < 0.){
1694  stringstream ss;
1695  ss<<"On line "<<__LINE__;
1696  infoPtr->errorMsg("Warning in "+__METHOD_NAME__
1697  +": kinematics map is broken.",ss.str());
1698  return false;
1699  }
1700  double bigRantMap = sqrt( bigRantMap2 );
1701  double p1dotpI = (sig2*(pow2(s02) - pow2(s02min))*
1702  (m2Ant + pow2(mass0) - pow2(mass1) - pow2(mass2) - s12)
1703  +8*rAntMap*(m2Ant + pow2(mass0) - pow2(mass1) - pow2(mass2) - s12)*gDet
1704  -bigRantMap*(pow2(s02) - pow2(s02min) + s01*s02-2*s12*pow2(mass0)))
1705  /(4*(4*gDet + m2Ant*(pow2(s02) - pow2(s02min))));
1706
1707  // Norm of the three-momentum and the energy of the first parent
1708  // particle (could also be obtained by explicitly boosting
1709  // incoming particles to CM).
1710  double apInum2 = pow2(m2Ant) + pow2(m2I) + pow2(m2K) - 2*m2Ant*m2I
1711  - 2*m2Ant*m2K - 2*m2I*m2K;
1712  if (apInum2 < 0.) {
1713  stringstream ss;
1714  ss<<"On line "<<__LINE__;
1715  infoPtr->errorMsg("Warning in "+__METHOD_NAME__
1716  +": kinematics map is broken.",ss.str());
1717  return false;
1718  }
1719  double apI = sqrt(apInum2)/(2*mAnt);
1720  double EI = sqrt( pow2(apI) + m2I );
1721
1722  // Protect against rounding errors before taking acos.
1723  double cospsi = ((E0*EI) - p1dotpI)/(ap0*apI);
1724  if (cospsi >= 1.0) {
1725  psi = 0.;
1726  } else if (cospsi <= -1.0) {
1727  psi = M_PI;
1728  } else if(isnan(cospsi)){
1729  psi= 0.;
1730  stringstream ss;
1731  ss << "ap0 = " << ap0;
1732  ss << " apI = " << apI;
1733  ss << " E0 = " << E0;
1734  ss << " mass0 = " << mass0;
1735  ss << " mAnt = " << mAnt;
1736  ss << " sum1 = " << pow2(sAnt) + pow2(m2I) + pow2(m2K);
1737  ss << " sum2 = " << 2*sAnt*m2I + 2*sAnt*m2K + 2*m2I*m2K;
1738  infoPtr->errorMsg("Error in "+__METHOD_NAME__
1739  +": cos(psi) = nan.",ss.str());
1740  return false;
1741  }
1742  else{
1743  psi = acos( cospsi );
1744  }
1745  }
1746
1747  // Perform global rotations.
1748  p1.rot(psi,phi);
1749  p2.rot(psi,phi);
1750  p3.rot(psi,phi);
1751
1752  // Verbose output.
1753  if (verbose >= 7) {
1754  Vec4 tot = p1+p2+p3;
1755  printOut("VinciaCommon::map2to3FFmassive:", "phi = " + num2str(phi,6)
1756  + " cospsi = " + num2str(cos(psi),6) +" psi = " + num2str(psi,6));
1757  printOut("VinciaCommon::map2to3FFmassive:",
1758  "mapType = " + num2str(kMapType));
1759  printOut("VinciaCommon::map2to3FFmassive:","final momenta in CM");
1760  cout << " k1cm = " << p1 << " k2cm = " << p2 << " k3cm = " << p3
1761  << " total = " << tot << endl;
1762  }
1763
1764  // Rotate and boost to lab frame.
1765  RotBstMatrix M;
1766  M.fromCMframe(pTwo[0],pTwo[1]);
1767  if (verbose >= superdebug) {
1768  Vec4 tot = pTwo[0]+pTwo[1];
1769  cout << " (VinciaCommon::map2to3FFmassive) boosting to LAB frame "
1770  << "defined by\n";
1771  cout << " p1 = " << pTwo[0] << " p2 = " << pTwo[1]
1772  << " total = " << tot << endl;
1773  }
1774  p1.rotbst(M);
1775  p2.rotbst(M);
1776  p3.rotbst(M);
1777  if (verbose >= 7) {
1778  Vec4 tot = p1+p2+p3;
1779  cout << " (VinciaCommon::map2to3FFmassive:) final momenta in LAB\n";
1780  cout << " k1 = " << p1 << " k2 = " << p2 << " k3 = " << p3
1781  << " total = " << tot << endl;
1782  }
1783
1784  // Save momenta.
1785  pThree.resize(0);
1786  pThree.push_back(p1);
1787  pThree.push_back(p2);
1788  pThree.push_back(p3);
1789
1790  Vec4 total = pTwo[0] + pTwo[1];
1791  total -= (p1+p2+p3);
1792  if (abs(total.e()) > SMALL || abs(total.px()) > SMALL
1793  || abs(total.py()) > SMALL || abs(total.pz()) > SMALL ){
1794  infoPtr->errorMsg("Error in "+__METHOD_NAME__+
1795  ": Failed momentum conservation test. Aborting.");
1796  infoPtr->setAbortPartonLevel(true);
1797  return false;
1798  }
1799  if (isnan(total)) {
1800  infoPtr->errorMsg("Error in "+__METHOD_NAME__+": (E,p) = nan.");
1801  return false;
1802  }
1803
1804  // Return.
1805  return true;
1806
1807 }
1808
1809 //--------------------------------------------------------------------------
1810
1811 // 2->3 branching kinematics, massless with output arguments first,
1812 // then input arguments. The output is p3, while the inputs are
1813 // kMapType, invariants(sIK, s01, s12), and phi.
1814
1815 bool VinciaCommon::map2to3FFmassless(vector<Vec4>& pThree,
1816  const vector<Vec4>& pTwo, int kMapType, const vector<double>& invariants,
1817  double phi) {
1818
1819  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
1820
1821  // Antenna invariant mass.
1822  double m2Ant = m2(pTwo[0],pTwo[1]);
1823  double mAnt = sqrt(m2Ant);
1824
1825  // Compute invariants (massless relation).
1826  double s01 = invariants[1];
1827  double s12 = invariants[2];
1828  double s02 = m2Ant - s01 - s12;
1829
1830  // Can check alternative hadronization vetos here.
1831  if (verbose >= 7) {
1832  cout << " (VinciaCommon::map2to3FFmassless:) m = " << num2str(mAnt)
1833  << " m12 =" << num2str(sqrt(s01))
1834  << " m23 =" << num2str(sqrt(s12))
1835  << " m13 =" << num2str(sqrt(s02)) << endl;
1836  RotBstMatrix M;
1837  Vec4 p1cm = pTwo[0];
1838  Vec4 p2cm = pTwo[1];
1839  M.toCMframe(p1cm,p2cm);
1840  p1cm.rotbst(M);
1841  p2cm.rotbst(M);
1842  Vec4 tot = p1cm+p2cm;
1843  cout << " (VinciaCommon::map2to3FFmassless:) starting dipole in CM\n"
1844  << " p1cm = " << p1cm << " p2cm = " << p2cm
1845  << " total = " << tot<<endl;
1846  }
1847
1848  // Set up kinematics in rest frame.
1849  double E0 = 1/mAnt*(s01/2 + s02/2);
1850  double E1 = 1/mAnt*(s01/2 + s12/2);
1851  double E2 = 1/mAnt*(s02/2 + s12/2);
1852  double ap0 = E0;
1853  double ap1 = E1;
1854  double ap2 = E2;
1855  double cos01 = (E0*E1 - s01/2)/(ap0*ap1);
1856  double cos02 = (E0*E2 - s02/2)/(ap0*ap2);
1857
1858  // Protection: num. precision loss for small (ultracollinear) invariants.
1859  if ( 1-abs(cos01) < 1e-15 ) cos01 = cos01 > 0 ? 1. : -1.;
1860  if ( 1-abs(cos02) < 1e-15 ) cos02 = cos02 > 0 ? 1. : -1.;
1861  double sin01 = (abs(cos01) < 1) ? sqrt(abs(1.0 - pow2(cos01))) : 0.0;
1862  double sin02 = (abs(cos02) < 1) ? sqrt(abs(1.0 - pow2(cos02))) : 0.0;
1863
1864  // Set momenta in CMz frame (with 1 oriented along positive z axis
1865  // and event in (x,z) plane).
1866  Vec4 p1(0.0,0.0,ap0,E0);
1867  Vec4 p2(-ap1*sin01,0.0,ap1*cos01,E1);
1868  Vec4 p3(ap2*sin02,0.0,ap2*cos02,E2);
1869
1870  // Verbose output.
1871  if (verbose >= superdebug) {
1872  Vec4 tot = p1+p2+p3;
1873  cout << " (map2to3FFmassless:) configuration in CM* (def: 1 along z)\n"
1874  << " k1* = " << p1 << " k2* = " << p2 << " k3* = " << p3
1875  << " total = " << tot << endl;
1876  }
1877
1878  // Choose global rotation around axis perpendicular to event plane.
1879  double psi;
1880
1881  // Force A to be longitudinal recoiler.
1882  if (kMapType == -2) {
1883  psi = 0.0;
1884
1885  // Force B to be longitudinal recoiler.
1886  } else if (kMapType == -1) {
1887  psi = M_PI - acos(max(-1.,min(1.,cos02)));
1888
1890  } else if (kMapType == 1) {
1891  psi = E2*E2/(E0*E0+E2*E2)*(M_PI-acos(cos02));
1892
1893  // psi PYTHIA-like. "Recoiler" remains along z-axis.
1894  } else if (kMapType == 2) {
1895  psi = 0.;
1896  if (s01 < s12 || (s01 == s12 && rndmPtr->flat() > 0.5) )
1897  psi = M_PI-acos(cos02);
1898
1899  // Kosower's map. Similar to ARIADNE.
1900  } else {
1901  double rMap(1);
1902  if (kMapType == 3) rMap = s12/(s01+s12);
1903  double rho=sqrt(1.0+4.0*rMap*(1.0-rMap)*s01*s12/s02/m2Ant);
1904  double s00=-( (1.0-rho)*m2Ant*s02 + 2.0*rMap*s01*s12 ) / 2.0 /
1905  (m2Ant - s01);
1906  psi=acos(1.0+2.0*s00/(m2Ant-s12));
1907  }
1908
1909  // Perform global rotations.
1910  p1.rot(psi,phi);
1911  p2.rot(psi,phi);
1912  p3.rot(psi,phi);
1913
1914  // Verbose output.
1915  if (verbose >= 7) {
1916  Vec4 tot = p1+p2+p3;
1917  printOut(__METHOD_NAME__, "phi = " + num2str(phi,6) + "psi = " +
1918  num2str(psi,6));
1919  printOut(__METHOD_NAME__, "final momenta in CM");
1920  cout << " k1cm = " << p1 << " k2cm = " << p2 << " k3cm = " << p3
1921  << " total = " << tot << endl;
1922  }
1923
1924  // Rotate and boost to lab frame.
1925  RotBstMatrix M;
1926  M.fromCMframe(pTwo[0],pTwo[1]);
1927  Vec4 total = pTwo[0] + pTwo[1];
1928  if (verbose >= superdebug) {
1929  cout << " (VinciaCommon::map2to3FFmassless:) boosting to LAB frame "
1930  << "defined by\n" << " p1 = " << pTwo[0] << " p2 = " << pTwo[1]
1931  << " total = " << total << endl;
1932  }
1933  p1.rotbst(M);
1934  p2.rotbst(M);
1935  p3.rotbst(M);
1936  if (verbose >= 7) {
1937  Vec4 tot = p1 + p2 + p3 ;
1938  printOut(__METHOD_NAME__,"final momenta in LAB");
1939  cout <<" k1 = "<<p1<<" k2 = "<<p2<<" k3 = "<<p3
1940  <<" total = "<<tot<<endl;
1941  }
1942
1943  // Save momenta.
1944  pThree.resize(0);
1945  pThree.push_back(p1);
1946  pThree.push_back(p2);
1947  pThree.push_back(p3);
1948
1949  // Check momentum conservation.
1950  Vec4 diff = total - (p1+p2+p3);
1951  if(abs(diff.e()) / abs(total.e()) > SMALL ||
1952  abs(diff.px()) / abs(total.e()) > SMALL ||
1953  abs(diff.py()) / abs(total.e()) > SMALL ||
1954  abs(diff.pz()) / abs(total.e()) > SMALL) {
1955  infoPtr->errorMsg("Error in "+__METHOD_NAME__
1956  +": (E,p) not conserved.","Aborting.");
1957  cout << setprecision(10) << " difference = " << total.px() << " "
1958  << total.py() << " " << total.pz() << " " << total.e() << endl;
1959  infoPtr->setAbortPartonLevel(true);
1960  return false;
1961  }
1962
1963  // Return.
1964  return true;
1965
1966 }
1967
1968 //--------------------------------------------------------------------------
1969
1970 // Implementations of RF clustering maps for massive partons.
1971
1972 bool VinciaCommon::map2to3RFmassive(vector<Vec4>& pThree, vector<Vec4> pTwo,
1973  vector<double> invariants,double phi,
1974  vector<double> masses) {
1975
1976  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
1977
1978  // Get momenta and boost to lab frame.
1979  if(pTwo.size() != 2){
1980  infoPtr->errorMsg("Error in "+__METHOD_NAME__
1981  +": Wrong number of momenta provided.");
1982  return false;
1983  }
1984
1985  // Momentum of recoiler(s), final state parton, and (modified) resonance.
1986  Vec4 pAKBefore = pTwo.at(0);
1987  Vec4 pKBefore = pTwo.at(1);
1988  Vec4 pABefore = pKBefore + pAKBefore;
1989  Vec4 pACoM = pABefore;
1990
1991  // Boost to CoM frame of (modified) resonance.
1992  pKBefore.bstback(pABefore);
1993  pAKBefore.bstback(pABefore);
1994  pACoM.bstback(pABefore);
1995
1996  // Get the polar and phi angle in CoM frame of K.
1997  double thetaK = pKBefore.theta();
1998  double phiK = pKBefore.phi();
1999
2000
2001  // Construct recoiled momenta in (modified) resonance CoM
2002  // frame. Recover masses and unscaled invariants.
2003  double saj = invariants[1];
2004  double sjk = invariants[2];
2005  double sak = invariants[3];
2006  double mA = masses[0];
2007  double mj = masses[1];
2008  double mk = masses[2];
2009  double mAK = masses[3];
2010  double invDiff = mA*mA + mj*mj + mk*mk -saj-sak+sjk;
2011  invDiff -= mAK*mAK;
2012
2013  // Get energies.
2014  double EjAfter = saj/(2.0*mA);
2015  double EkAfter = sak/(2.0*mA);
2016  if (EkAfter < mk) return false;
2017  else if (EjAfter < mj) return false;
2018  else if (invDiff > SMALL) return false;
2019
2020  // Get cosTheta.
2021  double cosTheta = getCosTheta(EjAfter,EkAfter, mj,mk, sjk);
2022  if (abs(cosTheta) > 1.0) return false;
2023  double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
2024  double pk = sqrt(EkAfter*EkAfter-mk*mk);
2025  double pj = sqrt(EjAfter*EjAfter-mj*mj);
2026
2027  // Construct three momenta, assuming k was along z.
2028  Vec4 pkAfter(0.,0.,pk, EkAfter);
2029  Vec4 pjAfter(pj*sinTheta,0.,pj*cosTheta, EjAfter);
2030  Vec4 pajkAfter = pACoM - pkAfter - pjAfter;
2031
2032  // Give some transverse recoil to k.
2033  double thetaEff = -(M_PI-pajkAfter.theta());
2034  pkAfter.rot(thetaEff,0.);
2035  pjAfter.rot(thetaEff,0.);
2036  pajkAfter.rot(thetaEff,0.);
2037
2038  // Rotate by arbitrary phi.
2039  pkAfter.rot(0.,phi);
2040  pjAfter.rot(0.,phi);
2041  pajkAfter.rot(0.,phi);
2042
2043  // Rotate to recover original orientation of frame.
2044  pkAfter.rot(thetaK,phiK);
2045  pjAfter.rot(thetaK,phiK);
2046  pajkAfter.rot(thetaK,phiK);
2047
2048  // Boost to lab frame.
2049  pkAfter.bst(pABefore);
2050  pjAfter.bst(pABefore);
2051  pajkAfter.bst(pABefore);
2052  pThree.clear();
2053  pThree.push_back(pajkAfter);
2054  pThree.push_back(pjAfter);
2055  pThree.push_back(pkAfter);
2056
2057  // Return.
2058  return true;
2059
2060 }
2061
2062 //--------------------------------------------------------------------------
2063
2064 // Implementations of resonance kineatic maps for massive partons. Inputs
2065 // are as follows:
2066 // pBefore = momenta of resonance and all downstream recoilers
2067 // before emission.
2068 // posF = position in pBefore of the momentum of the F end of the
2069 // antenna.
2070 // invariants = yaj and yjk scaled invariants.
2071 // phi = azimuthal angle of gluon emission.
2072 // mui = masses of a, j, k.
2073 // The output is as follows:
2074 // pAfter = momenta of resonance, emission and all downstream recoilers
2075 // after emission.
2076 // [0] = pa - will be unchanged
2077 // [1] = pj
2078 // [2] = pk
2079 // [i>3] = recoilers
2080
2081 bool VinciaCommon::map2toNRFmassive(vector<Vec4>& pAfter, vector<Vec4> pBefore,
2082  unsigned int posR, unsigned int posF, vector<double> invariants,double phi,
2083  vector<double> masses) {
2084
2085  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
2086  pAfter.clear();
2087
2088  // Momentum of "R", "F" end of antenna, and sum of downstream recoilers.
2089  Vec4 pR = pBefore.at(posR);
2090  Vec4 pF = pBefore.at(posF);
2091  Vec4 pSum(0.,0.,0.,0.);
2092  vector<Vec4> pRec;
2093  for(unsigned int imom = 0; imom < pBefore.size(); imom++){
2094  if (imom==posF || imom==posR) {
2095  continue;
2096  } else {
2097  pSum+= pBefore.at(imom);
2098  pRec.push_back(pBefore.at(imom));
2099  }
2100  }
2101  vector<Vec4> pTwo;
2102  vector<Vec4> pThree;
2103  pTwo.push_back(pSum);
2104  pTwo.push_back(pF);
2105
2106  // Recoil AK system.
2107  if (!map2to3RFmassive(pThree,pTwo,invariants,phi,masses)) {
2108  return false;
2109  } else if (pThree.size() != 3) {
2110  return false;
2111  }
2112
2113  // Add pa, pj, and k. Check mass.
2114  pAfter.push_back(pR);
2115  pAfter.push_back(pThree.at(1));
2116  pAfter.push_back(pThree.at(2));
2117  Vec4 pSumAfter = pThree.at(0);
2118  if (abs(pSumAfter.mCalc() - pSum.mCalc()) > SMALL) {
2119  infoPtr->errorMsg("Error in "+__METHOD_NAME__
2120  +": Failed to conserve mass of system.");
2121  return false;
2122  }
2123
2124  // If only a single recoiler, it just takes the remaining momentum.
2125  if (pRec.size() == 1) {
2126  pAfter.push_back(pSumAfter);
2127
2128  // Boost the downstream recoilers appropriately
2129  } else {
2130  for(unsigned int imom = 0; imom < pRec.size(); imom++) {
2131  double mRecBef = pRec.at(imom).mCalc();
2132  pRec.at(imom).bstback(pSum,pSum.mCalc());
2133  pRec.at(imom).bst(pSumAfter,pSum.mCalc());
2134  double mRecAfter = pRec.at(imom).mCalc();
2135
2136  // Check mass.
2137  if (abs(mRecAfter- mRecBef) > SMALL) {
2138  infoPtr->errorMsg("Error in "+__METHOD_NAME__
2139  +": Failed to conserve mass of recoilers.");
2140  return false;
2141  }
2142  pAfter.push_back(pRec.at(imom));
2143  }
2144  }
2145
2146  // Done.
2147  return true;
2148 }
2149
2150 //--------------------------------------------------------------------------
2151
2152 // 2 -> 3 kinematics map for initial-initial antennae, for general mj.
2153
2154 bool VinciaCommon::map2to3II(vector<Vec4>& pNew, vector<Vec4>& pRec,
2155  vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
2156  double phi, double mj2) {
2157
2158  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
2159
2160  // Hand off to massless map if mj2 = 0.
2161  if (mj2 == 0.0)
2162  return map2to3IImassless(pNew, pRec, pOld, sAB, saj, sjb, sab, phi);
2163
2164  // Do massive mapping.
2165  pNew.clear();
2166  pNew.resize(3);
2167
2168  // Force incoming momenta on shell (massless) with mass squared = sAB.
2169  pOld[0].py(0.);
2170  pOld[0].px(0.);
2171  pOld[1].py(0.);
2172  pOld[1].px(0.);
2173  double sCM = m2(pOld[0] + pOld[1]);
2174  double fac = sqrt(sAB/sCM);
2175  double e0 = pOld[0].e();
2176  double e1 = pOld[1].e();
2177  if (abs(1. - fac) > TINY) {
2178  if (verbose >= 3 && abs(1.-fac) > 1.01)
2179  printOut("VinClu::map2to3II", "Warning: scaling AB so m2(AB) = sAB");
2180  e0 *= fac;
2181  e1 *= fac;
2182  }
2183  int sign = pOld[0].pz() > 0 ? 1 : -1;
2184  pOld[0].e(e0);
2185  pOld[0].pz(sign * e0);
2186  pOld[1].e(e1);
2187  pOld[1].pz(-sign * e1);
2188
2189  // Initialise new momenta.
2190  pNew[0] = pOld[0];
2191  pNew[2] = pOld[1];
2192
2193  // Check if we're inside massive phase space.
2194  double G = saj*sjb*sab - mj2*sab*sab;
2195  if (G < 0. || sab < 0.) return false;
2196  if ((sab <= sjb) || (sab <= saj)) {
2197  infoPtr->errorMsg("Error in "+__METHOD_NAME__
2198  +": Incompatible invariants.");
2199  return false;
2200  }
2201
2202  // Incoming momenta.
2203  double rescaleFacA = sqrt(sab/sAB * (sab - saj)/(sab - sjb));
2204  double rescaleFacB = sqrt(sab/sAB * (sab - sjb)/(sab - saj));
2205  pNew[0].rescale4(rescaleFacA);
2206  pNew[2].rescale4(rescaleFacB);
2207
2208  // Emission.
2209  double preFacA = sjb*sqrt((sab - saj)/(sab - sjb)/sab/sAB);
2210  double preFacB = saj*sqrt((sab - sjb)/(sab - saj)/sab/sAB);
2211  double preFacT = sqrt(saj*sjb/sab - mj2);
2212  Vec4 pTrans(cos(phi), sin(phi), 0.0, 0.0);
2213  pNew[1] = preFacA*pOld[0] + preFacB*pOld[1] + preFacT*pTrans;
2214
2215  if (verbose >= superdebug) {
2216  printOut("VinClu::map2to3II","Invariants are");
2217  cout << scientific << " sAB = " << sAB << " saj = " << saj
2218  << " sjb = " << sjb << " sab = " << sab << endl
2219  << " Given momenta are" << endl;
2220  for (int i = 0; i < 2; i++) cout << " " << pOld[i];
2221  cout << " New momenta are" << endl;
2222  for (int i = 0; i < 3; i++) cout << " " << pNew[i];
2223  }
2224
2225  // Check the invariants, allow 0.0001% difference.
2226  double check = 1e-6;
2227  double sajNew = 2*pNew[0]*pNew[1];
2228  double sjbNew = 2*pNew[1]*pNew[2];
2229  double sabNew = 2*pNew[0]*pNew[2];
2230  if (abs(sabNew - sab)/sab > check) {
2231  if (verbose >= 5) {
2232  printOut("VinClu::map2to3II","ERROR! Invariants differ!");
2233  cout << scientific << " sab (" << sab << ") fracdiff = ydiff = "
2234  << abs(sabNew-sab)/sab << endl << " Old momenta are" << endl;
2235  for (int i=0; i<2; i++) cout << " " << pOld[i];
2236  cout << " New momenta are" << endl;
2237  for (int i=0; i<3; i++) cout << " " << pNew[i];
2238  }
2239  return false;
2240  } else if (abs(sajNew - saj)/sab > check) {
2241  if (verbose >= 5) {
2242  printOut("VinClu::map2to3II","ERROR! Invariants differ!");
2243  cout << scientific << " saj (" << saj << ") fracdiff = "
2244  << abs(sajNew-saj)/saj << " ydiff = "
2245  << abs(sajNew-saj)/sab << endl << " Old momenta are" << endl;
2246  for (int i=0; i<2; i++) cout << " " << pOld[i];
2247  cout << " New momenta are" << endl;
2248  for (int i=0; i<3; i++) cout << " " << pNew[i];
2249  }
2250  return false;
2251  } else if (abs(sjbNew - sjb)/sab > check) {
2252  if (verbose >= 5) {
2253  printOut("VinClu::map2to3II","ERROR! Invariants differ!");
2254  cout << scientific << " sjb (" << sjb << ") fracdiff = "
2255  << abs(sjbNew-sjb)/sjb << " ydiff = " << abs(sjbNew-sjb)/sab
2256  << endl << " Old momenta are" << endl;
2257  for (int i=0; i<2; i++) cout << " " << pOld[i];
2258  cout << " New momenta are" << endl;
2259  for (int i=0; i<3; i++) cout << " " << pNew[i];
2260  }
2261  return false;
2262  }
2263
2264  // Change the final state recoiler. The recoiler is currently sum
2265  // of initial guys => E,0,0,pz. Boost in center-of-mass frame AB
2266  // E,0,0,0.
2267  Vec4 pSum = pOld[0] + pOld[1];
2268  Vec4 pRecSumBefore(0.,0.,0.,0.);
2269  Vec4 pRecSumAfter(0.,0.,0.,0.);
2270  for (int i=0; i<(int)pRec.size(); i++){
2271  pRecSumBefore+=pRec[i];
2272  pRec[i].bstback(pSum);
2273  }
2274
2275  // Now boost from E,0,0,0 to E',px',py',pz'.
2276  Vec4 pPrime = pNew[0] + pNew[2] - pNew[1];
2277  for (int i=0; i<(int)pRec.size(); i++) {
2278  pRec[i].bst(pPrime, pPrime.mCalc());
2279  pRecSumAfter+=pRec[i];
2280  }
2281  if(verbose >= superdebug) {
2282  Vec4 total= pOld[0]+pOld[1];
2283  cout << " Total In before" << total << endl
2284  << " Total Out before" << pRecSumBefore << endl;
2285  total = pNew[0] + pNew[2] - pNew[1];
2286  cout << " Total In After" << total << endl
2287  << " Total Out After" << pRecSumAfter << endl
2288  << " Total diff After" << total-pRecSumAfter << endl;
2289  }
2290  return true;
2291
2292 }
2293
2294 //--------------------------------------------------------------------------
2295
2296 // 2 -> 3 kinematics map for initial-initial antennae, for mj= 0.
2297
2298 bool VinciaCommon::map2to3IImassless(vector<Vec4>& pNew, vector<Vec4>& pRec,
2299  vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
2300  double phi) {
2301
2302  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
2303  pNew.clear();
2304  pNew.resize(3);
2305
2306  // Force incoming momenta on shell (massless) with mass squared = sAB.
2307  pOld[0].py(0.);
2308  pOld[0].px(0.);
2309  pOld[1].py(0.);
2310  pOld[1].px(0.);
2311  double sCM = m2(pOld[0] + pOld[1]);
2312  double fac = sqrt(sAB/sCM);
2313  double e0 = pOld[0].e();
2314  double e1 = pOld[1].e();
2315  if (abs(1. - fac) > TINY) {
2316  if (verbose >= 3 && abs(1. - fac) > 1.01)
2317  printOut("VinciaCommon::map2to3IImassless",
2318  "Warning: scaling AB so m2(AB) = sAB");
2319  e0 *= fac;
2320  e1 *= fac;
2321  }
2322  int sign = pOld[0].pz() > 0 ? 1 : -1;
2323  pOld[0].e(e0);
2324  pOld[0].pz(sign * e0);
2325  pOld[1].e(e1);
2326  pOld[1].pz(-sign * e1);
2327
2328  // Initialise new momenta.
2329  pNew[0] = pOld[0];
2330  pNew[2] = pOld[1];
2331
2332  // Incoming momenta.
2333  double rescaleFacA = sqrt(sab/(sAB+saj) * (1. + sjb/sAB));
2334  double rescaleFacB = sqrt(sab/(sAB+sjb) * (1. + saj/sAB));
2335  pNew[0].rescale4(rescaleFacA);
2336  pNew[2].rescale4(rescaleFacB);
2337
2338  // Emission.
2339  double preFacA = sjb*sqrt((sAB+sjb)/(sAB+saj)/sab/sAB);
2340  double preFacB = saj*sqrt((sAB+saj)/(sAB+sjb)/sab/sAB);
2341  double preFacT = sqrt(saj*sjb/sab);
2342  Vec4 pTrans(cos(phi), sin(phi), 0.0, 0.0);
2343  pNew[1] = preFacA*pOld[0] + preFacB*pOld[1] + preFacT*pTrans;
2344
2345  // Debugging info.
2346  if (verbose >= superdebug) {
2347  stringstream ss;
2348  ss << "Invariants are: "
2349  << scientific << " sAB = " << sAB << " saj = " << saj
2350  << " sjb = " << sjb << " sab = " << sab;
2351  printOut(__METHOD_NAME__, ss.str());
2352  printOut(__METHOD_NAME__, "Given momenta are");
2353  for (int i = 0; i < 2; i++) cout << " " << pOld[i];
2354  printOut(__METHOD_NAME__, "New momenta are");
2355  for (int i = 0; i < 3; i++) cout << " " << pNew[i];
2356  }
2357
2358  // Check the invariants, allow 0.0001% difference.
2359  double check = 1e-6;
2360  double sajNew = 2*pNew[0]*pNew[1];
2361  double sjbNew = 2*pNew[1]*pNew[2];
2362  double sabNew = 2*pNew[0]*pNew[2];
2363  if (abs(sabNew - sab)/sab > check) {
2364  if (verbose >= 5) {
2365  printOut("VinciaCommon::map2to3IImassless","ERROR! Invariants differ!");
2366  cout << scientific << " sab (" << sab << ") fracdiff = ydiff = "
2367  << abs(sabNew - sab)/sab << endl
2368  << " Old momenta are" << endl;
2369  for (int i=0; i<2; i++) cout << " " << pOld[i];
2370  cout << " New momenta are" << endl;
2371  for (int i=0; i<3; i++) cout << " " << pNew[i];
2372  }
2373  return false;
2374  } else if (abs(sajNew - saj)/sab > check) {
2375  if (verbose >= 5) {
2376  printOut("VinciaCommon::map2to3IImassless","ERROR! Invariants differ!");
2377  cout << scientific << " saj (" << saj << ") fracdiff = "
2378  << abs(sajNew-saj)/saj << " ydiff = "<< abs(sajNew - saj)/sab
2379  << endl << " Old momenta are" << endl;
2380  for (int i=0; i<2; i++) cout << " " << pOld[i];
2381  cout << " New momenta are" << endl;
2382  for (int i=0; i<3; i++) cout << " " << pNew[i];
2383  }
2384  return false;
2385  } else if ( abs(sjbNew-sjb)/sab > check ) {
2386  if (verbose >= 5) {
2387  printOut("VinciaCommon::map2to3IImassless","ERROR! Invariants differ!");
2388  cout << scientific << " sjb (" << sjb << ") fracdiff = "
2389  << abs(sjbNew-sjb)/sjb << " ydiff = "<< abs(sjbNew - sjb)/sab
2390  << endl << " Old momenta are" << endl;
2391  for (int i=0; i<2; i++) cout << " " << pOld[i];
2392  cout << " New momenta are" << endl;
2393  for (int i=0; i<3; i++) cout << " " << pNew[i];
2394  }
2395  return false;
2396  }
2397
2398  // Change the final state recoiler. The recoiler is currently sum
2399  // of initial guys => E,0,0,pz. Boost in center-of-mass frame AB
2400  // E,0,0,0.
2401  Vec4 pSum = pOld[0] + pOld[1];
2402  for (int i=0; i<(int)pRec.size(); i++) pRec[i].bstback(pSum);
2403
2404  // Now boost from E,0,0,0 to E',px',py',pz' and return.
2405  Vec4 pPrime = pNew[0] + pNew[2] - pNew[1];
2406  for (int i=0; i<(int)pRec.size(); i++) pRec[i].bst(pPrime);
2407  return true;
2408
2409 }
2410
2411 //--------------------------------------------------------------------
2412
2413 // 2->3 kinematics map for local recoils, for general mj,mk. Assumes
2414 // partons from proton explicitly massless
2415
2416 bool VinciaCommon::map2to3IFlocal(vector<Vec4>& pNew, vector<Vec4>& pOld,
2417  double sAK, double saj, double sjk, double sak, double phi,
2418  double mK2, double mj2, double mk2) {
2419
2420  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
2421  pNew.clear();
2422  pNew.resize(3);
2423  if (verbose >= superdebug) {
2424  printOut("VinClu::map2to3IFlocal","Invariants are");
2425  cout << " sAK = " << sAK << " saj = " << saj
2426  << " sjk = " << sjk << " sak = " << sak << endl
2427  << " mK = " << sqrt(mK2) << " mj = " << sqrt(mj2)
2428  << " mk = " << sqrt(mk2) << endl
2429  << " Given momenta are" << endl;
2430  for (int i=0; i<2; i++) cout << " " << pOld[i];
2431  }
2432
2433  // Check invariants.
2434  double check = 1.e-3;
2435  double inv1Norm = (saj + sak)/(sAK + sjk);
2436  double inv2Norm = 1.0 + (mj2 + mk2 - mK2)/(sAK + sjk);
2437  double diff = abs(inv1Norm-inv2Norm);
2438  if(diff > check) {
2439  if (verbose >= 2) {
2440  infoPtr->errorMsg("Error in "+__METHOD_NAME__
2441  +": Inconsistent invariants.","Aborting.");
2442  cout <<" yaj + yak = " << inv1Norm
2443  << " 1 + muj2 + muk2 - muK2 = "<< inv2Norm
2444  << " Diff = " << diff << endl;
2445  }
2446  infoPtr->setAbortPartonLevel(true);
2447  return false;
2448  }
2449
2450  // Check if we're inside massive phase space.
2451  double G = saj*sjk*sak - mj2*sak*sak - mk2*saj*saj;
2452  if (G < 0. || sak < 0.) return false;
2453
2454  // Set up some variables for boosting pT.
2455  Vec4 pSum = pOld[0] + pOld[1];
2456  Vec4 pOldBst = pOld[0];
2457  pOldBst.bstback(pSum);
2458  double thetaRot = pOldBst.theta();
2459  double phiRot = pOldBst.phi();
2460  Vec4 pTrans(cos(phi), sin(phi), 0.0, 0.0);
2461
2462  // Rotate and boost.
2463  pTrans.rot(thetaRot, phiRot);
2464  pTrans.bst(pSum);
2465
2466  // Check if pT was properly boosted, allow 0.1% difference.
2467  if (pTrans*pOld[0] > pOld[0][0]*1e-3 || pTrans*pOld[1] > pOld[1][0]*1e-3) {
2468  if (verbose >= normal) {
2469  infoPtr->errorMsg("Error in "+__METHOD_NAME__
2470  +": The transverse momentum is not transverse after boosting");
2471  }
2472  return false;
2473  }
2474  double sig = sak + saj;
2475  double cj = (sig*(sak + mj2 - mk2) + mK2*(sak - saj) - sAK*sak)/(sAK*sig);
2476  double ck = (sig*(saj - mj2 + mk2) + mK2*(saj - sak) - sAK*saj)/(sAK*sig);
2477  double dj = saj/sig;
2478  double dk = sak/sig;
2479
2480  // Construct new momenta p_a, p_j, and p_k.
2481  pNew[0] = pOld[0];
2482  pNew[0].rescale4(sig/sAK);
2483  pNew[1] = cj*pOld[0] + dj*pOld[1] + (sqrt(G)/sig)*pTrans;
2484  pNew[2] = ck*pOld[0] + dk*pOld[1] - (sqrt(G)/sig)*pTrans;
2485
2486  // Check the invariants, allow 0.1% difference (due to boost).
2487  double sakNew = pNew[0]*pNew[2]*2;
2488  double sajNew = pNew[0]*pNew[1]*2;
2489  double sjkNew = pNew[1]*pNew[2]*2;
2490  if (abs(sakNew - sak)/sak > check) {
2491  if (verbose >= 2) {
2492  printOut("VinClu::map2to3IFlocal","ERROR! sak is inconsistent!");
2493  cout << scientific << " sak (" << sak << ") diff = "
2494  << abs(sakNew-sak)/sak << endl
2495  << " Old momenta are" << endl;
2496  for (int i=0; i<2; i++) cout << " " << pOld[i];
2497  cout << " New momenta are" << endl;
2498  for (int i=0; i<3; i++) cout << " " << pNew[i];
2499  cout <<"Masses: mK2 = " << mK2 << " mj2 = " << mj2
2500  << " mk2 = " << mk2 << endl;
2501  }
2502  return false;
2503  } else if (abs(sajNew - saj)/saj > check) {
2504  if (verbose >= 2 ) {
2505  printOut("VinClu::map2to3IFlocal","ERROR! saj is inconsistent!");
2506  cout << scientific << " saj (" << saj << ") diff = ";
2507  cout << abs(sajNew-saj)/saj << endl;
2508  cout << " Old momenta are" << endl;
2509  for (int i=0; i<2; i++) cout << " " << pOld[i];
2510  cout << " New momenta are" << endl;
2511  for (int i=0; i<3; i++) cout << " " << pNew[i];
2512  cout <<"Masses: mK2 = " << mK2 << " mj2 = " << mj2
2513  << " mk2 = " << mk2 << endl;
2514  }
2515  return false;
2516  } else if ( abs(sjkNew-sjk)/sjk > check ) {
2517  if (verbose >= 2 ){
2518  printOut("VinClu::map2to3IFlocal","ERROR! sjk is inconsistent!");
2519  cout << scientific << " sjk (" << sjk << ") diff = ";
2520  cout << abs(sjkNew-sjk)/sjk << endl;
2521  cout << " Old momenta are" << endl;
2522  for (int i=0; i<2; i++) cout << " " << pOld[i];
2523  cout << " New momenta are" << endl;
2524  for (int i=0; i<3; i++) cout << " " << pNew[i];
2525  cout <<"Masses: mK2 = " << mK2 << " mj2 = " << mj2
2526  << " mk2 = " << mk2 << endl;
2527  }
2528  infoPtr->setAbortPartonLevel(true);
2529  return false;
2530  }
2531  return true;
2532 }
2533
2534 //--------------------------------------------------------------------------
2535
2536 // 2->3 kinematics map for global recoils, for general mj,mk. Assumes
2537 // partons from proton explicitly massless.
2538
2539 bool VinciaCommon::map2to3IFglobal(vector<Vec4>& pNew,
2540  vector<Vec4>& pRec, vector<Vec4>& pOld, Vec4 &pB,
2541  double sAK, double saj, double sjk, double sak, double phi,
2542  double mK2, double mj2, double mk2) {
2543
2544  if (verbose >= superdebug) printOut(__METHOD_NAME__, "begin --------------");
2545  pNew.clear();
2546  pNew.resize(3);
2547
2548  // Set up some variables for boosting pT.
2549  Vec4 pSum = pOld[0] + pOld[1];
2550  Vec4 pAcm = pOld[0];
2551  pAcm.bstback(pSum);
2552  double thetaRot = pAcm.theta();
2553  double phiRot = pAcm.phi();
2554  Vec4 pTrans(cos(phi), sin(phi), 0.0, 0.0);
2555
2556  // Rotate and boost.
2557  pTrans.rot(thetaRot, phiRot);
2558  pTrans.bst(pSum);
2559
2560  // Check if pT was properly boosted, allow 0.1% difference.
2561  if (pTrans*pOld[0] > pOld[0].e()*1e-3 || pTrans*pOld[1] > pOld[1].e()*1e-3) {
2562  if (verbose >= superdebug) {
2563  printOut("VinciaCommon:map2to3IFglobal",
2564  "The transverse momentum is not transverse after boosting");
2565  }
2566  return false;
2567  }
2568
2569  // Solution vector, start from massless values.
2570  vector<double> v; v.resize(5);
2571  double sig = sAK - saj;
2572  v[0] = sak/sig;
2573  v[1] = (saj*sjk)/(sAK*sig);
2574  v[2] = sjk/sig;
2575  v[3] = (saj*sak)/(sAK*sig);
2576  v[4] = sjk*saj*sak/pow2(sig);
2577
2578  // Root finding with Newton-Raphson in 5D.
2579  int nCount = 0;
2580  while (true) {
2581  nCount ++;
2582
2583  // Construct function.
2584  vector<double> f(5, 0);
2585  f[0] = v[0]*v[1]*sAK + pow2(v[1])*mK2 - v[4];
2586  f[1] = v[2]*v[3]*sAK + pow2(v[3])*mK2 - v[4] - mj2;
2587  f[2] = (v[0] - v[2] - 1)*(v[1] - v[3] + 1)*sAK
2588  + pow2(v[1] - v[3] + 1)*mK2 - mk2;
2589  f[3] = (v[0]*v[3] + v[1]*v[2])*sAK + 2*v[1]*v[3]*mK2 - 2*v[4] - saj;
2590  f[4] = (v[2]*(v[1] - v[3] + 1) + v[3]*(v[0] - v[2] - 1))*sAK
2591  + 2*v[3]*(v[1] - v[3] + 1)*mK2 - sjk;
2592
2593  // Construct Jacobian.
2594  vector<vector<double> > A(5, vector<double>(5, 0));
2595  A[0][0] = v[1]*sAK;
2596  A[0][1] = v[0]*sAK + 2*v[1]*mK2;
2597  A[0][2] = 0;
2598  A[0][3] = 0;
2599  A[0][4] = -1;
2600  A[1][0] = 0;
2601  A[1][1] = 0;
2602  A[1][2] = v[3]*sAK;
2603  A[1][3] = v[2]*sAK + 2*v[3]*mK2;
2604  A[1][4] = -1;
2605  A[2][0] = (v[1] - v[3] + 1)*sAK;
2606  A[2][1] = (v[0] - v[2] - 1)*sAK + 2*(v[1] - v[3] + 1)*mK2;
2607  A[2][2] = -(v[1] - v[3] + 1)*sAK;
2608  A[2][3] = -( (v[0] - v[2] - 1)*sAK + 2*(v[1] - v[3] + 1)*mK2 );
2609  A[2][4] = 0;
2610  A[3][0] = v[3]*sAK;
2611  A[3][1] = v[2]*sAK + 2*v[3]*mK2;
2612  A[3][2] = v[1]*sAK;
2613  A[3][3] = v[0]*sAK + 2*v[1]*mK2;
2614  A[3][4] = -2;
2615  A[4][0] = v[3]*sAK;
2616  A[4][1] = v[2]*sAK + 2*v[3]*mK2;
2617  A[4][2] = (v[1] - 2*v[3] + 1)*sAK;
2618  A[4][3] = (v[0] - 2*v[2] - 1)*sAK + 2*(v[1] - 2*v[3] + 1)*mK2;
2619  A[4][4] = 0;
2620
2621  // Invert Jacobian and append identity.
2622  int n = 5;
2623  for (int i = 0; i < n; i++) {
2624  A[i].resize(2*n);
2625  A[i][n+i] = 1;
2626  }
2627
2628  for (int i = 0; i < n; i++) {
2629  // Find column max.
2630  double eleMax = abs(A[i][i]);
2631  int rowMax = i;
2632  for (int k=i+1; k<n; k++) {
2633  if (abs(A[k][i]) > eleMax) {
2634  eleMax = A[k][i];
2635  rowMax = k;
2636  }
2637  }
2638
2639  // Swap maximum row with current row.
2640  A[rowMax].swap(A[i]);
2641
2642  // Reduce current column.
2643  for (int k = i+1; k < n; k++) {
2644  double c = -A[k][i]/A[i][i];
2645  for (int j = i; j < 2*n; j++) {
2646  if (i == j) {
2647  A[k][j] = 0;
2648  } else {
2649  A[k][j] += c * A[i][j];
2650  }
2651  }
2652  }
2653  }
2654
2655  // Solve equation Ax = b for upper triangular matrix A.
2656  for (int i = n-1; i >= 0; i--) {
2657  for (int k = n; k < 2*n; k++) {
2658  A[i][k] /= A[i][i];
2659  }
2660  for (int j = i-1; j >= 0; j--) {
2661  for (int k = n; k < 2*n; k++) {
2662  A[j][k] -= A[i][k] * A[j][i];
2663  }
2664  }
2665  }
2666
2667  // Remove remaining identity.
2668  for (int i = 0; i < n; i++) {
2669  A[i].erase(A[i].begin(), A[i].begin()+n);
2670  }
2671
2672  // Find next iteration.
2673  vector<double> vNew(5, 0);
2674  for (int i = 0; i < n; i++) {
2675  vNew[i] = v[i];
2676  for (int j=0; j<n; j++) {
2677  vNew[i] -= A[i][j]*f[j];
2678  }
2679  }
2680
2681  // Check if done.
2682  double eps = 0;
2683  for (int i=0; i<n; i++) {
2684  if (abs(vNew[i] - v[i])/abs(v[i]) > eps) {
2685  eps = abs(vNew[i] - v[i])/abs(v[i]);
2686  }
2687  }
2688  v = vNew;
2689
2690  // Break if iterations/precision reached or randomly vary the variables.
2691  if (nCount == 1000) {return false;}
2692  if ((nCount%100) == 0) {
2693  for (int i=0; i<n; i++) {
2694  v[i] *= (2*rndmPtr->flat() - 1);
2695  }
2696  }
2697  if (eps < 1E-6) {break;}
2698  }
2699
2700  // Construct post-branching momenta.
2701  pNew[0] = v[0]*pOld[0] + v[1]*pOld[1] - sqrt(v[4])*pTrans;
2702  pNew[1] = v[2]*pOld[0] + v[3]*pOld[1] - sqrt(v[4])*pTrans;
2703  pNew[2] = (v[0] - v[2] - 1)*pOld[0] + (v[1] - v[3] + 1)*pOld[1];
2704
2705  // Set up the boost.
2706  Vec4 pa = pNew[0];
2707  Vec4 pA = pOld[0];
2708  double qaB = pa*pB;
2709  double qAB = pA*pB;
2710  double qAa = pA*pa;
2711
2712  // Perform boost.
2713  for (int i=0; i<3; i++) {
2714  Vec4 p = pNew[i];
2715  pNew[i] += pB*((pa*p)/qaB) - pa*((pB*p)/qaB) + pA*((pB*p)/qAB)
2716  - pB*((pA*p)/qAB) + pB*(qAa*(pB*p)/(qAB*qaB));
2717
2718  // Force the initial state to be on the beam axis.
2719  if (i==0) {
2720  double ea = pNew[i].e();
2721  double sign = (pNew[0].pz() > 0) ? 1 : -1;
2722  pNew[0] = Vec4(0, 0, sign*ea, ea);
2723  }
2724  }
2725
2726  // Perform boost on the rest of the system and return.
2727  for (int i=0; i<(int)pRec.size(); i++) {
2728  Vec4 p = pRec[i];
2729  pRec[i] += pB*((pa*p)/qaB) - pa*((pB*p)/qaB) + pA*((pB*p)/qAB)
2730  - pB*((pA*p)/qAB) + pB*(qAa*(pB*p)/(qAB*qaB));
2731  }
2732  double sajNew = 2*pNew[0]*pNew[1];
2733  double sjkNew = 2*pNew[1]*pNew[2];
2734  double sakNew = 2*pNew[0]*pNew[2];
2735  if (verbose >= 5) {
2736  if (abs(sajNew - saj)/saj > 1E-3)
2737  printOut("VinciaCommon:map2to3IFglobal", "saj not quite correct");
2738  if (abs(sjkNew - sjk)/sjk > 1E-3)
2739  printOut("VinciaCommon:map2to3IFglobal", "sjk not quite correct");
2740  if (abs(sakNew - sak)/sak > 1E-3)
2741  printOut("VinciaCommon:map2to3IFglobal", "sak not quite correct");
2742  }
2743  return true;
2744 }
2745
2746 //--------------------------------------------------------------------------
2747
2748 // pA + pK -> pa + pj + pk
2749
2750 bool VinciaCommon::map2to3RFmassive(vector<Vec4>& pNew, vector<Vec4>& pRec,
2751  vector<Vec4> pOld, double saj, double sjk, double phi,
2752  double mA2 = 0.0, double mj2 = 0.0, double mK2 = 0.0) {
2753
2754  // Get momenta and boost to lab frame.
2755  if (pOld.size() != 2) {return false;}
2756  Vec4 pABefore = pOld[0];
2757  Vec4 pKBefore = pOld[1];
2758  Vec4 pAKBefore = pABefore - pKBefore;
2759  Vec4 pACoM = pABefore;
2760  double sAK = 2.0*pABefore*pKBefore;
2761  double sak = sAK - saj + sjk;
2762
2763  // Check if inside phase space boundaries.
2764  if (sak < 0) {return false;}
2765  if (saj*sjk*sak - mA2*sjk*sjk - mj2*sak*sak - mK2*saj*saj < 0) return false;
2766
2767  // Boost to CoM frame.
2768  pKBefore.bstback(pABefore);
2769  pAKBefore.bstback(pABefore);
2770  pACoM.bstback(pABefore);
2771
2772  // Get energies, cos(theta), and sin(theta).
2773  double EjAfter = saj/(2.0*sqrt(mA2));
2774  double pVecjAfter = sqrt(pow2(EjAfter) - mj2);
2775  double EkAfter = sak/(2.0*sqrt(mA2));
2776  double pVeckAfter = sqrt(pow2(EkAfter) - mK2);
2777  double cosTheta = (2.*EjAfter*EkAfter - sjk)/2./pVecjAfter/pVeckAfter;
2778  if (abs(cosTheta) > 1) {return false;}
2779  double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
2780
2781  // Construct three momenta.
2782  Vec4 pkAfter(0., 0., pVeckAfter, EkAfter);
2783  Vec4 pjAfter(pVecjAfter*sinTheta*cos(phi), pVecjAfter*sinTheta*sin(phi),
2784  pVecjAfter*cosTheta, EjAfter);
2785  Vec4 pajkAfter = pACoM - pkAfter - pjAfter;
2786
2787  // Boost to lab frame.
2788  pkAfter.bst(pABefore);
2789  pjAfter.bst(pABefore);
2790  pajkAfter.bst(pABefore, sqrt(mA2));
2791  pNew.clear();
2792  pNew.push_back(pABefore);
2793  pNew.push_back(pjAfter);
2794  pNew.push_back(pkAfter);
2795
2796  // In case of a single recoiler, it just takes the remaining momentum.
2797  if (pRec.size() == 1) {
2798  pRec[0] = pajkAfter;
2799
2800  // Perform Lorentz boost for multiple recoilers, fails for a
2801  // single massless recoiler
2802  } else {
2803  for(int i=0; i<(int)pRec.size(); i++) {
2804  pRec[i].bstback(pAKBefore);
2805  pRec[i].bst(pajkAfter);
2806  }
2807  }
2808  return true;
2809
2810 }
2811
2812 //--------------------------------------------------------------------------
2813
2814 // Check if 2-particle system is on-shell and rescale if not.
2815
2816 bool VinciaCommon::onShellCM(Vec4& p1, Vec4& p2, double m1, double m2,
2817  double tol) {
2818
2819  double s1 = pow2(m1);
2820  double s2 = pow2(m2);
2821  double s01 = Vec4(p1+p2).m2Calc();
2822  double s1Calc = p1.m2Calc();
2823  double s2Calc = p2.m2Calc();
2824  if (abs(s1Calc-s1)/s01 > tol || abs(s2Calc-s2)/s01 > tol) {
2825  if (verbose >= 3)
2826  printOut("VinClu::onShellCM","forcing particles on mass shell");
2827  RotBstMatrix M;
2828  M.fromCMframe(p1,p2);
2829
2830  // Define massive on-shell momenta.
2831  double E0 = (s01 + s1 - s2)/(2*sqrt(s01));
2832  double E1 = (s01 - s1 + s2)/(2*sqrt(s01));
2833  double pz = pow2(E0)-s1;
2834  Vec4 p1new = Vec4(0.0,0.0,-pz,E0);
2835  Vec4 p2new = Vec4(0.0,0.0,pz,E1);
2836  p1new.rotbst(M);
2837  p2new.rotbst(M);
2838  double s1Test = p1new.m2Calc();
2839  double s2Test = p2new.m2Calc();
2840  if (verbose >= 3) {
2841  cout << " p1 : " << p1 << " p1new: " << p1new
2842  << " p2 : " << p1 << " p2new: " << p1new;
2843  }
2844
2845  // If this got them closer to mass shell, replace momenta.
2846  if (abs(s1Test-s1)/s01 <= abs(s1Calc-s1)/s01
2847  && abs(s2Test-s2)/s01 <= abs(s2Calc-s2)/s01) {
2848  p1 = p1new;
2849  p2 = p2new;
2850  }
2851  return false;
2852  }
2853  else return true;
2854
2855 }
2856
2857 //--------------------------------------------------------------------------
2858
2859 // Map partons partonSystems[iSys] to equivalent massless ones.
2860 // Return true if succeeded. Note, a method using only Particles or
2861 // Vec4 as input could in principle be split off from this, if needed,
2862 // but has not been required so far.
2863
2864 bool VinciaCommon::mapToMassless(int iSys, Event& event,
2865  PartonSystems* partonSystemsPtr, bool makeNewCopies) {
2866
2867  // Start by making new copies, if requested to do so.
2868  if (makeNewCopies) {
2869  int iOld, iNew;
2870
2871  // Copy incoming partons, interpret the copying operation as the
2872  // two incoming partons recoiling off each other. Assign status
2873  // code -42, incoming copy of recoiler (as mother).
2874  if (partonSystemsPtr->hasInAB(iSys)) {
2875  iOld = partonSystemsPtr->getInA(iSys);
2876  iNew = event.copy(iOld,-42);
2877  partonSystemsPtr->replace(iSys, iOld, iNew);
2878  iOld = partonSystemsPtr->getInB(iSys);
2879  iNew = event.copy(iOld,-42);
2880  partonSystemsPtr->replace(iSys, iOld, iNew);
2881  }
2882  // Note, a decaying resonance is not copied to preserve structure
2883  // of production and decay. Copy outgoing partons (use status code
2884  // 52).
2885  for (int i=0; i<partonSystemsPtr->sizeOut(iSys); ++i) {
2886  iOld = partonSystemsPtr->getOut(iSys,i);
2887  iNew = event.copy(iOld, 52);
2888  } // End loop to make new copies.
2889  } // End if new copies requested.
2890
2891  // Initial-state partons, always assumed massless in VINCIA.
2892  if (partonSystemsPtr->hasInAB(iSys) ) {
2893  int iA = partonSystemsPtr->getInA(iSys);
2894  int iB = partonSystemsPtr->getInB(iSys);
2895  if (event[iA].mCalc() != 0.0 || event[iB].mCalc() != 0.0) {
2896
2897  // Below we assume iA is the one with pz > 0; swap if opposite case.
2898  if (event[iA].pz() < 0 || event[iB].pz() > 0) {
2899  iA = partonSystemsPtr->getInB(iSys);
2900  iB = partonSystemsPtr->getInA(iSys);
2901  }
2902
2903  // Transverse components assumed zero: check.
2904  if (event[iA].pT() > 1.e-6 || event[iB].pT() > 1.e-6) {
2905  cout<<"Vincia::VinciaCommon::mapToMassless) Error: incoming partons "
2906  "have nonvanishing transverse momenta for system iSys = "
2907  <<iSys<<"; giving up"<<endl;
2908  return false;
2909  }
2910
2911  // Verbose output.
2912  if (verbose > superdebug) {
2913  stringstream ss;
2914  ss<<"Warning: forcing initial"
2915  "-state partons to be massless for system "<<iSys;
2916  printOut(__METHOD_NAME__,ss.str());
2917  }
2918
2919  // Define explicitly massless momenta (same as in Pythia::PartonLevel).
2920  double pPos = event[iA].pPos() + event[iB].pPos();
2921  double pNeg = event[iA].pNeg() + event[iB].pNeg();
2922  event[iA].pz( 0.5 * pPos);
2923  event[iA].e ( 0.5 * pPos);
2924  event[iA].m ( 0.);
2925  event[iB].pz(-0.5 * pNeg);
2926  event[iB].e ( 0.5 * pNeg);
2927  event[iB].m ( 0.);
2928  }
2929  } // End make initial-state partons massless.
2930
2931  // Final-state partons.
2932  if (partonSystemsPtr->sizeOut(iSys) >= 2) {
2933  vector<Vec4> momenta;
2934  vector<double> massOrg;
2935  bool makeMassless = false;
2936  Vec4 pSysOrg;
2937  for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
2938  momenta.push_back(event[partonSystemsPtr->getOut(iSys,i)].p());
2939  massOrg.push_back(event[partonSystemsPtr->getOut(iSys,i)].m());
2940  if (massOrg[i] > 0. && event[partonSystemsPtr->getOut(iSys,i)].idAbs()
2941  <= nFlavZeroMass) makeMassless = true;
2942  pSysOrg += momenta[i];
2943  }
2944  // Return if nothing needs to be done.
2945  if (!makeMassless) return true;
2946
2947  // Create copy of original momenta (in case of failure).
2948  vector<Vec4> momentaOrg = momenta;
2949
2950  // Boost to CM if original system not at rest.
2951  double sCM = m2(pSysOrg);
2952  bool isInCM = ( pow2(pSysOrg.pAbs())/sCM < 1e-10 );
2953  if (!isInCM)
2954  for (int i=0; i<(int)momenta.size(); ++i) momenta[i].bstback(pSysOrg);
2955
2956  // Define vector for computing CM energy of modified system.
2957  Vec4 pSysNew;
2958
2959  // Identify particles to be made massless (by ID code) and rescale
2960  // their momenta along direction of motion.
2961  for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
2962  int ip = partonSystemsPtr->getOut(iSys,i);
2963  if (event[ip].idAbs() <= nFlavZeroMass && event[ip].m() != 0.) {
2964  double facInv = momenta[i].pAbs()/momenta[i].e();
2965  // Sanity check.
2966  if (facInv <= 0.) {
2967  if (verbose >= 2)
2968  printOut("VinciaCommon:mapToMassless",
2969  "Remap failed. Particle is spacelike or at rest.");
2970  // Restore masses in case any were already changed.
2971  for (int j=0; j < partonSystemsPtr->sizeOut(iSys); ++j)
2972  event[partonSystemsPtr->getOut(iSys,j)].m(massOrg[j]);
2973  // Failed.
2974  return false;
2975  }
2976  momenta[i].rescale3(1./facInv);
2977  event[ip].m(0.);
2978  // Check new 4-vector.
2979  double mNew = momenta[i].mCalc();
2980  if (pow2(mNew/momenta[i].e()) > TINY) {
2981  printOut("VinciaCommon:mapToMassless","Warning: rounding problem.");
2982  if (verbose >= 7) {
2983  cout<<scientific << "(p,e) = "<<momenta[i].pAbs() << " "
2984  << momenta[i].e() << " facInv = " << facInv
2985  << " 1/facInv = " << 1./facInv << " mNew = " << mNew << endl;
2986  }
2987  } // End check new 4-vector.
2988  } // End if massless flavour with mass > 0.
2989
2990  // Add to new CM momentum.
2991  pSysNew += momenta[i];
2992  } // End loop over FS particles.
2993
2994  // New system generally has smaller invariant mass and some
2995  // motion. Determine if additional scalings or boosts are needed.
2996  Vec4 delta = pSysOrg - pSysNew;
2997  if (delta.e()/sqrt(sCM) < TINY && delta.pAbs()/sqrt(sCM) < TINY) {
2998  // Update event record (masses already updated above).
2999  for (int i = 0; i < (int)momenta.size(); ++i)
3000  event[partonSystemsPtr->getOut(iSys,i)].p(momenta[i]);
3001  if (verbose >= superdebug)
3002  printOut("VinciaCommon:mapToMassless","No further rescalings needed.");
3003
3004  //Check momentum conservation.
3005  if (!checkCoM(iSys,event,partonSystemsPtr)) {
3006  infoPtr->setAbortPartonLevel(true);
3007  infoPtr->errorMsg("Error in "+__METHOD_NAME__+
3008  ": Failed (E,p) conservation check.","Aborting.");
3009  return false;
3010  }
3011  return true;
3012  }
3013
3014  // If the new system has a different CM energy, rescale all
3015  // energies and momenta to restore the same CM energy as before.
3016  double sCMnew = m2(pSysNew);
3017  double scaleFac = sqrt(sCM/sCMnew);
3018  if (verbose >= 7 && pow2(scaleFac-1.0) > TINY)
3019  printOut("VinciaCommon:mapToMassless",
3020  "Rescaling 4-vectors to restore eCM");
3021  Vec4 pSysNewB;
3022  for (int i = 0; i < (int)momenta.size(); ++i) {
3023  momenta[i].rescale4(scaleFac);
3024  pSysNewB += momenta[i];
3025  }
3026  double sCMnewB = m2(pSysNewB);
3027  if (verbose >= superdebug)
3028  cout << "old CM energy = " << sqrt(sCM) << " intermediate CM energy = "
3029  << sqrt(sCMnew) << " new CM energy = " << sqrt(sCMnewB) << endl;
3030  // Then boost to CM frame (preserves CM energy).
3031  if (verbose >= 7)
3032  printOut("VinciaCommon:mapToMassless","Boosting back to CM frame");
3033  for (int i=0; i<(int)momenta.size(); ++i) {
3034  // Boost to new CM frame
3035  momenta[i].bstback(pSysNewB);
3036  // If required, also boost back to frame of input system
3037  if (!isInCM) momenta[i].bst(pSysOrg);
3038
3039  // Update event record (masses already updated above).
3040  event[partonSystemsPtr->getOut(iSys,i)].p(momenta[i]);
3041
3042  } // End do boosts.
3043
3044  // Verbose output: final configuration.
3045  if (verbose >= 7) {
3046  cout << "Final configuration:" << endl;
3047  for (int i = 0; i < (int)momenta.size(); ++i)
3048  cout << " " << i << " " << momenta[i];
3049  }
3050
3051  } // End make final-state momenta massless.
3052
3053  // Check momentum conservation.
3054  if(!checkCoM(iSys, event,partonSystemsPtr)){
3055  infoPtr->setAbortPartonLevel(true);
3056  infoPtr->errorMsg("Error in "+__METHOD_NAME__
3057  +": Failed (E,p) conservation check.","Aborting.");
3058  return false;
3059  }
3060
3062  return true;
3063
3064 }
3065
3066 //--------------------------------------------------------------------------
3067
3068 // More lightweight function to check conservation of momentum.
3069
3070 bool VinciaCommon::checkCoM(int iSys, Event& event,
3071  PartonSystems* partonSystemsPtr){
3072  Vec4 total(0.,0.,0.,0.);
3073  if (!partonSystemsPtr->hasInRes(iSys)){
3074  if (partonSystemsPtr->getInA(iSys) > 0)
3075  total+= event[partonSystemsPtr->getInA(iSys)].p();
3076  if (partonSystemsPtr->getInB(iSys) > 0)
3077  total+= event[partonSystemsPtr->getInB(iSys)].p();
3078  } else total+= event[partonSystemsPtr->getInRes(iSys)].p();
3079  double sysMass = total.mCalc();
3080
3081  // Loop over members of current system.
3082  for( int iPart=0; iPart<partonSystemsPtr->sizeOut(iSys); iPart++){
3083  int iOut = partonSystemsPtr->getOut(iSys,iPart);
3084
3085  // Sum total FS momentum.
3086  if (event[iOut].isFinal()) total -= event[iOut].p();
3087  else {
3088  stringstream ss;
3089  ss << "iSys = " << iSys << " iOut = " << iOut;
3090  infoPtr->errorMsg("Error in "+__METHOD_NAME__
3091  +": isFinal()=false for outgoing parton.",ss.str());
3092  partonSystemsPtr->list();
3093  event.list();
3094  return false;
3095  }
3096  }
3097  total/=sysMass;
3098  if(abs(total.e()) > SMALL || abs(total.px()) > SMALL
3099  || abs(total.py()) > SMALL || abs(total.pz()) > SMALL) {
3100  event.list();
3101  cout << "total = " << setprecision(10) << total.e() << " " << total.px()
3102  << " " << total.py() << " " << total.pz() << endl;
3103  infoPtr->errorMsg("Error in "+__METHOD_NAME__
3104  +" Failed (E,p) conservation check.");
3105  return false;
3106  } else if(isnan(total)){
3107  event.list();
3108  infoPtr->errorMsg("Error in "+__METHOD_NAME__
3109  +" Failed (E,p) isnan check.");
3110  return false;
3111  } else return true;
3112
3113 }
3114
3115 //==========================================================================
3116
3117 // TXiFunctor helper class.
3118
3119 //--------------------------------------------------------------------------
3120
3121 // Constructor.
3122
3123 TXiFunctor::TXiFunctor(vector<double> mIn, vector<double> energiesIn) {
3124  if (mIn.size() != energiesIn.size()) {
3125  m = vector<double>(0); energies = vector<double>(0);
3126  } else {m = mIn; energies = energiesIn;}
3127 }
3128
3129 double TXiFunctor::operator() (double xi) {
3130  double retval = 0.;
3131  for (vector<double>::size_type i = 0; i < m.size(); i++)
3132  retval += sqrt( pow2(m[i]) + pow2(xi)*pow2(energies[i]));
3133  return retval;
3134 }
3135
3136 //==========================================================================
3137
3138 // VINCIA Auxiliary helper functions.
3139
3140 //--------------------------------------------------------------------------
3141
3142 // External auxiliaries, extra four-products.
3143
3144 double m(const Vec4& v) {
3145  double s = m2(v); return (s >= 0) ? sqrt(s) : -sqrt(-s);}
3146
3147 double m2(const Vec4& v) {
3148  return pow2(v.e()) - pow2(v.px()) - pow2(v.py()) - pow2(v.pz());}
3149
3150 double m2(const Vec4& v1, const Vec4& v2, const Vec4& v3) {
3151  return pow2(v1.e() + v2.e() + v3.e())
3152  - pow2(v1.px() + v2.px() + v3.px())
3153  - pow2(v1.py() + v2.py() + v3.py())
3154  - pow2(v1.pz() + v2.pz() + v3.pz());}
3155
3156 double m2(const Vec4& v1, const Vec4& v2, const Vec4& v3, const Vec4& v4) {
3157  return pow2(v1.e() + v2.e() + v3.e() + v4.e())
3158  - pow2(v1.px() + v2.px() + v3.px() + v4.px())
3159  - pow2(v1.py() + v2.py() + v3.py() + v4.py())
3160  - pow2(v1.pz() + v2.pz() + v3.pz() + v4.pz());}
3161
3162 double m2(const Particle& p1, const Particle& p2, const Particle& p3) {
3163  return m2(p1.p(), p2.p(), p3.p());}
3164
3165 double dot4(const Particle& p1, const Particle& p2) {return p1.p()*p2.p();}
3166
3167 double getCosTheta(double E1, double E2, double m1, double m2, double s12){
3168  return (2.0*E1*E2 - s12)/(2.0*sqrt(E1*E1 - m1*m1)*sqrt(E2*E2 - m2*m2));}
3169
3170 //--------------------------------------------------------------------------
3171
3172 // External auxiliaries, string manipulation.
3173
3174 string num2str(int i, int width) {
3175  ostringstream tmp;
3176  if (width <= 1) tmp << i;
3177  else if (abs(i) < pow(10.0, width - 1) || ( i > 0 && i < pow(10.0, width)))
3178  tmp << fixed << setw(width) << i;
3179  else {
3180  string ab = "k";
3181  double r = i;
3182  if (abs(i) < 1e5) {r/=1e3;}
3183  else if (abs(i) < 1e8) {r/=1e6; ab = "M";}
3184  else if (abs(i) < 1e11) {r/=1e9; ab = "G";}
3185  else if (abs(i) < 1e14) {r/=1e12; ab = "T";}
3186  tmp << fixed << setw(width - 1)
3187  << (r > 10 ? setprecision(width-4) : setprecision(width-3)) << r << ab;
3188  }
3189  return tmp.str();
3190 }
3191
3192 string num2str(double r, int width) {
3193  ostringstream tmp;
3194  if (width <= 0) tmp << r;
3195  else if (r == 0.0 || (abs(r) > 0.1 && abs(r) < pow(10., max(width-3,1)))
3196  || width <= 8) tmp << fixed << setw(max(width,3))
3197  << setprecision(min(3, max(1, width - 2))) << r;
3198  else tmp << scientific << setprecision(max(2, width - 7))
3199  << setw(max(9, width)) << r;
3200  return tmp.str();
3201 }
3202
3203 string bool2str(bool b, int width) {
3204  string tmp = b ? "on" : "off";
3205  int nPad = width - tmp.length();
3206  for (int i = 1; i <= nPad; ++i) tmp = " " + tmp;
3207  return tmp;
3208 }
3209
3210 void printOut(string place, string message) {
3211  cout.setf(ios::internal);
3212  cout << " (" << (place + ") ") << message << "\n";
3213 }
3214
3215 //--------------------------------------------------------------------------
3216
3217 // Gram determinant, invariants used in the argument = 2*pi*pj.
3218
3219 double gramDet( double s01tilde, double s12tilde, double s02tilde,
3220  double m0, double m1, double m2) {
3221  return ((s01tilde*s12tilde*s02tilde - pow2(s01tilde)*pow2(m2)
3222  - pow2(s02tilde)*pow2(m1) - pow2(s12tilde)*pow2(m0))/4
3223  + pow2(m0)*pow2(m1)*pow2(m2));
3224 }
3225
3226 double gramDet(Vec4 p0, Vec4 p1, Vec4 p2) {
3227  return gramDet(2*p0*p1, 2*p1*p2, 2*p0*p2, p0.mCalc(), p1.mCalc(),
3228  p2.mCalc());
3229 }
3230
3231 //--------------------------------------------------------------------------
3232
3233 // Math support auxiliaries.
3234
3235 // Dilogarithm.
3236 double Li2(const double x, const double kmax, const double xerr) {
3237  if (x < 0.0) return 0.5*Li2(x*x) - Li2(-x);
3238  if (x <= 0.5) {
3239  double sum(x), term(x);
3240  for (int k = 2; k < kmax; k++) {
3241  double rk = (k-1.0)/k;
3242  term *= x*rk*rk;
3243  sum += term;
3244  if (abs(term/sum) < xerr) return sum;
3245  }
3246  cout << "Maximum number of iterations exceeded in Li2" << endl;
3247  return sum;
3248  }
3249  if (x < 1.0) return M_PI*M_PI/6.0 - Li2(1.0 - x) - log(x)*log(1.0 - x);
3250  if (x == 1.0) return M_PI*M_PI/6.0;
3251  if (x <= 1.01) {
3252  const double eps(x - 1.0), lne(log(eps)),
3253  c0(M_PI*M_PI/6.0), c1( 1.0 - lne),
3254  c2(-(1.0 - 2.0*lne)/4.0), c3( (1.0 - 3.0*lne)/9.0),
3255  c4(-(1.0 - 4.0*lne)/16.0), c5( (1.0 - 5.0*lne)/25.0),
3256  c6(-(1.0 - 6.0*lne)/36.0), c7( (1.0 - 7.0*lne)/49.0),
3257  c8(-(1.0 - 8.0*lne)/64.0);
3258  return c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*(c5 + eps*(
3259  c6 + eps*(c7 + eps*c8)))))));
3260  }
3261  double logx = log(x);
3262  if (x<=2.0) return M_PI*M_PI/6.0 + Li2(1.0 - 1.0/x) -
3263  logx*(log(1.0 - 1.0/x) + 0.5*logx);
3264  return M_PI*M_PI/3.0 - Li2(1.0/x) - 0.5*logx*logx;
3265 }
3266
3267
3268 // Standard factorial.
3269 double factorial(const int n) {
3270  double fac = 1;
3271  for (int i = 2; i <= n; i++) fac *= i;
3272  return fac;}
3273
3274 // Binomial coefficient.
3275 int binomial(const int n, const int m) {
3276  if (m < 0 || m > n) return 0;
3277  else if (m == n || m == 0) return 1;
3278  else if (m == 1 || m == n - 1) return n;
3279  else return factorial(n)/factorial(m)/factorial(n - m) + 0.01;
3280 }
3281
3282 // Lambert W function using the rational fit from Darko Veberic's
3283 // paper, arXiv:1209.0735v2. Should give 5 digits of precision for
3284 // positive arguments x not too large (fit region was 0.3, 2e, but
3285 // still has 5-digit accuracy at zero). Precision quickly drops for
3286 // negative values, but he has extra functions that can be implemented
3287 // if those are needed, and for very large values the asymptotic
3288 // log(x), log(log(x)) form could be used if precise solutions for
3289 // large values are needed. For now just write a warning if we are
3290 // ever asked for a value far outside region of validity.
3291 double LambertW(const double x) {
3292  if (x == 0.) return 0.;
3293  if (x < -0.2) {
3294  cout << "Warning in "<<__METHOD_NAME__
3295  << ": Accuracy less than three decimal places for x < -0.2";
3296  } else if (x > 10.) {
3297  cout << "Warning in "<<__METHOD_NAME__
3298  <<": Accuracy less than three decimal places for x > 10.";
3299  }
3300  return x*(1. + x*(2.445053 + x*(1.343664 + x*(0.14844 + 0.000804*x))))
3301  /(1. + x*(3.444708 + x*(3.292489 + x*(0.916460 + x*(0.053068)))));
3302 }
3303
3304 // Version of zbrent using function pointers, solve fun(x) - r = 0 for x.
3305 double zbrent(TFunctor& fun, double r, double x1, double x2, double tol) {
3306  int iter;
3307  double a(x1), b(x2), c(x2), d(x2-x1), e(x2-x1), min1, min2;
3308  double fa(fun(a) - r), fb(fun(b) - r), fc, p, q, r1, s, tol1, xm;
3309  double REALTINY = min(TINY, 1e-12);
3310  tol = max(tol, REALTINY);
3311
3312  // Check if there is a single zero in range.
3313  if (fa*fb > 0) return 0.0;
3314
3315  // Start search.
3316  fc = fb;
3317  for (iter = 1; iter < max(1000, int(1.0/sqrt(tol))); iter++) {
3318  if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
3319  c = a; fc = fa; e = d = b-a;
3320  }
3321  if (abs(fc) < abs(fb)) {
3322  a = b; b = c; c= a; fa = fb; fb = fc; fc = fa;
3323  }
3324  tol1 = 2.0*REALTINY*abs(b) + 0.5*tol;
3325  xm = 0.5*(c-b);
3326  if (abs(xm) <= tol1 || fb == 0.0) return b;
3327  if (abs(e) >= tol1 && abs(fa) > abs(fb)) {
3328  s = fb/fa;
3329  if (a == c) {p = 2.0*xm*s; q = 1.0-s;}
3330  else {
3331  q = fa/fc; r1 = fb/fc;
3332  p = s*(2.0*xm*q*(q - r1) - (b - a)*(r1 - 1.0));
3333  q = (q - 1.0)*(r1 - 1.0)*(s - 1.0);
3334  }
3335  if (p > 0.0) q = -q;
3336  p = abs(p);
3337  min1 = 3.0*xm*q - abs(tol1*q);
3338  min2 = abs(e*q);
3339  if (2.0*p < (min1 < min2 ? min1 : min2)) {e = d; d= p/q;}
3340  else {d = xm; e = d;}
3341  } else {d = xm; e = d;}
3342  a = b;
3343  fa = fb;
3344  b += abs(d) > tol1 ? d : (xm > 1) ? tol1 : - tol1;
3345  fb = fun(b) - r;
3346  }
3347  cerr << "(brent:) -> Maximum number of iterations exceeded" << endl;
3348  return 0.0;
3349
3350 }
3351
3352 //==========================================================================
3353
3354 } // end namespace Pythia8
Definition: AgUStep.h:26