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