]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8170/src/PartonLevel.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Stupid bug fix in new superlight mode (from Zurich airport)
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / PartonLevel.cxx
CommitLineData
63ba5337 1// PartonLevel.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2012 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 "PartonLevel.h"
9
10namespace 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.
22const int PartonLevel::NTRY = 10;
23
24//--------------------------------------------------------------------------
25
26// Main routine to initialize the parton-level generation process.
27
28bool 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 bool doMB = settings.flag("SoftQCD:minBias");
59 bool doSD = settings.flag("SoftQCD:singleDiffractive");
60 bool doDD = settings.flag("SoftQCD:doubleDiffractive");
61 bool doCD = settings.flag("SoftQCD:centralDiffractive");
62 doMinBias = doSQ || doMB;
63 doDiffraction = doSQ || doSD || doDD || doCD;
64
65 // Separate low-mass (unresolved) and high-mass (perturbative) diffraction.
66 mMinDiff = settings.parm("Diffraction:mMinPert");
67 mWidthDiff = settings.parm("Diffraction:mWidthPert");
68 pMaxDiff = settings.parm("Diffraction:probMaxPert");
69 if (mMinDiff > infoPtr->eCM()) doDiffraction = false;
70
71 // Initialise trial shower switch
72 doTrial = useAsTrial;
73
74 // Need MPI initialization for soft QCD processes, even if only first MPI.
75 // But no need to initialize MPI if never going to use it.
76 doMPI = settings.flag("PartonLevel:MPI");
77 doMPIMB = doMPI;
78 doMPISDA = doMPI;
79 doMPISDB = doMPI;
80 doMPICD = doMPI;
81 doMPIinit = doMPI;
82 if (doMinBias || doDiffraction) doMPIinit = true;
83 if (!settings.flag("PartonLevel:all")) doMPIinit = false;
84
85 // Merging initialization.
86 bool hasMergingHooks = (mergingHooksPtr != 0);
87 doMerging = settings.flag("Merging:doUserMerging")
88 || settings.flag("Merging:doMGMerging")
89 || settings.flag("Merging:doKTMerging")
90 || settings.flag("Merging:doPTLundMerging")
91 || settings.flag("Merging:doCutBasedMerging");
92 doMerging = doMerging && hasMergingHooks;
93
94 doMergeFirstEmm = doMerging && !doTrial;
95 nTrialEmissions = 1;
96 pTLastBranch = 0.0;
97 typeLastBranch = 0;
98
99 // Flags for showers: ISR and FSR.
100 doISR = settings.flag("PartonLevel:ISR");
101 bool FSR = settings.flag("PartonLevel:FSR");
102 bool FSRinProcess = settings.flag("PartonLevel:FSRinProcess");
103 bool interleaveFSR = settings.flag("TimeShower:interleave");
104 doFSRduringProcess = FSR && FSRinProcess && interleaveFSR;
105 doFSRafterProcess = FSR && FSRinProcess && !interleaveFSR;
106 doFSRinResonances = FSR && settings.flag("PartonLevel:FSRinResonances");
107
108 // Some other flags.
109 doRemnants = settings.flag("PartonLevel:Remnants");
110 doSecondHard = settings.flag("SecondHard:generate");
111
112 // Allow R-hadron formation.
113 allowRH = settings.flag("RHadrons:allow");
114
115 // Possibility to allow user veto during evolution.
116 canVetoPT = (userHooksPtr != 0)
117 ? userHooksPtr->canVetoPT() : false;
118 pTvetoPT = (canVetoPT)
119 ? userHooksPtr->scaleVetoPT() : -1.;
120 canVetoStep = (userHooksPtr != 0)
121 ? userHooksPtr->canVetoStep() : false;
122 nVetoStep = (canVetoStep)
123 ? userHooksPtr->numberVetoStep() : -1;
124 canVetoMPIStep = (userHooksPtr != 0)
125 ? userHooksPtr->canVetoMPIStep() : false;
126 nVetoMPIStep = (canVetoMPIStep)
127 ? userHooksPtr->numberVetoMPIStep() : -1;
128 canVetoEarly = (userHooksPtr != 0)
129 ? userHooksPtr->canVetoPartonLevelEarly() : false;
130
131 // Possibility to set maximal shower scale in resonance decays.
132 canSetScale = (userHooksPtr != 0)
133 ? userHooksPtr->canSetResonanceScale() : false;
134
135 // Done with initialization only for FSR in resonance decays.
136 if (beamAPtr == 0 || beamBPtr == 0) return true;
137
138 // Flag if lepton beams, and if non-resolved ones. May change main flags.
139 hasLeptonBeams = ( beamAPtr->isLepton() || beamBPtr->isLepton() );
140 hasPointLeptons = ( hasLeptonBeams
141 && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) );
142 if (hasLeptonBeams) {
143 doMPIMB = false;
144 doMPISDA = false;
145 doMPISDB = false;
146 doMPICD = false;
147 doMPIinit = false;
148 }
149 if (hasPointLeptons) {
150 doISR = false;
151 doRemnants = false;
152 }
153
154 // Set info and initialize the respective program elements.
155 timesPtr->init( beamAPtr, beamBPtr);
156 if (doISR) spacePtr->init( beamAPtr, beamBPtr);
157 doMPIMB = multiMB.init( doMPIinit, 0, infoPtr, settings, particleDataPtr,
158 rndmPtr, beamAPtr, beamBPtr, couplingsPtr, partonSystemsPtr, sigmaTotPtr,
159 userHooksPtr);
160 if (doSD || doDD || doSQ) doMPISDA = multiSDA.init( doMPIinit, 1, infoPtr,
161 settings, particleDataPtr, rndmPtr, beamAPtr, beamPomBPtr, couplingsPtr,
162 partonSystemsPtr, sigmaTotPtr, userHooksPtr);
163 if (doSD || doDD || doSQ) doMPISDB = multiSDB.init( doMPIinit, 2, infoPtr,
164 settings, particleDataPtr, rndmPtr, beamPomAPtr, beamBPtr, couplingsPtr,
165 partonSystemsPtr, sigmaTotPtr, userHooksPtr);
166 if (doCD || doSQ) doMPICD = multiCD.init( doMPIinit, 3, infoPtr, settings,
167 particleDataPtr, rndmPtr, beamPomAPtr, beamPomBPtr, couplingsPtr,
168 partonSystemsPtr, sigmaTotPtr, userHooksPtr);
169 remnants.init( infoPtr, settings, rndmPtr, beamAPtr, beamBPtr,
170 partonSystemsPtr);
171
172 // Succeeded, or not.
173 multiPtr = &multiMB;
174 if (doMPIinit && !doMPIMB) return false;
175 if (doMPIinit && (doSD || doDD || doSQ) && (!doMPISDA || !doMPISDB))
176 return false;
177 if (doMPIinit && (doCD || doSQ) && !doMPICD) return false;
178 if (!doMPIMB || !doMPISDA || !doMPISDB || !doMPICD) doMPI = false;
179 return true;
180}
181
182//--------------------------------------------------------------------------
183
184// Function to reset PartonLevel object for trial shower usage.
185
186void PartonLevel::resetTrial() {
187
188 // Clear input pointers.
189 partonSystemsPtr->clear();
190 beamAPtr->clear();
191 beamBPtr->clear();
192 beamHadAPtr->clear();
193 beamHadBPtr->clear();
194 beamPomAPtr->clear();
195 beamPomBPtr->clear();
196
197 // Clear last branching return values.
198 pTLastBranch = 0.0;
199 typeLastBranch = 0;
200
201}
202
203//--------------------------------------------------------------------------
204
205// Main routine to do the parton-level evolution.
206
207bool PartonLevel::next( Event& process, Event& event) {
208
209 // Current event classification.
210 isResolved = infoPtr->isResolved();
211 isResolvedA = isResolved;
212 isResolvedB = isResolved;
213 isResolvedC = isResolved;
214 isDiffA = infoPtr->isDiffractiveA();
215 isDiffB = infoPtr->isDiffractiveB();
216 isDiffC = infoPtr->isDiffractiveC();
217 isDiff = isDiffA || isDiffB || isDiffC;
218 isCentralDiff = isDiffC;
219 isDoubleDiff = isDiffA && isDiffB;
220 isSingleDiff = isDiff && !isDoubleDiff && !isCentralDiff;
221 isMinBias = infoPtr->isMinBias();
222
223 // nHardLoop counts how many hard-scattering subsystems are to be processed.
224 // Almost always 1, but elastic and low-mass diffraction gives 0, while
225 // double diffraction can give up to 2. Not to be confused with SecondHard.
226 int nHardLoop = 1;
227 if (!isResolved) nHardLoop = (isDiff) ? decideResolvedDiff( process) : 0;
228
229 // Handle unresolved subsystems. Done if no resolved ones.
230 sizeProcess = 0;
231 sizeEvent = 0;
232 if (!isResolvedA || !isResolvedB || !isResolvedC) {
233 bool physical = setupUnresolvedSys( process, event);
234 if (!physical || nHardLoop == 0) return physical;
235 sizeProcess = process.size();
236 sizeEvent = event.size();
237 }
238
239 // Number of actual branchings.
240 int nBranch = 0;
241 // Number of desired branchings, negative value means no restriction.
242 int nBranchMax = (doTrial) ? nTrialEmissions : -1;
243
244 // Big outer loop to handle up to two systems (in double diffraction),
245 // but normally one. (Not indented in following, but end clearly marked.)
246 for (int iHardLoop = 1; iHardLoop <= nHardLoop; ++iHardLoop) {
247 infoPtr->setCounter(20, iHardLoop);
248 infoPtr->setCounter(21);
249
250 // Classification of diffractive system: 1 = A, 2 = B, 3 = central.
251 iDS = 0;
252 if (isDiffA || isDiffB) iDS = (iHardLoop == 2 || !isResolvedA) ? 2 : 1;
253 if (isDiffC) iDS = 3;
254
255 // Process and event records can be out of step for diffraction.
256 if (iHardLoop == 2) {
257 sizeProcess = process.size();
258 sizeEvent = event.size();
259 partonSystemsPtr->clear();
260 }
261
262 // If you need to restore then do not throw existing diffractive system.
263 if (isDiff) {
264 event.saveSize();
265 event.saveJunctionSize();
266
267 // Allow special treatment of diffractive systems.
268 setupResolvedDiff( process);
269 }
270
271 // Prepare to do multiparton interactions; at new mass for diffraction.
272 if (doMPIinit) multiPtr->reset();
273
274 // Special case if minimum bias: do hardest interaction.
275 if (isMinBias || isDiff) {
276 multiPtr->pTfirst();
277 multiPtr->setupFirstSys( process);
278 }
279
280 // Allow up to ten tries; failure possible for beam remnants.
281 // Main cause: inconsistent colour flow at the end of the day.
282 bool physical = true;
283 int nRad = 0;
284 for (int iTry = 0; iTry < NTRY; ++ iTry) {
285 infoPtr->addCounter(21);
286 for (int i = 22; i < 32; ++i) infoPtr->setCounter(i);
287
288 // Reset flag, counters and max scales.
289 physical = true;
290 nMPI = (doSecondHard) ? 2 : 1;
291 nISR = 0;
292 nFSRinProc = 0;
293 nFSRinRes = 0;
294 nISRhard = 0;
295 nFSRhard = 0;
296 pTsaveMPI = 0.;
297 pTsaveISR = 0.;
298 pTsaveFSR = 0.;
299
300 // Identify hard interaction system for showers.
301 setupHardSys( process, event);
302
303 // Optionally check for a veto after the hardest interaction.
304 if (canVetoMPIStep) {
305 doVeto = userHooksPtr->doVetoMPIStep( 1, event);
306 // Abort event if vetoed.
307 if (doVeto) {
308 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
309 return false;
310 }
311 }
312
313 // Check matching of process scale to maximum ISR/FSR/MPI scales.
314 double Q2Fac = infoPtr->Q2Fac();
315 double Q2Ren = infoPtr->Q2Ren();
316 bool limitPTmaxISR = (doISR)
317 ? spacePtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
318 bool limitPTmaxFSR = (doFSRduringProcess)
319 ? timesPtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
320 bool limitPTmaxMPI = (doMPI) ? multiPtr->limitPTmax( event) : false;
321
322 // Set hard scale, maximum for showers and multiparton interactions,
323 double pTscale = process.scale();
324 if (doSecondHard) pTscale = max( pTscale, process.scaleSecond() );
325 double pTmaxMPI = (limitPTmaxMPI) ? pTscale : infoPtr->eCM();
326 double pTmaxISR = (limitPTmaxISR) ? spacePtr->enhancePTmax() * pTscale
327 : infoPtr->eCM();
328 double pTmaxFSR = (limitPTmaxFSR) ? timesPtr->enhancePTmax() * pTscale
329 : infoPtr->eCM();
330 double pTmax = max( pTmaxMPI, max( pTmaxISR, pTmaxFSR) );
331 pTsaveMPI = pTmaxMPI;
332 pTsaveISR = pTmaxISR;
333 pTsaveFSR = pTmaxFSR;
334
335 // Prepare the classes to begin the generation.
336 if (doMPI) multiPtr->prepare( event, pTmaxMPI);
337 if (doISR) spacePtr->prepare( 0, event, limitPTmaxISR);
338 if (doFSRduringProcess) timesPtr->prepare( 0, event, limitPTmaxFSR);
339 if (doSecondHard && doISR) spacePtr->prepare( 1, event, limitPTmaxISR);
340 if (doSecondHard && doFSRduringProcess)
341 timesPtr->prepare( 1, event, limitPTmaxFSR);
342
343 // Impact parameter has now been chosen, except for diffraction.
344 if (!isDiff) infoPtr->setImpact( multiPtr->bMPI(),
345 multiPtr->enhanceMPI(), true);
346 // Set up initial veto scale.
347 doVeto = false;
348 double pTveto = pTvetoPT;
349 typeLatest = 0;
350
351 // Begin evolution down in pT from hard pT scale.
352 do {
353 infoPtr->addCounter(22);
354 typeVetoStep = 0;
355 nRad = nISR + nFSRinProc;
356
357 // Find next pT value for FSR, MPI and ISR.
358 // Order calls to minimize time expenditure.
359 double pTgen = 0.;
360 double pTtimes = (doFSRduringProcess)
361 ? timesPtr->pTnext( event, pTmaxFSR, pTgen) : -1.;
362 pTgen = max( pTgen, pTtimes);
363 double pTmulti = (doMPI)
364 ? multiPtr->pTnext( pTmaxMPI, pTgen, event) : -1.;
365 pTgen = max( pTgen, pTmulti);
366 double pTspace = (doISR)
367 ? spacePtr->pTnext( event, pTmaxISR, pTgen, nRad) : -1.;
368 double pTnow = max( pTtimes, max( pTmulti, pTspace));
369 infoPtr->setPTnow( pTnow);
370
371 // Allow a user veto. Only do it once, so remember to change pTveto.
372 if (pTveto > 0. && pTveto > pTnow) {
373 pTveto = -1.;
374 doVeto = userHooksPtr->doVetoPT( typeLatest, event);
375 // Abort event if vetoed.
376 if (doVeto) {
377 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
378 return false;
379 }
380 }
381
382 // Do a multiparton interaction (if allowed).
383 if (pTmulti > 0. && pTmulti > pTspace && pTmulti > pTtimes) {
384 infoPtr->addCounter(23);
385 if (multiPtr->scatter( event)) {
386 typeLatest = 1;
387 ++nMPI;
388 if (canVetoMPIStep && nMPI <= nVetoMPIStep) typeVetoStep = 1;
389
390 // Update ISR and FSR dipoles.
391 if (doISR) spacePtr->prepare( nMPI - 1, event);
392 if (doFSRduringProcess) timesPtr->prepare( nMPI - 1, event);
393 }
394
395 // Set maximal scales for next pT to pick.
396 pTmaxMPI = pTmulti;
397 pTmaxISR = min( pTmulti, pTmaxISR);
398 pTmaxFSR = min( pTmulti, pTmaxFSR);
399 pTmax = pTmulti;
400 nBranch++;
401 pTLastBranch = pTmulti;
402 typeLastBranch = 1;
403 }
404
405 // Do an initial-state emission (if allowed).
406 else if (pTspace > 0. && pTspace > pTtimes) {
407 infoPtr->addCounter(24);
408 if (spacePtr->branch( event)) {
409 typeLatest = 2;
410 iSysNow = spacePtr->system();
411 ++nISR;
412 if (iSysNow == 0) ++nISRhard;
413 if (canVetoStep && iSysNow == 0 && nISRhard <= nVetoStep)
414 typeVetoStep = 2;
415
416 // Update FSR dipoles.
417 if (doFSRduringProcess) timesPtr->update( iSysNow, event);
418 nBranch++;
419 pTLastBranch = pTspace;
420 typeLastBranch = 2;
421
422 // Rescatter: it is possible for kinematics to fail, in which
423 // case we need to restart the parton level processing.
424 } else if (spacePtr->doRestart()) {
425 physical = false;
426 break;
427 }
428
429 // Set maximal scales for next pT to pick.
430 pTmaxMPI = min( pTspace, pTmaxMPI);
431 pTmaxISR = pTspace;
432 pTmaxFSR = min( pTspace, pTmaxFSR);
433 pTmax = pTspace;
434 }
435
436 // Do a final-state emission (if allowed).
437 else if (pTtimes > 0.) {
438 infoPtr->addCounter(25);
439 if (timesPtr->branch( event, true)) {
440 typeLatest = 3;
441 iSysNow = timesPtr->system();
442 ++nFSRinProc;
443 if (iSysNow == 0) ++nFSRhard;
444 if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
445 typeVetoStep = 3;
446
447 // Update ISR dipoles.
448 if (doISR) spacePtr->update( iSysNow, event);
449 nBranch++;
450 pTLastBranch = pTtimes;
451 typeLastBranch = 3;
452
453 }
454
455 // Set maximal scales for next pT to pick.
456 pTmaxMPI = min( pTtimes, pTmaxMPI);
457 pTmaxISR = min( pTtimes, pTmaxISR);
458 pTmaxFSR = pTtimes;
459 pTmax = pTtimes;
460 }
461
462 // If no pT scales above zero then nothing to be done.
463 else pTmax = 0.;
464
465 // Optionally check for a veto after the first few interactions,
466 // or after the first few emissions, ISR or FSR, in the hardest system.
467 if (typeVetoStep == 1) {
468 doVeto = userHooksPtr->doVetoMPIStep( nMPI, event);
469 } else if (typeVetoStep > 1) {
470 doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
471 nFSRhard, event);
472 }
473
474 // Abort event if vetoed.
475 if (doVeto) {
476 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
477 return false;
478 }
479
480 // Keep on evolving until nothing is left to be done.
481 if (typeLatest > 0 && typeLatest < 4)
482 infoPtr->addCounter(25 + typeLatest);
483 infoPtr->setPartEvolved( nMPI, nISR);
484
485 // Handle potential merging veto.
486 if ( doMergeFirstEmm && (nISRhard + nFSRhard == 1) ){
487 // Get number of clustering steps.
488 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(process);
489 // Get maximal number of additional jets.
490 int nJetMax = mergingHooksPtr->nMaxJets();
491 // Get merging scale value.
492 double tms = mergingHooksPtr->tms();
493 // Get merging scale in current event.
494 double tnow = 0.;
495 // Use KT/Durham merging scale definition.
496 if ( mergingHooksPtr->doKTMerging() || mergingHooksPtr->doMGMerging() )
497 tnow = mergingHooksPtr->kTms(event);
498 // Use Lund PT merging scale definition.
499 else if ( mergingHooksPtr->doPTLundMerging() )
500 tnow = mergingHooksPtr->rhoms(event, false);
501 // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition.
502 else if ( mergingHooksPtr->doCutBasedMerging() )
503 tnow = mergingHooksPtr->cutbasedms(event);
504 // Use user-defined merging scale.
505 else
506 tnow = mergingHooksPtr->tmsDefinition(event);
507
508 // Check merging veto condition.
509 if ( nSteps < nJetMax && tnow > tms ) {
510 // Reset stored weights.
511 mergingHooksPtr->setWeightCKKWL(0.);
512 infoPtr->setWeightCKKWL(0.);
513 doVeto = true;
514 }
515
516 }
517
518 // Abort event if vetoed.
519 if (doVeto) {
520 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
521 return false;
522 }
523
524 // End loop evolution down in pT from hard pT scale.
525 } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
526
527 // Do all final-state emissions if not already considered above.
528 if (doFSRafterProcess && (nBranchMax <= 0 || nBranch < nBranchMax) ) {
529
530 // Find largest scale for final partons.
531 pTmax = 0.;
532 for (int i = 0; i < event.size(); ++i)
533 if (event[i].isFinal() && event[i].scale() > pTmax)
534 pTmax = event[i].scale();
535 pTsaveFSR = pTmax;
536
537 // Prepare all subsystems for evolution.
538 for (int iSys = 0; iSys < partonSystemsPtr->sizeSys(); ++iSys)
539 timesPtr->prepare( iSys, event);
540
541 // Set up initial veto scale.
542 doVeto = false;
543 pTveto = pTvetoPT;
544
545 // Begin evolution down in pT from hard pT scale.
546 do {
547 infoPtr->addCounter(29);
548 typeVetoStep = 0;
549 double pTtimes = timesPtr->pTnext( event, pTmax, 0.);
550 infoPtr->setPTnow( pTtimes);
551
552 // Allow a user veto. Only do it once, so remember to change pTveto.
553 if (pTveto > 0. && pTveto > pTtimes) {
554 pTveto = -1.;
555 doVeto = userHooksPtr->doVetoPT( 4, event);
556 // Abort event if vetoed.
557 if (doVeto) {
558 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
559 return false;
560 }
561 }
562
563 // Do a final-state emission (if allowed).
564 if (pTtimes > 0.) {
565 infoPtr->addCounter(30);
566 if (timesPtr->branch( event, true)) {
567 iSysNow = timesPtr->system();
568 ++nFSRinProc;
569 if (iSysNow == 0) ++nFSRhard;
570 if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
571 typeVetoStep = 4;
572
573 nBranch++;
574 pTLastBranch = pTtimes;
575 typeLastBranch = 4;
576
577 }
578 pTmax = pTtimes;
579 }
580
581 // If no pT scales above zero then nothing to be done.
582 else pTmax = 0.;
583
584 // Optionally check for a veto after the first few emissions.
585 if (typeVetoStep > 0) {
586 doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
587 nFSRhard, event);
588 // Abort event if vetoed.
589 if (doVeto) {
590 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
591 return false;
592 }
593 }
594
595 // Handle potential merging veto.
596 if ( doMergeFirstEmm && (nISRhard + nFSRhard == 1) ){
597 // Get number of clustering steps.
598 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(process);
599 // Get maximal number of additional jets.
600 int nJetMax = mergingHooksPtr->nMaxJets();
601 // Get merging scale value.
602 double tms = mergingHooksPtr->tms();
603 // Get merging scale in current event.
604 double tnow = 0.;
605 // Use KT/Durham merging scale definition.
606 if ( mergingHooksPtr->doKTMerging()
607 || mergingHooksPtr->doMGMerging() )
608 tnow = mergingHooksPtr->kTms(event);
609 // Use Lund PT merging scale definition.
610 else if ( mergingHooksPtr->doPTLundMerging() )
611 tnow = mergingHooksPtr->rhoms(event, false);
612 // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition.
613 else if ( mergingHooksPtr->doCutBasedMerging() )
614 tnow = mergingHooksPtr->cutbasedms(event);
615 // Use user-defined merging scale.
616 else
617 tnow = mergingHooksPtr->tmsDefinition(event);
618
619 // Check merging veto condition.
620 if ( nSteps < nJetMax && tnow > tms ) {
621 // Reset stored weights.
622 mergingHooksPtr->setWeightCKKWL(0.);
623 infoPtr->setWeightCKKWL(0.);
624 doVeto = true;
625 }
626
627 }
628
629 // Abort event if vetoed.
630 if (doVeto) {
631 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
632 return false;
633 }
634
635 // Keep on evolving until nothing is left to be done.
636 infoPtr->addCounter(31);
637
638 } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
639 }
640
641 // Handle veto after ISR + FSR + MPI, but before beam remnants
642 // and resonance decays, e.g. for MLM matching.
643 if (canVetoEarly && userHooksPtr->doVetoPartonLevelEarly( event)) {
644 doVeto = true;
645 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
646 return false;
647 }
648
649 // Add beam remnants, including primordial kT kick and colour tracing.
650 if (!doTrial && physical && doRemnants && !remnants.add( event))
651 physical = false;
652
653 // If no problems then done.
654 if (physical) break;
655
656 // Else restore and loop, but do not throw existing diffractive system.
657 if (!isDiff) event.clear();
658 else {
659 event.restoreSize();
660 event.restoreJunctionSize();
661 }
662 beamAPtr->clear();
663 beamBPtr->clear();
664 partonSystemsPtr->clear();
665
666 // End loop over ten tries. Restore from diffraction. Hopefully it worked.
667 }
668 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
669 if (!physical) return false;
670
671 // End big outer loop to handle two systems in double diffraction.
672 }
673
674 // Perform showers in resonance decay chains.
675 if(nBranchMax <= 0 || nBranch < nBranchMax)
676 doVeto = !resonanceShowers( process, event, true);
677 // Abort event if vetoed.
678 if (doVeto) return false;
679
680 // Store event properties. Not available for diffraction.
681 if (!isDiff) infoPtr->setEvolution( pTsaveMPI, pTsaveISR, pTsaveFSR,
682 nMPI, nISR, nFSRinProc, nFSRinRes);
683 if (isDiff) {
684 multiPtr->setEmpty();
685 infoPtr->setImpact( multiPtr->bMPI(), multiPtr->enhanceMPI(), false);
686 }
687
688 // Done.
689 return true;
690
691}
692
693//--------------------------------------------------------------------------
694
695// Decide which diffractive subsystems are resolved (= perturbative).
696
697int PartonLevel::decideResolvedDiff( Event& process) {
698
699 // Loop over two systems.
700 int nHighMass = 0;
701 int iDSmin = (isDiffC) ? 3 : 1;
702 int iDSmax = (isDiffC) ? 3 : 2;
703 for (int iDSnow = iDSmin; iDSnow <= iDSmax; ++iDSnow) {
704 int iDiffMot = iDSnow + 2;
705
706 // Only high-mass diffractive systems should be resolved.
707 double mDiff = process[iDiffMot].m();
708 bool isHighMass = ( mDiff > mMinDiff && rndmPtr->flat()
709 < pMaxDiff * ( 1. - exp( -(mDiff - mMinDiff) / mWidthDiff ) ) );
710
711 // Set outcome and done.
712 if (isHighMass) ++nHighMass;
713 if (iDSnow == 1) isResolvedA = isHighMass;
714 if (iDSnow == 2) isResolvedB = isHighMass;
715 if (iDSnow == 3) isResolvedC = isHighMass;
716 }
717 return nHighMass;
718
719}
720
721//--------------------------------------------------------------------------
722
723// Set up an unresolved process, i.e. elastic or diffractive.
724
725bool PartonLevel::setupUnresolvedSys( Event& process, Event& event) {
726
727 // No hard scale in event.
728 process.scale( 0.);
729
730 // Copy particles from process to event.
731 for (int i = 0; i < process.size(); ++ i) event.append( process[i]);
732
733 // Loop to find diffractively excited beams.
734 for (iDS = 1; iDS < 4; ++iDS)
735 if ( (iDS == 1 && isDiffA && !isResolvedA)
736 || (iDS == 2 && isDiffB && !isResolvedB)
737 || (iDS == 3 && isDiffC && !isResolvedC) ) {
738 int iBeam = iDS + 2;
739
740 // Diffractive mass. Boost and rotation from diffractive system
741 // rest frame, aligned along z axis, to event cm frame.
742 double mDiff = process[iBeam].m();
743 double m2Diff = mDiff * mDiff;
744 Vec4 pDiffA = (iDS == 1) ? process[1].p()
745 : process[1].p() - process[3].p();
746 Vec4 pDiffB = (iDS == 2) ? process[2].p()
747 : process[2].p() - process[4].p();
748 RotBstMatrix MtoCM;
749 MtoCM.fromCMframe( pDiffA, pDiffB);
750
751 // Beam Particle used for flavour content kicked out by Pomeron.
752 // Randomize for central diffraction; misses closed gluon loop case.
753 bool beamSideA = (iDS == 1 || (iDS == 3 && rndmPtr->flat() < 0.5));
754 BeamParticle* beamPtr = (beamSideA) ? beamAPtr : beamBPtr;
755 if (iDS == 3) beamPtr = (beamSideA) ? beamPomAPtr : beamPomBPtr;
756
757 // Pick quark or gluon kicked out and flavour subdivision.
758 beamPtr->newValenceContent();
759 bool gluonIsKicked = beamPtr->pickGluon(mDiff);
760 int id1 = beamPtr->pickValence();
761 int id2 = beamPtr->pickRemnant();
762
763 // Find flavour masses. Scale them down if too big.
764 double m1 = particleDataPtr->constituentMass(id1);
765 double m2 = particleDataPtr->constituentMass(id2);
766 if (m1 + m2 > 0.5 * mDiff) {
767 double reduce = 0.5 * mDiff / (m1 + m2);
768 m1 *= reduce;
769 m2 *= reduce;
770 }
771
772 // If quark is kicked out, then trivial kinematics in rest frame.
773 if (!gluonIsKicked) {
774 double pAbs = sqrt( pow2(m2Diff - m1*m1 - m2*m2)
775 - pow2(2. * m1 * m2) ) / (2. * mDiff);
776 if (!beamSideA) pAbs = -pAbs;
777 double e1 = (m2Diff + m1*m1 - m2*m2) / (2. * mDiff);
778 double e2 = (m2Diff + m2*m2 - m1*m1) / (2. * mDiff);
779 Vec4 p1( 0., 0., -pAbs, e1);
780 Vec4 p2( 0., 0., pAbs, e2);
781
782 // Boost and rotate to event cm frame.
783 p1.rotbst( MtoCM);
784 p2.rotbst( MtoCM);
785
786 // Set colours.
787 int col1, acol1, col2, acol2;
788 if (particleDataPtr->colType(id1) == 1) {
789 col1 = event.nextColTag();
790 acol1 = 0;
791 col2 = 0;
792 acol2 = col1;
793 } else {
794 col1 = 0;
795 acol1 = event.nextColTag();
796 col2 = acol1;
797 acol2 = 0;
798 }
799 // Update process colours to stay in step.
800 process.nextColTag();
801
802 // Store partons of diffractive system and mark system decayed.
803 int iDauBeg = event.append( id1, 23, iBeam, 0, 0, 0, col1, acol1,
804 p1, m1);
805 int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
806 p2, m2);
807 event[iBeam].statusNeg();
808 event[iBeam].daughters(iDauBeg, iDauEnd);
809
810 // If gluon is kicked out: share momentum between two remnants.
811 } else {
812 double m2Sys, zSys, pxSys, pySys, mTS1, mTS2;
813 zSys = beamPtr->zShare(mDiff, m1, m2);
814
815 // Provide relative pT kick in remnant. Construct (transverse) masses.
816 pxSys = beamPtr->pxShare();
817 pySys = beamPtr->pyShare();
818 mTS1 = m1*m1 + pxSys*pxSys + pySys*pySys;
819 mTS2 = m2*m2 + pxSys*pxSys + pySys*pySys;
820 m2Sys = mTS1 / zSys + mTS2 / (1. - zSys);
821
822 // Momentum of kicked-out massless gluon in diffractive rest frame.
823 double pAbs = (m2Diff - m2Sys) / (2. * mDiff);
824 double pLRem = (beamSideA) ? pAbs : -pAbs;
825 Vec4 pG( 0., 0., -pLRem, pAbs);
826 Vec4 pRem(0., 0., pLRem, mDiff - pAbs);
827
828 // Momenta of the two beam remnant flavours. (Lightcone p+ = m_diff!)
829 double e1 = 0.5 * (zSys * mDiff + mTS1 / (zSys * mDiff));
830 double pL1 = 0.5 * (zSys * mDiff - mTS1 / (zSys * mDiff));
831 if (!beamSideA) pL1 = -pL1;
832 Vec4 p1(pxSys, pySys, pL1, e1);
833 Vec4 p2 = pRem - p1;
834
835 // Boost and rotate to event cm frame.
836 pG.rotbst( MtoCM);
837 p1.rotbst( MtoCM);
838 p2.rotbst( MtoCM);
839
840 // Set colours.
841 int colG, acolG, col1, acol1, col2, acol2;
842 if (particleDataPtr->colType(id1) == 1) {
843 col1 = event.nextColTag();
844 acol1 = 0;
845 colG = event.nextColTag();
846 acolG = col1;
847 col2 = 0;
848 acol2 = colG;
849 } else {
850 col1 = 0;
851 acol1 = event.nextColTag();
852 colG = acol1;
853 acolG = event.nextColTag();
854 col2 = acolG;
855 acol2 = 0;
856 }
857 // Update process colours to stay in step.
858 process.nextColTag();
859 process.nextColTag();
860
861 // Store partons of diffractive system and mark system decayed.
862 int iDauBeg = event.append( 21, 23, iBeam, 0, 0, 0, colG, acolG,
863 pG, m1);
864 event.append( id1, 63, iBeam, 0, 0, 0, col1, acol1, p1, m1);
865 int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
866 p2, m2);
867 event[iBeam].statusNeg();
868 event[iBeam].daughters(iDauBeg, iDauEnd);
869 }
870
871 // End loop over beams. Done.
872 }
873 return true;
874
875}
876
877//--------------------------------------------------------------------------
878
879// Set up the hard process(es), excluding subsequent resonance decays.
880
881void PartonLevel::setupHardSys( Event& process, Event& event) {
882
883 // Incoming partons to hard process are stored in slots 3 and 4.
884 int inS = 0;
885 int inP = 3;
886 int inM = 4;
887
888 // Mother and last entry of diffractive system. Offset.
889 int iDiffMot = iDS + 2;
890 int iDiffDau = process.size() - 1;
891 int nOffset = sizeEvent - sizeProcess;
892
893 // Resolved diffraction means more entries.
894 if (isDiff) {
895 inS = iDiffMot;
896 inP = iDiffDau - 3;
897 inM = iDiffDau - 2;
898
899 // Diffractively excited particle described as Pomeron-hadron beams.
900 event[inS].statusNeg();
901 event[inS].daughters( inP - 2 + nOffset, inM - 2 + nOffset);
902 }
903
904 // If two hard interactions then find where second begins.
905 int iBeginSecond = process.size();
906 if (doSecondHard) {
907 iBeginSecond = 5;
908 while (process[iBeginSecond].status() != -21) ++iBeginSecond;
909 }
910
911 // If incoming partons are massive then recalculate to put them massless.
912 if (process[inP].m() != 0. || process[inM].m() != 0.) {
913 double pPos = process[inP].pPos() + process[inM].pPos();
914 double pNeg = process[inP].pNeg() + process[inM].pNeg();
915 process[inP].pz( 0.5 * pPos);
916 process[inP].e( 0.5 * pPos);
917 process[inP].m( 0.);
918 process[inM].pz(-0.5 * pNeg);
919 process[inM].e( 0.5 * pNeg);
920 process[inM].m( 0.);
921 }
922
923 // Add incoming hard-scattering partons to list in beam remnants.
924 double x1 = process[inP].pPos() / process[inS].m();
925 beamAPtr->append( inP + nOffset, process[inP].id(), x1);
926 double x2 = process[inM].pNeg() / process[inS].m();
927 beamBPtr->append( inM + nOffset, process[inM].id(), x2);
928
929 // Scale. Find whether incoming partons are valence or sea. Store.
930 // When an x-dependent matter profile is used with minBias,
931 // trial interactions mean that the valence/sea choice has already
932 // been made and should be restored here.
933 double scale = process.scale();
934 int vsc1, vsc2;
935 beamAPtr->xfISR( 0, process[inP].id(), x1, scale*scale);
936 if (isMinBias && (vsc1 = multiPtr->getVSC1()) != 0)
937 (*beamAPtr)[0].companion(vsc1);
938 else vsc1 = beamAPtr->pickValSeaComp();
939 beamBPtr->xfISR( 0, process[inM].id(), x2, scale*scale);
940 if (isMinBias && (vsc2 = multiPtr->getVSC2()) != 0)
941 (*beamBPtr)[0].companion(vsc2);
942 else vsc2 = beamBPtr->pickValSeaComp();
943 bool isVal1 = (vsc1 == -3);
944 bool isVal2 = (vsc2 == -3);
945 infoPtr->setValence( isVal1, isVal2);
946
947 // Initialize info needed for subsequent sequential decays + showers.
948 nHardDone = sizeProcess;
949 iPosBefShow.resize( process.size() );
950 fill( iPosBefShow.begin(), iPosBefShow.end(), 0);
951
952 // Add the beam and hard subprocess partons to the event record.
953 for (int i = sizeProcess; i < iBeginSecond; ++i) {
954 if (process[i].mother1() > inM) break;
955 int j = event.append(process[i]);
956 iPosBefShow[i] = i;
957
958 // Offset history if process and event not in step.
959 if (nOffset != 0) {
960 int iOrd = i - iBeginSecond + 7;
961 if (iOrd == 1 || iOrd == 2)
962 event[j].offsetHistory( 0, 0, 0, nOffset);
963 else if (iOrd == 3 || iOrd == 4)
964 event[j].offsetHistory( 0, nOffset, 0, nOffset);
965 else if (iOrd == 5 || iOrd == 6)
966 event[j].offsetHistory( 0, nOffset, 0, 0);
967 }
968
969 // Currently outgoing ones should not count as decayed.
970 if (event[j].status() == -22) {
971 event[j].statusPos();
972 event[j].daughters(0, 0);
973 }
974
975 // Complete task of copying hard subsystem into event record.
976 ++nHardDone;
977 }
978
979 // Store participating partons as first set in list of all systems.
980 partonSystemsPtr->addSys();
981 partonSystemsPtr->setInA(0, inP + nOffset);
982 partonSystemsPtr->setInB(0, inM + nOffset);
983 for (int i = inM + 1; i < nHardDone; ++i)
984 partonSystemsPtr->addOut(0, i + nOffset);
985 partonSystemsPtr->setSHat( 0,
986 (event[inP + nOffset].p() + event[inM + nOffset].p()).m2Calc() );
987 partonSystemsPtr->setPTHat( 0, scale);
988
989 // Identify second hard process where applicable.
990 // Since internally generated incoming partons are guaranteed massless.
991 if (doSecondHard) {
992 int inP2 = iBeginSecond;
993 int inM2 = iBeginSecond + 1;
994
995 // Add incoming hard-scattering partons to list in beam remnants.
996 // Not valid if not in rest frame??
997 x1 = process[inP2].pPos() / process[0].e();
998 beamAPtr->append( inP2, process[inP2].id(), x1);
999 x2 = process[inM2].pNeg() / process[0].e();
1000 beamBPtr->append( inM2, process[inM2].id(), x2);
1001
1002 // Find whether incoming partons are valence or sea.
1003 scale = process.scaleSecond();
1004 beamAPtr->xfISR( 1, process[inP2].id(), x1, scale*scale);
1005 beamAPtr->pickValSeaComp();
1006 beamBPtr->xfISR( 1, process[inM2].id(), x2, scale*scale);
1007 beamBPtr->pickValSeaComp();
1008
1009 // Add the beam and hard subprocess partons to the event record.
1010 for (int i = inP2; i < process.size(); ++ i) {
1011 int mother = process[i].mother1();
1012 if ( (mother > 2 && mother < inP2) || mother > inM2 ) break;
1013 event.append(process[i]);
1014 iPosBefShow[i] = i;
1015
1016 // Currently outgoing ones should not count as decayed.
1017 if (event[i].status() == -22) {
1018 event[i].statusPos();
1019 event[i].daughters(0, 0);
1020 }
1021
1022 // Complete task of copying hard subsystem into event record.
1023 ++nHardDone;
1024 }
1025
1026 // Store participating partons as second set in list of all systems.
1027 partonSystemsPtr->addSys();
1028 partonSystemsPtr->setInA(1, inP2);
1029 partonSystemsPtr->setInB(1, inM2);
1030 for (int i = inM2 + 1; i < nHardDone; ++i)
1031 partonSystemsPtr->addOut(1, i);
1032 partonSystemsPtr->setSHat( 1,
1033 (event[inP2].p() + event[inM2].p()).m2Calc() );
1034 partonSystemsPtr->setPTHat( 1, scale);
1035
1036 // End code for second hard process.
1037 }
1038
1039 // Update event colour tag to maximum in whole process.
1040 int maxColTag = 0;
1041 for (int i = 0; i < process.size(); ++ i) {
1042 if (process[i].col() > maxColTag) maxColTag = process[i].col();
1043 if (process[i].acol() > maxColTag) maxColTag = process[i].acol();
1044 }
1045 event.initColTag(maxColTag);
1046
1047 // Copy junctions from process to event.
1048 for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1049 // Resonance decay products may not have been copied from process to
1050 // event yet. If so, do not add junctions associated with decays yet.
1051 int kindJunction = process.kindJunction(iJun);
1052 bool doCopy = true;
1053 // For junction types <= 4, check if final-state legs were copied.
1054 if (kindJunction <= 4) {
1055 int iLegF1 = (kindJunction - 1) / 2;
1056 for (int iLeg = iLegF1; iLeg <= 2; ++iLeg) {
1057 bool colFound = false;
1058 for (int i = inM + 1; i < event.size(); ++i) {
1059 int col = (kindJunction % 2 == 1) ? event[i].col() : event[i].acol();
1060 if (col == process.colJunction(iJun,iLeg)) colFound = true;
1061 }
1062 if (!colFound) doCopy = false;
1063 }
1064 }
1065 if (doCopy) event.appendJunction( process.getJunction(iJun));
1066 }
1067
1068 // Done.
1069}
1070
1071//--------------------------------------------------------------------------
1072
1073// Set up the event for subsequent resonance decays and showers.
1074
1075void PartonLevel::setupShowerSys( Event& process, Event& event) {
1076
1077 // Reset event record to only contain line 0.
1078 event.clear();
1079 event.append( process[0]);
1080
1081 // Initialize info needed for subsequent sequential decays + showers.
1082 nHardDone = 1;
1083 iPosBefShow.resize( process.size());
1084 fill( iPosBefShow.begin(), iPosBefShow.end(), 0);
1085
1086 // Add the hard subprocess partons to the event record.
1087 for (int i = 1; i < process.size(); ++i) {
1088 if (process[i].mother1() > 0) break;
1089 int j = event.append(process[i]);
1090 iPosBefShow[i] = i;
1091
1092 // Currently outgoing ones should not count as decayed.
1093 if (event[j].status() == -22) {
1094 event[j].statusPos();
1095 event[j].daughters(0, 0);
1096 }
1097
1098 // Complete task of copying hard subsystem into event record.
1099 ++nHardDone;
1100 }
1101
1102 // Store participating partons as first set in list of all systems.
1103 partonSystemsPtr->clear();
1104 partonSystemsPtr->addSys();
1105 for (int i = 1; i < nHardDone; ++i) partonSystemsPtr->addOut(0, i);
1106 partonSystemsPtr->setSHat( 0, pow2(process[0].m()) );
1107 partonSystemsPtr->setPTHat( 0, 0.5 * process[0].m() );
1108
1109 // Copy junctions from process to event.
1110 for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1111 // Resonance decay products may not have been copied from process to
1112 // event yet. If so, do not add junctions associated with decays yet.
1113 int kindJunction = process.kindJunction(iJun);
1114 bool doCopy = true;
1115 // For junction types <= 4, check if final-state legs were copied.
1116 if (kindJunction <= 4) {
1117 int iLegF1 = (kindJunction - 1) / 2;
1118 for (int iLeg = iLegF1; iLeg <= 2; ++iLeg) {
1119 bool colFound = false;
1120 for (int i = 1; i < event.size(); ++i) {
1121 int col = (kindJunction % 2 == 1) ? event[i].col() : event[i].acol();
1122 if (col == process.colJunction(iJun,iLeg)) colFound = true;
1123 }
1124 if (!colFound) doCopy = false;
1125 }
1126 }
1127 if (doCopy) event.appendJunction( process.getJunction(iJun));
1128 }
1129
1130 // Done.
1131}
1132
1133//--------------------------------------------------------------------------
1134
1135// Resolved diffraction: replace full event with diffractive subsystem.
1136
1137void PartonLevel::setupResolvedDiff( Event& process) {
1138
1139 // Mother and last entry of diffractive system.
1140 int iDiffMot = iDS + 2;
1141 int iDiffDau = process.size() - 1;
1142
1143 // Diffractively excited particle to be replaced by Pomeron-hadron beams
1144 // (or Pomeron-Pomeron beams for central diffraction).
1145 process[iDiffMot].statusNeg();
1146 process[iDiffMot].daughters( iDiffDau + 1, iDiffDau + 2);
1147
1148 // Diffractive system mass.
1149 double mDiff = process[iDiffMot].m();
1150 double m2Diff = mDiff * mDiff;
1151
1152 // Set up Pomeron-proton or Pomeron-Pomeron system as if it were
1153 // the complete collision. Set Pomeron "beam particle" massless.
1154 int idDiffA = (iDS == 1) ? process[1].id() : 990;
1155 int idDiffB = (iDS == 2) ? process[2].id() : 990;
1156 double mDiffA = (iDS == 1) ? process[1].m() : 0.;
1157 double mDiffB = (iDS == 2) ? process[2].m() : 0.;
1158 double m2DiffA = mDiffA * mDiffA;
1159 double m2DiffB = mDiffB * mDiffB;
1160 double eDiffA = 0.5 * (m2Diff + m2DiffA - m2DiffB) / mDiff;
1161 double eDiffB = 0.5 * (m2Diff + m2DiffB - m2DiffA) / mDiff;
1162 double pzDiff = 0.5 * sqrtpos( pow2(m2Diff - m2DiffA - m2DiffB)
1163 - 4. * m2DiffA * m2DiffB ) / mDiff;
1164 process.append( idDiffA, 13, iDiffMot, 0, 0, 0, 0, 0,
1165 0., 0., pzDiff, eDiffA, mDiffA);
1166 process.append( idDiffB, 13, iDiffMot, 0, 0, 0, 0, 0,
1167 0., 0., -pzDiff, eDiffB, mDiffB);
1168
1169 // Reassign beam pointers to refer to subsystem effective beams.
1170 beamAPtr = (iDS == 1) ? beamHadAPtr : beamPomAPtr;
1171 beamBPtr = (iDS == 2) ? beamHadBPtr : beamPomBPtr;
1172
1173 // Pretend that the diffractive system is the whole collision.
1174 eCMsave = infoPtr->eCM();
1175 infoPtr->setECM( mDiff);
1176 beamAPtr->newPzE( pzDiff, eDiffA);
1177 beamBPtr->newPzE( -pzDiff, eDiffB);
1178
1179 // Beams not found in normal slots 1 and 2.
1180 int beamOffset = (sizeEvent > 0) ? sizeEvent - 1 : 4;
1181
1182 // Reassign beam pointers in other classes.
1183 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1184 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1185 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, iDS);
1186
1187 // Reassign multiparton interactions pointer to right object.
1188 if (iDS == 1) multiPtr = &multiSDA;
1189 else if (iDS == 2) multiPtr = &multiSDB;
1190 else multiPtr = &multiCD;
1191
1192}
1193
1194//--------------------------------------------------------------------------
1195
1196// Resolved diffraction: restore to original behaviour.
1197
1198void PartonLevel::leaveResolvedDiff( int iHardLoop, Event& process,
1199 Event& event) {
1200
1201 // Reconstruct boost and rotation to event cm frame.
1202 Vec4 pDiffA = (iDS == 1) ? process[1].p()
1203 : process[1].p() - process[3].p();
1204 Vec4 pDiffB = (iDS == 2) ? process[2].p()
1205 : process[2].p() - process[4].p();
1206 RotBstMatrix MtoCM;
1207 MtoCM.fromCMframe( pDiffA, pDiffB);
1208
1209 // Perform rotation and boost on diffractive system.
1210 for (int i = sizeProcess; i < process.size(); ++i)
1211 process[i].rotbst( MtoCM);
1212 int iFirst = (iHardLoop == 1) ? 5 + sizeEvent - sizeProcess : sizeEvent;
1213 if(isDiffC) iFirst = 6 + sizeEvent - sizeProcess;
1214 for (int i = iFirst; i < event.size(); ++i)
1215 event[i].rotbst( MtoCM);
1216
1217 // Restore cm energy.
1218 infoPtr->setECM( eCMsave);
1219 beamAPtr->newPzE( event[1].pz(), event[1].e());
1220 beamBPtr->newPzE( event[2].pz(), event[2].e());
1221
1222 // Restore beam pointers to incoming hadrons.
1223 beamAPtr = beamHadAPtr;
1224 beamBPtr = beamHadBPtr;
1225
1226 // Reassign beam pointers in other classes.
1227 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1228 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1229 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1230
1231 // Restore multiparton interactions pointer to default object.
1232 multiPtr = &multiMB;
1233
1234}
1235
1236//--------------------------------------------------------------------------
1237
1238// Handle showers in successive resonance decays.
1239
1240bool PartonLevel::resonanceShowers( Event& process, Event& event,
1241 bool skipForR) {
1242
1243 // Prepare to start over from beginning for R-hadron decays.
1244 if (allowRH) {
1245 if (skipForR) {
1246 nHardDoneRHad = nHardDone;
1247 inRHadDecay.resize(0);
1248 for (int i = 0; i < process.size(); ++i)
1249 inRHadDecay.push_back( false);
1250 } else nHardDone = nHardDoneRHad;
1251 }
1252
1253 // Isolate next system to be processed, if anything remains.
1254 int nRes = 0;
1255 int nFSRres = 0;
1256 // Number of desired branchings, negative value means no restriction.
1257 int nBranchMax = (doTrial) ? nTrialEmissions : -1;
1258
1259 while (nHardDone < process.size()) {
1260 ++nRes;
1261 int iBegin = nHardDone;
1262
1263 // In first call (skipForR = true) skip over resonances
1264 // that should form R-hadrons, and their daughters.
1265 if (allowRH) {
1266 if (skipForR) {
1267 bool comesFromR = false;
1268 int iTraceUp = iBegin;
1269 do {
1270 if ( rHadronsPtr->givesRHadron(process[iTraceUp].id()) )
1271 comesFromR = true;
1272 iTraceUp = process[iTraceUp].mother1();
1273 } while (iTraceUp > 0 && !comesFromR);
1274 if (comesFromR) {
1275 inRHadDecay[iBegin] = true;
1276 ++nHardDone;
1277 continue;
1278 }
1279
1280 // In optional second call (skipForR = false) process decay chains
1281 // inside R-hadrons.
1282 } else if (!inRHadDecay[iBegin]) {
1283 ++nHardDone;
1284 continue;
1285 }
1286 }
1287
1288 // Mother in hard process and in complete event (after shower).
1289 int iHardMother = process[iBegin].mother1();
1290 Particle& hardMother = process[iHardMother];
1291 int iBefMother = iPosBefShow[iHardMother];
1292 int iAftMother = event.iBotCopyId(iBefMother);
1293 // Possibly trace across intermediate R-hadron state.
1294 if (allowRH) {
1295 int iTraceRHadron = rHadronsPtr->trace( iAftMother);
1296 if (iTraceRHadron > 0) iAftMother = iTraceRHadron;
1297 }
1298 Particle& aftMother = event[iAftMother];
1299
1300 // From now on mother counts as decayed.
1301 aftMother.statusNeg();
1302
1303 // Mother can have been moved by showering (in any of previous steps),
1304 // so prepare to update colour and momentum information for system.
1305 int colBef = hardMother.col();
1306 int acolBef = hardMother.acol();
1307 int colAft = aftMother.col();
1308 int acolAft = aftMother.acol();
1309 RotBstMatrix M;
1310 M.bst( hardMother.p(), aftMother.p());
1311
1312 // Extract next partons from hard event into normal event record.
1313 vector<bool> iJunCopy;
1314 for (int i = iBegin; i < process.size(); ++i) {
1315 if (process[i].mother1() != iHardMother) break;
1316 int iNow = event.append( process[i] );
1317 iPosBefShow[i] = iNow;
1318 Particle& now = event.back();
1319
1320 // Currently outgoing ones should not count as decayed.
1321 if (now.status() == -22) {
1322 now.statusPos();
1323 now.daughters(0, 0);
1324 }
1325
1326 // Update daughter and mother information.
1327 if (i == iBegin) event[iAftMother].daughter1( iNow);
1328 else event[iAftMother].daughter2( iNow);
1329 now.mother1(iAftMother);
1330
1331 // Check if this parton carries a junction color in hard event.
1332 for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1333 if (iJun >= int(iJunCopy.size())) iJunCopy.push_back(false);
1334 int kindJunction = process.kindJunction(iJun);
1335 // Only consider junctions that can appear in decays.
1336 if (kindJunction >= 5) continue;
1337 int col = (kindJunction % 2 == 1) ? now.col() : now.acol();
1338 int iLegF1 = (kindJunction - 1) / 2;
1339 for (int iLeg = iLegF1; iLeg <= 2; ++iLeg)
1340 if (col == process.colJunction(iJun,iLeg)) iJunCopy[iJun] = true;
1341 }
1342
1343 // Update colour and momentum information.
1344 if (now.col() == colBef) now.col( colAft);
1345 if (now.acol() == acolBef) now.acol( acolAft);
1346 now.rotbst( M);
1347
1348 // Update vertex information.
1349 if (now.hasVertex()) now.vProd( aftMother.vDec() );
1350
1351 // Complete task of copying next subsystem into event record.
1352 ++nHardDone;
1353 }
1354 int iEnd = nHardDone - 1;
1355
1356 // Copy down junctions from hard event into normal event record.
1357 for (int iJun = 0; iJun < int(iJunCopy.size()); ++iJun) {
1358 if (!iJunCopy[iJun]) continue;
1359 Junction junCopy = process.getJunction(iJun);
1360 // Check for changed colors and update as necessary.
1361 for (int iLeg = 0; iLeg <= 2; ++iLeg) {
1362 int colLeg = junCopy.col(iLeg);
1363 if (colLeg == colBef) junCopy.col(iLeg, colAft);
1364 if (colLeg == acolBef) junCopy.col(iLeg, acolAft);
1365 }
1366 event.appendJunction(junCopy);
1367 }
1368
1369 // Reset pT of last branching
1370 pTLastBranch = 0.0;
1371
1372 // Do parton showers inside subsystem: maximum scale by mother mass.
1373 if (doFSRinResonances) {
1374 double pTmax = 0.5 * hardMother.m();
1375 if (canSetScale) pTmax
1376 = userHooksPtr->scaleResonance( iAftMother, event);
1377 nFSRhard = 0;
1378
1379 // Add new system, automatically with two empty beam slots.
1380 int iSys = partonSystemsPtr->addSys();
1381 partonSystemsPtr->setSHat(iSys, pow2(hardMother.m()) );
1382 partonSystemsPtr->setPTHat(iSys, 0.5 * hardMother.m() );
1383
1384 // Loop over allowed range to find all final-state particles.
1385 for (int i = iPosBefShow[iBegin]; i <= iPosBefShow[iEnd]; ++i)
1386 if (event[i].isFinal()) partonSystemsPtr->addOut( iSys, i);
1387
1388 // Let prepare routine do the setup.
1389 timesDecPtr->prepare( iSys, event);
1390
1391 // Number of actual branchings
1392 int nBranch = 0;
1393
1394 // Set up initial veto scale.
1395 doVeto = false;
1396 double pTveto = pTvetoPT;
1397
1398 // Begin evolution down in pT from hard pT scale.
1399 do {
1400 typeVetoStep = 0;
1401 double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
1402
1403 // Allow a user veto. Only do it once, so remember to change pTveto.
1404 if (pTveto > 0. && pTveto > pTtimes) {
1405 pTveto = -1.;
1406 doVeto = userHooksPtr->doVetoPT( 5, event);
1407 // Abort event if vetoed.
1408 if (doVeto) return false;
1409 }
1410
1411 // Do a final-state emission (if allowed).
1412 if (pTtimes > 0.) {
1413 if (timesDecPtr->branch( event)) {
1414 ++nFSRres;
1415 ++nFSRhard;
1416 if (canVetoStep && nFSRhard <= nVetoStep) typeVetoStep = 5;
1417
1418 nBranch++;
1419 pTLastBranch = pTtimes;
1420 typeLastBranch = 5;
1421
1422 }
1423 pTmax = pTtimes;
1424 }
1425
1426 // If no pT scales above zero then nothing to be done.
1427 else pTmax = 0.;
1428
1429 // Optionally check for a veto after the first few emissions.
1430 if (typeVetoStep > 0) {
1431 doVeto = userHooksPtr->doVetoStep( typeVetoStep, 0, nFSRhard,
1432 event);
1433 // Abort event if vetoed.
1434 if (doVeto) return false;
1435 }
1436
1437 // Handle potential merging veto.
1438 if( doMergeFirstEmm && nFSRhard == 1){
1439 // Get number of clustering steps.
1440 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(process);
1441 // Get maximal number of additional jets.
1442 int nJetMax = mergingHooksPtr->nMaxJets();
1443 // Get merging scale value.
1444 double tms = mergingHooksPtr->tms();
1445 // Get merging scale in current event.
1446 double tnow = 0.;
1447 // Use KT/Durham merging scale definition.
1448 if(mergingHooksPtr->doKTMerging() || mergingHooksPtr->doMGMerging())
1449 tnow = mergingHooksPtr->kTms(event);
1450 // Use Lund PT merging scale definition.
1451 else if(mergingHooksPtr->doPTLundMerging())
1452 tnow = mergingHooksPtr->rhoms(event, false);
1453 // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale
1454 // definition.
1455 else if(mergingHooksPtr->doCutBasedMerging())
1456 tnow = mergingHooksPtr->cutbasedms(event);
1457 // Use user-defined merging scale.
1458 else
1459 tnow = mergingHooksPtr->tmsDefinition(event);
1460
1461 // Check veto condition.
1462 if( nSteps < nJetMax && tnow > tms) {
1463 // Set stored weights to zero.
1464 mergingHooksPtr->setWeightCKKWL(0.);
1465 infoPtr->setWeightCKKWL(0.);
1466 // Now allow veto.
1467 doVeto = true;
1468 }
1469
1470 infoPtr->setWeightCKKWL(mergingHooksPtr->getWeightCKKWL());
1471 // Abort event if vetoed.
1472 if (doVeto) return false;
1473 }
1474
1475 // Keep on evolving until nothing is left to be done.
1476 } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
1477
1478 }
1479
1480 // No more systems to be processed. Set total number of emissions.
1481 }
1482 if (skipForR) nFSRinRes = nFSRres;
1483 return true;
1484
1485}
1486
1487//==========================================================================
1488
1489} // end namespace Pythia8