EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lxsect.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lxsect.F
1 
2 C **********************************************************************
3 
4  SUBROUTINE lxsect
5 
6  IMPLICIT NONE
7 
8 C...Integrate differential cross-section using GADAP, RIWIAD or DIVONNE
9 
10  COMMON /linteg/ ntot,npass
11  INTEGER ntot,npass
12  SAVE /linteg/
13 
14 *
15 * to avoid variable conflictions, a second keep element is necessary
16 * with the same common block name (see LPTOU2)
17 *
18  COMMON /leptou/ cut(14),lst(40),parl(30),
19  & x,y,w2,q2,u
20  REAL cut,parl,x,y,w2,q2,u
21  INTEGER lst
22  SAVE /leptou/
23 
24  COMMON /lintrl/ psave(3,4,5),ksave(4),xmin,xmax,ymin,ymax,
25  &q2min,q2max,w2min,w2max,ilep,inu,ig,iz
26  REAL psave,xmin,xmax,ymin,ymax,q2min,q2max,w2min,w2max
27  INTEGER ksave,ilep,inu,ig,iz
28  SAVE /lintrl/
29 
30 * this Common block appears only in Subroutine LXSECT
31  COMMON /params/ acc,ndim,nsub,iter
32  DOUBLE PRECISION acc
33  INTEGER ndim,nsub,iter
34  SAVE /params/
35 
36 * this Common block appears only in Subroutine LXSECT
37  COMMON /answer/ value,erriw
38  DOUBLE PRECISION value,erriw
39  SAVE /answer/
40 
41 * this Common block appears only in Subroutine LXSECT
42  COMMON /bndlmt/ flow,fhigh
43  DOUBLE PRECISION flow,fhigh
44  SAVE /bndlmt/
45 
46 * this Common block appears only in Subroutine LXSECT
47  COMMON /sample/ npoint
48  INTEGER npoint
49  SAVE /sample/
50 
51 
52  INTEGER ncall,i,maxnum,maxpts,npt,it,jdeg
53  REAL xminus,xplus,ti1,sigma,eps,sprdmx,accur,ti2,errest
54  dimension xminus(2),xplus(2)
55  EXTERNAL dcross,dlower,dupper,riwfun
56  DATA ncall/0/
57 
58  ncall=ncall+1
59  CALL ltimex(ti1)
60  ntot=0
61  npass=0
62  sigma=0.
63  errest=0.
64  ndim=2
65 C...Parameters for RIWIAD integration.
66  acc=parl(15)
67  nsub=100
68  iter=100
69 C...Parameters for DIVON integration.
70  DO 20 i=1,2
71  xminus(i)=0.
72  20 xplus(i)=1.
73  eps=parl(15)
74  maxnum=50000
75  flow=-1.d0
76  fhigh=1.d+20
77  npoint=100
78 C...Additional parameters for detailed DIVON integration.
79  sprdmx=2.
80  maxpts=50000
81  jdeg=0
82  npt=1000
83 
84  sigma=0.
85  errest=0.
86  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,1000)
87  IF(lst(10).EQ.1) THEN
88 C...Integration using GADAP.
89  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)) WRITE(6,1100)
90  accur=parl(15)
91  it=0
92  100 it=it+1
93  errest=accur
94  CALL gadap2(xmin,xmax,dlower,dupper,dcross,errest,sigma)
95  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
96  & WRITE(6,1110) it,ntot,npass,sigma
97  IF(sigma.GT.1.) THEN
98  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
99  & WRITE(6,1120) accur
100  ELSE
101  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1))
102  & WRITE(6,1130) accur,accur/max(1.e-22,sigma),parl(15)
103  accur=max(1.e-22,sigma*parl(15))
104  IF(it.LT.2) goto 100
105  ENDIF
106  ELSEIF(lst(10).EQ.2) THEN
107 C...Integration using RIWIAD. When RIWIAD cannot be loaded:
108 C...activate next two lines and deactivate RIWIAD call.
109  WRITE(6,*) ' RIWIAD not available, execution stopped.'
110  stop
111 C.HI IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1))
112 C.HI & WRITE(6,1200) SNGL(ACC),NSUB,ITER
113 C CALL RIWIAD(RIWFUN)
114 C.HI SIGMA=SNGL(VALUE)
115 C.HI ERREST=SNGL(ERRIW)
116  ELSEIF(lst(10).EQ.3) THEN
117 C...Integration using simple DIVONNE. When DIVONNE cannot be loaded:
118 C...activate next two lines and deactivate DIVONNE call.
119  WRITE(6,*) ' DIVONNE not available, execution stopped.'
120  stop
121 C.HI IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1))
122 C.HI & WRITE(6,1300) EPS,MAXNUM,SNGL(FLOW),SNGL(FHIGH),NPOINT
123 C CALL DIVON(NDIM,XMINUS,XPLUS,EPS,MAXNUM,SIGMA,ERREST)
124  ELSEIF(lst(10).EQ.4) THEN
125 C...Integration using detailed DIVONNE. When DIVONNE cannot be loaded:
126 C...activate next two lines and deactivate PARTN and INTGRL calls.
127  WRITE(6,*) ' DIVONNE not available, execution stopped.'
128  stop
129 C.HI IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1))
130 C.HI & WRITE(6,1400) EPS,MAXNUM,
131 C.HI & SNGL(FLOW),SNGL(FHIGH),NPOINT,SPRDMX,MAXPTS,JDEG,NPT
132 C CALL PARTN(NDIM,XMINUS,XPLUS,SPRDMX,MAXPTS)
133 C CALL INTGRL(NDIM,JDEG,NPT,SIGMA,ERREST)
134  ELSE
135  IF(lst(3).GE.1) WRITE(6,*) ' Warning: LST(10) = ',lst(10),
136  & ' not allowed.'
137  ENDIF
138  CALL ltimex(ti2)
139  IF(lst(3).GE.4.OR.(lst(3).EQ.3.AND.ncall.EQ.1)
140  &.OR.(lst(3).GE.1.AND.npass.EQ.0)) THEN
141  WRITE(6,1500) sigma,errest,ntot,npass,ti2-ti1
142  IF(lst(3).GE.1.AND.npass.EQ.0) WRITE(6,1600)
143  ENDIF
144  parl(23)=sigma
145 
146  RETURN
147  1000 FORMAT(/,' Integration of cross section:',/,1x,28('-'))
148  1100 FORMAT(5x,'using GADAP = adaptive Gaussian integration')
149  1110 FORMAT(5x,'Iteration #',i3,/,
150  &10x,'# function evaluations; total & non-zero =',2i8,/,
151  &10x,'sigma =',g10.2,' pb')
152  1120 FORMAT(10x,'required relative error = ',g10.2)
153  1130 FORMAT(10x,'effective absolute error = ',g10.2,/,
154  & 10x,'effective relative error = ',g10.2,/,
155  & 10x,'required relative error = ',g10.2)
156  1200 FORMAT(5x,'using RIWIAD with parameters: rel. acc. = ',f10.4,
157  &/,5x,'# of subvolumes = ',i5,5x,'max # iterations = ',i5)
158  1300 FORMAT(5x,'using automatic DIVONNE with parameters: ',
159  &'rel. acc. = ',f10.4,/,5x,'max # function calls = ',i5,
160  &/,5x,'lower and upper bound on integrand =',2e12.4,
161  &/,5x,'# sample points/region =',i5)
162  1400 FORMAT(5x,'using detailed DIVONNE with parameters: ',
163  &'rel. acc. = ',f10.4,/,5x,'max # function calls = ',i5,
164  &/,5x,'lower and upper bound on integrand =',2e12.4,
165  &/,5x,'# sample points/region =',i5,
166  &/,5x,'SPRDMX, MAXPTS, JDEG, NPT =',f5.2,3i10)
167  1500 FORMAT(/,' ===> Cross-section =',1p,g12.3,
168  &' pb, error estimate = ',g12.3,/,
169  &6x,'# of integrand evaluations; total & non-zero =',2i8,/,
170  &6x,'cpu time for integration =',g12.3,' seconds',/)
171  1600 FORMAT(' Warning: integrand always zero, probably no allowed',
172  &' phase space due to cuts',/,
173  &10x,'check, in particular, CUT(11) to CUT(14)')
174  END