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