cepgen is hosted by Hepforge, IPPP Durham
CepGen 1.2.5
Central exclusive processes event generator
Loading...
Searching...
No Matches
nucl_to_ff.f
Go to the documentation of this file.
1 function nucl_to_ff()
2 implicit none
3 double precision nucl_to_ff
4c =================================================================
5c CepGen common blocks for kinematics definition
6c =================================================================
7 include 'CepGen/Process/Fortran/KTBlocks.inc'
8
9c =================================================================
10c input parameters
11c =================================================================
12 double precision am_l,q_l
13 integer imethod,pdg_l
14 integer iterm11,iterm22,iterm12,itermtt
15 integer imat1,imat2
16
17c =================================================================
18c local variables
19c =================================================================
20 double precision s
21 double precision alpha1,alpha2,amt1,amt2
22 double precision pt1,pt2,eta1,eta2,dely
23 double precision pt1x,pt1y,pt2x,pt2y
24 double precision ak10,ak1x,ak1y,ak1z,ak20,ak2x,ak2y,ak2z
25 double precision beta1,beta2,x1,x2
26 double precision z1p,z1m,z2p,z2m
27 double precision q1tx,q1ty,q1z,q10,q1t2
28 double precision q2tx,q2ty,q2z,q20,q2t2
29 double precision ptsum
30 double precision that1,that2,that,uhat1,uhat2,uhat
31 double precision term1,term2,term3,term4,term5,term6,term7
32 double precision term8,term9,term10,amat2
33 double precision ak1_x,ak1_y,ak2_x,ak2_y
34 double precision eps12,eps22
35 double precision aux2_1,aux2_2
36 double precision phi10,phi102,phi11_x,phi11_y,phi112
37 double precision phi20,phi202,phi21_x,phi21_y,phi212
38 double precision phi11_dot_e,phi11_cross_e
39 double precision phi21_dot_e,phi21_cross_e
40 double precision aintegral
41
42 double precision px_plus,px_minus,py_plus,py_minus
43 double precision r1,r2
44 double precision ptdiffx,ptsumx,ptdiffy,ptsumy
45 double precision invm,invm2,s1_eff,s2_eff
46 double precision t1abs,t2abs
47 double precision inp_a,inp_b,am_a,am_b
48 double precision p1_plus, p2_minus
49 double precision amat2_1,amat2_2
50
51 double precision coupling
52
53 logical first_init
54 data first_init/.true./
55 save first_init,imethod
56 save iterm11,iterm22,iterm12,itermtt
57 save imat1,imat2
58 save pdg_l,am_l,q_l
59
60c =================================================================
61c quarks production
62c =================================================================
63 double precision t_max,amu2
64
65c =================================================================
66c at initialisation, retrieve a few user-defined parameters
67c and initialise the alpha(s) algorithm if necessary
68c =================================================================
69
70 if(first_init) then
71 call cepgen_print
72 imethod = cepgen_param_int('method', 1) ! kinematics mode
73 pdg_l = cepgen_param_int('pair', 13) ! central particles PDG
74c polarisation terms to consider in the matrix element
75 iterm11 = cepgen_param_int('term11', 1) ! LL
76 iterm22 = cepgen_param_int('term22', 1) ! TT
77 iterm12 = cepgen_param_int('term12', 1) ! LT
78 itermtt = cepgen_param_int('termtt', 1) ! TTprime
79c two terms in Wolfgang formula for off-shell gamma gamma --> l^+ l^-
80 imat1 = cepgen_param_int('mat1', 1)
81 imat2 = cepgen_param_int('mat2', 1)
82c central particles properties
83 am_l = cepgen_particle_mass(pdg_l) ! central particles mass
84 q_l = cepgen_particle_charge(pdg_l) ! central particles charge
85 if(iflux1.eq.21) then
86 if(icontri.eq.3.or.icontri.eq.4) then
87 print *,'Invalid process mode for gluon emission!'
88 stop
89 endif
90 endif
91 first_init = .false.
92 endif
93
94c =================================================================
95c start by initialising a few variables
96c =================================================================
97
98 nucl_to_ff = 0.d0
99 amat2 = 0.d0
100 eps12 = 0.d0
101 eps22 = 0.d0
102 q10 = 0.d0
103 q1z = 0.d0
104 q20 = 0.d0
105 q2z = 0.d0
106
107c =================================================================
108c go to energy of Nucleus = A*inp in order that generator puts out
109c proper momenta in LAB frame
110c =================================================================
111
112 inp_a = inp1*a_nuc1
113 am_a = am_p*a_nuc1
114 inp_b = inp2*a_nuc2
115 am_b = am_p*a_nuc2
116
117c =================================================================
118c four-momenta for incoming beams in LAB !!!!
119c =================================================================
120
121 r1 = dsqrt(1.d0+am_a**2/inp_a**2)
122 r2 = dsqrt(1.d0+am_b**2/inp_b**2)
123
124 ak10 = inp_a*r1
125 ak1x = 0.d0
126 ak1y = 0.d0
127 ak1z = inp_a
128
129 ak20 = inp_b*r2
130 ak2x = 0.d0
131 ak2y = 0.d0
132 ak2z = -inp_b
133
134 s = 4.*inp_a*inp_b*(1.d0 + r1*r2)/2.d0+(am_a**2+am_b**2)
135
136 p1_plus = (ak10+ak1z)/dsqrt(2.d0)
137 p2_minus = (ak20-ak2z)/dsqrt(2.d0)
138
139c =================================================================
140c Outgoing proton final state mass
141c =================================================================
142 if((icontri.eq.1).or.(icontri.eq.2)) am_x = am_a
143 if((icontri.eq.1).or.(icontri.eq.3)) am_y = am_b
144
145 q1tx = q1t*cos(phiq1t)
146 q1ty = q1t*sin(phiq1t)
147
148 q2tx = q2t*cos(phiq2t)
149 q2ty = q2t*sin(phiq2t)
150
151 ptsumx = q1tx+q2tx
152 ptsumy = q1ty+q2ty
153
154 ptsum = sqrt(ptsumx**2+ptsumy**2)
155
156c =================================================================
157c a window in final state transverse momentum
158c =================================================================
159
160 if(iptsum) then
161 if(ptsum.lt.ptsum_min.or.ptsum.gt.ptsum_max) return
162 endif
163
164c =================================================================
165c compute the individual central particles momentum
166c =================================================================
167
168 ptdiffx = ptdiff*cos(phiptdiff)
169 ptdiffy = ptdiff*sin(phiptdiff)
170
171 pt1x = 0.5*(ptsumx+ptdiffx)
172 pt1y = 0.5*(ptsumy+ptdiffy)
173
174 pt2x = 0.5*(ptsumx-ptdiffx)
175 pt2y = 0.5*(ptsumy-ptdiffy)
176
177 pt1 = sqrt(pt1x**2+pt1y**2)
178 pt2 = sqrt(pt2x**2+pt2y**2)
179
180 if(ipt) then
181 if(pt1.lt.pt_min.or.pt2.lt.pt_min) return
182 if(pt_max.gt.0d0.and.(pt2.gt.pt_max.or.pt2.gt.pt_max)) return
183 endif
184
185 amt1 = dsqrt(pt1**2+am_l**2)
186 amt2 = dsqrt(pt2**2+am_l**2)
187
188 invm2 = amt1**2 + amt2**2 + 2.d0*amt1*amt2*dcosh(y1-y2)
189 & -ptsum**2
190
191 invm = dsqrt(invm2)
192
193c =================================================================
194c a window in final state invariant mass
195c =================================================================
196
197 if(iinvm) then
198 if(invm.lt.invm_min) return
199 if(invm_max.gt.0d0.and.invm.gt.invm_max) return
200 endif
201
202c =================================================================
203c a window in rapidity distance
204c =================================================================
205
206 dely = dabs(y1-y2)
207 if(idely) then
208 if(dely.lt.dely_min.or.dely.gt.dely_max) return
209 endif
210
211c =================================================================
212c auxiliary quantities
213c =================================================================
214
215 alpha1 = amt1/(dsqrt(2.d0)*p1_plus)*dexp( y1)
216 alpha2 = amt2/(dsqrt(2.d0)*p1_plus)*dexp( y2)
217 beta1 = amt1/(dsqrt(2.d0)*p2_minus)*dexp(-y1)
218 beta2 = amt2/(dsqrt(2.d0)*p2_minus)*dexp(-y2)
219
220 q1t2 = q1tx**2 + q1ty**2
221 q2t2 = q2tx**2 + q2ty**2
222
223 x1 = alpha1 + alpha2
224 x2 = beta1 + beta2
225
226 z1p = alpha1/x1
227 z1m = alpha2/x1
228 z2p = beta1/x2
229 z2m = beta2/x2
230
231c -----------------------------------------------------------------
232 if(x1.gt.1.0.or.x2.gt.1.0) then
233c call CepGen_warning('Unphysical x1/x2')
234 return
235 endif
236c -----------------------------------------------------------------
237
238 s1_eff = x1*s - q1t**2
239 s2_eff = x2*s - q2t**2
240
241c -----------------------------------------------------------------
242c Additional conditions for energy-momentum conservation
243c -----------------------------------------------------------------
244 if(((icontri.eq.2).or.(icontri.eq.4))
245 1 .and.(dsqrt(s1_eff).le.(am_y+invm))) return
246 if(((icontri.eq.3).or.(icontri.eq.4))
247 1 .and.(dsqrt(s2_eff).le.(am_x+invm))) return
248
249c =================================================================
250c >>> TO THE OUTPUT COMMON BLOCK
251c =================================================================
252
253c =================================================================
254c four-momenta of the outgoing protons (or remnants)
255c =================================================================
256
257 px_plus = (1.d0-x1) * p1_plus
258 px_minus = ((a_nuc1*am_x)**2 + q1tx**2 + q1ty**2)/2.d0/px_plus
259
260 px(1) = -q1tx
261 px(2) = -q1ty
262 px(3) = (px_plus - px_minus)/dsqrt(2.d0)
263 px(4) = (px_plus + px_minus)/dsqrt(2.d0)
264
265 py_minus = (1.d0-x2) * p2_minus
266 py_plus = ((a_nuc2*am_y)**2 + q2tx**2 + q2ty**2)/2.d0/py_minus
267
268 py(1) = -q2tx
269 py(2) = -q2ty
270 py(3) = (py_plus - py_minus)/dsqrt(2.d0)
271 py(4) = (py_plus + py_minus)/dsqrt(2.d0)
272
273 q1t = dsqrt(q1t2)
274 q2t = dsqrt(q2t2)
275
276c =================================================================
277c four-momenta of the outgoing central particles
278c =================================================================
279
280 nout = 2
281
282 ipdg(1) = pdg_l
283 pc(1,1) = pt1x
284 pc(2,1) = pt1y
285 pc(3,1) = alpha1*ak1z + beta1*ak2z
286 pc(4,1) = alpha1*ak10 + beta1*ak20
287
288 ipdg(2) = -pdg_l
289 pc(1,2) = pt2x
290 pc(2,2) = pt2y
291 pc(3,2) = alpha2*ak1z + beta2*ak2z
292 pc(4,2) = alpha2*ak10 + beta2*ak20
293
294 eta1 = 0.5d0*dlog((dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2)
295 2 + amt1*dsinh(y1))/(dsqrt(amt1**2*(dcosh(y1))**2 - am_l**2)
296 3 - amt1*dsinh(y1)))
297
298 eta2 = 0.5d0*dlog((dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2)
299 2 + amt2*dsinh(y2))/(dsqrt(amt2**2*(dcosh(y2))**2 - am_l**2)
300 3 - amt2*dsinh(y2)))
301
302 if(ieta) then
303 if(eta1.lt.eta_min.or.eta1.gt.eta_max) return
304 if(eta2.lt.eta_min.or.eta2.gt.eta_max) return
305 endif
306
307c =================================================================
308c matrix element squared
309c averaged over initial spin polarizations
310c and summed over final spin polarizations
311c (--> see Wolfgang notes
312c =================================================================
313
314c =================================================================
315c Mendelstam variables
316c =================================================================
317
318 that1 = (q10-pc(4,1))**2
319 & -(q1tx-pc(1,1))**2-(q1ty-pc(2,1))**2-(q1z-pc(3,1))**2
320 uhat1 = (q10-pc(4,2))**2
321 & -(q1tx-pc(1,2))**2-(q1ty-pc(2,2))**2-(q1z-pc(3,2))**2
322 that2 = (q20-pc(4,2))**2
323 & -(q2tx-pc(1,2))**2-(q2ty-pc(2,2))**2-(q2z-pc(3,2))**2
324 uhat2 = (q20-pc(4,1))**2
325 & -(q2tx-pc(1,1))**2-(q2ty-pc(2,1))**2-(q2z-pc(3,1))**2
326
327 that = (that1+that2)/2.d0
328 uhat = (uhat1+uhat2)/2.d0
329
330c =================================================================
331c matrix elements
332c =================================================================
333c How matrix element is calculated
334c imethod = 0: on-shell formula
335c imethod = 1: off-shell formula
336c =================================================================
337 if(imethod.eq.0) then
338c =================================================================
339c on-shell formula for M^2
340c =================================================================
341 term1 = 6.d0*am_l**8
342 term2 = -3.d0*am_l**4*that**2
343 term3 = -14.d0*am_l**4*that*uhat
344 term4 = -3.d0*am_l**4*uhat**2
345 term5 = am_l**2*that**3
346 term6 = 7.d0*am_l**2*that**2*uhat
347 term7 = 7.d0*am_l**2*that*uhat**2
348 term8 = am_l**2*uhat**3
349 term9 = -that**3*uhat
350 term10 = -that*uhat**3
351
352 amat2 = -2.d0*(term1+term2+term3+term4+term5
353 2 +term6+term7+term8+term9+term10)
354 3 / ( (am_l**2-that)**2 * (am_l**2-uhat)**2 )
355
356 elseif(imethod.eq.1)then
357c =================================================================
358c Wolfgang formulae
359c =================================================================
360
361 ak1_x = z1m*pt1x-z1p*pt2x
362 ak1_y = z1m*pt1y-z1p*pt2y
363
364 ak2_x = z2m*pt1x-z2p*pt2x
365 ak2_y = z2m*pt1y-z2p*pt2y
366
367 !FIXME FIXME FIXME am_p or am_A/B???
368 t1abs = (q1t2+x1*(am_x**2-am_p**2)+x1**2*am_p**2)/(1.d0-x1)
369 t2abs = (q2t2+x2*(am_y**2-am_p**2)+x2**2*am_p**2)/(1.d0-x2)
370
371 eps12 = am_l**2 + z1p*z1m*t1abs
372 eps22 = am_l**2 + z2p*z2m*t2abs
373
374 phi10 = 1.d0/((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
375 2 - 1.d0/((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
376 phi11_x = (ak1_x+z1p*q2tx)/
377 2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
378 3 - (ak1_x-z1m*q2tx)/
379 4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
380 phi11_y = (ak1_y+z1p*q2ty)/
381 2 ((ak1_x+z1p*q2tx)**2+(ak1_y+z1p*q2ty)**2+eps12)
382 3 - (ak1_y-z1m*q2ty)/
383 4 ((ak1_x-z1m*q2tx)**2+(ak1_y-z1m*q2ty)**2+eps12)
384
385 phi102 = phi10*phi10
386 phi112 = phi11_x**2+phi11_y**2
387
388 phi20 = 1.d0/((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
389 2 - 1.d0/((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
390
391 phi21_x = (ak2_x+z2p*q1tx)/
392 2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
393 3 - (ak2_x-z2m*q1tx)/
394 4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
395 phi21_y = (ak2_y+z2p*q1ty)/
396 2 ((ak2_x+z2p*q1tx)**2+(ak2_y+z2p*q1ty)**2+eps22)
397 3 - (ak2_y-z2m*q1ty)/
398 4 ((ak2_x-z2m*q1tx)**2+(ak2_y-z2m*q1ty)**2+eps22)
399
400 phi202 = phi20*phi20
401 phi212 = phi21_x**2+phi21_y**2
402
403 phi11_dot_e = (phi11_x*q1tx + phi11_y*q1ty)/dsqrt(q1t2)
404 phi11_cross_e = (phi11_x*q1ty - phi11_y*q1tx)/dsqrt(q1t2)
405
406 phi21_dot_e = (phi21_x*q2tx +phi21_y*q2ty)/dsqrt(q2t2)
407 phi21_cross_e = (phi21_x*q2ty -phi21_y*q2tx)/dsqrt(q2t2)
408
409 aux2_1 = iterm11*(am_l**2+4.d0*z1p**2*z1m**2*t1abs)*phi102
410 1 +iterm22*( (z1p**2 + z1m**2)*(phi11_dot_e**2 +
411 2 phi11_cross_e**2) )
412 3 + itermtt*( phi11_cross_e**2 - phi11_dot_e**2)
413 4 - iterm12*4.d0*z1p*z1m*(z1p-z1m)*phi10
414 5 *(q1tx*phi11_x+q1ty*phi11_y)
415
416 aux2_2 = iterm11*(am_l**2+4.d0*z2p**2*z2m**2*t2abs)*phi202
417 1 +iterm22*( (z2p**2 + z2m**2)*(phi21_dot_e**2 +
418 2 phi21_cross_e**2) )
419 3 + itermtt*( phi21_cross_e**2 - phi21_dot_e**2)
420 4 - iterm12*4.d0*z2p*z2m*(z2p-z2m)*phi20
421 5 *(q2tx*phi21_x+q2ty*phi21_y)
422
423c =================================================================
424c convention of matrix element as in our kt-factorization
425c for heavy flavours
426c =================================================================
427
428 amat2_1 = (x1*x2*s)**2 * aux2_1 * 2.*z1p*z1m*q1t2 / (q1t2*q2t2)
429 amat2_2 = (x1*x2*s)**2 * aux2_2 * 2.*z2p*z2m*q2t2 / (q1t2*q2t2)
430
431c =================================================================
432c symmetrization
433c =================================================================
434
435 amat2 = (imat1*amat2_1 + imat2*amat2_2)/2.d0
436
437 endif
438
439 coupling = 1.d0
440c =================================================================
441c first parton coupling
442c =================================================================
443 t_max = max(amt1,amt2)**2
444 amu2 = max(eps12,t_max)
445 if(iflux1.eq.21) then ! at least one gluon exchanged
446 coupling = coupling * 4.d0*pi*cepgen_alphas(dsqrt(amu2))/2.d0
447 else ! photon exchange
448 coupling = coupling * 4.d0*pi*cepgen_alphaem(dsqrt(amu2))*q_l**2
449 endif
450c =================================================================
451c second parton coupling
452c =================================================================
453 coupling = coupling * 4.d0*pi*cepgen_alphaem(dsqrt(amu2))
454 & * q_l**2 ! photon exchange
455
456c =================================================================
457c factor 2.*pi below from integration over phi_sum
458c factor 1/4 below from jacobian of transformations
459c factors 1/pi and 1/pi due to integration
460c over d^2 kappa_1 d^2 kappa_2 instead d kappa_1^2 d kappa_2^2
461c =================================================================
462
463 aintegral = (2.d0*pi)*1.d0/(16.d0*pi**2*(x1*x2*s)**2) * amat2
464 & * coupling
465 & * (1.d0/4.d0) * units
466 & * 0.5d0 * 4.0d0 / (4.d0*pi)
467
468c *****************************************************************
469c =================================================================
470 nucl_to_ff = aintegral*q1t*q2t*ptdiff
471c =================================================================
472c *****************************************************************
473
474 return
475 end
cepgen_print
subroutine cepgen_print
Definition cepgen_print.f:18
nucl_to_ff
double precision function nucl_to_ff()
Definition nucl_to_ff.f:2