Line data Source code
1 : // HardDiffraction.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 : // Author: Christine O. Rasmussen.
7 :
8 : // Function definitions (not found in the header) for the
9 : // HardDiffraction class.
10 :
11 : #include "Pythia8/HardDiffraction.h"
12 : namespace Pythia8 {
13 :
14 : //==========================================================================
15 :
16 : // Constants: could be changed here if desired, but normally should not.
17 : // These are of technical nature, as described for each.
18 :
19 : // Lower limit on PDF value in order to avoid division by zero.
20 : const double HardDiffraction::TINYPDF = 1e-10;
21 :
22 : // Ficticious Pomeron mass to leave room for beam remnant
23 : const double HardDiffraction::POMERONMASS = 1.;
24 :
25 : // Proton mass
26 : const double HardDiffraction::PROTONMASS = 0.938;
27 :
28 : //--------------------------------------------------------------------------
29 :
30 : void HardDiffraction::init(Info* infoPtrIn, Settings& settingsPtrIn,
31 : Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
32 : BeamParticle* beamPomAPtrIn, BeamParticle* beamPomBPtrIn) {
33 :
34 : // Store pointers.
35 0 : infoPtr = infoPtrIn;
36 0 : settings = settingsPtrIn;
37 0 : rndmPtr = rndmPtrIn;
38 0 : beamAPtr = beamAPtrIn;
39 0 : beamBPtr = beamBPtrIn;
40 0 : beamPomAPtr = beamPomAPtrIn;
41 0 : beamPomBPtr = beamPomBPtrIn;
42 :
43 : // Set diffraction parameters.
44 0 : pomSet = settings.mode("PDF:PomSet");
45 0 : pomFlux = settings.mode("Diffraction:PomFlux");
46 :
47 : // Read out some properties of beams to allow shorthand.
48 0 : idA = (beamAPtr != 0) ? beamAPtr->id() : 0;
49 0 : idB = (beamBPtr != 0) ? beamBPtr->id() : 0;
50 0 : mA = (beamAPtr != 0) ? beamAPtr->m() : 0.;
51 0 : mB = (beamBPtr != 0) ? beamBPtr->m() : 0.;
52 :
53 : // Set up Pomeron flux constants.
54 0 : if (pomFlux == 1) {
55 0 : double sigmaRefPomP = settings.parm("Diffraction:sigmaRefPomP");
56 0 : normPom = pow2(sigmaRefPomP) * 0.02;
57 0 : b0 = 2.3;
58 0 : ap = 0.25;
59 0 : } else if (pomFlux == 2) {
60 0 : normPom = 1/2.3;
61 0 : A1 = 6.38;
62 0 : A2 = 0.424;
63 0 : a1 = 8.;
64 0 : a2 = 3.;
65 0 : } else if (pomFlux == 3) {
66 0 : normPom = 1.99;
67 0 : a0 = 1.085;
68 0 : ap = 0.25;
69 0 : a1 = 4.7;
70 0 : } else if (pomFlux == 4) {
71 0 : normPom = 0.74;
72 0 : ap = 0.06;
73 0 : a0 = 1.1182;
74 0 : A1 = 0.27;
75 0 : a1 = 8.38;
76 0 : A2 = 0.56;
77 0 : a2 = 3.78;
78 0 : A3 = 0.18;
79 0 : a3 = 1.36;
80 0 : } else if (pomFlux == 5) {
81 0 : normPom = 0.858;
82 0 : A1 = 0.9;
83 0 : a1 = 4.6;
84 0 : A2 = 0.1;
85 0 : a2 = 0.6;
86 0 : a0 = 1.104;
87 0 : ap = 0.25;
88 0 : } else if (pomFlux == 6 || pomFlux == 7) {
89 0 : normPom = 1.57285;
90 0 : ap = 0.06;
91 0 : b0 = 5.5;
92 0 : if (pomFlux == 6) a0 = 1.1182;
93 0 : else a0 = 1.1110;
94 : }
95 :
96 : // Initialise Pomeron values to zero.
97 0 : xPomA = tPomA = thetaPomA = 0.;
98 0 : xPomB = tPomB = thetaPomB = 0.;
99 :
100 : // Done.
101 0 : }
102 :
103 : //--------------------------------------------------------------------------
104 :
105 : bool HardDiffraction::isDiffractive( int iBeamIn, int partonIn, double xIn,
106 : double Q2In, double xfIncIn) {
107 :
108 : // Store incoming values.
109 0 : iBeam = iBeamIn;
110 : int parton = partonIn;
111 : double x = xIn;
112 : double Q2 = Q2In;
113 : double xfInc = xfIncIn;
114 0 : tmpPDFPtr = (iBeam == 1) ? beamPomAPtr : beamPomBPtr;
115 :
116 : // Return false if value of inclusive PDF is zero.
117 0 : if (xfInc < TINYPDF) {
118 0 : infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive: "
119 : "inclusive PDF is zero");
120 0 : return false;
121 : }
122 :
123 : // Generate an xNow according to 1 / x.
124 0 : double xNow = pow(xIn, rndmPtr->flat());
125 :
126 : // Overestimated function:
127 : // g(xP) = c/xP * x*f_{P/p}(x) * x*f_{i/P}(x, Q^2_{max})
128 : // G(x) = int_x^1 dxP g(xP)
129 : // = c * x*f_{P/p}(x) * x*f_{i/P}(x, Q^2_{max}) * log( 1/x )
130 : // True function:
131 : // f(xP) = 1/xP * xP *f_{P/p}(xP) * x/xP * f_{i/P}(x/xP, Q^2)
132 : // F(x) = int_x^1 dxP f(xP)
133 : // We have diffraction if G(x)/xfInc > R and if f(xP)/g(xP) > R =>
134 : // R < (G(x) * f(xP)) / (g(xP) * xfInc) = (log(1/x) * f(xP)) / (xP * xfInc)
135 0 : double over = log(1./x) * xfPom(xNow) * tmpPDFPtr->xf(parton, x/xNow, Q2);
136 0 : if (over > xfInc) {
137 0 : stringstream msg;
138 0 : msg << " Weight above unity with parton " << parton;
139 0 : infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive:", msg.str());
140 0 : }
141 :
142 : // Discard if overestimate/inclusive PDF is less than random number.
143 0 : if (over/xfInc < rndmPtr->flat()) return false;
144 :
145 : // Make sure there is momentum left for beam remnant
146 0 : double m2Diff = xNow * pow2( infoPtr->eCM());
147 0 : double mDiff = sqrt(m2Diff);
148 0 : double mDiffA = (iBeam == 1) ? 0. : PROTONMASS;
149 0 : double mDiffB = (iBeam == 2) ? 0. : PROTONMASS;
150 0 : double m2DiffA = mDiffA * mDiffA;
151 0 : double m2DiffB = mDiffB * mDiffB;
152 0 : double eDiff = (iBeam == 1) ? 0.5 * (m2Diff + m2DiffA - m2DiffB) / mDiff :
153 0 : 0.5 * (m2Diff + m2DiffB - m2DiffA) / mDiff;
154 0 : if ( 1. - x / xNow < POMERONMASS / eDiff) {
155 0 : infoPtr->errorMsg("Warning in HardDiffraction::isDiffractive: "
156 : "No momentum left for beam remnant.");
157 0 : return false;
158 : }
159 :
160 : // The chosen xNow is accepted, now find t and theta.
161 0 : double tNow = pickTNow(xNow);
162 0 : double thetaNow = getThetaNow(xNow, tNow);
163 :
164 : // Set the chosen diffractive values, to be able to use them later.
165 0 : if (iBeam == 1) {
166 0 : xPomA = xNow;
167 0 : tPomA = tNow;
168 0 : thetaPomA = thetaNow;
169 0 : } else {
170 0 : xPomB = xNow;
171 0 : tPomB = tNow;
172 0 : thetaPomB = thetaNow;
173 : }
174 :
175 : // Done.
176 : return true;
177 0 : }
178 :
179 : //--------------------------------------------------------------------------
180 :
181 : // Return x*f_P/p( x), ie. Pomeron flux inside proton, integrated over t.
182 :
183 : double HardDiffraction::xfPom(double xIn) {
184 :
185 : // Setup t range.
186 0 : pair<double, double> tLim = tRange(xIn);
187 : double tMin = tLim.first;
188 : double tMax = tLim.second;
189 : double x = xIn;
190 : double xFlux = 0.;
191 :
192 : // Schuler-Sjöstrand Pomeron flux, see Phys. Rev. D.49 (1994) 2259.
193 : // flux = normPom * 1/x * exp(2t(2.3 + 0.25 * log(1/x)))
194 : // => x * flux = normPom * exp(2t(2.3 + 0.25*log(1/x)))
195 0 : if (pomFlux == 1) {
196 0 : double b = b0 + ap * log(1./x);
197 0 : xFlux = normPom/(2.*b) * ( exp(2.*b*tMax) - exp(2.*b*tMin));
198 0 : }
199 :
200 : // Bruni-Ingelman Pomeron flux, see Phys. Lett. B311 (1993) 317.
201 : // flux = normPom * (1/x) * (6.38 *exp(8*t)+ 0.424 * exp(3*t))
202 : // => x * flux = normPom * (6.38 *exp(8*t)+ 0.424 * exp(3*t))
203 0 : else if (pomFlux == 2) {
204 0 : xFlux = normPom * (A1/a1 * (exp(a1*tMax) - exp(a1*tMin))
205 0 : + A2/a2 * (exp(a2*tMax) - exp(a2*tMin)));
206 0 : }
207 :
208 : // Streng-Berger Pomeron flux, see Comp. Phys. Comm. 86 (1995) 147.
209 : // flux = normPom * x^(1 - 2*alpha(t)) * exp(-R_N^2 * t)
210 : // => x * flux = normPom * x^(2 - 2*alpha(t)) * exp(-R_N^2 * t)
211 0 : else if (pomFlux == 3) {
212 0 : double b = (a1 + 2. * ap * log(1./x));
213 0 : xFlux = normPom * exp(log(1./x) * (2.*a0 - 2.));
214 0 : xFlux *= (exp(b*tMax) - exp(b*tMin))/b;
215 0 : }
216 :
217 : // Donnachie-Landshoff Pomeron flux, see Phys. Lett. B 191 (1987) 309.
218 : // flux = beta^2(0)/(16 pi) x^(1 - 2*alpha(t)) F_1(t)^2 with
219 : // F_1(t)^2 = 0.27 exp(8.38 t) + 0.56 exp(3.78 t) + 0.18 exp(1.36 t)
220 : // = (4m_p^2-2.8t)^2/(4m_p^2-t)^2*(1/(1-t/0.7))^4
221 : // => x * flux = beta^2(0)/(16 pi) * x^(2 - 2*\alpha(t)) F_1(t)^2
222 0 : else if (pomFlux == 4) {
223 0 : double Q = 2. * ap * log(1./x);
224 0 : xFlux = normPom * exp(log(1./x) * (2.*a0 - 2.));
225 0 : xFlux *= (A1/(Q + a1) * (exp((Q + a1)*tMax) - exp((Q + a1)*tMin))
226 0 : + A2/(Q + a2) * (exp((Q + a2)*tMax) - exp((Q + a2)*tMin))
227 0 : + A3/(Q + a3) * (exp((Q + a3)*tMax) - exp((Q + a3)*tMin)));
228 0 : }
229 :
230 : // MBR Pomeron flux, see arXiv:1205.1446v2 [hep-ph] 2012.
231 : // flux = normPom * F_1(t)^2 * exp((2 * alpha(t) - 1)*log(1/x))
232 : // F_1(t)^2 = 0.9 exp(4.6 t) + 0.1 exp(0.6 t)
233 : // = (4m_p^2-2.8t)^2/(4m_p^2-t)^2*(1/(1-t/0.7))^4
234 : // => x * flux = normPom * F_1(t)^2 * exp( 2*(alpha(t) -1)*log(1/x))
235 0 : else if (pomFlux == 5) {
236 0 : double Q = 2. * ap * log(1./x);
237 0 : xFlux = normPom * exp(log(1./x) * ( 2.*a0 - 2.));
238 0 : xFlux *= (A1/(Q + a1) * (exp((Q + a1)*tMax) - exp((Q + a1)*tMin))
239 0 : + A2/(Q + a2) * (exp((Q + a2)*tMax) - exp((Q + a2)*tMin)));
240 0 : }
241 :
242 : // H1 Pomeron flux, see Eur. Phys. J. C48 (2006) 715, ibid. 749
243 : // flux = normPom * exp(B_Pom*t)/x^(2*\alpha(t)-1)
244 : // => x * flux = normPom * exp(B_Pom * t) / x^(2*\alpha(t)-2)
245 0 : else if (pomFlux == 6 || pomFlux == 7) {
246 0 : double b = b0 + 2. * ap * log(1./x);
247 0 : xFlux = normPom * exp(log(1./x) * ( 2.*a0 - 2.));
248 0 : xFlux *= (exp(b*tMax) - exp(b*tMin))/b;
249 0 : }
250 :
251 : // Done
252 0 : return xFlux;
253 : }
254 :
255 : //--------------------------------------------------------------------------
256 :
257 : // Pick a t value according to different Pomeron flux parametrizations.
258 :
259 : double HardDiffraction::pickTNow(double xIn) {
260 :
261 : // Get kinematical limits for t. initial values.
262 0 : pair<double, double> tLim = HardDiffraction::tRange(xIn);
263 : double tMin = tLim.first;
264 : double tMax = tLim.second;
265 : double tTmp = 0.;
266 0 : double rndm = rndmPtr->flat();
267 :
268 : // Schuler-Sjöstrndm Pomeron flux, see Phys. Rev. D.49 (1994) 2259.
269 0 : if (pomFlux == 1) {
270 0 : double b = b0 + ap * log(1./xIn);
271 0 : tTmp = log( rndm*exp(2.*b*tMin) + (1. - rndm)*exp(2.*b*tMax))/(2.*b);
272 0 : }
273 :
274 : // Bruni-Ingelman Pomeron flux, see Phys. Lett. B311 (1993) 317.
275 0 : else if (pomFlux == 2) {
276 0 : double prob1 = A1/a1 * (exp(a1*tMax) - exp(a1*tMin));
277 0 : double prob2 = A2/a2 * (exp(a2*tMax) - exp(a2*tMin));
278 0 : prob1 /= (prob1 + prob2);
279 0 : tTmp = (prob1 > rndmPtr->flat())
280 0 : ? log( rndm * exp(a1*tMin) + (1. - rndm) * exp(a1*tMax))/a1
281 0 : : log( rndm * exp(a2*tMin) + (1. - rndm) * exp(a2*tMax))/a2;
282 0 : }
283 :
284 : // Streng-Berger Pomeron flux, see Comp. Phys. Comm. 86 (1995) 147.
285 0 : else if (pomFlux == 3) {
286 0 : double b = (2. * ap * log(1./xIn) + a1);
287 0 : tTmp = log( rndm * exp(b*tMin) + (1. - rndm) * exp(b*tMax))/b;
288 0 : }
289 :
290 : // Donnachie-Landshoff Pomeron flux, see Phys. Lett. B 191 (1987) 309.
291 0 : else if (pomFlux == 4) {
292 0 : double b1 = 2. * ap * log(1./xIn) + a1;
293 0 : double b2 = 2. * ap * log(1./xIn) + a2;
294 0 : double b3 = 2. * ap * log(1./xIn) + a3;
295 0 : double prob1 = A1/b1 * ( exp(b1*tMax) - exp(b1*tMin));
296 0 : double prob2 = A2/b2 * ( exp(b2*tMax) - exp(b2*tMin));
297 0 : double prob3 = A3/b3 * ( exp(b3*tMax) - exp(b3*tMin));
298 0 : double rndmProb = (prob1 + prob2 + prob3) * rndmPtr->flat() ;
299 0 : if (rndmProb < prob1)
300 0 : tTmp = log( rndm * exp(b1*tMin) + (1. - rndm) * exp(b1*tMax))/b1;
301 0 : else if ( rndmProb < prob1 + prob2)
302 0 : tTmp = log( rndm * exp(b2*tMin) + (1. - rndm) * exp(b2*tMax))/b2;
303 : else
304 0 : tTmp = log( rndm * exp(b3*tMin) + (1. - rndm) * exp(b3*tMax))/b3;
305 0 : }
306 :
307 : // MBR Pomeron flux, see arXiv:1205.1446v2 [hep-ph] 2012.
308 0 : else if (pomFlux == 5) {
309 0 : double b1 = a1 + 2. * ap * log(1./xIn);
310 0 : double b2 = a2 + 2. * ap * log(1./xIn);
311 0 : double prob1 = A1/b1 * (exp(b1*tMax) - exp(b1*tMin));
312 0 : double prob2 = A2/b2 * (exp(b2*tMax) - exp(b2*tMin));
313 0 : prob1 /= (prob1 + prob2);
314 0 : tTmp = (prob1 > rndmPtr->flat())
315 0 : ? log( rndm * exp(b1*tMin) + (1. - rndm) * exp(b1*tMax))/b1
316 0 : : log( rndm * exp(b2*tMin) + (1. - rndm) * exp(b2*tMax))/b2;
317 0 : }
318 :
319 : // H1 Pomeron flux, see Eur. Phys. J. C48 (2006) 715, ibid. 749
320 0 : else if (pomFlux == 6 || pomFlux == 7){
321 0 : double b = b0 + 2. * ap * log(1./xIn);
322 0 : tTmp = log( rndm * exp(b*tMin) + (1. - rndm) * exp(b*tMax))/b;
323 0 : }
324 :
325 : // Done.
326 0 : return tTmp;
327 : }
328 :
329 : //--------------------------------------------------------------------------
330 :
331 : // Return x*f_P/p( x, t), ie. Pomeron flux inside proton differential in t.
332 :
333 : double HardDiffraction::xfPomWithT(double xIn, double tIn) {
334 :
335 : // Initial values.
336 : double x = xIn;
337 : double t = tIn;
338 : double xFlux = 0.;
339 :
340 : // Schuler-Sjöstrand Pomeron flux, see Phys. Rev. D.49 (1994) 2259.
341 0 : if (pomFlux == 1) {
342 0 : double b = b0 + ap * log(1./x);
343 0 : xFlux = normPom * exp( 2.*b*t);
344 0 : }
345 :
346 : // Bruni-Ingelman Pomeron flux, see Phys. Lett. B311 (1993) 317.
347 0 : else if (pomFlux == 2)
348 0 : xFlux = normPom * (A1 * exp(a1*t) + A2 * exp(a2*t));
349 :
350 : // Streng-Berger Pomeron flux, see Comp. Phys. Comm. 86 (1995) 147.
351 0 : else if (pomFlux == 3) {
352 0 : xFlux = normPom * exp(log(1./x) * (2.*a0 - 2.))
353 0 : * exp(t * (a1 + 2.*ap*log(1./x)));
354 0 : }
355 :
356 : // Donnachie-Landshoff Pomeron flux, see Phys. Lett. B 191 (1987) 309.
357 0 : else if (pomFlux == 4){
358 0 : double sqrF1 = A1 * exp(a1*t) + A2 * exp(a2*t) + A3 * exp(a3*t);
359 0 : xFlux = normPom * pow(x, 2. + 2. * (a0 + ap*t)) * sqrF1;
360 0 : }
361 :
362 : // MBR Pomeron flux, see arXiv:1205.1446v2 [hep-ph] 2012.
363 0 : else if (pomFlux == 5) {
364 0 : double sqrF1 = A1 * exp(a1*t) + A2 * exp(a2*t);
365 0 : xFlux = normPom * sqrF1 * exp(log(1./x) * (-2. + a0 + ap*t));
366 0 : }
367 :
368 : // H1 Pomeron flux, see Eur. Phys. J. C48 (2006) 715, ibid. 749
369 0 : else if (pomFlux == 6 || pomFlux == 7)
370 0 : xFlux = normPom * exp(b0*t)/pow(x, 2. * (a0 + ap*t) - 2.);
371 :
372 : // Done
373 0 : return xFlux;
374 : }
375 :
376 : //--------------------------------------------------------------------------
377 :
378 : // Set up t range. See p. 113 of 6.4 manual.
379 :
380 : pair<double, double> HardDiffraction::tRange(double xIn) {
381 :
382 : // Set up diffractive masses.
383 0 : double eCM = infoPtr->eCM();
384 0 : s = eCM * eCM;
385 0 : s1 = pow2(mA);
386 0 : s2 = pow2(mB);
387 0 : s3 = (iBeam == 1) ? xIn * s : s1;
388 0 : s4 = (iBeam == 2) ? xIn * s : s2;
389 :
390 : // Calculate kinematics.
391 0 : double lambda12 = sqrtpos(pow2(s - s1 - s2) - 4. * s1 * s2);
392 0 : double lambda34 = sqrtpos(pow2(s - s3 - s4) - 4. * s3 * s4);
393 0 : double tmp1 = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4) / s;
394 0 : double tmp2 = lambda12 * lambda34 / s;
395 0 : double tmp3 = (s1 + s4 - s2 - s3) * (s1 * s4 - s2 * s3) / s
396 0 : + (s3 - s1) * (s4 - s2);
397 0 : double tMin = -0.5 * (tmp1 + tmp2);
398 0 : double tMax = tmp3 / tMin;
399 :
400 : // Done.
401 0 : return make_pair(tMin, tMax);
402 0 : }
403 :
404 : //--------------------------------------------------------------------------
405 :
406 : // Get the scattering angle from the chosen t.
407 :
408 : double HardDiffraction::getThetaNow( double xIn, double tIn) {
409 :
410 : // Set up diffractive masses.
411 0 : double eCM = infoPtr->eCM();
412 0 : s = eCM * eCM;
413 0 : s1 = pow2( mA);
414 0 : s2 = pow2( mB);
415 0 : s3 = (iBeam == 1) ? xIn * s : s1;
416 0 : s4 = (iBeam == 2) ? xIn * s : s2;
417 :
418 : // Find theta from the chosen t.
419 0 : double lambda12 = sqrtpos(pow2(s - s1 - s2) - 4. * s1 * s2);
420 0 : double lambda34 = sqrtpos(pow2(s - s3 - s4) - 4. * s3 * s4);
421 0 : double tmp1 = s - (s1 + s2 + s3 + s4) + (s1 - s2) * (s3 - s4)/s;
422 0 : double tmp2 = lambda12 * lambda34 / s;
423 0 : double tmp3 = (s1 + s4 - s2 - s3) * (s1 * s4 - s2 * s3) / s
424 0 : + (s3 - s1) * (s4 - s2);
425 0 : double cosTheta = min(1., max(-1., (tmp1 + 2. * tIn) / tmp2));
426 0 : double sinTheta = 2. * sqrtpos( -(tmp3 + tmp1 * tIn + tIn * tIn) ) / tmp2;
427 0 : double theta = asin( min(1., sinTheta));
428 0 : if (cosTheta < 0.) theta = M_PI - theta;
429 :
430 : // Done.
431 0 : return theta;
432 0 : }
433 :
434 : //==========================================================================
435 :
436 : } // end namespace Pythia8
|
