Line data Source code
1 : // PhaseSpace.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 :
6 : // Function definitions (not found in the header) for the
7 : // PhaseSpace and PhaseSpace2to2tauyz classes.
8 :
9 : #include "Pythia8/PhaseSpace.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // The PhaseSpace class.
16 : // Base class for phase space generators.
17 :
18 : //--------------------------------------------------------------------------
19 :
20 : // Constants: could be changed here if desired, but normally should not.
21 : // These are of technical nature, as described for each.
22 :
23 : // Number of trial maxima around which maximum search is performed.
24 : const int PhaseSpace::NMAXTRY = 2;
25 :
26 : // Number of three-body trials in phase space optimization.
27 : const int PhaseSpace::NTRY3BODY = 20;
28 :
29 : // Maximum cross section increase, just in case true maximum not found.
30 : const double PhaseSpace::SAFETYMARGIN = 1.05;
31 :
32 : // Small number to avoid division by zero.
33 : const double PhaseSpace::TINY = 1e-20;
34 :
35 : // Fraction of total weight that is shared evenly between all shapes.
36 : const double PhaseSpace::EVENFRAC = 0.4;
37 :
38 : // Two cross sections with a small relative error are assumed same.
39 : const double PhaseSpace::SAMESIGMA = 1e-6;
40 :
41 : // Do not include resonances peaked too far outside allowed mass region.
42 : const double PhaseSpace::WIDTHMARGIN = 20.;
43 :
44 : // Special optimization treatment when two resonances at almost same mass.
45 : const double PhaseSpace::SAMEMASS = 0.01;
46 :
47 : // Minimum phase space left when kinematics constraints are combined.
48 : const double PhaseSpace::MASSMARGIN = 0.01;
49 :
50 : // When using Breit-Wigners in 2 -> 2 raise maximum weight estimate.
51 : const double PhaseSpace::EXTRABWWTMAX = 1.25;
52 :
53 : // Size of Breit-Wigner threshold region, for mass selection biasing.
54 : const double PhaseSpace::THRESHOLDSIZE = 3.;
55 :
56 : // Step size in optimal-mass search, for mass selection biasing.
57 : const double PhaseSpace::THRESHOLDSTEP = 0.2;
58 :
59 : // Minimal rapidity range for allowed open range (in 2 -> 3).
60 : const double PhaseSpace::YRANGEMARGIN = 1e-6;
61 :
62 : // Cutoff for f_e^e at x < 1 - 10^{-10} to be used in phase space selection.
63 : // Note: the ...MIN quantities come from 1 - x_max or 1 - tau_max.
64 : const double PhaseSpace::LEPTONXMIN = 1e-10;
65 : const double PhaseSpace::LEPTONXMAX = 1. - 1e-10;
66 6 : const double PhaseSpace::LEPTONXLOGMIN = log(1e-10);
67 6 : const double PhaseSpace::LEPTONXLOGMAX = log(1. - 1e-10);
68 : const double PhaseSpace::LEPTONTAUMIN = 2e-10;
69 :
70 : // Safety to avoid division with unreasonably small value for z selection.
71 : const double PhaseSpace::SHATMINZ = 1.;
72 :
73 : // Regularization for small pT2min in z = cos(theta) selection.
74 : const double PhaseSpace::PT2RATMINZ = 0.0001;
75 :
76 : // These numbers are hardwired empirical parameters,
77 : // intended to speed up the M-generator.
78 : const double PhaseSpace::WTCORRECTION[11] = { 1., 1., 1.,
79 : 2., 5., 15., 60., 250., 1250., 7000., 50000. };
80 :
81 : // MBR: form factor appoximation with two exponents, [FFB1,FFB2] = GeV-2.
82 : const double PhaseSpace::FFA1 = 0.9;
83 : const double PhaseSpace::FFA2 = 0.1;
84 : const double PhaseSpace::FFB1 = 4.6;
85 : const double PhaseSpace::FFB2 = 0.6;
86 :
87 : //--------------------------------------------------------------------------
88 :
89 : // Perform simple initialization and store pointers.
90 :
91 : void PhaseSpace::init(bool isFirst, SigmaProcess* sigmaProcessPtrIn,
92 : Info* infoPtrIn, Settings* settingsPtrIn, ParticleData* particleDataPtrIn,
93 : Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
94 : Couplings* couplingsPtrIn, SigmaTotal* sigmaTotPtrIn,
95 : UserHooks* userHooksPtrIn) {
96 :
97 : // Store input pointers for future use.
98 0 : sigmaProcessPtr = sigmaProcessPtrIn;
99 0 : infoPtr = infoPtrIn;
100 0 : settingsPtr = settingsPtrIn;
101 0 : particleDataPtr = particleDataPtrIn;
102 0 : rndmPtr = rndmPtrIn;
103 0 : beamAPtr = beamAPtrIn;
104 0 : beamBPtr = beamBPtrIn;
105 0 : couplingsPtr = couplingsPtrIn;
106 0 : sigmaTotPtr = sigmaTotPtrIn;
107 0 : userHooksPtr = userHooksPtrIn;
108 :
109 : // Some commonly used beam information.
110 0 : idA = beamAPtr->id();
111 0 : idB = beamBPtr->id();
112 0 : mA = beamAPtr->m();
113 0 : mB = beamBPtr->m();
114 0 : eCM = infoPtr->eCM();
115 0 : s = eCM * eCM;
116 :
117 : // Flag if lepton beams, and if non-resolved ones.
118 0 : hasLeptonBeamA = beamAPtr->isLepton();
119 0 : hasLeptonBeamB = beamBPtr->isLepton();
120 0 : hasTwoLeptonBeams = hasLeptonBeamA && hasLeptonBeamB;
121 0 : hasOneLeptonBeam = (hasLeptonBeamA || hasLeptonBeamB) && !hasTwoLeptonBeams;
122 0 : bool hasPointLepton = (hasLeptonBeamA && beamAPtr->isUnresolved())
123 0 : || (hasLeptonBeamB && beamBPtr->isUnresolved());
124 0 : hasOnePointLepton = hasOneLeptonBeam && hasPointLepton;
125 0 : hasTwoPointLeptons = hasTwoLeptonBeams && hasPointLepton;
126 :
127 : // Standard phase space cuts.
128 0 : if (isFirst || settingsPtr->flag("PhaseSpace:sameForSecond")) {
129 0 : mHatGlobalMin = settingsPtr->parm("PhaseSpace:mHatMin");
130 0 : mHatGlobalMax = settingsPtr->parm("PhaseSpace:mHatMax");
131 0 : pTHatGlobalMin = settingsPtr->parm("PhaseSpace:pTHatMin");
132 0 : pTHatGlobalMax = settingsPtr->parm("PhaseSpace:pTHatMax");
133 :
134 : // Optionally separate phase space cuts for second hard process.
135 0 : } else {
136 0 : mHatGlobalMin = settingsPtr->parm("PhaseSpace:mHatMinSecond");
137 0 : mHatGlobalMax = settingsPtr->parm("PhaseSpace:mHatMaxSecond");
138 0 : pTHatGlobalMin = settingsPtr->parm("PhaseSpace:pTHatMinSecond");
139 0 : pTHatGlobalMax = settingsPtr->parm("PhaseSpace:pTHatMaxSecond");
140 : }
141 :
142 : // Cutoff against divergences at pT -> 0.
143 0 : pTHatMinDiverge = settingsPtr->parm("PhaseSpace:pTHatMinDiverge");
144 :
145 : // When to use Breit-Wigners.
146 0 : useBreitWigners = settingsPtr->flag("PhaseSpace:useBreitWigners");
147 0 : minWidthBreitWigners = settingsPtr->parm("PhaseSpace:minWidthBreitWigners");
148 :
149 : // Whether generation is with variable energy.
150 0 : doEnergySpread = settingsPtr->flag("Beams:allowMomentumSpread");
151 :
152 : // Flags for maximization information and violation handling.
153 0 : showSearch = settingsPtr->flag("PhaseSpace:showSearch");
154 0 : showViolation = settingsPtr->flag("PhaseSpace:showViolation");
155 0 : increaseMaximum = settingsPtr->flag("PhaseSpace:increaseMaximum");
156 :
157 : // Know whether a Z0 is pure Z0 or admixed with gamma*.
158 0 : gmZmodeGlobal = settingsPtr->mode("WeakZ0:gmZmode");
159 :
160 : // Flags if user should be allowed to reweight cross section.
161 0 : canModifySigma = (userHooksPtr != 0)
162 0 : ? userHooksPtr->canModifySigma() : false;
163 0 : canBiasSelection = (userHooksPtr != 0)
164 0 : ? userHooksPtr->canBiasSelection() : false;
165 :
166 : // Parameters for simplified reweighting of 2 -> 2 processes.
167 0 : canBias2Sel = settingsPtr->flag("PhaseSpace:bias2Selection");
168 0 : bias2SelPow = settingsPtr->parm("PhaseSpace:bias2SelectionPow");
169 0 : bias2SelRef = settingsPtr->parm("PhaseSpace:bias2SelectionRef");
170 :
171 : // Default event-specific kinematics properties.
172 0 : x1H = 1.;
173 0 : x2H = 1.;
174 0 : m3 = 0.;
175 0 : m4 = 0.;
176 0 : m5 = 0.;
177 0 : s3 = m3 * m3;
178 0 : s4 = m4 * m4;
179 0 : s5 = m5 * m5;
180 0 : mHat = eCM;
181 0 : sH = s;
182 0 : tH = 0.;
183 0 : uH = 0.;
184 0 : pTH = 0.;
185 0 : theta = 0.;
186 0 : phi = 0.;
187 0 : runBW3H = 1.;
188 0 : runBW4H = 1.;
189 0 : runBW5H = 1.;
190 :
191 : // Default cross section information.
192 0 : sigmaNw = 0.;
193 0 : sigmaMx = 0.;
194 0 : sigmaPos = 0.;
195 0 : sigmaNeg = 0.;
196 0 : newSigmaMx = false;
197 0 : biasWt = 1.;
198 :
199 0 : }
200 :
201 : //--------------------------------------------------------------------------
202 :
203 : // Allow for nonisotropic decays when ME's available.
204 :
205 : void PhaseSpace::decayKinematics( Event& process) {
206 :
207 : // Identify sets of sister partons.
208 : int iResEnd = 4;
209 0 : for (int iResBeg = 5; iResBeg < process.size(); ++iResBeg) {
210 0 : if (iResBeg <= iResEnd) continue;
211 : iResEnd = iResBeg;
212 0 : while ( iResEnd < process.size() - 1
213 0 : && process[iResEnd + 1].mother1() == process[iResBeg].mother1()
214 0 : && process[iResEnd + 1].mother2() == process[iResBeg].mother2() )
215 : ++iResEnd;
216 :
217 : // Check that at least one of them is a resonance.
218 : bool hasRes = false;
219 0 : for (int iRes = iResBeg; iRes <= iResEnd; ++iRes)
220 0 : if ( !process[iRes].isFinal() ) hasRes = true;
221 0 : if ( !hasRes ) continue;
222 :
223 : // Evaluate matrix element and decide whether to keep kinematics.
224 0 : double decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd);
225 0 : if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
226 : "Kinematics: negative angular weight");
227 0 : if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
228 : "Kinematics: angular weight above unity");
229 0 : while (decWt < rndmPtr->flat() ) {
230 :
231 : // Find resonances for which to redo decay angles.
232 0 : for (int iRes = iResBeg; iRes < process.size(); ++iRes) {
233 0 : if ( process[iRes].isFinal() ) continue;
234 : int iResMother = iRes;
235 0 : while (iResMother > iResEnd)
236 0 : iResMother = process[iResMother].mother1();
237 0 : if (iResMother < iResBeg) continue;
238 :
239 : // Do decay of this mother isotropically in phase space.
240 0 : decayKinematicsStep( process, iRes);
241 :
242 : // End loop over resonance decay chains.
243 0 : }
244 :
245 : // Ready to allow new test of matrix element.
246 0 : decWt = sigmaProcessPtr->weightDecay( process, iResBeg, iResEnd);
247 0 : if (decWt < 0.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
248 : "Kinematics: negative angular weight");
249 0 : if (decWt > 1.) infoPtr->errorMsg("Warning in PhaseSpace::decay"
250 : "Kinematics: angular weight above unity");
251 : }
252 :
253 : // End loop over sets of sister resonances/partons.
254 0 : }
255 :
256 0 : }
257 :
258 : //--------------------------------------------------------------------------
259 :
260 : // Reselect decay products momenta isotropically in phase space.
261 : // Does not redo secondary vertex position!
262 :
263 : void PhaseSpace::decayKinematicsStep( Event& process, int iRes) {
264 :
265 : // Multiplicity and mother mass and four-momentum.
266 0 : int i1 = process[iRes].daughter1();
267 0 : int mult = process[iRes].daughter2() + 1 - i1;
268 0 : double m0 = process[iRes].m();
269 0 : Vec4 pRes = process[iRes].p();
270 :
271 : // Description of two-body decays as simple special case.
272 0 : if (mult == 2) {
273 :
274 : // Products and product masses.
275 0 : int i2 = i1 + 1;
276 0 : double m1t = process[i1].m();
277 0 : double m2t = process[i2].m();
278 :
279 : // Energies and absolute momentum in the rest frame.
280 0 : double e1 = 0.5 * (m0*m0 + m1t*m1t - m2t*m2t) / m0;
281 0 : double e2 = 0.5 * (m0*m0 + m2t*m2t - m1t*m1t) / m0;
282 0 : double p12 = 0.5 * sqrtpos( (m0 - m1t - m2t) * (m0 + m1t + m2t)
283 0 : * (m0 + m1t - m2t) * (m0 - m1t + m2t) ) / m0;
284 :
285 : // Pick isotropic angles to give three-momentum.
286 0 : double cosTheta = 2. * rndmPtr->flat() - 1.;
287 0 : double sinTheta = sqrt(1. - cosTheta*cosTheta);
288 0 : double phi12 = 2. * M_PI * rndmPtr->flat();
289 0 : double pX = p12 * sinTheta * cos(phi12);
290 0 : double pY = p12 * sinTheta * sin(phi12);
291 0 : double pZ = p12 * cosTheta;
292 :
293 : // Fill four-momenta in mother rest frame and then boost to lab frame.
294 0 : Vec4 p1( pX, pY, pZ, e1);
295 0 : Vec4 p2( -pX, -pY, -pZ, e2);
296 0 : p1.bst( pRes );
297 0 : p2.bst( pRes );
298 :
299 : // Done for two-body decay.
300 0 : process[i1].p( p1 );
301 0 : process[i2].p( p2 );
302 : return;
303 0 : }
304 :
305 : // Description of three-body decays as semi-simple special case.
306 0 : if (mult == 3) {
307 :
308 : // Products and product masses.
309 0 : int i2 = i1 + 1;
310 0 : int i3 = i2 + 1;
311 0 : double m1t = process[i1].m();
312 0 : double m2t = process[i2].m();
313 0 : double m3t = process[i3].m();
314 0 : double mDiff = m0 - (m1t + m2t + m3t);
315 :
316 : // Kinematical limits for 2+3 mass. Maximum phase-space weight.
317 0 : double m23Min = m2t + m3t;
318 0 : double m23Max = m0 - m1t;
319 0 : double p1Max = 0.5 * sqrtpos( (m0 - m1t - m23Min)
320 0 : * (m0 + m1t + m23Min) * (m0 + m1t - m23Min)
321 0 : * (m0 - m1t + m23Min) ) / m0;
322 0 : double p23Max = 0.5 * sqrtpos( (m23Max - m2t - m3t)
323 0 : * (m23Max + m2t + m3t) * (m23Max + m2t - m3t)
324 0 : * (m23Max - m2t + m3t) ) / m23Max;
325 0 : double wtPSmax = 0.5 * p1Max * p23Max;
326 :
327 : // Pick an intermediate mass m23 flat in the allowed range.
328 : double wtPS, m23, p1Abs, p23Abs;
329 0 : do {
330 0 : m23 = m23Min + rndmPtr->flat() * mDiff;
331 :
332 : // Translate into relative momenta and find phase-space weight.
333 0 : p1Abs = 0.5 * sqrtpos( (m0 - m1t - m23) * (m0 + m1t + m23)
334 0 : * (m0 + m1t - m23) * (m0 - m1t + m23) ) / m0;
335 0 : p23Abs = 0.5 * sqrtpos( (m23 - m2t - m3t) * (m23 + m2t + m3t)
336 0 : * (m23 + m2t - m3t) * (m23 - m2t + m3t) ) / m23;
337 0 : wtPS = p1Abs * p23Abs;
338 :
339 : // If rejected, try again with new invariant masses.
340 0 : } while ( wtPS < rndmPtr->flat() * wtPSmax );
341 :
342 : // Set up m23 -> m2 + m3 isotropic in its rest frame.
343 0 : double cosTheta = 2. * rndmPtr->flat() - 1.;
344 0 : double sinTheta = sqrt(1. - cosTheta*cosTheta);
345 0 : double phi23 = 2. * M_PI * rndmPtr->flat();
346 0 : double pX = p23Abs * sinTheta * cos(phi23);
347 0 : double pY = p23Abs * sinTheta * sin(phi23);
348 0 : double pZ = p23Abs * cosTheta;
349 0 : double e2 = sqrt( m2t*m2t + p23Abs*p23Abs);
350 0 : double e3 = sqrt( m3t*m3t + p23Abs*p23Abs);
351 0 : Vec4 p2( pX, pY, pZ, e2);
352 0 : Vec4 p3( -pX, -pY, -pZ, e3);
353 :
354 : // Set up 0 -> 1 + 23 isotropic in its rest frame.
355 0 : cosTheta = 2. * rndmPtr->flat() - 1.;
356 0 : sinTheta = sqrt(1. - cosTheta*cosTheta);
357 0 : phi23 = 2. * M_PI * rndmPtr->flat();
358 0 : pX = p1Abs * sinTheta * cos(phi23);
359 0 : pY = p1Abs * sinTheta * sin(phi23);
360 0 : pZ = p1Abs * cosTheta;
361 0 : double e1 = sqrt( m1t*m1t + p1Abs*p1Abs);
362 0 : double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
363 0 : Vec4 p1( pX, pY, pZ, e1);
364 :
365 : // Boost 2 + 3 to the 0 rest frame and then boost to lab frame.
366 0 : Vec4 p23( -pX, -pY, -pZ, e23);
367 0 : p2.bst( p23 );
368 0 : p3.bst( p23 );
369 0 : p1.bst( pRes );
370 0 : p2.bst( pRes );
371 0 : p3.bst( pRes );
372 :
373 : // Done for three-body decay.
374 0 : process[i1].p( p1 );
375 0 : process[i2].p( p2 );
376 0 : process[i3].p( p3 );
377 : return;
378 0 : }
379 :
380 : // Do a multibody decay using the M-generator algorithm.
381 :
382 : // Set up masses and four-momenta in a vector, with mother in slot 0.
383 0 : vector<double> mProd;
384 0 : mProd.push_back( m0);
385 0 : for (int i = i1; i <= process[iRes].daughter2(); ++i)
386 0 : mProd.push_back( process[i].m() );
387 0 : vector<Vec4> pProd;
388 0 : pProd.push_back( pRes);
389 :
390 : // Sum of daughter masses.
391 0 : double mSum = mProd[1];
392 0 : for (int i = 2; i <= mult; ++i) mSum += mProd[i];
393 0 : double mDiff = m0 - mSum;
394 :
395 : // Begin setup of intermediate invariant masses.
396 0 : vector<double> mInv;
397 0 : for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
398 :
399 : // Calculate the maximum weight in the decay.
400 0 : double wtPSmax = 1. / WTCORRECTION[mult];
401 0 : double mMaxWT = mDiff + mProd[mult];
402 : double mMinWT = 0.;
403 0 : for (int i = mult - 1; i > 0; --i) {
404 0 : mMaxWT += mProd[i];
405 0 : mMinWT += mProd[i+1];
406 0 : double mNow = mProd[i];
407 0 : wtPSmax *= 0.5 * sqrtpos( (mMaxWT - mMinWT - mNow)
408 0 : * (mMaxWT + mMinWT + mNow) * (mMaxWT + mMinWT - mNow)
409 0 : * (mMaxWT - mMinWT + mNow) ) / mMaxWT;
410 : }
411 :
412 : // Begin loop to find the set of intermediate invariant masses.
413 0 : vector<double> rndmOrd;
414 : double wtPS;
415 0 : do {
416 : wtPS = 1.;
417 :
418 : // Find and order random numbers in descending order.
419 0 : rndmOrd.resize(0);
420 0 : rndmOrd.push_back(1.);
421 0 : for (int i = 1; i < mult - 1; ++i) {
422 0 : double rndm = rndmPtr->flat();
423 0 : rndmOrd.push_back(rndm);
424 0 : for (int j = i - 1; j > 0; --j) {
425 0 : if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
426 0 : else break;
427 : }
428 0 : }
429 0 : rndmOrd.push_back(0.);
430 :
431 : // Translate into intermediate masses and find weight.
432 0 : for (int i = mult - 1; i > 0; --i) {
433 0 : mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
434 0 : wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
435 0 : * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
436 0 : * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
437 : }
438 :
439 : // If rejected, try again with new invariant masses.
440 0 : } while ( wtPS < rndmPtr->flat() * wtPSmax );
441 :
442 : // Perform two-particle decays in the respective rest frame.
443 0 : vector<Vec4> pInv;
444 0 : pInv.resize(mult + 1);
445 0 : for (int i = 1; i < mult; ++i) {
446 0 : double p12 = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
447 0 : * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
448 0 : * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
449 :
450 : // Isotropic angles give three-momentum.
451 0 : double cosTheta = 2. * rndmPtr->flat() - 1.;
452 0 : double sinTheta = sqrt(1. - cosTheta*cosTheta);
453 0 : double phiLoc = 2. * M_PI * rndmPtr->flat();
454 0 : double pX = p12 * sinTheta * cos(phiLoc);
455 0 : double pY = p12 * sinTheta * sin(phiLoc);
456 0 : double pZ = p12 * cosTheta;
457 :
458 : // Calculate energies, fill four-momenta.
459 0 : double eHad = sqrt( mProd[i]*mProd[i] + p12*p12);
460 0 : double eInv = sqrt( mInv[i+1]*mInv[i+1] + p12*p12);
461 0 : pProd.push_back( Vec4( pX, pY, pZ, eHad) );
462 0 : pInv[i+1].p( -pX, -pY, -pZ, eInv);
463 : }
464 0 : pProd.push_back( pInv[mult] );
465 :
466 : // Boost decay products to the mother rest frame and on to lab frame.
467 0 : pInv[1] = pProd[0];
468 0 : for (int iFrame = mult - 1; iFrame > 0; --iFrame)
469 0 : for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);
470 :
471 : // Done for multibody decay.
472 0 : for (int i = 1; i < mult; ++i)
473 0 : process[i1 + i - 1].p( pProd[i] );
474 : return;
475 :
476 0 : }
477 :
478 : //--------------------------------------------------------------------------
479 :
480 : // Determine how 3-body phase space should be sampled.
481 :
482 : void PhaseSpace::setup3Body() {
483 :
484 : // Check for massive t-channel propagator particles.
485 0 : int idTchan1 = abs( sigmaProcessPtr->idTchan1() );
486 0 : int idTchan2 = abs( sigmaProcessPtr->idTchan2() );
487 0 : mTchan1 = (idTchan1 == 0) ? pTHatMinDiverge
488 0 : : particleDataPtr->m0(idTchan1);
489 0 : mTchan2 = (idTchan2 == 0) ? pTHatMinDiverge
490 0 : : particleDataPtr->m0(idTchan2);
491 0 : sTchan1 = mTchan1 * mTchan1;
492 0 : sTchan2 = mTchan2 * mTchan2;
493 :
494 : // Find coefficients of different pT2 selection terms. Mirror choice.
495 0 : frac3Pow1 = sigmaProcessPtr->tChanFracPow1();
496 0 : frac3Pow2 = sigmaProcessPtr->tChanFracPow2();
497 0 : frac3Flat = 1. - frac3Pow1 - frac3Pow2;
498 0 : useMirrorWeight = sigmaProcessPtr->useMirrorWeight();
499 :
500 0 : }
501 :
502 : //--------------------------------------------------------------------------
503 :
504 : // Determine how phase space should be sampled.
505 :
506 : bool PhaseSpace::setupSampling123(bool is2, bool is3, ostream& os) {
507 :
508 : // Optional printout.
509 0 : if (showSearch) os << "\n PYTHIA Optimization printout for "
510 0 : << sigmaProcessPtr->name() << "\n \n" << scientific << setprecision(3);
511 :
512 : // Check that open range in tau (+ set tauMin, tauMax).
513 0 : if (!limitTau(is2, is3)) return false;
514 :
515 : // Reset coefficients and matrices of equation system to solve.
516 0 : int binTau[8], binY[8], binZ[8];
517 0 : double vecTau[8], matTau[8][8], vecY[8], matY[8][8], vecZ[8], matZ[8][8];
518 0 : for (int i = 0; i < 8; ++i) {
519 0 : tauCoef[i] = 0.;
520 0 : yCoef[i] = 0.;
521 0 : zCoef[i] = 0.;
522 0 : binTau[i] = 0;
523 0 : binY[i] = 0;
524 0 : binZ[i] = 0;
525 0 : vecTau[i] = 0.;
526 0 : vecY[i] = 0.;
527 0 : vecZ[i] = 0.;
528 0 : for (int j = 0; j < 8; ++j) {
529 0 : matTau[i][j] = 0.;
530 0 : matY[i][j] = 0.;
531 0 : matZ[i][j] = 0.;
532 : }
533 : }
534 0 : sigmaMx = 0.;
535 0 : sigmaNeg = 0.;
536 :
537 : // Number of used coefficients/points for each dimension: tau, y, c.
538 0 : nTau = (hasTwoPointLeptons) ? 1 : 2;
539 0 : nY = (hasOnePointLepton || hasTwoPointLeptons) ? 1 : 5;
540 0 : nZ = (is2) ? 5 : 1;
541 :
542 : // Identify if any resonances contribute in s-channel.
543 0 : idResA = sigmaProcessPtr->resonanceA();
544 0 : if (idResA != 0) {
545 0 : mResA = particleDataPtr->m0(idResA);
546 0 : GammaResA = particleDataPtr->mWidth(idResA);
547 0 : if (mHatMin > mResA + WIDTHMARGIN * GammaResA || (mHatMax > 0.
548 0 : && mHatMax < mResA - WIDTHMARGIN * GammaResA) ) idResA = 0;
549 : }
550 0 : idResB = sigmaProcessPtr->resonanceB();
551 0 : if (idResB != 0) {
552 0 : mResB = particleDataPtr->m0(idResB);
553 0 : GammaResB = particleDataPtr->mWidth(idResB);
554 0 : if (mHatMin > mResB + WIDTHMARGIN * GammaResB || (mHatMax > 0.
555 0 : && mHatMax < mResB - WIDTHMARGIN * GammaResB) ) idResB = 0;
556 : }
557 0 : if (idResA == 0 && idResB != 0) {
558 0 : idResA = idResB;
559 0 : mResA = mResB;
560 0 : GammaResA = GammaResB;
561 0 : idResB = 0;
562 0 : }
563 :
564 : // More sampling in tau if resonances in s-channel.
565 0 : if (idResA !=0 && !hasTwoPointLeptons) {
566 0 : nTau += 2;
567 0 : tauResA = mResA * mResA / s;
568 0 : widResA = mResA * GammaResA / s;
569 0 : }
570 0 : if (idResB != 0 && !hasTwoPointLeptons) {
571 0 : nTau += 2;
572 0 : tauResB = mResB * mResB / s;
573 0 : widResB = mResB * GammaResB / s;
574 0 : }
575 :
576 : // More sampling in tau (and different in y) if incoming lepton beams.
577 0 : if (hasTwoLeptonBeams && !hasTwoPointLeptons) ++nTau;
578 :
579 : // Special case when both resonances have same mass.
580 0 : sameResMass = false;
581 0 : if (idResB != 0 && abs(mResA - mResB) < SAMEMASS * (GammaResA + GammaResB))
582 0 : sameResMass = true;
583 :
584 : // Default z value and weight required for 2 -> 1. Number of dimensions.
585 0 : z = 0.;
586 0 : wtZ = 1.;
587 0 : int nVar = (is2) ? 3 : 2;
588 :
589 : // Initial values, to be modified later.
590 0 : tauCoef[0] = 1.;
591 0 : yCoef[1] = 0.5;
592 0 : yCoef[2] = 0.5;
593 0 : zCoef[0] = 1.;
594 :
595 : // Step through grid in tau. Set limits on y and z generation.
596 0 : for (int iTau = 0; iTau < nTau; ++iTau) {
597 : double posTau = 0.5;
598 0 : if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6;
599 0 : selectTau( iTau, posTau, is2);
600 0 : if (!limitY()) continue;
601 0 : if (is2 && !limitZ()) continue;
602 :
603 : // Step through grids in y and z.
604 0 : for (int iY = 0; iY < nY; ++iY) {
605 0 : selectY( iY, 0.5);
606 0 : for (int iZ = 0; iZ < nZ; ++iZ) {
607 0 : if (is2) selectZ( iZ, 0.5);
608 : double sigmaTmp = 0.;
609 :
610 : // 2 -> 1: calculate cross section, weighted by phase-space volume.
611 0 : if (!is2 && !is3) {
612 0 : sigmaProcessPtr->set1Kin( x1H, x2H, sH);
613 0 : sigmaTmp = sigmaProcessPtr->sigmaPDF();
614 0 : sigmaTmp *= wtTau * wtY;
615 :
616 : // 2 -> 2: calculate cross section, weighted by phase-space volume
617 : // and Breit-Wigners for masses
618 0 : } else if (is2) {
619 0 : sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
620 0 : runBW3H, runBW4H);
621 0 : sigmaTmp = sigmaProcessPtr->sigmaPDF();
622 0 : sigmaTmp *= wtTau * wtY * wtZ * wtBW;
623 :
624 : // 2 -> 3: repeat internal 3-body phase space several times and
625 : // keep maximal cross section, weighted by phase-space volume
626 : // and Breit-Wigners for masses
627 0 : } else if (is3) {
628 0 : for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
629 0 : if (!select3Body()) continue;
630 0 : sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
631 0 : m3, m4, m5, runBW3H, runBW4H, runBW5H);
632 0 : double sigmaTry = sigmaProcessPtr->sigmaPDF();
633 0 : sigmaTry *= wtTau * wtY * wt3Body * wtBW;
634 0 : if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
635 0 : }
636 0 : }
637 :
638 : // Allow possibility for user to modify cross section. (3body??)
639 0 : if (canModifySigma) sigmaTmp
640 0 : *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
641 0 : if (canBiasSelection) sigmaTmp
642 0 : *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, false);
643 0 : if (canBias2Sel) sigmaTmp *= pow( pTH / bias2SelRef, bias2SelPow);
644 :
645 : // Check if current maximum exceeded.
646 0 : if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp;
647 :
648 : // Optional printout. Protect against negative cross sections.
649 0 : if (showSearch) os << " tau =" << setw(11) << tau << " y ="
650 0 : << setw(11) << y << " z =" << setw(11) << z
651 0 : << " sigma =" << setw(11) << sigmaTmp << "\n";
652 0 : if (sigmaTmp < 0.) sigmaTmp = 0.;
653 :
654 : // Sum up tau cross-section pieces in points used.
655 0 : if (!hasTwoPointLeptons) {
656 0 : binTau[iTau] += 1;
657 0 : vecTau[iTau] += sigmaTmp;
658 0 : matTau[iTau][0] += 1. / intTau0;
659 0 : matTau[iTau][1] += (1. / intTau1) / tau;
660 0 : if (idResA != 0) {
661 0 : matTau[iTau][2] += (1. / intTau2) / (tau + tauResA);
662 0 : matTau[iTau][3] += (1. / intTau3)
663 0 : * tau / ( pow2(tau - tauResA) + pow2(widResA) );
664 0 : }
665 0 : if (idResB != 0) {
666 0 : matTau[iTau][4] += (1. / intTau4) / (tau + tauResB);
667 0 : matTau[iTau][5] += (1. / intTau5)
668 0 : * tau / ( pow2(tau - tauResB) + pow2(widResB) );
669 0 : }
670 0 : if (hasTwoLeptonBeams) matTau[iTau][nTau - 1] += (1. / intTau6)
671 0 : * tau / max( LEPTONTAUMIN, 1. - tau);
672 : }
673 :
674 : // Sum up y cross-section pieces in points used.
675 0 : if (!hasOnePointLepton && !hasTwoPointLeptons) {
676 0 : binY[iY] += 1;
677 0 : vecY[iY] += sigmaTmp;
678 0 : matY[iY][0] += (yMax / intY0) / cosh(y);
679 0 : matY[iY][1] += (yMax / intY12) * (y + yMax);
680 0 : matY[iY][2] += (yMax / intY12) * (yMax - y);
681 0 : if (!hasTwoLeptonBeams) {
682 0 : matY[iY][3] += (yMax / intY34) * exp(y);
683 0 : matY[iY][4] += (yMax / intY34) * exp(-y);
684 0 : } else {
685 0 : matY[iY][3] += (yMax / intY56)
686 0 : / max( LEPTONXMIN, 1. - exp( y - yMax) );
687 0 : matY[iY][4] += (yMax / intY56)
688 0 : / max( LEPTONXMIN, 1. - exp(-y - yMax) );
689 : }
690 : }
691 :
692 : // Integrals over z expressions at tauMax, to be used below.
693 0 : if (is2) {
694 0 : double p2AbsMax = 0.25 * (pow2(tauMax * s - s3 - s4)
695 0 : - 4. * s3 * s4) / (tauMax * s);
696 0 : double zMaxMax = sqrtpos( 1. - pT2HatMin / p2AbsMax );
697 0 : double zPosMaxMax = max(ratio34, unity34 + zMaxMax);
698 0 : double zNegMaxMax = max(ratio34, unity34 - zMaxMax);
699 0 : double intZ0Max = 2. * zMaxMax;
700 0 : double intZ12Max = log( zPosMaxMax / zNegMaxMax);
701 0 : double intZ34Max = 1. / zNegMaxMax - 1. / zPosMaxMax;
702 :
703 : // Sum up z cross-section pieces in points used.
704 0 : binZ[iZ] += 1;
705 0 : vecZ[iZ] += sigmaTmp;
706 0 : matZ[iZ][0] += 1.;
707 0 : matZ[iZ][1] += (intZ0Max / intZ12Max) / zNeg;
708 0 : matZ[iZ][2] += (intZ0Max / intZ12Max) / zPos;
709 0 : matZ[iZ][3] += (intZ0Max / intZ34Max) / pow2(zNeg);
710 0 : matZ[iZ][4] += (intZ0Max / intZ34Max) / pow2(zPos);
711 0 : }
712 :
713 : // End of loops over phase space points.
714 : }
715 : }
716 0 : }
717 :
718 : // Fail if no non-vanishing cross sections.
719 0 : if (sigmaMx <= 0.) {
720 0 : sigmaMx = 0.;
721 0 : return false;
722 : }
723 :
724 : // Solve respective equation system for better phase space coefficients.
725 0 : if (!hasTwoPointLeptons) solveSys( nTau, binTau, vecTau, matTau, tauCoef);
726 0 : if (!hasOnePointLepton && !hasTwoPointLeptons)
727 0 : solveSys( nY, binY, vecY, matY, yCoef);
728 0 : if (is2) solveSys( nZ, binZ, vecZ, matZ, zCoef);
729 0 : if (showSearch) os << "\n";
730 :
731 : // Provide cumulative sum of coefficients.
732 0 : tauCoefSum[0] = tauCoef[0];
733 0 : yCoefSum[0] = yCoef[0];
734 0 : zCoefSum[0] = zCoef[0];
735 0 : for (int i = 1; i < 8; ++ i) {
736 0 : tauCoefSum[i] = tauCoefSum[i - 1] + tauCoef[i];
737 0 : yCoefSum[i] = yCoefSum[i - 1] + yCoef[i];
738 0 : zCoefSum[i] = zCoefSum[i - 1] + zCoef[i];
739 : }
740 : // The last element should be > 1 to be on safe side in selection below.
741 0 : tauCoefSum[nTau - 1] = 2.;
742 0 : yCoefSum[nY - 1] = 2.;
743 0 : zCoefSum[nZ - 1] = 2.;
744 :
745 :
746 : // Begin find two most promising maxima among same points as before.
747 0 : int iMaxTau[NMAXTRY + 2], iMaxY[NMAXTRY + 2], iMaxZ[NMAXTRY + 2];
748 0 : double sigMax[NMAXTRY + 2];
749 : int nMax = 0;
750 :
751 : // Scan same grid as before in tau, y, z.
752 0 : for (int iTau = 0; iTau < nTau; ++iTau) {
753 : double posTau = 0.5;
754 0 : if (sameResMass && iTau > 1 && iTau < 6) posTau = (iTau < 4) ? 0.4 : 0.6;
755 0 : selectTau( iTau, posTau, is2);
756 0 : if (!limitY()) continue;
757 0 : if (is2 && !limitZ()) continue;
758 0 : for (int iY = 0; iY < nY; ++iY) {
759 0 : selectY( iY, 0.5);
760 0 : for (int iZ = 0; iZ < nZ; ++iZ) {
761 0 : if (is2) selectZ( iZ, 0.5);
762 : double sigmaTmp = 0.;
763 :
764 : // 2 -> 1: calculate cross section, weighted by phase-space volume.
765 0 : if (!is2 && !is3) {
766 0 : sigmaProcessPtr->set1Kin( x1H, x2H, sH);
767 0 : sigmaTmp = sigmaProcessPtr->sigmaPDF();
768 0 : sigmaTmp *= wtTau * wtY;
769 :
770 : // 2 -> 2: calculate cross section, weighted by phase-space volume
771 : // and Breit-Wigners for masses
772 0 : } else if (is2) {
773 0 : sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
774 0 : runBW3H, runBW4H);
775 0 : sigmaTmp = sigmaProcessPtr->sigmaPDF();
776 0 : sigmaTmp *= wtTau * wtY * wtZ * wtBW;
777 :
778 : // 2 -> 3: repeat internal 3-body phase space several times and
779 : // keep maximal cross section, weighted by phase-space volume
780 : // and Breit-Wigners for masses
781 0 : } else if (is3) {
782 0 : for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
783 0 : if (!select3Body()) continue;
784 0 : sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
785 0 : m3, m4, m5, runBW3H, runBW4H, runBW5H);
786 0 : double sigmaTry = sigmaProcessPtr->sigmaPDF();
787 0 : sigmaTry *= wtTau * wtY * wt3Body * wtBW;
788 0 : if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
789 0 : }
790 0 : }
791 :
792 : // Allow possibility for user to modify cross section. (3body??)
793 0 : if (canModifySigma) sigmaTmp
794 0 : *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
795 0 : if (canBiasSelection) sigmaTmp
796 0 : *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, false);
797 0 : if (canBias2Sel) sigmaTmp *= pow( pTH / bias2SelRef, bias2SelPow);
798 :
799 : // Optional printout. Protect against negative cross section.
800 0 : if (showSearch) os << " tau =" << setw(11) << tau << " y ="
801 0 : << setw(11) << y << " z =" << setw(11) << z
802 0 : << " sigma =" << setw(11) << sigmaTmp << "\n";
803 0 : if (sigmaTmp < 0.) sigmaTmp = 0.;
804 :
805 : // Check that point is not simply mirror of already found one.
806 : bool mirrorPoint = false;
807 0 : for (int iMove = 0; iMove < nMax; ++iMove)
808 0 : if (abs(sigmaTmp - sigMax[iMove]) < SAMESIGMA
809 0 : * (sigmaTmp + sigMax[iMove])) mirrorPoint = true;
810 :
811 : // Add to or insert in maximum list. Only first two count.
812 0 : if (!mirrorPoint) {
813 : int iInsert = 0;
814 0 : for (int iMove = nMax - 1; iMove >= -1; --iMove) {
815 0 : iInsert = iMove + 1;
816 0 : if (iInsert == 0 || sigmaTmp < sigMax[iMove]) break;
817 0 : iMaxTau[iMove + 1] = iMaxTau[iMove];
818 0 : iMaxY[iMove + 1] = iMaxY[iMove];
819 0 : iMaxZ[iMove + 1] = iMaxZ[iMove];
820 0 : sigMax[iMove + 1] = sigMax[iMove];
821 : }
822 0 : iMaxTau[iInsert] = iTau;
823 0 : iMaxY[iInsert] = iY;
824 0 : iMaxZ[iInsert] = iZ;
825 0 : sigMax[iInsert] = sigmaTmp;
826 0 : if (nMax < NMAXTRY) ++nMax;
827 0 : }
828 :
829 : // Found two most promising maxima.
830 : }
831 : }
832 0 : }
833 0 : if (showSearch) os << "\n";
834 :
835 : // Read out starting position for search.
836 0 : sigmaMx = sigMax[0];
837 0 : int beginVar = (hasTwoPointLeptons) ? 2 : 0;
838 0 : if (hasOnePointLepton) beginVar = 1;
839 0 : for (int iMax = 0; iMax < nMax; ++iMax) {
840 0 : int iTau = iMaxTau[iMax];
841 0 : int iY = iMaxY[iMax];
842 0 : int iZ = iMaxZ[iMax];
843 : double tauVal = 0.5;
844 : double yVal = 0.5;
845 : double zVal = 0.5;
846 : int iGrid;
847 0 : double varVal, varNew, deltaVar, marginVar, sigGrid[3];
848 :
849 : // Starting point and step size in parameter space.
850 0 : for (int iRepeat = 0; iRepeat < 2; ++iRepeat) {
851 : // Run through (possibly a subset of) tau, y and z.
852 0 : for (int iVar = beginVar; iVar < nVar; ++iVar) {
853 0 : bool isTauVar = iVar == 0 || (beginVar == 1 && iVar == 1);
854 0 : if (isTauVar) varVal = tauVal;
855 0 : else if (iVar == 1) varVal = yVal;
856 : else varVal = zVal;
857 0 : deltaVar = (iRepeat == 0) ? 0.1
858 0 : : max( 0.01, min( 0.05, min( varVal - 0.02, 0.98 - varVal) ) );
859 0 : marginVar = (iRepeat == 0) ? 0.02 : 0.002;
860 0 : int moveStart = (iRepeat == 0 && isTauVar) ? 0 : 1;
861 0 : for (int move = moveStart; move < 9; ++move) {
862 :
863 : // Define new parameter-space point by step in one dimension.
864 0 : if (move == 0) {
865 : iGrid = 1;
866 : varNew = varVal;
867 0 : } else if (move == 1) {
868 : iGrid = 2;
869 0 : varNew = varVal + deltaVar;
870 0 : } else if (move == 2) {
871 : iGrid = 0;
872 0 : varNew = varVal - deltaVar;
873 0 : } else if (sigGrid[2] >= max( sigGrid[0], sigGrid[1])
874 0 : && varVal + 2. * deltaVar < 1. - marginVar) {
875 0 : varVal += deltaVar;
876 0 : sigGrid[0] = sigGrid[1];
877 0 : sigGrid[1] = sigGrid[2];
878 : iGrid = 2;
879 0 : varNew = varVal + deltaVar;
880 0 : } else if (sigGrid[0] >= max( sigGrid[1], sigGrid[2])
881 0 : && varVal - 2. * deltaVar > marginVar) {
882 0 : varVal -= deltaVar;
883 0 : sigGrid[2] = sigGrid[1];
884 0 : sigGrid[1] = sigGrid[0];
885 : iGrid = 0;
886 0 : varNew = varVal - deltaVar;
887 0 : } else if (sigGrid[2] >= sigGrid[0]) {
888 0 : deltaVar *= 0.5;
889 0 : varVal += deltaVar;
890 0 : sigGrid[0] = sigGrid[1];
891 : iGrid = 1;
892 : varNew = varVal;
893 0 : } else {
894 : deltaVar *= 0.5;
895 0 : varVal -= deltaVar;
896 0 : sigGrid[2] = sigGrid[1];
897 : iGrid = 1;
898 : varNew = varVal;
899 : }
900 :
901 : // Convert to relevant variables and find derived new limits.
902 : bool insideLimits = true;
903 0 : if (isTauVar) {
904 : tauVal = varNew;
905 0 : selectTau( iTau, tauVal, is2);
906 0 : if (!limitY()) insideLimits = false;
907 0 : if (is2 && !limitZ()) insideLimits = false;
908 0 : if (insideLimits) {
909 0 : selectY( iY, yVal);
910 0 : if (is2) selectZ( iZ, zVal);
911 : }
912 0 : } else if (iVar == 1) {
913 : yVal = varNew;
914 0 : selectY( iY, yVal);
915 0 : } else if (iVar == 2) {
916 : zVal = varNew;
917 0 : selectZ( iZ, zVal);
918 0 : }
919 :
920 : // Evaluate cross-section.
921 : double sigmaTmp = 0.;
922 0 : if (insideLimits) {
923 :
924 : // 2 -> 1: calculate cross section, weighted by phase-space volume.
925 0 : if (!is2 && !is3) {
926 0 : sigmaProcessPtr->set1Kin( x1H, x2H, sH);
927 0 : sigmaTmp = sigmaProcessPtr->sigmaPDF();
928 0 : sigmaTmp *= wtTau * wtY;
929 :
930 : // 2 -> 2: calculate cross section, weighted by phase-space volume
931 : // and Breit-Wigners for masses
932 0 : } else if (is2) {
933 0 : sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4,
934 0 : runBW3H, runBW4H);
935 0 : sigmaTmp = sigmaProcessPtr->sigmaPDF();
936 0 : sigmaTmp *= wtTau * wtY * wtZ * wtBW;
937 :
938 : // 2 -> 3: repeat internal 3-body phase space several times and
939 : // keep maximal cross section, weighted by phase-space volume
940 : // and Breit-Wigners for masses
941 0 : } else if (is3) {
942 0 : for (int iTry3 = 0; iTry3 < NTRY3BODY; ++iTry3) {
943 0 : if (!select3Body()) continue;
944 0 : sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
945 0 : m3, m4, m5, runBW3H, runBW4H, runBW5H);
946 0 : double sigmaTry = sigmaProcessPtr->sigmaPDF();
947 0 : sigmaTry *= wtTau * wtY * wt3Body * wtBW;
948 0 : if (sigmaTry > sigmaTmp) sigmaTmp = sigmaTry;
949 0 : }
950 0 : }
951 :
952 : // Allow possibility for user to modify cross section.
953 0 : if (canModifySigma) sigmaTmp
954 0 : *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, false);
955 0 : if (canBiasSelection) sigmaTmp
956 0 : *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, false);
957 0 : if (canBias2Sel) sigmaTmp *= pow( pTH / bias2SelRef, bias2SelPow);
958 :
959 : // Optional printout. Protect against negative cross section.
960 0 : if (showSearch) os << " tau =" << setw(11) << tau << " y ="
961 0 : << setw(11) << y << " z =" << setw(11) << z
962 0 : << " sigma =" << setw(11) << sigmaTmp << "\n";
963 0 : if (sigmaTmp < 0.) sigmaTmp = 0.;
964 : }
965 :
966 : // Save new maximum. Final maximum.
967 0 : sigGrid[iGrid] = sigmaTmp;
968 0 : if (sigmaTmp > sigmaMx) sigmaMx = sigmaTmp;
969 : }
970 : }
971 : }
972 0 : }
973 0 : sigmaMx *= SAFETYMARGIN;
974 0 : sigmaPos = sigmaMx;
975 :
976 : // Optional printout.
977 0 : if (showSearch) os << "\n Final maximum = " << setw(11) << sigmaMx << endl;
978 :
979 : // Done.
980 : return true;
981 0 : }
982 :
983 : //--------------------------------------------------------------------------
984 :
985 : // Select a trial kinematics phase space point.
986 : // Note: by In is meant the integral over the quantity multiplying
987 : // coefficient cn. The sum of cn is normalized to unity.
988 :
989 : bool PhaseSpace::trialKin123(bool is2, bool is3, bool inEvent, ostream& os) {
990 :
991 : // Allow for possibility that energy varies from event to event.
992 0 : if (doEnergySpread) {
993 0 : eCM = infoPtr->eCM();
994 0 : s = eCM * eCM;
995 :
996 : // Find shifted tauRes values.
997 0 : if (idResA !=0 && !hasTwoPointLeptons) {
998 0 : tauResA = mResA * mResA / s;
999 0 : widResA = mResA * GammaResA / s;
1000 0 : }
1001 0 : if (idResB != 0 && !hasTwoPointLeptons) {
1002 0 : tauResB = mResB * mResB / s;
1003 0 : widResB = mResB * GammaResB / s;
1004 0 : }
1005 : }
1006 :
1007 : // Choose tau according to h1(tau)/tau, where
1008 : // h1(tau) = c0/I0 + (c1/I1) * 1/tau
1009 : // + (c2/I2) / (tau + tauResA)
1010 : // + (c3/I3) * tau / ((tau - tauResA)^2 + widResA^2)
1011 : // + (c4/I4) / (tau + tauResB)
1012 : // + (c5/I5) * tau / ((tau - tauResB)^2 + widResB^2)
1013 : // + (c6/I6) * tau / (1 - tau).
1014 0 : if (!limitTau(is2, is3)) return false;
1015 : int iTau = 0;
1016 0 : if (!hasTwoPointLeptons) {
1017 0 : double rTau = rndmPtr->flat();
1018 0 : while (rTau > tauCoefSum[iTau]) ++iTau;
1019 0 : }
1020 0 : selectTau( iTau, rndmPtr->flat(), is2);
1021 :
1022 : // Choose y according to h2(y), where
1023 : // h2(y) = (c0/I0) * 1/cosh(y)
1024 : // + (c1/I1) * (y-ymin) + (c2/I2) * (ymax-y)
1025 : // + (c3/I3) * exp(y) + (c4/i4) * exp(-y) (for hadron; for lepton instead)
1026 : // + (c5/I5) * 1 / (1 - exp(y-ymax)) + (c6/I6) * 1 / (1 - exp(ymin-y)).
1027 0 : if (!limitY()) return false;
1028 : int iY = 0;
1029 0 : if (!hasOnePointLepton && !hasTwoPointLeptons) {
1030 0 : double rY = rndmPtr->flat();
1031 0 : while (rY > yCoefSum[iY]) ++iY;
1032 0 : }
1033 0 : selectY( iY, rndmPtr->flat());
1034 :
1035 : // Choose z = cos(thetaHat) according to h3(z), where
1036 : // h3(z) = c0/I0 + (c1/I1) * 1/(A - z) + (c2/I2) * 1/(A + z)
1037 : // + (c3/I3) * 1/(A - z)^2 + (c4/I4) * 1/(A + z)^2,
1038 : // where A = 1 + 2*(m3*m4/sH)^2 (= 1 for massless products).
1039 0 : if (is2) {
1040 0 : if (!limitZ()) return false;
1041 : int iZ = 0;
1042 0 : double rZ = rndmPtr->flat();
1043 0 : while (rZ > zCoefSum[iZ]) ++iZ;
1044 0 : selectZ( iZ, rndmPtr->flat());
1045 0 : }
1046 :
1047 : // 2 -> 1: calculate cross section, weighted by phase-space volume.
1048 0 : if (!is2 && !is3) {
1049 0 : sigmaProcessPtr->set1Kin( x1H, x2H, sH);
1050 0 : sigmaNw = sigmaProcessPtr->sigmaPDF();
1051 0 : sigmaNw *= wtTau * wtY;
1052 :
1053 : // 2 -> 2: calculate cross section, weighted by phase-space volume
1054 : // and Breit-Wigners for masses
1055 0 : } else if (is2) {
1056 0 : sigmaProcessPtr->set2Kin( x1H, x2H, sH, tH, m3, m4, runBW3H, runBW4H);
1057 0 : sigmaNw = sigmaProcessPtr->sigmaPDF();
1058 0 : sigmaNw *= wtTau * wtY * wtZ * wtBW;
1059 :
1060 : // 2 -> 3: also sample internal 3-body phase, weighted by
1061 : // 2 -> 1 phase-space volume and Breit-Wigners for masses
1062 0 : } else if (is3) {
1063 0 : if (!select3Body()) sigmaNw = 0.;
1064 : else {
1065 0 : sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
1066 0 : m3, m4, m5, runBW3H, runBW4H, runBW5H);
1067 0 : sigmaNw = sigmaProcessPtr->sigmaPDF();
1068 0 : sigmaNw *= wtTau * wtY * wt3Body * wtBW;
1069 : }
1070 : }
1071 :
1072 : // Allow possibility for user to modify cross section.
1073 0 : if (canModifySigma) sigmaNw
1074 0 : *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, inEvent);
1075 0 : if (canBiasSelection) sigmaNw
1076 0 : *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, inEvent);
1077 0 : if (canBias2Sel) sigmaNw *= pow( pTH / bias2SelRef, bias2SelPow);
1078 :
1079 : // Check if maximum violated.
1080 0 : newSigmaMx = false;
1081 0 : if (sigmaNw > sigmaMx) {
1082 0 : infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin: "
1083 : "maximum for cross section violated");
1084 :
1085 : // Violation strategy 1: increase maximum (always during initialization).
1086 0 : if (increaseMaximum || !inEvent) {
1087 0 : double violFact = SAFETYMARGIN * sigmaNw / sigmaMx;
1088 0 : sigmaMx = SAFETYMARGIN * sigmaNw;
1089 0 : newSigmaMx = true;
1090 0 : if (showViolation) {
1091 0 : if (violFact < 9.99) os << fixed;
1092 0 : else os << scientific;
1093 0 : os << " PYTHIA Maximum for " << sigmaProcessPtr->name()
1094 0 : << " increased by factor " << setprecision(3) << violFact
1095 0 : << " to " << scientific << sigmaMx << endl;
1096 0 : }
1097 :
1098 : // Violation strategy 2: weight event (done in ProcessContainer).
1099 0 : } else if (showViolation && sigmaNw > sigmaPos) {
1100 0 : double violFact = sigmaNw / sigmaMx;
1101 0 : if (violFact < 9.99) os << fixed;
1102 0 : else os << scientific;
1103 0 : os << " PYTHIA Maximum for " << sigmaProcessPtr->name()
1104 0 : << " exceeded by factor " << setprecision(3) << violFact << endl;
1105 0 : sigmaPos = sigmaNw;
1106 0 : }
1107 : }
1108 :
1109 : // Check if negative cross section.
1110 0 : if (sigmaNw < sigmaNeg) {
1111 0 : infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::trialKin:"
1112 0 : " negative cross section set 0", "for " + sigmaProcessPtr->name() );
1113 0 : sigmaNeg = sigmaNw;
1114 :
1115 : // Optional printout of (all) violations.
1116 0 : if (showViolation) os << " PYTHIA Negative minimum for "
1117 0 : << sigmaProcessPtr->name() << " changed to " << scientific
1118 0 : << setprecision(3) << sigmaNeg << endl;
1119 : }
1120 0 : if (sigmaNw < 0.) sigmaNw = 0.;
1121 :
1122 : // Set event weight, where relevant.
1123 0 : biasWt = (canBiasSelection) ? userHooksPtr->biasedSelectionWeight() : 1.;
1124 0 : if (canBias2Sel) biasWt /= pow( pTH / bias2SelRef, bias2SelPow);
1125 :
1126 : // Done.
1127 0 : return true;
1128 0 : }
1129 :
1130 : //--------------------------------------------------------------------------
1131 :
1132 : // Find range of allowed tau values.
1133 :
1134 : bool PhaseSpace::limitTau(bool is2, bool is3) {
1135 :
1136 : // Trivial reply for unresolved lepton beams.
1137 0 : if (hasTwoPointLeptons) {
1138 0 : tauMin = 1.;
1139 0 : tauMax = 1.;
1140 0 : return true;
1141 : }
1142 :
1143 : // Requirements from allowed mHat range.
1144 0 : tauMin = sHatMin / s;
1145 0 : tauMax = (mHatMax < mHatMin) ? 1. : min( 1., sHatMax / s);
1146 :
1147 : // Requirements from allowed pT range and masses.
1148 0 : if (is2 || is3) {
1149 0 : double mT3Min = sqrt(s3 + pT2HatMin);
1150 0 : double mT4Min = sqrt(s4 + pT2HatMin);
1151 0 : double mT5Min = (is3) ? sqrt(s5 + pT2HatMin) : 0.;
1152 0 : tauMin = max( tauMin, pow2(mT3Min + mT4Min + mT5Min) / s);
1153 0 : }
1154 :
1155 : // Check that there is an open range.
1156 0 : return (tauMax > tauMin);
1157 0 : }
1158 :
1159 : //--------------------------------------------------------------------------
1160 :
1161 : // Find range of allowed y values.
1162 :
1163 : bool PhaseSpace::limitY() {
1164 :
1165 : // Trivial reply for unresolved lepton beams.
1166 0 : if (hasTwoPointLeptons) {
1167 0 : yMax = 1.;
1168 0 : return true;
1169 : }
1170 :
1171 : // Requirements from selected tau value. Trivial for one unresolved beam.
1172 0 : yMax = -0.5 * log(tau);
1173 0 : if (hasOnePointLepton) return true;
1174 :
1175 : // For lepton beams requirements from cutoff for f_e^e.
1176 0 : double yMaxMargin = (hasTwoLeptonBeams) ? yMax + LEPTONXLOGMAX : yMax;
1177 :
1178 : // Check that there is an open range.
1179 0 : return (yMaxMargin > 0.);
1180 0 : }
1181 :
1182 : //--------------------------------------------------------------------------
1183 :
1184 : // Find range of allowed z = cos(theta) values.
1185 :
1186 : bool PhaseSpace::limitZ() {
1187 :
1188 : // Default limits.
1189 0 : zMin = 0.;
1190 0 : zMax = 1.;
1191 :
1192 : // Requirements from pTHat limits.
1193 0 : zMax = sqrtpos( 1. - pT2HatMin / p2Abs );
1194 0 : if (pTHatMax > pTHatMin) zMin = sqrtpos( 1. - pT2HatMax / p2Abs );
1195 :
1196 : // Check that there is an open range.
1197 0 : return (zMax > zMin);
1198 : }
1199 :
1200 : //--------------------------------------------------------------------------
1201 :
1202 : // Select tau according to a choice of shapes.
1203 :
1204 : void PhaseSpace::selectTau(int iTau, double tauVal, bool is2) {
1205 :
1206 : // Trivial reply for unresolved lepton beams.
1207 0 : if (hasTwoPointLeptons) {
1208 0 : tau = 1.;
1209 0 : wtTau = 1.;
1210 0 : sH = s;
1211 0 : mHat = sqrt(sH);
1212 0 : if (is2) {
1213 0 : p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1214 0 : pAbs = sqrtpos( p2Abs );
1215 0 : }
1216 : return;
1217 : }
1218 :
1219 : // Contributions from s-channel resonances.
1220 : double tRatA = 0.;
1221 : double aLowA = 0.;
1222 : double aUppA = 0.;
1223 0 : if (idResA !=0) {
1224 0 : tRatA = ((tauResA + tauMax) / (tauResA + tauMin)) * (tauMin / tauMax);
1225 0 : aLowA = atan( (tauMin - tauResA) / widResA);
1226 0 : aUppA = atan( (tauMax - tauResA) / widResA);
1227 0 : }
1228 : double tRatB = 0.;
1229 : double aLowB = 0.;
1230 : double aUppB = 0.;
1231 0 : if (idResB != 0) {
1232 0 : tRatB = ((tauResB + tauMax) / (tauResB + tauMin)) * (tauMin / tauMax);
1233 0 : aLowB = atan( (tauMin - tauResB) / widResB);
1234 0 : aUppB = atan( (tauMax - tauResB) / widResB);
1235 0 : }
1236 :
1237 : // Contributions from 1 / (1 - tau) for lepton beams.
1238 : double aLowT = 0.;
1239 : double aUppT = 0.;
1240 0 : if (hasTwoLeptonBeams) {
1241 0 : aLowT = log( max( LEPTONTAUMIN, 1. - tauMin) );
1242 0 : aUppT = log( max( LEPTONTAUMIN, 1. - tauMax) );
1243 0 : intTau6 = aLowT - aUppT;
1244 0 : }
1245 :
1246 : // Select according to 1/tau or 1/tau^2.
1247 0 : if (iTau == 0) tau = tauMin * pow( tauMax / tauMin, tauVal);
1248 0 : else if (iTau == 1) tau = tauMax * tauMin
1249 0 : / (tauMin + (tauMax - tauMin) * tauVal);
1250 :
1251 : // Select according to 1 / (1 - tau) for lepton beams.
1252 0 : else if (hasTwoLeptonBeams && iTau == nTau - 1)
1253 0 : tau = 1. - exp( aUppT + intTau6 * tauVal );
1254 :
1255 : // Select according to 1 / (tau * (tau + tauRes)) or
1256 : // 1 / ((tau - tauRes)^2 + widRes^2) for resonances A and B.
1257 0 : else if (iTau == 2) tau = tauResA * tauMin
1258 0 : / ((tauResA + tauMin) * pow( tRatA, tauVal) - tauMin);
1259 0 : else if (iTau == 3) tau = tauResA + widResA
1260 0 : * tan( aLowA + (aUppA - aLowA) * tauVal);
1261 0 : else if (iTau == 4) tau = tauResB * tauMin
1262 0 : / ((tauResB + tauMin) * pow( tRatB, tauVal) - tauMin);
1263 0 : else if (iTau == 5) tau = tauResB + widResB
1264 0 : * tan( aLowB + (aUppB - aLowB) * tauVal);
1265 :
1266 : // Phase space weight in tau.
1267 0 : intTau0 = log( tauMax / tauMin);
1268 0 : intTau1 = (tauMax - tauMin) / (tauMax * tauMin);
1269 0 : double invWtTau = (tauCoef[0] / intTau0) + (tauCoef[1] / intTau1) / tau;
1270 0 : if (idResA != 0) {
1271 0 : intTau2 = -log(tRatA) / tauResA;
1272 0 : intTau3 = (aUppA - aLowA) / widResA;
1273 0 : invWtTau += (tauCoef[2] / intTau2) / (tau + tauResA)
1274 0 : + (tauCoef[3] / intTau3) * tau / ( pow2(tau - tauResA) + pow2(widResA) );
1275 0 : }
1276 0 : if (idResB != 0) {
1277 0 : intTau4 = -log(tRatB) / tauResB;
1278 0 : intTau5 = (aUppB - aLowB) / widResB;
1279 0 : invWtTau += (tauCoef[4] / intTau4) / (tau + tauResB)
1280 0 : + (tauCoef[5] / intTau5) * tau / ( pow2(tau - tauResB) + pow2(widResB) );
1281 0 : }
1282 0 : if (hasTwoLeptonBeams)
1283 0 : invWtTau += (tauCoef[nTau - 1] / intTau6)
1284 0 : * tau / max( LEPTONTAUMIN, 1. - tau);
1285 0 : wtTau = 1. / invWtTau;
1286 :
1287 : // Calculate sHat and absolute momentum of outgoing partons.
1288 0 : sH = tau * s;
1289 0 : mHat = sqrt(sH);
1290 0 : if (is2) {
1291 0 : p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
1292 0 : pAbs = sqrtpos( p2Abs );
1293 0 : }
1294 :
1295 0 : }
1296 :
1297 : //--------------------------------------------------------------------------
1298 :
1299 : // Select y according to a choice of shapes.
1300 :
1301 : void PhaseSpace::selectY(int iY, double yVal) {
1302 :
1303 : // Trivial reply for two unresolved lepton beams.
1304 0 : if (hasTwoPointLeptons) {
1305 0 : y = 0.;
1306 0 : wtY = 1.;
1307 0 : x1H = 1.;
1308 0 : x2H = 1.;
1309 0 : return;
1310 : }
1311 :
1312 : // Trivial replies for one unresolved lepton beam.
1313 0 : if (hasOnePointLepton) {
1314 0 : if (hasLeptonBeamA) {
1315 0 : y = yMax;
1316 0 : x1H = 1.;
1317 0 : x2H = tau;
1318 0 : } else {
1319 0 : y = -yMax;
1320 0 : x1H = tau;
1321 0 : x2H = 1.;
1322 : }
1323 0 : wtY = 1.;
1324 0 : return;
1325 : }
1326 :
1327 : // For lepton beams skip options 3&4 and go straight to 5&6.
1328 0 : if (hasTwoLeptonBeams && iY > 2) iY += 2;
1329 :
1330 : // Standard expressions used below.
1331 0 : double expYMax = exp( yMax );
1332 0 : double expYMin = exp(-yMax );
1333 0 : double atanMax = atan( expYMax );
1334 0 : double atanMin = atan( expYMin );
1335 0 : double aUppY = (hasTwoLeptonBeams)
1336 0 : ? log( max( LEPTONXMIN, LEPTONXMAX / tau - 1. ) ) : 0.;
1337 0 : double aLowY = LEPTONXLOGMIN;
1338 :
1339 : // 1 / cosh(y).
1340 0 : if (iY == 0) y = log( tan( atanMin + (atanMax - atanMin) * yVal ) );
1341 :
1342 : // y - y_min or mirrored y_max - y.
1343 0 : else if (iY <= 2) y = yMax * (2. * sqrt(yVal) - 1.);
1344 :
1345 : // exp(y) or mirrored exp(-y).
1346 0 : else if (iY <= 4) y = log( expYMin + (expYMax - expYMin) * yVal );
1347 :
1348 : // 1 / (1 - exp(y - y_max)) or mirrored 1 / (1 - exp(y_min - y)).
1349 0 : else y = yMax - log( 1. + exp(aLowY + (aUppY - aLowY) * yVal) );
1350 :
1351 : // Mirror two cases.
1352 0 : if (iY == 2 || iY == 4 || iY == 6) y = -y;
1353 :
1354 : // Phase space integral in y.
1355 0 : intY0 = 2. * (atanMax - atanMin);
1356 0 : intY12 = 0.5 * pow2(2. * yMax);
1357 0 : intY34 = expYMax - expYMin;
1358 0 : intY56 = aUppY - aLowY;
1359 0 : double invWtY = (yCoef[0] / intY0) / cosh(y)
1360 0 : + (yCoef[1] / intY12) * (y + yMax) + (yCoef[2] / intY12) * (yMax - y);
1361 0 : if (!hasTwoLeptonBeams) invWtY
1362 0 : += (yCoef[3] / intY34) * exp(y) + (yCoef[4] / intY34) * exp(-y);
1363 : else invWtY
1364 0 : += (yCoef[3] / intY56) / max( LEPTONXMIN, 1. - exp( y - yMax) )
1365 0 : + (yCoef[4] / intY56) / max( LEPTONXMIN, 1. - exp(-y - yMax) );
1366 0 : wtY = 1. / invWtY;
1367 :
1368 : // Calculate x1 and x2.
1369 0 : x1H = sqrt(tau) * exp(y);
1370 0 : x2H = sqrt(tau) * exp(-y);
1371 0 : }
1372 :
1373 : //--------------------------------------------------------------------------
1374 :
1375 : // Select z = cos(theta) according to a choice of shapes.
1376 : // The selection is split in the positive- and negative-z regions,
1377 : // since a pTmax cut can remove the region around z = 0.
1378 :
1379 : void PhaseSpace::selectZ(int iZ, double zVal) {
1380 :
1381 : // Mass-dependent dampening of pT -> 0 limit.
1382 0 : ratio34 = max(TINY, 2. * s3 * s4 / pow2(sH));
1383 0 : unity34 = 1. + ratio34;
1384 0 : double ratiopT2 = 2. * pT2HatMin / max( SHATMINZ, sH);
1385 0 : if (ratiopT2 < PT2RATMINZ) ratio34 = max( ratio34, ratiopT2);
1386 :
1387 : // Common expressions in z limits.
1388 0 : double zPosMax = max(ratio34, unity34 + zMax);
1389 0 : double zNegMax = max(ratio34, unity34 - zMax);
1390 0 : double zPosMin = max(ratio34, unity34 + zMin);
1391 0 : double zNegMin = max(ratio34, unity34 - zMin);
1392 :
1393 : // Flat in z.
1394 0 : if (iZ == 0) {
1395 0 : if (zVal < 0.5) z = -(zMax + (zMin - zMax) * 2. * zVal);
1396 0 : else z = zMin + (zMax - zMin) * (2. * zVal - 1.);
1397 :
1398 : // 1 / (unity34 - z).
1399 0 : } else if (iZ == 1) {
1400 0 : double areaNeg = log(zPosMax / zPosMin);
1401 0 : double areaPos = log(zNegMin / zNegMax);
1402 0 : double area = areaNeg + areaPos;
1403 0 : if (zVal * area < areaNeg) {
1404 0 : double zValMod = zVal * area / areaNeg;
1405 0 : z = unity34 - zPosMax * pow(zPosMin / zPosMax, zValMod);
1406 0 : } else {
1407 0 : double zValMod = (zVal * area - areaNeg)/ areaPos;
1408 0 : z = unity34 - zNegMin * pow(zNegMax / zNegMin, zValMod);
1409 : }
1410 :
1411 : // 1 / (unity34 + z).
1412 0 : } else if (iZ == 2) {
1413 0 : double areaNeg = log(zNegMin / zNegMax);
1414 0 : double areaPos = log(zPosMax / zPosMin);
1415 0 : double area = areaNeg + areaPos;
1416 0 : if (zVal * area < areaNeg) {
1417 0 : double zValMod = zVal * area / areaNeg;
1418 0 : z = zNegMax * pow(zNegMin / zNegMax, zValMod) - unity34;
1419 0 : } else {
1420 0 : double zValMod = (zVal * area - areaNeg)/ areaPos;
1421 0 : z = zPosMin * pow(zPosMax / zPosMin, zValMod) - unity34;
1422 : }
1423 :
1424 : // 1 / (unity34 - z)^2.
1425 0 : } else if (iZ == 3) {
1426 0 : double areaNeg = 1. / zPosMin - 1. / zPosMax;
1427 0 : double areaPos = 1. / zNegMax - 1. / zNegMin;
1428 0 : double area = areaNeg + areaPos;
1429 0 : if (zVal * area < areaNeg) {
1430 0 : double zValMod = zVal * area / areaNeg;
1431 0 : z = unity34 - 1. / (1./zPosMax + areaNeg * zValMod);
1432 0 : } else {
1433 0 : double zValMod = (zVal * area - areaNeg)/ areaPos;
1434 0 : z = unity34 - 1. / (1./zNegMin + areaPos * zValMod);
1435 : }
1436 :
1437 : // 1 / (unity34 + z)^2.
1438 0 : } else if (iZ == 4) {
1439 0 : double areaNeg = 1. / zNegMax - 1. / zNegMin;
1440 0 : double areaPos = 1. / zPosMin - 1. / zPosMax;
1441 0 : double area = areaNeg + areaPos;
1442 0 : if (zVal * area < areaNeg) {
1443 0 : double zValMod = zVal * area / areaNeg;
1444 0 : z = 1. / (1./zNegMax - areaNeg * zValMod) - unity34;
1445 0 : } else {
1446 0 : double zValMod = (zVal * area - areaNeg)/ areaPos;
1447 0 : z = 1. / (1./zPosMin - areaPos * zValMod) - unity34;
1448 : }
1449 0 : }
1450 :
1451 : // Safety check for roundoff errors. Combinations with z.
1452 0 : if (z < 0.) z = min(-zMin, max(-zMax, z));
1453 0 : else z = min(zMax, max(zMin, z));
1454 0 : zNeg = max(ratio34, unity34 - z);
1455 0 : zPos = max(ratio34, unity34 + z);
1456 :
1457 : // Phase space integral in z.
1458 0 : double intZ0 = 2. * (zMax - zMin);
1459 0 : double intZ12 = log( (zPosMax * zNegMin) / (zPosMin * zNegMax) );
1460 0 : double intZ34 = 1. / zPosMin - 1. / zPosMax + 1. / zNegMax
1461 0 : - 1. / zNegMin;
1462 0 : wtZ = mHat * pAbs / ( (zCoef[0] / intZ0) + (zCoef[1] / intZ12) / zNeg
1463 0 : + (zCoef[2] / intZ12) / zPos + (zCoef[3] / intZ34) / pow2(zNeg)
1464 0 : + (zCoef[4] / intZ34) / pow2(zPos) );
1465 :
1466 : // Calculate tHat and uHat. Also gives pTHat.
1467 0 : double sH34 = -0.5 * (sH - s3 - s4);
1468 0 : tH = sH34 + mHat * pAbs * z;
1469 0 : uH = sH34 - mHat * pAbs * z;
1470 0 : pTH = sqrtpos( (tH * uH - s3 * s4) / sH);
1471 :
1472 0 : }
1473 :
1474 : //--------------------------------------------------------------------------
1475 :
1476 : // Select three-body phase space according to a cylindrically based form
1477 : // that can be chosen to favour low pT based on the form of propagators.
1478 :
1479 : bool PhaseSpace::select3Body() {
1480 :
1481 : // Upper and lower limits of pT choice for 4 and 5.
1482 0 : double m35S = pow2(m3 + m5);
1483 0 : double pT4Smax = 0.25 * ( pow2(sH - s4 - m35S) - 4. * s4 * m35S ) / sH;
1484 0 : if (pTHatMax > pTHatMin) pT4Smax = min( pT2HatMax, pT4Smax);
1485 0 : double pT4Smin = pT2HatMin;
1486 0 : double m34S = pow2(m3 + m4);
1487 0 : double pT5Smax = 0.25 * ( pow2(sH - s5 - m34S) - 4. * s5 * m34S ) / sH;
1488 0 : if (pTHatMax > pTHatMin) pT5Smax = min( pT2HatMax, pT5Smax);
1489 0 : double pT5Smin = pT2HatMin;
1490 :
1491 : // Check that pT ranges not closed.
1492 0 : if ( pT4Smax < pow2(pTHatMin + MASSMARGIN) ) return false;
1493 0 : if ( pT5Smax < pow2(pTHatMin + MASSMARGIN) ) return false;
1494 :
1495 : // Select pT4S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2.
1496 0 : double pTSmaxProp = pT4Smax + sTchan1;
1497 0 : double pTSminProp = pT4Smin + sTchan1;
1498 0 : double pTSratProp = pTSmaxProp / pTSminProp;
1499 0 : double pTSdiff = pT4Smax - pT4Smin;
1500 0 : double rShape = rndmPtr->flat();
1501 : double pT4S = 0.;
1502 0 : if (rShape < frac3Flat) pT4S = pT4Smin + rndmPtr->flat() * pTSdiff;
1503 0 : else if (rShape < frac3Flat + frac3Pow1) pT4S = max( pT2HatMin,
1504 0 : pTSminProp * pow( pTSratProp, rndmPtr->flat() ) - sTchan1 );
1505 0 : else pT4S = max( pT2HatMin, pTSminProp * pTSmaxProp
1506 0 : / (pTSminProp + rndmPtr->flat()* pTSdiff) - sTchan1 );
1507 0 : double wt4 = pTSdiff / ( frac3Flat
1508 0 : + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT4S + sTchan1))
1509 0 : + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT4S + sTchan1) );
1510 :
1511 : // Select pT5S according to c0 + c1/(M^2 + pT^2) + c2/(M^2 + pT^2)^2.
1512 0 : pTSmaxProp = pT5Smax + sTchan2;
1513 0 : pTSminProp = pT5Smin + sTchan2;
1514 0 : pTSratProp = pTSmaxProp / pTSminProp;
1515 0 : pTSdiff = pT5Smax - pT5Smin;
1516 0 : rShape = rndmPtr->flat();
1517 : double pT5S = 0.;
1518 0 : if (rShape < frac3Flat) pT5S = pT5Smin + rndmPtr->flat() * pTSdiff;
1519 0 : else if (rShape < frac3Flat + frac3Pow1) pT5S = max( pT2HatMin,
1520 0 : pTSminProp * pow( pTSratProp, rndmPtr->flat() ) - sTchan2 );
1521 0 : else pT5S = max( pT2HatMin, pTSminProp * pTSmaxProp
1522 0 : / (pTSminProp + rndmPtr->flat()* pTSdiff) - sTchan2 );
1523 0 : double wt5 = pTSdiff / ( frac3Flat
1524 0 : + frac3Pow1 * pTSdiff / (log(pTSratProp) * (pT5S + sTchan2))
1525 0 : + frac3Pow2 * pTSminProp * pTSmaxProp / pow2(pT5S + sTchan2) );
1526 :
1527 : // Select azimuthal angles and check that third pT in range.
1528 0 : double phi4 = 2. * M_PI * rndmPtr->flat();
1529 0 : double phi5 = 2. * M_PI * rndmPtr->flat();
1530 0 : double pT3S = max( 0., pT4S + pT5S + 2. * sqrt(pT4S * pT5S)
1531 0 : * cos(phi4 - phi5) );
1532 0 : if ( pT3S < pT2HatMin || (pTHatMax > pTHatMin && pT3S > pT2HatMax) )
1533 0 : return false;
1534 :
1535 : // Calculate transverse masses and check that phase space not closed.
1536 0 : double sT3 = s3 + pT3S;
1537 0 : double sT4 = s4 + pT4S;
1538 0 : double sT5 = s5 + pT5S;
1539 0 : double mT3 = sqrt(sT3);
1540 0 : double mT4 = sqrt(sT4);
1541 0 : double mT5 = sqrt(sT5);
1542 0 : if ( mT3 + mT4 + mT5 + MASSMARGIN > mHat ) return false;
1543 :
1544 : // Select rapidity for particle 3 and check that phase space not closed.
1545 0 : double m45S = pow2(mT4 + mT5);
1546 0 : double y3max = log( ( sH + sT3 - m45S + sqrtpos( pow2(sH - sT3 - m45S)
1547 0 : - 4 * sT3 * m45S ) ) / (2. * mHat * mT3) );
1548 0 : if (y3max < YRANGEMARGIN) return false;
1549 0 : double y3 = (2. * rndmPtr->flat() - 1.) * (1. - YRANGEMARGIN) * y3max;
1550 0 : double pz3 = mT3 * sinh(y3);
1551 0 : double e3 = mT3 * cosh(y3);
1552 :
1553 : // Find momentum transfers in the two mirror solutions (in 4-5 frame).
1554 0 : double pz45 = -pz3;
1555 0 : double e45 = mHat - e3;
1556 0 : double sT45 = e45 * e45 - pz45 * pz45;
1557 0 : double lam45 = sqrtpos( pow2(sT45 - sT4 - sT5) - 4. * sT4 * sT5 );
1558 0 : if (lam45 < YRANGEMARGIN * sH) return false;
1559 0 : double lam4e = sT45 + sT4 - sT5;
1560 0 : double lam5e = sT45 + sT5 - sT4;
1561 0 : double tFac = -0.5 * mHat / sT45;
1562 0 : double t1Pos = tFac * (e45 - pz45) * (lam4e - lam45);
1563 0 : double t1Neg = tFac * (e45 - pz45) * (lam4e + lam45);
1564 0 : double t2Pos = tFac * (e45 + pz45) * (lam5e - lam45);
1565 0 : double t2Neg = tFac * (e45 + pz45) * (lam5e + lam45);
1566 :
1567 : // Construct relative mirror weights and make choice.
1568 : double wtPosUnnorm = 1.;
1569 : double wtNegUnnorm = 1.;
1570 0 : if (useMirrorWeight) {
1571 0 : wtPosUnnorm = 1./ pow2( (t1Pos - sTchan1) * (t2Pos - sTchan2) );
1572 0 : wtNegUnnorm = 1./ pow2( (t1Neg - sTchan1) * (t2Neg - sTchan2) );
1573 0 : }
1574 0 : double wtPos = wtPosUnnorm / (wtPosUnnorm + wtNegUnnorm);
1575 0 : double wtNeg = wtNegUnnorm / (wtPosUnnorm + wtNegUnnorm);
1576 0 : double epsilon = (rndmPtr->flat() < wtPos) ? 1. : -1.;
1577 :
1578 : // Construct four-vectors in rest frame of subprocess.
1579 0 : double px4 = sqrt(pT4S) * cos(phi4);
1580 0 : double py4 = sqrt(pT4S) * sin(phi4);
1581 0 : double px5 = sqrt(pT5S) * cos(phi5);
1582 0 : double py5 = sqrt(pT5S) * sin(phi5);
1583 0 : double pz4 = 0.5 * (pz45 * lam4e + epsilon * e45 * lam45) / sT45;
1584 0 : double pz5 = pz45 - pz4;
1585 0 : double e4 = sqrt(sT4 + pz4 * pz4);
1586 0 : double e5 = sqrt(sT5 + pz5 * pz5);
1587 0 : p3cm = Vec4( -(px4 + px5), -(py4 + py5), pz3, e3);
1588 0 : p4cm = Vec4( px4, py4, pz4, e4);
1589 0 : p5cm = Vec4( px5, py5, pz5, e5);
1590 :
1591 : // Total weight to associate with kinematics choice.
1592 0 : wt3Body = wt4 * wt5 * (2. * y3max) / (128. * pow3(M_PI) * lam45);
1593 0 : wt3Body *= (epsilon > 0.) ? 1. / wtPos : 1. / wtNeg;
1594 :
1595 : // Cross section of form |M|^2/(2 sHat) dPS_3 so need 1/(2 sHat).
1596 0 : wt3Body /= (2. * sH);
1597 :
1598 : // Done.
1599 : return true;
1600 :
1601 0 : }
1602 :
1603 : //--------------------------------------------------------------------------
1604 :
1605 : // Solve linear equation system for better phase space coefficients.
1606 :
1607 : void PhaseSpace::solveSys( int n, int bin[8], double vec[8],
1608 : double mat[8][8], double coef[8], ostream& os) {
1609 :
1610 : // Optional printout.
1611 0 : if (showSearch) {
1612 0 : os << "\n Equation system: " << setw(5) << bin[0];
1613 0 : for (int j = 0; j < n; ++j) os << setw(12) << mat[0][j];
1614 0 : os << setw(12) << vec[0] << "\n";
1615 0 : for (int i = 1; i < n; ++i) {
1616 0 : os << " " << setw(5) << bin[i];
1617 0 : for (int j = 0; j < n; ++j) os << setw(12) << mat[i][j];
1618 0 : os << setw(12) << vec[i] << "\n";
1619 : }
1620 0 : }
1621 :
1622 : // Local variables.
1623 0 : double vecNor[8], coefTmp[8];
1624 0 : for (int i = 0; i < n; ++i) coefTmp[i] = 0.;
1625 :
1626 : // Check if equation system solvable.
1627 : bool canSolve = true;
1628 0 : for (int i = 0; i < n; ++i) if (bin[i] == 0) canSolve = false;
1629 : double vecSum = 0.;
1630 0 : for (int i = 0; i < n; ++i) vecSum += vec[i];
1631 0 : if (abs(vecSum) < TINY) canSolve = false;
1632 :
1633 : // Solve to find relative importance of cross-section pieces.
1634 0 : if (canSolve) {
1635 0 : for (int i = 0; i < n; ++i) vecNor[i] = max( 0.1, vec[i] / vecSum);
1636 0 : for (int k = 0; k < n - 1; ++k) {
1637 0 : for (int i = k + 1; i < n; ++i) {
1638 0 : if (abs(mat[k][k]) < TINY) {canSolve = false; break;}
1639 0 : double ratio = mat[i][k] / mat[k][k];
1640 0 : vec[i] -= ratio * vec[k];
1641 0 : for (int j = k; j < n; ++j) mat[i][j] -= ratio * mat[k][j];
1642 : }
1643 0 : if (!canSolve) break;
1644 : }
1645 0 : if (canSolve) {
1646 0 : for (int k = n - 1; k >= 0; --k) {
1647 0 : for (int j = k + 1; j < n; ++j) vec[k] -= mat[k][j] * coefTmp[j];
1648 0 : coefTmp[k] = vec[k] / mat[k][k];
1649 : }
1650 0 : }
1651 : }
1652 :
1653 : // Share evenly if failure.
1654 0 : if (!canSolve) for (int i = 0; i < n; ++i) {
1655 0 : coefTmp[i] = 1.;
1656 0 : vecNor[i] = 0.1;
1657 0 : if (vecSum > TINY) vecNor[i] = max(0.1, vec[i] / vecSum);
1658 0 : }
1659 :
1660 : // Normalize coefficients, with piece shared democratically.
1661 : double coefSum = 0.;
1662 : vecSum = 0.;
1663 0 : for (int i = 0; i < n; ++i) {
1664 0 : coefTmp[i] = max( 0., coefTmp[i]);
1665 0 : coefSum += coefTmp[i];
1666 0 : vecSum += vecNor[i];
1667 : }
1668 0 : if (coefSum > 0.) for (int i = 0; i < n; ++i) coef[i] = EVENFRAC / n
1669 0 : + (1. - EVENFRAC) * 0.5 * (coefTmp[i] / coefSum + vecNor[i] / vecSum);
1670 0 : else for (int i = 0; i < n; ++i) coef[i] = 1. / n;
1671 :
1672 : // Optional printout.
1673 0 : if (showSearch) {
1674 0 : os << " Solution: ";
1675 0 : for (int i = 0; i < n; ++i) os << setw(12) << coef[i];
1676 0 : os << "\n";
1677 0 : }
1678 0 : }
1679 :
1680 : //--------------------------------------------------------------------------
1681 :
1682 : // Setup mass selection for one resonance at a time - part 1.
1683 :
1684 : void PhaseSpace::setupMass1(int iM) {
1685 :
1686 : // Identity for mass seletion; is 0 also for light quarks (not yet selected).
1687 0 : if (iM == 3) idMass[iM] = abs(sigmaProcessPtr->id3Mass());
1688 0 : if (iM == 4) idMass[iM] = abs(sigmaProcessPtr->id4Mass());
1689 0 : if (iM == 5) idMass[iM] = abs(sigmaProcessPtr->id5Mass());
1690 :
1691 : // Masses and widths of resonances.
1692 0 : if (idMass[iM] == 0) {
1693 0 : mPeak[iM] = 0.;
1694 0 : mWidth[iM] = 0.;
1695 0 : mMin[iM] = 0.;
1696 0 : mMax[iM] = 0.;
1697 0 : } else {
1698 0 : mPeak[iM] = particleDataPtr->m0(idMass[iM]);
1699 0 : mWidth[iM] = particleDataPtr->mWidth(idMass[iM]);
1700 0 : mMin[iM] = particleDataPtr->mMin(idMass[iM]);
1701 0 : mMax[iM] = particleDataPtr->mMax(idMass[iM]);
1702 : // gmZmode == 1 means pure photon propagator; set at lower mass limit.
1703 0 : if (idMass[iM] == 23 && gmZmode == 1) mPeak[iM] = mMin[iM];
1704 : }
1705 :
1706 : // Mass and width combinations for Breit-Wigners.
1707 0 : sPeak[iM] = mPeak[iM] * mPeak[iM];
1708 0 : useBW[iM] = useBreitWigners && (mWidth[iM] > minWidthBreitWigners);
1709 0 : if (!useBW[iM]) mWidth[iM] = 0.;
1710 0 : mw[iM] = mPeak[iM] * mWidth[iM];
1711 0 : wmRat[iM] = (idMass[iM] == 0 || mPeak[iM] == 0.)
1712 0 : ? 0. : mWidth[iM] / mPeak[iM];
1713 :
1714 : // Simple Breit-Wigner range, upper edge to be corrected subsequently.
1715 0 : if (useBW[iM]) {
1716 0 : mLower[iM] = mMin[iM];
1717 0 : mUpper[iM] = mHatMax;
1718 0 : }
1719 :
1720 0 : }
1721 :
1722 : //--------------------------------------------------------------------------
1723 :
1724 : // Setup mass selection for one resonance at a time - part 2.
1725 :
1726 : void PhaseSpace::setupMass2(int iM, double distToThresh) {
1727 :
1728 : // Store reduced Breit-Wigner range.
1729 0 : if (mMax[iM] > mMin[iM]) mUpper[iM] = min( mUpper[iM], mMax[iM]);
1730 0 : sLower[iM] = mLower[iM] * mLower[iM];
1731 0 : sUpper[iM] = mUpper[iM] * mUpper[iM];
1732 :
1733 : // Prepare to select m3 by BW + flat + 1/s_3.
1734 : // Determine relative coefficients by allowed mass range.
1735 0 : if (distToThresh > THRESHOLDSIZE) {
1736 0 : fracFlat[iM] = 0.1;
1737 0 : fracInv[iM] = 0.1;
1738 0 : } else if (distToThresh > - THRESHOLDSIZE) {
1739 0 : fracFlat[iM] = 0.25 - 0.15 * distToThresh / THRESHOLDSIZE;
1740 0 : fracInv [iM] = 0.15 - 0.05 * distToThresh / THRESHOLDSIZE;
1741 0 : } else {
1742 0 : fracFlat[iM] = 0.4;
1743 0 : fracInv[iM] = 0.2;
1744 : }
1745 :
1746 : // For gamma*/Z0: increase 1/s_i part and introduce 1/s_i^2 part.
1747 0 : fracInv2[iM] = 0.;
1748 0 : if (idMass[iM] == 23 && gmZmode == 0) {
1749 0 : fracFlat[iM] *= 0.5;
1750 0 : fracInv[iM] = 0.5 * fracInv[iM] + 0.25;
1751 0 : fracInv2[iM] = 0.25;
1752 0 : } else if (idMass[iM] == 23 && gmZmode == 1) {
1753 0 : fracFlat[iM] = 0.1;
1754 0 : fracInv[iM] = 0.4;
1755 0 : fracInv2[iM] = 0.4;
1756 0 : }
1757 :
1758 : // Normalization integrals for the respective contribution.
1759 0 : atanLower[iM] = atan( (sLower[iM] - sPeak[iM])/ mw[iM] );
1760 0 : atanUpper[iM] = atan( (sUpper[iM] - sPeak[iM])/ mw[iM] );
1761 0 : intBW[iM] = atanUpper[iM] - atanLower[iM];
1762 0 : intFlat[iM] = sUpper[iM] - sLower[iM];
1763 0 : intInv[iM] = log( sUpper[iM] / sLower[iM] );
1764 0 : intInv2[iM] = 1./sLower[iM] - 1./sUpper[iM];
1765 :
1766 0 : }
1767 :
1768 : //--------------------------------------------------------------------------
1769 :
1770 : // Select Breit-Wigner-distributed or fixed masses.
1771 :
1772 : void PhaseSpace::trialMass(int iM) {
1773 :
1774 : // References to masses to be set.
1775 0 : double& mSet = (iM == 3) ? m3 : ( (iM == 4) ? m4 : m5 );
1776 0 : double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 );
1777 :
1778 : // Distribution for m_i is BW + flat + 1/s_i + 1/s_i^2.
1779 0 : if (useBW[iM]) {
1780 0 : double pickForm = rndmPtr->flat();
1781 0 : if (pickForm > fracFlat[iM] + fracInv[iM] + fracInv2[iM])
1782 0 : sSet = sPeak[iM] + mw[iM] * tan( atanLower[iM]
1783 0 : + rndmPtr->flat() * intBW[iM] );
1784 0 : else if (pickForm > fracInv[iM] + fracInv2[iM])
1785 0 : sSet = sLower[iM] + rndmPtr->flat() * (sUpper[iM] - sLower[iM]);
1786 0 : else if (pickForm > fracInv2[iM])
1787 0 : sSet = sLower[iM] * pow( sUpper[iM] / sLower[iM], rndmPtr->flat() );
1788 0 : else sSet = sLower[iM] * sUpper[iM]
1789 0 : / (sLower[iM] + rndmPtr->flat() * (sUpper[iM] - sLower[iM]));
1790 0 : mSet = sqrt(sSet);
1791 :
1792 : // Else m_i is fixed at peak value.
1793 0 : } else {
1794 0 : mSet = mPeak[iM];
1795 0 : sSet = sPeak[iM];
1796 : }
1797 :
1798 0 : }
1799 :
1800 : //--------------------------------------------------------------------------
1801 :
1802 : // Naively a fixed-width Breit-Wigner is used to pick the mass.
1803 : // Here come the correction factors for
1804 : // (i) preselection according to BW + flat in s_i + 1/s_i + 1/s_i^2,
1805 : // (ii) reduced allowed mass range,
1806 : // (iii) running width, i.e. m0*Gamma0 -> s*Gamma0/m0.
1807 : // In the end, the weighted distribution is a running-width BW.
1808 :
1809 : double PhaseSpace::weightMass(int iM) {
1810 :
1811 : // Reference to mass and to Breit-Wigner weight to be set.
1812 0 : double& sSet = (iM == 3) ? s3 : ( (iM == 4) ? s4 : s5 );
1813 0 : double& runBWH = (iM == 3) ? runBW3H : ( (iM == 4) ? runBW4H : runBW5H );
1814 :
1815 : // Default weight if no Breit-Wigner.
1816 0 : runBWH = 1.;
1817 0 : if (!useBW[iM]) return 1.;
1818 :
1819 : // Weight of generated distribution.
1820 0 : double genBW = (1. - fracFlat[iM] - fracInv[iM] - fracInv2[iM])
1821 0 : * mw[iM] / ( (pow2(sSet - sPeak[iM]) + pow2(mw[iM])) * intBW[iM])
1822 0 : + fracFlat[iM] / intFlat[iM] + fracInv[iM] / (sSet * intInv[iM])
1823 0 : + fracInv2[iM] / (sSet*sSet * intInv2[iM]);
1824 :
1825 : // Weight of distribution with running width in Breit-Wigner.
1826 0 : double mwRun = sSet * wmRat[iM];
1827 0 : runBWH = mwRun / (pow2(sSet - sPeak[iM]) + pow2(mwRun)) / M_PI;
1828 :
1829 : // Done.
1830 0 : return (runBWH / genBW);
1831 :
1832 0 : }
1833 :
1834 : //==========================================================================
1835 :
1836 : // PhaseSpace2to1tauy class.
1837 : // 2 -> 1 kinematics for normal subprocesses.
1838 :
1839 : //--------------------------------------------------------------------------
1840 :
1841 : // Set limits for resonance mass selection.
1842 :
1843 : bool PhaseSpace2to1tauy::setupMass() {
1844 :
1845 : // Treat Z0 as such or as gamma*/Z0
1846 0 : gmZmode = gmZmodeGlobal;
1847 0 : int gmZmodeProc = sigmaProcessPtr->gmZmode();
1848 0 : if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
1849 :
1850 : // Mass limits for current resonance.
1851 0 : int idRes = abs(sigmaProcessPtr->resonanceA());
1852 0 : int idTmp = abs(sigmaProcessPtr->resonanceB());
1853 0 : if (idTmp > 0) idRes = idTmp;
1854 0 : double mResMin = (idRes == 0) ? 0. : particleDataPtr->mMin(idRes);
1855 0 : double mResMax = (idRes == 0) ? 0. : particleDataPtr->mMax(idRes);
1856 :
1857 : // Compare with global mass limits and pick tighter of them.
1858 0 : mHatMin = max( mResMin, mHatGlobalMin);
1859 0 : sHatMin = mHatMin*mHatMin;
1860 0 : mHatMax = eCM;
1861 0 : if (mResMax > mResMin) mHatMax = min( mHatMax, mResMax);
1862 0 : if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( mHatMax, mHatGlobalMax);
1863 0 : sHatMax = mHatMax*mHatMax;
1864 :
1865 : // Default Breit-Wigner weight.
1866 0 : wtBW = 1.;
1867 :
1868 : // Fail if mass window (almost) closed.
1869 0 : return (mHatMax > mHatMin + MASSMARGIN);
1870 :
1871 0 : }
1872 :
1873 : //--------------------------------------------------------------------------
1874 :
1875 : // Construct the four-vector kinematics from the trial values.
1876 :
1877 : bool PhaseSpace2to1tauy::finalKin() {
1878 :
1879 : // Particle masses; incoming always on mass shell.
1880 0 : mH[1] = 0.;
1881 0 : mH[2] = 0.;
1882 0 : mH[3] = mHat;
1883 :
1884 : // Incoming partons along beam axes. Outgoing has sum of momenta.
1885 0 : pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
1886 0 : pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
1887 0 : pH[3] = pH[1] + pH[2];
1888 :
1889 : // Done.
1890 0 : return true;
1891 : }
1892 :
1893 : //==========================================================================
1894 :
1895 : // PhaseSpace2to2tauyz class.
1896 : // 2 -> 2 kinematics for normal subprocesses.
1897 :
1898 : //--------------------------------------------------------------------------
1899 :
1900 : // Set up for fixed or Breit-Wigner mass selection.
1901 :
1902 : bool PhaseSpace2to2tauyz::setupMasses() {
1903 :
1904 : // Treat Z0 as such or as gamma*/Z0
1905 0 : gmZmode = gmZmodeGlobal;
1906 0 : int gmZmodeProc = sigmaProcessPtr->gmZmode();
1907 0 : if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
1908 :
1909 : // Set sHat limits - based on global limits only.
1910 0 : mHatMin = mHatGlobalMin;
1911 0 : sHatMin = mHatMin*mHatMin;
1912 0 : mHatMax = eCM;
1913 0 : if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax);
1914 0 : sHatMax = mHatMax*mHatMax;
1915 :
1916 : // Masses and widths of resonances.
1917 0 : setupMass1(3);
1918 0 : setupMass1(4);
1919 :
1920 : // Reduced mass range when two massive particles.
1921 0 : if (useBW[3]) mUpper[3] -= (useBW[4]) ? mMin[4] : mPeak[4];
1922 0 : if (useBW[4]) mUpper[4] -= (useBW[3]) ? mMin[3] : mPeak[3];
1923 :
1924 : // If closed phase space then unallowed process.
1925 : bool physical = true;
1926 0 : if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false;
1927 0 : if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false;
1928 0 : if (!useBW[3] && !useBW[4] && mHatMax < mPeak[3] + mPeak[4] + MASSMARGIN)
1929 0 : physical = false;
1930 0 : if (!physical) return false;
1931 :
1932 : // If either particle is massless then need extra pTHat cut.
1933 0 : pTHatMin = pTHatGlobalMin;
1934 0 : if (mPeak[3] < pTHatMinDiverge || mPeak[4] < pTHatMinDiverge)
1935 0 : pTHatMin = max( pTHatMin, pTHatMinDiverge);
1936 0 : pT2HatMin = pTHatMin * pTHatMin;
1937 0 : pTHatMax = pTHatGlobalMax;
1938 0 : pT2HatMax = pTHatMax * pTHatMax;
1939 :
1940 : // Prepare to select m3 by BW + flat + 1/s_3.
1941 0 : if (useBW[3]) {
1942 0 : double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[3]
1943 0 : / (pow2(mWidth[3]) + pow2(mWidth[4]));
1944 0 : double distToThreshB = (mHatMax - mPeak[3] - mMin[4]) / mWidth[3];
1945 0 : double distToThresh = min( distToThreshA, distToThreshB);
1946 0 : setupMass2(3, distToThresh);
1947 0 : }
1948 :
1949 : // Prepare to select m4 by BW + flat + 1/s_4.
1950 0 : if (useBW[4]) {
1951 0 : double distToThreshA = (mHatMax - mPeak[3] - mPeak[4]) * mWidth[4]
1952 0 : / (pow2(mWidth[3]) + pow2(mWidth[4]));
1953 0 : double distToThreshB = (mHatMax - mMin[3] - mPeak[4]) / mWidth[4];
1954 0 : double distToThresh = min( distToThreshA, distToThreshB);
1955 0 : setupMass2(4, distToThresh);
1956 0 : }
1957 :
1958 : // Initialization masses. Special cases when constrained phase space.
1959 0 : m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3];
1960 0 : m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4];
1961 0 : if (m3 + m4 + THRESHOLDSIZE * (mWidth[3] + mWidth[4]) + MASSMARGIN
1962 0 : > mHatMax) {
1963 0 : if (useBW[3] && useBW[4]) physical = constrainedM3M4();
1964 0 : else if (useBW[3]) physical = constrainedM3();
1965 0 : else if (useBW[4]) physical = constrainedM4();
1966 : }
1967 0 : s3 = m3*m3;
1968 0 : s4 = m4*m4;
1969 :
1970 : // Correct selected mass-spectrum to running-width Breit-Wigner.
1971 : // Extra safety margin for maximum search.
1972 0 : wtBW = 1.;
1973 0 : if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX;
1974 0 : if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX;
1975 :
1976 : // Done.
1977 0 : return physical;
1978 :
1979 0 : }
1980 :
1981 :
1982 : //--------------------------------------------------------------------------
1983 :
1984 : // Select Breit-Wigner-distributed or fixed masses.
1985 :
1986 : bool PhaseSpace2to2tauyz::trialMasses() {
1987 :
1988 : // By default vanishing cross section.
1989 0 : sigmaNw = 0.;
1990 0 : wtBW = 1.;
1991 :
1992 : // Pick m3 and m4 independently.
1993 0 : trialMass(3);
1994 0 : trialMass(4);
1995 :
1996 : // If outside phase space then reject event.
1997 0 : if (m3 + m4 + MASSMARGIN > mHatMax) return false;
1998 :
1999 : // Correct selected mass-spectrum to running-width Breit-Wigner.
2000 0 : if (useBW[3]) wtBW *= weightMass(3);
2001 0 : if (useBW[4]) wtBW *= weightMass(4);
2002 :
2003 : // Done.
2004 0 : return true;
2005 0 : }
2006 :
2007 : //--------------------------------------------------------------------------
2008 :
2009 : // Construct the four-vector kinematics from the trial values.
2010 :
2011 : bool PhaseSpace2to2tauyz::finalKin() {
2012 :
2013 : // Assign masses to particles assumed massless in matrix elements.
2014 0 : int id3 = sigmaProcessPtr->id(3);
2015 0 : int id4 = sigmaProcessPtr->id(4);
2016 0 : if (idMass[3] == 0) { m3 = particleDataPtr->m0(id3); s3 = m3*m3; }
2017 0 : if (idMass[4] == 0) { m4 = particleDataPtr->m0(id4); s4 = m4*m4; }
2018 :
2019 : // Sometimes swap tHat <-> uHat to reflect chosen final-state order.
2020 0 : if (sigmaProcessPtr->swappedTU()) {
2021 0 : swap(tH, uH);
2022 0 : z = -z;
2023 0 : }
2024 :
2025 : // Check that phase space still open after new mass assignment.
2026 0 : if (m3 + m4 + MASSMARGIN > mHat) {
2027 0 : infoPtr->errorMsg("Warning in PhaseSpace2to2tauyz::finalKin: "
2028 : "failed after mass assignment");
2029 0 : return false;
2030 : }
2031 0 : p2Abs = 0.25 * (pow2(sH - s3 - s4) - 4. * s3 * s4) / sH;
2032 0 : pAbs = sqrtpos( p2Abs );
2033 :
2034 : // Particle masses; incoming always on mass shell.
2035 0 : mH[1] = 0.;
2036 0 : mH[2] = 0.;
2037 0 : mH[3] = m3;
2038 0 : mH[4] = m4;
2039 :
2040 : // Incoming partons along beam axes.
2041 0 : pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
2042 0 : pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
2043 :
2044 : // Outgoing partons initially in collision CM frame along beam axes.
2045 0 : pH[3] = Vec4( 0., 0., pAbs, 0.5 * (sH + s3 - s4) / mHat);
2046 0 : pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (sH + s4 - s3) / mHat);
2047 :
2048 : // Then rotate and boost them to overall CM frame.
2049 0 : theta = acos(z);
2050 0 : phi = 2. * M_PI * rndmPtr->flat();
2051 0 : betaZ = (x1H - x2H)/(x1H + x2H);
2052 0 : pH[3].rot( theta, phi);
2053 0 : pH[4].rot( theta, phi);
2054 0 : pH[3].bst( 0., 0., betaZ);
2055 0 : pH[4].bst( 0., 0., betaZ);
2056 0 : pTH = pAbs * sin(theta);
2057 :
2058 : // Done.
2059 0 : return true;
2060 0 : }
2061 :
2062 : //--------------------------------------------------------------------------
2063 :
2064 : // Special choice of m3 and m4 when mHatMax push them off mass shell.
2065 : // Vary x in expression m3 + m4 = mHatMax - x * (Gamma3 + Gamma4).
2066 : // For each x try to put either 3 or 4 as close to mass shell as possible.
2067 : // Maximize BW_3 * BW_4 * beta_34, where latter approximate phase space.
2068 :
2069 : bool PhaseSpace2to2tauyz::constrainedM3M4() {
2070 :
2071 : // Initial values.
2072 : bool foundNonZero = false;
2073 : double wtMassMax = 0.;
2074 : double m3WtMax = 0.;
2075 : double m4WtMax = 0.;
2076 0 : double xMax = (mHatMax - mLower[3] - mLower[4]) / (mWidth[3] + mWidth[4]);
2077 0 : double xStep = THRESHOLDSTEP * min(1., xMax);
2078 : double xNow = 0.;
2079 : double wtMassXbin, wtMassMaxOld, m34, mT34Min, wtMassNow,
2080 : wtBW3Now, wtBW4Now, beta34Now;
2081 :
2082 : // Step through increasing x values.
2083 0 : do {
2084 0 : xNow += xStep;
2085 : wtMassXbin = 0.;
2086 : wtMassMaxOld = wtMassMax;
2087 0 : m34 = mHatMax - xNow * (mWidth[3] + mWidth[4]);
2088 :
2089 : // Study point where m3 as close as possible to on-shell.
2090 0 : m3 = min( mUpper[3], m34 - mLower[4]);
2091 0 : if (m3 > mPeak[3]) m3 = max( mLower[3], mPeak[3]);
2092 0 : m4 = m34 - m3;
2093 0 : if (m4 < mLower[4]) {m4 = mLower[4]; m3 = m34 - m4;}
2094 :
2095 : // Check that inside phase space limit set by pTmin.
2096 0 : mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin);
2097 0 : if (mT34Min < mHatMax) {
2098 :
2099 : // Breit-Wigners and beta factor give total weight.
2100 : wtMassNow = 0.;
2101 0 : if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4]
2102 0 : && m4 < mUpper[4]) {
2103 0 : wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2104 0 : wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2105 0 : beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2106 0 : - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2107 0 : wtMassNow = wtBW3Now * wtBW4Now * beta34Now;
2108 0 : }
2109 :
2110 : // Store new maximum, if any.
2111 0 : if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow;
2112 0 : if (wtMassNow > wtMassMax) {
2113 : foundNonZero = true;
2114 : wtMassMax = wtMassNow;
2115 0 : m3WtMax = m3;
2116 0 : m4WtMax = m4;
2117 0 : }
2118 : }
2119 :
2120 : // Study point where m4 as close as possible to on-shell.
2121 0 : m4 = min( mUpper[4], m34 - mLower[3]);
2122 0 : if (m4 > mPeak[4]) m4 = max( mLower[4], mPeak[4]);
2123 0 : m3 = m34 - m4;
2124 0 : if (m3 < mLower[3]) {m3 = mLower[3]; m4 = m34 - m3;}
2125 :
2126 : // Check that inside phase space limit set by pTmin.
2127 0 : mT34Min = sqrt(m3*m3 + pT2HatMin) + sqrt(m4*m4 + pT2HatMin);
2128 0 : if (mT34Min < mHatMax) {
2129 :
2130 : // Breit-Wigners and beta factor give total weight.
2131 : wtMassNow = 0.;
2132 0 : if (m3 > mLower[3] && m3 < mUpper[3] && m4 > mLower[4]
2133 0 : && m4 < mUpper[4]) {
2134 0 : wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2135 0 : wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2136 0 : beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2137 0 : - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2138 0 : wtMassNow = wtBW3Now * wtBW4Now * beta34Now;
2139 0 : }
2140 :
2141 : // Store new maximum, if any.
2142 0 : if (wtMassNow > wtMassXbin) wtMassXbin = wtMassNow;
2143 0 : if (wtMassNow > wtMassMax) {
2144 : foundNonZero = true;
2145 : wtMassMax = wtMassNow;
2146 0 : m3WtMax = m3;
2147 0 : m4WtMax = m4;
2148 0 : }
2149 : }
2150 :
2151 : // Continue stepping if increasing trend and more x range available.
2152 0 : } while ( (!foundNonZero || wtMassXbin > wtMassMaxOld)
2153 0 : && xNow < xMax - xStep);
2154 :
2155 : // Restore best values for subsequent maximization. Return.
2156 0 : m3 = m3WtMax;
2157 0 : m4 = m4WtMax;
2158 0 : return foundNonZero;
2159 :
2160 0 : }
2161 :
2162 : //--------------------------------------------------------------------------
2163 :
2164 : // Special choice of m3 when mHatMax pushes it off mass shell.
2165 : // Vary x in expression m3 = mHatMax - m4 - x * Gamma3.
2166 : // Maximize BW_3 * beta_34, where latter approximate phase space.
2167 :
2168 : bool PhaseSpace2to2tauyz::constrainedM3() {
2169 :
2170 : // Initial values.
2171 : bool foundNonZero = false;
2172 : double wtMassMax = 0.;
2173 : double m3WtMax = 0.;
2174 0 : double mT4Min = sqrt(m4*m4 + pT2HatMin);
2175 0 : double xMax = (mHatMax - mLower[3] - m4) / mWidth[3];
2176 0 : double xStep = THRESHOLDSTEP * min(1., xMax);
2177 : double xNow = 0.;
2178 : double wtMassNow, mT34Min, wtBW3Now, beta34Now;
2179 :
2180 : // Step through increasing x values; gives m3 unambiguously.
2181 0 : do {
2182 0 : xNow += xStep;
2183 : wtMassNow = 0.;
2184 0 : m3 = mHatMax - m4 - xNow * mWidth[3];
2185 :
2186 : // Check that inside phase space limit set by pTmin.
2187 0 : mT34Min = sqrt(m3*m3 + pT2HatMin) + mT4Min;
2188 0 : if (mT34Min < mHatMax) {
2189 :
2190 : // Breit-Wigner and beta factor give total weight.
2191 0 : wtBW3Now = mw[3] / ( pow2(m3*m3 - sPeak[3]) + pow2(mw[3]) );
2192 0 : beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2193 0 : - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2194 0 : wtMassNow = wtBW3Now * beta34Now;
2195 :
2196 : // Store new maximum, if any.
2197 0 : if (wtMassNow > wtMassMax) {
2198 : foundNonZero = true;
2199 : wtMassMax = wtMassNow;
2200 0 : m3WtMax = m3;
2201 0 : }
2202 : }
2203 :
2204 : // Continue stepping if increasing trend and more x range available.
2205 0 : } while ( (!foundNonZero || wtMassNow > wtMassMax)
2206 0 : && xNow < xMax - xStep);
2207 :
2208 : // Restore best value for subsequent maximization. Return.
2209 0 : m3 = m3WtMax;
2210 0 : return foundNonZero;
2211 :
2212 0 : }
2213 :
2214 : //--------------------------------------------------------------------------
2215 :
2216 : // Special choice of m4 when mHatMax pushes it off mass shell.
2217 : // Vary x in expression m4 = mHatMax - m3 - x * Gamma4.
2218 : // Maximize BW_4 * beta_34, where latter approximate phase space.
2219 :
2220 : bool PhaseSpace2to2tauyz::constrainedM4() {
2221 :
2222 : // Initial values.
2223 : bool foundNonZero = false;
2224 : double wtMassMax = 0.;
2225 : double m4WtMax = 0.;
2226 0 : double mT3Min = sqrt(m3*m3 + pT2HatMin);
2227 0 : double xMax = (mHatMax - mLower[4] - m3) / mWidth[4];
2228 0 : double xStep = THRESHOLDSTEP * min(1., xMax);
2229 : double xNow = 0.;
2230 : double wtMassNow, mT34Min, wtBW4Now, beta34Now;
2231 :
2232 : // Step through increasing x values; gives m4 unambiguously.
2233 0 : do {
2234 0 : xNow += xStep;
2235 : wtMassNow = 0.;
2236 0 : m4 = mHatMax - m3 - xNow * mWidth[4];
2237 :
2238 : // Check that inside phase space limit set by pTmin.
2239 0 : mT34Min = mT3Min + sqrt(m4*m4 + pT2HatMin);
2240 0 : if (mT34Min < mHatMax) {
2241 :
2242 : // Breit-Wigner and beta factor give total weight.
2243 0 : wtBW4Now = mw[4] / ( pow2(m4*m4 - sPeak[4]) + pow2(mw[4]) );
2244 0 : beta34Now = sqrt( pow2(mHatMax*mHatMax - m3*m3 - m4*m4)
2245 0 : - pow2(2. * m3 * m4) ) / (mHatMax*mHatMax);
2246 0 : wtMassNow = wtBW4Now * beta34Now;
2247 :
2248 : // Store new maximum, if any.
2249 0 : if (wtMassNow > wtMassMax) {
2250 : foundNonZero = true;
2251 : wtMassMax = wtMassNow;
2252 0 : m4WtMax = m4;
2253 0 : }
2254 : }
2255 :
2256 : // Continue stepping if increasing trend and more x range available.
2257 0 : } while ( (!foundNonZero || wtMassNow > wtMassMax)
2258 0 : && xNow < xMax - xStep);
2259 :
2260 : // Restore best value for subsequent maximization.
2261 0 : m4 = m4WtMax;
2262 0 : return foundNonZero;
2263 :
2264 0 : }
2265 :
2266 : //==========================================================================
2267 :
2268 : // PhaseSpace2to2elastic class.
2269 : // 2 -> 2 kinematics set up for elastic scattering.
2270 :
2271 : //--------------------------------------------------------------------------
2272 :
2273 : // Constants: could be changed here if desired, but normally should not.
2274 : // These are of technical nature, as described for each.
2275 :
2276 : // Maximum positive/negative argument for exponentiation.
2277 : const double PhaseSpace2to2elastic::EXPMAX = 50.;
2278 :
2279 : // Conversion coefficients = 1/(16pi) * (mb <-> GeV^2).
2280 : const double PhaseSpace2to2elastic::CONVERTEL = 0.0510925;
2281 :
2282 : //--------------------------------------------------------------------------
2283 :
2284 : // Form of phase space sampling already fixed, so no optimization.
2285 : // However, need to read out relevant parameters from SigmaTotal.
2286 :
2287 : bool PhaseSpace2to2elastic::setupSampling() {
2288 :
2289 : // Find maximum = value of cross section.
2290 0 : sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2291 0 : sigmaMx = sigmaNw;
2292 :
2293 : // Squared and outgoing masses of particles.
2294 0 : s1 = mA * mA;
2295 0 : s2 = mB * mB;
2296 0 : m3 = mA;
2297 0 : m4 = mB;
2298 :
2299 : // Elastic slope.
2300 0 : bSlope = sigmaTotPtr->bSlopeEl();
2301 :
2302 : // Determine maximum possible t range.
2303 0 : lambda12S = pow2(s - s1 - s2) - 4. * s1 * s2 ;
2304 0 : tLow = - lambda12S / s;
2305 0 : tUpp = 0.;
2306 :
2307 : // Production model with Coulomb corrections need more parameters.
2308 0 : useCoulomb = settingsPtr->flag("SigmaTotal:setOwn")
2309 0 : && settingsPtr->flag("SigmaElastic:setOwn");
2310 0 : if (useCoulomb) {
2311 0 : sigmaTot = sigmaTotPtr->sigmaTot();
2312 0 : rho = settingsPtr->parm("SigmaElastic:rho");
2313 0 : lambda = settingsPtr->parm("SigmaElastic:lambda");
2314 0 : tAbsMin = settingsPtr->parm("SigmaElastic:tAbsMin");
2315 0 : phaseCst = settingsPtr->parm("SigmaElastic:phaseConst");
2316 0 : alphaEM0 = settingsPtr->parm("StandardModel:alphaEM0");
2317 :
2318 : // Relative rate of nuclear and Coulombic parts in trials.
2319 0 : tUpp = -tAbsMin;
2320 0 : sigmaNuc = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho) / bSlope
2321 0 : * exp(-bSlope * tAbsMin);
2322 0 : sigmaCou = (useCoulomb) ?
2323 0 : pow2(alphaEM0) / (4. * CONVERTEL * tAbsMin) : 0.;
2324 0 : signCou = (idA == idB) ? 1. : -1.;
2325 :
2326 : // Dummy values.
2327 0 : } else {
2328 0 : sigmaNuc = sigmaNw;
2329 0 : sigmaCou = 0.;
2330 : }
2331 :
2332 : // Calculate coefficient of generation.
2333 0 : tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2334 :
2335 0 : return true;
2336 :
2337 0 : }
2338 :
2339 : //--------------------------------------------------------------------------
2340 :
2341 : // Select a trial kinematics phase space point. Perform full
2342 : // Monte Carlo acceptance/rejection at this stage.
2343 :
2344 : bool PhaseSpace2to2elastic::trialKin( bool, bool ) {
2345 :
2346 : // Allow for possibility that energy varies from event to event.
2347 0 : if (doEnergySpread) {
2348 0 : eCM = infoPtr->eCM();
2349 0 : s = eCM * eCM;
2350 0 : lambda12S = pow2(s - s1 - s2) - 4. * s1 * s2 ;
2351 0 : tLow = - lambda12S / s;
2352 0 : tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2353 0 : }
2354 :
2355 : // Select t according to exp(bSlope*t).
2356 0 : if (!useCoulomb || sigmaNuc > rndmPtr->flat() * (sigmaNuc + sigmaCou))
2357 0 : tH = tUpp + log(1. + tAux * rndmPtr->flat()) / bSlope;
2358 :
2359 : // Select t according to 1/t^2.
2360 0 : else tH = tLow * tUpp / (tUpp + rndmPtr->flat() * (tLow - tUpp));
2361 :
2362 : // Correction factor for ratio full/simulated.
2363 0 : if (useCoulomb) {
2364 0 : double sigmaN = CONVERTEL * pow2(sigmaTot) * (1. + rho*rho)
2365 0 : * exp(bSlope * tH);
2366 0 : double alpEM = couplingsPtr->alphaEM(-tH);
2367 0 : double sigmaC = pow2(alpEM) / (4. * CONVERTEL * tH*tH);
2368 0 : double sigmaGen = 2. * (sigmaN + sigmaC);
2369 0 : double form2 = pow4(lambda/(lambda - tH));
2370 0 : double phase = signCou * alpEM
2371 0 : * (-phaseCst - log(-0.5 * bSlope * tH));
2372 0 : double sigmaCor = sigmaN + pow2(form2) * sigmaC
2373 0 : - signCou * alpEM * sigmaTot * (form2 / (-tH))
2374 0 : * exp(0.5 * bSlope * tH) * (rho * cos(phase) + sin(phase));
2375 0 : sigmaNw = sigmaMx * sigmaCor / sigmaGen;
2376 0 : }
2377 :
2378 : // Careful reconstruction of scattering angle.
2379 0 : double tRat = s * tH / lambda12S;
2380 0 : double cosTheta = min(1., max(-1., 1. + 2. * tRat ) );
2381 0 : double sinTheta = 2. * sqrtpos( -tRat * (1. + tRat) );
2382 0 : theta = asin( min(1., sinTheta));
2383 0 : if (cosTheta < 0.) theta = M_PI - theta;
2384 :
2385 0 : return true;
2386 :
2387 0 : }
2388 :
2389 : //--------------------------------------------------------------------------
2390 :
2391 : // Construct the four-vector kinematics from the trial values.
2392 :
2393 : bool PhaseSpace2to2elastic::finalKin() {
2394 :
2395 : // Particle masses.
2396 0 : mH[1] = mA;
2397 0 : mH[2] = mB;
2398 0 : mH[3] = m3;
2399 0 : mH[4] = m4;
2400 :
2401 : // Incoming particles along beam axes.
2402 0 : pAbs = 0.5 * sqrtpos(lambda12S) / eCM;
2403 0 : pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2404 0 : pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2405 :
2406 : // Outgoing particles initially along beam axes.
2407 0 : pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2408 0 : pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2409 :
2410 : // Then rotate them
2411 0 : phi = 2. * M_PI * rndmPtr->flat();
2412 0 : pH[3].rot( theta, phi);
2413 0 : pH[4].rot( theta, phi);
2414 :
2415 : // Set some further info for completeness.
2416 0 : x1H = 1.;
2417 0 : x2H = 1.;
2418 0 : sH = s;
2419 0 : uH = 2. * (s1 + s2) - sH - tH;
2420 0 : mHat = eCM;
2421 0 : p2Abs = pAbs * pAbs;
2422 0 : betaZ = 0.;
2423 0 : pTH = pAbs * sin(theta);
2424 :
2425 : // Done.
2426 0 : return true;
2427 :
2428 : }
2429 :
2430 : //==========================================================================
2431 :
2432 : // PhaseSpace2to2diffractive class.
2433 : // 2 -> 2 kinematics set up for diffractive scattering.
2434 :
2435 : //--------------------------------------------------------------------------
2436 :
2437 : // Constants: could be changed here if desired, but normally should not.
2438 : // These are of technical nature, as described for each.
2439 :
2440 : // Number of tries to find acceptable (m^2, t) set.
2441 : const int PhaseSpace2to2diffractive::NTRY = 500;
2442 :
2443 : // Maximum positive/negative argument for exponentiation.
2444 : const double PhaseSpace2to2diffractive::EXPMAX = 50.;
2445 :
2446 : // Safety margin so sum of diffractive masses not too close to eCM.
2447 : const double PhaseSpace2to2diffractive::DIFFMASSMARGIN = 0.2;
2448 :
2449 : //--------------------------------------------------------------------------
2450 :
2451 : // Form of phase space sampling already fixed, so no optimization.
2452 : // However, need to read out relevant parameters from SigmaTotal.
2453 :
2454 : bool PhaseSpace2to2diffractive::setupSampling() {
2455 :
2456 : // Pomeron flux parametrization, and parameters of some options.
2457 0 : PomFlux = settingsPtr->mode("Diffraction:PomFlux");
2458 0 : epsilonPF = settingsPtr->parm("Diffraction:PomFluxEpsilon");
2459 0 : alphaPrimePF = settingsPtr->parm("Diffraction:PomFluxAlphaPrime");
2460 :
2461 : // Find maximum = value of cross section.
2462 0 : sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2463 0 : sigmaMx = sigmaNw;
2464 :
2465 : // Masses of particles and minimal masses of diffractive states.
2466 0 : m3ElDiff = (isDiffA) ? sigmaTotPtr->mMinXB() : mA;
2467 0 : m4ElDiff = (isDiffB) ? sigmaTotPtr->mMinAX() : mB;
2468 0 : s1 = mA * mA;
2469 0 : s2 = mB * mB;
2470 0 : s3 = pow2( m3ElDiff);
2471 0 : s4 = pow2( m4ElDiff);
2472 :
2473 : // Determine maximum possible t range and coefficient of generation.
2474 0 : lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
2475 0 : lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2476 0 : double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2477 0 : double tempB = lambda12 * lambda34 / s;
2478 0 : double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2479 0 : * (s1 * s4 - s2 * s3) / s;
2480 0 : tLow = -0.5 * (tempA + tempB);
2481 0 : tUpp = tempC / tLow;
2482 :
2483 : // Default for all parametrization-specific parameters.
2484 0 : cRes = sResXB = sResAX = sProton = bMin = bSlope = bSlope1 = bSlope2
2485 0 : = probSlope1 = xIntPF = xtCorPF = mp24DL = coefDL = tAux
2486 0 : = tAux1 = tAux2 = 0.;
2487 :
2488 : // Schuler&Sjostrand: parameters of low-mass-resonance enhancement.
2489 0 : if (PomFlux == 1) {
2490 0 : cRes = sigmaTotPtr->cRes();
2491 0 : sResXB = pow2( sigmaTotPtr->mResXB());
2492 0 : sResAX = pow2( sigmaTotPtr->mResAX());
2493 0 : sProton = sigmaTotPtr->sProton();
2494 :
2495 : // Schuler&Sjostrand: lower limit diffractive slope.
2496 0 : if (!isDiffB) bMin = sigmaTotPtr->bMinSlopeXB();
2497 0 : else if (!isDiffA) bMin = sigmaTotPtr->bMinSlopeAX();
2498 0 : else bMin = sigmaTotPtr->bMinSlopeXX();
2499 0 : tAux = exp( max(-EXPMAX, bMin * (tLow - tUpp)) ) - 1.;
2500 :
2501 : // Bruni&Ingelman: relative weight of two diffractive slopes.
2502 0 : } else if (PomFlux == 2) {
2503 0 : bSlope1 = 8.0;
2504 0 : probSlope1 = 6.38 * ( exp(max(-EXPMAX, bSlope1 * tUpp))
2505 0 : - exp(max(-EXPMAX, bSlope1 * tLow)) ) / bSlope1;
2506 0 : bSlope2 = 3.0;
2507 0 : double pS2 = 0.424 * ( exp(max(-EXPMAX, bSlope2 * tUpp))
2508 0 : - exp(max(-EXPMAX, bSlope2 * tLow)) ) / bSlope2;
2509 0 : probSlope1 /= probSlope1 + pS2;
2510 0 : tAux1 = exp( max(-EXPMAX, bSlope1 * (tLow - tUpp)) ) - 1.;
2511 0 : tAux2 = exp( max(-EXPMAX, bSlope2 * (tLow - tUpp)) ) - 1.;
2512 :
2513 : // Streng&Berger (RapGap): diffractive slope, power of mass spectrum.
2514 0 : } else if (PomFlux == 3) {
2515 0 : bSlope = 4.7;
2516 0 : double xPowPF = 1. - 2. * (1. + epsilonPF);
2517 0 : xIntPF = 2. * (1. + xPowPF);
2518 0 : xtCorPF = 2. * alphaPrimePF;
2519 0 : tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2520 :
2521 : // Donnachie&Landshoff (RapGap): power of mass spectrum.
2522 0 : } else if (PomFlux == 4) {
2523 0 : mp24DL = 4. * pow2(particleDataPtr->m0(2212));
2524 0 : double xPowPF = 1. - 2. * (1. + epsilonPF);
2525 0 : xIntPF = 2. * (1. + xPowPF);
2526 0 : xtCorPF = 2. * alphaPrimePF;
2527 : // Upper estimate of t dependence, for preliminary choice.
2528 0 : coefDL = 0.85;
2529 0 : tAux1 = 1. / pow3(1. - coefDL * tLow);
2530 0 : tAux2 = 1. / pow3(1. - coefDL * tUpp);
2531 :
2532 : // MBR model.
2533 0 : } else if (PomFlux == 5) {
2534 0 : eps = settingsPtr->parm("Diffraction:MBRepsilon");
2535 0 : alph = settingsPtr->parm("Diffraction:MBRalpha");
2536 0 : alph2 = alph * alph;
2537 0 : m2min = settingsPtr->parm("Diffraction:MBRm2Min");
2538 0 : dyminSD = settingsPtr->parm("Diffraction:MBRdyminSD");
2539 0 : dyminDD = settingsPtr->parm("Diffraction:MBRdyminDD");
2540 0 : dyminSigSD = settingsPtr->parm("Diffraction:MBRdyminSigSD");
2541 0 : dyminSigDD = settingsPtr->parm("Diffraction:MBRdyminSigDD");
2542 :
2543 : // Max f(dy) for Von Neumann method, from SigmaTot.
2544 0 : sdpmax= sigmaTotPtr->sdpMax();
2545 0 : ddpmax= sigmaTotPtr->ddpMax();
2546 :
2547 : // H1 Fit A/B.
2548 0 : } else if (PomFlux == 6 || PomFlux == 7) {
2549 0 : bSlope = 5.5;
2550 0 : epsilonPF = (PomFlux == 6) ? 0.1182 : 0.1110;
2551 0 : alphaPrimePF = 0.06;
2552 0 : double xPowPF = 1. - 2. * (1. + epsilonPF);
2553 0 : xIntPF = 2. * (1. + xPowPF);
2554 0 : xtCorPF = 2. * alphaPrimePF;
2555 0 : tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2556 0 : }
2557 :
2558 : // Done.
2559 0 : return true;
2560 :
2561 0 : }
2562 :
2563 : //--------------------------------------------------------------------------
2564 :
2565 : // Select a trial kinematics phase space point. Perform full
2566 : // Monte Carlo acceptance/rejection at this stage.
2567 :
2568 : bool PhaseSpace2to2diffractive::trialKin( bool, bool ) {
2569 :
2570 : // Allow for possibility that energy varies from event to event.
2571 0 : if (doEnergySpread) {
2572 0 : eCM = infoPtr->eCM();
2573 0 : s = eCM * eCM;
2574 0 : lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
2575 0 : lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2576 0 : double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2577 0 : double tempB = lambda12 * lambda34 / s;
2578 0 : double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2579 0 : * (s1 * s4 - s2 * s3) / s;
2580 0 : tLow = -0.5 * (tempA + tempB);
2581 0 : tUpp = tempC / tLow;
2582 0 : if (PomFlux == 1) {
2583 0 : tAux = exp( max(-EXPMAX, bMin * (tLow - tUpp)) ) - 1.;
2584 0 : } else if (PomFlux == 2) {
2585 0 : tAux1 = exp( max(-EXPMAX, bSlope1 * (tLow - tUpp)) ) - 1.;
2586 0 : tAux2 = exp( max(-EXPMAX, bSlope2 * (tLow - tUpp)) ) - 1.;
2587 0 : } else if (PomFlux == 3 || PomFlux == 6 || PomFlux == 7) {
2588 0 : tAux = exp( max(-EXPMAX, bSlope * (tLow - tUpp)) ) - 1.;
2589 0 : } else if (PomFlux == 4) {
2590 0 : tAux1 = 1. / pow3(1. - coefDL * tLow);
2591 0 : tAux2 = 1. / pow3(1. - coefDL * tUpp);
2592 0 : }
2593 0 : }
2594 :
2595 : // Loop over attempts to set up masses and t consistently.
2596 0 : for (int loop = 0; ; ++loop) {
2597 0 : if (loop == NTRY) {
2598 0 : infoPtr->errorMsg("Error in PhaseSpace2to2diffractive::trialKin: "
2599 : " quit after repeated tries");
2600 0 : return false;
2601 : }
2602 :
2603 : // Schuler and Sjostrand:
2604 0 : if (PomFlux == 1) {
2605 :
2606 : // Select diffractive mass(es) according to dm^2/m^2.
2607 0 : m3 = (isDiffA) ? m3ElDiff * pow( max(mA, eCM - m4ElDiff) / m3ElDiff,
2608 0 : rndmPtr->flat()) : m3ElDiff;
2609 0 : m4 = (isDiffB) ? m4ElDiff * pow( max(mB, eCM - m3ElDiff) / m4ElDiff,
2610 0 : rndmPtr->flat()) : m4ElDiff;
2611 0 : if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue;
2612 0 : s3 = m3 * m3;
2613 0 : s4 = m4 * m4;
2614 :
2615 : // Additional mass factors, including resonance enhancement.
2616 0 : if (isDiffA && !isDiffB) {
2617 0 : double facXB = (1. - s3 / s)
2618 0 : * (1. + cRes * sResXB / (sResXB + s3));
2619 0 : if (facXB < rndmPtr->flat() * (1. + cRes)) continue;
2620 0 : } else if (isDiffB && !isDiffA) {
2621 0 : double facAX = (1. - s4 / s)
2622 0 : * (1. + cRes * sResAX / (sResAX + s4));
2623 0 : if (facAX < rndmPtr->flat() * (1. + cRes)) continue;
2624 0 : } else {
2625 0 : double facXX = (1. - pow2(m3 + m4) / s)
2626 0 : * (s * sProton / (s * sProton + s3 * s4))
2627 0 : * (1. + cRes * sResXB / (sResXB + s3))
2628 0 : * (1. + cRes * sResAX / (sResAX + s4));
2629 0 : if (facXX < rndmPtr->flat() * pow2(1. + cRes)) continue;
2630 0 : }
2631 :
2632 : // Select t according to exp(bMin*t) and correct to right slope.
2633 0 : tH = tUpp + log(1. + tAux * rndmPtr->flat()) / bMin;
2634 0 : double bDiff = 0.;
2635 0 : if (isDiffA && !isDiffB) bDiff = sigmaTotPtr->bSlopeXB(s3) - bMin;
2636 0 : else if (!isDiffA) bDiff = sigmaTotPtr->bSlopeAX(s4) - bMin;
2637 0 : else bDiff = sigmaTotPtr->bSlopeXX(s3, s4) - bMin;
2638 0 : bDiff = max(0., bDiff);
2639 0 : if (exp( max(-EXPMAX, bDiff * (tH - tUpp)) ) < rndmPtr->flat()) continue;
2640 :
2641 : // Bruni and Ingelman:
2642 0 : } else if (PomFlux == 2) {
2643 :
2644 : // Select diffractive mass(es) according to dm^2/m^2.
2645 0 : m3 = (isDiffA) ? m3ElDiff * pow( max(mA, eCM - m4ElDiff) / m3ElDiff,
2646 0 : rndmPtr->flat()) : m3ElDiff;
2647 0 : m4 = (isDiffB) ? m4ElDiff * pow( max(mB, eCM - m3ElDiff) / m4ElDiff,
2648 0 : rndmPtr->flat()) : m4ElDiff;
2649 0 : if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue;
2650 0 : s3 = m3 * m3;
2651 0 : s4 = m4 * m4;
2652 :
2653 : // Select t according to exp(bSlope*t) with two possible slopes.
2654 0 : tH = (rndmPtr->flat() < probSlope1)
2655 0 : ? tUpp + log(1. + tAux1 * rndmPtr->flat()) / bSlope1
2656 0 : : tUpp + log(1. + tAux2 * rndmPtr->flat()) / bSlope2;
2657 :
2658 : // Streng and Berger et al. (RapGap) & H1 Fit A/B:
2659 0 : } else if (PomFlux == 3 || PomFlux == 6 || PomFlux == 7) {
2660 :
2661 : // Select diffractive mass(es) according to dm^2/(m^2)^(1 + 2 epsilon).
2662 0 : m3 = m3ElDiff;
2663 0 : m4 = m4ElDiff;
2664 0 : if (isDiffA) {
2665 0 : double s3MinPow = pow( m3ElDiff, xIntPF );
2666 0 : double s3MaxPow = pow( max(mA, eCM - m4ElDiff), xIntPF );
2667 0 : m3 = pow( s3MinPow + rndmPtr->flat() * (s3MaxPow - s3MinPow),
2668 0 : 1. / xIntPF );
2669 0 : }
2670 0 : if (isDiffB) {
2671 0 : double s4MinPow = pow( m4ElDiff, xIntPF );
2672 0 : double s4MaxPow = pow( max(mB, eCM - m3ElDiff), xIntPF );
2673 0 : m4 = pow( s4MinPow + rndmPtr->flat() * (s4MaxPow - s4MinPow),
2674 0 : 1. / xIntPF );
2675 0 : }
2676 0 : if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue;
2677 0 : s3 = m3 * m3;
2678 0 : s4 = m4 * m4;
2679 :
2680 : // Select t according to exponential and weigh by x_P^(2 alpha' |t|).
2681 0 : tH = tUpp + log(1. + tAux * rndmPtr->flat()) / bSlope;
2682 0 : if ( isDiffA && pow( s3 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() )
2683 : continue;
2684 0 : if ( isDiffB && pow( s4 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() )
2685 : continue;
2686 :
2687 : // Donnachie and Landshoff (RapGap):
2688 0 : } else if (PomFlux == 4) {
2689 :
2690 : // Select diffractive mass(es) according to dm^2/(m^2)^(1 + 2 epsilon).
2691 0 : m3 = m3ElDiff;
2692 0 : m4 = m4ElDiff;
2693 0 : if (isDiffA) {
2694 0 : double s3MinPow = pow( m3ElDiff, xIntPF );
2695 0 : double s3MaxPow = pow( max(mA, eCM - m4ElDiff), xIntPF );
2696 0 : m3 = pow( s3MinPow + rndmPtr->flat() * (s3MaxPow - s3MinPow),
2697 0 : 1. / xIntPF );
2698 0 : }
2699 0 : if (isDiffB) {
2700 0 : double s4MinPow = pow( m4ElDiff, xIntPF );
2701 0 : double s4MaxPow = pow( max(mB, eCM - m3ElDiff), xIntPF );
2702 0 : m4 = pow( s4MinPow + rndmPtr->flat() * (s4MaxPow - s4MinPow),
2703 0 : 1. / xIntPF );
2704 0 : }
2705 0 : if (m3 + m4 + DIFFMASSMARGIN >= eCM) continue;
2706 0 : s3 = m3 * m3;
2707 0 : s4 = m4 * m4;
2708 :
2709 : // Select t according to power and weigh by x_P^(2 alpha' |t|).
2710 0 : tH = - (1. / pow( tAux1 + rndmPtr->flat() * (tAux2 - tAux1), 1./3.)
2711 0 : - 1.) / coefDL;
2712 0 : double wDL = pow2( (mp24DL - 2.8 * tH) / (mp24DL - tH) )
2713 0 : / pow4( 1. - tH / 0.7);
2714 0 : double wMX = 1. / pow4( 1. - coefDL * tH);
2715 0 : if (wDL < rndmPtr->flat() * wMX) continue;
2716 0 : if ( isDiffA && pow( s3 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() )
2717 0 : continue;
2718 0 : if ( isDiffB && pow( s4 / s, xtCorPF * abs(tH) ) < rndmPtr->flat() )
2719 0 : continue;
2720 :
2721 : // MBR model:
2722 0 : } else if (PomFlux == 5) {
2723 0 : m3 = mA;
2724 0 : m4 = mB;
2725 : double xi, P, yRnd, dy;
2726 :
2727 : // MBR double diffractive.
2728 0 : if (isDiffA && isDiffB) {
2729 0 : dymin0 = 0.;
2730 0 : dymax = log(s/pow2(m2min));
2731 :
2732 : // Von Neumann method to generate dy, uses ddpmax from SigmaTot.
2733 0 : do {
2734 0 : dy = dymin0 + (dymax - dymin0) * rndmPtr->flat();
2735 0 : P = (dymax - dy) * exp(eps*dy) * ( exp(-2. * alph * dy * exp(-dy))
2736 0 : - exp(-2. * alph * dy * exp(dy)) ) / dy;
2737 : // Suppress smaller gap, smooth transition to non-diffractive.
2738 0 : P *= 0.5 * (1 + erf( ( dy - dyminDD) / dyminSigDD ) );
2739 0 : if (P > ddpmax) {
2740 0 : ostringstream osWarn;
2741 0 : osWarn << "ddpmax = " << scientific << setprecision(3)
2742 0 : << ddpmax << " " << P << " " << dy << endl;
2743 0 : infoPtr->errorMsg("Warning in PhaseSpace2to2diffractive::"
2744 0 : "trialKin for double diffraction:", osWarn.str());
2745 0 : }
2746 0 : yRnd = ddpmax * rndmPtr->flat();
2747 0 : } while (yRnd > P);
2748 :
2749 0 : double y0max = (dymax - dy)/2.;
2750 0 : double y0min = -y0max;
2751 0 : double y0 = y0min + (y0max - y0min) * rndmPtr->flat();
2752 0 : am1 = sqrt( eCM * exp( -y0 - dy/2. ) );
2753 0 : am2 = sqrt( eCM * exp( y0 - dy/2. ) );
2754 :
2755 : // Generate 4-momentum transfer, t from exp.
2756 0 : double b = 2. * alph * dy;
2757 0 : tUpp = -exp( -dy );
2758 0 : tLow = -exp( dy );
2759 0 : tAux = exp( b * (tLow - tUpp) ) - 1.;
2760 0 : t = tUpp + log(1. + tAux * rndmPtr->flat()) / b;
2761 0 : m3 = am1;
2762 0 : m4 = am2;
2763 0 : tH = t;
2764 :
2765 : // MBR single diffractive.
2766 0 : } else if (isDiffA || isDiffB) {
2767 0 : dymin0 = 0.;
2768 0 : dymax = log(s/m2min);
2769 :
2770 : // Von Neumann method to generate dy, uses sdpmax from SigmaTot.
2771 0 : do {
2772 0 : dy = dymin0 + (dymax - dymin0) * rndmPtr->flat();
2773 0 : P = exp(eps * dy) * ( (FFA1 / (FFB1 + 2. * alph * dy) )
2774 0 : + (FFA2 / (FFB2 + 2. * alph * dy) ) );
2775 : // Suppress smaller gap.
2776 0 : P *= 0.5 * (1. + erf( (dy - dyminSD) / dyminSigSD) );
2777 0 : if (P > sdpmax) {
2778 0 : ostringstream osWarn;
2779 0 : osWarn << "sdpmax = " << scientific << setprecision(3)
2780 0 : << sdpmax << " " << P << " " << dy << endl;
2781 0 : infoPtr->errorMsg("Warning in PhaseSpace2to2diffractive::"
2782 0 : "trialKin for single diffraction:", osWarn.str());
2783 0 : }
2784 0 : yRnd = sdpmax * rndmPtr->flat();
2785 0 : } while (yRnd > P);
2786 0 : xi = exp( -dy );
2787 0 : amx = sqrt( xi * s );
2788 :
2789 : // Generate 4-momentum transfer, t. First exponent, then FF*exp.
2790 0 : double tmin = -s1 * xi * xi / (1 - xi);
2791 0 : do {
2792 0 : t = tmin + log(1. - rndmPtr->flat());
2793 0 : double pFF = (4. * s1 - 2.8 * t) / ( (4. * s1 - t)
2794 0 : * pow2(1. - t / 0.71) );
2795 0 : P = pow2(pFF) * exp(2. * alph * dy * t);
2796 0 : yRnd = exp(t) * rndmPtr->flat();
2797 0 : } while (yRnd > P);
2798 0 : if(isDiffA) m3 = amx;
2799 0 : if(isDiffB) m4 = amx;
2800 0 : tH = t;
2801 0 : }
2802 :
2803 : // End of MBR model code.
2804 0 : s3 = m3 * m3;
2805 0 : s4 = m4 * m4;
2806 0 : }
2807 :
2808 : // Check whether m^2 and t choices are consistent.
2809 0 : lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2810 0 : double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2811 0 : double tempB = lambda12 * lambda34 / s;
2812 0 : double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2813 0 : * (s1 * s4 - s2 * s3) / s;
2814 0 : double tLowNow = -0.5 * (tempA + tempB);
2815 0 : double tUppNow = tempC / tLowNow;
2816 0 : if (tH < tLowNow || tH > tUppNow) continue;
2817 :
2818 : // Careful reconstruction of scattering angle.
2819 0 : double cosTheta = min(1., max(-1., (tempA + 2. * tH) / tempB));
2820 0 : double sinTheta = 2. * sqrtpos( -(tempC + tempA * tH + tH * tH) )
2821 0 : / tempB;
2822 0 : theta = asin( min(1., sinTheta));
2823 0 : if (cosTheta < 0.) theta = M_PI - theta;
2824 :
2825 : // Found acceptable kinematics, so no more looping. Done
2826 : break;
2827 0 : }
2828 0 : return true;
2829 :
2830 0 : }
2831 :
2832 : //--------------------------------------------------------------------------
2833 :
2834 : // Construct the four-vector kinematics from the trial values.
2835 :
2836 : bool PhaseSpace2to2diffractive::finalKin() {
2837 :
2838 : // Particle masses; incoming always on mass shell.
2839 0 : mH[1] = mA;
2840 0 : mH[2] = mB;
2841 0 : mH[3] = m3;
2842 0 : mH[4] = m4;
2843 :
2844 : // Incoming particles along beam axes.
2845 0 : pAbs = 0.5 * lambda12 / eCM;
2846 0 : pH[1] = Vec4( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
2847 0 : pH[2] = Vec4( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
2848 :
2849 : // Outgoing particles initially along beam axes.
2850 0 : pAbs = 0.5 * lambda34 / eCM;
2851 0 : pH[3] = Vec4( 0., 0., pAbs, 0.5 * (s + s3 - s4) / eCM);
2852 0 : pH[4] = Vec4( 0., 0., -pAbs, 0.5 * (s + s4 - s3) / eCM);
2853 :
2854 : // Then rotate them
2855 0 : phi = 2. * M_PI * rndmPtr->flat();
2856 0 : pH[3].rot( theta, phi);
2857 0 : pH[4].rot( theta, phi);
2858 :
2859 : // Set some further info for completeness (but Info can be for subprocess).
2860 0 : x1H = 1.;
2861 0 : x2H = 1.;
2862 0 : sH = s;
2863 0 : uH = s1 + s2 + s3 + s4 - sH - tH;
2864 0 : mHat = eCM;
2865 0 : p2Abs = pAbs * pAbs;
2866 0 : betaZ = 0.;
2867 0 : pTH = pAbs * sin(theta);
2868 :
2869 : // Done.
2870 0 : return true;
2871 :
2872 : }
2873 :
2874 : //==========================================================================
2875 :
2876 : // PhaseSpace2to3diffractive class.
2877 : // 2 -> 3 kinematics set up for central diffractive scattering.
2878 :
2879 : //--------------------------------------------------------------------------
2880 :
2881 : // Constants: could be changed here if desired, but normally should not.
2882 : // These are of technical nature, as described for each.
2883 :
2884 : // Number of tries to find acceptable (m^2, t1, t2) set.
2885 : const int PhaseSpace2to3diffractive::NTRY = 500;
2886 : const int PhaseSpace2to3diffractive::NINTEG2 = 40;
2887 :
2888 : // Maximum positive/negative argument for exponentiation.
2889 : const double PhaseSpace2to3diffractive::EXPMAX = 50.;
2890 :
2891 : // Minimal mass of central diffractive system.
2892 : const double PhaseSpace2to3diffractive::DIFFMASSMIN = 0.8;
2893 :
2894 : // Safety margin so sum of diffractive masses not too close to eCM.
2895 : const double PhaseSpace2to3diffractive::DIFFMASSMARGIN = 0.2;
2896 :
2897 : //--------------------------------------------------------------------------
2898 :
2899 : // Set up for phase space sampling.
2900 :
2901 : bool PhaseSpace2to3diffractive::setupSampling() {
2902 :
2903 : // Pomeron flux parametrization, and parameters of some options.
2904 0 : PomFlux = settingsPtr->mode("Diffraction:PomFlux");
2905 0 : epsilonPF = settingsPtr->parm("Diffraction:PomFluxEpsilon");
2906 0 : alphaPrimePF = settingsPtr->parm("Diffraction:PomFluxAlphaPrime");
2907 :
2908 : // Find maximum = value of cross section.
2909 0 : sigmaNw = sigmaProcessPtr->sigmaHatWrap();
2910 0 : sigmaMx = sigmaNw;
2911 :
2912 : // Squared masses of particles and minimal mass of diffractive states.
2913 0 : s1 = mA * mA;
2914 0 : s2 = mB * mB;
2915 0 : m5min = sigmaTotPtr->mMinAXB();
2916 0 : s5min = m5min * m5min;
2917 :
2918 : // Loop over two cases: s4 = (X + B)^2 and s3 = (A + X)^2.
2919 0 : for (int i = 0; i < 2; ++i) {
2920 0 : s3 = (i == 0) ? s1 : pow2(mA + m5min);
2921 0 : s4 = (i == 0) ? pow2(mB + m5min) : s2;
2922 :
2923 : // Determine maximum possible t range and coefficient of generation.
2924 0 : double lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
2925 0 : double lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
2926 0 : double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
2927 0 : double tempB = lambda12 * lambda34 / s;
2928 0 : double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
2929 0 : * (s1 * s4 - s2 * s3) / s;
2930 0 : tLow[i] = -0.5 * (tempA + tempB);
2931 0 : tUpp[i] = tempC / tLow[i];
2932 : }
2933 0 : s3 = s1;
2934 0 : s4 = s2;
2935 :
2936 : // Default for all parametrization-specific parameters.
2937 0 : bSlope1 = bSlope2 = bSlope = xIntPF = xIntInvPF = xtCorPF = mp24DL
2938 0 : = coefDL = 0.;
2939 0 : for (int i = 0; i < 2; ++i)
2940 0 : bMin[i] = tAux[i] = probSlope1[i] = tAux1[i] = tAux2[i] = 0.;
2941 :
2942 : // Schuler&Sjostrand: lower limit diffractive slope.
2943 0 : if (PomFlux == 1) {
2944 0 : bMin[0] = sigmaTotPtr->bMinSlopeAX();
2945 0 : tAux[0] = exp( max(-EXPMAX, bMin[0] * (tLow[0] - tUpp[0])) ) - 1.;
2946 0 : bMin[1] = sigmaTotPtr->bMinSlopeXB();
2947 0 : tAux[1] = exp( max(-EXPMAX, bMin[1] * (tLow[1] - tUpp[1])) ) - 1.;
2948 :
2949 : // Bruni&Ingelman: relative weight of two diffractive slopes.
2950 0 : } else if (PomFlux == 2) {
2951 0 : bSlope1 = 8.0;
2952 0 : bSlope2 = 3.0;
2953 0 : for (int i = 0; i < 2; ++i) {
2954 0 : probSlope1[i] = 6.38 * ( exp(max(-EXPMAX, bSlope1 * tUpp[i]))
2955 0 : - exp(max(-EXPMAX, bSlope1 * tLow[i])) ) / bSlope1;
2956 0 : double pS2 = 0.424 * ( exp(max(-EXPMAX, bSlope2 * tUpp[i]))
2957 0 : - exp(max(-EXPMAX, bSlope2 * tLow[i])) ) / bSlope2;
2958 0 : probSlope1[i] /= probSlope1[i] + pS2;
2959 0 : tAux1[i] = exp( max(-EXPMAX, bSlope1 * (tLow[i] - tUpp[i])) ) - 1.;
2960 0 : tAux2[i] = exp( max(-EXPMAX, bSlope2 * (tLow[i] - tUpp[i])) ) - 1.;
2961 : }
2962 :
2963 : // Streng&Berger (RapGap): diffractive slope, power of mass spectrum.
2964 0 : } else if (PomFlux == 3) {
2965 0 : bSlope = 4.7;
2966 0 : double xPowPF = 1. - 2. * (1. + epsilonPF);
2967 0 : xIntPF = 1. + xPowPF;
2968 0 : xIntInvPF = 1. / xIntPF;
2969 0 : xtCorPF = 2. * alphaPrimePF;
2970 0 : tAux[0] = exp( max(-EXPMAX, bSlope * (tLow[0] - tUpp[0])) ) - 1.;
2971 0 : tAux[1] = exp( max(-EXPMAX, bSlope * (tLow[1] - tUpp[1])) ) - 1.;
2972 :
2973 : // Donnachie&Landshoff (RapGap): power of mass spectrum.
2974 0 : } else if (PomFlux == 4) {
2975 0 : mp24DL = 4. * pow2(particleDataPtr->m0(2212));
2976 0 : double xPowPF = 1. - 2. * (1. + epsilonPF);
2977 0 : xIntPF = 1. + xPowPF;
2978 0 : xIntInvPF = 1. / xIntPF;
2979 0 : xtCorPF = 2. * alphaPrimePF;
2980 : // Upper estimate of t dependence, for preliminary choice.
2981 0 : coefDL = 0.85;
2982 0 : tAux1[0] = 1. / pow3(1. - coefDL * tLow[0]);
2983 0 : tAux2[0] = 1. / pow3(1. - coefDL * tUpp[0]);
2984 0 : tAux1[1] = 1. / pow3(1. - coefDL * tLow[1]);
2985 0 : tAux2[1] = 1. / pow3(1. - coefDL * tUpp[1]);
2986 :
2987 : // Setup for the MBR model.
2988 0 : } else if (PomFlux == 5) {
2989 0 : epsMBR = settingsPtr->parm("Diffraction:MBRepsilon");
2990 0 : alphMBR = settingsPtr->parm("Diffraction:MBRalpha");
2991 0 : m2minMBR = settingsPtr->parm("Diffraction:MBRm2Min");
2992 0 : dyminMBR = settingsPtr->parm("Diffraction:MBRdyminCD");
2993 0 : dyminSigMBR = settingsPtr->parm("Diffraction:MBRdyminSigCD");
2994 0 : dyminInvMBR = sqrt(2.) / dyminSigMBR;
2995 : // Max f(dy) for Von Neumann method, dpepmax from SigmaTot.
2996 0 : dpepmax = sigmaTotPtr->dpepMax();
2997 :
2998 : // H1 Fit A/B.
2999 0 : } else if (PomFlux == 6 || PomFlux == 7) {
3000 0 : bSlope = 5.5;
3001 0 : epsilonPF = (PomFlux == 6) ? 0.1182 : 0.1110;
3002 0 : alphaPrimePF = 0.06;
3003 0 : double xPowPF = 1. - 2. * (1. + epsilonPF);
3004 0 : xIntPF = 1. + xPowPF;
3005 0 : xIntInvPF = 1. / xIntPF;
3006 0 : xtCorPF = 2. * alphaPrimePF;
3007 0 : tAux[0] = exp( max(-EXPMAX, bSlope * (tLow[0] - tUpp[0])) ) - 1.;
3008 0 : tAux[1] = exp( max(-EXPMAX, bSlope * (tLow[1] - tUpp[1])) ) - 1.;
3009 0 : }
3010 :
3011 : // Done.
3012 0 : return true;
3013 :
3014 0 : }
3015 :
3016 : //--------------------------------------------------------------------------
3017 :
3018 : // Select a trial kinematics phase space point. Perform full
3019 : // Monte Carlo acceptance/rejection at this stage.
3020 :
3021 : bool PhaseSpace2to3diffractive::trialKin( bool, bool ) {
3022 :
3023 : // Allow for possibility that energy varies from event to event.
3024 0 : if (doEnergySpread) {
3025 0 : eCM = infoPtr->eCM();
3026 0 : s = eCM * eCM;
3027 0 : for (int i = 0; i < 2; ++i) {
3028 0 : s3 = (i == 0) ? s1 : pow2(mA + m5min);
3029 0 : s4 = (i == 0) ? pow2(mB + m5min) : s2;
3030 0 : double lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
3031 0 : double lambda34 = sqrtpos( pow2( s - s3 - s4) - 4. * s3 * s4 );
3032 0 : double tempA = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
3033 0 : double tempB = lambda12 * lambda34 / s;
3034 0 : double tempC = (s3 - s1) * (s4 - s2) + (s1 + s4 - s2 - s3)
3035 0 : * (s1 * s4 - s2 * s3) / s;
3036 0 : tLow[i] = -0.5 * (tempA + tempB);
3037 0 : tUpp[i] = tempC / tLow[i];
3038 : }
3039 0 : s3 = s1;
3040 0 : s4 = s2;
3041 0 : if (PomFlux == 1) {
3042 0 : tAux[0] = exp( max(-EXPMAX, bMin[0] * (tLow[0] - tUpp[0])) ) - 1.;
3043 0 : tAux[1] = exp( max(-EXPMAX, bMin[1] * (tLow[1] - tUpp[1])) ) - 1.;
3044 0 : } else if (PomFlux == 2) {
3045 0 : for (int i = 0; i < 2; ++i) {
3046 0 : tAux1[i] = exp( max(-EXPMAX, bSlope1 * (tLow[i] - tUpp[i])) ) - 1.;
3047 0 : tAux2[i] = exp( max(-EXPMAX, bSlope2 * (tLow[i] - tUpp[i])) ) - 1.;
3048 : }
3049 0 : } else if (PomFlux == 3 || PomFlux == 6 || PomFlux == 7) {
3050 0 : tAux[0] = exp( max(-EXPMAX, bSlope * (tLow[0] - tUpp[0])) ) - 1.;
3051 0 : tAux[1] = exp( max(-EXPMAX, bSlope * (tLow[1] - tUpp[1])) ) - 1.;
3052 0 : } else if (PomFlux == 4) {
3053 0 : tAux1[0] = 1. / pow3(1. - coefDL * tLow[0]);
3054 0 : tAux2[0] = 1. / pow3(1. - coefDL * tUpp[0]);
3055 0 : tAux1[1] = 1. / pow3(1. - coefDL * tLow[1]);
3056 0 : tAux2[1] = 1. / pow3(1. - coefDL * tUpp[1]);
3057 0 : }
3058 : }
3059 :
3060 : // Trivial kinematics of incoming hadrons.
3061 0 : double lambda12 = sqrtpos( pow2( s - s1 - s2) - 4. * s1 * s2 );
3062 0 : pAbs = 0.5 * lambda12 / eCM;
3063 0 : p1.p( 0., 0., pAbs, 0.5 * (s + s1 - s2) / eCM);
3064 0 : p2.p( 0., 0., -pAbs, 0.5 * (s + s2 - s1) / eCM);
3065 :
3066 : // Loop over attempts to set up mass, t1, t2 consistently.
3067 0 : for (int loop = 0; ; ++loop) {
3068 0 : if (loop == NTRY) {
3069 0 : infoPtr->errorMsg("Error in PhaseSpace2to3diffractive::trialKin: "
3070 : " quit after repeated tries");
3071 0 : return false;
3072 : }
3073 0 : double xi1 = 0.;
3074 0 : double xi2 = 0.;
3075 0 : double tVal[2] = { 0., 0.};
3076 :
3077 : // Schuler and Sjostrand:
3078 0 : if (PomFlux == 1) {
3079 :
3080 : // Select mass according to dxi_1/xi_1 * dxi_2/xi_2 * (1 - m^2/s).
3081 : do {
3082 0 : xi1 = pow( s5min / s, rndmPtr->flat());
3083 0 : xi2 = pow( s5min / s, rndmPtr->flat());
3084 0 : s5 = xi1 * xi2 * s;
3085 0 : } while (s5 < s5min || xi1 * xi2 > rndmPtr->flat());
3086 0 : if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue;
3087 :
3088 : // Select t according to exp(bMin*t) and correct to right slope.
3089 : bool tryAgain = false;
3090 0 : for (int i = 0; i < 2; ++i) {
3091 0 : tVal[i] = tUpp[i] + log(1. + tAux[i] * rndmPtr->flat()) / bMin[i];
3092 0 : double bDiff = (i == 0) ? sigmaTotPtr->bSlopeAX(s2 + xi1 * s)
3093 0 : : sigmaTotPtr->bSlopeXB(s1 + xi2 * s);
3094 0 : bDiff = max(0., bDiff - bMin[i]);
3095 0 : if (exp( max(-EXPMAX, bDiff * (tVal[i] - tUpp[i])) )
3096 0 : < rndmPtr->flat()) tryAgain = true;
3097 : }
3098 0 : if (tryAgain) continue;
3099 :
3100 : // Bruni and Ingelman:
3101 0 : } else if (PomFlux == 2) {
3102 :
3103 : // Select mass according to dxi_1/xi_1 * dxi_2/xi_2.
3104 : do {
3105 0 : xi1 = pow( s5min / s, rndmPtr->flat());
3106 0 : xi2 = pow( s5min / s, rndmPtr->flat());
3107 0 : s5 = xi1 * xi2 * s;
3108 0 : } while (s5 < s5min);
3109 0 : if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue;
3110 :
3111 : // Select t according to exp(bSlope*t) with two possible slopes.
3112 0 : for (int i = 0; i < 2; ++i)
3113 0 : tVal[i] = (rndmPtr->flat() < probSlope1[i])
3114 0 : ? tUpp[i] + log(1. + tAux1[i] * rndmPtr->flat()) / bSlope1
3115 0 : : tUpp[i] + log(1. + tAux2[i] * rndmPtr->flat()) / bSlope2;
3116 :
3117 : // Streng and Berger et al. (RapGap) and H1 Fit A/B:
3118 0 : } else if (PomFlux == 3 || PomFlux == 6 || PomFlux == 7) {
3119 :
3120 : // Select mass by dxi_1 * dxi_2 / (xi_1 * xi_2)^(1 + 2 epsilon).
3121 0 : double sMinPow = pow( s5min / s, xIntPF);
3122 0 : do {
3123 0 : xi1 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF );
3124 0 : xi2 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF );
3125 0 : s5 = xi1 * xi2 * s;
3126 0 : } while (s5 < s5min);
3127 0 : if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue;
3128 :
3129 : // Select t according to exponential and weigh by x_P^(2 alpha' |t|).
3130 : bool tryAgain = false;
3131 0 : for (int i = 0; i < 2; ++i) {
3132 0 : tVal[i] = tUpp[i] + log(1. + tAux[i] * rndmPtr->flat()) / bSlope;
3133 0 : double xi = (i == 0) ? xi1 : xi2;
3134 0 : if ( pow( xi, xtCorPF * abs(tVal[i]) ) < rndmPtr->flat() )
3135 0 : tryAgain = true;
3136 : }
3137 0 : if (tryAgain) continue;
3138 :
3139 : // Donnachie and Landshoff (RapGap):
3140 0 : } else if (PomFlux == 4) {
3141 :
3142 : // Select mass by dxi_1 * dxi_2 / (xi_1 * xi_2)^(1 + 2 epsilon).
3143 0 : double sMinPow = pow( s5min / s, xIntPF);
3144 0 : do {
3145 0 : xi1 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF );
3146 0 : xi2 = pow( sMinPow + rndmPtr->flat() * (1. - sMinPow), xIntInvPF );
3147 0 : s5 = xi1 * xi2 * s;
3148 0 : } while (s5 < s5min);
3149 0 : if (mA + mB + sqrt(s5) + DIFFMASSMARGIN >= eCM) continue;
3150 :
3151 : // Select t according to power and weigh by x_P^(2 alpha' |t|).
3152 : bool tryAgain = false;
3153 0 : for (int i = 0; i < 2; ++i) {
3154 0 : tVal[i] = - (1. / pow( tAux1[i] + rndmPtr->flat()
3155 0 : * (tAux2[i] - tAux1[i]), 1./3.) - 1.) / coefDL;
3156 0 : double wDL = pow2( (mp24DL - 2.8 * tVal[i]) / (mp24DL - tVal[i]) )
3157 0 : / pow4( 1. - tVal[i] / 0.7);
3158 0 : double wMX = 1. / pow4( 1. - coefDL * tVal[i]);
3159 0 : if (wDL < rndmPtr->flat() * wMX) tryAgain = true;
3160 0 : double xi = (i == 0) ? xi1 : xi2;
3161 0 : if ( pow( xi, xtCorPF * abs(tVal[i]) ) < rndmPtr->flat() )
3162 0 : tryAgain = true;
3163 : }
3164 0 : if (tryAgain) continue;
3165 :
3166 : // The MBR model (PomFlux == 5).
3167 0 : } else if (PomFlux == 5) {
3168 : double dymin0 = 0.;
3169 0 : double dymax = log(s/m2minMBR);
3170 : double f1, f2, step2, dy, yc, ycmin, ycmax, dy1, dy2,
3171 : P, P1, P2, yRnd, yRnd1, yRnd2;
3172 :
3173 : // Von Neumann method to generate dy, uses dpepmax from SigmaTot.
3174 0 : do {
3175 0 : dy = dymin0 + (dymax - dymin0) * rndmPtr->flat();
3176 : P = 0.;
3177 0 : step2 = (dy - dymin0) / NINTEG2;
3178 0 : for (int j = 0; j < NINTEG2 ; ++j) {
3179 0 : yc = -(dy - dymin0) / 2. + (j + 0.5) * step2;
3180 0 : dy1 = 0.5 * dy - yc;
3181 0 : dy2 = 0.5 * dy + yc;
3182 0 : f1 = exp(epsMBR * dy1) * ( (FFA1 / (FFB1 + 2. * alphMBR * dy1) )
3183 0 : + (FFA2 / (FFB2 + 2. * alphMBR * dy1) ) );
3184 0 : f2 = exp(epsMBR * dy2) * ( (FFA1 / (FFB1 + 2. * alphMBR * dy2) )
3185 0 : + (FFA2 / (FFB2 + 2. * alphMBR * dy2) ) );
3186 0 : f1 *= 0.5 * (1. + erf( (dy1 - 0.5 * dyminMBR) * dyminInvMBR ));
3187 0 : f2 *= 0.5 * (1. + erf( (dy2 - 0.5 * dyminMBR) * dyminInvMBR ));
3188 0 : P += f1 * f2 * step2;
3189 : }
3190 0 : if (P > dpepmax) {
3191 0 : ostringstream osWarn;
3192 0 : osWarn << "dpepmax = " << scientific << setprecision(3)
3193 0 : << dpepmax << " " << P << " " << dy << endl;
3194 0 : infoPtr->errorMsg("Warning in PhaseSpace2to2diffractive::"
3195 0 : "trialKin for central diffraction:", osWarn.str());
3196 0 : }
3197 0 : yRnd = dpepmax * rndmPtr->flat();
3198 :
3199 : // Generate dyc.
3200 0 : ycmax = (dy - dymin0) / 2.;
3201 0 : ycmin = -ycmax;
3202 0 : yc = ycmin + (ycmax - ycmin) * rndmPtr->flat();
3203 :
3204 : // xi1, xi2 from dy and dy0.
3205 0 : dy1 = 0.5 * dy + yc;
3206 0 : dy2 = 0.5 * dy - yc;
3207 0 : P1 = 0.5 * (1. + erf( (dy1 - 0.5 * dyminMBR) * dyminInvMBR ));
3208 0 : P2 = 0.5 * (1. + erf( (dy2 - 0.5 * dyminMBR) * dyminInvMBR ));
3209 0 : yRnd1 = rndmPtr->flat();
3210 0 : yRnd2 = rndmPtr->flat();
3211 0 : } while( !(yRnd < P && yRnd1 < P1 && yRnd2 < P2) );
3212 0 : xi1 = exp( -dy1 );
3213 0 : xi2 = exp( -dy2 );
3214 :
3215 : // Generate t1 at vertex1. First exponent, then FF*exp.
3216 0 : double tmin = -s1 * xi1 * xi1 / (1. - xi1);
3217 0 : do {
3218 0 : t1 = tmin + log(1. - rndmPtr->flat());
3219 0 : double pFF = (4. * s1 - 2.8 * t1) / ( (4. * s1 - t1)
3220 0 : * pow2(1. - t1 / 0.71));
3221 0 : P = pow2(pFF) * exp(2. * alphMBR * dy1 * t1);
3222 0 : yRnd = exp(t1) * rndmPtr->flat();
3223 0 : } while (yRnd > P);
3224 :
3225 : // Generate t2 at vertex2. First exponent, then FF*exp.
3226 0 : tmin = -s2 * xi2 * xi2 / (1. - xi2);
3227 0 : do {
3228 0 : t2 = tmin + log(1. - rndmPtr->flat());
3229 0 : double pFF = (4. * s2 - 2.8 * t2) / ((4. * s2 - t2)
3230 0 : * pow2(1. - t2 / 0.71));
3231 0 : P = pow2(pFF) * exp(2. * alphMBR * dy2 * t2);
3232 0 : yRnd = exp(t2) * rndmPtr->flat();
3233 0 : } while (yRnd > P);
3234 :
3235 0 : }
3236 :
3237 : // Checks and kinematics construction four first options.
3238 : double pz3 = 0.;
3239 : double pz4 = 0.;
3240 : double pT3 = 0.;
3241 : double pT4 = 0.;
3242 0 : if (PomFlux != 5) {
3243 :
3244 : // Check whether m^2 (i.e. xi) and t choices are consistent.
3245 : bool tryAgain = false;
3246 0 : for (int i = 0; i < 2; ++i) {
3247 0 : double sx1 = (i == 0) ? s1 : s2;
3248 0 : double sx2 = (i == 0) ? s2 : s1;
3249 : double sx3 = sx1;
3250 0 : double sx4 = (i == 0) ? s2 + xi1 * s : s1 + xi2 * s;
3251 0 : if (sqrt(sx3) + sqrt(sx4) + DIFFMASSMARGIN > eCM) tryAgain = true;
3252 0 : double lambda34 = sqrtpos( pow2( s - sx3 - sx4) - 4. * sx3 * sx4 );
3253 0 : double tempA = s - (sx1 + sx2 + sx3 + sx4)
3254 0 : + (sx1 - sx2) * (sx3 - sx4) / s;
3255 0 : double tempB = lambda12 * lambda34 / s;
3256 0 : double tempC = (sx3 - sx1) * (sx4 - sx2) + (sx1 + sx4 - sx2 - sx3)
3257 0 : * (sx1 * sx4 - sx2 * sx3) / s;
3258 0 : double tLowNow = -0.5 * (tempA + tempB);
3259 0 : double tUppNow = tempC / tLowNow;
3260 0 : if (tVal[i] < tLowNow || tVal[i] > tUppNow) tryAgain = true;
3261 0 : if (tryAgain) break;
3262 :
3263 : // Careful reconstruction of scattering angle.
3264 0 : double cosTheta = min(1., max(-1., (tempA + 2. * tVal[i]) / tempB));
3265 0 : double sinTheta = 2. * sqrtpos( -(tempC + tempA * tVal[i]
3266 0 : + tVal[i] * tVal[i]) ) / tempB;
3267 0 : theta = asin( min(1., sinTheta));
3268 0 : if (cosTheta < 0.) theta = M_PI - theta;
3269 0 : double pAbs34 = 0.5 * lambda34 / eCM;
3270 0 : if (i == 0) {
3271 0 : pz3 = pAbs34 * cos(theta);
3272 0 : pT3 = pAbs34 * sin(theta);
3273 0 : } else {
3274 0 : pz4 = -pAbs34 * cos(theta);
3275 0 : pT4 = pAbs34 * sin(theta);
3276 : }
3277 0 : }
3278 0 : if (tryAgain) continue;
3279 0 : t1 = tVal[0];
3280 0 : t2 = tVal[1];
3281 :
3282 : // Kinematics construction in the MBR model.
3283 0 : } else {
3284 0 : pz3 = pAbs * (1. - xi1);
3285 0 : pz4 = -pAbs * (1. - xi2);
3286 0 : pT3 = sqrt( (1. - xi1) * abs(t1) - s1 * pow2(xi1) );
3287 0 : pT4 = sqrt( (1. - xi2) * abs(t2) - s2 * pow2(xi2) );
3288 : }
3289 :
3290 : // Common final steps of kinematics.
3291 0 : double phi3 = 2. * M_PI * rndmPtr->flat();
3292 0 : double phi4 = 2. * M_PI * rndmPtr->flat();
3293 0 : p3.p( pT3 * cos(phi3), pT3 * sin(phi3), pz3,
3294 0 : sqrt(pz3 * pz3 + pT3 * pT3 + s1) );
3295 0 : p4.p( pT4 * cos(phi4), pT4 * sin(phi4), pz4,
3296 0 : sqrt(pz4 * pz4 + pT4 * pT4 + s2) );
3297 :
3298 : // Central dissociated system, from Pomeron-Pomeron 4 vectors.
3299 0 : p5 = (p1 - p3) + (p2 - p4);
3300 0 : mHat = p5.mCalc();
3301 :
3302 : // If acceptable diffractive mass then no more looping.
3303 0 : if (mHat > DIFFMASSMIN) break;
3304 0 : }
3305 0 : return true;
3306 :
3307 0 : }
3308 :
3309 : //--------------------------------------------------------------------------
3310 :
3311 : // Construct the four-vector kinematics from the trial values.
3312 :
3313 : bool PhaseSpace2to3diffractive::finalKin() {
3314 :
3315 : // Particle four-momenta and masses.
3316 0 : pH[1] = p1;
3317 0 : pH[2] = p2;
3318 0 : pH[3] = p3;
3319 0 : pH[4] = p4;
3320 0 : pH[5] = p5;
3321 0 : mH[1] = mA;
3322 0 : mH[2] = mB;
3323 0 : mH[3] = mA;
3324 0 : mH[4] = mB;
3325 0 : mH[5] = mHat;
3326 :
3327 : // Set some further info for completeness (but Info can be for subprocess).
3328 0 : x1H = 1.;
3329 0 : x2H = 1.;
3330 0 : sH = s;
3331 0 : tH = (p1 - p3).m2Calc();
3332 0 : uH = (p2 - p4).m2Calc();
3333 0 : mHat = eCM;
3334 0 : p2Abs = pAbs * pAbs;
3335 0 : betaZ = 0.;
3336 : // Store average pT of three final particles for documentation.
3337 0 : pTH = (p3.pT() + p4.pT() + p5.pT()) / 3.;
3338 :
3339 : // Done.
3340 0 : return true;
3341 :
3342 : }
3343 :
3344 : //==========================================================================
3345 :
3346 : // PhaseSpace2to3tauycyl class.
3347 : // 2 -> 3 kinematics for normal subprocesses.
3348 :
3349 : //--------------------------------------------------------------------------
3350 :
3351 : // Constants: could be changed here if desired, but normally should not.
3352 : // These are of technical nature, as described for each.
3353 :
3354 : // Number of Newton-Raphson iterations of kinematics when masses introduced.
3355 : const int PhaseSpace2to3tauycyl::NITERNR = 5;
3356 :
3357 : //--------------------------------------------------------------------------
3358 :
3359 : // Set up for fixed or Breit-Wigner mass selection.
3360 :
3361 : bool PhaseSpace2to3tauycyl::setupMasses() {
3362 :
3363 : // Treat Z0 as such or as gamma*/Z0
3364 0 : gmZmode = gmZmodeGlobal;
3365 0 : int gmZmodeProc = sigmaProcessPtr->gmZmode();
3366 0 : if (gmZmodeProc >= 0) gmZmode = gmZmodeProc;
3367 :
3368 : // Set sHat limits - based on global limits only.
3369 0 : mHatMin = mHatGlobalMin;
3370 0 : sHatMin = mHatMin*mHatMin;
3371 0 : mHatMax = eCM;
3372 0 : if (mHatGlobalMax > mHatGlobalMin) mHatMax = min( eCM, mHatGlobalMax);
3373 0 : sHatMax = mHatMax*mHatMax;
3374 :
3375 : // Masses and widths of resonances.
3376 0 : setupMass1(3);
3377 0 : setupMass1(4);
3378 0 : setupMass1(5);
3379 :
3380 : // Reduced mass range - do not make it as fancy as in two-body case.
3381 0 : if (useBW[3]) mUpper[3] -= (mPeak[4] + mPeak[5]);
3382 0 : if (useBW[4]) mUpper[4] -= (mPeak[3] + mPeak[5]);
3383 0 : if (useBW[5]) mUpper[5] -= (mPeak[3] + mPeak[4]);
3384 :
3385 : // If closed phase space then unallowed process.
3386 : bool physical = true;
3387 0 : if (useBW[3] && mUpper[3] < mLower[3] + MASSMARGIN) physical = false;
3388 0 : if (useBW[4] && mUpper[4] < mLower[4] + MASSMARGIN) physical = false;
3389 0 : if (useBW[5] && mUpper[5] < mLower[5] + MASSMARGIN) physical = false;
3390 0 : if (!useBW[3] && !useBW[4] && !useBW[5] && mHatMax < mPeak[3]
3391 0 : + mPeak[4] + mPeak[5] + MASSMARGIN) physical = false;
3392 0 : if (!physical) return false;
3393 :
3394 : // No extra pT precautions in massless limit - assumed fixed by ME's.
3395 0 : pTHatMin = pTHatGlobalMin;
3396 0 : pT2HatMin = pTHatMin * pTHatMin;
3397 0 : pTHatMax = pTHatGlobalMax;
3398 0 : pT2HatMax = pTHatMax * pTHatMax;
3399 :
3400 : // Prepare to select m3 by BW + flat + 1/s_3.
3401 0 : if (useBW[3]) {
3402 0 : double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
3403 0 : * mWidth[3] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
3404 0 : double distToThreshB = (mHatMax - mPeak[3] - mMin[4] - mMin[5])
3405 0 : / mWidth[3];
3406 0 : double distToThresh = min( distToThreshA, distToThreshB);
3407 0 : setupMass2(3, distToThresh);
3408 0 : }
3409 :
3410 : // Prepare to select m4 by BW + flat + 1/s_3.
3411 0 : if (useBW[4]) {
3412 0 : double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
3413 0 : * mWidth[4] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
3414 0 : double distToThreshB = (mHatMax - mPeak[4] - mMin[3] - mMin[5])
3415 0 : / mWidth[4];
3416 0 : double distToThresh = min( distToThreshA, distToThreshB);
3417 0 : setupMass2(4, distToThresh);
3418 0 : }
3419 :
3420 : // Prepare to select m5 by BW + flat + 1/s_3.
3421 0 : if (useBW[5]) {
3422 0 : double distToThreshA = (mHatMax - mPeak[3] - mPeak[4] - mPeak[5])
3423 0 : * mWidth[5] / (pow2(mWidth[3]) + pow2(mWidth[4]) + pow2(mWidth[5]));
3424 0 : double distToThreshB = (mHatMax - mPeak[5] - mMin[3] - mMin[4])
3425 0 : / mWidth[5];
3426 0 : double distToThresh = min( distToThreshA, distToThreshB);
3427 0 : setupMass2(5, distToThresh);
3428 0 : }
3429 :
3430 : // Initialization masses. For now give up when constrained phase space.
3431 0 : m3 = (useBW[3]) ? min(mPeak[3], mUpper[3]) : mPeak[3];
3432 0 : m4 = (useBW[4]) ? min(mPeak[4], mUpper[4]) : mPeak[4];
3433 0 : m5 = (useBW[5]) ? min(mPeak[5], mUpper[5]) : mPeak[5];
3434 0 : if (m3 + m4 + m5 + MASSMARGIN > mHatMax) physical = false;
3435 0 : s3 = m3*m3;
3436 0 : s4 = m4*m4;
3437 0 : s5 = m5*m5;
3438 :
3439 : // Correct selected mass-spectrum to running-width Breit-Wigner.
3440 : // Extra safety margin for maximum search.
3441 0 : wtBW = 1.;
3442 0 : if (useBW[3]) wtBW *= weightMass(3) * EXTRABWWTMAX;
3443 0 : if (useBW[4]) wtBW *= weightMass(4) * EXTRABWWTMAX;
3444 0 : if (useBW[5]) wtBW *= weightMass(5) * EXTRABWWTMAX;
3445 :
3446 : // Done.
3447 0 : return physical;
3448 :
3449 0 : }
3450 :
3451 : //--------------------------------------------------------------------------
3452 :
3453 : // Select Breit-Wigner-distributed or fixed masses.
3454 :
3455 : bool PhaseSpace2to3tauycyl::trialMasses() {
3456 :
3457 : // By default vanishing cross section.
3458 0 : sigmaNw = 0.;
3459 0 : wtBW = 1.;
3460 :
3461 : // Pick m3, m4 and m5 independently.
3462 0 : trialMass(3);
3463 0 : trialMass(4);
3464 0 : trialMass(5);
3465 :
3466 : // If outside phase space then reject event.
3467 0 : if (m3 + m4 + m5 + MASSMARGIN > mHatMax) return false;
3468 :
3469 : // Correct selected mass-spectrum to running-width Breit-Wigner.
3470 0 : if (useBW[3]) wtBW *= weightMass(3);
3471 0 : if (useBW[4]) wtBW *= weightMass(4);
3472 0 : if (useBW[5]) wtBW *= weightMass(5);
3473 :
3474 : // Done.
3475 0 : return true;
3476 :
3477 0 : }
3478 :
3479 : //--------------------------------------------------------------------------
3480 :
3481 : // Construct the four-vector kinematics from the trial values.
3482 :
3483 : bool PhaseSpace2to3tauycyl::finalKin() {
3484 :
3485 : // Assign masses to particles assumed massless in matrix elements.
3486 0 : int id3 = sigmaProcessPtr->id(3);
3487 0 : int id4 = sigmaProcessPtr->id(4);
3488 0 : int id5 = sigmaProcessPtr->id(5);
3489 0 : if (idMass[3] == 0) { m3 = particleDataPtr->m0(id3); s3 = m3*m3; }
3490 0 : if (idMass[4] == 0) { m4 = particleDataPtr->m0(id4); s4 = m4*m4; }
3491 0 : if (idMass[5] == 0) { m5 = particleDataPtr->m0(id5); s5 = m5*m5; }
3492 :
3493 : // Check that phase space still open after new mass assignment.
3494 0 : if (m3 + m4 + m5 + MASSMARGIN > mHat) {
3495 0 : infoPtr->errorMsg("Warning in PhaseSpace2to3tauycyl::finalKin: "
3496 : "failed after mass assignment");
3497 0 : return false;
3498 : }
3499 :
3500 : // Particle masses; incoming always on mass shell.
3501 0 : mH[1] = 0.;
3502 0 : mH[2] = 0.;
3503 0 : mH[3] = m3;
3504 0 : mH[4] = m4;
3505 0 : mH[5] = m5;
3506 :
3507 : // Incoming partons along beam axes.
3508 0 : pH[1] = Vec4( 0., 0., 0.5 * eCM * x1H, 0.5 * eCM * x1H);
3509 0 : pH[2] = Vec4( 0., 0., -0.5 * eCM * x2H, 0.5 * eCM * x2H);
3510 :
3511 : // Begin three-momentum rescaling to compensate for masses.
3512 0 : if (idMass[3] == 0 || idMass[4] == 0 || idMass[5] == 0) {
3513 0 : double p3S = p3cm.pAbs2();
3514 0 : double p4S = p4cm.pAbs2();
3515 0 : double p5S = p5cm.pAbs2();
3516 : double fac = 1.;
3517 : double e3, e4, e5, value, deriv;
3518 :
3519 : // Iterate rescaling solution five times, using Newton-Raphson.
3520 0 : for (int i = 0; i < NITERNR; ++i) {
3521 0 : e3 = sqrt(s3 + fac * p3S);
3522 0 : e4 = sqrt(s4 + fac * p4S);
3523 0 : e5 = sqrt(s5 + fac * p5S);
3524 0 : value = e3 + e4 + e5 - mHat;
3525 0 : deriv = 0.5 * (p3S / e3 + p4S / e4 + p5S / e5);
3526 0 : fac -= value / deriv;
3527 : }
3528 :
3529 : // Rescale momenta appropriately.
3530 0 : double facRoot = sqrt(fac);
3531 0 : p3cm.rescale3( facRoot );
3532 0 : p4cm.rescale3( facRoot );
3533 0 : p5cm.rescale3( facRoot );
3534 0 : p3cm.e( sqrt(s3 + fac * p3S) );
3535 0 : p4cm.e( sqrt(s4 + fac * p4S) );
3536 0 : p5cm.e( sqrt(s5 + fac * p5S) );
3537 0 : }
3538 :
3539 : // Outgoing partons initially in collision CM frame along beam axes.
3540 0 : pH[3] = p3cm;
3541 0 : pH[4] = p4cm;
3542 0 : pH[5] = p5cm;
3543 :
3544 : // Then boost them to overall CM frame
3545 0 : betaZ = (x1H - x2H)/(x1H + x2H);
3546 0 : pH[3].rot( theta, phi);
3547 0 : pH[4].rot( theta, phi);
3548 0 : pH[3].bst( 0., 0., betaZ);
3549 0 : pH[4].bst( 0., 0., betaZ);
3550 0 : pH[5].bst( 0., 0., betaZ);
3551 :
3552 : // Store average pT of three final particles for documentation.
3553 0 : pTH = (p3cm.pT() + p4cm.pT() + p5cm.pT()) / 3.;
3554 :
3555 : // Done.
3556 0 : return true;
3557 0 : }
3558 :
3559 : //==========================================================================
3560 :
3561 : // The PhaseSpace2to3yyycyl class.
3562 : // Phase space for 2 -> 3 QCD processes, 1 + 2 -> 3 + 4 + 5 set up in
3563 : // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut.
3564 : // Note: here cout is used for output, not os. Change??
3565 :
3566 : //--------------------------------------------------------------------------
3567 :
3568 : // Sample the phase space of the process.
3569 :
3570 : bool PhaseSpace2to3yyycyl::setupSampling() {
3571 :
3572 : // Phase space cuts specifically for 2 -> 3 QCD processes.
3573 0 : pTHat3Min = settingsPtr->parm("PhaseSpace:pTHat3Min");
3574 0 : pTHat3Max = settingsPtr->parm("PhaseSpace:pTHat3Max");
3575 0 : pTHat5Min = settingsPtr->parm("PhaseSpace:pTHat5Min");
3576 0 : pTHat5Max = settingsPtr->parm("PhaseSpace:pTHat5Max");
3577 0 : RsepMin = settingsPtr->parm("PhaseSpace:RsepMin");
3578 0 : R2sepMin = pow2(RsepMin);
3579 :
3580 : // If both beams are baryons then softer PDF's than for mesons/Pomerons.
3581 0 : hasBaryonBeams = ( beamAPtr->isBaryon() && beamBPtr->isBaryon() );
3582 :
3583 : // Work with massless partons.
3584 0 : for (int i = 0; i < 6; ++i) mH[i] = 0.;
3585 :
3586 : // Constrain to possible cuts at current CM energy and check consistency.
3587 0 : pT3Min = pTHat3Min;
3588 0 : pT3Max = pTHat3Max;
3589 0 : if (pT3Max < pT3Min) pT3Max = 0.5 * eCM;
3590 0 : pT5Min = pTHat5Min;
3591 0 : pT5Max = pTHat5Max;
3592 0 : if (pT5Max < pT5Min) pT5Max = 0.5 * eCM;
3593 0 : if (pT5Max > pT3Max || pT5Min > pT3Min || pT3Min + 2. * pT5Min > eCM) {
3594 0 : infoPtr->errorMsg("Error in PhaseSpace2to3yyycyl::setupSampling: "
3595 : "inconsistent pT limits in 3-body phase space");
3596 0 : return false;
3597 : }
3598 :
3599 : // Loop over some configurations where cross section could be maximal.
3600 : // In all cases all sum p_z = 0, for maximal PDF weights.
3601 : // Also pT3 and R45 are minimal, while pT5 may vary.
3602 0 : sigmaMx = 0.;
3603 0 : double pT5EffMax = min( pT5Max, 0.5 * pT3Min / cos(0.5 * RsepMin) );
3604 0 : double pT3EffMin = max( pT3Min, 2.0 * pT5Min * cos(0.5 * RsepMin) ) ;
3605 0 : double sinhR = sinh(0.5 * RsepMin);
3606 0 : double coshR = cosh(0.5 * RsepMin);
3607 0 : for (int iStep = 0; iStep < 120; ++iStep) {
3608 :
3609 : // First kind: |phi4 - phi5| = R, all p_z = 0, i.e. separation in phi.
3610 0 : if (iStep < 10) {
3611 0 : pT3 = pT3EffMin;
3612 0 : pT5 = pT5Min * pow( pT5EffMax / pT5Min, iStep / 9.);
3613 0 : double pTRat = pT5 / pT3;
3614 0 : double sin2Rsep = pow2( sin(RsepMin) );
3615 0 : double cosPhi35 = - cos(RsepMin) * sqrtpos(1. - sin2Rsep
3616 0 : * pow2(pTRat)) - sin2Rsep * pTRat;
3617 0 : cosPhi35 = max( cosPhi35, cos(M_PI - 0.5 * RsepMin) );
3618 0 : double sinPhi35 = sqrt(1. - pow2(cosPhi35));
3619 0 : pT4 = sqrt( pow2(pT3) + pow2(pT5) + 2. * pT3 * pT5 * cosPhi35);
3620 0 : p3cm = pT3 * Vec4( 1., 0., 0., 1.);
3621 0 : p4cm = Vec4(-pT3 - pT5 * cosPhi35, -pT5 * sinPhi35, 0., pT4);
3622 0 : p5cm = pT5 * Vec4( cosPhi35, sinPhi35, 0., 1.);
3623 :
3624 : // Second kind: |y4 - y5| = R, phi4 = phi5, i.e. separation in y.
3625 0 : } else {
3626 0 : pT5 = pT5Min * pow( pT5Max / pT5Min, iStep%10 / 9. );
3627 0 : pT3 = max( pT3Min, 2. * pT5);
3628 0 : pT4 = pT3 - pT5;
3629 0 : p4cm = pT4 * Vec4( -1., 0., sinhR, coshR );
3630 0 : p5cm = pT5 * Vec4( -1., 0., -sinhR, coshR );
3631 0 : y3 = -1.2 + 0.2 * (iStep/10);
3632 0 : p3cm = pT3 * Vec4( 1., 0., sinh(y3), cosh(y3));
3633 0 : betaZ = (p3cm.pz() + p4cm.pz() + p5cm.pz())
3634 0 : / (p3cm.e() + p4cm.e() + p5cm.e());
3635 0 : p3cm.bst( 0., 0., -betaZ);
3636 0 : p4cm.bst( 0., 0., -betaZ);
3637 0 : p5cm.bst( 0., 0., -betaZ);
3638 : }
3639 :
3640 : // Find cross section in chosen phase space point.
3641 0 : pInSum = p3cm + p4cm + p5cm;
3642 0 : x1H = pInSum.e() / eCM;
3643 0 : x2H = x1H;
3644 0 : sH = pInSum.m2Calc();
3645 0 : sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
3646 : 0., 0., 0., 1., 1., 1.);
3647 0 : sigmaNw = sigmaProcessPtr->sigmaPDF();
3648 :
3649 : // Multiply by Jacobian.
3650 0 : double flux = 1. /(8. * pow2(sH) * pow5(2. * M_PI));
3651 0 : double pTRng = pow2(M_PI)
3652 0 : * pow4(pT3) * (1./pow2(pT3Min) - 1./pow2(pT3Max))
3653 0 : * pow2(pT5) * 2.* log(pT5Max/pT5Min);
3654 0 : double yRng = 8. * log(eCM / pT3) * log(eCM / pT4) * log(eCM / pT5);
3655 0 : sigmaNw *= SAFETYMARGIN * flux * pTRng * yRng;
3656 :
3657 : // Update to largest maximum.
3658 0 : if (showSearch && sigmaNw > sigmaMx) cout << "\n New sigmamax is "
3659 0 : << scientific << setprecision(3) << sigmaNw << " for x1 = " << x1H
3660 0 : << " x2 = " << x2H << " sH = " << sH << endl << " p3 = " << p3cm
3661 0 : << " p4 = " << p4cm << " p5 = " << p5cm;
3662 0 : if (sigmaNw > sigmaMx) sigmaMx = sigmaNw;
3663 : }
3664 0 : sigmaPos = sigmaMx;
3665 :
3666 : // Done.
3667 : return true;
3668 0 : }
3669 :
3670 : //--------------------------------------------------------------------------
3671 :
3672 : // Sample the phase space of the process.
3673 :
3674 : bool PhaseSpace2to3yyycyl::trialKin(bool inEvent, bool) {
3675 :
3676 : // Allow for possibility that energy varies from event to event.
3677 0 : if (doEnergySpread) {
3678 0 : eCM = infoPtr->eCM();
3679 0 : s = eCM * eCM;
3680 0 : }
3681 0 : sigmaNw = 0.;
3682 :
3683 : // Constrain to possible cuts at current CM energy and check consistency.
3684 0 : pT3Min = pTHat3Min;
3685 0 : pT3Max = pTHat3Max;
3686 0 : if (pT3Max < pT3Min) pT3Max = 0.5 * eCM;
3687 0 : pT5Min = pTHat5Min;
3688 0 : pT5Max = pTHat5Max;
3689 0 : if (pT5Max < pT5Min) pT5Max = 0.5 * eCM;
3690 0 : if (pT5Max > pT3Max || pT5Min > pT3Min || pT3Min + 2. * pT5Min > eCM) {
3691 0 : infoPtr->errorMsg("Error in PhaseSpace2to3yyycyl::trialKin: "
3692 : "inconsistent pT limits in 3-body phase space");
3693 0 : return false;
3694 : }
3695 :
3696 : // Pick pT3 according to d^2(pT3)/pT3^4 and pT5 to d^2(pT5)/pT5^2.
3697 0 : pT3 = pT3Min * pT3Max / sqrt( pow2(pT3Min) +
3698 0 : rndmPtr->flat() * (pow2(pT3Max) - pow2(pT3Min)) );
3699 0 : pT5Max = min(pT5Max, pT3);
3700 0 : if (pT5Max < pT5Min) return false;
3701 0 : pT5 = pT5Min * pow( pT5Max / pT5Min, rndmPtr->flat() );
3702 :
3703 : // Pick azimuthal angles flat and reconstruct pT4, between pT3 and pT5.
3704 0 : phi3 = 2. * M_PI * rndmPtr->flat();
3705 0 : phi5 = 2. * M_PI * rndmPtr->flat();
3706 0 : pT4 = sqrt( pow2(pT3) + pow2(pT5) + 2. * pT3 * pT5 * cos(phi3 - phi5) );
3707 0 : if (pT4 > pT3 || pT4 < pT5) return false;
3708 0 : phi4 = atan2( -(pT3 * sin(phi3) + pT5 * sin(phi5)),
3709 0 : -(pT3 * cos(phi3) + pT5 * cos(phi5)) );
3710 :
3711 : // Pick rapidities flat in allowed ranges.
3712 0 : y3Max = log(eCM / pT3);
3713 0 : y4Max = log(eCM / pT4);
3714 0 : y5Max = log(eCM / pT5);
3715 0 : y3 = y3Max * (2. * rndmPtr->flat() - 1.);
3716 0 : y4 = y4Max * (2. * rndmPtr->flat() - 1.);
3717 0 : y5 = y5Max * (2. * rndmPtr->flat() - 1.);
3718 :
3719 : // Reject some events at large rapidities to improve efficiency.
3720 : // (Works for baryons, not pions or Pomerons if they have hard PDF's.)
3721 0 : double WTy = (hasBaryonBeams) ? (1. - pow2(y3/y3Max))
3722 0 : * (1. - pow2(y4/y4Max)) * (1. - pow2(y5/y5Max)) : 1.;
3723 0 : if (WTy < rndmPtr->flat()) return false;
3724 :
3725 : // Check that any pair separated more then RsepMin in (y, phi) space.
3726 0 : dphi = abs(phi3 - phi4);
3727 0 : if (dphi > M_PI) dphi = 2. * M_PI - dphi;
3728 0 : if (pow2(y3 - y4) + pow2(dphi) < R2sepMin) return false;
3729 0 : dphi = abs(phi3 - phi5);
3730 0 : if (dphi > M_PI) dphi = 2. * M_PI - dphi;
3731 0 : if (pow2(y3 - y5) + pow2(dphi) < R2sepMin) return false;
3732 0 : dphi = abs(phi4 - phi5);
3733 0 : if (dphi > M_PI) dphi = 2. * M_PI - dphi;
3734 0 : if (pow2(y4 - y5) + pow2(dphi) < R2sepMin) return false;
3735 :
3736 : // Reconstruct all four-vectors.
3737 0 : pH[3] = pT3 * Vec4( cos(phi3), sin(phi3), sinh(y3), cosh(y3) );
3738 0 : pH[4] = pT4 * Vec4( cos(phi4), sin(phi4), sinh(y4), cosh(y4) );
3739 0 : pH[5] = pT5 * Vec4( cos(phi5), sin(phi5), sinh(y5), cosh(y5) );
3740 0 : pInSum = pH[3] + pH[4] + pH[5];
3741 :
3742 : // Check that x values physical and sHat in allowed range.
3743 0 : x1H = (pInSum.e() + pInSum.pz()) / eCM;
3744 0 : x2H = (pInSum.e() - pInSum.pz()) / eCM;
3745 0 : if (x1H >= 1. || x2H >= 1.) return false;
3746 0 : sH = pInSum.m2Calc();
3747 0 : if ( sH < pow2(mHatGlobalMin) ||
3748 0 : (mHatGlobalMax > mHatGlobalMin && sH > pow2(mHatGlobalMax)) )
3749 0 : return false;
3750 :
3751 : // Boost four-vectors to rest frame of collision.
3752 0 : betaZ = (x1H - x2H)/(x1H + x2H);
3753 0 : p3cm = pH[3]; p3cm.bst( 0., 0., -betaZ);
3754 0 : p4cm = pH[4]; p4cm.bst( 0., 0., -betaZ);
3755 0 : p5cm = pH[5]; p5cm.bst( 0., 0., -betaZ);
3756 :
3757 : // Find cross section in chosen phase space point.
3758 0 : sigmaProcessPtr->set3Kin( x1H, x2H, sH, p3cm, p4cm, p5cm,
3759 : 0., 0., 0., 1., 1., 1.);
3760 0 : sigmaNw = sigmaProcessPtr->sigmaPDF();
3761 :
3762 : // Multiply by Jacobian. Correct for rejection of large rapidities.
3763 0 : double flux = 1. /(8. * pow2(sH) * pow5(2. * M_PI));
3764 0 : double yRng = 8. * y3Max * y4Max * y5Max;
3765 0 : double pTRng = pow2(M_PI)
3766 0 : * pow4(pT3) * (1./pow2(pT3Min) - 1./pow2(pT3Max))
3767 0 : * pow2(pT5) * 2.* log(pT5Max/pT5Min);
3768 0 : sigmaNw *= flux * yRng * pTRng / WTy;
3769 :
3770 : // Allow possibility for user to modify cross section.
3771 0 : if (canModifySigma) sigmaNw
3772 0 : *= userHooksPtr->multiplySigmaBy( sigmaProcessPtr, this, inEvent);
3773 0 : if (canBiasSelection) sigmaNw
3774 0 : *= userHooksPtr->biasSelectionBy( sigmaProcessPtr, this, inEvent);
3775 0 : if (canBias2Sel) sigmaNw *= pow( pTH / bias2SelRef, bias2SelPow);
3776 :
3777 : // Check if maximum violated.
3778 0 : newSigmaMx = false;
3779 0 : if (sigmaNw > sigmaMx) {
3780 0 : infoPtr->errorMsg("Warning in PhaseSpace2to3yyycyl::trialKin: "
3781 : "maximum for cross section violated");
3782 :
3783 : // Violation strategy 1: increase maximum (always during initialization).
3784 0 : if (increaseMaximum || !inEvent) {
3785 0 : double violFact = SAFETYMARGIN * sigmaNw / sigmaMx;
3786 0 : sigmaMx = SAFETYMARGIN * sigmaNw;
3787 0 : newSigmaMx = true;
3788 0 : if (showViolation) {
3789 0 : if (violFact < 9.99) cout << fixed;
3790 0 : else cout << scientific;
3791 0 : cout << " PYTHIA Maximum for " << sigmaProcessPtr->name()
3792 0 : << " increased by factor " << setprecision(3) << violFact
3793 0 : << " to " << scientific << sigmaMx << endl;
3794 0 : }
3795 :
3796 : // Violation strategy 2: weight event (done in ProcessContainer).
3797 0 : } else if (showViolation && sigmaNw > sigmaPos) {
3798 0 : double violFact = sigmaNw / sigmaMx;
3799 0 : if (violFact < 9.99) cout << fixed;
3800 0 : else cout << scientific;
3801 0 : cout << " PYTHIA Maximum for " << sigmaProcessPtr->name()
3802 0 : << " exceeded by factor " << setprecision(3) << violFact << endl;
3803 0 : sigmaPos = sigmaNw;
3804 0 : }
3805 : }
3806 :
3807 : // Check if negative cross section.
3808 0 : if (sigmaNw < sigmaNeg) {
3809 0 : infoPtr->errorMsg("Warning in PhaseSpace2to3yyycyl::trialKin:"
3810 0 : " negative cross section set 0", "for " + sigmaProcessPtr->name() );
3811 0 : sigmaNeg = sigmaNw;
3812 :
3813 : // Optional printout of (all) violations.
3814 0 : if (showViolation) cout << " PYTHIA Negative minimum for "
3815 0 : << sigmaProcessPtr->name() << " changed to " << scientific
3816 0 : << setprecision(3) << sigmaNeg << endl;
3817 : }
3818 0 : if (sigmaNw < 0.) sigmaNw = 0.;
3819 :
3820 :
3821 : // Done.
3822 : return true;
3823 0 : }
3824 :
3825 : //--------------------------------------------------------------------------
3826 :
3827 : // Construct the final kinematics of the process: not much left
3828 :
3829 : bool PhaseSpace2to3yyycyl::finalKin() {
3830 :
3831 : // Work with massless partons.
3832 0 : for (int i = 0; i < 6; ++i) mH[i] = 0.;
3833 :
3834 : // Ibncoming partons to collision.
3835 0 : pH[1] = 0.5 * (pInSum.e() + pInSum.pz()) * Vec4( 0., 0., 1., 1.);
3836 0 : pH[2] = 0.5 * (pInSum.e() - pInSum.pz()) * Vec4( 0., 0., -1., 1.);
3837 :
3838 : // Some quantities meaningless for 2 -> 3. pT defined as average value.
3839 0 : tH = 0.;
3840 0 : uH = 0.;
3841 0 : pTH = (pH[3].pT() + pH[4].pT() + pH[5].pT()) / 3.;
3842 0 : theta = 0.;
3843 0 : phi = 0.;
3844 :
3845 0 : return true;
3846 : }
3847 :
3848 :
3849 : //==========================================================================
3850 :
3851 : // The PhaseSpaceLHA class.
3852 : // A derived class for Les Houches events.
3853 : // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
3854 :
3855 : //--------------------------------------------------------------------------
3856 :
3857 : // Constants: could be changed here if desired, but normally should not.
3858 : // These are of technical nature, as described for each.
3859 :
3860 : // LHA convention with cross section in pb forces conversion to mb.
3861 : const double PhaseSpaceLHA::CONVERTPB2MB = 1e-9;
3862 :
3863 : //--------------------------------------------------------------------------
3864 :
3865 : // Find maximal cross section for comparison with internal processes.
3866 :
3867 : bool PhaseSpaceLHA::setupSampling() {
3868 :
3869 : // Find which strategy Les Houches events are produced with.
3870 0 : strategy = lhaUpPtr->strategy();
3871 0 : stratAbs = abs(strategy);
3872 0 : if (strategy == 0 || stratAbs > 4) {
3873 0 : ostringstream stratCode;
3874 0 : stratCode << strategy;
3875 0 : infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: unknown "
3876 0 : "Les Houches Accord weighting stategy", stratCode.str());
3877 : return false;
3878 0 : }
3879 :
3880 : // Number of contributing processes.
3881 0 : nProc = lhaUpPtr->sizeProc();
3882 :
3883 : // Loop over all processes. Read out maximum and cross section.
3884 0 : xMaxAbsSum = 0.;
3885 0 : xSecSgnSum = 0.;
3886 0 : int idPr;
3887 0 : double xMax, xSec, xMaxAbs;
3888 0 : for (int iProc = 0 ; iProc < nProc; ++iProc) {
3889 0 : idPr = lhaUpPtr->idProcess(iProc);
3890 0 : xMax = lhaUpPtr->xMax(iProc);
3891 0 : xSec = lhaUpPtr->xSec(iProc);
3892 :
3893 : // Check for inconsistencies between strategy and stored values.
3894 0 : if ( (strategy == 1 || strategy == 2) && xMax < 0.) {
3895 0 : infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: "
3896 : "negative maximum not allowed");
3897 0 : return false;
3898 : }
3899 0 : if ( ( strategy == 2 || strategy == 3) && xSec < 0.) {
3900 0 : infoPtr->errorMsg("Error in PhaseSpaceLHA::setupSampling: "
3901 : "negative cross section not allowed");
3902 0 : return false;
3903 : }
3904 :
3905 : // Store maximal cross sections for later choice.
3906 0 : if (stratAbs == 1) xMaxAbs = abs(xMax);
3907 0 : else if (stratAbs < 4) xMaxAbs = abs(xSec);
3908 0 : else xMaxAbs = 1.;
3909 0 : idProc.push_back( idPr );
3910 0 : xMaxAbsProc.push_back( xMaxAbs );
3911 :
3912 : // Find sum and convert to mb.
3913 0 : xMaxAbsSum += xMaxAbs;
3914 0 : xSecSgnSum += xSec;
3915 : }
3916 0 : sigmaMx = xMaxAbsSum * CONVERTPB2MB;
3917 0 : sigmaSgn = xSecSgnSum * CONVERTPB2MB;
3918 :
3919 : // Done.
3920 0 : return true;
3921 :
3922 0 : }
3923 :
3924 : //--------------------------------------------------------------------------
3925 :
3926 : // Construct the next process, by interface to Les Houches class.
3927 :
3928 : bool PhaseSpaceLHA::trialKin( bool, bool repeatSame ) {
3929 :
3930 : // Must select process type in some cases.
3931 : int idProcNow = 0;
3932 0 : if (repeatSame) idProcNow = idProcSave;
3933 0 : else if (stratAbs <= 2) {
3934 0 : double xMaxAbsRndm = xMaxAbsSum * rndmPtr->flat();
3935 : int iProc = -1;
3936 0 : do xMaxAbsRndm -= xMaxAbsProc[++iProc];
3937 0 : while (xMaxAbsRndm > 0. && iProc < nProc - 1);
3938 0 : idProcNow = idProc[iProc];
3939 0 : }
3940 :
3941 : // Generate Les Houches event. Return if fail (= end of file).
3942 0 : bool physical = lhaUpPtr->setEvent(idProcNow);
3943 0 : if (!physical) return false;
3944 :
3945 : // Find which process was generated.
3946 0 : int idPr = lhaUpPtr->idProcess();
3947 : int iProc = 0;
3948 0 : for (int iP = 0; iP < int(idProc.size()); ++iP)
3949 0 : if (idProc[iP] == idPr) iProc = iP;
3950 0 : idProcSave = idPr;
3951 :
3952 : // Extract cross section and rescale according to strategy.
3953 0 : double wtPr = lhaUpPtr->weight();
3954 0 : if (stratAbs == 1) sigmaNw = wtPr * CONVERTPB2MB
3955 0 : * xMaxAbsSum / xMaxAbsProc[iProc];
3956 0 : else if (stratAbs == 2) sigmaNw = (wtPr / abs(lhaUpPtr->xMax(iProc)))
3957 0 : * sigmaMx;
3958 0 : else if (strategy == 3) sigmaNw = sigmaMx;
3959 0 : else if (strategy == -3 && wtPr > 0.) sigmaNw = sigmaMx;
3960 0 : else if (strategy == -3) sigmaNw = -sigmaMx;
3961 0 : else if (stratAbs == 4) sigmaNw = wtPr * CONVERTPB2MB;
3962 :
3963 : // Set x scales.
3964 0 : x1H = lhaUpPtr->x1();
3965 0 : x2H = lhaUpPtr->x2();
3966 :
3967 : // Done.
3968 : return true;
3969 :
3970 0 : }
3971 :
3972 : //==========================================================================
3973 :
3974 : } // end namespace Pythia8
|
