StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
PartonLevel.cc
1 // PartonLevel.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) for the PartonLevel class.
7 
8 #include "Pythia8/PartonLevel.h"
9 
10 namespace Pythia8 {
11 
12 //==========================================================================
13 
14 // The PartonLevel class.
15 
16 //--------------------------------------------------------------------------
17 
18 // Constants: could be changed here if desired, but normally should not.
19 // These are of technical nature, as described for each.
20 
21 // Maximum number of tries to produce parton level from given input.
22 const int PartonLevel::NTRY = 10;
23 
24 //--------------------------------------------------------------------------
25 
26 // Main routine to initialize the parton-level generation process.
27 
28 bool PartonLevel::init( Info* infoPtrIn, Settings& settings,
29  ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
30  BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
31  BeamParticle* beamPomAPtrIn, BeamParticle* beamPomBPtrIn,
32  Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn,
33  SigmaTotal* sigmaTotPtr, TimeShower* timesDecPtrIn, TimeShower* timesPtrIn,
34  SpaceShower* spacePtrIn, RHadrons* rHadronsPtrIn, UserHooks* userHooksPtrIn,
35  MergingHooks* mergingHooksPtrIn, bool useAsTrial ) {
36 
37  // Store input pointers and modes for future use.
38  infoPtr = infoPtrIn;
39  particleDataPtr = particleDataPtrIn;
40  rndmPtr = rndmPtrIn;
41  beamAPtr = beamAPtrIn;
42  beamBPtr = beamBPtrIn;
43  beamHadAPtr = beamAPtr;
44  beamHadBPtr = beamBPtr;
45  beamPomAPtr = beamPomAPtrIn;
46  beamPomBPtr = beamPomBPtrIn;
47  couplingsPtr = couplingsPtrIn;
48  partonSystemsPtr = partonSystemsPtrIn;
49  timesDecPtr = timesDecPtrIn;
50  timesPtr = timesPtrIn;
51  spacePtr = spacePtrIn;
52  rHadronsPtr = rHadronsPtrIn;
53  userHooksPtr = userHooksPtrIn;
54  mergingHooksPtr = mergingHooksPtrIn;
55 
56  // Min bias and diffraction processes need special treatment.
57  bool doSQ = settings.flag("SoftQCD:all")
58  || settings.flag("SoftQCD:inelastic");
59  bool doND = settings.flag("SoftQCD:nonDiffractive");
60  bool doSD = settings.flag("SoftQCD:singleDiffractive");
61  bool doDD = settings.flag("SoftQCD:doubleDiffractive");
62  bool doCD = settings.flag("SoftQCD:centralDiffractive");
63  doNonDiff = doSQ || doND;
64  doDiffraction = doSQ || doSD || doDD || doCD;
65 
66  // Separate low-mass (unresolved) and high-mass (perturbative) diffraction.
67  mMinDiff = settings.parm("Diffraction:mMinPert");
68  mWidthDiff = settings.parm("Diffraction:mWidthPert");
69  pMaxDiff = settings.parm("Diffraction:probMaxPert");
70  if (mMinDiff > infoPtr->eCM()) doDiffraction = false;
71 
72  // Need MPI initialization for soft QCD processes, even if only first MPI.
73  // But no need to initialize MPI if never going to use it.
74  doMPI = settings.flag("PartonLevel:MPI");
75  doMPIMB = doMPI;
76  doMPISDA = doMPI;
77  doMPISDB = doMPI;
78  doMPICD = doMPI;
79  doMPIinit = doMPI;
80  if (doNonDiff || doDiffraction) doMPIinit = true;
81  if (!settings.flag("PartonLevel:all")) doMPIinit = false;
82 
83  // Initialise trial shower switch.
84  doTrial = useAsTrial;
85  // Merging initialization.
86  bool hasMergingHooks = (mergingHooksPtr != 0);
87  canRemoveEvent = !doTrial && hasMergingHooks
88  && ( mergingHooksPtr->doCKKWLMerging() || mergingHooksPtr->doNL3Merging());
89  canRemoveEmission = !doTrial && hasMergingHooks
90  && ( mergingHooksPtr->doUMEPSMerging() || mergingHooksPtr->doNL3Merging()
91  || mergingHooksPtr->doUNLOPSMerging() );
92  nTrialEmissions = 1;
93  pTLastBranch = 0.0;
94  typeLastBranch = 0;
95 
96  // Flags for showers: ISR and FSR.
97  doISR = settings.flag("PartonLevel:ISR");
98  bool FSR = settings.flag("PartonLevel:FSR");
99  bool FSRinProcess = settings.flag("PartonLevel:FSRinProcess");
100  bool interleaveFSR = settings.flag("TimeShower:interleave");
101  doFSRduringProcess = FSR && FSRinProcess && interleaveFSR;
102  doFSRafterProcess = FSR && FSRinProcess && !interleaveFSR;
103  doFSRinResonances = FSR && settings.flag("PartonLevel:FSRinResonances");
104 
105  // Some other flags.
106  doRemnants = settings.flag("PartonLevel:Remnants");
107  doSecondHard = settings.flag("SecondHard:generate");
108  earlyResDec = settings.flag("PartonLevel:earlyResDec");
109 
110  // Allow R-hadron formation.
111  allowRH = settings.flag("RHadrons:allow");
112 
113  // Possibility to allow user veto during evolution.
114  canVetoPT = (userHooksPtr != 0)
115  ? userHooksPtr->canVetoPT() : false;
116  pTvetoPT = (canVetoPT)
117  ? userHooksPtr->scaleVetoPT() : -1.;
118  canVetoStep = (userHooksPtr != 0)
119  ? userHooksPtr->canVetoStep() : false;
120  nVetoStep = (canVetoStep)
121  ? userHooksPtr->numberVetoStep() : -1;
122  canVetoMPIStep = (userHooksPtr != 0)
123  ? userHooksPtr->canVetoMPIStep() : false;
124  nVetoMPIStep = (canVetoMPIStep)
125  ? userHooksPtr->numberVetoMPIStep() : -1;
126  canVetoEarly = (userHooksPtr != 0)
127  ? userHooksPtr->canVetoPartonLevelEarly() : false;
128 
129  // Settings for vetoing of QCD emission for Drell-Yan weak boson production.
130  vetoWeakJets = settings.flag("WeakShower:vetoQCDjets");
131  vetoWeakDeltaR2 = pow2(settings.parm("WeakShower:vetoWeakDeltaR"));
132 
133  // Possibility to set maximal shower scale in resonance decays.
134  canSetScale = (userHooksPtr != 0)
135  ? userHooksPtr->canSetResonanceScale() : false;
136 
137  // Possibility to reconnect specifically for resonance decays.
138  canReconResSys = (userHooksPtr != 0)
139  ? userHooksPtr->canReconnectResonanceSystems() : false;
140 
141  // Done with initialization only for FSR in resonance decays.
142  if (beamAPtr == 0 || beamBPtr == 0) return true;
143 
144  // Flag if lepton beams, and if non-resolved ones. May change main flags.
145  hasLeptonBeams = ( beamAPtr->isLepton() || beamBPtr->isLepton() );
146  hasPointLeptons = ( hasLeptonBeams
147  && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) );
148  if (hasLeptonBeams) {
149  doMPIMB = false;
150  doMPISDA = false;
151  doMPISDB = false;
152  doMPICD = false;
153  doMPIinit = false;
154  }
155  if (hasPointLeptons) {
156  doISR = false;
157  doRemnants = false;
158  }
159 
160  // Set info and initialize the respective program elements.
161  timesPtr->init( beamAPtr, beamBPtr);
162  if (doISR) spacePtr->init( beamAPtr, beamBPtr);
163  doMPIMB = multiMB.init( doMPIinit, 0, infoPtr, settings, particleDataPtr,
164  rndmPtr, beamAPtr, beamBPtr, couplingsPtr, partonSystemsPtr, sigmaTotPtr,
165  userHooksPtr);
166  if (doSD || doDD || doSQ) doMPISDA = multiSDA.init( doMPIinit, 1, infoPtr,
167  settings, particleDataPtr, rndmPtr, beamAPtr, beamPomBPtr, couplingsPtr,
168  partonSystemsPtr, sigmaTotPtr, userHooksPtr);
169  if (doSD || doDD || doSQ) doMPISDB = multiSDB.init( doMPIinit, 2, infoPtr,
170  settings, particleDataPtr, rndmPtr, beamPomAPtr, beamBPtr, couplingsPtr,
171  partonSystemsPtr, sigmaTotPtr, userHooksPtr);
172  if (doCD || doSQ) doMPICD = multiCD.init( doMPIinit, 3, infoPtr, settings,
173  particleDataPtr, rndmPtr, beamPomAPtr, beamPomBPtr, couplingsPtr,
174  partonSystemsPtr, sigmaTotPtr, userHooksPtr);
175  remnants.init( infoPtr, settings, rndmPtr, beamAPtr, beamBPtr,
176  partonSystemsPtr);
177  resonanceDecays.init( infoPtr, particleDataPtr, rndmPtr);
178 
179  // Succeeded, or not.
180  multiPtr = &multiMB;
181  if (doMPIinit && !doMPIMB) return false;
182  if (doMPIinit && (doSD || doDD || doSQ) && (!doMPISDA || !doMPISDB))
183  return false;
184  if (doMPIinit && (doCD || doSQ) && !doMPICD) return false;
185  if (!doMPIMB || !doMPISDA || !doMPISDB || !doMPICD) doMPI = false;
186  return true;
187 }
188 
189 //--------------------------------------------------------------------------
190 
191 // Function to reset PartonLevel object for trial shower usage.
192 
193 void PartonLevel::resetTrial() {
194 
195  // Clear input pointers.
196  partonSystemsPtr->clear();
197  beamAPtr->clear();
198  beamBPtr->clear();
199  beamHadAPtr->clear();
200  beamHadBPtr->clear();
201  beamPomAPtr->clear();
202  beamPomBPtr->clear();
203 
204  // Clear last branching return values.
205  pTLastBranch = 0.0;
206  typeLastBranch = 0;
207 
208 }
209 
210 //--------------------------------------------------------------------------
211 
212 // Main routine to do the parton-level evolution.
213 
214 bool PartonLevel::next( Event& process, Event& event) {
215 
216  // Current event classification.
217  isResolved = infoPtr->isResolved();
218  isResolvedA = isResolved;
219  isResolvedB = isResolved;
220  isResolvedC = isResolved;
221  isDiffA = infoPtr->isDiffractiveA();
222  isDiffB = infoPtr->isDiffractiveB();
223  isDiffC = infoPtr->isDiffractiveC();
224  isDiff = isDiffA || isDiffB || isDiffC;
225  isCentralDiff = isDiffC;
226  isDoubleDiff = isDiffA && isDiffB;
227  isSingleDiff = isDiff && !isDoubleDiff && !isCentralDiff;
228  isNonDiff = infoPtr->isNonDiffractive();
229  doVeto = false;
230 
231  // nHardLoop counts how many hard-scattering subsystems are to be processed.
232  // Almost always 1, but elastic and low-mass diffraction gives 0, while
233  // double diffraction can give up to 2. Not to be confused with SecondHard.
234  int nHardLoop = 1;
235  if (!isResolved) nHardLoop = (isDiff) ? decideResolvedDiff( process) : 0;
236 
237  // Handle unresolved subsystems. Done if no resolved ones.
238  sizeProcess = 0;
239  sizeEvent = 0;
240  if (!isResolvedA || !isResolvedB || !isResolvedC) {
241  bool physical = setupUnresolvedSys( process, event);
242  if (!physical || nHardLoop == 0) return physical;
243  sizeProcess = process.size();
244  sizeEvent = event.size();
245  }
246 
247  // Number of actual branchings.
248  int nBranch = 0;
249  // Number of desired branchings, negative value means no restriction.
250  int nBranchMax = (doTrial) ? nTrialEmissions : -1;
251 
252  // Store merging weight.
253  bool hasMergingHooks = (mergingHooksPtr != 0);
254  if ( hasMergingHooks && canRemoveEvent )
255  mergingHooksPtr->storeWeights(infoPtr->getWeightCKKWL());
256 
257  // Big outer loop to handle up to two systems (in double diffraction),
258  // but normally one. (Not indented in following, but end clearly marked.)
259  for (int iHardLoop = 1; iHardLoop <= nHardLoop; ++iHardLoop) {
260  infoPtr->setCounter(20, iHardLoop);
261  infoPtr->setCounter(21);
262 
263  // Classification of diffractive system: 1 = A, 2 = B, 3 = central.
264  iDS = 0;
265  if (isDiffA || isDiffB) iDS = (iHardLoop == 2 || !isResolvedA) ? 2 : 1;
266  if (isDiffC) iDS = 3;
267 
268  // Process and event records can be out of step for diffraction.
269  if (iHardLoop == 2) {
270  sizeProcess = process.size();
271  sizeEvent = event.size();
272  partonSystemsPtr->clear();
273  if (event.lastColTag() > process.lastColTag())
274  process.initColTag(event.lastColTag());
275  }
276 
277  // If you need to restore then do not throw existing diffractive system.
278  if (isDiff) {
279  event.saveSize();
280  event.saveJunctionSize();
281 
282  // Allow special treatment of diffractive systems.
283  setupResolvedDiff( process);
284  }
285 
286  // Prepare to do multiparton interactions; at new mass for diffraction.
287  if (doMPIinit) multiPtr->reset();
288 
289  // Special case if nondiffractive: do hardest interaction.
290  if (isNonDiff || isDiff) {
291  multiPtr->pTfirst();
292  multiPtr->setupFirstSys( process);
293  }
294 
295  // Allow up to ten tries; failure possible for beam remnants.
296  // Main cause: inconsistent colour flow at the end of the day.
297  bool physical = true;
298  int nRad = 0;
299  for (int iTry = 0; iTry < NTRY; ++ iTry) {
300  infoPtr->addCounter(21);
301  for (int i = 22; i < 32; ++i) infoPtr->setCounter(i);
302 
303  // Reset flag, counters and max scales.
304  physical = true;
305  nMPI = (doSecondHard) ? 2 : 1;
306  nISR = 0;
307  nFSRinProc = 0;
308  nFSRinRes = 0;
309  nISRhard = 0;
310  nFSRhard = 0;
311  pTsaveMPI = 0.;
312  pTsaveISR = 0.;
313  pTsaveFSR = 0.;
314 
315  // Identify hard interaction system for showers.
316  setupHardSys( process, event);
317 
318  // Optionally check for a veto after the hardest interaction.
319  if (canVetoMPIStep) {
320  doVeto = userHooksPtr->doVetoMPIStep( 1, event);
321  // Abort event if vetoed.
322  if (doVeto) {
323  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
324  return false;
325  }
326  }
327 
328  // Check matching of process scale to maximum ISR/FSR/MPI scales.
329  double Q2Fac = infoPtr->Q2Fac();
330  double Q2Ren = infoPtr->Q2Ren();
331  bool limitPTmaxISR = (doISR)
332  ? spacePtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
333  bool limitPTmaxFSR = (doFSRduringProcess)
334  ? timesPtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
335  bool limitPTmaxMPI = (doMPI) ? multiPtr->limitPTmax( event) : false;
336 
337  // Global recoil: reset counters and store locations of outgoing partons.
338  timesPtr->prepareGlobal( event);
339  bool isFirstTrial = true;
340 
341  // Set hard scale, maximum for showers and multiparton interactions.
342  double pTscaleRad = process.scale();
343  double pTscaleMPI = pTscaleRad;
344  if (doSecondHard) {
345  pTscaleRad = max( pTscaleRad, process.scaleSecond() );
346  pTscaleMPI = min( pTscaleMPI, process.scaleSecond() );
347  }
348  double pTmaxMPI = (limitPTmaxMPI) ? pTscaleMPI : infoPtr->eCM();
349  double pTmaxISR = (limitPTmaxISR) ? spacePtr->enhancePTmax() * pTscaleRad
350  : infoPtr->eCM();
351  double pTmaxFSR = (limitPTmaxFSR) ? timesPtr->enhancePTmax() * pTscaleRad
352  : infoPtr->eCM();
353 
354  // Potentially reset up starting scales for matrix element merging.
355  if ( hasMergingHooks && (doTrial || canRemoveEvent || canRemoveEmission) )
356  mergingHooksPtr->setShowerStartingScales( doTrial,
357  (canRemoveEvent || canRemoveEmission), pTscaleRad, process, pTmaxFSR,
358  limitPTmaxFSR, pTmaxISR, limitPTmaxISR, pTmaxMPI, limitPTmaxMPI );
359 
360  double pTmax = max( pTmaxMPI, max( pTmaxISR, pTmaxFSR) );
361  pTsaveMPI = pTmaxMPI;
362  pTsaveISR = pTmaxISR;
363  pTsaveFSR = pTmaxFSR;
364  // Prepare the classes to begin the generation.
365  if (doMPI) multiPtr->prepare( event, pTmaxMPI);
366  if (doISR) spacePtr->prepare( 0, event, limitPTmaxISR);
367  if (doFSRduringProcess) timesPtr->prepare( 0, event, limitPTmaxFSR);
368  if (doSecondHard && doISR) spacePtr->prepare( 1, event, limitPTmaxISR);
369  if (doSecondHard && doFSRduringProcess) timesPtr->prepare( 1, event,
370  limitPTmaxFSR);
371 
372  // Impact parameter has now been chosen, except for diffraction.
373  if (!isDiff) infoPtr->setImpact( multiPtr->bMPI(),
374  multiPtr->enhanceMPI(), true);
375  // Set up initial veto scale.
376  doVeto = false;
377  double pTveto = pTvetoPT;
378  typeLatest = 0;
379 
380  // Begin evolution down in pT from hard pT scale.
381  do {
382  infoPtr->addCounter(22);
383  typeVetoStep = 0;
384  nRad = nISR + nFSRinProc;
385 
386  // Find next pT value for FSR, MPI and ISR.
387  // Order calls to minimize time expenditure.
388  double pTgen = 0.;
389  double pTtimes = (doFSRduringProcess)
390  ? timesPtr->pTnext( event, pTmaxFSR, pTgen, isFirstTrial) : -1.;
391  pTgen = max( pTgen, pTtimes);
392  double pTmulti = (doMPI)
393  ? multiPtr->pTnext( pTmaxMPI, pTgen, event) : -1.;
394  pTgen = max( pTgen, pTmulti);
395  double pTspace = (doISR)
396  ? spacePtr->pTnext( event, pTmaxISR, pTgen, nRad) : -1.;
397  double pTnow = max( pTtimes, max( pTmulti, pTspace));
398 
399  // Update information.
400  infoPtr->setPTnow( pTnow);
401  isFirstTrial = false;
402 
403  // Allow a user veto. Only do it once, so remember to change pTveto.
404  if (pTveto > 0. && pTveto > pTnow) {
405  pTveto = -1.;
406  doVeto = userHooksPtr->doVetoPT( typeLatest, event);
407  // Abort event if vetoed.
408  if (doVeto) {
409  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
410  return false;
411  }
412  }
413 
414  // Do a multiparton interaction (if allowed).
415  if (pTmulti > 0. && pTmulti > pTspace && pTmulti > pTtimes) {
416  infoPtr->addCounter(23);
417  if (multiPtr->scatter( event)) {
418  typeLatest = 1;
419  ++nMPI;
420  if (canVetoMPIStep && nMPI <= nVetoMPIStep) typeVetoStep = 1;
421 
422  // Update ISR and FSR dipoles.
423  if (doISR) spacePtr->prepare( nMPI - 1, event);
424  if (doFSRduringProcess) timesPtr->prepare( nMPI - 1, event);
425  }
426 
427  // Set maximal scales for next pT to pick.
428  pTmaxMPI = pTmulti;
429  pTmaxISR = min(pTmulti, pTmaxISR);
430  pTmaxFSR = min(pTmulti, pTmaxFSR);
431  pTmax = pTmulti;
432  nBranch++;
433  pTLastBranch = pTmulti;
434  typeLastBranch = 1;
435  }
436 
437  // Do an initial-state emission (if allowed).
438  else if (pTspace > 0. && pTspace > pTtimes) {
439  infoPtr->addCounter(24);
440  if (spacePtr->branch( event)) {
441  typeLatest = 2;
442  iSysNow = spacePtr->system();
443  ++nISR;
444  if (iSysNow == 0) ++nISRhard;
445  if (canVetoStep && iSysNow == 0 && nISRhard <= nVetoStep)
446  typeVetoStep = 2;
447 
448  // Update FSR dipoles.
449  if (doFSRduringProcess) timesPtr->update( iSysNow, event,
450  spacePtr->getHasWeaklyRadiated());
451  nBranch++;
452  pTLastBranch = pTspace;
453  typeLastBranch = 2;
454 
455  // Rescatter: it is possible for kinematics to fail, in which
456  // case we need to restart the parton level processing.
457  } else if (spacePtr->doRestart()) {
458  physical = false;
459  break;
460  }
461 
462  // Set maximal scales for next pT to pick.
463  pTmaxMPI = min(pTspace, pTmaxMPI);
464  pTmaxISR = pTspace;
465  pTmaxFSR = min(pTspace, pTmaxFSR);
466  pTmax = pTspace;
467  }
468 
469  // Do a final-state emission (if allowed).
470  else if (pTtimes > 0.) {
471  infoPtr->addCounter(25);
472  if (timesPtr->branch( event, true)) {
473  typeLatest = 3;
474  iSysNow = timesPtr->system();
475  ++nFSRinProc;
476  if (iSysNow == 0) ++nFSRhard;
477  if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
478  typeVetoStep = 3;
479 
480  // Update ISR dipoles.
481  if (doISR) spacePtr->update( iSysNow, event,
482  timesPtr->getHasWeaklyRadiated());
483  nBranch++;
484  pTLastBranch = pTtimes;
485  typeLastBranch = 3;
486 
487  }
488 
489  // Set maximal scales for next pT to pick.
490  pTmaxMPI = min(pTtimes, pTmaxMPI);
491  pTmaxISR = min(pTtimes, pTmaxISR);
492  pTmaxFSR = pTtimes;
493  pTmax = pTtimes;
494  }
495 
496  // If no pT scales above zero then nothing to be done.
497  else pTmax = 0.;
498 
499  // Check for double counting for Drell-Yan weak production.
500  // Only look at the second emission.
501  if ( (infoPtr->code() == 221 || infoPtr->code() == 222) &&
502  nISRhard + nFSRhard == 2 && vetoWeakJets) {
503  int id1 = event[partonSystemsPtr->getOut(0,0)].id();
504  int id2 = event[partonSystemsPtr->getOut(0,1)].id();
505  int id3 = event[partonSystemsPtr->getOut(0,2)].id();
506  Vec4 p1 = event[partonSystemsPtr->getOut(0,0)].p();
507  Vec4 p2 = event[partonSystemsPtr->getOut(0,1)].p();
508  Vec4 p3 = event[partonSystemsPtr->getOut(0,2)].p();
509 
510  // Make sure id1 is weak boson, and check that there
511  // only is a single weak boson and no photons.
512  bool doubleCountEvent = true;
513  if (abs(id1) == 24 || abs(id1) == 23) {
514  if (abs(id2) > 21 || abs(id3) > 21)
515  doubleCountEvent = false;
516  } else if (abs(id2) == 24 || abs(id2) == 23) {
517  swap(id1,id2);
518  swap(p1,p2);
519  if (abs(id3) > 21)
520  doubleCountEvent = false;
521  } else if ( abs(id3) == 24 || abs(id3) == 23) {
522  swap(id1,id3);
523  swap(p1,p3);
524  }
525 
526  if (doubleCountEvent) {
527  double d = p1.pT2();
528  bool cut = true;
529  if (p2.pT2() < d) {d = p2.pT2(); cut = false;}
530  if (p3.pT2() < d) {d = p3.pT2(); cut = false;}
531 
532  // Check for angle between weak boson and quarks.
533  // (require final state particle to be a fermion)
534  if (abs(id2) < 20) {
535  double dij = min(p1.pT2(),p2.pT2())
536  * pow2(RRapPhi(p1,p2)) / vetoWeakDeltaR2;
537  if (dij < d) {
538  d = dij;
539  cut = true;
540  }
541  }
542 
543  if (abs(id3) < 20) {
544  double dij = min(p1.pT2(),p3.pT2())
545  * pow2(RRapPhi(p1,p3)) / vetoWeakDeltaR2;
546  if (dij < d) {
547  d = dij;
548  cut = true;
549  }
550  }
551 
552  // Check for angle between recoiler and radiator,
553  // if it is a quark anti-quark pair
554  // or if the recoiler is a gluon.
555  if (abs(id2) == 21 || abs(id3) == 21 || id2 == - id3) {
556  double dij = min(p2.pT2(),p3.pT2())
557  * pow2(RRapPhi(p2,p3)) / vetoWeakDeltaR2;
558  if (dij < d) {
559  d = dij;
560  cut = false;
561  }
562  }
563 
564  // Veto event if it does not belong to Drell-Yan production.
565  if (cut) return false;
566  }
567  }
568 
569  // Optionally check for a veto after the first few interactions,
570  // or after the first few emissions, ISR or FSR, in the hardest system.
571  if (typeVetoStep == 1) {
572  doVeto = userHooksPtr->doVetoMPIStep( nMPI, event);
573  } else if (typeVetoStep > 1) {
574  doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
575  nFSRhard, event);
576  }
577 
578  // Abort event if vetoed.
579  if (doVeto) {
580  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
581  return false;
582  }
583 
584  // Keep on evolving until nothing is left to be done.
585  if (typeLatest > 0 && typeLatest < 4)
586  infoPtr->addCounter(25 + typeLatest);
587  if (!isDiff) infoPtr->setPartEvolved( nMPI, nISR);
588 
589  // Handle potential merging veto.
590  if ( canRemoveEvent && nISRhard + nFSRhard == 1 ){
591  // Simply check, and possibly reset weights.
592  mergingHooksPtr->doVetoStep( process, event );
593  }
594 
595  // End loop evolution down in pT from hard pT scale.
596  } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
597 
598  // Do all final-state emissions if not already considered above.
599  if (doFSRafterProcess && (nBranchMax <= 0 || nBranch < nBranchMax) ) {
600 
601  // Find largest scale for final partons.
602  pTmax = 0.;
603  for (int i = 0; i < event.size(); ++i)
604  if (event[i].isFinal() && event[i].scale() > pTmax)
605  pTmax = event[i].scale();
606  pTsaveFSR = pTmax;
607 
608  // Prepare all subsystems for evolution.
609  for (int iSys = 0; iSys < partonSystemsPtr->sizeSys(); ++iSys)
610  timesPtr->prepare( iSys, event);
611 
612  // Set up initial veto scale.
613  doVeto = false;
614  pTveto = pTvetoPT;
615 
616  // Begin evolution down in pT from hard pT scale.
617  do {
618  infoPtr->addCounter(29);
619  typeVetoStep = 0;
620  double pTtimes = timesPtr->pTnext( event, pTmax, 0.);
621  infoPtr->setPTnow( pTtimes);
622 
623  // Allow a user veto. Only do it once, so remember to change pTveto.
624  if (pTveto > 0. && pTveto > pTtimes) {
625  pTveto = -1.;
626  doVeto = userHooksPtr->doVetoPT( 4, event);
627  // Abort event if vetoed.
628  if (doVeto) {
629  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
630  return false;
631  }
632  }
633 
634  // Do a final-state emission (if allowed).
635  if (pTtimes > 0.) {
636  infoPtr->addCounter(30);
637  if (timesPtr->branch( event, true)) {
638  iSysNow = timesPtr->system();
639  ++nFSRinProc;
640  if (iSysNow == 0) ++nFSRhard;
641  if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
642  typeVetoStep = 4;
643 
644  nBranch++;
645  pTLastBranch = pTtimes;
646  typeLastBranch = 4;
647 
648  }
649  pTmax = pTtimes;
650  }
651 
652  // If no pT scales above zero then nothing to be done.
653  else pTmax = 0.;
654 
655  // Optionally check for a veto after the first few emissions.
656  if (typeVetoStep > 0) {
657  doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
658  nFSRhard, event);
659  // Abort event if vetoed.
660  if (doVeto) {
661  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
662  return false;
663  }
664  }
665 
666  // Handle potential merging veto.
667  if ( canRemoveEvent && nISRhard + nFSRhard == 1 ){
668  // Simply check, and possibly reset weights.
669  mergingHooksPtr->doVetoStep( process, event );
670  }
671 
672  // Keep on evolving until nothing is left to be done.
673  infoPtr->addCounter(31);
674 
675  } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
676  }
677 
678  // Handle veto after ISR + FSR + MPI, but before beam remnants
679  // and resonance decays, e.g. for MLM matching.
680  if (canVetoEarly && userHooksPtr->doVetoPartonLevelEarly( event)) {
681  doVeto = true;
682  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
683  return false;
684  }
685 
686  // Perform showers in resonance decay chains before beams & reconnection.
687  if (earlyResDec) {
688  int oldSizeEvt = event.size();
689  int oldSizeSys = partonSystemsPtr->sizeSys();
690  if (nBranchMax <= 0 || nBranch < nBranchMax)
691  doVeto = !resonanceShowers( process, event, true);
692  // Abort event if vetoed.
693  if (doVeto) return false;
694 
695  // Reassign new decay products to original system.
696  for (int iSys = oldSizeSys; iSys < partonSystemsPtr->sizeSys(); ++iSys)
697  for (int iOut = 0; iOut < partonSystemsPtr->sizeOut(iSys); ++iOut)
698  partonSystemsPtr->addOut(0, partonSystemsPtr->getOut( iSys, iOut) );
699  partonSystemsPtr->setSizeSys( oldSizeSys);
700 
701  // Perform decays and showers of W and Z emitted in shower.
702  // To do:check if W/Z emission is on in ISR or FSR??
703  if (!wzDecayShowers( event)) return false;
704 
705  // User hook to reconnect colours specifically in resonance decays.
706  if (canReconResSys && !userHooksPtr->doReconnectResonanceSystems(
707  oldSizeEvt, event)) return false;
708  }
709 
710  // Add beam remnants, including primordial kT kick and colour tracing.
711  if (!doTrial && physical && doRemnants && !remnants.add( event))
712  physical = false;
713 
714  // If no problems then done.
715  if (physical) break;
716 
717  // Else restore and loop, but do not throw existing diffractive system.
718  if (!isDiff) event.clear();
719  else {
720  event.restoreSize();
721  event.restoreJunctionSize();
722  }
723  beamAPtr->clear();
724  beamBPtr->clear();
725  partonSystemsPtr->clear();
726 
727  // End loop over ten tries. Restore from diffraction. Hopefully it worked.
728  }
729  if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
730  if (!physical) return false;
731 
732  // End big outer loop to handle two systems in double diffraction.
733  }
734 
735  // Perform showers in resonance decay chains after beams & reconnection.
736  if (!earlyResDec) {
737  int oldSizeEvt = event.size();
738  if (nBranchMax <= 0 || nBranch < nBranchMax)
739  doVeto = !resonanceShowers( process, event, true);
740  // Abort event if vetoed.
741  if (doVeto) return false;
742 
743  // Perform decays and showers of W and Z emitted in shower.
744  // To do:check if W/Z emission is on in ISR or FSR??
745  if (!wzDecayShowers( event)) return false;
746 
747  // User hook to reconnect colours specifically in resonance decays.
748  if (canReconResSys && !userHooksPtr->doReconnectResonanceSystems(
749  oldSizeEvt, event)) return false;
750  }
751 
752  // Store event properties. Not available for diffraction.
753  if (!isDiff) infoPtr->setEvolution( pTsaveMPI, pTsaveISR, pTsaveFSR,
754  nMPI, nISR, nFSRinProc, nFSRinRes);
755  if (isDiff) {
756  multiPtr->setEmpty();
757  infoPtr->setImpact( multiPtr->bMPI(), multiPtr->enhanceMPI(), false);
758  }
759 
760  // Done.
761  return true;
762 
763 }
764 
765 //--------------------------------------------------------------------------
766 
767 // Decide which diffractive subsystems are resolved (= perturbative).
768 
769 int PartonLevel::decideResolvedDiff( Event& process) {
770 
771  // Loop over two systems.
772  int nHighMass = 0;
773  int iDSmin = (isDiffC) ? 3 : 1;
774  int iDSmax = (isDiffC) ? 3 : 2;
775  for (int iDSnow = iDSmin; iDSnow <= iDSmax; ++iDSnow) {
776  int iDiffMot = iDSnow + 2;
777 
778  // Only high-mass diffractive systems should be resolved.
779  double mDiff = process[iDiffMot].m();
780  bool isHighMass = ( mDiff > mMinDiff && rndmPtr->flat()
781  < pMaxDiff * ( 1. - exp( -(mDiff - mMinDiff) / mWidthDiff ) ) );
782 
783  // Set outcome and done.
784  if (isHighMass) ++nHighMass;
785  if (iDSnow == 1) isResolvedA = isHighMass;
786  if (iDSnow == 2) isResolvedB = isHighMass;
787  if (iDSnow == 3) isResolvedC = isHighMass;
788  }
789  return nHighMass;
790 
791 }
792 
793 //--------------------------------------------------------------------------
794 
795 // Set up an unresolved process, i.e. elastic or diffractive.
796 
797 bool PartonLevel::setupUnresolvedSys( Event& process, Event& event) {
798 
799  // No hard scale in event.
800  process.scale( 0.);
801 
802  // Copy particles from process to event.
803  for (int i = 0; i < process.size(); ++ i) event.append( process[i]);
804 
805  // Loop to find diffractively excited beams.
806  for (iDS = 1; iDS < 4; ++iDS)
807  if ( (iDS == 1 && isDiffA && !isResolvedA)
808  || (iDS == 2 && isDiffB && !isResolvedB)
809  || (iDS == 3 && isDiffC && !isResolvedC) ) {
810  int iBeam = iDS + 2;
811 
812  // Diffractive mass. Boost and rotation from diffractive system
813  // rest frame, aligned along z axis, to event cm frame.
814  double mDiff = process[iBeam].m();
815  double m2Diff = mDiff * mDiff;
816  Vec4 pDiffA = (iDS == 1) ? process[1].p()
817  : process[1].p() - process[3].p();
818  Vec4 pDiffB = (iDS == 2) ? process[2].p()
819  : process[2].p() - process[4].p();
820  RotBstMatrix MtoCM;
821  MtoCM.fromCMframe( pDiffA, pDiffB);
822 
823  // Beam Particle used for flavour content kicked out by Pomeron.
824  // Randomize for central diffraction; misses closed gluon loop case.
825  bool beamSideA = (iDS == 1 || (iDS == 3 && rndmPtr->flat() < 0.5));
826  BeamParticle* beamPtr = (beamSideA) ? beamAPtr : beamBPtr;
827  if (iDS == 3) beamPtr = (beamSideA) ? beamPomAPtr : beamPomBPtr;
828 
829  // Pick quark or gluon kicked out and flavour subdivision.
830  beamPtr->newValenceContent();
831  bool gluonIsKicked = beamPtr->pickGluon(mDiff);
832  int id1 = beamPtr->pickValence();
833  int id2 = beamPtr->pickRemnant();
834 
835  // Find flavour masses. Scale them down if too big.
836  double m1 = particleDataPtr->constituentMass(id1);
837  double m2 = particleDataPtr->constituentMass(id2);
838  if (m1 + m2 > 0.5 * mDiff) {
839  double reduce = 0.5 * mDiff / (m1 + m2);
840  m1 *= reduce;
841  m2 *= reduce;
842  }
843 
844  // If quark is kicked out, then trivial kinematics in rest frame.
845  if (!gluonIsKicked) {
846  double pAbs = sqrt( pow2(m2Diff - m1*m1 - m2*m2)
847  - pow2(2. * m1 * m2) ) / (2. * mDiff);
848  if (!beamSideA) pAbs = -pAbs;
849  double e1 = (m2Diff + m1*m1 - m2*m2) / (2. * mDiff);
850  double e2 = (m2Diff + m2*m2 - m1*m1) / (2. * mDiff);
851  Vec4 p1( 0., 0., -pAbs, e1);
852  Vec4 p2( 0., 0., pAbs, e2);
853 
854  // Boost and rotate to event cm frame.
855  p1.rotbst( MtoCM);
856  p2.rotbst( MtoCM);
857 
858  // Set colours.
859  int col1, acol1, col2, acol2;
860  if (particleDataPtr->colType(id1) == 1) {
861  col1 = event.nextColTag();
862  acol1 = 0;
863  col2 = 0;
864  acol2 = col1;
865  } else {
866  col1 = 0;
867  acol1 = event.nextColTag();
868  col2 = acol1;
869  acol2 = 0;
870  }
871  // Update process colours to stay in step.
872  process.nextColTag();
873 
874  // Store partons of diffractive system and mark system decayed.
875  int iDauBeg = event.append( id1, 24, iBeam, 0, 0, 0, col1, acol1,
876  p1, m1);
877  int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
878  p2, m2);
879  event[iBeam].statusNeg();
880  event[iBeam].daughters(iDauBeg, iDauEnd);
881 
882  // If gluon is kicked out: share momentum between two remnants.
883  } else {
884  double m2Sys, zSys, pxSys, pySys, mTS1, mTS2;
885  zSys = beamPtr->zShare(mDiff, m1, m2);
886 
887  // Provide relative pT kick in remnant. Construct (transverse) masses.
888  pxSys = beamPtr->pxShare();
889  pySys = beamPtr->pyShare();
890  mTS1 = m1*m1 + pxSys*pxSys + pySys*pySys;
891  mTS2 = m2*m2 + pxSys*pxSys + pySys*pySys;
892  m2Sys = mTS1 / zSys + mTS2 / (1. - zSys);
893 
894  // Momentum of kicked-out massless gluon in diffractive rest frame.
895  double pAbs = (m2Diff - m2Sys) / (2. * mDiff);
896  double pLRem = (beamSideA) ? pAbs : -pAbs;
897  Vec4 pG( 0., 0., -pLRem, pAbs);
898  Vec4 pRem(0., 0., pLRem, mDiff - pAbs);
899 
900  // Momenta of the two beam remnant flavours. (Lightcone p+ = m_diff!)
901  double e1 = 0.5 * (zSys * mDiff + mTS1 / (zSys * mDiff));
902  double pL1 = 0.5 * (zSys * mDiff - mTS1 / (zSys * mDiff));
903  if (!beamSideA) pL1 = -pL1;
904  Vec4 p1(pxSys, pySys, pL1, e1);
905  Vec4 p2 = pRem - p1;
906 
907  // Boost and rotate to event cm frame. Improve precision.
908  pG.rotbst( MtoCM);
909  p1.rotbst( MtoCM);
910  p2.rotbst( MtoCM);
911  pG.e( pG.pAbs());
912 
913  // Set colours.
914  int colG, acolG, col1, acol1, col2, acol2;
915  if (particleDataPtr->colType(id1) == 1) {
916  col1 = event.nextColTag();
917  acol1 = 0;
918  colG = event.nextColTag();
919  acolG = col1;
920  col2 = 0;
921  acol2 = colG;
922  } else {
923  col1 = 0;
924  acol1 = event.nextColTag();
925  colG = acol1;
926  acolG = event.nextColTag();
927  col2 = acolG;
928  acol2 = 0;
929  }
930  // Update process colours to stay in step.
931  process.nextColTag();
932  process.nextColTag();
933 
934  // Store partons of diffractive system and mark system decayed.
935  int iDauBeg = event.append( 21, 24, iBeam, 0, 0, 0, colG, acolG, pG, 0.);
936  event.append( id1, 63, iBeam, 0, 0, 0, col1, acol1, p1, m1);
937  int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
938  p2, m2);
939  event[iBeam].statusNeg();
940  event[iBeam].daughters(iDauBeg, iDauEnd);
941  }
942 
943  // End loop over beams. Done.
944  }
945  return true;
946 
947 }
948 
949 //--------------------------------------------------------------------------
950 
951 // Set up the hard process(es), excluding subsequent resonance decays.
952 
953 void PartonLevel::setupHardSys( Event& process, Event& event) {
954 
955  // Incoming partons to hard process are stored in slots 3 and 4.
956  int inS = 0;
957  int inP = 3;
958  int inM = 4;
959 
960  // Mother and last entry of diffractive system. Offset.
961  int iDiffMot = iDS + 2;
962  int iDiffDau = process.size() - 1;
963  int nOffset = sizeEvent - sizeProcess;
964 
965  // Resolved diffraction means more entries.
966  if (isDiff) {
967  inS = iDiffMot;
968  inP = iDiffDau - 3;
969  inM = iDiffDau - 2;
970 
971  // Diffractively excited particle described as Pomeron-hadron beams.
972  event[inS].statusNeg();
973  event[inS].daughters( inP - 2 + nOffset, inM - 2 + nOffset);
974  }
975 
976  // If two hard interactions then find where second begins.
977  int iBeginSecond = process.size();
978  if (doSecondHard) {
979  iBeginSecond = 5;
980  while (process[iBeginSecond].status() != -21) ++iBeginSecond;
981  }
982 
983  // If incoming partons are massive then recalculate to put them massless.
984  if (process[inP].m() != 0. || process[inM].m() != 0.) {
985  double pPos = process[inP].pPos() + process[inM].pPos();
986  double pNeg = process[inP].pNeg() + process[inM].pNeg();
987  process[inP].pz( 0.5 * pPos);
988  process[inP].e( 0.5 * pPos);
989  process[inP].m( 0.);
990  process[inM].pz(-0.5 * pNeg);
991  process[inM].e( 0.5 * pNeg);
992  process[inM].m( 0.);
993  }
994 
995  // Add incoming hard-scattering partons to list in beam remnants.
996  double x1 = process[inP].pPos() / process[inS].m();
997  beamAPtr->append( inP + nOffset, process[inP].id(), x1);
998  double x2 = process[inM].pNeg() / process[inS].m();
999  beamBPtr->append( inM + nOffset, process[inM].id(), x2);
1000 
1001  // Scale. Find whether incoming partons are valence or sea. Store.
1002  // When an x-dependent matter profile is used with nonDiffractive,
1003  // trial interactions mean that the valence/sea choice has already
1004  // been made and should be restored here.
1005  double scale = process.scale();
1006  int vsc1, vsc2;
1007  beamAPtr->xfISR( 0, process[inP].id(), x1, scale*scale);
1008  if (isNonDiff && (vsc1 = multiPtr->getVSC1()) != 0)
1009  (*beamAPtr)[0].companion(vsc1);
1010  else vsc1 = beamAPtr->pickValSeaComp();
1011  beamBPtr->xfISR( 0, process[inM].id(), x2, scale*scale);
1012  if (isNonDiff && (vsc2 = multiPtr->getVSC2()) != 0)
1013  (*beamBPtr)[0].companion(vsc2);
1014  else vsc2 = beamBPtr->pickValSeaComp();
1015  bool isVal1 = (vsc1 == -3);
1016  bool isVal2 = (vsc2 == -3);
1017  infoPtr->setValence( isVal1, isVal2);
1018 
1019  // Initialize info needed for subsequent sequential decays + showers.
1020  nHardDone = sizeProcess;
1021  iPosBefShow.resize( process.size() );
1022  fill( iPosBefShow.begin(), iPosBefShow.end(), 0);
1023 
1024  // Add the beam and hard subprocess partons to the event record.
1025  for (int i = sizeProcess; i < iBeginSecond; ++i) {
1026  if (process[i].mother1() > inM) break;
1027  int j = event.append(process[i]);
1028  iPosBefShow[i] = i;
1029 
1030  // Offset history if process and event not in step.
1031  if (nOffset != 0) {
1032  int iOrd = i - iBeginSecond + 7;
1033  if (iOrd == 1 || iOrd == 2)
1034  event[j].offsetHistory( 0, 0, 0, nOffset);
1035  else if (iOrd == 3 || iOrd == 4)
1036  event[j].offsetHistory( 0, nOffset, 0, nOffset);
1037  else if (iOrd == 5 || iOrd == 6)
1038  event[j].offsetHistory( 0, nOffset, 0, 0);
1039  }
1040 
1041  // Currently outgoing ones should not count as decayed.
1042  if (event[j].status() == -22) {
1043  event[j].statusPos();
1044  event[j].daughters(0, 0);
1045  }
1046 
1047  // Complete task of copying hard subsystem into event record.
1048  ++nHardDone;
1049  }
1050 
1051  // Store participating partons as first set in list of all systems.
1052  partonSystemsPtr->addSys();
1053  partonSystemsPtr->setInA(0, inP + nOffset);
1054  partonSystemsPtr->setInB(0, inM + nOffset);
1055  for (int i = inM + 1; i < nHardDone; ++i)
1056  partonSystemsPtr->addOut(0, i + nOffset);
1057  partonSystemsPtr->setSHat( 0,
1058  (event[inP + nOffset].p() + event[inM + nOffset].p()).m2Calc() );
1059  partonSystemsPtr->setPTHat( 0, scale);
1060 
1061  // Identify second hard process where applicable.
1062  // Since internally generated incoming partons are guaranteed massless.
1063  if (doSecondHard) {
1064  int inP2 = iBeginSecond;
1065  int inM2 = iBeginSecond + 1;
1066 
1067  // Add incoming hard-scattering partons to list in beam remnants.
1068  // Not valid if not in rest frame??
1069  x1 = process[inP2].pPos() / process[0].e();
1070  beamAPtr->append( inP2, process[inP2].id(), x1);
1071  x2 = process[inM2].pNeg() / process[0].e();
1072  beamBPtr->append( inM2, process[inM2].id(), x2);
1073 
1074  // Find whether incoming partons are valence or sea.
1075  scale = process.scaleSecond();
1076  beamAPtr->xfISR( 1, process[inP2].id(), x1, scale*scale);
1077  beamAPtr->pickValSeaComp();
1078  beamBPtr->xfISR( 1, process[inM2].id(), x2, scale*scale);
1079  beamBPtr->pickValSeaComp();
1080 
1081  // Add the beam and hard subprocess partons to the event record.
1082  for (int i = inP2; i < process.size(); ++ i) {
1083  int mother = process[i].mother1();
1084  if ( (mother > 2 && mother < inP2) || mother > inM2 ) break;
1085  event.append(process[i]);
1086  iPosBefShow[i] = i;
1087 
1088  // Currently outgoing ones should not count as decayed.
1089  if (event[i].status() == -22) {
1090  event[i].statusPos();
1091  event[i].daughters(0, 0);
1092  }
1093 
1094  // Complete task of copying hard subsystem into event record.
1095  ++nHardDone;
1096  }
1097 
1098  // Store participating partons as second set in list of all systems.
1099  partonSystemsPtr->addSys();
1100  partonSystemsPtr->setInA(1, inP2);
1101  partonSystemsPtr->setInB(1, inM2);
1102  for (int i = inM2 + 1; i < nHardDone; ++i)
1103  partonSystemsPtr->addOut(1, i);
1104  partonSystemsPtr->setSHat( 1,
1105  (event[inP2].p() + event[inM2].p()).m2Calc() );
1106  partonSystemsPtr->setPTHat( 1, scale);
1107 
1108  // End code for second hard process.
1109  }
1110 
1111  // Update event colour tag to maximum in whole process.
1112  int maxColTag = 0;
1113  for (int i = 0; i < process.size(); ++ i) {
1114  if (process[i].col() > maxColTag) maxColTag = process[i].col();
1115  if (process[i].acol() > maxColTag) maxColTag = process[i].acol();
1116  }
1117  event.initColTag(maxColTag);
1118 
1119  // Copy junctions from process to event.
1120  for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1121  // Resonance decay products may not have been copied from process to
1122  // event yet. If so, do not add junctions associated with decays yet.
1123  int kindJunction = process.kindJunction(iJun);
1124  bool doCopy = true;
1125  // For junction types <= 4, check if final-state legs were copied.
1126  if (kindJunction <= 4) {
1127  int iLegF1 = (kindJunction - 1) / 2;
1128  for (int iLeg = iLegF1; iLeg <= 2; ++iLeg) {
1129  bool colFound = false;
1130  for (int i = inM + 1; i < event.size(); ++i) {
1131  int col = (kindJunction % 2 == 1) ? event[i].col() : event[i].acol();
1132  if (col == process.colJunction(iJun,iLeg)) colFound = true;
1133  }
1134  if (!colFound) doCopy = false;
1135  }
1136  }
1137  if (doCopy) event.appendJunction( process.getJunction(iJun));
1138  }
1139 
1140  // Done.
1141 }
1142 
1143 //--------------------------------------------------------------------------
1144 
1145 // Set up the event for subsequent resonance decays and showers.
1146 
1147 void PartonLevel::setupShowerSys( Event& process, Event& event) {
1148 
1149  // Reset event record to only contain line 0.
1150  event.clear();
1151  event.append( process[0]);
1152 
1153  // Initialize info needed for subsequent sequential decays + showers.
1154  nHardDone = 1;
1155  iPosBefShow.resize( process.size());
1156  fill( iPosBefShow.begin(), iPosBefShow.end(), 0);
1157 
1158  // Add the hard subprocess partons to the event record.
1159  for (int i = 1; i < process.size(); ++i) {
1160  if (process[i].mother1() > 0) break;
1161  int j = event.append(process[i]);
1162  iPosBefShow[i] = i;
1163 
1164  // Currently outgoing ones should not count as decayed.
1165  if (event[j].status() == -22) {
1166  event[j].statusPos();
1167  event[j].daughters(0, 0);
1168  }
1169 
1170  // Complete task of copying hard subsystem into event record.
1171  ++nHardDone;
1172  }
1173 
1174  // Store participating partons as first set in list of all systems.
1175  partonSystemsPtr->clear();
1176  partonSystemsPtr->addSys();
1177  for (int i = 1; i < nHardDone; ++i) partonSystemsPtr->addOut(0, i);
1178  partonSystemsPtr->setSHat( 0, pow2(process[0].m()) );
1179  partonSystemsPtr->setPTHat( 0, 0.5 * process[0].m() );
1180 
1181  // Copy junctions from process to event.
1182  for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1183  // Resonance decay products may not have been copied from process to
1184  // event yet. If so, do not add junctions associated with decays yet.
1185  int kindJunction = process.kindJunction(iJun);
1186  bool doCopy = true;
1187  // For junction types <= 4, check if final-state legs were copied.
1188  if (kindJunction <= 4) {
1189  int iLegF1 = (kindJunction - 1) / 2;
1190  for (int iLeg = iLegF1; iLeg <= 2; ++iLeg) {
1191  bool colFound = false;
1192  for (int i = 1; i < event.size(); ++i) {
1193  int col = (kindJunction % 2 == 1) ? event[i].col() : event[i].acol();
1194  if (col == process.colJunction(iJun,iLeg)) colFound = true;
1195  }
1196  if (!colFound) doCopy = false;
1197  }
1198  }
1199  if (doCopy) event.appendJunction( process.getJunction(iJun));
1200  }
1201 
1202  // Done.
1203 }
1204 
1205 //--------------------------------------------------------------------------
1206 
1207 // Resolved diffraction: replace full event with diffractive subsystem.
1208 
1209 void PartonLevel::setupResolvedDiff( Event& process) {
1210 
1211  // Mother and last entry of diffractive system.
1212  int iDiffMot = iDS + 2;
1213  int iDiffDau = process.size() - 1;
1214 
1215  // Diffractively excited particle to be replaced by Pomeron-hadron beams
1216  // (or Pomeron-Pomeron beams for central diffraction).
1217  process[iDiffMot].statusNeg();
1218  process[iDiffMot].daughters( iDiffDau + 1, iDiffDau + 2);
1219 
1220  // Diffractive system mass.
1221  double mDiff = process[iDiffMot].m();
1222  double m2Diff = mDiff * mDiff;
1223 
1224  // Set up Pomeron-proton or Pomeron-Pomeron system as if it were
1225  // the complete collision. Set Pomeron "beam particle" massless.
1226  int idDiffA = (iDS == 1) ? process[1].id() : 990;
1227  int idDiffB = (iDS == 2) ? process[2].id() : 990;
1228  double mDiffA = (iDS == 1) ? process[1].m() : 0.;
1229  double mDiffB = (iDS == 2) ? process[2].m() : 0.;
1230  double m2DiffA = mDiffA * mDiffA;
1231  double m2DiffB = mDiffB * mDiffB;
1232  double eDiffA = 0.5 * (m2Diff + m2DiffA - m2DiffB) / mDiff;
1233  double eDiffB = 0.5 * (m2Diff + m2DiffB - m2DiffA) / mDiff;
1234  double pzDiff = 0.5 * sqrtpos( pow2(m2Diff - m2DiffA - m2DiffB)
1235  - 4. * m2DiffA * m2DiffB ) / mDiff;
1236  process.append( idDiffA, 13, iDiffMot, 0, 0, 0, 0, 0,
1237  0., 0., pzDiff, eDiffA, mDiffA);
1238  process.append( idDiffB, 13, iDiffMot, 0, 0, 0, 0, 0,
1239  0., 0., -pzDiff, eDiffB, mDiffB);
1240 
1241  // Reassign beam pointers to refer to subsystem effective beams.
1242  beamAPtr = (iDS == 1) ? beamHadAPtr : beamPomAPtr;
1243  beamBPtr = (iDS == 2) ? beamHadBPtr : beamPomBPtr;
1244 
1245  // Pretend that the diffractive system is the whole collision.
1246  eCMsave = infoPtr->eCM();
1247  infoPtr->setECM( mDiff);
1248  beamAPtr->newPzE( pzDiff, eDiffA);
1249  beamBPtr->newPzE( -pzDiff, eDiffB);
1250 
1251  // Beams not found in normal slots 1 and 2.
1252  int beamOffset = (sizeEvent > 0) ? sizeEvent - 1 : 4;
1253 
1254  // Reassign beam pointers in other classes.
1255  timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1256  spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1257  remnants.reassignBeamPtrs( beamAPtr, beamBPtr, iDS);
1258 
1259  // Reassign multiparton interactions pointer to right object.
1260  if (iDS == 1) multiPtr = &multiSDA;
1261  else if (iDS == 2) multiPtr = &multiSDB;
1262  else multiPtr = &multiCD;
1263 
1264 }
1265 
1266 //--------------------------------------------------------------------------
1267 
1268 // Resolved diffraction: restore to original behaviour.
1269 
1270 void PartonLevel::leaveResolvedDiff( int iHardLoop, Event& process,
1271  Event& event) {
1272 
1273  // Reconstruct boost and rotation to event cm frame.
1274  Vec4 pDiffA = (iDS == 1) ? process[1].p()
1275  : process[1].p() - process[3].p();
1276  Vec4 pDiffB = (iDS == 2) ? process[2].p()
1277  : process[2].p() - process[4].p();
1278  RotBstMatrix MtoCM;
1279  MtoCM.fromCMframe( pDiffA, pDiffB);
1280 
1281  // Perform rotation and boost on diffractive system.
1282  for (int i = sizeProcess; i < process.size(); ++i)
1283  process[i].rotbst( MtoCM);
1284  int iFirst = (iHardLoop == 1) ? 5 + sizeEvent - sizeProcess : sizeEvent;
1285  if(isDiffC) iFirst = 6 + sizeEvent - sizeProcess;
1286  for (int i = iFirst; i < event.size(); ++i)
1287  event[i].rotbst( MtoCM);
1288 
1289  // Restore cm energy.
1290  infoPtr->setECM( eCMsave);
1291  beamAPtr->newPzE( event[1].pz(), event[1].e());
1292  beamBPtr->newPzE( event[2].pz(), event[2].e());
1293 
1294  // Restore beam pointers to incoming hadrons.
1295  beamAPtr = beamHadAPtr;
1296  beamBPtr = beamHadBPtr;
1297 
1298  // Reassign beam pointers in other classes.
1299  timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1300  spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1301  remnants.reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1302 
1303  // Restore multiparton interactions pointer to default object.
1304  multiPtr = &multiMB;
1305 
1306 }
1307 
1308 //--------------------------------------------------------------------------
1309 
1310 // Handle showers in successive resonance decays.
1311 
1312 bool PartonLevel::resonanceShowers( Event& process, Event& event,
1313  bool skipForR) {
1314 
1315  // Prepare to start over from beginning for R-hadron decays.
1316  if (allowRH) {
1317  if (skipForR) {
1318  nHardDoneRHad = nHardDone;
1319  inRHadDecay.resize(0);
1320  for (int i = 0; i < process.size(); ++i)
1321  inRHadDecay.push_back( false);
1322  } else nHardDone = nHardDoneRHad;
1323  }
1324 
1325  // Isolate next system to be processed, if anything remains.
1326  int nRes = 0;
1327  int nFSRres = 0;
1328  // Number of desired branchings, negative value means no restriction.
1329  int nBranchMax = (doTrial) ? nTrialEmissions : -1;
1330  // Vector to tell which junctions have already been copied
1331  vector<int> iJunCopied;
1332 
1333  while (nHardDone < process.size()) {
1334  ++nRes;
1335  int iBegin = nHardDone;
1336 
1337  // In first call (skipForR = true) skip over resonances
1338  // that should form R-hadrons, and their daughters.
1339  if (allowRH) {
1340  if (skipForR) {
1341  bool comesFromR = false;
1342  int iTraceUp = iBegin;
1343  do {
1344  if ( rHadronsPtr->givesRHadron(process[iTraceUp].id()) )
1345  comesFromR = true;
1346  iTraceUp = process[iTraceUp].mother1();
1347  } while (iTraceUp > 0 && !comesFromR);
1348  if (comesFromR) {
1349  inRHadDecay[iBegin] = true;
1350  ++nHardDone;
1351  continue;
1352  }
1353 
1354  // In optional second call (skipForR = false) process decay chains
1355  // inside R-hadrons.
1356  } else if (!inRHadDecay[iBegin]) {
1357  ++nHardDone;
1358  continue;
1359  }
1360  }
1361 
1362  // Mother in hard process and in complete event (after shower).
1363  int iHardMother = process[iBegin].mother1();
1364  Particle& hardMother = process[iHardMother];
1365  int iBefMother = iPosBefShow[iHardMother];
1366  int iAftMother = event.iBotCopyId(iBefMother);
1367  // Possibly trace across intermediate R-hadron state.
1368  if (allowRH) {
1369  int iTraceRHadron = rHadronsPtr->trace( iAftMother);
1370  if (iTraceRHadron > 0) iAftMother = iTraceRHadron;
1371  }
1372  Particle& aftMother = event[iAftMother];
1373 
1374  // From now on mother counts as decayed.
1375  aftMother.statusNeg();
1376 
1377  // Mother can have been moved by showering (in any of previous steps),
1378  // so prepare to update colour and momentum information for system.
1379  int colBef = hardMother.col();
1380  int acolBef = hardMother.acol();
1381  int colAft = aftMother.col();
1382  int acolAft = aftMother.acol();
1383  RotBstMatrix M;
1384  M.bst( hardMother.p(), aftMother.p());
1385 
1386  // Extract next partons from hard event into normal event record.
1387  vector<bool> doCopyJun;
1388  for (int i = iBegin; i < process.size(); ++i) {
1389  if (process[i].mother1() != iHardMother) break;
1390  int iNow = event.append( process[i] );
1391  iPosBefShow[i] = iNow;
1392  Particle& now = event.back();
1393 
1394  // Currently outgoing ones should not count as decayed.
1395  if (now.status() == -22) {
1396  now.statusPos();
1397  now.daughters(0, 0);
1398  }
1399 
1400  // Update daughter and mother information.
1401  if (i == iBegin) event[iAftMother].daughter1( iNow);
1402  else event[iAftMother].daughter2( iNow);
1403  now.mother1(iAftMother);
1404 
1405  // Check if this parton carries a junction color in hard event.
1406  for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1407  if (iJun >= int(doCopyJun.size())) doCopyJun.push_back(false);
1408  // Only consider junctions that can appear in decays.
1409  int kindJunction = process.kindJunction(iJun);
1410  if (kindJunction >= 5) continue;
1411  int col = (kindJunction % 2 == 1) ? now.col() : now.acol();
1412  int iLegF1 = (kindJunction - 1) / 2;
1413  for (int iLeg = iLegF1; iLeg <= 2; ++iLeg)
1414  if (col == process.colJunction(iJun,iLeg)) doCopyJun[iJun] = true;
1415  }
1416 
1417  // Update colour and momentum information.
1418  if (now.col() == colBef) now.col( colAft);
1419  if (now.acol() == acolBef) now.acol( acolAft);
1420  now.rotbst( M);
1421 
1422  // Update vertex information.
1423  if (now.hasVertex()) now.vProd( event[iAftMother].vDec() );
1424 
1425  // Complete task of copying next subsystem into event record.
1426  ++nHardDone;
1427  }
1428  int iEnd = nHardDone - 1;
1429 
1430  // Copy down junctions from hard event into normal event record.
1431  for (int iJun = 0; iJun < int(doCopyJun.size()); ++iJun) {
1432  // Check if this junction was already copied
1433  for (int jJun = 0; jJun < int(iJunCopied.size()); ++jJun)
1434  if (iJunCopied[jJun] == iJun) doCopyJun[iJun] = false;
1435  // Skip if not doing anything
1436  if (!doCopyJun[iJun]) continue;
1437  // Check for changed colors and update as necessary.
1438  Junction junCopy = process.getJunction(iJun);
1439  for (int iLeg = 0; iLeg <= 2; ++iLeg) {
1440  int colLeg = junCopy.col(iLeg);
1441  if (colLeg == colBef) junCopy.col(iLeg, colAft);
1442  if (colLeg == acolBef) junCopy.col(iLeg, acolAft);
1443  }
1444  event.appendJunction(junCopy);
1445  // Mark junction as copied (to avoid later recopying)
1446  iJunCopied.push_back(iJun);
1447  }
1448 
1449  // Reset pT of last branching
1450  pTLastBranch = 0.0;
1451 
1452  // Add new system, automatically with two empty beam slots.
1453  int iSys = partonSystemsPtr->addSys();
1454  partonSystemsPtr->setSHat(iSys, pow2(hardMother.m()) );
1455  partonSystemsPtr->setPTHat(iSys, 0.5 * hardMother.m() );
1456 
1457  // Loop over allowed range to find all final-state particles.
1458  for (int i = iPosBefShow[iBegin]; i <= iPosBefShow[iEnd]; ++i)
1459  if (event[i].isFinal()) partonSystemsPtr->addOut( iSys, i);
1460 
1461  // Do parton showers inside subsystem: maximum scale by mother mass.
1462  if (doFSRinResonances) {
1463  double pTmax = 0.5 * hardMother.m();
1464  if (canSetScale) pTmax
1465  = userHooksPtr->scaleResonance( iAftMother, event);
1466  nFSRhard = 0;
1467 
1468  // Set correct scale for trial showers.
1469  if (doTrial) pTmax = process.scale();
1470 
1471  // Let prepare routine do the setup.
1472  timesDecPtr->prepare( iSys, event);
1473 
1474  // Number of actual branchings
1475  int nBranch = 0;
1476 
1477  // Set up initial veto scale.
1478  doVeto = false;
1479  double pTveto = pTvetoPT;
1480 
1481  // Begin evolution down in pT from hard pT scale.
1482  do {
1483  typeVetoStep = 0;
1484  double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
1485 
1486  // Allow a user veto. Only do it once, so remember to change pTveto.
1487  if (pTveto > 0. && pTveto > pTtimes) {
1488  pTveto = -1.;
1489  doVeto = userHooksPtr->doVetoPT( 5, event);
1490  // Abort event if vetoed.
1491  if (doVeto) return false;
1492  }
1493 
1494  // Do a final-state emission (if allowed).
1495  if (pTtimes > 0.) {
1496  if (timesDecPtr->branch( event)) {
1497  ++nFSRres;
1498  ++nFSRhard;
1499  if (canVetoStep && nFSRhard <= nVetoStep) typeVetoStep = 5;
1500 
1501  nBranch++;
1502  pTLastBranch = pTtimes;
1503  typeLastBranch = 5;
1504 
1505  }
1506  pTmax = pTtimes;
1507  }
1508 
1509  // If no pT scales above zero then nothing to be done.
1510  else pTmax = 0.;
1511 
1512  // Optionally check for a veto after the first few emissions.
1513  if (typeVetoStep > 0) {
1514  doVeto = userHooksPtr->doVetoStep( typeVetoStep, 0, nFSRhard,
1515  event);
1516  // Abort event if vetoed.
1517  if (doVeto) return false;
1518  }
1519 
1520  // Handle potential merging veto.
1521  if ( canRemoveEvent && nFSRhard == 1 ) {
1522  // Simply check, and possibly reset weights.
1523  mergingHooksPtr->doVetoStep( process, event, true );
1524  }
1525 
1526  // Keep on evolving until nothing is left to be done.
1527  } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
1528 
1529  }
1530 
1531  // No more systems to be processed. Set total number of emissions.
1532  }
1533  if (skipForR) nFSRinRes = nFSRres;
1534  return true;
1535 
1536 }
1537 
1538 //--------------------------------------------------------------------------
1539 
1540 // Perform decays and showers of W and Z emitted in shower.
1541 
1542 bool PartonLevel::wzDecayShowers( Event& event) {
1543 
1544  // Identify W/Z produced by a parton shower.
1545  for (int iWZ = 0; iWZ < event.size(); ++iWZ)
1546  if (event[iWZ].isFinal()
1547  && (event[iWZ].id() == 23 || event[iWZ].idAbs() == 24) ) {
1548  int iWZtop = event[iWZ].iTopCopy();
1549  int typeWZ = 0;
1550  if (event[iWZtop].statusAbs() == 56) typeWZ = 1;
1551  if (event[iWZtop].statusAbs() == 47) typeWZ = 2;
1552  int sizeSave = event.size();
1553 
1554  // Map id_Z = 23 -> 93 and id_W = 24 -> 94, for separate decay settings.
1555  // Let W/Z resonance decay. Restore correct identity and status codes.
1556  if (typeWZ > 0) {
1557  int idSave = event[iWZ].id();
1558  event[iWZ].id( (idSave > 0) ? idSave + 70 : idSave - 70);
1559  int statusSave = event[iWZ].status();
1560  resonanceDecays.next( event, iWZ);
1561  event[iWZ].id( idSave);
1562  if (event.size() - sizeSave != 2) return true;
1563  event[iWZ].status( -statusSave);
1564  }
1565 
1566  // FSR decays.
1567  if (typeWZ == 1) {
1568 
1569  // Identify fermion after W/Z emission.
1570  vector<int> iSisters = event[iWZtop].sisterList();
1571  if (iSisters.size() != 1) {
1572  infoPtr->errorMsg("Error in PartonLevel::wzDecayShowers: "
1573  "Not able to find a single sister particle");
1574  return false;
1575  }
1576  int iEmitter = iSisters[0];
1577 
1578  // Boosts to study decay in W/Z rest frame.
1579  RotBstMatrix MtoNew, MtoRest, MtoCM;
1580  MtoNew.bst( event[iWZtop].p(), event[iWZ].p());
1581  MtoRest.bstback( event[iWZ].p());
1582  MtoCM.bst( event[iWZ].p());
1583 
1584  // Emitter and recoiler in W/Z rest frame.
1585  Vec4 pEmitter = event[iEmitter].p();
1586  pEmitter.rotbst( MtoNew);
1587  pEmitter.rotbst( MtoRest);
1588  if (event[iWZtop + 1].statusAbs() != 52) {
1589  infoPtr->errorMsg("Error in PartonLevel::wzDecayShowers: "
1590  "Found wrong recoiler");
1591  return false;
1592  }
1593  Vec4 pRecoiler = event[iWZtop + 1].p();
1594  pRecoiler.rotbst( MtoNew);
1595  pRecoiler.rotbst( MtoRest);
1596  Vec4 pWZRest = event[iWZ].p();
1597  pWZRest.rotbst( MtoRest);
1598 
1599  // Always choose p4 as the particle and p5 as the anti-particle.
1600  Vec4 p4 = pEmitter;
1601  Vec4 p5 = pRecoiler;
1602  if (event[iEmitter].id() < 0) swap( p4, p5);
1603 
1604  // Decay daughters in W/Z rest frame.
1605  // Always choose pDec1 as the particle and p2Dec as the anti-particle.
1606  Vec4 pDec1 = event[sizeSave].p();
1607  Vec4 pDec2 = event[sizeSave + 1].p();
1608  if (event[sizeSave].id() < 0) swap( pDec1, pDec2);
1609  pDec1.rotbst( MtoRest);
1610  pDec2.rotbst( MtoRest);
1611 
1612  // Couplings.
1613  double li2, ri2, lf2, rf2;
1614  // Z couplings: make use of initial fermion polarization if set.
1615  if (event[iWZ].id() == 23) {
1616  li2 = pow2(couplingsPtr->lf( event[iEmitter].idAbs() ));
1617  ri2 = pow2(couplingsPtr->rf( event[iEmitter].idAbs() ));
1618  lf2 = pow2(couplingsPtr->lf( event[sizeSave].idAbs() ));
1619  rf2 = pow2(couplingsPtr->rf( event[sizeSave].idAbs() ));
1620  if ( abs( event[iEmitter].pol() + 1.) < 0.1) ri2 = 0.;
1621  if ( abs( event[iEmitter].pol() - 1.) < 0.1) li2 = 0.;
1622  // W couplings.
1623  } else {
1624  li2 = 1.;
1625  ri2 = 0.;
1626  lf2 = 1.;
1627  rf2 = 0.;
1628  }
1629 
1630  // Different needed kinematic variables.
1631  double sWZER = (p4 + pWZRest + p5).m2Calc();
1632  double x1 = 2. * p4 * (p4 + pWZRest + p5) / sWZER;
1633  double x2 = 2. * p5 * (p4 + pWZRest + p5) / sWZER;
1634  double x1s = x1 * x1;
1635  double x2s = x2 * x2;
1636  double m2Sel = pWZRest.m2Calc();
1637  double rWZER = m2Sel / sWZER;
1638 
1639  // Calculate constants needed in correction.
1640  double con[9];
1641  con[0] = 2. * m2Sel * (1.-x1) * ((x2s+1.-x1-x2) - rWZER * (1.-x2))
1642  * (li2 * lf2 + ri2 * rf2);
1643  con[1] = 2. * m2Sel * (1.-x2) * ((x1s+1.-x1-x2) - rWZER * (1.-x1))
1644  * (li2 * rf2 + ri2 * lf2);
1645  con[2] = 2. * m2Sel * (1.-x1) * ((x2s+1.-x1-x2) - rWZER * (1.-x2))
1646  * (li2 * rf2 + ri2 * lf2);
1647  con[3] = 2. * m2Sel * (1.-x2) * ((x1s+1.-x1-x2) - rWZER * (1.-x1))
1648  * (li2 * lf2 + ri2 * rf2);
1649  con[4] = m2Sel * sWZER * (1.-x1) * (1.-x2) * ((x1+x2-1.) + rWZER)
1650  * (li2 + ri2) * (lf2 + rf2);
1651  con[5] = -4. * (1.-x1) * (1.-x2) * (li2 + ri2) * (lf2 + rf2);
1652  con[6] = -4. * (1.-x1) * (1.-x2) * (li2 + ri2) * (lf2 + rf2);
1653  con[7] = 4. * (1.-x1) * ((1.-x1) - rWZER * (1.-x2))
1654  * (li2 + ri2) * (lf2 + rf2);
1655  con[8] = 4. * (1.-x2) * ((1.-x2) - rWZER * (1.-x1))
1656  * (li2 + ri2) * (lf2 + rf2);
1657 
1658  // Find maximum value: try pDec1 and pDec2 = -pDec1 along +-x, +-y, +-z.
1659  double wtMax = 0.;
1660  double pAbs12 = pDec1.pAbs();
1661  for (int j = 0; j < 6; ++j) {
1662  Vec4 pDec1Test( 0., 0., 0., pDec1.e());
1663  Vec4 pDec2Test( 0., 0., 0., pDec2.e());
1664  if (j == 0) { pDec1Test.px( pAbs12); pDec1Test.px( -pAbs12);}
1665  else if (j == 1) { pDec1Test.px( -pAbs12); pDec1Test.px( pAbs12);}
1666  else if (j == 2) { pDec1Test.py( pAbs12); pDec1Test.py( -pAbs12);}
1667  else if (j == 3) { pDec1Test.py( -pAbs12); pDec1Test.py( pAbs12);}
1668  else if (j == 4) { pDec1Test.pz( pAbs12); pDec1Test.pz( -pAbs12);}
1669  else if (j == 5) { pDec1Test.pz( -pAbs12); pDec1Test.pz( pAbs12);}
1670 
1671  // Evaluate matrix element and compare with current maximum.
1672  double p2p4Test = p4 * pDec1Test;
1673  double p3p4Test = p4 * pDec2Test;
1674  double p2p5Test = p5 * pDec1Test;
1675  double p3p5Test = p5 * pDec2Test;
1676  double testValues[9] = { p2p4Test, p2p5Test, p3p4Test, p3p5Test, 1.,
1677  p2p5Test * p3p4Test, p2p4Test * p3p5Test, p2p4Test * p3p4Test,
1678  p2p5Test * p3p5Test};
1679  double wtTest = 0.;
1680  for (int i = 0; i < 9; ++i) wtTest += con[i] * testValues[i];
1681  if (wtTest > wtMax) wtMax = wtTest;
1682  }
1683 
1684  // Multiply by four to ensure maximum is an overestimate.
1685  wtMax *= 4.;
1686 
1687  // Iterate with new angles until weighting succeeds.
1688  int nRot = -1;
1689  double wt = 0.;
1690  do {
1691  ++nRot;
1692  if (nRot > 0) {
1693  RotBstMatrix MrndmRot;
1694  MrndmRot.rot( acos(2. * rndmPtr->flat() - 1.),
1695  2. * M_PI * rndmPtr->flat());
1696  pDec1.rotbst(MrndmRot);
1697  pDec2.rotbst(MrndmRot);
1698  }
1699 
1700  // p2 is decay product, p3 is anti decay product,
1701  // p4 is dipole particle, p5 is dipole anti particle.
1702  // So far assumed that we always have qQ-dipole.
1703  double p2p4 = p4 * pDec1;
1704  double p3p4 = p4 * pDec2;
1705  double p2p5 = p5 * pDec1;
1706  double p3p5 = p5 * pDec2;
1707 
1708  // Calculate weight and compare with maximum weight.
1709  double wtValues[9] = { p2p4, p2p5, p3p4, p3p5, 1., p2p5 * p3p4,
1710  p2p4 * p3p5, p2p4 * p3p4, p2p5 * p3p5};
1711  wt = 0.;
1712  for (int i = 0; i < 9; ++i) wt += con[i] * wtValues[i];
1713  if (wt > wtMax || wt < 0.) {
1714  infoPtr->errorMsg("Error in PartonLevel::wzDecayShowers: "
1715  "wt bigger than wtMax or less than zero");
1716  return false;
1717  }
1718  } while (wt < wtMax * rndmPtr->flat());
1719 
1720  // If momenta rotated then store new ones.
1721  if (nRot > 0) {
1722  pDec1.rotbst( MtoCM);
1723  pDec2.rotbst( MtoCM);
1724  if (event[sizeSave].id() > 0) {
1725  event[sizeSave].p( pDec1);
1726  event[sizeSave + 1].p( pDec2);
1727  }
1728  else {
1729  event[sizeSave].p( pDec2);
1730  event[sizeSave + 1].p( pDec1);
1731  }
1732  }
1733  }
1734 
1735  // ISR decays.
1736  if (typeWZ == 2) {
1737 
1738  // Identify mother of W/Z emission.
1739  double iMother = event[iWZtop].mother1();
1740 
1741  // Boosts to study decay in W/Z rest frame.
1742  RotBstMatrix MtoNew, MtoRest, MtoCM;
1743  MtoNew.bst( event[iWZtop].p(), event[iWZ].p());
1744  MtoRest.bstback( event[iWZ].p());
1745  MtoCM.bst( event[iWZ].p());
1746 
1747  // Find recoiler.
1748  int iRecoiler;
1749  if (event[iMother + 1].statusAbs() == 42) iRecoiler = iMother + 1;
1750  else if (event[iMother - 1].statusAbs() == 42) iRecoiler = iMother - 1;
1751  else {
1752  infoPtr->errorMsg("Error in PartonLevel::wzDecayShowers: "
1753  "Found wrong recoiler");
1754  return false;
1755  }
1756 
1757  // Emitter and recoiler in W/Z rest frame.
1758  Vec4 pMother = event[iMother].p();
1759  pMother.rotbst( MtoNew);
1760  pMother.rotbst( MtoRest);
1761  Vec4 pRecoiler = event[iRecoiler].p();
1762  pRecoiler.rotbst( MtoNew);
1763  pRecoiler.rotbst( MtoRest);
1764  Vec4 pWZRest = event[iWZ].p();
1765  pWZRest.rotbst( MtoRest);
1766 
1767  // Always choose p1 as the particle and p2 as the anti-particle.
1768  // If no anti-particles just continue.
1769  Vec4 p1 = pMother;
1770  Vec4 p2 = pRecoiler;
1771  if (event[iMother].id() < 0) swap( p1, p2);
1772 
1773  // Decay daughters in W/Z rest frame.
1774  // Always choose pDec1 as the particle and p2Dec as the anti-particle.
1775  Vec4 pDec1 = event[sizeSave].p();
1776  Vec4 pDec2 = event[sizeSave + 1].p();
1777  if (event[sizeSave].id() < 0) swap( pDec1, pDec2);
1778  pDec1.rotbst( MtoRest);
1779  pDec2.rotbst( MtoRest);
1780 
1781  // Couplings.
1782  double li2, ri2, lf2, rf2;
1783  // Z couplings: make use of initial fermion polarization if set.
1784  if (event[iWZ].id() == 23) {
1785  li2 = pow2(couplingsPtr->lf( event[iMother].idAbs() ));
1786  ri2 = pow2(couplingsPtr->rf( event[iMother].idAbs() ));
1787  lf2 = pow2(couplingsPtr->lf( event[sizeSave].idAbs() ));
1788  rf2 = pow2(couplingsPtr->rf( event[sizeSave].idAbs() ));
1789  if ( abs( event[iMother].pol() + 1.) < 0.1) ri2 = 0.;
1790  if ( abs( event[iMother].pol() - 1.) < 0.1) li2 = 0.;
1791  // W couplings.
1792  } else {
1793  li2 = 1.;
1794  ri2 = 0.;
1795  lf2 = 1.;
1796  rf2 = 0.;
1797  }
1798 
1799  // Different needed kinematic variables.
1800  double s = (p1 + p2).m2Calc();
1801  double t = (p1-pWZRest).m2Calc();
1802  double u = (p2-pWZRest).m2Calc();
1803  double m2Sel = pWZRest.m2Calc();
1804  double m2Final = t + u + s - m2Sel;
1805 
1806  // Calculate constants needed in correction.
1807  double con[9];
1808  con[0] = 2. * m2Sel * (u * (1. - m2Final / t) + s)
1809  * (li2 * rf2 + ri2 * lf2);
1810  con[1] = 2. * m2Sel * (u * (1. - m2Final / t) + s)
1811  * (li2 * lf2 + ri2 * rf2);
1812  con[2] = 2. * m2Sel * (t * (1. - m2Final / u) + s)
1813  * (li2 * lf2 + ri2 * rf2);
1814  con[3] = 2. * m2Sel * (t * (1. - m2Final / u) + s)
1815  * (li2 * rf2 + ri2 * lf2);
1816  con[4] = m2Sel * m2Final * s * (li2 + ri2) * (lf2 + rf2);
1817  con[5] = -4. * m2Final * (li2 + ri2) * (lf2 + rf2);
1818  con[6] = -4. * m2Final * (li2 + ri2) * (lf2 + rf2);
1819  con[7] = 4. * (m2Final * u / t - m2Sel) * (li2 + ri2) * (lf2 + rf2);
1820  con[8] = 4. * (m2Final * t / u - m2Sel) * (li2 + ri2) * (lf2 + rf2);
1821 
1822  // Find maximum value: try pDec1 and pDec2 = -pDec1 along +-x, +-y, +-z.
1823  double wtMax = 0.;
1824  double pAbs12 = pDec1.pAbs();
1825  for (int j = 0; j < 6; ++j) {
1826  Vec4 pDec1Test( 0., 0., 0., pDec1.e());
1827  Vec4 pDec2Test( 0., 0., 0., pDec2.e());
1828  if (j == 0) { pDec1Test.px( pAbs12); pDec1Test.px( -pAbs12);}
1829  else if (j == 1) { pDec1Test.px( -pAbs12); pDec1Test.px( pAbs12);}
1830  else if (j == 2) { pDec1Test.py( pAbs12); pDec1Test.py( -pAbs12);}
1831  else if (j == 3) { pDec1Test.py( -pAbs12); pDec1Test.py( pAbs12);}
1832  else if (j == 4) { pDec1Test.pz( pAbs12); pDec1Test.pz( -pAbs12);}
1833  else if (j == 5) { pDec1Test.pz( -pAbs12); pDec1Test.pz( pAbs12);}
1834 
1835  // Evaluate matrix element and compare with current maximum.
1836  double p1p4Test = p1 * pDec1Test;
1837  double p1p5Test = p1 * pDec2Test;
1838  double p2p4Test = p2 * pDec1Test;
1839  double p2p5Test = p2 * pDec2Test;
1840  double testValues[9] = { p1p4Test, p1p5Test, p2p4Test, p2p5Test, 1.,
1841  p1p4Test * p2p5Test, p1p5Test * p2p4Test, p1p4Test * p1p5Test,
1842  p2p4Test * p2p5Test};
1843  double wtTest = 0.;
1844  for (int i = 0; i < 9; ++i) wtTest += con[i] * testValues[i];
1845  if (wtTest > wtMax) wtMax = wtTest;
1846  }
1847 
1848  // Multiply by four to ensure maximum is an overestimate.
1849  wtMax *= 4.;
1850 
1851  // Iterate with new angles until weighting succeeds.
1852  int nRot = -1;
1853  double wt = 0.;
1854  do {
1855  ++nRot;
1856  if (nRot > 0) {
1857  RotBstMatrix MrndmRot;
1858  MrndmRot.rot( acos(2. * rndmPtr->flat() - 1.),
1859  2. * M_PI * rndmPtr->flat());
1860  pDec1.rotbst(MrndmRot);
1861  pDec2.rotbst(MrndmRot);
1862  }
1863 
1864  // p2 is decay product, p3 is anti decay product,
1865  // p4 is dipole particle, p5 is dipole anti particle.
1866  // So far assumed that we always have qQ-dipole.
1867  double p1p4 = p1 * pDec1;
1868  double p1p5 = p1 * pDec2;
1869  double p2p4 = p2 * pDec1;
1870  double p2p5 = p2 * pDec2;
1871 
1872  // Calculate weight and compare with maximum weight.
1873  double wtValues[9] = { p1p4, p1p5, p2p4, p2p5, 1., p1p4 * p2p5,
1874  p1p5 * p2p4, p1p4 * p1p5, p2p4 * p2p5};
1875  wt = 0.;
1876  for (int i = 0; i < 9; ++i) wt += con[i] * wtValues[i];
1877  if (wt > wtMax || wt < 0.) {
1878  infoPtr->errorMsg("Error in PartonLevel::wzDecayShowers: "
1879  "wt bigger than wtMax or less than zero");
1880  return false;
1881  }
1882  } while (wt < wtMax * rndmPtr->flat());
1883 
1884  // If momenta rotated then store new ones.
1885  if (nRot > 0) {
1886  pDec1.rotbst( MtoCM);
1887  pDec2.rotbst( MtoCM);
1888  if (event[sizeSave].id() > 0) {
1889  event[sizeSave].p( pDec1);
1890  event[sizeSave + 1].p( pDec2);
1891  }
1892  else {
1893  event[sizeSave].p( pDec2);
1894  event[sizeSave + 1].p( pDec1);
1895  }
1896  }
1897  }
1898 
1899  // Add new system, automatically with two empty beam slots.
1900  if (typeWZ > 0) {
1901  // Maximum shower scale set by mother mass.
1902  double pTmax = 0.5 * event[iWZ].m();
1903  int iSys = partonSystemsPtr->addSys();
1904  partonSystemsPtr->setSHat(iSys, pow2(event[iWZ].m()) );
1905  partonSystemsPtr->setPTHat(iSys, pTmax );
1906  for (int i = sizeSave; i < event.size(); ++i)
1907  partonSystemsPtr->addOut( iSys, i);
1908 
1909  // Do parton showers inside subsystem.
1910  if (doFSRinResonances) {
1911 
1912  // Let prepare routine do the setup.
1913  timesDecPtr->prepare( iSys, event);
1914 
1915  // Begin evolution down in pT from hard pT scale.
1916  do {
1917  double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
1918 
1919  // Do a final-state emission (if allowed).
1920  if (pTtimes > 0.) {
1921  timesDecPtr->branch( event);
1922  pTmax = pTtimes;
1923  }
1924 
1925  // If no pT scales above zero then nothing to be done.
1926  else pTmax = 0.;
1927 
1928  // Keep on evolving until nothing is left to be done.
1929  } while (pTmax > 0.);
1930  }
1931  }
1932 
1933  // End loop over event to find W/Z gauge bosons.
1934  }
1935 
1936  // Done.
1937  return true;
1938 
1939 }
1940 
1941 //==========================================================================
1942 
1943 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5