11 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
16 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
17 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
22 DOUBLE PRECISION minq2, rccorr,
pyth_xsec, temp, etheta
23 DOUBLE PRECISION ppt, ppx, ppy
24 include
"mcRadCor.inc"
33 dimension pms(2),
xmin(2),
xmax(2),q2min(2),q2max(2),pmc(3),
35 SAVE pms,
xmin,
xmax,q2min,q2max,pmc,
x,q2,theta,
phi,
pt,w2min,
49 pmc(3)=
vint(302)-pms(1)-pms(2)
50 w2min=
max(ckin(77),2d0*ckin(3),2d0*ckin(5))**2
54 pmc(
i)=
vint(302)+pms(
i)-pms(3-
i)
55 IF(
mint(140+
i).NE.0)
THEN
57 xmax(
i)=min(ckin(60+2*
i),1d0-2d0*
vint(301)*sqrt(pms(
i))/
62 & (
ymin*pmc(3)-ckin(64+2*
i))/pmc(
i))
64 themin=
max(ckin(67+2*
i),0d0)
65 themax=min(ckin(68+2*
i),paru(1))
66 IF(ckin(68+2*
i).LT.0d0) themax=paru(1)
69 & 2d0*pms(
i)/(1d0-
xmax(
i)))*sin(themin/2d0)**2,0d0)
72 & 2d0*pms(
i)/(1d0-
xmin(
i)))*sin(themax/2d0)**2
73 IF(ckin(64+2*
i).GT.0d0) q2max(
i)=min(ckin(64+2*
i),q2max(
i))
75 IF(
mint(143-
i).EQ.0)
THEN
78 & (ckin(78)**2-pms(3-
i))/pmc(
i))
84 IF(
mint(141).NE.0.AND.
mint(142).NE.0)
THEN
85 IF(ckin(78).GT.0d0)
xmax(1)=min(
xmax(1),
86 & (ckin(78)**2+pmc(3)-pmc(2)*
xmin(2))/pmc(1))
87 IF(ckin(78).GT.0d0)
xmax(2)=min(
xmax(2),
88 & (ckin(78)**2+pmc(3)-pmc(1)*
xmin(1))/pmc(2))
89 IF(iabs(
mint(141)).NE.iabs(
mint(142)))
THEN
91 & pms(1)-pms(2))/(pmc(2)*
xmax(2)+pms(1)-pms(2)))/pmc(1))
93 & pms(1)-pms(2))/(pmc(1)*
xmax(1)+pms(2)-pms(1)))/pmc(2))
101 ELSEIF(igaga.EQ.2)
THEN
108 IF(
mint(141).NE.0.AND.
mint(142).EQ.0)
THEN
110 ELSEIF(
mint(141).EQ.0.AND.
mint(142).NE.0)
THEN
112 ELSEIF(isub.GE.137.AND.isub.LE.140)
THEN
113 vint(2)=
max(ckin(77)**2,12d0*
max(ckin(3),ckin(5))**2)
114 IF(ckin(78).GT.0d0)
vint(2)=min(
vint(2),ckin(78)**2)
117 IF(ckin(78).GT.0d0)
vint(2)=min(
vint(2),ckin(78)**2)
125 IF(
mint(140+
i).NE.0)
THEN
126 IF(
mstp(199).EQ.1)
then
127 wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
128 & (
log(
real(mcset_ymax/mcset_ymin)))*
129 & (
log(
real(mcset_q2max/mcset_q2min)))
131 wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
134 IF(isub.EQ.99.AND.
mint(106+
i).EQ.4.AND.
mint(109-
i).EQ.3)
136 q2init=5d0+q2min(3-
i)
137 ELSEIF(isub.EQ.99.AND.
mint(106+
i).EQ.4)
THEN
138 q2init=pmas(
pycomp(113),1)**2+q2min(3-
i)
139 ELSEIF(isub.EQ.132.OR.isub.EQ.134.OR.isub.EQ.136)
THEN
140 q2init=
max(ckin(1),2d0*ckin(3),2d0*ckin(5))**2/3d0
141 ELSEIF((isub.EQ.138.AND.
i.EQ.2).OR.
142 & (isub.EQ.139.AND.
i.EQ.1))
THEN
144 ELSEIF(isub.EQ.140)
THEN
149 vint(2+
i)=-sqrt(
max(q2min(
i),min(q2max(
i),q2init)))
150 IF(
mstp(14).EQ.0.OR.(isub.GE.131.AND.isub.LE.140))
158 IF(
mstp(82).LE.1)
THEN
170 ELSEIF(igaga.EQ.3)
THEN
179 IF(
mint(140+
i).NE.0)
THEN
182 q2(
i)=q2min(
i)*(q2max(
i)/q2min(
i))**
pyr(0)
184 IF(q2(
i).LT.
x(
i)**2*pms(
i)/(1d0-
x(
i))) goto 120
185 IF(q2(
i).GT.(1d0-
x(
i))*(
vint(302)-2d0*pms(3-
i))-
186 & (2d0-
x(
i)**2)*pms(
i)/(1d0-
x(
i))) goto 120
188 y(
i)=(pmc(
i)*
x(
i)+q2(
i))/pmc(3)
189 IF(
y(
i).LT.ckin(71+2*
i).OR.
y(
i).GT.ckin(72+2*
i)) goto 120
190 rat=((1d0-
x(
i))*q2(
i)-
x(
i)**2*pms(
i))/
191 & ((1d0-
x(
i))**2*(
vint(302)-2d0*pms(3-
i)-2d0*pms(
i)))
192 theta(
i)=2d0*asin(sqrt(
max(0d0,min(1d0,rat))))
193 IF(theta(
i).LT.ckin(67+2*
i)) goto 120
194 IF(ckin(68+2*
i).GT.0d0.AND.theta(
i).GT.ckin(68+2*
i))
199 pt(
i)=sqrt(((1d0-
x(
i))*pmc(
i))**2/(4d0*
vint(302))-
200 & pms(
i))*sin(theta(
i))
219 IF(
mint(141).NE.0.AND.
mint(142).NE.0)
THEN
220 w2=-q2(1)-q2(2)+0.5d0*
x(1)*pmc(1)*
x(2)*pmc(2)/
vint(302)-
222 & sqrt((0.5d0*
x(1)*pmc(1)/
vint(301))**2+q2(1)-
pt(1)**2)*
223 & sqrt((0.5d0*
x(2)*pmc(2)/
vint(301))**2+q2(2)-
pt(2)**2)
224 IF(w2.LT.w2min) goto 120
225 IF(ckin(78).GT.0d0.AND.w2.GT.ckin(78)**2) goto 120
228 ELSEIF(
mint(141).NE.0)
THEN
229 w2=(
vint(302)+pms(1))*
x(1)+pms(2)*(1d0-
x(1))
232 ELSEIF(
mint(142).NE.0)
THEN
233 w2=(
vint(302)+pms(2))*
x(2)+pms(1)*(1d0-
x(2))
243 vint(293)=0.5d0*sqrt((w2-pms1-pms2)**2-4d0*pms1*pms2)/
vint(1)
244 vint(294)=0.5d0*(w2+pms1-pms2)/
vint(1)
245 vint(295)=
sign(sqrt(abs(pms1)),pms1)
249 vint(299)=0.5d0*(w2+pms2-pms1)/
vint(1)
250 vint(300)=
sign(sqrt(abs(pms2)),pms2)
256 IF(
mint(140+
i).NE.0)
THEN
257 wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
259 IF(
mstp(16).EQ.0)
THEN
262 wtgaga=wtgaga*
x(
i)/
y(
i)
265 IF(isub.EQ.132.OR.isub.EQ.134.OR.isub.EQ.136)
THEN
266 IF((
mint(11).EQ.22).and.
267 & (
mint(12).EQ.2212.or.
mint(12).EQ.2112))
THEN
268 wtgaga=wtgaga*(1d0/(1d0+(q2(
i)/xy**2/ebeamenucl**2))*
269 & (1d0-xy-(q2(
i)/4d0/ebeamenucl**2)))/
270 & q2(
i)/xy**2/ebeamenucl*
271 & (ebeamenucl*xy-q2(
i)/2d0/massp)*xy*q2(
i)
273 wtgaga=wtgaga*(1d0-xy)
275 ELSEIF(
i.EQ.1.AND.(isub.EQ.139.OR.isub.EQ.140))
THEN
276 wtgaga=wtgaga*(1d0-xy)
277 ELSEIF(
i.EQ.2.AND.(isub.EQ.138.OR.isub.EQ.140))
THEN
278 wtgaga=wtgaga*(1d0-xy)
279 ELSEIF((
mint(11).EQ.22).and.
280 & (
mint(12).EQ.2212.or.
mint(12).EQ.2112))
THEN
281 wtgaga=wtgaga*(0.5d0*((ebeamenucl*xy-q2(
i)/2d0/
282 & massp)/q2(
i)/xy**2/ebeamenucl*
283 & (xy**2*(1d0-(2d0*masse**2/q2(
i)))+
284 & (2d0/(1d0+(q2(
i)/xy**2/ebeamenucl**2)))*
285 & (1d0-xy-(q2(
i)/4d0/ebeamenucl**2))))*
288 wtgaga=wtgaga*(0.5d0*(1d0+(1d0-xy)**2)-
289 & pms(
i)*xy**2/q2(
i))
298 IF(
mstp(82).LE.1)
THEN
310 ELSEIF(igaga.EQ.5)
THEN
311 write(*,*).EQ.
"In pygaga with IGAGA5"
320 IF(
mint(140+
i).NE.0)
THEN
323 geny=mcset_ymin*(mcset_ymax/mcset_ymin)**
pyr(0)
324 genq2=mcset_q2min*(mcset_q2max/mcset_q2min)**
pyr(0)
325 genx = genq2/geny/(4.*ebeam*pbeam)
326 genw2 = massp**2.+(genq2*(1./genx-1.))
327 gennu=geny*ebeamenucl
328 geneprim = ebeamenucl - gennu
329 write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
331 minq2 = pms(1) * geny**2. / (1.- geny)
332 if (genq2.lt.minq2)
then
336 if (genq2.gt.(2.*gennu*massp))
then
340 temp=(genq2/(2.*ebeamenucl*geneprim))-1.
341 if ((temp.le.1.00).and.(temp.ge.-1.0))
then
343 genthe =
sngl(etheta)
344 genphi=paru(2)*
pyr(0)
353 if ((genw2.lt.ckin(77)**2).or.
354 & (ckin(78).gt.0.and.genw2.gt.ckin(78)**2))
then
357 write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
360 122
if (qedrad.eq.1)
then
361 write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
364 if (qedrad.eq.0)
then
367 elseif ((mcradcor_ebrems.eq.mcradcor_ebrems).and.
368 & (mcradcor_thetabrems.eq.mcradcor_thetabrems))
then
369 y(
i)=dble(mcradcor_nutrue)/dble(mcset_enebeam)
370 q2(
i)=dble(mcradcor_q2true)
371 write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
372 write(*,*) mcradcor_nutrue, mcset_enebeam, mcradcor_q2true
374 write(*,*)
"I go to 122 again"
375 write(*,*) mcradcor_thetabrems,mcradcor_ebrems,
381 x(
i)=((pmc(3)*
y(
i))-q2(
i))/pmc(
i)
385 IF (qedrad.ne.0)
then
392 IF(dble(mcradcor_w2true).LT.
393 & (ckin(77)**2-1.d-4*abs(ckin(77)**2)))
THEN
399 IF(ntries.ge.20) goto 121
405 if (qedrad.eq.1)
then
406 call
mkf2(genq2,genx,
407 + mcset_tara,mcset_tarz,py6f2,py6f1)
408 sigobs=
pyth_xsec(geny,genx,genq2,py6f1, py6f2)
409 IF(mcradcor_ebrems.eq.0)
THEN
410 IF (sig1g.gt.0.d0)
then
411 rccorr=(tbor+tine)/sig1g/(dble(
mint(199))+1.0d0)
415 ELSEIF(mcradcor_ebrems.gt.0)
THEN
416 call
mkf2(dble(mcradcor_q2true),dble(mcradcor_xtrue),
417 + mcset_tara,mcset_tarz,py6f2,py6f1)
418 sigtrue=
pyth_xsec(dble( mcradcor_ytrue), dble(mcradcor_xtrue),
419 + dble(mcradcor_q2true), py6f1, py6f2)
420 IF ((sig1g.gt.0.d0).and.(sigtrue.gt.0.d0))
then
421 rccorr=(tbor+tine)/sig1g*sigobs/sigtrue/
422 & (dble(
mint(199))+1.0d0)
432 IF(q2(
i).LT.
x(
i)**2*pms(
i)/(1d0-
x(
i)))
then
435 IF(q2(
i).GT.(1d0-
x(
i))*(
vint(302)-2d0*pms(3-
i))-
436 & (2d0-
x(
i)**2)*pms(
i)/(1d0-
x(
i)))
THEN
440 IF(
y(
i).LT.ckin(71+2*
i).OR.
y(
i).GT.ckin(72+2*
i))
THEN
443 rat=((1d0-
x(
i))*q2(
i)-
x(
i)**2*pms(
i))/
444 & ((1d0-
x(
i))**2*(
vint(302)-2d0*pms(3-
i)-2d0*pms(
i)))
445 theta(
i)=2d0*asin(sqrt(
max(0d0,min(1d0,rat))))
446 IF(theta(
i).LT.ckin(67+2*
i))
THEN
449 IF(ckin(68+2*
i).GT.0d0.AND.theta(
i).GT.ckin(68+2*
i))
455 temp=((1d0-
x(
i))*pmc(
i))**2/(4d0*
vint(302))-pms(
i)
456 pt(
i)=(sqrt(temp))*sin(theta(
i))
458 IF ((qedrad.ne.0).and.(mcradcor_ebrems.gt.0))
then
459 emom=sqrt(geneprim**2-masse**2)
460 phi(
i)=atan2((emom*ppy+dplabg(2)),(emom*ppx+dplabg(1)))
461 IF (
phi(
i).lt.0)
THEN
482 IF(
mint(141).NE.0.AND.
mint(142).NE.0)
THEN
483 w2=-q2(1)-q2(2)+0.5d0*
x(1)*pmc(1)*
x(2)*pmc(2)/
vint(302)-
485 & sqrt((0.5d0*
x(1)*pmc(1)/
vint(301))**2+q2(1)-
pt(1)**2)*
486 & sqrt((0.5d0*
x(2)*pmc(2)/
vint(301))**2+q2(2)-
pt(2)**2)
490 IF(ckin(78).GT.0d0.AND.w2.GT.ckin(78)**2) goto 121
493 ELSEIF(
mint(141).NE.0)
THEN
494 w2=(
vint(302)+pms(1))*
x(1)+pms(2)*(1d0-
x(1))
497 ELSEIF(
mint(142).NE.0)
THEN
498 w2=(
vint(302)+pms(2))*
x(2)+pms(1)*(1d0-
x(2))
508 vint(293)=0.5d0*sqrt((w2-pms1-pms2)**2-4d0*pms1*pms2)/
vint(1)
509 vint(294)=0.5d0*(w2+pms1-pms2)/
vint(1)
510 vint(295)=
sign(sqrt(abs(pms1)),pms1)
514 vint(299)=0.5d0*(w2+pms2-pms1)/
vint(1)
515 vint(300)=
sign(sqrt(abs(pms2)),pms2)
521 IF(
mint(140+
i).NE.0)
THEN
522 wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
523 & (
log(
real(mcset_ymax))-
log(
real(mcset_ymin)))*
524 & (
log(
real(mcset_q2max))-
log(
real(mcset_q2min)))
526 IF(isub.EQ.132.OR.isub.EQ.134.OR.isub.EQ.136)
THEN
527 IF((
mint(11).EQ.22).and.
528 & (
mint(12).EQ.2212.or.
mint(12).EQ.2112))
THEN
529 xxy=xy*ebeamenucl/ebeamenucl
530 wtgaga=wtgaga*(1d0/(1d0+(q2(
i)/xxy**2/ebeamenucl**2))*
531 & (1d0-xxy-(q2(
i)/4d0/ebeamenucl**2)))/
532 & q2(
i)/xxy**2/ebeamenucl*
533 & (ebeamenucl*xxy-q2(
i)/2d0/massp)*xxy*q2(
i)
535 wtgaga=wtgaga*(1d0-xy)
537 ELSEIF(
i.EQ.1.AND.(isub.EQ.139.OR.isub.EQ.140))
THEN
538 wtgaga=wtgaga*(1d0-xy)
539 ELSEIF(
i.EQ.2.AND.(isub.EQ.138.OR.isub.EQ.140))
THEN
540 wtgaga=wtgaga*(1d0-xy)
541 ELSEIF((
mint(11).EQ.22).and.
542 & (
mint(12).EQ.2212.or.
mint(12).EQ.2112))
THEN
543 xxy=xy*ebeamenucl/ebeamenucl
544 wtgaga=wtgaga*(0.5d0*((ebeamenucl*xxy-q2(
i)/2d0/
545 & massp)/q2(
i)/xxy**2/ebeamenucl*
546 & (xxy**2*(1d0-(2d0*masse**2/q2(
i)))+
547 & (2d0/(1d0+(q2(
i)/xxy**2/ebeamenucl**2)))*
548 & (1d0-xxy-(q2(
i)/4d0/ebeamenucl**2))))*xxy*q2(
i))
550 wtgaga=wtgaga*(0.5d0*(1d0+(1d0-xy)**2)-
551 & pms(
i)*xy**2/q2(
i))
560 IF(
mstp(82).LE.1)
THEN
570 ELSEIF(igaga.EQ.4)
THEN
574 IF(
mint(141).NE.0.AND.
mint(142).NE.0) move=4
577 IF(
k(
i,1).EQ.21)
THEN
584 &
k(
i+move,3)=
k(
i,3)+move
586 &
k(
i+move,4)=
k(
i,4)+move
588 &
k(
i+move,5)=
k(
i,5)+move
606 IF(
mint(140+
i).NE.0)
THEN
613 p(
mint(83)+
i,3)=0.5d0*sqrt((pmc(3)**2-4d0*pms(1)*pms(2))/
614 &
vint(302))*(-1d0)**(
i+1)
619 IF(
mint(141).NE.0.AND.
mint(142).NE.0)
THEN
628 ELSEIF(
mint(141).NE.0)
THEN
635 ELSEIF(
mint(142).NE.0)
THEN
646 IF(
mint(140+
i).NE.0)
THEN
647 lsc=
mint(83)+min(
i+2,move)
653 p(lsc,3)=sqrt(
p(lsc,4)**2-pms(
i))*cos(theta(
i))*
661 IF(
mint(141).NE.0)
THEN
671 IF(
mint(142).NE.0)
THEN
688 CALL
pyrobo(
n+1,
n+2,0d0,-bphi,0d0,0d0,0d0)
695 IF(
mint(140+
i).NE.0)
THEN
696 lsc=
mint(83)+min(
i+2,move)