1 jones 1.1 real*8 function epc_func(EBEAM,Z1,N1,FERMI1,MPI,PART,P,THP)
2
3 CGAW PROGRAM EPC
4
5 C VAX VERSION
6 C ELECTROPRODUCTION YIELDS OF NUCLEONS AND PIONS
7 C WRITTEN BY J.S. OCONNELL AND J.W. LIGHTBODY, JR.
8 C NATIONAL BUREAU OF STANDARDS
9 C APRIL 1988
10 C TRANSVERSE SCALING REGION ADDED
11
12 C Modified by OARA to plot like older EPC's
13 C Modified slightly by Glen Warren to compile under Linux (g77) Sept. 02
14
15 IMPLICIT REAL*8 (A-H,O-Z)
16 REAL*8 MP,MPI0,MN
17 CHARACTER*1 MPI,SCAL,FERMI,FERMI1
18 CHARACTER*3 PART
19
20 COMMON/SG/IA
21 COMMON/QD/QDF
22 jones 1.1 COMMON/DEL/IP
23 COMMON /FER/ FERMI,SCAL
24 COMMON/M/Z,N
25 COMMON/KF/E1
26 DATA AM/938.28/,AMP/139.6/,MPI0/135.9/,MN/939.56/
27
28 SCAL = 'N'
29
30 FERMI = FERMI1
31 Z = Z1
32 N = N1
33 E1 = EBEAM
34
35 PI=ACOS(-1.0D0)
36 IA=Z+N
37
38 C 'AN' is effective number of nucleons for one pion production
39
40 IF(PART.EQ.'P')THEN
41 AN=N/3.d0+2.d0*Z/3.d0
42 MP = AM
43 jones 1.1 IP=1
44 ELSEIF(PART.EQ.'N')THEN
45 AN=Z/3.d0+2.d0*N/3.d0
46 MP = MN
47 IP=-1
48 ELSEIF(PART.EQ.'PI+')THEN
49 AN=Z/3.d0
50 MP = AMP
51 IP=2
52 ELSEIF(PART.EQ.'PI-')THEN
53 AN=N/3.d0
54 MP = AMP
55 IP=2
56 ELSEIF(PART.EQ.'PI0')THEN
57 AN=2.d0*(N+Z)/3.d0
58 MP = MPI0
59 IP=0
60 ELSE
61 CALL EXIT
62 ENDIF
63
64 jones 1.1 IF(ABS(IP).EQ.1)THEN
65 IF(IA.EQ.1)THEN
66 DLF = 0
67 ELSE
68 IF(IA.GT.1.AND.IA.LT.5)THEN
69 DLF=IA
70 ELSE
71 DLF=7.
72 ENDIF
73
74 AL=DLF ! use default setting for Levinger factor
75 QDF=AL*N*Z/FLOAT(IA)
76 END IF
77 ENDIF
78
79 TH=THP*PI/180.
80
81 IF(ABS(IP).EQ.1)THEN
82 E=SQRT(P**2+AM**2)
83 TP=P**2/(E+AM)
84 AJ=P/E ! Converts cross section from 1/MeV to 1/MeV/c
85 jones 1.1 IF(IA.EQ.1)THEN
86 D2QD=0.
87 D2QF=0.
88 ELSEIF(IA.GT.1)THEN
89 CALL DEP_OR(E1,TP,TH,IP,D2QD)
90 IF(SCAL.NE.'Y') THEN
91 D2QD=D2QD*AJ
92 END IF
93 IF(SCAL.NE.'Y') THEN
94 CALL EP_OR(E1,TP,TH,D2QF)
95 D2QF=D2QF*AJ
96 END IF
97 ENDIF
98 ELSEIF(ABS(IP).EQ.2.OR.IP.EQ.0)THEN
99 E=SQRT(P**2+AMP**2)
100 TP=P**2/(E+AMP)
101 AJ=P/E
102 D2QD=0.
103 D2QF=0.
104 ELSE
105 CALL EXIT
106 jones 1.1 ENDIF
107
108 CALL DELTA_OR(E1,TP,TH,D2DEL)
109 IF(SCAL.EQ.'Y') THEN
110 D2DEL=AN*D2DEL
111 ELSE
112 D2DEL=AN*D2DEL*AJ
113 END IF
114
115 D2SC = 0.d0
116 IF(MPI.EQ.'Y') THEN
117 IF(ABS(IP).EQ.2.OR.IP.EQ.0)THEN
118 CALL S2PI(2,E1,TP,TH,D2SC1)
119 CALL S2PI(-2,E1,TP,TH,D2SC2)
120 IF(PART.EQ.'PI+')THEN
121 D2SC=Z*D2SC1+N*D2SC2
122 ELSEIF(PART.EQ.'PI-')THEN
123 D2SC=N*D2SC1+Z*D2SC2
124 ELSEIF(PART.EQ.'PI0')THEN
125 D2SC=IA*D2SC1
126 ELSE
127 jones 1.1 CALL EXIT
128 ENDIF
129 ELSEIF(ABS(IP).EQ.1)THEN
130 CALL S2PI(1,E1,TP,TH,D2SC1)
131 D2SC=IA*D2SC1
132 ENDIF
133 END IF
134
135 TOTAL=D2QD+D2QF+D2DEL+D2SC
136
137 IF(SCAL.EQ.'Y') THEN
138 TOTALE=TOTAL*P*1.0D-6
139 ELSE
140 TOTALE=TOTAL/AJ
141 END IF
142
143 epc_func = total*1.d6 ! nanobarns / GeV/c / sr
144 return
145 END
146
147 C------------------------------------------------------------------------------
148 jones 1.1
149 *VTP_OR
150 SUBROUTINE VTP_OR(AMT,AM1,EI,W0,TP,TH,GN)
151 C TIATOR-WRIGHT VIRTUAL PHOTON SPECTRUM
152 C PHYS. REV. C26,2349(1982) AND NUC. PHYS. A379,407(1982)
153 IMPLICIT REAL*8 (A-H,O-Z)
154 DATA AME/.511/
155 PI=ACOS(-1.0D0)
156 EF0=EI-W0
157 AKI=SQRT(EI**2-AME**2)
158 AKF0=SQRT(EF0**2-AME**2)
159 AKP=SQRT(TP**2+2.*AM1*TP)
160 EP=TP+AM1
161 AR=EI+AMT-EP
162 BR=EF0*(AKP*COS(TH)-AKI)/AKF0
163 C BRP=(AKF0/EF0)**2*BR
164 A=AME**2-EI*EF0
165 B=AKI*AKF0
166 D=-AME**2*BR*(EI/EF0-1.)/AR
167 AP=A-D
168 BP=B+D
169 jones 1.1 AN1=1./137./2./PI*W0**2/AKI**2
170 ! APB=-AME**2*(AKI-AKF0)**2/(AME**2+EI*EF0+AKI*AKF0)
171 ! PRINT *,'APB',APB
172 APB = A + B ! New
173 ! PRINT *,'A+B',APB
174 AN1=AN1*B/BP*(AR+BR)/(AR-AP/BP*BR)
175 AN2=1.-2.*A/W0**2
176 AN4=((AP-BP)*(AR+BR)/APB/(AR-BR))
177 IF(AN4.LE.0.)GO TO 1
178 AN2=AN2*LOG(AN4)
179 AN3=-4.*B/W0**2
180 ANE=AN1*(AN2+AN3)
181 D0=AMT+EI-EP+EF0/AKF0*(AKP*COS(TH)-AKI)
182 R=(AMT+W0-EP/AKP*W0*COS(TH))/D0
183 GN=ANE*R/W0
184 IF(GN.LT.0.)GN=0.
185 ! PRINT *,'TP,W0,R,ANE,GN',TP,W0,R,ANE,GN
186 RETURN
187 1 GN=0.
188 RETURN
189 END
190 jones 1.1 *DEP
191 SUBROUTINE DEP_OR(E1,TP,TH,IP,D2QD)
192 C QUASI-DEUTERON CROSS SECTION
193 IMPLICIT REAL*8 (A-H,O-Z)
194 CHARACTER*1 SCAL,FERMI
195 COMMON/SG/IA
196 COMMON/QD/QDF
197 COMMON /FER/ FERMI,SCAL
198 DATA AM/939./,AMD/1876./
199 IF(IA.EQ.1)GOTO 1
200 PN=SQRT(TP**2+2.*AM*TP)
201 EP = TP + AM
202 CALL KINE_OR(AMD,AM,AM,PN,TH,W0,THC)
203 IF(W0.GE.E1)GO TO 1
204 IF(W0.LE.0.)GO TO 1
205 W0G=W0/1000.
206 CALL SIGD_OR(W0G,THC,IP,DSQD)
207 CALL PART_OR(AMD,AM,AM,PN,TH,AJT,AJW)
208 IF(SCAL.NE.'Y') THEN
209 CALL VTP_OR(AMD,AM,E1,W0,TP,TH,PHI)
210 C CROSS SECTION IN UB/MEV-SR
211 jones 1.1 C OARA: AJW gives (w+-w-)/(p+-p-) so it is dw/dp, not dw/dt!
212 C Cross section is then in ub/((MeV/c).sr) directly
213 ! PRINT *,'QDF,PHI,DSQD,AJW,AJT',QDF,PHI,DSQD,AJW,AJT
214 D2QD=QDF*PHI*DSQD*AJW*AJT*EP/PN
215 ELSE
216 BREM=1./W0
217 C CROSS SECTION IN UB/MEV-SR-Q
218 C Cross section was then in ub/(Q.(MeV/c).sr) directly
219 D2QD = QDF*BREM*DSQD*AJW*AJT*EP/PN ! From GPC
220 ! PRINT *,'QDF,BREM,DSQD,AJW,AJT',QDF,BREM,DSQD,AJW,AJT
221 ! D2QD = QDF*BREM*DSQD*AJW*AJT One step conversion: next line
222 ! D2QD = D2QD*EP*1.0D6/PN**2 ! Convert to ub.c/(Q.(GeV/c)^2.sr)
223 D2QD = D2QD*1.0D6/PN ! Convert to units ub.c/(Q.(GeV/c)^2.sr)
224 END IF
225 RETURN
226 1 D2QD=0.
227 RETURN
228 END
229 *SIGD_OR
230 SUBROUTINE SIGD_OR(E,TH,IP,DSQD)
231 C DEUTERON CROSS SECTION
232 jones 1.1 C BASED ON FIT OF THORLACIUS & FEARING
233 C PHYS. REV. C33,1830(1986)
234 C PHOTON ENERGY RANGE 10 - 625 MEV
235 C E[GEV] IN LAB SYSTEM
236 C TH[RAD] & DSQD[UB/SR] IN CENTER-OF-MOMENTUM SYSTEM
237 IMPLICIT REAL*8 (A-H,O-Z)
238 DIMENSION C0(8),C1(4),C2(4),C3(4),C4(4)
239 DIMENSION A(0:4),B(4,4)
240 DATA C0/2.61E2,-1.10E2,2.46E1,-1.71E1,5.76E0,
241 # -2.05E0,2.67E-1,1.13E2/
242 DATA C1/1.68E1,-4.66E1,2.56E0,-4.72E0/
243 DATA C2/-2.03E2,-8.12E1,-4.05E0,-5.99E0/
244 DATA C3/-1.77E1,-3.74E1,-5.07E-1,-5.40E0/
245 DATA C4/-2.05E0,-7.05E0,9.40E-1,-2.05E0/
246 X=COS(TH)
247 IF(E.LE.625.)THEN
248 C TEST FOR NEUTRON
249 X=IP*X
250 C COEFICIENTS
251 A(0)=C0(1)*EXP(C0(2)*E)+C0(3)*EXP(C0(4)*E)
252 A(0)=A(0)+(C0(5)+C0(6)*E)/(1.+C0(8)*(E-C0(7))**2)
253 jones 1.1 DSQD=A(0)*PL(0,X)
254 DO 2 L=1,4
255 B(1,L)=C1(L)
256 B(2,L)=C2(L)
257 B(3,L)=C3(L)
258 2 B(4,L)=C4(L)
259 DO 1 L=1,4
260 A(L)=B(L,1)*EXP(B(L,2)*E)+ B(L,3)*EXP(B(L,4)*E)
261 1 DSQD=DSQD+A(L)*PL(L,X)
262 ELSEIF(E.LT..700)THEN
263 DSQD=.3
264 ELSEIF(E.LT..800)THEN
265 DSQD=.15
266 ELSEIF(E.LT..900)THEN
267 DSQD=.1
268 ELSE
269 DSQD=55./(E-.350)
270 ENDIF
271 RETURN
272 END
273 *LEG
274 jones 1.1 REAL*8 FUNCTION PL(L,X)
275 C LEGENDRE POLYNOMIALS
276 IMPLICIT REAL*8 (A-H,O-Z)
277 IF(L.EQ.0)THEN
278 PL=1.
279 ELSEIF(L.EQ.1)THEN
280 PL=X
281 ELSEIF(L.EQ.2)THEN
282 PL=.5*(3.*X**2-1.)
283 ELSEIF(L.EQ.3)THEN
284 PL=.5*(5.*X**3-3.*X)
285 ELSEIF(L.EQ.4)THEN
286 PL=1./8.*(35.*X**4-30.*X**2+3.)
287 ELSE
288 PL=0.
289 ENDIF
290 RETURN
291 END
292 *DELTA_OR
293 SUBROUTINE DELTA_OR(E1,TP,TH,D2DEL)
294 C PHOTOPRODUCTION OF NUCLEONS AND PIONS VIA DELTA
295 jones 1.1 IMPLICIT REAL*8 (A-H,O-Z)
296 COMMON/DEL/IP
297 CHARACTER*1 SCAL,FERMI
298 COMMON /FER/ FERMI,SCAL
299 DATA AM/939./,AMP/139./
300 IF(ABS(IP).EQ.1)THEN
301 AM1=AM
302 AM2=AMP
303 ELSE
304 AM1=AMP
305 AM2=AM
306 ENDIF
307 EP=TP+AM1
308 PN=SQRT(EP**2-AM1**2)
309 ! PRINT *
310 ! PRINT *,'DELTA',PN
311 IF(FERMI.NE.'Y') THEN
312 CALL KINE_OR(AM,AM1,AM2,PN,TH,W,TC)
313 CPF = 1.
314 ELSE
315 CALL KINDEL(AM,AM1,AM2,PN,TH,W,TC,CPF)
316 jones 1.1 END IF
317 ! PRINT *,'PN,W,CPF',PN,W,CPF
318 IF(W.LE.0.)GO TO 1
319 ! IF(W.GE.E1)GO TO 1
320 IF(W.GT.E1)GO TO 1
321 IF(FERMI.NE.'Y') THEN
322 CALL PART_OR(AM,AM1,AM2,PN,TH,AJT,AJW)
323 ELSE
324 CALL PARTDEL(AM,AM1,AM2,PN,TH,AJT,AJW)
325 END IF
326 CALL SIGMA_OR(W,TC,DSIGG)
327 IF(SCAL.NE.'Y') THEN
328 C CROSS SECTION IN UB/MEV-SR
329 CALL VTP_OR(AM,AM1,E1,W,TP,TH,PHI)
330 ! WRITE(*,100) AM,AM1,E1,W,TH,PHI
331 100 FORMAT(2X,2F8.1,5G12.3)
332 D2DEL = PHI*DSIGG*AJT*CPF
333 ! PRINT *,'PHI,AJT,DSIGG,D2DEL,CPF',PHI,AJT,DSIGG,D2DEL,CPF
334 ! PRINT *
335 ELSE
336 C CROSS SECTION IN UB/MEV-SR-Q
337 jones 1.1 BREM = 1./W
338 DSIG = BREM*DSIGG*AJT*AJW
339 ! D2DEL = DSIG*CPF
340 ! D2DEL = D2DEL*EP*1.0D6/PN**2 ! Convert to ub.c/(Q.(GeV/c)^2.sr)
341 D2DEL = DSIG*EP/PN
342 D2DEL = D2DEL*1.0D6/PN ! Convert to units ub.c/(Q.(GeV/c)^2.sr)
343 D2DEL = D2DEL*CPF
344 END IF
345 RETURN
346 1 D2DEL=0.
347 RETURN
348 END
349 *PART_OR
350 SUBROUTINE PART_OR(AMT,AM1,AM2,PN,TN,AJT,AJW)
351 C PARTIAL DERIVATIVES
352 IMPLICIT REAL*8 (A-H,O-Z)
353 PI=ACOS(-1.0D0)
354 ! DT=PI/180.0D0
355 ! DP=10.0D0
356 DT=PI/720.0D0
357 DP=2.0D0
358 jones 1.1 C ANGLE
359 TNP=TN+DT
360 TNM=TN-DT
361 TN0=TN
362 CALL KINE_OR(AMT,AM1,AM2,PN,TNP,WP,TCP)
363 CALL KINE_OR(AMT,AM1,AM2,PN,TNM,WM,TCM)
364 CALL KINE_OR(AMT,AM1,AM2,PN,TN0,W,TC0)
365 DEN=COS(TNP)-COS(TNM)
366 DEN=ABS(DEN)
367 ! PRINT *,'DEN,TCP,TCM,TC0',DEN,TCP,TCM,TC0
368 IF(DEN.GT.1.0D-3.AND.(W*WP*WM.GT.0.))THEN
369 AJT=(COS(TCP)-COS(TCM))/DEN
370 AJT=ABS(AJT)
371 ELSE
372 AJT=(COS(TC0)-COS(TCM))/(COS(TN0)-COS(TNM))
373 AJT=ABS(AJT)
374 ENDIF
375 C ENERGY
376 PNP=PN+DP
377 PNM=PN-DP
378 CALL KINE_OR(AMT,AM1,AM2,PNP,TN,WP,TC)
379 jones 1.1 CALL KINE_OR(AMT,AM1,AM2,PNM,TN,WM,TC)
380 AM12 = AM1**2
381 TP = SQRT(PNP**2 + AM12) - AM1
382 TM = SQRT(PNM**2 + AM12) - AM1
383 AJW=(WP-WM)/(PNP-PNM)
384 ! AJW=(WP-WM)/(TP-TM)
385 AJW=ABS(AJW)
386 RETURN
387 END
388 *KINE_OR
389 SUBROUTINE KINE_OR(AMT,AM1,AM2,PN,TH,W,TC)
390 C COMPUTES CM VARIABLES FROM LAB VARIABLES
391 IMPLICIT REAL*8 (A-H,O-Z)
392 COMMON /KF/ E1
393 AMT2 = AMT**2
394 AM12 = AM1**2
395 EP=SQRT(PN**2+AM12)
396 PNT=PN*SIN(TH)
397 PNL=PN*COS(TH)
398 ! ANUM=PN**2+AM2**2-(AMT-EP)**2
399 ANUM = AM2**2 - AM12 - AMT2 +2*EP*AMT
400 jones 1.1 DEN=2.*(PNL+AMT-EP)
401 W=ANUM/DEN
402 IF(W.LE.0.)W=0.
403 C INVARIANT MASS
404 WW2 = AMT**2+2.*W*AMT
405 WW = SQRT(WW2)
406 ! PRINT *,' P,W,E_thr ',PN,WW,W,ANUM/DEN
407 C CM VARIABLES
408 PCT=PNT
409 B=W/(AMT+W)
410 G=(W+AMT)/WW
411 PCL=G*(PNL-B*EP)
412 PCS=PCL**2+PCT**2
413 PC=SQRT(PCS)
414 CTHC=PCL/PC
415 TC=ACOS(CTHC)
416 RETURN
417 END
418 *SIGMA_OR
419 SUBROUTINE SIGMA_OR(E,THRCM,SIGCM)
420 C REAL PHOTON CROSS SECTION IN DELTA REGION
421 jones 1.1 C MICROBARNS PER STERADIAN
422 IMPLICIT REAL*8 (A-H,O-Z)
423 GAM=100.
424 PI=ACOS(-1.0D0)
425 IF(E.GT.420.)THEN
426 SIGCM=(1.+420./E)*90./4./PI
427 ELSE
428 SIGCM=360.*(5.-3.*COS(THRCM)**2)
429 SIGCM=SIGCM/16./PI/(1.+(E-320)**2/GAM**2)
430 ENDIF
431 RETURN
432 END
433 *EP_OR
434 SUBROUTINE EP_OR(E1,TP,THP,DSEP)
435 C ELECTRO PROTON PRODUCTION CROSS SECTIONS
436 IMPLICIT REAL*8 (A-H,O-Z)
437 COMMON/Pcm/PH(10),WPH(10)
438 DATA AML/.511/
439 PI=ACOS(-1.0D0)
440 CALL GAUSAB_OR(10,PH,WPH,0.D0,2.*PI,PI)
441 AK=SQRT(E1**2-AML**2)
442 jones 1.1 CALL SEP_OR(AK,TP,THP,DSEP)
443 DSEP=DSEP*1.E4
444 C CROSS SECTION IN UB/MEV-SR
445 END
446 *DOT_OR
447 REAL*8 FUNCTION DOT_OR(V,U)
448 IMPLICIT REAL*8 (A-H,O-Z)
449 DIMENSION V(3),U(3)
450 DOT_OR=0.
451 DO 1 I=1,3
452 1 DOT_OR=DOT_OR+V(I)*U(I)
453 RETURN
454 END
455 *CROSS
456 SUBROUTINE CROSS_OR(V,U,W)
457 IMPLICIT REAL*8 (A-H,O-Z)
458 DIMENSION V(3),U(3),W(3)
459 W(1)=V(2)*U(3)-V(3)*U(2)
460 W(2)=V(3)*U(1)-V(1)*U(3)
461 W(3)=V(1)*U(2)-V(2)*U(1)
462 RETURN
463 jones 1.1 END
464 *GAUSAB_OR
465 SUBROUTINE GAUSAB_OR(N,E,W,A,B,C)
466 IMPLICIT REAL*8 (A-H,O-Z)
467 DIMENSION E(24),W(24)
468 DATA EPS/1.D-16/
469 IF(A.GE.C.OR.C.GE.B)STOP
470 C STOPS PROGRAM IF A, C, B ARE OUT OF SEQUENCE
471 PI=ACOS(-1.0D0)
472 AL=(C*(A+B)-2*A*B)/(B-A)
473 BE=(A+B-2*C)/(B-A)
474 M=(N+1)/2
475 DN=N
476 DO 5 I=1,M
477 DI=I
478 X=PI*(4.D0*(DN-DI)+3.D0)/(4.D0*DN+2.D0)
479 XN=(1.D0-(DN-1.D0)/(8.D0*DN*DN*DN))*COS(X)
480 IF(I.GT.N/2) XN=0
481 DO 3 ITER=1,10
482 X=XN
483 Y1=1.D0
484 jones 1.1 Y=X
485 IF(N.LT.2) GO TO 2
486 DO 1 J=2,N
487 DJ=J
488 Y2=Y1
489 Y1=Y
490 1 Y=((2.D0*DJ-1.D0)*X*Y1-(DJ-1.D0)*Y2)/DJ
491 2 CONTINUE
492 YS=DN*(X*Y-Y1)/(X*X-1.D0)
493 H=-Y/YS
494 XN=X+H
495 IF(ABS(H).LT.EPS) GO TO 4
496 3 CONTINUE
497 4 E(I)=(C+AL*X)/(1.D0-BE*X)
498 E(N-I+1)=(C-AL*X)/(1.D0+BE*X)
499 GEW=2.D0/((1.D0-X*X)*YS*YS)
500 W(I)=GEW*(AL+BE*C)/(1.D0-BE*X)**2
501 W(N-I+1)=GEW*(AL+BE*C)/(1.D0+BE*X)**2
502 5 CONTINUE
503 RETURN
504 END
505 jones 1.1 *VECT_OR
506 SUBROUTINE VECT_OR(THP,THE,PHI,P,AK1,AK2)
507 C CARTESIAN COMPONENTS OF ELECTRON AND PROTON VECTORS
508 IMPLICIT REAL*8 (A-H,O-Z)
509 COMMON/V/AK1V(3),AK2V(3),QV(3),PV(3),PP(3)
510 PV(1)=P*SIN(THP)
511 PV(2)=0.
512 PV(3)=P*COS(THP)
513 AK1V(1)=0.
514 AK1V(2)=0.
515 AK1V(3)=AK1
516 AK2V(1)=AK2*SIN(THE)*COS(PHI)
517 AK2V(2)=AK2*SIN(THE)*SIN(PHI)
518 AK2V(3)=AK2*COS(THE)
519 QV(1)=AK1V(1)-AK2V(1)
520 QV(2)=AK1V(2)-AK2V(2)
521 QV(3)=AK1V(3)-AK2V(3)
522 PP(1)=PV(1)-QV(1)
523 PP(2)=PV(2)-QV(2)
524 PP(3)=PV(3)-QV(3)
525 RETURN
526 jones 1.1 END
527 *AMAG_OR
528 REAL*8 FUNCTION AMAG_OR(V)
529 IMPLICIT REAL*8 (A-H,O-Z)
530 DIMENSION V(3)
531 AMAG_OR=0.
532 DO 1 I=1,3
533 1 AMAG_OR=AMAG_OR+V(I)**2
534 AMAG_OR=SQRT(AMAG_OR)
535 RETURN
536 END
537 *LEPT_OR
538 SUBROUTINE LEPT_OR(E1,E2,AK1,AK2,AML,QS,QUS,THE,V)
539 C LEPTON FACTORS FOR COINCIDENCE CROSS SECTION
540 IMPLICIT REAL*8 (A-H,O-Z)
541 DIMENSION V(5)
542 V(1)=(QUS/QS)**2*(E1*E2+AK1*AK2*COS(THE)+AML**2)
543 X=AK1*AK2*SIN(THE)
544 V(2)=X**2/QS+QUS/2.
545 V(3)=QUS/QS*X/SQRT(QS)*(E1+E2)
546 V(4)=X**2/QS
547 jones 1.1 V(5)=0.
548 RETURN
549 END
550 *D4S_OR
551 SUBROUTINE D4S_OR(AK1,AK2,THE,P,PP,THQP,CPHIP,DSIG)
552 C FULLY DIFFERENTIAL CROSS SECTION
553 IMPLICIT REAL*8 (A-H,O-Z)
554 DIMENSION V(5),W(5)
555 DATA AM/939./,AML/.511/,A/855./
556 QS=AK1**2+AK2**2-2.*AK1*AK2*COS(THE)
557 E1=SQRT(AK1**2+AML**2)
558 E2=SQRT(AK2**2+AML**2)
559 QUS=2.*(E1*E2-AK1*AK2*COS(THE)-AML**2)
560 SM=2.*(1.44)**2/QUS**2*AK2/AK1
561 EP=SQRT(AM**2+P**2)
562 PS=EP*P
563 FNS=1./(1.+QUS/A**2)**4
564 CALL LEPT_OR(E1,E2,AK1,AK2,AML,QS,QUS,THE,V)
565 CALL FORM_OR(QS,P,THQP,CPHIP,W)
566 SUM=0.
567 DO 1 I=1,5
568 jones 1.1 1 SUM=SUM+V(I)*W(I)
569 DSIG=SM*PS*FNS*SUM*SGSL_OR(PP)
570 RETURN
571 END
572 *STHE_OR
573 SUBROUTINE STHE_OR(D2S)
574 C INTEGRAL OVER ELECTRON POLAR ANGLE
575 IMPLICIT REAL*8 (A-H,O-Z)
576 COMMON/S/ AK1,AK2,THE,P,THP
577 COMMON/E/TH1(12),WT1(12),TH2(12),WT2(12)
578 COMMON/E1/TH3(24),WT3(24)
579 D2S1=0.
580 DO 1 I=1,12
581 THE=TH1(I)
582 CALL SPHI_OR(D3S)
583 1 D2S1=D2S1+D3S*WT1(I)*SIN(THE)
584 D2S2=0.
585 DO 2 I=1,12
586 THE=TH2(I)
587 CALL SPHI_OR(D3S)
588 2 D2S2=D2S2+D3S*WT2(I)*SIN(THE)
589 jones 1.1 D2S3=0.
590 DO 3 I=1,24
591 THE=TH3(I)
592 CALL SPHI_OR(D3S)
593 3 D2S3=D2S3+D3S*WT3(I)*SIN(THE)
594 D2S=D2S1+D2S2+D2S3
595 RETURN
596 END
597 *SPHI_OR
598 SUBROUTINE SPHI_OR(D3S)
599 C INTEGRATE OVER ELECTRON AZIMUTHAL ANGLE
600 IMPLICIT REAL*8 (A-H,O-Z)
601 COMMON/S/ AK1,AK2,THE,P,THP
602 COMMON/V/AK1V(3),AK2V(3),QV(3),PV(3),PP(3)
603 COMMON/Pcm/PH(10),WPH(10)
604 DIMENSION QXP(3),AK1X2(3)
605 D3S=0.
606 DO 1 I=1,10
607 PHI=PH(I)
608 CALL VECT_OR(THP,THE,PHI,P,AK1,AK2)
609 CALL CROSS_OR(QV,PV,QXP)
610 jones 1.1 CALL CROSS_OR(AK1V,AK2V,AK1X2)
611 C PROTON THETA
612 CTHEP=DOT_OR(PV,QV)/AMAG_OR(PV)/AMAG_OR(QV)
613 THQP=ACOS(CTHEP)
614 C PROTON PHI
615 CPHIP=DOT_OR(QXP,AK1X2)
616 IF (CPHIP.EQ.0.) THEN
617 CPHIP=1.
618 ELSE
619 CPHIP=CPHIP/AMAG_OR(QXP)/AMAG_OR(AK1X2)
620 ENDIF
621 PPM=AMAG_OR(PP)
622 CALL D4S_OR(AK1,AK2,THE,P,PPM,THQP,CPHIP,DSIG)
623 1 D3S=D3S+DSIG*WPH(I)
624 RETURN
625 END
626 *FORM_OR
627 SUBROUTINE FORM_OR(QS,P,THQP,CPHIP,W)
628 C NUCLEAR FORM FACTORS
629 IMPLICIT REAL*8 (A-H,O-Z)
630 COMMON/M/ZZ,NN
631 jones 1.1 COMMON/DEL/IP
632 DIMENSION W(5)
633 DATA AM/939./,UP/2.79/,UN/-1.91/
634 IF(IP.EQ.1)THEN
635 Z=ZZ
636 N=0.
637 ELSEIF(IP.EQ.-1)THEN
638 Z=0.
639 N=NN
640 ELSE
641 Z=0.
642 N=0.
643 ENDIF
644 Y=P/AM*SIN(THQP)
645 W(1)=Z
646 W(2)=Z*Y**2
647 W(2)=W(2)+(Z*UP**2+N*UN**2)*QS/2./AM**2
648 W(3)=-2.*Z*Y*CPHIP
649 W(4)=Z*Y**2*(2.*CPHIP**2-1.)
650 W(5)=0.
651 RETURN
652 jones 1.1 END
653 *SEP_OR
654 SUBROUTINE SEP_OR(AK,TP,THPP,D2S)
655 IMPLICIT REAL*8 (A-H,O-Z)
656 COMMON/S/AK1,AK2,THE,P,THP
657 COMMON/E/TH1(12),WT1(12),TH2(12),WT2(12)
658 COMMON/E1/TH3(24),WT3(24)
659 DATA AM/939./,AML/.511/,BE/16./
660
661 PI=ACOS(-1.0D0)
662 THP=THPP
663 AK1=AK
664 AK2=AK1-TP-BE
665 C GAUSSIAN POINTS FOR THE
666 THEMAX=AML*(AK1-AK2)/AK1/AK2
667 CALL GAUSAB_OR(12,TH1,WT1,0.D0,2.*THEMAX,THEMAX)
668 CALL GAUSAB_OR(12,TH2,WT2,2.*THEMAX,100.*THEMAX,10.*THEMAX)
669 A3=100.*THEMAX
670 C3=A3+(PI-A3)/10.
671 CALL GAUSAB_OR(24,TH3,WT3,A3,PI,C3)
672 P=SQRT(TP**2+2.*TP*AM)
673 jones 1.1 CALL STHE_OR(D2S)
674 IF(AK2.LE.0.)D2S=0.
675 RETURN
676 END
677 *SGSL_OR
678 REAL*8 FUNCTION SGSL_OR(P)
679 IMPLICIT REAL*8 (A-H,O-Z)
680 C P INTEGRAL OVER SGSL_OR NORMALIZED TO 1/4PI
681 COMMON/SG/IA
682 COMMON/QD/QDF
683 IF(IA.EQ.2)THEN
684 C BEGIN 2-H
685 PP=P/197.3
686 SGS=3.697-7.428*PP-2.257*PP**2
687 SGS=SGS+3.618*PP**3-1.377*PP**4+.221*PP**5-.013*PP**6
688 IF(SGS.LT.-293.)GO TO 1
689 SGS=EXP(SGS)
690 SGS=SGS/.18825/4./3.1416/(197.3)**3
691 ! SGSL_OR=SGS/1.
692 ELSEIF(IA.EQ.3)THEN
693 C BEGIN 3-HE
694 jones 1.1 IF(-(P/33)**2.LT.-293.)GO TO 1
695 SGS=2.4101E-6*EXP(-P/33)
696 SGS=SGS-1.4461E-6*EXP(-(P/33)**2)
697 SGS=SGS+1.6871E-10*EXP(-(P/493)**2)
698 SGSL_OR=SGS/2.0D0
699 ELSEIF(IA.EQ.4)THEN
700 C BEGIN 4-HE
701 IF(-(P/113.24)**2.LT.-293.)GO TO 1
702 SGS=1.39066E-6*EXP(-(P/113.24)**2)
703 SGS=SGS+3.96476E-9*EXP(-(P/390.75)**2)
704 SGSL_OR=SGS/2.0D0
705 SGSL_OR=SGSL_OR/2.0D0/3.1416
706 ELSEIF(IA.GT.4.AND.IA.LT.12)THEN
707 IF(-(P/127)**2.LT.-293.)GO TO 1
708 SGS=1.7052E-7*(1.+(P/127)**2)*EXP(-(P/127)**2)
709 SGS=SGS+1.7052E-9*EXP(-(P/493)**2)
710 SGSL_OR=SGS/(FLOAT(IA)/2.0D0)
711 ELSEIF(IA.EQ.12)THEN
712 C BEGIN 12-C
713 IF(-(P/127)**2.LT.-293.)GO TO 1
714 SGS=1.7052E-7*(1.+(P/127)**2)*EXP(-(P/127)**2)
715 jones 1.1 SGS=SGS+1.7052E-9*EXP(-(P/493)**2)
716 SGSL_OR=SGS/6.
717 ELSE
718 C BEGIN 16-O
719 IF(-(P/120)**2.LT.-293.)GO TO 1
720 SGS=3.0124E-7*(1.+(P/120)**2)*EXP(-(P/120)**2)
721 SGS=SGS+1.1296E-9*EXP(-(P/493)**2)
722 SGSL_OR=SGS/(FLOAT(IA)/2.0D0)
723 ENDIF
724 RETURN
725 1 SGSL_OR=0.
726 RETURN
727 END
728 *S2PI
729 SUBROUTINE S2PI(IP,E1,TP,TH,D2SC)
730 C INTEGRAL OVER SCALING CROSS SECTION
731 IMPLICIT REAL*8 (A-H,O-Z)
732 CHARACTER*1 SCAL,FERMI
733 COMMON /FER/ FERMI,SCAL
734 DATA AM/939./,AMP/139./
735 IF(ABS(IP).EQ.1)THEN
736 jones 1.1 C ONE PION THR
737 AP=AM
738 AM2=AMP
739 ELSEIF(IP.EQ.2.OR.IP.EQ.0)THEN
740 C ONE PION THR
741 AP=AMP
742 AM2=AM
743 ELSEIF(IP.EQ.-2)THEN
744 C TWO PION THR
745 AP=AMP
746 AM2=AM+AMP
747 ELSE
748 STOP
749 ENDIF
750 P=SQRT(TP**2+2.*AP*TP)
751 E=TP+AP
752 ! PRINT *,'S2PI',P
753 IF(FERMI.NE.'Y') THEN
754 CALL KINE_OR(AM,AP,AM2,P,TH,THR,TC)
755 CPF = 1.
756 ELSE
757 jones 1.1 CALL KINDEL(AM,AP,AM2,P,TH,THR,TC,CPF)
758 END IF
759 ! PRINT *,'P,THR',P,THR
760 IF(THR.LE.0.)GOTO 2
761 ! IF(E1.LE.THR)GOTO 2
762 IF(E1.LT.THR)GOTO 2
763 DW=(E1-THR)/20.0D0
764 SUM=0.0D0
765 SUMBR=0.0D0
766 ! DO 1 I=1,20
767 DO I=1,20
768 W=THR+(FLOAT(I)-.5D0)*DW
769 IF(W.LT.(E1-0.511)) THEN
770 CALL VTP_OR(AM,AMP,E1,W,TP,TH,GN)
771 ! WRITE(*,100) I,AM,AMP,E1,W,DW,THR,GN
772 100 FORMAT(I4,2F8.1,5G12.3)
773 CALL WISER(W/1.E3,P/1.E3,TH,F)
774 ! PRINT *,'P,W,GN,F',P,W,GN,F
775 SUM=SUM+GN*F*DW
776 SUMBR = SUMBR + F*DW/W
777 ! 1 CONTINUE
778 jones 1.1 END IF
779 END DO
780 IF(SCAL.EQ.'Y') THEN
781 D2SC=SUMBR
782 ELSE
783 D2SC=SUM*P**2/E*1.E-6
784 END IF
785 D2SC = D2SC*CPF
786 ! PRINT *,'D2PI',D2SC,CPF
787 RETURN
788 2 D2SC=0.
789 RETURN
790 END
791 *WISER
792 SUBROUTINE WISER(W,P,TH,F)
793 C INVARIANT INCLUSIVE CROSS SECTION
794 C UNITS IN GEV
795 IMPLICIT REAL*8 (A-H,O-Z)
796 COMMON/DEL/IP
797 DIMENSION A(7),B(7),C(7)
798 DATA A/5.66E2,8.29E2,1.79,2.10,-5.49,-1.73,0./
799 jones 1.1 DATA B/4.86E2,1.15E2,1.77,2.18,-5.23,-1.82,0./
800 DATA C/1.33E5,5.69E4,1.41,0.72,-6.77,1.90,-1.17E-2/
801 DATA AM/.939/,AMP/.139/
802 IF(ABS(IP).EQ.1)THEN
803 AP=AM
804 ELSEIF(ABS(IP).EQ.2.OR.IP.EQ.0)THEN
805 AP=AMP
806 ELSE
807 STOP
808 ENDIF
809 E=SQRT(P**2+AP**2)
810 C MANDELSTAM VARIABLES
811 S=2.*W*AM+AM**2
812 C T=-2.*W*E+2.*W*P*COS(TH) +AP**2
813 U=AM**2-2.*AM*E+AP**2
814 C FITTING VARIABLES
815 PT=P*SIN(TH)
816 AML=SQRT(PT**2+AP**2)
817 CALL FXR(W,P,E,TH,XR)
818 C FITTED
819 IF(IP.EQ.2.OR.IP.EQ.0)THEN
820 jones 1.1 X1=A(1)+A(2)/SQRT(S)
821 X2=(1.-XR+A(3)**2/S)**A(4)
822 X3=EXP(A(5)*AML)
823 X4=EXP(A(6)*PT**2/E)
824 F=X1*X2*X3*X4
825 ELSEIF(IP.EQ.-2)THEN
826 X1=B(1)+B(2)/SQRT(S)
827 X2=(1.-XR+B(3)**2/S)**B(4)
828 X3=EXP(B(5)*AML)
829 X4=EXP(B(6)*PT**2/E)
830 F=X1*X2*X3*X4
831 ELSEIF(ABS(IP).EQ.1)THEN
832 X1=C(1)+C(2)/SQRT(S)
833 X2=(1.-XR+C(3)**2/S)**C(4)
834 X3=EXP(C(5)*AML)
835 X4=1./(1.+ABS(U))**(C(6)+C(7)*S)
836 F=X1*X2*X3*X4
837 ELSE
838 STOP
839 ENDIF
840 RETURN
841 jones 1.1 END
842 *FXR
843 SUBROUTINE FXR(W,P,E,TH,XR)
844 C COMPUTES RATIO OF CM PARTICLE MOMENTUM TO PHOTON MOMENTUM
845 C GEV UNITS
846 IMPLICIT REAL*8 (A-H,O-Z)
847 DATA AM/.939/
848 PT=P*SIN(TH)
849 PL=P*COS(TH)
850 C LORENTZ TRANSFORMATION
851 B=W/(W+AM)
852 D=SQRT(2.*W*AM+AM**2)
853 G=(W+AM)/D
854 BG=B*G
855 C CM VARIABLES
856 WC=G*W-BG*W
857 PLC=G*PL-BG*E
858 PC=SQRT(PT**2+PLC**2)
859 XR=PC/WC
860 RETURN
861 END
862 jones 1.1 *KINDEL
863 SUBROUTINE KINDEL(AMT,AM1,AM2,PN,TH,W,TC,CPF)
864 C COMPUTES CM VARIABLES FROM LAB VARIABLES
865 IMPLICIT REAL*8 (A-H,O-Z)
866 COMMON /KF/ E1
867 DATA COSAV/0.63662/
868 AMT2 = AMT**2
869 AM12 = AM1**2
870 EP=SQRT(PN**2+AM12)
871 PNT=PN*SIN(TH)
872 PNL=PN*COS(TH)
873 ! ANUM=PN**2+AM2**2-(AMT-EP)**2
874 ANUM = AM2**2 - AM12 - AMT2 +2*EP*AMT
875 DEN=2.*(PNL+AMT-EP)
876 W=ANUM/DEN
877 IF(W.LE.0.)W=0.
878 C INVARIANT MASS
879 WW2 = AMT**2+2.*W*AMT
880 WW = SQRT(WW2)
881 ! PRINT *,' P,W,E_thr ',PN,WW,W,ANUM/DEN
882
883 jones 1.1 C Section on Fermi motion
884
885 IF(CPF.GE.0.) THEN ! CPF = -1. for Jacobian computation
886 CPF = 1.0D0
887 IF(W.GE.E1) THEN
888 DM2 = WW2 - AMT2
889 PFMIN = E1*((DM2/2./E1)**2 - AMT2)/DM2
890 CPF0 = SGSL_OR(0.d0)
891 CPF = SGSL_OR(PFMIN)
892 IF(CPF.GT.0.0D0) THEN
893 DPF = 10.0D0
894 PF = PFMIN
895 DSPF = CPF*PF
896 SUMPF = DSPF
897 SUMCF = CPF
898 DO WHILE((DSPF/SUMPF).GT.1.0D-4)
899 PF = PF + DPF
900 CPF = SGSL_OR(PF)
901 DSPF = CPF*PF
902 SUMPF = DSPF + SUMPF
903 SUMCF = CPF + SUMCF
904 jones 1.1 END DO
905 ! PRINT *,'PF,DSPF,SUMPF',PF,DSPF,SUMPF
906 PFMEAN = SUMPF/SUMCF
907 PFMEAN = PFMEAN*COSAV
908 CPF = SGSL_OR(PFMEAN)/CPF0
909 CPF = 0.5*CPF
910 EF = SQRT(PFMEAN**2 + AMT**2)
911 WEFF = DM2/(2.0D0*(EF + PFMEAN))
912 ! WEFF = DM2/(2.0D0*EF)
913 ! PRINT *,'W,WW,PFMIN,CPF0',W,WW,PFMIN,CPF0
914 ! PRINT *,'PFMEAN,CPF,WEFF',PFMEAN,CPF,WEFF
915 W = WEFF
916 ELSE
917 W = 0.
918 END IF
919 END IF
920 END IF
921
922 C CM VARIABLES
923 PCT=PNT
924 B=W/(AMT+W)
925 jones 1.1 G=(W+AMT)/WW
926 PCL=G*(PNL-B*EP)
927 PCS=PCL**2+PCT**2
928 PC=SQRT(PCS)
929 CTHC=PCL/PC
930 TC=ACOS(CTHC)
931 RETURN
932 END
933 *PARTDEL
934 SUBROUTINE PARTDEL(AMT,AM1,AM2,PN,TN,AJT,AJW)
935 C PARTIAL DERIVATIVES
936 IMPLICIT REAL*8 (A-H,O-Z)
937 PI=ACOS(-1.0D0)
938 ! DT=PI/180.0D0
939 ! DP=10.0D0
940 DT=PI/720.0D0
941 DP=2.0D0
942 C ANGLE
943 TNP=TN+DT
944 TNM=TN-DT
945 TN0=TN
946 jones 1.1 CPF = -1.
947 CALL KINDEL(AMT,AM1,AM2,PN,TNP,WP,TCP,CPF)
948 CPF = -1.
949 CALL KINDEL(AMT,AM1,AM2,PN,TNM,WM,TCM,CPF)
950 CPF = -1.
951 CALL KINDEL(AMT,AM1,AM2,PN,TN0,W,TC0,CPF)
952 DEN=COS(TNP)-COS(TNM)
953 DEN=ABS(DEN)
954 ! PRINT *,'DEN,TCP,TCM,TC0',DEN,TCP,TCM,TC0
955 IF(DEN.GT.1.0D-3.AND.(W*WP*WM.GT.0.))THEN
956 AJT=(COS(TCP)-COS(TCM))/DEN
957 AJT=ABS(AJT)
958 ELSE
959 AJT=(COS(TC0)-COS(TCM))/(COS(TN0)-COS(TNM))
960 AJT=ABS(AJT)
961 ENDIF
962 C ENERGY
963 PNP=PN+DP
964 PNM=PN-DP
965 CALL KINE_OR(AMT,AM1,AM2,PNP,TN,WP,TC)
966 CALL KINE_OR(AMT,AM1,AM2,PNM,TN,WM,TC)
967 jones 1.1 AM12 = AM1**2
968 TP = SQRT(PNP**2 + AM12) - AM1
969 TM = SQRT(PNM**2 + AM12) - AM1
970 AJW=(WP-WM)/(PNP-PNM)
971 ! AJW=(WP-WM)/(TP-TM)
972 AJW=ABS(AJW)
973 RETURN
974 END
975 *OARA
976 SUBROUTINE OARA(W,P,TH,DSDWDP)
977 IMPLICIT REAL*8 (A-H,M-Z)
978 INTEGER N
979 COMMON/SG/IA
980 COMMON/QD/QDF
981 COMMON/DEL/IP
982 CHARACTER*1 SCAL,FERMI
983 COMMON /FER/ FERMI,SCAL
984 COMMON/M/Z,N
985 COMMON/KF/E1
986
987 DATA MPI/0.140/
988 jones 1.1 DATA B/7.0D0/,C/9.0D0/,STOT/125.0D0/
989 DATA CP/8.0D3/,T0/124.0D-3/
990
991 IF(ABS(IP).EQ.1) THEN
992 DSDWDP = STOT*C*EXP(-B*P**2)
993 ELSE IF (ABS(IP).EQ.2.OR.IP.EQ.0) THEN
994 T = SQRT(P**2+MPI**2) - MPI
995 DSDWDP = CP*EXP(-T/T0)
996 END IF
997 RETURN
998 END
999
|
Return to
epc_or.f
CVS log
Up to [HallC] / sane_geant_mc / SRC