9 #include "Pythia8/PartonLevel.h"
23 const int PartonLevel::NTRY = 10;
29 bool PartonLevel::init( Info* infoPtrIn, Settings& settings,
30 ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
31 BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
32 BeamParticle* beamPomAPtrIn, BeamParticle* beamPomBPtrIn,
33 BeamParticle* beamGamAPtrIn, BeamParticle* beamGamBPtrIn,
34 BeamParticle* beamVMDAPtrIn, BeamParticle* beamVMDBPtrIn,
35 Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn,
36 SigmaTotal* sigmaTotPtr, TimeShower* timesDecPtrIn, TimeShower* timesPtrIn,
37 SpaceShower* spacePtrIn, RHadrons* rHadronsPtrIn, UserHooks* userHooksPtrIn,
38 MergingHooks* mergingHooksPtrIn, PartonVertex* partonVertexPtrIn,
43 particleDataPtr = particleDataPtrIn;
45 beamAPtr = beamAPtrIn;
46 beamBPtr = beamBPtrIn;
47 beamHadAPtr = beamAPtr;
48 beamHadBPtr = beamBPtr;
49 beamPomAPtr = beamPomAPtrIn;
50 beamPomBPtr = beamPomBPtrIn;
51 beamGamAPtr = beamGamAPtrIn;
52 beamGamBPtr = beamGamBPtrIn;
53 beamVMDAPtr = beamVMDAPtrIn;
54 beamVMDBPtr = beamVMDBPtrIn;
55 couplingsPtr = couplingsPtrIn;
56 partonSystemsPtr = partonSystemsPtrIn;
57 timesDecPtr = timesDecPtrIn;
58 timesPtr = timesPtrIn;
59 spacePtr = spacePtrIn;
60 rHadronsPtr = rHadronsPtrIn;
61 userHooksPtr = userHooksPtrIn;
62 mergingHooksPtr = mergingHooksPtrIn;
63 partonVertexPtr = partonVertexPtrIn;
66 bool doSQ = settings.flag(
"SoftQCD:all")
67 || settings.flag(
"SoftQCD:inelastic");
68 bool doND = settings.flag(
"SoftQCD:nonDiffractive");
69 bool doSD = settings.flag(
"SoftQCD:singleDiffractive");
70 bool doDD = settings.flag(
"SoftQCD:doubleDiffractive");
71 bool doCD = settings.flag(
"SoftQCD:centralDiffractive");
72 doNonDiff = doSQ || doND;
73 doDiffraction = doSQ || doSD || doDD || doCD;
74 doHardDiff = settings.flag(
"Diffraction:doHard");
75 sampleTypeDiff = (doHardDiff) ? settings.mode(
"Diffraction:sampleType")
79 mMinDiff = settings.parm(
"Diffraction:mMinPert");
80 mWidthDiff = settings.parm(
"Diffraction:mWidthPert");
81 pMaxDiff = settings.parm(
"Diffraction:probMaxPert");
82 if (mMinDiff > infoPtr->eCM()) doDiffraction =
false;
85 gammaMode = settings.mode(
"Photon:ProcessType");
87 beamHasGamma = settings.flag(
"PDF:lepton2gamma")
88 && (beamAPtr != 0) && (beamBPtr != 0);
91 beamAisGamma = (beamAPtr != 0) ? beamAPtr->isGamma() :
false;
92 beamBisGamma = (beamBPtr != 0) ? beamBPtr->isGamma() :
false;
93 beamAhasGamma = (beamHasGamma && beamAPtr->isLepton());
94 beamBhasGamma = (beamHasGamma && beamBPtr->isLepton());
95 beamAhasResGamma = (beamAPtr != 0) ? beamAPtr->hasResGamma() :
false;
96 beamBhasResGamma = (beamBPtr != 0) ? beamBPtr->hasResGamma() :
false;
97 beamHasResGamma = (gammaMode < 4) && beamHasGamma;
98 isGammaHadronDir =
false;
99 bool isAgamBhad = (beamAisGamma || beamAhasGamma) && beamBPtr->isHadron();
100 bool isBgamAhad = (beamBisGamma || beamBhasGamma) && beamAPtr->isHadron();
101 bool isAgamBgam = (beamAisGamma || beamAhasGamma)
102 && (beamBisGamma || beamBhasGamma);
103 bool onlyDirGamma =
false;
104 if ( gammaMode == 4 || ( gammaMode == 3 && isAgamBhad )
105 || ( gammaMode == 2 && isBgamAhad ) || ( gammaMode > 1 && isAgamBgam) )
109 showUnresGamma = settings.flag(
"Photon:showUnres");
113 doMPI = settings.flag(
"PartonLevel:MPI");
120 if (doNonDiff || doDiffraction) doMPIinit =
true;
121 if (!settings.flag(
"PartonLevel:all")) doMPIinit =
false;
124 pTmaxMatchMPI = settings.mode(
"MultipartonInteractions:pTmaxMatch");
127 doTrial = useAsTrial;
129 bool hasMergingHooks = (mergingHooksPtr != 0);
130 canRemoveEvent = !doTrial && hasMergingHooks
131 && ( mergingHooksPtr->doCKKWLMerging() || mergingHooksPtr->doNL3Merging());
132 canRemoveEmission = !doTrial && hasMergingHooks
133 && ( mergingHooksPtr->doUMEPSMerging() || mergingHooksPtr->doNL3Merging()
134 || mergingHooksPtr->doUNLOPSMerging() );
140 doISR = settings.flag(
"PartonLevel:ISR");
141 bool FSR = settings.flag(
"PartonLevel:FSR");
142 bool FSRinProcess = settings.flag(
"PartonLevel:FSRinProcess");
143 bool interleaveFSR = settings.flag(
"TimeShower:interleave");
144 doFSRduringProcess = FSR && FSRinProcess && interleaveFSR;
145 doFSRafterProcess = FSR && FSRinProcess && !interleaveFSR;
146 doFSRinResonances = FSR && settings.flag(
"PartonLevel:FSRinResonances");
149 doReconnect = settings.flag(
"ColourReconnection:reconnect");
150 reconnectMode = settings.mode(
"ColourReconnection:mode");
151 forceResonanceCR = settings.flag(
"ColourReconnection:forceResonance");
154 doRemnants = settings.flag(
"PartonLevel:Remnants");
155 doSecondHard = settings.flag(
"SecondHard:generate");
156 earlyResDec = settings.flag(
"PartonLevel:earlyResDec");
159 allowRH = settings.flag(
"RHadrons:allow");
162 canVetoPT = (userHooksPtr != 0)
163 ? userHooksPtr->canVetoPT() :
false;
164 pTvetoPT = (canVetoPT)
165 ? userHooksPtr->scaleVetoPT() : -1.;
166 canVetoStep = (userHooksPtr != 0)
167 ? userHooksPtr->canVetoStep() :
false;
168 nVetoStep = (canVetoStep)
169 ? userHooksPtr->numberVetoStep() : -1;
170 canVetoMPIStep = (userHooksPtr != 0)
171 ? userHooksPtr->canVetoMPIStep() :
false;
172 nVetoMPIStep = (canVetoMPIStep)
173 ? userHooksPtr->numberVetoMPIStep() : -1;
174 canVetoEarly = (userHooksPtr != 0)
175 ? userHooksPtr->canVetoPartonLevelEarly() :
false;
178 vetoWeakJets = settings.flag(
"WeakShower:vetoQCDjets");
179 vetoWeakDeltaR2 = pow2(settings.parm(
"WeakShower:vetoWeakDeltaR"));
182 canSetScale = (userHooksPtr != 0)
183 ? userHooksPtr->canSetResonanceScale() :
false;
186 canReconResSys = (userHooksPtr != 0)
187 ? userHooksPtr->canReconnectResonanceSystems() :
false;
190 if (beamAPtr == 0 || beamBPtr == 0)
return true;
194 if ( (beamAisGamma || beamAhasGamma) && gammaMode == 0) {
195 beamAPtr->setGammaMode(1);
196 if (beamAhasGamma) beamGamAPtr->setGammaMode(1);
198 if ( (beamBisGamma || beamBhasGamma) && gammaMode == 0) {
199 beamBPtr->setGammaMode(1);
200 if (beamBhasGamma) beamGamBPtr->setGammaMode(1);
204 if ( ( (gammaMode == 3) &&
205 ( (beamAPtr->isGamma() || beamAhasGamma) && beamBPtr->isHadron() ) )
206 || ( (gammaMode == 2) &&
207 ( (beamBPtr->isGamma() || beamBhasGamma) && beamAPtr->isHadron() ) )
208 || ( (gammaMode > 1) &&
209 ( (beamAPtr->isGamma() || beamAhasGamma)
210 && (beamBPtr->isGamma() || beamBhasGamma) ) ) )
214 hasTwoLeptonBeams = beamAPtr->isLepton() && beamBPtr->isLepton();
215 hasOneLeptonBeam = (beamAPtr->isLepton() || beamBPtr->isLepton())
216 && !hasTwoLeptonBeams;
217 hasPointLeptons = (hasOneLeptonBeam || hasTwoLeptonBeams)
218 && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved());
219 if ( (hasOneLeptonBeam || hasTwoLeptonBeams) && !beamHasResGamma ) {
227 if (hasTwoLeptonBeams && hasPointLeptons) {
234 if (beamHasResGamma) doMPIinit =
false;
235 if (beamHasResGamma && doND) doNDgamma =
true;
238 timesPtr->init( beamAPtr, beamBPtr);
239 if (doISR) spacePtr->init( beamAPtr, beamBPtr);
241 doMPIMB = multiMB.init( doMPIinit, 0, infoPtr, settings, particleDataPtr,
242 rndmPtr, beamAPtr, beamBPtr, couplingsPtr,
243 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr);
247 if (doSD || doDD || doSQ || (doHardDiff && beamBPtr->getGammaMode() < 2)) {
248 BeamParticle* tmpBeamA = (beamAhasGamma) ? beamGamAPtr : beamAPtr;
249 if (infoPtr->isVMDstateA()) tmpBeamA = beamVMDAPtr;
250 doMPISDA = multiSDA.init( !onlyDirGamma, 1, infoPtr, settings,
251 particleDataPtr, rndmPtr, tmpBeamA, beamPomBPtr, couplingsPtr,
252 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr,
253 (beamAisGamma || beamAhasGamma) );
255 if (doSD || doDD || doSQ || (doHardDiff && beamAPtr->getGammaMode() < 2)) {
256 BeamParticle* tmpBeamB = (beamBhasGamma) ? beamGamBPtr : beamBPtr;
257 if (infoPtr->isVMDstateB()) tmpBeamB = beamVMDBPtr;
258 doMPISDB = multiSDB.init( !onlyDirGamma, 2, infoPtr, settings,
259 particleDataPtr, rndmPtr, beamPomAPtr, tmpBeamB, couplingsPtr,
260 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr,
261 (beamBisGamma || beamBhasGamma) );
263 if (doCD || doSQ) doMPICD = multiCD.init( doMPIinit, 3, infoPtr, settings,
264 particleDataPtr, rndmPtr, beamPomAPtr, beamPomBPtr, couplingsPtr,
265 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr);
266 if (!remnants.init( infoPtr, settings, rndmPtr, beamAPtr, beamBPtr,
267 partonSystemsPtr, partonVertexPtr, particleDataPtr, &colourReconnection))
269 resonanceDecays.init( infoPtr, particleDataPtr, rndmPtr);
270 colourReconnection.init( infoPtr, settings, rndmPtr, particleDataPtr,
271 beamAPtr, beamBPtr, partonSystemsPtr);
272 junctionSplitting.init(infoPtr, settings, rndmPtr, particleDataPtr);
275 if (doHardDiff && (gammaMode !=4) )
276 hardDiffraction.init(infoPtr, settings, rndmPtr, beamAhasGamma ?
277 beamGamAPtr : beamAPtr, beamBhasGamma ? beamGamBPtr : beamBPtr,
278 beamPomAPtr, beamPomBPtr, sigmaTotPtr);
281 if ( beamHasResGamma && (doMPI || doNDgamma) ) {
284 if (beamAPtr->isLepton() && beamBPtr->isHadron() ) {
285 doMPIgmgm = multiGmGm.init( doMPIinit, 0, infoPtr, settings,
286 particleDataPtr, rndmPtr, beamGamAPtr, beamBPtr, couplingsPtr,
287 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr,
true);
289 }
else if (beamBPtr->isLepton() && beamAPtr->isHadron() ) {
290 doMPIgmgm = multiGmGm.init( doMPIinit, 0, infoPtr, settings,
291 particleDataPtr, rndmPtr, beamAPtr, beamGamBPtr, couplingsPtr,
292 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr,
true);
295 doMPIgmgm = multiGmGm.init( doMPIinit, 0, infoPtr, settings,
296 particleDataPtr, rndmPtr, beamGamAPtr, beamGamBPtr, couplingsPtr,
297 partonSystemsPtr, sigmaTotPtr, userHooksPtr, partonVertexPtr,
true);
304 if (doMPIinit && !doMPIMB)
return false;
305 if (doMPIinit && (doSD || doDD || doSQ) && (!doMPISDA || !doMPISDB))
307 if (doMPIinit && (doCD || doSQ) && !doMPICD)
return false;
308 if (!doMPIMB || !doMPISDA || !doMPISDB || !doMPICD) doMPI =
false;
317 void PartonLevel::resetTrial() {
320 partonSystemsPtr->clear();
323 beamHadAPtr->clear();
324 beamHadBPtr->clear();
325 beamPomAPtr->clear();
326 beamPomBPtr->clear();
327 beamGamAPtr->clear();
328 beamGamBPtr->clear();
329 beamVMDAPtr->clear();
330 beamVMDBPtr->clear();
342 bool PartonLevel::next(
Event& process,
Event& event) {
345 isResolved = infoPtr->isResolved();
346 isResolvedA = isResolved;
347 isResolvedB = isResolved;
348 isResolvedC = isResolved;
349 isDiffA = infoPtr->isDiffractiveA();
350 isDiffB = infoPtr->isDiffractiveB();
351 isDiffC = infoPtr->isDiffractiveC();
352 isDiff = isDiffA || isDiffB || isDiffC;
353 isCentralDiff = isDiffC;
354 isDoubleDiff = isDiffA && isDiffB;
355 isSingleDiff = isDiff && !isDoubleDiff && !isCentralDiff;
356 isNonDiff = infoPtr->isNonDiffractive();
364 bool doDiffCR =
false;
366 int nHardDiffLoop = 1;
374 infoPtr->setAbortPartonLevel(
false);
377 beamAhasResGamma = beamAPtr->hasResGamma();
378 beamBhasResGamma = beamBPtr->hasResGamma();
379 beamHasResGamma = beamAhasResGamma || beamBhasResGamma;
382 hasGammaA = (beamAPtr != 0) ? beamAPtr->getGammaMode() > 0 :
false;
383 hasGammaB = (beamBPtr != 0) ? beamBPtr->getGammaMode() > 0 :
false;
386 gammaModeEvent = gammaMode;
387 if ( hasGammaA || hasGammaB ) {
388 if (beamAPtr->getGammaMode() < 2 && beamBPtr->getGammaMode() < 2)
390 if (beamAPtr->getGammaMode() < 2 && beamBPtr->getGammaMode() == 2)
392 if (beamAPtr->getGammaMode() == 2 && beamBPtr->getGammaMode() < 2)
394 if (beamAPtr->getGammaMode() == 2 && beamBPtr->getGammaMode() == 2)
399 isGammaHadronDir = !(hasGammaA || hasGammaB) ?
false :
400 ( (beamAPtr->getGammaMode() == 2) && (beamBPtr->getGammaMode() == 0) )
401 || ( (beamAPtr->getGammaMode() == 0) && (beamBPtr->getGammaMode() == 2) );
405 if ( !setupResolvedLeptonGamma( process) )
return false;
417 if ( beamBPtr->isHadron() || ((beamBisGamma || beamBhasGamma)
418 && beamBPtr->getGammaMode() == 1)) {
419 double xGammaB = beamBhasGamma ? beamHadBPtr->xGamma() : 1.;
420 double xBrescaled = infoPtr->x2pdf() / xGammaB;
421 double pdfB = beamBhasGamma ? beamGamBPtr->xf(infoPtr->id2pdf(),
422 xBrescaled, infoPtr->Q2Fac()) : infoPtr->pdf2();
423 isHardDiffA = hardDiffraction.isDiffractive(2, infoPtr->id2pdf(),
424 xBrescaled, infoPtr->Q2Fac(), pdfB);
426 if ( beamAPtr->isHadron() || ((beamAisGamma || beamAhasGamma)
427 && beamAPtr->getGammaMode() == 1)) {
428 double xGammaA = beamAhasGamma ? beamHadAPtr->xGamma() : 1.;
429 double xArescaled = infoPtr->x1pdf() / xGammaA;
430 double pdfA = beamAhasGamma ? beamGamAPtr->xf(infoPtr->id1pdf(),
431 xArescaled, infoPtr->Q2Fac()) : infoPtr->pdf1();
432 isHardDiffB = hardDiffraction.isDiffractive(1, infoPtr->id1pdf(),
433 xArescaled, infoPtr->Q2Fac(), pdfA);
437 if (isHardDiffA && isHardDiffB) {
438 if (rndmPtr->flat() < 0.5) isHardDiffA =
false;
439 else isHardDiffB =
false;
441 isHardDiff = isHardDiffA || isHardDiffB;
444 double xPomA = (isHardDiffB) ? hardDiffraction.getXPomeronA() : 0.;
445 double xPomB = (isHardDiffA) ? hardDiffraction.getXPomeronB() : 0.;
446 double tPomA = (isHardDiffB) ? hardDiffraction.getTPomeronA() : 0.;
447 double tPomB = (isHardDiffA) ? hardDiffraction.getTPomeronB() : 0.;
448 infoPtr->setHardDiff(
false,
false, isHardDiffA, isHardDiffB,
449 xPomA, xPomB, tPomA, tPomB);
452 if (!isHardDiff && sampleTypeDiff > 2) {
455 if (beamHasGamma) leaveResolvedLeptonGamma( process, event,
false);
461 if (sampleTypeDiff%2 == 1) setupHardDiff( process);
463 else nHardDiffLoop = 2;
471 if (!isResolved) nHardLoop = (isDiff) ? decideResolvedDiff( process) : 0;
476 if (!isResolvedA || !isResolvedB || !isResolvedC) {
477 bool physical = setupUnresolvedSys( process, event);
478 if (!physical || nHardLoop == 0)
return physical;
479 sizeProcess = process.size();
480 sizeEvent =
event.size();
486 int nBranchMax = (doTrial) ? nTrialEmissions : -1;
489 bool hasMergingHooks = (mergingHooksPtr != 0);
490 if ( hasMergingHooks && canRemoveEvent )
491 mergingHooksPtr->storeWeights(infoPtr->getWeightCKKWL());
494 if (userHooksPtr != 0) userHooksPtr->setEnhancedEventWeight(1.);
497 for (
int iHardDiffLoop = 1; iHardDiffLoop <= nHardDiffLoop;
502 for (
int iHardLoop = 1; iHardLoop <= nHardLoop; ++iHardLoop) {
503 infoPtr->setCounter(20, iHardLoop);
504 infoPtr->setCounter(21);
508 if (isDiffA || isDiffB) iDS = (iHardLoop == 2 || !isResolvedA) ? 2 : 1;
509 if (isDiffC) iDS = 3;
512 if (iHardLoop == 2) {
513 sizeProcess = process.size();
514 sizeEvent =
event.size();
515 partonSystemsPtr->clear();
516 if (event.lastColTag() > process.lastColTag())
517 process.initColTag(event.lastColTag());
523 event.saveJunctionSize();
526 setupResolvedDiff( process);
530 if (doMPIinit) multiPtr->reset();
533 if (isNonDiff || isDiff) {
535 multiPtr->setupFirstSys( process);
540 bool physical =
true;
542 for (
int iTry = 0; iTry < NTRY; ++ iTry) {
543 infoPtr->addCounter(21);
544 for (
int i = 22; i < 32; ++i) infoPtr->setCounter(i);
548 nMPI = (doSecondHard) ? 2 : 1;
559 infoPtr->setPartEvolved(nMPI, nISR);
563 if (beamAPtr->isGamma() && ( doMPI || doISR ) ) beamAPtr->resetGamma();
564 if (beamBPtr->isGamma() && ( doMPI || doISR ) ) beamBPtr->resetGamma();
567 setupHardSys( process, event);
570 if (canVetoMPIStep) {
571 doVeto = userHooksPtr->doVetoMPIStep( 1, event);
574 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
575 if (beamHasResGamma) leaveResolvedLeptonGamma( process, event,
false);
581 double Q2Fac = infoPtr->Q2Fac();
582 double Q2Ren = infoPtr->Q2Ren();
583 bool limitPTmaxISR = (doISR)
584 ? spacePtr->limitPTmax( event, Q2Fac, Q2Ren) :
false;
585 bool limitPTmaxFSR = (doFSRduringProcess)
586 ? timesPtr->limitPTmax( event, Q2Fac, Q2Ren) :
false;
587 bool limitPTmaxMPI = (doMPI) ? multiPtr->limitPTmax( event) :
false;
590 timesPtr->prepareGlobal( event);
591 bool isFirstTrial =
true;
594 double pTscaleRad = process.scale();
595 double pTscaleMPI = (doMPI && pTmaxMatchMPI == 3)
596 ? multiPtr->scaleLimitPT() : pTscaleRad;
598 pTscaleRad = max( pTscaleRad, process.scaleSecond() );
599 pTscaleMPI = min( pTscaleMPI, process.scaleSecond() );
601 double pTmaxMPI = (limitPTmaxMPI) ? pTscaleMPI : infoPtr->eCM();
602 double pTmaxISR = (limitPTmaxISR) ? spacePtr->enhancePTmax() * pTscaleRad
604 double pTmaxFSR = (limitPTmaxFSR) ? timesPtr->enhancePTmax() * pTscaleRad
609 beamAPtr->pTMPI( process.scale() );
610 beamBPtr->pTMPI( process.scale() );
613 if ( hasMergingHooks && (doTrial || canRemoveEvent || canRemoveEmission) )
614 mergingHooksPtr->setShowerStartingScales( doTrial,
615 (canRemoveEvent || canRemoveEmission), pTscaleRad, process, pTmaxFSR,
616 limitPTmaxFSR, pTmaxISR, limitPTmaxISR, pTmaxMPI, limitPTmaxMPI );
617 double pTmax = max( pTmaxMPI, max( pTmaxISR, pTmaxFSR) );
618 pTsaveMPI = pTmaxMPI;
619 pTsaveISR = pTmaxISR;
620 pTsaveFSR = pTmaxFSR;
623 if (doMPI) multiPtr->prepare( event, pTmaxMPI, (iHardDiffLoop == 2) );
624 if (doISR) spacePtr->prepare( 0, event, limitPTmaxISR);
625 if (doFSRduringProcess) timesPtr->prepare( 0, event, limitPTmaxFSR);
626 if (doSecondHard && doISR) spacePtr->prepare( 1, event, limitPTmaxISR);
627 if (doSecondHard && doFSRduringProcess) timesPtr->prepare( 1, event,
632 if (!isDiff) infoPtr->setImpact( multiPtr->bMPI(), multiPtr->enhanceMPI(),
633 multiPtr->enhanceMPIavg(),
true, (iHardDiffLoop == 2) );
637 double pTveto = pTvetoPT;
642 infoPtr->addCounter(22);
644 nRad = nISR + nFSRinProc;
647 bool unresolvedGammaA = (beamAPtr->isGamma()
648 && !(beamAPtr->resolvedGamma()) );
649 bool unresolvedGammaB = (beamBPtr->isGamma()
650 && !(beamBPtr->resolvedGamma()) );
651 bool unresolvedGamma = (unresolvedGammaA || unresolvedGammaB)
652 || gammaModeEvent == 4;
660 if ( hasMergingHooks && doTrial)
661 pTgen = max( pTgen, mergingHooksPtr->getShowerStoppingScale() );
663 double pTtimes = (doFSRduringProcess)
664 ? timesPtr->pTnext( event, pTmaxFSR, pTgen, isFirstTrial, doTrial)
666 pTgen = max( pTgen, pTtimes);
668 double pTmulti = (doMPI && !unresolvedGamma)
669 ? multiPtr->pTnext( pTmaxMPI, pTgen, event) : -1.;
670 pTgen = max( pTgen, pTmulti);
671 double pTspace = (doISR)
672 ? spacePtr->pTnext( event, pTmaxISR, pTgen, nRad, doTrial) : -1.;
673 double pTnow = max( pTtimes, max( pTmulti, pTspace));
676 infoPtr->setPTnow( pTnow);
677 isFirstTrial =
false;
680 if (pTveto > 0. && pTveto > pTnow) {
682 doVeto = userHooksPtr->doVetoPT( typeLatest, event);
685 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
687 leaveResolvedLeptonGamma( process, event,
false);
693 if (pTmulti > 0. && pTmulti > pTspace && pTmulti > pTtimes) {
694 infoPtr->addCounter(23);
695 if (multiPtr->scatter( event)) {
698 if (canVetoMPIStep && nMPI <= nVetoMPIStep) typeVetoStep = 1;
701 if (isHardDiff && sampleTypeDiff == 4 && iHardDiffLoop == 1) {
702 infoPtr->setHardDiff(
false,
false,
false,
false, 0., 0., 0., 0.);
706 leaveResolvedLeptonGamma( process, event,
false);
711 if (doISR) spacePtr->prepare( nMPI - 1, event);
712 if (doFSRduringProcess) timesPtr->prepare( nMPI - 1, event);
714 pTLastBranch = pTmulti;
720 pTmaxISR = min(pTmulti, pTmaxISR);
721 pTmaxFSR = min(pTmulti, pTmaxFSR);
726 else if (pTspace > 0. && pTspace > pTtimes) {
727 infoPtr->addCounter(24);
730 if (spacePtr->branch( event)
731 && ( !(nMPI > 1 && spacePtr->wasGamma2qqbar()) ) ) {
733 iSysNow = spacePtr->system();
735 if (iSysNow == 0) ++nISRhard;
736 if (canVetoStep && iSysNow == 0 && nISRhard <= nVetoStep)
740 if (doFSRduringProcess) timesPtr->update( iSysNow, event,
741 spacePtr->getHasWeaklyRadiated());
743 pTLastBranch = pTspace;
748 }
else if (spacePtr->doRestart()) {
754 pTmaxMPI = min( min(pTspace,pTmaxISR), pTmaxMPI);
755 pTmaxISR = min(pTspace,pTmaxISR);
756 pTmaxFSR = min( min(pTspace,pTmaxISR), pTmaxFSR);
761 else if (pTtimes > 0.) {
762 infoPtr->addCounter(25);
763 if (timesPtr->branch( event,
true)) {
765 iSysNow = timesPtr->system();
767 if (iSysNow == 0) ++nFSRhard;
768 if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
772 if (doISR) spacePtr->update( iSysNow, event,
773 timesPtr->getHasWeaklyRadiated());
775 pTLastBranch = pTtimes;
781 pTmaxMPI = min( min(pTtimes,pTmaxFSR), pTmaxMPI);
782 pTmaxISR = min( min(pTtimes,pTmaxFSR), pTmaxISR);
783 pTmaxFSR = min(pTtimes, pTmaxFSR);
792 if ( (infoPtr->code() == 221 || infoPtr->code() == 222) &&
793 nISRhard + nFSRhard == 2 && vetoWeakJets) {
794 int id1 =
event[partonSystemsPtr->getOut(0,0)].id();
795 int id2 =
event[partonSystemsPtr->getOut(0,1)].id();
796 int id3 =
event[partonSystemsPtr->getOut(0,2)].id();
797 Vec4 p1 =
event[partonSystemsPtr->getOut(0,0)].p();
798 Vec4 p2 =
event[partonSystemsPtr->getOut(0,1)].p();
799 Vec4 p3 =
event[partonSystemsPtr->getOut(0,2)].p();
803 bool doubleCountEvent =
true;
804 if (abs(id1) == 24 || abs(id1) == 23) {
805 if (abs(id2) > 21 || abs(id3) > 21)
806 doubleCountEvent =
false;
807 }
else if (abs(id2) == 24 || abs(id2) == 23) {
811 doubleCountEvent =
false;
812 }
else if ( abs(id3) == 24 || abs(id3) == 23) {
817 if (doubleCountEvent) {
820 if (p2.pT2() < d) {d = p2.pT2(); cut =
false;}
821 if (p3.pT2() < d) {d = p3.pT2(); cut =
false;}
826 double dij = min(p1.pT2(),p2.pT2())
827 * pow2(RRapPhi(p1,p2)) / vetoWeakDeltaR2;
835 double dij = min(p1.pT2(),p3.pT2())
836 * pow2(RRapPhi(p1,p3)) / vetoWeakDeltaR2;
846 if (abs(id2) == 21 || abs(id3) == 21 || id2 == - id3) {
847 double dij = min(p2.pT2(),p3.pT2())
848 * pow2(RRapPhi(p2,p3)) / vetoWeakDeltaR2;
856 if (cut)
return false;
862 if (typeVetoStep == 1) {
863 doVeto = userHooksPtr->doVetoMPIStep( nMPI, event);
864 }
else if (typeVetoStep > 1) {
865 doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
871 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
872 if (beamHasResGamma) leaveResolvedLeptonGamma( process, event);
877 if (typeLatest > 0 && typeLatest < 4)
878 infoPtr->addCounter(25 + typeLatest);
879 if (!isDiff) infoPtr->setPartEvolved( nMPI, nISR);
882 if ( canRemoveEvent && nISRhard + nFSRhard == 1 ) {
884 mergingHooksPtr->doVetoStep( process, event );
888 }
while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
891 if (doFSRafterProcess && (nBranchMax <= 0 || nBranch < nBranchMax) ) {
895 for (
int i = 0; i <
event.size(); ++i)
896 if (event[i].isFinal() &&
event[i].scale() > pTmax)
897 pTmax = event[i].scale();
901 for (
int iSys = 0; iSys < partonSystemsPtr->sizeSys(); ++iSys)
902 timesPtr->prepare( iSys, event);
910 infoPtr->addCounter(29);
912 double pTtimes = timesPtr->pTnext( event, pTmax, 0.);
913 infoPtr->setPTnow( pTtimes);
916 if (pTveto > 0. && pTveto > pTtimes) {
918 doVeto = userHooksPtr->doVetoPT( 4, event);
921 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
922 if (beamHasResGamma) leaveResolvedLeptonGamma( process, event);
929 infoPtr->addCounter(30);
930 if (timesPtr->branch( event,
true)) {
931 iSysNow = timesPtr->system();
933 if (iSysNow == 0) ++nFSRhard;
934 if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
938 pTLastBranch = pTtimes;
949 if (typeVetoStep > 0) {
950 doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
954 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
955 if (beamHasResGamma) leaveResolvedLeptonGamma( process, event);
961 if ( canRemoveEvent && nISRhard + nFSRhard == 1 ) {
963 mergingHooksPtr->doVetoStep( process, event );
967 infoPtr->addCounter(31);
969 }
while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
974 if (canVetoEarly && userHooksPtr->doVetoPartonLevelEarly( event)) {
976 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
977 if (beamHasResGamma) leaveResolvedLeptonGamma( process, event);
983 int oldSizeEvt =
event.size();
984 int oldSizeSys = partonSystemsPtr->sizeSys();
985 if (nBranchMax <= 0 || nBranch < nBranchMax)
986 doVeto = !resonanceShowers( process, event,
true);
988 if (doVeto)
return false;
991 for (
int iSys = oldSizeSys; iSys < partonSystemsPtr->sizeSys(); ++iSys)
992 for (
int iOut = 0; iOut < partonSystemsPtr->sizeOut(iSys); ++iOut)
993 partonSystemsPtr->addOut(0, partonSystemsPtr->getOut( iSys, iOut) );
994 partonSystemsPtr->setSizeSys( oldSizeSys);
998 if (!wzDecayShowers( event))
return false;
1001 if (canReconResSys && !userHooksPtr->doReconnectResonanceSystems(
1002 oldSizeEvt, event))
return false;
1009 iFirst = (iHardLoop == 1) ? 5 + sizeEvent - sizeProcess : sizeEvent;
1010 if (isDiffC) iFirst = 6 + sizeEvent - sizeProcess;
1014 if (infoPtr->hasPomPsystem()) {
1020 if (!doTrial && physical && doRemnants
1021 && (!beamHasGamma || gammaModeEvent != 4)
1022 && !remnants.add( event, iFirst, doDiffCR)) physical =
false;
1025 if (physical)
break;
1028 if (!isDiff)
event.clear();
1030 event.restoreSize();
1031 event.restoreJunctionSize();
1035 partonSystemsPtr->clear();
1038 if (beamAhasResGamma) beamHadAPtr->clear();
1039 if (beamBhasResGamma) beamHadBPtr->clear();
1040 if (infoPtr->isVMDstateA()) beamVMDAPtr->clear();
1041 if (infoPtr->isVMDstateB()) beamVMDBPtr->clear();
1045 if (isDiff) leaveResolvedDiff( iHardLoop, process, event);
1049 if (infoPtr->hasPomPsystem()) leaveHardDiff( process, event);
1051 if (beamHasGamma) leaveResolvedLeptonGamma( process, event,
false);
1060 if (isHardDiff && sampleTypeDiff%2 == 0 && iHardDiffLoop == 1 && nMPI == 1) {
1064 partonSystemsPtr->clear();
1065 setupHardDiff( process);
1070 if (doReconnect && !doDiffCR && reconnectMode > 0) {
1071 Event eventSave = event;
1072 bool colCorrect =
false;
1073 for (
int i = 0; i < 10; ++i) {
1074 colourReconnection.next(event, 0);
1075 if (junctionSplitting.checkColours(event)) {
1079 else event = eventSave;
1082 infoPtr->errorMsg(
"Error in PartonLevel::next: "
1083 "Colour reconnection failed.");
1089 int oldSizeEvt =
event.size();
1091 if (nBranchMax <= 0 || nBranch < nBranchMax)
1092 doVeto = !resonanceShowers( process, event,
true);
1094 if (doVeto)
return false;
1098 if (!wzDecayShowers( event))
return false;
1101 if (canReconResSys && !userHooksPtr->doReconnectResonanceSystems(
1102 oldSizeEvt, event))
return false;
1106 if (!isDiff) infoPtr->setEvolution( pTsaveMPI, pTsaveISR, pTsaveFSR,
1107 nMPI, nISR, nFSRinProc, nFSRinRes);
1109 multiPtr->setEmpty();
1110 infoPtr->setImpact( multiPtr->bMPI(), multiPtr->enhanceMPI(),
1111 multiPtr->enhanceMPIavg(),
false);
1115 if (!earlyResDec && forceResonanceCR && doReconnect &&
1116 !doDiffCR && reconnectMode != 0) {
1117 Event eventSave = event;
1118 bool colCorrect =
false;
1119 for (
int i = 0; i < 10; ++i) {
1120 colourReconnection.next(event, oldSizeEvt);
1121 if (junctionSplitting.checkColours(event)) {
1125 else event = eventSave;
1128 infoPtr->errorMsg(
"Error in PartonLevel::next: "
1129 "Colour reconnection failed.");
1139 if (sampleTypeDiff == 2 && iHardDiffLoop == 1 && nMPI > 1) {
1140 infoPtr->setHardDiff(
false,
false,
false,
false, 0., 0., 0., 0.);
1145 if (infoPtr->hasPomPsystem()) leaveHardDiff( process, event);
1153 if ( ( beamAPtr->isGamma() || beamBPtr->isGamma() )
1154 && ( !beamAPtr->resolvedGamma() || !beamBPtr->resolvedGamma() ) ) {
1155 if (!showUnresGamma && (gammaModeEvent != 4) ) cleanEventFromGamma( event);
1159 if (beamHasGamma) leaveResolvedLeptonGamma( process, event,
true);
1170 int PartonLevel::decideResolvedDiff(
Event& process) {
1174 int iDSmin = (isDiffC) ? 3 : 1;
1175 int iDSmax = (isDiffC) ? 3 : 2;
1176 for (
int iDSnow = iDSmin; iDSnow <= iDSmax; ++iDSnow) {
1177 int iDiffMot = iDSnow + 2;
1180 double mDiff = process[iDiffMot].m();
1181 bool isHighMass = ( mDiff > mMinDiff && rndmPtr->flat()
1182 < pMaxDiff * ( 1. - exp( -(mDiff - mMinDiff) / mWidthDiff ) ) );
1185 if (isHighMass) ++nHighMass;
1186 if (iDSnow == 1) isResolvedA = isHighMass;
1187 if (iDSnow == 2) isResolvedB = isHighMass;
1188 if (iDSnow == 3) isResolvedC = isHighMass;
1198 bool PartonLevel::setupUnresolvedSys(
Event& process,
Event& event) {
1204 for (
int i = 0; i < process.size(); ++ i) event.append( process[i]);
1207 for (iDS = 1; iDS < 4; ++iDS)
1208 if ( (iDS == 1 && isDiffA && !isResolvedA)
1209 || (iDS == 2 && isDiffB && !isResolvedB)
1210 || (iDS == 3 && isDiffC && !isResolvedC) ) {
1211 int iBeam = iDS + 2;
1215 double mDiff = process[iBeam].m();
1216 double m2Diff = mDiff * mDiff;
1217 Vec4 pDiffA = (iDS == 1) ? process[1].p()
1218 : process[1].p() - process[3].p();
1219 Vec4 pDiffB = (iDS == 2) ? process[2].p()
1220 : process[2].p() - process[4].p();
1222 MtoCM.fromCMframe( pDiffA, pDiffB);
1226 bool beamSideA = (iDS == 1 || (iDS == 3 && rndmPtr->flat() < 0.5));
1227 BeamParticle* beamPtr = (beamSideA) ? beamAPtr : beamBPtr;
1228 if (iDS == 3) beamPtr = (beamSideA) ? beamPomAPtr : beamPomBPtr;
1231 beamPtr->newValenceContent();
1232 bool gluonIsKicked = beamPtr->pickGluon(mDiff);
1233 int id1 = beamPtr->pickValence();
1234 int id2 = beamPtr->pickRemnant();
1237 double m1 = particleDataPtr->constituentMass(id1);
1238 double m2 = particleDataPtr->constituentMass(id2);
1239 if (m1 + m2 > 0.5 * mDiff) {
1240 double reduce = 0.5 * mDiff / (m1 + m2);
1246 if (!gluonIsKicked) {
1247 double pAbs = sqrt( pow2(m2Diff - m1*m1 - m2*m2)
1248 - pow2(2. * m1 * m2) ) / (2. * mDiff);
1249 if (!beamSideA) pAbs = -pAbs;
1250 double e1 = (m2Diff + m1*m1 - m2*m2) / (2. * mDiff);
1251 double e2 = (m2Diff + m2*m2 - m1*m1) / (2. * mDiff);
1252 Vec4 p1( 0., 0., -pAbs, e1);
1253 Vec4 p2( 0., 0., pAbs, e2);
1260 int col1, acol1, col2, acol2;
1261 if (particleDataPtr->colType(id1) == 1) {
1262 col1 =
event.nextColTag();
1268 acol1 =
event.nextColTag();
1273 process.nextColTag();
1276 int iDauBeg =
event.append( id1, 24, iBeam, 0, 0, 0, col1, acol1,
1278 int iDauEnd =
event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
1280 event[iBeam].statusNeg();
1281 event[iBeam].daughters(iDauBeg, iDauEnd);
1285 double m2Sys, zSys, pxSys, pySys, mTS1, mTS2;
1286 zSys = beamPtr->zShare(mDiff, m1, m2);
1289 pxSys = beamPtr->pxShare();
1290 pySys = beamPtr->pyShare();
1291 mTS1 = m1*m1 + pxSys*pxSys + pySys*pySys;
1292 mTS2 = m2*m2 + pxSys*pxSys + pySys*pySys;
1293 m2Sys = mTS1 / zSys + mTS2 / (1. - zSys);
1296 double pAbs = (m2Diff - m2Sys) / (2. * mDiff);
1297 double pLRem = (beamSideA) ? pAbs : -pAbs;
1298 Vec4 pG( 0., 0., -pLRem, pAbs);
1299 Vec4 pRem(0., 0., pLRem, mDiff - pAbs);
1302 double e1 = 0.5 * (zSys * mDiff + mTS1 / (zSys * mDiff));
1303 double pL1 = 0.5 * (zSys * mDiff - mTS1 / (zSys * mDiff));
1304 if (!beamSideA) pL1 = -pL1;
1305 Vec4 p1(pxSys, pySys, pL1, e1);
1306 Vec4 p2 = pRem - p1;
1315 int colG, acolG, col1, acol1, col2, acol2;
1316 if (particleDataPtr->colType(id1) == 1) {
1317 col1 =
event.nextColTag();
1319 colG =
event.nextColTag();
1325 acol1 =
event.nextColTag();
1327 acolG =
event.nextColTag();
1332 process.nextColTag();
1333 process.nextColTag();
1336 int iDauBeg =
event.append( 21, 24, iBeam, 0, 0, 0, colG, acolG, pG, 0.);
1337 event.append( id1, 63, iBeam, 0, 0, 0, col1, acol1, p1, m1);
1338 int iDauEnd =
event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
1340 event[iBeam].statusNeg();
1341 event[iBeam].daughters(iDauBeg, iDauEnd);
1354 void PartonLevel::setupHardSys(
Event& process,
Event& event) {
1362 int iDiffMot = iDS + 2;
1363 int iDiffDau = process.size() - 1;
1364 int nOffset = sizeEvent - sizeProcess;
1367 if (infoPtr->hasPomPsystem()) {
1368 iDiffMot = (isHardDiffB) ? 4 : 3;
1381 event[inS].statusNeg();
1382 event[inS].daughters( inP - 2 + nOffset, inM - 2 + nOffset);
1387 int nGammaOffset = 0;
1388 int nGammaDiffOffset = 0;
1389 if ( beamHasResGamma || (beamHasGamma && isGammaHadronDir) ) {
1391 inP += nGammaOffset;
1392 inM += nGammaOffset;
1393 if ( (isHardDiffA || isHardDiffB) && hardDiffSet ) nGammaDiffOffset = 4;
1397 int iBeginSecond = process.size();
1400 while (process[iBeginSecond].status() != -21) ++iBeginSecond;
1404 if (process[inP].m() != 0. || process[inM].m() != 0.) {
1405 double pPos = process[inP].pPos() + process[inM].pPos();
1406 double pNeg = process[inP].pNeg() + process[inM].pNeg();
1407 process[inP].pz( 0.5 * pPos);
1408 process[inP].e( 0.5 * pPos);
1409 process[inP].m( 0.);
1410 process[inM].pz(-0.5 * pNeg);
1411 process[inM].e( 0.5 * pNeg);
1412 process[inM].m( 0.);
1416 double x1 = process[inP].pPos() / process[inS].m();
1417 double x2 = process[inM].pNeg() / process[inS].m();
1420 if ( beamHasResGamma || (beamHasGamma && isGammaHadronDir) ) {
1422 beamHadAPtr->append( 3, process[3].
id(), beamHadAPtr->xGamma() );
1424 beamHadBPtr->append( 4, process[4].
id(), beamHadBPtr->xGamma() );
1425 x1 = process[inP].pPos() / ( process[3 + nGammaDiffOffset].p()
1426 + process[4 + nGammaDiffOffset].p() ).mCalc();
1427 x2 = process[inM].pNeg() / ( process[3 + nGammaDiffOffset].p()
1428 + process[4 + nGammaDiffOffset].p() ).mCalc();
1430 beamAPtr->append( inP + nOffset, process[inP].
id(), x1);
1431 beamBPtr->append( inM + nOffset, process[inM].
id(), x2);
1437 double scale = process.scale();
1439 beamAPtr->xfISR( 0, process[inP].
id(), x1, scale*scale);
1440 if (isNonDiff && (vsc1 = multiPtr->getVSC1()) != 0)
1441 (*beamAPtr)[0].companion(vsc1);
1442 else if (beamAPtr->isGamma()) vsc1 = beamAPtr->gammaValSeaComp(0);
1443 else vsc1 = beamAPtr->pickValSeaComp();
1444 beamBPtr->xfISR( 0, process[inM].
id(), x2, scale*scale);
1445 if (isNonDiff && (vsc2 = multiPtr->getVSC2()) != 0)
1446 (*beamBPtr)[0].companion(vsc2);
1447 else if (beamBPtr->isGamma()) vsc2 = beamBPtr->gammaValSeaComp(0);
1448 else vsc2 = beamBPtr->pickValSeaComp();
1449 bool isVal1 = (vsc1 == -3);
1450 bool isVal2 = (vsc2 == -3);
1451 infoPtr->setValence( isVal1, isVal2);
1454 nHardDone = sizeProcess;
1455 iPosBefShow.resize( process.size() );
1456 fill( iPosBefShow.begin(), iPosBefShow.end(), 0);
1459 for (
int i = sizeProcess; i < iBeginSecond; ++i) {
1460 if (process[i].mother1() > inM)
break;
1461 int j =
event.append(process[i]);
1466 int iOrd = i - iBeginSecond + 7;
1467 if (iOrd == 1 || iOrd == 2)
1468 event[j].offsetHistory( 0, 0, 0, nOffset);
1469 else if (iOrd == 3 || iOrd == 4)
1470 event[j].offsetHistory( 0, nOffset, 0, nOffset);
1471 else if (iOrd == 5 || iOrd == 6)
1472 event[j].offsetHistory( 0, nOffset, 0, 0);
1476 if (event[j].status() == -22) {
1477 event[j].statusPos();
1478 event[j].daughters(0, 0);
1486 partonSystemsPtr->addSys();
1487 partonSystemsPtr->setInA(0, inP + nOffset);
1488 partonSystemsPtr->setInB(0, inM + nOffset);
1489 for (
int i = inM + 1; i < nHardDone; ++i)
1490 partonSystemsPtr->addOut(0, i + nOffset);
1491 partonSystemsPtr->setSHat( 0,
1492 (event[inP + nOffset].p() + event[inM + nOffset].p()).m2Calc() );
1493 partonSystemsPtr->setPTHat( 0, scale);
1498 int inP2 = iBeginSecond;
1499 int inM2 = iBeginSecond + 1;
1503 x1 = process[inP2].pPos() / process[0].e();
1504 beamAPtr->append( inP2, process[inP2].
id(), x1);
1505 x2 = process[inM2].pNeg() / process[0].e();
1506 beamBPtr->append( inM2, process[inM2].
id(), x2);
1509 scale = process.scaleSecond();
1510 beamAPtr->xfISR( 1, process[inP2].
id(), x1, scale*scale);
1511 beamAPtr->pickValSeaComp();
1512 beamBPtr->xfISR( 1, process[inM2].
id(), x2, scale*scale);
1513 beamBPtr->pickValSeaComp();
1516 for (
int i = inP2; i < process.size(); ++ i) {
1517 int mother = process[i].mother1();
1518 if ( (mother > 2 && mother < inP2) || mother > inM2 )
break;
1519 event.append(process[i]);
1523 if (event[i].status() == -22) {
1524 event[i].statusPos();
1525 event[i].daughters(0, 0);
1533 partonSystemsPtr->addSys();
1534 partonSystemsPtr->setInA(1, inP2);
1535 partonSystemsPtr->setInB(1, inM2);
1536 for (
int i = inM2 + 1; i < nHardDone; ++i)
1537 partonSystemsPtr->addOut(1, i);
1538 partonSystemsPtr->setSHat( 1,
1539 (event[inP2].p() + event[inM2].p()).m2Calc() );
1540 partonSystemsPtr->setPTHat( 1, scale);
1547 for (
int i = 0; i < process.size(); ++ i) {
1548 if (process[i].col() > maxColTag) maxColTag = process[i].col();
1549 if (process[i].acol() > maxColTag) maxColTag = process[i].acol();
1551 event.initColTag(maxColTag);
1554 for (
int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1557 int kindJunction = process.kindJunction(iJun);
1560 if (kindJunction <= 4) {
1561 int iLegF1 = (kindJunction - 1) / 2;
1562 for (
int iLeg = iLegF1; iLeg <= 2; ++iLeg) {
1563 bool colFound =
false;
1564 for (
int i = inM + 1; i <
event.size(); ++i) {
1565 int col = (kindJunction % 2 == 1) ? event[i].col() :
event[i].acol();
1566 if (col == process.colJunction(iJun,iLeg)) colFound =
true;
1568 if (!colFound) doCopy =
false;
1572 event.appendJunction( process.getJunction(iJun));
1583 void PartonLevel::setupShowerSys(
Event& process,
Event& event) {
1587 event.append( process[0]);
1591 iPosBefShow.resize( process.size());
1592 fill( iPosBefShow.begin(), iPosBefShow.end(), 0);
1595 for (
int i = 1; i < process.size(); ++i) {
1596 if (process[i].mother1() > 0)
break;
1597 int j =
event.append(process[i]);
1601 if (event[j].status() == -22) {
1602 event[j].statusPos();
1603 event[j].daughters(0, 0);
1611 partonSystemsPtr->clear();
1612 partonSystemsPtr->addSys();
1613 for (
int i = 1; i < nHardDone; ++i) partonSystemsPtr->addOut(0, i);
1614 partonSystemsPtr->setSHat( 0, pow2(process[0].m()) );
1615 partonSystemsPtr->setPTHat( 0, 0.5 * process[0].m() );
1618 for (
int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1621 int kindJunction = process.kindJunction(iJun);
1624 if (kindJunction <= 4) {
1625 int iLegF1 = (kindJunction - 1) / 2;
1626 for (
int iLeg = iLegF1; iLeg <= 2; ++iLeg) {
1627 bool colFound =
false;
1628 for (
int i = 1; i <
event.size(); ++i) {
1629 int col = (kindJunction % 2 == 1) ? event[i].col() :
event[i].acol();
1630 if (col == process.colJunction(iJun,iLeg)) colFound =
true;
1632 if (!colFound) doCopy =
false;
1636 event.appendJunction( process.getJunction(iJun));
1647 void PartonLevel::setupResolvedDiff(
Event& process) {
1650 int iDiffMot = iDS + 2;
1651 int iDiffDau = process.size() - 1;
1655 process[iDiffMot].statusNeg();
1656 process[iDiffMot].daughters( iDiffDau + 1, iDiffDau + 2);
1659 double mDiff = process[iDiffMot].m();
1660 double m2Diff = mDiff * mDiff;
1664 int idDiffA = (iDS == 1) ? process[1].
id() : 990;
1665 int idDiffB = (iDS == 2) ? process[2].
id() : 990;
1666 double mDiffA = (iDS == 1) ? process[1].m() : 0.;
1667 double mDiffB = (iDS == 2) ? process[2].m() : 0.;
1668 if (idDiffA == 22 && infoPtr->isVMDstateA()) {
1669 idDiffA = (iDS == 1) ? infoPtr->idVMDA() : 990;
1670 mDiffA = (iDS == 1) ? infoPtr->mVMDA() : 0.;
1672 if (idDiffB == 22 && infoPtr->isVMDstateB()) {
1673 idDiffB = (iDS == 2) ? infoPtr->idVMDB() : 990;
1674 mDiffB = (iDS == 2) ? infoPtr->mVMDB() : 0.;
1676 double m2DiffA = mDiffA * mDiffA;
1677 double m2DiffB = mDiffB * mDiffB;
1678 double eDiffA = 0.5 * (m2Diff + m2DiffA - m2DiffB) / mDiff;
1679 double eDiffB = 0.5 * (m2Diff + m2DiffB - m2DiffA) / mDiff;
1680 double pzDiff = 0.5 * sqrtpos( pow2(m2Diff - m2DiffA - m2DiffB)
1681 - 4. * m2DiffA * m2DiffB ) / mDiff;
1682 process.append( idDiffA, 13, iDiffMot, 0, 0, 0, 0, 0,
1683 0., 0., pzDiff, eDiffA, mDiffA);
1684 process.append( idDiffB, 13, iDiffMot, 0, 0, 0, 0, 0,
1685 0., 0., -pzDiff, eDiffB, mDiffB);
1689 beamAPtr = (iDS == 1) ? beamHadAPtr : beamPomAPtr;
1690 beamBPtr = (iDS == 2) ? beamHadBPtr : beamPomBPtr;
1691 if (infoPtr->isVMDstateA())
1692 beamAPtr = (iDS == 1) ? beamVMDAPtr : beamPomAPtr;
1693 if (infoPtr->isVMDstateB())
1694 beamBPtr = (iDS == 2) ? beamVMDBPtr : beamPomBPtr;
1697 eCMsave = infoPtr->eCM();
1698 infoPtr->setECM( mDiff);
1699 beamAPtr->newPzE( pzDiff, eDiffA);
1700 beamBPtr->newPzE( -pzDiff, eDiffB);
1703 if ( beamAPtr->id() == 990 )
1704 beamAPtr->xPom(pow2(mDiff/eCMsave));
1705 if ( beamBPtr->id() == 990 )
1706 beamBPtr->xPom(pow2(mDiff/eCMsave));
1709 int beamOffset = (sizeEvent > 0) ? sizeEvent - 1 : 4;
1710 int beamRemnOffset = iDS;
1711 if (beamAPtr->isGamma() || beamBPtr->isGamma()) beamRemnOffset = 4;
1714 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1715 timesDecPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1716 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1717 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, beamRemnOffset);
1718 colourReconnection.reassignBeamPtrs( beamAPtr, beamBPtr);
1722 if (iDS == 1) multiPtr = &multiSDA;
1723 else if (iDS == 2) multiPtr = &multiSDB;
1724 else multiPtr = &multiCD;
1732 void PartonLevel::leaveResolvedDiff(
int iHardLoop,
Event& process,
1736 Vec4 pDiffA = (iDS == 1) ? process[1].p()
1737 : process[1].p() - process[3].p();
1738 Vec4 pDiffB = (iDS == 2) ? process[2].p()
1739 : process[2].p() - process[4].p();
1741 MtoCM.fromCMframe( pDiffA, pDiffB);
1744 for (
int i = sizeProcess; i < process.size(); ++i)
1745 process[i].rotbst( MtoCM);
1746 int iFirst = (iHardLoop == 1) ? 5 + sizeEvent - sizeProcess : sizeEvent;
1747 if (isDiffC) iFirst = 6 + sizeEvent - sizeProcess;
1748 for (
int i = iFirst; i <
event.size(); ++i)
1749 event[i].rotbst( MtoCM);
1752 infoPtr->setECM( eCMsave);
1753 beamAPtr->newPzE( event[1].pz(), event[1].e());
1754 beamBPtr->newPzE( event[2].pz(), event[2].e());
1760 beamAPtr = beamHadAPtr;
1761 beamBPtr = beamHadBPtr;
1764 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1765 timesDecPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1766 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1767 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, 0);
1768 colourReconnection.reassignBeamPtrs( beamAPtr, beamBPtr);
1771 multiPtr = &multiMB;
1779 void PartonLevel::setupHardDiff(
Event& process) {
1782 Event tmpProcess = process;
1784 process.scale(tmpProcess.scale());
1787 for (
int iEntry = 0; iEntry < 3; ++iEntry)
1788 process.append( tmpProcess[iEntry]);
1791 double eCM = infoPtr->eCM();
1792 double sNow = eCM * eCM;
1793 double xPom = (isHardDiffB) ? infoPtr->xPomeronA()
1794 : infoPtr->xPomeronB();
1795 double phiPom = 2. * M_PI * rndmPtr->flat();
1796 double thetaPom = (isHardDiffB) ? hardDiffraction.getThetaPomeronA()
1797 : hardDiffraction.getThetaPomeronB();
1798 double m2Diff = xPom * sNow;
1799 double mDiff = sqrt(m2Diff);
1802 if (beamAhasGamma || beamBhasGamma) {
1803 process.append( tmpProcess[3]);
1804 process.append( tmpProcess[4]);
1811 int id1 = process[1 + gammaOffset].id();
1812 int id2 = process[2 + gammaOffset].id();
1813 int idEl = (isHardDiffB) ? id1 : id2;
1814 if (idEl == 22) idEl = 113;
1815 int idX = (isHardDiffB) ? id2 : id1;
1816 double m1 = process[1 + gammaOffset].m();
1817 double m2 = process[2 + gammaOffset].m();
1818 double mEl = (isHardDiffB) ? m1 : m2;
1823 double m3, m4, s3, s4, lambda34;
1827 mEl = particleDataPtr->mSel(113);
1828 m3 = (isHardDiffB) ? mEl : mDiff;
1829 m4 = (isHardDiffA) ? mEl : mDiff;
1832 lambda34 = sqrtpos( pow2( sNow - s3 - s4) - 4. * s3 * s4 );
1834 }
while (lambda34 <= 0. && nTry < 10);
1836 m3 = (isHardDiffB) ? mEl : mDiff;
1837 m4 = (isHardDiffA) ? mEl : mDiff;
1840 lambda34 = sqrtpos( pow2( sNow - s3 - s4) - 4. * s3 * s4 );
1842 double pAbs = 0.5 * lambda34 / eCM;
1843 Vec4 p3 = Vec4( 0., 0., pAbs, 0.5 * (sNow + s3 - s4) / eCM);
1844 Vec4 p4 = Vec4( 0., 0., -pAbs, 0.5 * (sNow + s4 - s3) / eCM);
1849 p3.rot( thetaPom, phiPom);
1850 p4.rot( thetaPom, phiPom);
1853 int status3 = (isHardDiffB) ? 14 : 15;
1854 int status4 = (isHardDiffA) ? 14 : 15;
1855 int sign = (idX > 0) ? 1 : -1;
1856 int idDiff = (idX == 22) ? 9900020 : sign * 9902210;
1857 int id3 = (isHardDiffB) ? idEl : idDiff;
1858 int id4 = (isHardDiffA) ? idEl : idDiff;
1859 process.append( id3, status3, 1 + gammaOffset, 0, 0, 0, 0, 0, p3, m3);
1860 process.append( id4, status4, 2 + gammaOffset, 0, 0, 0, 0, 0, p4, m4);
1863 process[1 + gammaOffset].daughters(3 + gammaOffset, 0);
1864 process[2 + gammaOffset].daughters(4 + gammaOffset, 0);
1865 int iDiffMot = (isHardDiffB) ? 4 + gammaOffset: 3 + gammaOffset;
1866 int iDiffRad = process.size() - 1;
1867 process[iDiffMot].statusNeg();
1868 process[iDiffMot].daughters( iDiffRad + 1, iDiffRad + 2);
1872 int idDiffA = (isHardDiffB) ? 990 : process[1 + gammaOffset].
id();
1873 int idDiffB = (isHardDiffA) ? 990 : process[2 + gammaOffset].
id();
1874 double mDiffA = (isHardDiffB) ? 0. : process[1 + gammaOffset].m();
1875 double mDiffB = (isHardDiffA) ? 0. : process[2 + gammaOffset].m();
1876 double m2DiffA = mDiffA * mDiffA;
1877 double m2DiffB = mDiffB * mDiffB;
1878 double eDiffA = 0.5 * (m2Diff + m2DiffA - m2DiffB) / mDiff;
1879 double eDiffB = 0.5 * (m2Diff + m2DiffB - m2DiffA) / mDiff;
1880 double pzDiff = 0.5 * sqrtpos( pow2(m2Diff - m2DiffA - m2DiffB)
1881 - 4. * m2DiffA * m2DiffB ) / mDiff;
1882 process.append( idDiffA, 13, iDiffMot, 0, 0, 0, 0, 0,
1883 0., 0., pzDiff, eDiffA, mDiffA);
1884 process.append( idDiffB, 13, iDiffMot, 0, 0, 0, 0, 0,
1885 0., 0., -pzDiff, eDiffB, mDiffB);
1888 vector<int> hardParton;
1889 for (
int iHard = 3 + gammaOffset; iHard < tmpProcess.size(); ++iHard)
1890 hardParton.push_back( process.append(tmpProcess[iHard]) );
1893 Vec4 pDiffA = (isHardDiffA) ? process[1 + gammaOffset].p()
1894 : process[1 + gammaOffset].p() - pD3;
1895 Vec4 pDiffB = (isHardDiffB) ? process[2 + gammaOffset].p()
1896 : process[2 + gammaOffset].p() - pD4;
1898 MtoCM.toCMframe( pDiffA, pDiffB);
1899 for (
unsigned int i = 0; i < hardParton.size(); ++i)
1900 process[hardParton[i]].rotbst(MtoCM);
1903 for (
unsigned int j = 0; j < hardParton.size(); ++j) {
1904 int mother1 = (tmpProcess[j + 3 + gammaOffset].mother1() == 0)
1905 ? 0 : tmpProcess[j + 3 + gammaOffset].mother1() + 4;
1906 int mother2 = (tmpProcess[j + 3 + gammaOffset].mother2() == 0)
1907 ? 0 : tmpProcess[j + 3 + gammaOffset].mother2() + 4;
1908 int daughter1 = (tmpProcess[j + 3 + gammaOffset].daughter1() == 0)
1909 ? 0 : tmpProcess[j + 3 + gammaOffset].daughter1() + 4;
1910 int daughter2 = (tmpProcess[j + 3 + gammaOffset].daughter2() == 0)
1911 ? 0 : tmpProcess[j + 3 + gammaOffset].daughter2() + 4;
1912 process[hardParton[j]].mothers( mother1,mother2);
1913 process[hardParton[j]].daughters( daughter1, daughter2);
1919 for (
int i = 0; i < process.size(); ++i) {
1920 if (process[i].
id() == 990 && process[i].status() == 13) iPomeron = i;
1921 if (process[i].idAbs() == idX && process[i].status() == 13) iPrtcl = i;
1925 process[iPomeron].daughters(hardParton[0], 0);
1926 process[iPrtcl].daughters(hardParton[1],0);
1927 process[hardParton[0]].mothers(iPomeron,0);
1928 process[hardParton[1]].mothers(iPrtcl, 0);
1930 process[iPomeron].daughters(hardParton[1], 0);
1931 process[iPrtcl].daughters(hardParton[0],0);
1932 process[hardParton[1]].mothers(iPomeron,0);
1933 process[hardParton[0]].mothers(iPrtcl, 0);
1937 process[iPomeron].statusNeg();
1938 process[iPrtcl].statusNeg();
1941 infoPtr->setHasUnresolvedBeams(
true);
1945 beamAPtr = (isHardDiffB) ? beamPomAPtr
1946 : (beamAhasGamma ? beamGamAPtr : beamHadAPtr);
1947 beamBPtr = (isHardDiffA) ? beamPomBPtr
1948 : (beamBhasGamma ? beamGamBPtr : beamHadBPtr);
1951 eCMsave = infoPtr->eCM();
1952 infoPtr->setECM( mDiff);
1953 beamAPtr->newPzE( pzDiff, eDiffA);
1954 beamBPtr->newPzE( -pzDiff, eDiffB);
1957 int beamOffset = 4 + gammaOffset;
1958 int beamRemnOffset = (isHardDiffB) ? 2 : 1;
1959 if (beamAPtr->isGamma() || beamBPtr->isGamma())
1960 beamRemnOffset = 4 + gammaOffset;
1963 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1964 timesDecPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1965 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
1966 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, beamRemnOffset);
1967 colourReconnection.reassignBeamPtrs( beamAPtr, beamBPtr);
1970 if (isHardDiffA) multiPtr = &multiSDA;
1971 else if (isHardDiffB) multiPtr = &multiSDB;
1974 multiPtr->setBeamOffset(beamOffset);
1980 infoPtr->setHasPomPsystem(
true);
1988 void PartonLevel::leaveHardDiff(
Event& process,
Event& event) {
1991 Vec4 pDiffA = (isHardDiffA) ? process[1 + gammaOffset].p()
1992 : process[1 + gammaOffset].p() - process[3 + gammaOffset].p();
1993 Vec4 pDiffB = (isHardDiffB) ? process[2 + gammaOffset].p()
1994 : process[2 + gammaOffset].p() - process[4 + gammaOffset].p();
1996 MtoCM.fromCMframe( pDiffA, pDiffB);
1999 for (
int i = 5 + gammaOffset; i < process.size(); ++i)
2000 process[i].rotbst( MtoCM);
2001 for (
int i = 5 + gammaOffset; i <
event.size(); ++i)
2002 event[i].rotbst( MtoCM);
2005 isHardDiffA = isHardDiffB = isHardDiff =
false;
2008 infoPtr->setECM( eCMsave);
2009 beamAPtr->newPzE( event[1 + gammaOffset].pz(), event[1 + gammaOffset].e());
2010 beamBPtr->newPzE( event[2 + gammaOffset].pz(), event[2 + gammaOffset].e());
2014 beamAPtr = beamAhasGamma ? beamGamAPtr : beamHadAPtr;
2015 beamBPtr = beamBhasGamma ? beamGamBPtr : beamHadBPtr;
2018 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2019 timesDecPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2020 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2021 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2022 colourReconnection.reassignBeamPtrs( beamAPtr, beamBPtr);
2025 multiPtr->setBeamOffset(0);
2028 multiPtr = &multiMB;
2037 bool PartonLevel::setupResolvedLeptonGamma(
Event& process) {
2040 eCMsaveGamma = infoPtr->eCM();
2046 gammaOffset = beamOffset;
2049 double mGmGm = (infoPtr->nFinal() > 1) || (gammaModeEvent != 4)
2050 ? infoPtr->eCMsub() : sqrt( infoPtr->sHat());
2051 double m2GmGm = pow2(mGmGm);
2054 double m2Gamma1 = hasGammaA ? 0. : pow2(beamAPtr->m());
2055 double m2Gamma2 = hasGammaB ? 0. : pow2(beamBPtr->m());
2058 double eGamA = 0.5 * (m2GmGm + m2Gamma1 - m2Gamma2) / mGmGm;
2059 double eGamB = 0.5 * (m2GmGm + m2Gamma2 - m2Gamma1) / mGmGm;
2060 double pzGam = 0.5 * sqrtpos( pow2(m2GmGm - m2Gamma1 - m2Gamma2)
2061 - 4. * m2Gamma1 * m2Gamma2 ) / mGmGm;
2062 Vec4 pGammaANew( 0, 0, pzGam, eGamA);
2063 Vec4 pGammaBNew( 0, 0, -pzGam, eGamB);
2066 beamGamAPtr->newPzE( pzGam, eGamA);
2067 beamGamBPtr->newPzE( -pzGam, eGamB);
2070 Vec4 pGammaA = process[iBeamA].p();
2071 Vec4 pGammaB = process[iBeamB].p();
2074 RotBstMatrix MtoGammaGamma;
2075 MtoGammaGamma.toCMframe( pGammaA, pGammaB);
2076 process.rotbst(MtoGammaGamma);
2079 process[iBeamA].p(pGammaANew);
2080 process[iBeamB].p(pGammaBNew);
2083 if ( !hasGammaA && !(beamBPtr->getGammaMode() == 2) )
2084 process[iBeamA].m( sqrt(m2Gamma1) );
2085 if ( !hasGammaB && !(beamAPtr->getGammaMode() == 2) )
2086 process[iBeamB].m( sqrt(m2Gamma2) );
2089 if ( gammaModeEvent == 4)
return true;
2092 if ( hasGammaA) beamAPtr = beamGamAPtr;
2093 else beamAPtr->newPzE( pzGam, eGamA);
2094 if ( hasGammaB) beamBPtr = beamGamBPtr;
2095 else beamBPtr->newPzE( -pzGam, eGamB);
2098 if ( (beamAhasResGamma && (!beamBhasResGamma && hasGammaB))
2099 || ( (!beamAhasResGamma && hasGammaA) && beamBhasResGamma) )
2100 infoPtr->setHasUnresolvedBeams(
true);
2103 infoPtr->setECM( mGmGm);
2106 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
2107 timesDecPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
2108 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
2109 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
2110 colourReconnection.reassignBeamPtrs( beamAPtr, beamBPtr);
2113 multiPtr = &multiGmGm;
2114 multiPtr->setBeamOffset(beamOffset);
2126 void PartonLevel::leaveResolvedLeptonGamma(
Event& process,
Event& event,
2130 Vec4 pLeptonA = process[1].p();
2131 Vec4 pLeptonB = process[2].p();
2138 int iSkipHardDiff = -1;
2141 RotBstMatrix MtoLeptonLepton;
2142 MtoLeptonLepton.toCMframe( pLeptonA, pLeptonB);
2145 process.rotbst(MtoLeptonLepton);
2146 event.rotbst(MtoLeptonLepton);
2149 infoPtr->setECM( eCMsaveGamma);
2152 if ( hasGammaA) beamAPtr = beamHadAPtr;
2153 if ( hasGammaB) beamBPtr = beamHadBPtr;
2156 double m2BeamA = pow2( beamAPtr->m());
2157 double m2BeamB = pow2( beamBPtr->m());
2160 double sCM = infoPtr->s();
2161 double eCM2A = 0.25 * pow2(sCM + m2BeamA - m2BeamB) / sCM;
2162 double eCM2B = 0.25 * pow2(sCM - m2BeamA + m2BeamB) / sCM;
2165 Vec4 pGamma1Orig = process[3].p();
2166 Vec4 pGamma2Orig = process[4].p();
2167 Vec4 pGamma1 = pGamma1Orig;
2168 Vec4 pGamma2 = pGamma2Orig;
2169 double mGamma1 = sqrt(m2BeamA);
2170 double mGamma2 = sqrt(m2BeamB);
2174 if ( (infoPtr->nFinal() == 1) && (gammaModeEvent == 4)
2175 && hasGammaA && hasGammaB ) {
2178 mGamma1 = -sqrt( beamAPtr->Q2Gamma());
2179 mGamma2 = -sqrt( beamBPtr->Q2Gamma());
2180 double kTxA = beamAPtr->gammaKTx();
2181 double kTyA = beamAPtr->gammaKTy();
2182 double kTxB = beamBPtr->gammaKTx();
2183 double kTyB = beamBPtr->gammaKTy();
2186 double mT2A = pow2( beamAPtr->gammaKT() ) - pow2(mGamma1);
2187 double mT2B = pow2( beamBPtr->gammaKT() ) - pow2(mGamma2);
2188 double sHatBef = infoPtr->sHat();
2189 double sHatAft = infoPtr->sHat() + pow2(kTxA + kTxB) + pow2(kTyA + kTyB);
2190 double lambda = pow2(sHatAft) + pow2(mT2A) + pow2(mT2B)
2191 - 2. * ( sHatAft * (mT2A + mT2B) + mT2A * mT2B );
2192 double kz2New = 0.25 * lambda / sHatAft;
2195 Vec4 pGamma1New(kTxA, kTyA, sqrt(kz2New), sqrt(kz2New + mT2A) );
2196 Vec4 pGamma2New(kTxB, kTyB, -sqrt(kz2New), sqrt(kz2New + mT2B) );
2199 RotBstMatrix MfromGmGmOrig;
2200 MfromGmGmOrig.toCMframe( pGamma1Orig, pGamma2Orig);
2201 MfromGmGmOrig.invert();
2202 pGamma1New.rotbst(MfromGmGmOrig);
2203 pGamma2New.rotbst(MfromGmGmOrig);
2206 pGamma1 = pGamma1New;
2207 pGamma2 = pGamma2New;
2210 event[3].p( pGamma1);
2211 event[3].m( mGamma1);
2212 event[4].p( pGamma2);
2213 event[4].m( mGamma2);
2216 double rescale = sqrt(sHatAft / sHatBef);
2217 double wPosBef = beamAPtr->xGamma() * infoPtr->eCM();
2218 double wNegBef = beamBPtr->xGamma() * infoPtr->eCM();
2221 double w2BeamA = pow2(beamAPtr->gammaKT()) + m2BeamA;
2222 double w2BeamB = pow2(beamBPtr->gammaKT()) + m2BeamB;
2223 double wPosRem = infoPtr->eCM() - rescale * wPosBef;
2224 double wNegRem = infoPtr->eCM() - rescale * wNegBef;
2225 double w2Rem = wPosRem * wNegRem;
2226 double xLepA = 1. - beamAPtr->xGamma();
2227 double xLepB = 1. - beamBPtr->xGamma();
2230 double lambdaRoot = sqrtpos( pow2(w2Rem - w2BeamA - w2BeamB)
2231 - 4. * w2BeamA * w2BeamB );
2232 double rescaleA = (w2Rem + w2BeamA - w2BeamB + lambdaRoot)
2233 / (2. * w2Rem * xLepA);
2234 double rescaleB = (w2Rem + w2BeamB - w2BeamA + lambdaRoot)
2235 / (2. * w2Rem * xLepB);
2238 double pPosA = rescaleA * xLepA * wPosRem;
2239 double pNegA = w2BeamA / pPosA;
2240 double eLepA = 0.5 * (pPosA + pNegA);
2241 double pzLepA = 0.5 * (pPosA - pNegA);
2242 double pNegB = rescaleB * xLepB * wNegRem;
2243 double pPosB = w2BeamB / pNegB;
2244 double eLepB = 0.5 * (pPosB + pNegB);
2245 double pzLepB = 0.5 * (pPosB - pNegB);
2246 Vec4 pLeptonANew( -kTxA, -kTyA, pzLepA, eLepA);
2247 Vec4 pLeptonBNew( -kTxB, -kTyB, pzLepB, eLepB);
2250 pLepton1scat = pLeptonANew;
2251 pLepton2scat = pLeptonBNew;
2254 double theta1 = pLepton1scat.theta();
2255 double theta2 = M_PI - pLepton2scat.theta();
2256 infoPtr->setTheta1(theta1);
2257 infoPtr->setTheta2(theta2);
2258 infoPtr->setECMsub(sqrt(sHatBef));
2259 infoPtr->setsHatNew(sHatBef);
2268 double xGamma1 = beamAPtr->xGamma();
2269 double Q2gamma1 = beamAPtr->Q2Gamma();
2270 mGamma1 = -sqrt( Q2gamma1);
2271 beamGamAPtr->newM( mGamma1);
2274 double eGamma1 = xGamma1 * sqrt( eCM2A);
2275 double kz1 = (eCM2A * xGamma1 + 0.5 * Q2gamma1)
2276 / sqrt(eCM2A - m2BeamA);
2279 pGamma1 = Vec4( beamAPtr->gammaKTx(), beamAPtr->gammaKTy(),
2281 event[3].p( pGamma1);
2282 event[3].m( mGamma1);
2285 if (infoPtr->isHardDiffractiveA()) {
2286 event[7].p( pGamma1);
2287 event[7].m( mGamma1);
2296 double xGamma2 = beamBPtr->xGamma();
2297 double Q2gamma2 = beamBPtr->Q2Gamma();
2298 mGamma2 = -sqrt( Q2gamma2);
2299 beamGamBPtr->newM( mGamma2);
2302 double eGamma2 = xGamma2 * sqrt( eCM2B);
2303 double kz2 = (eCM2B * xGamma2 + 0.5 * Q2gamma2)
2304 / sqrt(eCM2B - m2BeamB);
2307 pGamma2 = Vec4( beamBPtr->gammaKTx(), beamBPtr->gammaKTy(),
2309 event[4].p( pGamma2);
2310 event[4].m( mGamma2);
2313 if (infoPtr->isHardDiffractiveB()) {
2314 event[8].p( pGamma2);
2315 event[8].m( mGamma2);
2321 Vec4 pLepton1 = process[1].p();
2322 Vec4 pLepton2 = process[2].p();
2323 pLepton1scat = pLepton1 - pGamma1;
2324 pLepton2scat = pLepton2 - pGamma2;
2329 RotBstMatrix MfromGmGm;
2330 MfromGmGm.toCMframe( pGamma1, pGamma2);
2331 MfromGmGm.fromCMframe( pGamma1Orig, pGamma2Orig);
2336 int iSkipGamma = -1;
2337 if ( gammaModeEvent == 3 ) {
2339 event[iSkipGamma].m( mGamma1);
2340 event[iSkipGamma].p( pGamma1);
2341 }
else if ( gammaModeEvent == 2 ) {
2343 event[iSkipGamma].m( mGamma2);
2344 event[iSkipGamma].p( pGamma2);
2348 for (
int i = 5; i <
event.size(); ++i) {
2349 if ( (i != iSkipGamma) && (i != iSkipHardDiff) )
2350 event[i].rotbst( MfromGmGm);
2354 if ( doRemnants && physical) {
2358 int iPosLepton1 =
event.append( beamAPtr->id(), 63, 1, 0, 0, 0, 0, 0,
2359 pLepton1scat, beamAPtr->m());
2360 event[1].daughter2( event[1].daughter1());
2361 event[1].daughter1( iPosLepton1);
2364 int iPosLepton2 =
event.append( beamBPtr->id(), 63, 2, 0, 0, 0, 0, 0,
2365 pLepton2scat, beamBPtr->m());
2366 event[2].daughter2( event[2].daughter1());
2367 event[2].daughter1( iPosLepton2);
2372 if ( gammaModeEvent == 4)
return;
2375 timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2376 timesDecPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2377 spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2378 remnants.reassignBeamPtrs( beamAPtr, beamBPtr, 0);
2379 colourReconnection.reassignBeamPtrs( beamAPtr, beamBPtr);
2382 multiPtr = &multiMB;
2383 multiPtr->setBeamOffset(0);
2391 void PartonLevel::cleanEventFromGamma(
Event& event) {
2396 if (infoPtr->isHardDiffractiveA() || infoPtr->isHardDiffractiveB())
2398 int iPosBeam1 = 1 + beamOffset;
2399 int iPosBeam2 = 2 + beamOffset;
2404 for (
int i = event.size() - 1; i > 0; --i) {
2405 if ( (event[i].
id() == 22 && event[i].mother1() == iPosBeam1)
2406 && beamAhasResGamma ) iPosGamma1 = i;
2407 if ( (event[i].
id() == 22 && event[i].mother1() == iPosBeam2)
2408 && beamBhasResGamma ) iPosGamma2 = i;
2413 if (iPosGamma1 > 0) ++nGamma;
2414 if (iPosGamma2 > 0) ++nGamma;
2417 if ( nGamma == 0 )
return;
2420 for (
int i = 0; i < nGamma; ++i) {
2423 int iPosGamma = (iPosGamma1 > 0 && i == 0) ? iPosGamma1 : iPosGamma2;
2424 int iPosBeam = (iPosGamma1 > 0 && i == 0) ? iPosBeam1 : iPosBeam2;
2427 while ( iPosGamma > iPosBeam ) {
2428 int iDaughter1 =
event[iPosGamma].daughter1();
2429 int iDaughter2 =
event[iPosGamma].daughter2();
2430 int iMother1 =
event[iPosGamma].mother1();
2431 int iMother2 =
event[iPosGamma].mother2();
2434 if ( iDaughter1 == iDaughter2 ) {
2435 event[iDaughter1].mothers( iMother1, iMother2 );
2436 event.remove( iPosGamma, iPosGamma);
2437 iPosGamma = iDaughter1;
2441 event[iMother1].daughters( iDaughter1, iDaughter2 );
2442 event[iDaughter1].mother1( iMother1 );
2443 event[iDaughter2].mother1( iMother1 );
2444 event.remove( iPosGamma, iPosGamma);
2445 iPosGamma = iMother1;
2449 if ( (i == 0 && nGamma > 1) && iPosGamma2 > iPosGamma ) --iPosGamma2;
2459 bool PartonLevel::resonanceShowers(
Event& process,
Event& event,
2465 nHardDoneRHad = nHardDone;
2466 inRHadDecay.resize(0);
2467 for (
int i = 0; i < process.size(); ++i)
2468 inRHadDecay.push_back(
false);
2469 }
else nHardDone = nHardDoneRHad;
2476 int nBranchMax = (doTrial) ? nTrialEmissions : -1;
2478 vector<int> iJunCopied;
2480 while (nHardDone < process.size()) {
2482 int iBegin = nHardDone;
2488 bool comesFromR =
false;
2489 int iTraceUp = process[iBegin].mother1();
2491 if ( rHadronsPtr->givesRHadron(process[iTraceUp].id()) )
2493 iTraceUp = process[iTraceUp].mother1();
2494 }
while (iTraceUp > 0 && !comesFromR);
2496 inRHadDecay[iBegin] =
true;
2503 }
else if (!inRHadDecay[iBegin]) {
2510 int iHardMother = process[iBegin].mother1();
2511 Particle& hardMother = process[iHardMother];
2512 int iBefMother = iPosBefShow[iHardMother];
2513 int iAftMother =
event[iBefMother].iBotCopyId();
2516 int iTraceRHadron = rHadronsPtr->trace( iAftMother);
2517 if (iTraceRHadron > 0) iAftMother = iTraceRHadron;
2519 Particle& aftMother =
event[iAftMother];
2522 aftMother.statusNeg();
2526 int colBef = hardMother.col();
2527 int acolBef = hardMother.acol();
2528 int colAft = aftMother.col();
2529 int acolAft = aftMother.acol();
2531 M.bst( hardMother.p(), aftMother.p());
2535 if ( (colBef != 0 || acolBef != 0) && doReconnect && reconnectMode == 1
2536 && forceResonanceCR && !earlyResDec) {
2537 infoPtr->errorMsg(
"Abort from PartonLevel::resonanceShower: "
2538 "new CR can't handle separate CR for coloured resonance decays");
2539 infoPtr->setAbortPartonLevel(
true);
2544 vector<bool> doCopyJun;
2545 for (
int i = iBegin; i < process.size(); ++i) {
2546 if (process[i].mother1() != iHardMother)
break;
2547 int iNow =
event.append( process[i] );
2548 iPosBefShow[i] = iNow;
2549 Particle& now =
event.back();
2552 if (now.status() == -22) {
2554 now.daughters(0, 0);
2558 if (i == iBegin)
event[iAftMother].daughter1( iNow);
2559 else event[iAftMother].daughter2( iNow);
2560 now.mother1(iAftMother);
2563 for (
int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
2564 if (iJun >=
int(doCopyJun.size())) doCopyJun.push_back(
false);
2566 if (doCopyJun[iJun])
continue;
2568 int kindJunction = process.kindJunction(iJun);
2569 if (kindJunction <= 2) {
2573 int iDau1 = hardMother.daughter1();
2574 int iDau2 = hardMother.daughter2();
2576 if (iDau1 == 0 || iDau2 - iDau1 < 2)
continue;
2577 for (
int iDau=iDau1; iDau<=iDau2; ++iDau) {
2578 int colDau = (kindJunction == 1 ? process[iDau].col()
2579 : process[iDau].acol());
2580 for (
int iLeg = 0; iLeg <= 2; ++iLeg)
2581 if ( process.colJunction(iJun,iLeg) == colDau ) nMatch += 1;
2584 if (nMatch == 3) doCopyJun[iJun] =
true;
2585 }
else if (kindJunction <= 4) {
2588 int col = (kindJunction == 3 ? hardMother.acol() : hardMother.col());
2589 if ( process.colJunction(iJun,0) == col ) doCopyJun[iJun] =
true;
2593 for (
int kJun = 0; kJun <
event.sizeJunction(); ++kJun) {
2595 for (
int iLeg = 0; iLeg <= 2; ++iLeg)
2596 if (event.colJunction(kJun,iLeg) == process.colJunction(iJun,iLeg))
2598 if (nMatch == 3) doCopyJun[iJun] =
false;
2603 if (now.col() == colBef) now.col( colAft);
2604 if (now.acol() == acolBef) now.acol( acolAft);
2606 if (now.col() == -acolBef) now.col( -acolAft);
2607 if (now.acol() == -colBef) now.acol( -colAft);
2611 if (now.hasVertex()) now.vProd( event[iAftMother].vDec() );
2616 int iEnd = nHardDone - 1;
2619 for (
int iJun = 0; iJun < int(doCopyJun.size()); ++iJun) {
2621 for (
int jJun = 0; jJun < int(iJunCopied.size()); ++jJun)
2622 if (iJunCopied[jJun] == iJun) doCopyJun[iJun] =
false;
2624 if (!doCopyJun[iJun])
continue;
2626 Junction junCopy = process.getJunction(iJun);
2627 for (
int iLeg = 0; iLeg <= 2; ++iLeg) {
2628 int colLeg = junCopy.col(iLeg);
2629 if (colLeg == colBef) junCopy.col(iLeg, colAft);
2630 if (colLeg == acolBef) junCopy.col(iLeg, acolAft);
2632 event.appendJunction(junCopy);
2634 iJunCopied.push_back(iJun);
2641 int iSys = partonSystemsPtr->addSys();
2642 partonSystemsPtr->setSHat(iSys, pow2(hardMother.m()) );
2643 partonSystemsPtr->setPTHat(iSys, 0.5 * hardMother.m() );
2646 for (
int i = iPosBefShow[iBegin]; i <= iPosBefShow[iEnd]; ++i)
2647 if (event[i].isFinal()) partonSystemsPtr->addOut( iSys, i);
2650 if (doFSRinResonances) {
2651 double pTmax = 0.5 * hardMother.m();
2652 if (canSetScale) pTmax
2653 = userHooksPtr->scaleResonance( iAftMother, event);
2657 if (doTrial) pTmax = process.scale();
2661 int iMother1 = hardMother.mother1();
2662 int iMother2 = hardMother.mother2();
2664 && event[iMother1].colType() == 0 && event[iMother2].colType() == 0)
2665 pTmax = process.scale();
2668 timesDecPtr->prepare( iSys, event);
2675 double pTveto = pTvetoPT;
2680 double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
2683 if (pTveto > 0. && pTveto > pTtimes) {
2685 doVeto = userHooksPtr->doVetoPT( 5, event);
2687 if (doVeto)
return false;
2692 if (timesDecPtr->branch( event)) {
2695 if (canVetoStep && nFSRhard <= nVetoStep) typeVetoStep = 5;
2698 pTLastBranch = pTtimes;
2709 if (typeVetoStep > 0) {
2710 doVeto = userHooksPtr->doVetoStep( typeVetoStep, 0, nFSRhard,
2713 if (doVeto)
return false;
2717 if ( canRemoveEvent && nFSRhard == 1 ) {
2719 mergingHooksPtr->doVetoStep( process, event,
true );
2723 }
while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax) );
2729 if (skipForR) nFSRinRes = nFSRres;
2738 bool PartonLevel::wzDecayShowers(
Event& event) {
2741 for (
int iWZ = 0; iWZ <
event.size(); ++iWZ)
2742 if (event[iWZ].isFinal()
2743 && (
event[iWZ].id() == 23 ||
event[iWZ].idAbs() == 24) ) {
2746 if ( event[iWZ].canDecay() &&
event[iWZ].mayDecay()
2747 &&
event[iWZ].isResonance() ) ;
2750 int iWZtop =
event[iWZ].iTopCopy();
2752 if (event[iWZtop].statusAbs() == 56) typeWZ = 1;
2753 if (event[iWZtop].statusAbs() == 47) typeWZ = 2;
2754 int sizeSave =
event.size();
2759 int idSave =
event[iWZ].id();
2760 event[iWZ].id( (idSave > 0) ? idSave + 70 : idSave - 70);
2761 int statusSave =
event[iWZ].status();
2762 resonanceDecays.next( event, iWZ);
2763 event[iWZ].id( idSave);
2764 if (event.size() - sizeSave != 2)
return true;
2765 event[iWZ].status( -statusSave);
2772 vector<int> iSisters =
event[iWZtop].sisterList();
2773 if (iSisters.size() != 1) {
2774 infoPtr->errorMsg(
"Error in PartonLevel::wzDecayShowers: "
2775 "Not able to find a single sister particle");
2778 int iEmitter = iSisters[0];
2781 RotBstMatrix MtoNew, MtoRest, MtoCM;
2782 MtoNew.bst( event[iWZtop].p(), event[iWZ].p());
2783 MtoRest.bstback( event[iWZ].p());
2784 MtoCM.bst( event[iWZ].p());
2787 Vec4 pEmitter =
event[iEmitter].p();
2788 pEmitter.rotbst( MtoNew);
2789 pEmitter.rotbst( MtoRest);
2790 if (event[iWZtop + 1].statusAbs() != 52) {
2791 infoPtr->errorMsg(
"Error in PartonLevel::wzDecayShowers: "
2792 "Found wrong recoiler");
2795 Vec4 pRecoiler =
event[iWZtop + 1].p();
2796 pRecoiler.rotbst( MtoNew);
2797 pRecoiler.rotbst( MtoRest);
2798 Vec4 pWZRest =
event[iWZ].p();
2799 pWZRest.rotbst( MtoRest);
2803 Vec4 p5 = pRecoiler;
2804 if (event[iEmitter].
id() < 0) swap( p4, p5);
2808 Vec4 pDec1 =
event[sizeSave].p();
2809 Vec4 pDec2 =
event[sizeSave + 1].p();
2810 if (event[sizeSave].
id() < 0) swap( pDec1, pDec2);
2811 pDec1.rotbst( MtoRest);
2812 pDec2.rotbst( MtoRest);
2815 double li2, ri2, lf2, rf2;
2817 if (event[iWZ].
id() == 23) {
2818 li2 = pow2(couplingsPtr->lf( event[iEmitter].idAbs() ));
2819 ri2 = pow2(couplingsPtr->rf( event[iEmitter].idAbs() ));
2820 lf2 = pow2(couplingsPtr->lf( event[sizeSave].idAbs() ));
2821 rf2 = pow2(couplingsPtr->rf( event[sizeSave].idAbs() ));
2822 if ( abs( event[iEmitter].pol() + 1.) < 0.1) ri2 = 0.;
2823 if ( abs( event[iEmitter].pol() - 1.) < 0.1) li2 = 0.;
2833 double sWZER = (p4 + pWZRest + p5).m2Calc();
2834 double x1 = 2. * p4 * (p4 + pWZRest + p5) / sWZER;
2835 double x2 = 2. * p5 * (p4 + pWZRest + p5) / sWZER;
2836 double x1s = x1 * x1;
2837 double x2s = x2 * x2;
2838 double m2Sel = pWZRest.m2Calc();
2839 double rWZER = m2Sel / sWZER;
2843 con[0] = 2. * m2Sel * (1.-x1) * ((x2s+1.-x1-x2) - rWZER * (1.-x2))
2844 * (li2 * lf2 + ri2 * rf2);
2845 con[1] = 2. * m2Sel * (1.-x2) * ((x1s+1.-x1-x2) - rWZER * (1.-x1))
2846 * (li2 * rf2 + ri2 * lf2);
2847 con[2] = 2. * m2Sel * (1.-x1) * ((x2s+1.-x1-x2) - rWZER * (1.-x2))
2848 * (li2 * rf2 + ri2 * lf2);
2849 con[3] = 2. * m2Sel * (1.-x2) * ((x1s+1.-x1-x2) - rWZER * (1.-x1))
2850 * (li2 * lf2 + ri2 * rf2);
2851 con[4] = m2Sel * sWZER * (1.-x1) * (1.-x2) * ((x1+x2-1.) + rWZER)
2852 * (li2 + ri2) * (lf2 + rf2);
2853 con[5] = -4. * (1.-x1) * (1.-x2) * (li2 + ri2) * (lf2 + rf2);
2854 con[6] = -4. * (1.-x1) * (1.-x2) * (li2 + ri2) * (lf2 + rf2);
2855 con[7] = 4. * (1.-x1) * ((1.-x1) - rWZER * (1.-x2))
2856 * (li2 + ri2) * (lf2 + rf2);
2857 con[8] = 4. * (1.-x2) * ((1.-x2) - rWZER * (1.-x1))
2858 * (li2 + ri2) * (lf2 + rf2);
2862 double pAbs12 = pDec1.pAbs();
2863 for (
int j = 0; j < 6; ++j) {
2864 Vec4 pDec1Test( 0., 0., 0., pDec1.e());
2865 Vec4 pDec2Test( 0., 0., 0., pDec2.e());
2866 if (j == 0) { pDec1Test.px( pAbs12); pDec1Test.px( -pAbs12);}
2867 else if (j == 1) { pDec1Test.px( -pAbs12); pDec1Test.px( pAbs12);}
2868 else if (j == 2) { pDec1Test.py( pAbs12); pDec1Test.py( -pAbs12);}
2869 else if (j == 3) { pDec1Test.py( -pAbs12); pDec1Test.py( pAbs12);}
2870 else if (j == 4) { pDec1Test.pz( pAbs12); pDec1Test.pz( -pAbs12);}
2871 else if (j == 5) { pDec1Test.pz( -pAbs12); pDec1Test.pz( pAbs12);}
2874 double p2p4Test = p4 * pDec1Test;
2875 double p3p4Test = p4 * pDec2Test;
2876 double p2p5Test = p5 * pDec1Test;
2877 double p3p5Test = p5 * pDec2Test;
2878 double testValues[9] = { p2p4Test, p2p5Test, p3p4Test, p3p5Test, 1.,
2879 p2p5Test * p3p4Test, p2p4Test * p3p5Test, p2p4Test * p3p4Test,
2880 p2p5Test * p3p5Test};
2882 for (
int i = 0; i < 9; ++i) wtTest += con[i] * testValues[i];
2883 if (wtTest > wtMax) wtMax = wtTest;
2895 RotBstMatrix MrndmRot;
2896 MrndmRot.rot( acos(2. * rndmPtr->flat() - 1.),
2897 2. * M_PI * rndmPtr->flat());
2898 pDec1.rotbst(MrndmRot);
2899 pDec2.rotbst(MrndmRot);
2905 double p2p4 = p4 * pDec1;
2906 double p3p4 = p4 * pDec2;
2907 double p2p5 = p5 * pDec1;
2908 double p3p5 = p5 * pDec2;
2911 double wtValues[9] = { p2p4, p2p5, p3p4, p3p5, 1., p2p5 * p3p4,
2912 p2p4 * p3p5, p2p4 * p3p4, p2p5 * p3p5};
2914 for (
int i = 0; i < 9; ++i) wt += con[i] * wtValues[i];
2915 if (wt > wtMax || wt < 0.) {
2916 infoPtr->errorMsg(
"Error in PartonLevel::wzDecayShowers: "
2917 "wt bigger than wtMax or less than zero");
2920 }
while (wt < wtMax * rndmPtr->flat());
2924 pDec1.rotbst( MtoCM);
2925 pDec2.rotbst( MtoCM);
2926 if (event[sizeSave].
id() > 0) {
2927 event[sizeSave].p( pDec1);
2928 event[sizeSave + 1].p( pDec2);
2931 event[sizeSave].p( pDec2);
2932 event[sizeSave + 1].p( pDec1);
2941 double iMother =
event[iWZtop].mother1();
2944 RotBstMatrix MtoNew, MtoRest, MtoCM;
2945 MtoNew.bst( event[iWZtop].p(), event[iWZ].p());
2946 MtoRest.bstback( event[iWZ].p());
2947 MtoCM.bst( event[iWZ].p());
2951 if (event[iMother + 1].statusAbs() == 42) iRecoiler = iMother + 1;
2952 else if (event[iMother - 1].statusAbs() == 42) iRecoiler = iMother - 1;
2954 infoPtr->errorMsg(
"Error in PartonLevel::wzDecayShowers: "
2955 "Found wrong recoiler");
2960 Vec4 pMother =
event[iMother].p();
2961 pMother.rotbst( MtoNew);
2962 pMother.rotbst( MtoRest);
2963 Vec4 pRecoiler =
event[iRecoiler].p();
2964 pRecoiler.rotbst( MtoNew);
2965 pRecoiler.rotbst( MtoRest);
2966 Vec4 pWZRest =
event[iWZ].p();
2967 pWZRest.rotbst( MtoRest);
2972 Vec4 p2 = pRecoiler;
2973 if (event[iMother].
id() < 0) swap( p1, p2);
2977 Vec4 pDec1 =
event[sizeSave].p();
2978 Vec4 pDec2 =
event[sizeSave + 1].p();
2979 if (event[sizeSave].
id() < 0) swap( pDec1, pDec2);
2980 pDec1.rotbst( MtoRest);
2981 pDec2.rotbst( MtoRest);
2984 double li2, ri2, lf2, rf2;
2986 if (event[iWZ].
id() == 23) {
2987 li2 = pow2(couplingsPtr->lf( event[iMother].idAbs() ));
2988 ri2 = pow2(couplingsPtr->rf( event[iMother].idAbs() ));
2989 lf2 = pow2(couplingsPtr->lf( event[sizeSave].idAbs() ));
2990 rf2 = pow2(couplingsPtr->rf( event[sizeSave].idAbs() ));
2991 if ( abs( event[iMother].pol() + 1.) < 0.1) ri2 = 0.;
2992 if ( abs( event[iMother].pol() - 1.) < 0.1) li2 = 0.;
3002 double s = (p1 + p2).m2Calc();
3003 double t = (p1-pWZRest).m2Calc();
3004 double u = (p2-pWZRest).m2Calc();
3005 double m2Sel = pWZRest.m2Calc();
3006 double m2Final = t + u + s - m2Sel;
3010 con[0] = 2. * m2Sel * (u * (1. - m2Final / t) + s)
3011 * (li2 * rf2 + ri2 * lf2);
3012 con[1] = 2. * m2Sel * (u * (1. - m2Final / t) + s)
3013 * (li2 * lf2 + ri2 * rf2);
3014 con[2] = 2. * m2Sel * (t * (1. - m2Final / u) + s)
3015 * (li2 * lf2 + ri2 * rf2);
3016 con[3] = 2. * m2Sel * (t * (1. - m2Final / u) + s)
3017 * (li2 * rf2 + ri2 * lf2);
3018 con[4] = m2Sel * m2Final * s * (li2 + ri2) * (lf2 + rf2);
3019 con[5] = -4. * m2Final * (li2 + ri2) * (lf2 + rf2);
3020 con[6] = -4. * m2Final * (li2 + ri2) * (lf2 + rf2);
3021 con[7] = 4. * (m2Final * u / t - m2Sel) * (li2 + ri2) * (lf2 + rf2);
3022 con[8] = 4. * (m2Final * t / u - m2Sel) * (li2 + ri2) * (lf2 + rf2);
3026 double pAbs12 = pDec1.pAbs();
3027 for (
int j = 0; j < 6; ++j) {
3028 Vec4 pDec1Test( 0., 0., 0., pDec1.e());
3029 Vec4 pDec2Test( 0., 0., 0., pDec2.e());
3030 if (j == 0) { pDec1Test.px( pAbs12); pDec1Test.px( -pAbs12);}
3031 else if (j == 1) { pDec1Test.px( -pAbs12); pDec1Test.px( pAbs12);}
3032 else if (j == 2) { pDec1Test.py( pAbs12); pDec1Test.py( -pAbs12);}
3033 else if (j == 3) { pDec1Test.py( -pAbs12); pDec1Test.py( pAbs12);}
3034 else if (j == 4) { pDec1Test.pz( pAbs12); pDec1Test.pz( -pAbs12);}
3035 else if (j == 5) { pDec1Test.pz( -pAbs12); pDec1Test.pz( pAbs12);}
3038 double p1p4Test = p1 * pDec1Test;
3039 double p1p5Test = p1 * pDec2Test;
3040 double p2p4Test = p2 * pDec1Test;
3041 double p2p5Test = p2 * pDec2Test;
3042 double testValues[9] = { p1p4Test, p1p5Test, p2p4Test, p2p5Test, 1.,
3043 p1p4Test * p2p5Test, p1p5Test * p2p4Test, p1p4Test * p1p5Test,
3044 p2p4Test * p2p5Test};
3046 for (
int i = 0; i < 9; ++i) wtTest += con[i] * testValues[i];
3047 if (wtTest > wtMax) wtMax = wtTest;
3059 RotBstMatrix MrndmRot;
3060 MrndmRot.rot( acos(2. * rndmPtr->flat() - 1.),
3061 2. * M_PI * rndmPtr->flat());
3062 pDec1.rotbst(MrndmRot);
3063 pDec2.rotbst(MrndmRot);
3069 double p1p4 = p1 * pDec1;
3070 double p1p5 = p1 * pDec2;
3071 double p2p4 = p2 * pDec1;
3072 double p2p5 = p2 * pDec2;
3075 double wtValues[9] = { p1p4, p1p5, p2p4, p2p5, 1., p1p4 * p2p5,
3076 p1p5 * p2p4, p1p4 * p1p5, p2p4 * p2p5};
3078 for (
int i = 0; i < 9; ++i) wt += con[i] * wtValues[i];
3079 if (wt > wtMax || wt < 0.) {
3080 infoPtr->errorMsg(
"Error in PartonLevel::wzDecayShowers: "
3081 "wt bigger than wtMax or less than zero");
3084 }
while (wt < wtMax * rndmPtr->flat());
3088 pDec1.rotbst( MtoCM);
3089 pDec2.rotbst( MtoCM);
3090 if (event[sizeSave].
id() > 0) {
3091 event[sizeSave].p( pDec1);
3092 event[sizeSave + 1].p( pDec2);
3095 event[sizeSave].p( pDec2);
3096 event[sizeSave + 1].p( pDec1);
3104 double pTmax = 0.5 *
event[iWZ].m();
3105 int iSys = partonSystemsPtr->addSys();
3106 partonSystemsPtr->setSHat(iSys, pow2(event[iWZ].m()) );
3107 partonSystemsPtr->setPTHat(iSys, pTmax );
3108 for (
int i = sizeSave; i <
event.size(); ++i)
3109 partonSystemsPtr->addOut( iSys, i);
3112 if (doFSRinResonances) {
3115 timesDecPtr->prepare( iSys, event);
3119 double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
3123 timesDecPtr->branch( event);
3131 }
while (pTmax > 0.);