]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TUHKMgen/PYQUEN/progs_fortran.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixed warnings, coding conventions, updates (Ionut, Ludmila)
[u/mrichter/AliRoot.git] / TUHKMgen / PYQUEN / progs_fortran.f
CommitLineData
b1c2e580 1c*******************************************
2
3 SUBROUTINE myini
4
5 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
6 double precision nbcol,npart,npart0
7 real p,v,plu
8 external ludata,pydata
9 double precision numpar,npar0,nbco0
7b7936e9 10c common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
b1c2e580 11 common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
12 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
13 common /hyflow/ ytfl,ylfl,Tf,fpart
14 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
15 common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
16 common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
17 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
18 common /pyqpar/ T0,tau0,nf,ienglu,ianglu
19 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
20 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
21
22 COMMON/PYDATR/MRPY(6),RRPY(100)
23 common/SERVICE/iseed_fromC,iPythDecay,parPYTH(100)
24
25
26 common /hypart/ppart(10,50000),bmin,bmax,njp
7b7936e9 27c save /lujets/,/hyjets/,/hyflow/,/hyjpar/,/hyfpar/,/pyqpar/,
28c > /pysubs/,/pypars/,/pydat1/
b1c2e580 29
30* ----------- INITIALIZATION OF RUNDOM GENERATOR
31 MRPY(1)=iseed_fromC
32 write(*,*)'--initialization of high-pt part-- : '
33* write(*,*)'evnt generator phase ',iseed_fromC
34
35* initialize HYINIT with the input parameters
36 write(*,*)' AW= ', AW,' energy= ', ene, ' ptmin= ',ptmin
37 write(*,*)' bminh= ',bminh,' bmaxh= ', bmaxh
38 write(*,*)' ifb= ', ifb, ' bfix= ', bfix
39 write(*,*)' ishad= ',ishad, ' nhsel= ',nhsel
40 write(*,*)' ienglu= ',ienglu, ' ianglu= ',ianglu
41 write(*,*)' T0= ',T0,' tau0= ',tau0,' nf= ', nf
42
43* set input PYTHIA parameters from SERVICE common:
44 CKIN(3)=ptmin ! minimum pt in initial hard scattering, GeV
45 MSEL=0 ! QCD-dijet production
46 MSTU(21)=1 ! avoid stopping run
47 PARU(14)=1. ! tolerance parameter to adjust fragmentation
48 MSTP(81)=1 ! pp multiple scattering off
49 MSTJ(21)=0 ! hadron decays off
50 MSTP(2)=1 !which order running alphaS
51 MSTP(33)=0 !unclusion of k factor in cross section
52 MSTP(51)=7 !structure function chosen - CTEQ5M pdf
53 MSTP(82)=4 !Defines the multi-parton model
54 PARP(67)=1 !amount of initial-state radiation
55 PARP(82)=2. !pt cutoff for multiparton interactions
56 MSTJ(11)=3 !Choice of the fragmentation function
57
58* flags to select dijet production channels (active if msel=0)
59 MSUB(11)=1 ! q+q->q+q
60 MSUB(12)=1 ! q+qbar->q+qbar
61 MSUB(13)=1 ! q+qbar->g+g
62 MSUB(28)=1 ! q+g->q+g
63 MSUB(53)=1 ! g+g->q+qbar
64 MSUB(68)=1 ! g+g->g+g
65
66* flags to select prompt photon production channels (active if msel=0)
67 MSUB(14)=1 ! q+qbar->g+gamma
68 MSUB(18)=1 ! q+qbar->gamma+gamma
69 MSUB(29)=1 ! q+g->q+gamma
70 MSUB(114)=1 ! g+g->gamma+gamma
71 MSUB(115)=1 ! g+g->g+gamma
72
73
74 write(*,*)'--initialization of PYTHIA parameters-- : '
75
76 write(*,*)'CKIN(3) MSEL MSTU(21) PARU(14)'
77 write(*,*)CKIN(3), MSEL, MSTU(21), PARU(14)
78
79 write(*,*)'MSTP(81) MSTJ(21) MSTP(2) MSTP(33)'
80 write(*,*)MSTP(81), MSTJ(21), MSTP(2), MSTP(33)
81
82 write(*,*)'MSTP(51) MSTP(82) PARP(67) PARP(82) MSTJ(11)'
83 write(*,*)MSTP(51),MSTP(82),PARP(67),PARP(82),MSTJ(11)
84
85 write(*,*)'MSUB(11) MSUB(12) MSUB(13) MSUB(28) MSUB(53) MSUB(68)'
86 write(*,*)MSUB(11),MSUB(12),MSUB(13),MSUB(28),MSUB(53),MSUB(68)
87
88 write(*,*)'MSUB(14) MSUB(18) MSUB(29) MSUB(114) MSUB(115)'
89 write(*,*)MSUB(14),MSUB(18),MSUB(29),MSUB(114),MSUB(115)
90
91 call hyinit ! (ene,AW,ifb,bmin,bmax,bfix1,nh)
92
93 end
94****************end myini***************************************************************
95
96
97
98********************************* HYINIT ***************************
99 SUBROUTINE HYINIT
100
101*ml-- (energy,A,ifb1,bmin,bmax,bfix1,nh1)
102* PYTHIA inizialization, calculation of total and hard cross sections and
103* # of participants and binary sub-collisions at "reference point" (Pb,b=0),
104* tabulation of nuclear thickness function and nuclear overlap function,
105* test of input parametes
106 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
107 double precision numpar,npar0,nbco0
108 external numpar,rhoaa,hfunc3
109 common /pyint7/ sigt(0:6,0:6,0:5)
110 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
111 common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
112 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
113 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
114 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
115 common /hyflow/ ytfl,ylfl,Tf,fpart
116 common /hygeom/ BC
117 common /hythic/ BAB(110),TAB(110),TAAB(110)
118 common /hynup1/ bp,x
119
120*
121
7b7936e9 122c-- save /pyint7/,/pypars/,/pysubs/,/hyjpar/,/hyipar/,/hypyin/,
123c-- > /hyflow/,/hygeom/,/hythic/
b1c2e580 124
125* start HYDJET initialization
126 init=1
127 ipr=1
128
129* set beam paramters
130c-ml ene=energy ! c.m.s. energy per nucleon
131c-ml AW=A ! atomic weight
132
133
134 RA=1.15d0*AW**0.333333d0 ! nuclear radius
135 rzta=1.d0/(1.98d0+0.015d0*A**0.666667d0) ! fraction of protons in nucleus
136 rnta=1.d0-rzta ! fraction of neutrons in nucleus
137c-ml ifb=ifb1 ! centrality flag
138c-ml bfix=bfix1 ! fixed impact parameter
139c-ml nh=nh1 ! mean soft mult. in central PbPb
140
141c-ml ptmin=ckin(3) ! minimum pt of hard scattering
142
143c* printout of HYDJET versioning
144c write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145c >%%%%%%%%%%%%%%%%%%%%%%%%'
146c write(6,*) '%% This is HYDJET version 1.5
147c > %%'
148c write(6,*) '%% Corresponding author: Igor Lokhtin (Igor.Lo
149c >khtin@cern.ch) %%'
150c write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151c >%%%%%%%%%%%%%%%%%%%%%%%%'
152
153* Pythia inizialization for pp collisions
154 call pyinit('cms','p','p',ene)
155
156* calculation of hard scattering cross section
157 if(nhsel.ne.0) then
158
159 mstp(111)=0
160* no printout of Pythia initialization information hereinafter
161 mstp(122)=0
162
163
164* initialize HYINIT with the input parameters
165 write(*,*)' in hyinit AW= ', AW,' energy= ', ene, ' ptmin= ',ptmin
166 write(*,*)' bminh= ',bminh,' bmaxh= ', bmaxh
167 write(*,*)' ifb= ', ifb, ' bfix= ', bfix
168
169
170* Pythia test pp event run
171 do i=1,1000
172 call pyevnt
173 end do
174* hard scattering pp cross section
175 sjpp=pari(1)
176
177* Pythia inizialization for pn collisions
178 call pyinit('cms','p','n',ene)
179* Pythia test pn event run
180 do i=1,1000
181 call pyevnt
182 end do
183* hard scattering pn cross section
184 sjpn=pari(1)
185
186* Pythia inizialization for nn collisions
187 call pyinit('cms','n','n',ene)
188* Pythia test nn event run
189 do i=1,1000
190 call pyevnt
191 end do
192* hard scattering nn cross section
193 sjnn=pari(1)
194
195 sigjet=rzta*rzta*sjpp+rnta*rnta*sjnn+2.d0*rzta*rnta*sjpn
196
197 end if
198
199* total inelastic cross section
200 if(sigin.lt.10.d0.or.sigin.gt.200.d0)
201 > sigin=sigt(0,0,0)-sigt(0,0,1)
202
203* # of nucelons-participants and NN sub-collisions at "reference point" (Pb,b=0)
204 Apb=207.d0
205 Rpb=1.15d0*Apb**0.333333d0
206 EPS=0.005d0
207 Z2=4.d0*Rpb
208 Z1=-1.d0*Z2
209 H=0.01d0*(Z2-Z1)
210 do ib=1,110
211 BC=3.d0*Rpb*(ib-1)/109.d0
212 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS)
213 BAB(ib)=BC
214 TAB(ib)=Apb*RES
215 end do
216 Z1=0.d0
217 Z2=6.28318d0
218 H=0.01d0*(Z2-Z1)
219 bp=0.d0
220 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,TAAPB0,AIH,AIABS)
221
222 npar0=numpar(0.d0) ! Npart(Pb,b=0)
223 nbco0=0.1d0*sigin*TAAPB0 ! Nsub(Pb,b=0)
224
225 init=0
226
227* creation of arrays for tabulation of beam/target nuclear thickness function
228 Z2=4.d0*RA
229 Z1=-1.d0*Z2
230 H=0.01d0*(Z2-Z1)
231 do ib=1,110
232 BC=3.d0*RA*(ib-1)/109.d0
233 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS)
234 BAB(ib)=BC
235 TAB(ib)=AW*RES
236 end do
237
238* creation of arrays for tabulation of nuclear overlap function
239 Z1=0.d0
240 Z2=6.28318d0
241 H=0.01d0*(Z2-Z1)
242 do ib=1,110
243 bp=BAB(ib)
244 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,RES,AIH,AIABS)
245 TAAB(ib)=RES
246 end do
247
248* test of centrality selection
249 if(ifb.eq.0) then
250 if(bfix.lt.0.d0) then
251 write(6,*) 'Impact parameter less than zero!'
252 bfix=0.d0
253 end if
254 if (bfix.gt.3.d0) then
255 write(6,*) 'Impact parameter larger than three nuclear radius!'
256 bfix=3.d0
257 end if
258 else
259 if(bmin.lt.0.d0) then
260 write(6,*) 'Impact parameter less than zero!'
261 bmin=0.d0
262 end if
263 if(bmax.gt.3.d0) then
264 write(6,*) 'Impact parameter larger than three nuclear radius!'
265 bmax=3.d0
266 end if
267c-ml bminh=bmin
268c-ml bmaxh=bmax
269 end if
270
271* test of flow parameter selection
272 if (Tf.lt.0.08d0.or.Tf.gt.0.2d0) Tf=0.1d0 ! freeze-out temperature
273 if (ylfl.lt.0.01d0.or.ylfl.gt.7.d0) ylfl=4.d0 ! longitudinal flow rapidity
274 if (ytfl.lt.0.01d0.or.ytfl.gt.3.d0) ytfl=1.5d0 ! transverse flow rapidity
275 if (fpart.le.0.d0.or.fpart.gt.1.d0) fpart=1.d0 ! fraction of soft multiplicity
276 ! proport. to # of participants
277* test of 'nhsel' selection
278 if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4)
279 > nhsel=0
280
281 return
282 end
283****************************** END HYINIT ***************************
284 SUBROUTINE MYDELTA
285 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
286 common/SERVICEEV/delta,KC,ipdg
287 COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
288 integer KC, PYCOMP
289 real delta
290* allows to transmit delta -limits for BW from PYTHIA to c++ part
291 KC=PYCOMP(ipdg)
292 delta=PMAS(int(KC),3)
293c write(*,*)" ipdg ", ipdg, " KC ",KC," delta ",delta
294 return
295 end
296
297
298********************************* HYEVNT ****************************
299 SUBROUTINE HYEVNT
300* generation of single HYDJET event (soft+hard parts) at given parameters
301 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
302 double precision numpar,npar0,nbco0,npart,nbcol,npart0
303 real p,v, delta
304 external hsin,hftaa,numpar,hyhard,hipsear,pyr,pymass,PYCOMP
305 external ludata
7b7936e9 306c common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
b1c2e580 307 common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
308 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
309 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
310 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
311 common /hyflow/ ytfl,ylfl,Tf,fpart
312 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
313c-ml
314 common /hypart/ppart(10,50000),bmin,bmax,njp
315 COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
316
7b7936e9 317c-- save/lujets/,/hyjets/,/hyipar/,/hyfpar/,/hyflow/,/hyjpar/,/hypyin/
b1c2e580 318
319
320* reset lujets and hyjets arrays before event generation
321
7b7936e9 322 write(*,*)'in hyevnt 0'
b1c2e580 323
324 n=0
325 nhj=0
7b7936e9 326
b1c2e580 327 do ncl=1,150000
328 do j=1,5
7b7936e9 329c p(ncl,j)=0.
b1c2e580 330 phj(ncl,j)=0.d0
7b7936e9 331c v(ncl,j)=0.
b1c2e580 332 vhj(ncl,j)=0.d0
7b7936e9 333c k(ncl,j)=0
b1c2e580 334 khj(ncl,j)=0
335 enddo
336 end do
337*
338 pi=3.14159d0
339
340* generate impact parameter of A-A collision
341 if(ifb.eq.0) then
342 b1=bfix*RA
343 bgen=bfix
344 else
345 call hipsear(fmax1,xmin1)
346 fmax=fmax1
347 xmin=xmin1
348 3 bb1=xmin*pyr(0)+bminh*RA
349 ff1=fmax*pyr(0)
350 fb=hsin(bb1)
351 if(ff1.gt.fb) goto 3
352 b1=bb1
353 bgen=bb1/RA
354 end if
355
7b7936e9 356 write(*,*)'in hyevnt bminh bgen sigin',RA,bminh,bgen,sigin
b1c2e580 357
358* calculate # of nucelons-participants and binary NN sub-collisions
359 npart=numpar(b1) ! Npart(b)
360 npart0=numpar(0.d0) ! Npart(b)
361
362c write(*,*)'in hyevnt npart',npart
363 nbcol=0.1d0*sigin*hftaa(b1) ! Nsub(b)
364c write(*,*)'in hyevnt nbcol',nbcol
365* generate hard parton-parton scatterings (Q>ptmin) 'njet' times
366 njet=0
367 if(nhsel.ne.0) then
368 pjet=sigjet/sigin
369c write(*,*)'in hyevnt pjet',pjet
370 do i=1,int(nbcol)
371 if(pyr(0).lt.pjet) njet=njet+1
372 end do
373c write(*,*)'before hyhard'
374 call hyhard
375 end if
376
377 npyt=nhj
378
379c write(*,*)'in hyevnt pjet njet=', pjet,njet
380c-ml
381* fill array 'ppart'
382
383 do ih=1,nhj
384 if(i.le.50000)then
385 ppart(1,ih)= khj(ih,1) ! pdg ?
386 ppart(2,ih)= khj(ih,2) ! pdg -mother
387 ppart(3,ih)=phj(ih,1) ! px
388 ppart(4,ih)=phj(ih,2) ! py
389 ppart(5,ih)=phj(ih,3) ! pz
390 ppart(6,ih)=phj(ih,4) ! E
391 ppart(7,ih)= vhj(ih,1)
392 ppart(8,ih)= vhj(ih,2)
393 ppart(9,ih)= vhj(ih,3)
394 ppart(10,ih)= vhj(ih,4)
395
396 endif
397 end do
398
399 njp=nhj ! fill number of jet particles for InitialStateBjorken
400c--
401
402* fill array 'lujets'
403
7b7936e9 404c do ih=1,nhj
405c do jh=1,5
406c p(ih,jh)=phj(ih,jh)
407c k(ih,jh)=khj(ih,jh)
408c v(ih,jh)=vhj(ih,jh)
409c end do
410c end do
b1c2e580 411
412c-ml
413c- call luedit(2) ! remove unstable particles and partons
414
415
416
417 return
418 end
419****************************** END HYEVNT **************************
420********************************* HYHARD ***************************
421 SUBROUTINE HYHARD
422* generate 'njet' number of hard parton-parton scatterings with PYTHIA
423 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
424 double precision npar0,nbco0,npart,nbcol,npart0
425 double precision gaus
426
427 INTEGER PYK,PYCHGE,PYCOMP
428 CHARACTER beam*2,targ*2
429 external pydata
430 external pyp,pyr,pyk,pyquen,shad1,gaus
431
432 common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
433 common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
434 COMMON /PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
435 COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
436 COMMON /PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
437 COMMON /PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
438 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
439 common /parimp/ b1, psib1, r0, rb1, rb2, noquen
440 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
441 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
442 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
443 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
7b7936e9 444c-- save /pyjets/,/pypars/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,
445c-- + /hyjets/,/parimp/,/hyjpar/,/hyipar/,/hyfpar/,/hypyin/
b1c2e580 446
447* generate 'njet' PYTHIA events and fill arrays for partons and hadrons
448
449c write(*,*)'in hyhard 0'
450
451 nshad=0
452 noquen=0
453 if(nhsel.eq.1.or.nhsel.eq.3) noquen=1
454 if(njet.ge.1) then
455 mdcy(pycomp(111),1)=0 ! no pi0 decay
456 mdcy(pycomp(310),1)=0 ! no K_S0 decay
457
458 ifbp=0 ! fix impact parameter
459 bfixp=RA*bgen
460 Ap=AW ! atomic weight
461
462 do ihard=1,njet
463 mstp(111)=0
464
465* generate type of NN sub-collision (pp, pn or nn)
466 rand1=pyr(0)
467 if(rand1.lt.rzta) then
468 beam='p'
469 else
470 beam='n'
471 end if
472 rand2=pyr(0)
473 if(rand2.lt.rzta) then
474 targ='p'
475 else
476 targ='n'
477 end if
478 call pyinit('cms',beam,targ,ene)
479c mstj(41)=0 ! vacuum showering off
480
481
482c write(*,*)'in hyhard 1'
483
484 call pyevnt ! generate hard scattering
485
486c write(*,*)'in hyhard 2'
487
488
489* PYQUEN: quenched jets if noquen=0 or non-quenched jets if noquen=1 and ishad=1
490 if(ishad.eq.1.or.nhsel.eq.2.or.nhsel.eq.4)
491 > call pyquen(Ap,ifbp,bfixp)
492
493c-ml coordinate info if we need of it
494
495c write(*,*)'in hyhard pari(21)',pari(21)
496
497 Q=pari(21)
498 x=r0*cos(psib1) !fm
499 y=r0*sin(psib1) !fm
500c write(*,*)'in hyhard ',x,y
501
502 if(Q.gt.0)then
503 tau=1./Q ! 1/GeV
504 tau=tau*0.197 !fm/c
505 rm=0.0d0
506 sig=1.0d0
507
508c write(*,*)'in hyhard tau',tau
509
510 etaLj=gaus(rm,sig) !fm
511
512c write(*,*)'in hyhard after gauss',tau
513 z=tau*sinh(etaLj)
514 t=tau*cosh(etaLj)
515c write(*,*)'in hyhard z t',z,t
516 endif
517
518c write(*,*)'x y z t', x,y,z,t,etaLj,Q
519
520* treatment of "nuclear shadowing" (for Pb, Au, Pd or Ca beams only)
521
522
523 if(ishad.eq.1) then
524 kfh1=abs(k(3,2))
525 kfh2=abs(k(4,2))
526 xh1=pari(33)
527 xh2=pari(34)
528 Q2=pari(22)
529 shad=shad1(kfh1,xh1,Q2,rb1)*shad1(kfh2,xh2,Q2,rb2)
530 if(pyr(0).gt.shad) then
531 nshad=nshad+1
532 goto 53
533 end if
534 end if
535
536
537c write(*,*)'in hyhard 4'
538
539 call pyexec ! hadronization done
540
541c write(*,*)'in hyhard 5'
542
543c-ml
544 call pyedit(2) ! remove partons & leave hadrons
545
546c write(*,*)'in hyhard 6'
547
548
549* fill array of final particles
550 nu=nhj+n
551 if(nu.gt.150000-np) then
552 write(6,*) 'Warning, multiplicity too large! Cut hard part.'
553 goto 52
554 end if
555 nhj=nu
556
557
558 do i=nhj-n+1,nhj
559 ip=i+n-nhj
560c write(*,*)'in hyhard 7'
561
562 do j=1,5
563 phj(i,j)=p(ip,j)
564 end do
565
566 vhj(i,1)=x
567 vhj(i,2)=y
568 vhj(i,3)=z
569 vhj(i,4)=t
570
571 do j=1,5
572 khj(i,j)=k(ip,j)
573 end do
574 do j=3,5
575 kk=khj(i,j)
576 if(kk.gt.0) then
577 khj(i,j)=kk+nhj-n
578 end if
579 end do
580 end do
581
582
583 53 continue
584
585c write(*,*)'in hyhard 8'
586
587 end do
588 52 njet=ihard-1
589
590c write(*,*)'in hyhard 9', njet,nshad
591
592 end if
593 njet=njet-nshad
594
595
596c write(*,*)'in hyhard OUT'
597
598
599 return
600 end
601****************************** END HYHARD **************************
602********************************* HIPSEAR ***************************
603 SUBROUTINE HIPSEAR (fmax,xmin)
604* find maximum and 'sufficient minimum' of differential inelasic AA cross
605* section as a function of impact paramater (xm, fm are outputs)
606 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
607 double precision npar0,nbco0
608 external hsin
609 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
7b7936e9 610c-- save /hyipar/
b1c2e580 611 xmin=(bmaxh-bminh)*RA
612
613c write(*,*)'bmaxh bminh',bmaxh,bminh
614
615 fmax=0.d0
616 do j=1,1000
617 x=bminh*RA+xmin*(j-1)/999.d0
618 f=hsin(x)
619 if(f.gt.fmax) then
620 fmax=f
621 endif
622 end do
623 return
624 end
625****************************** END HIPSEAR **************************
626
627* differential inelastic AA cross section
628 double precision function hsin(x)
629 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
630 external hftaa
631 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
7b7936e9 632c-- save /hyjpar/
b1c2e580 633 br=x
634 hsin=br*(1.d0-dexp(-0.1d0*hftaa(br)*sigin))
635 return
636 end
637
638* number of nucleons-participants at impact parameter b
639 double precision function numpar(c)
640 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
641 external HFUNC1
642 common /hynup1/ bp,x
643 EPS=0.005d0
644 A=0.d0
645 B=6.28318d0
646 H=0.01d0*(B-A)
647 bp=c
648 CALL SIMPA(A,B,H,EPS,1.d-8,HFUNC1,X,RES,AIH,AIABS)
649 numpar=2.d0*RES
650 return
651 end
652*
653 double precision function HFUNC1(x)
654 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
655 double precision npar0,nbco0
656 external HFUNC2
657 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
658 common /hynup1/ bp,xx
7b7936e9 659c-- save /hyipar/
b1c2e580 660 if(init.eq.1) then
661 Rl=Rpb
662 else
663 Rl=RA
664 end if
665 xx=x
666 EPS=0.005d0
667 A=0.d0
668 B=3.d0*Rl
669 H=0.01d0*(B-A)
670 CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC2,Y,RES,AIH,AIABS)
671 HFUNC1=RES
672 return
673 end
674*
675 double precision function HFUNC2(y)
676 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
677 double precision npar0,nbco0
678 external hythik
679 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
680 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
681 common /hynup1/ bp,x
7b7936e9 682c-- save /hyipar/,/hyjpar/
b1c2e580 683 r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x)))
684 r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x)))
685 s=1.d0-dexp(-0.1d0*sigin*hythik(r2))
686 HFUNC2=y*hythik(r1)*s
687 return
688 end
689
690* nuclear overlap function at impact parameter b
691 double precision function hftaa(c)
692 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
693 common /hythic/ BAB(110),TAB(110),TAAB(110)
7b7936e9 694c-- save /hythic/
b1c2e580 695 call parinv(c,BAB,TAAB,110,RES)
696 hftaa=RES
697 return
698 end
699*
700 double precision function HFUNC3(x)
701 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
702 double precision npar0,nbco0
703 external HFUNC4
704 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
705 common /hynup1/ bp,xx
7b7936e9 706c-- save /hyipar/
b1c2e580 707 if(init.eq.1) then
708 Rl=Rpb
709 else
710 Rl=RA
711 end if
712 xx=x
713 EPS=0.005d0
714 A=0.d0
715 B=3.d0*Rl
716 H=0.01d0*(B-A)
717 CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC4,Y,RES,AIH,AIABS)
718 HFUNC3=RES
719 return
720 end
721*
722 double precision function HFUNC4(y)
723 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
724 double precision npar0,nbco0
725 external hythik
726 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
727 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
728 common /hynup1/ bp,x
7b7936e9 729c-- save /hyipar/,/hyjpar/
b1c2e580 730 r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x)))
731 r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x)))
732 HFUNC4=y*hythik(r1)*hythik(r2)
733 return
734 end
735
736* nuclear thickness function
737 double precision function hythik(r)
738 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
739 common /hythic/ BAB(110),TAB(110),TAAB(110)
7b7936e9 740c-- save /hythic/
b1c2e580 741 call parinv(r,BAB,TAB,110,RES)
742 hythik=RES
743 return
744 end
745
746* Wood-Saxon nucleon distrubution
747 double precision function rhoaa(z)
748 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
749 double precision npar0,nbco0
750 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
751 common /hygeom/ BC
7b7936e9 752c-- save /hyipar/,/hygeom/
b1c2e580 753 if(init.eq.1) then
754 Rl=Rpb
755 else
756 Rl=RA
757 end if
758 pi=3.14159d0
759 df=0.54d0
760 r=sqrt(bc*bc+z*z)
761 rho0=3.d0/(4.d0*pi*Rl**3)/(1.d0+(pi*df/Rl)**2)
762 rhoaa=rho0/(1.d0+dexp((r-Rl)/df))
763 return
764 end
765
766* function to calculate nuclear shadowing factor (for Pb, Au, Pd or Ca beams)
767 double precision function shad1(kfh,xbjh,Q2h,r)
768 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
769 double precision npar0,nbco0,nbcol,npart,npart0
770 external ggshad
771 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
772 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
773 common /hyshad/ bbmin,bbmax,inuc
7b7936e9 774c-- save /hyipar/,/hyfpar/,/hyshad/
b1c2e580 775 dimension res(2)
776 kf=kfh
777 xbj=xbjh
778 Q2=Q2h
779 bb=r
780 inuc=0
781 if(AW.gt.205.d0.and.AW.lt.209.d0) inuc=4 ! Pb-206, 207 or 208
782 if(AW.gt.196.d0.and.AW.lt.198.d0) inuc=3 ! Au-197
783 if(AW.gt.109.d0.and.AW.lt.111.d0) inuc=2 ! Pd-110
784 if(AW.gt.39.d0.and.AW.lt.41.d0) inuc=1 ! Ca-40
785 if(inuc.eq.0.and.ipr.eq.1) then
786 write(6,*)
787 > 'Warning! Shadowing is not foreseen for atomic weigth A ='
788 > ,AW
789 write(6,*)'******************************************************
790 >************************'
791 ipr=0
792 end if
793 if(inuc.ne.0) then
794 xbj=max(5.d-5,xbj)
795 xbj=min(0.95d0,xbj)
796 Q2=max(4.d0,Q2)
797 Q2=min(520.d0,Q2)
798 bb=max(0.d0,bb)
799 call ggshad(inuc,xbj,Q2,bb,res,ta)
800 if(kf.eq.21) then
801 shad1=res(1)
802 elseif(kf.eq.1.or.kf.eq.2.or.kf.eq.3) then
803 shad1=res(2)
804 else
805 shad1=1.d0
806 end if
807 else
808 shad1=1.d0
809 end if
810 return
811 end
812
813******************************************************************************
814* The part of the code which follows below, includes nuclear shadowing model *
815* and has been written by Konrad Tywoniuk (Oslo University, Norway) *
816******************************************************************************
817c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$
818c$$$
819c$$$
820c$$$ Shadowing from Glauber-Gribov theory for
821c$$$ gluons and light quarks (u,d,s and their antiquarks).
822c$$$ All Pomeron tree diagrams have been summed (Schwimmer model).
823c$$$ More details about the model in
824c$$$ K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and E.
825c$$$ Zabrodin, Phys. Lett. B 657 (2007) 170
826c$$$
827c$$$ We use FIT B from the H1 parameterization published in
828c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 715
829c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 749
830c$$$
831c$$$ Main routine is GGSHAD which provides the user with the
832c$$$ ratio of nuclear and nucleon parton distribution function
833c$$$ normalized by the atomic number for a given
834c$$$ INUCL: nucleus (1=Ca,2=Pd,3=Au,4=Pb)
835c$$$ X: Bjorken x
836c$$$ Q2: q**2, scale squared
837c$$$ B: impact parameter
838c$$$
839c$$$ Then RES(1): gluon shadowing
840c$$$ RES(2): sea quark shadowing
841c$$$ TA: nuclear profile function
842c$$$
843c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$
844 SUBROUTINE GGSHAD(INUCL,X,Q2,B,RES,TAF)
845 IMPLICIT NONE
846 DOUBLE PRECISION TA(31,4),IMPAR(31),ANUCL(4)
847 DOUBLE PRECISION XB(36),Q2V(13)
848 DOUBLE PRECISION XMAX,XMAXX,Q2MIN,Q2MAX,BMAX
849 DOUBLE PRECISION G(36,13,4),LQ(36,13,4)
850 DOUBLE PRECISION TATMP(31),SHAD(2)
851 DOUBLE PRECISION C(100),D(100),E(100)
852 DOUBLE PRECISION X,Q2,B
853 DOUBLE PRECISION RES(2)
854 DOUBLE PRECISION TAF
855 DOUBLE PRECISION SEVAL
856 INTEGER INUCL,IREAD,TMAX
857 INTEGER I,J,K,IK
858
859 PARAMETER(TMAX=31)
860 PARAMETER(XMAX=0.1)
861 PARAMETER(XMAXX=0.95)
862 PARAMETER(Q2MIN=4.)
863 PARAMETER(Q2MAX=520.)
864 PARAMETER(BMAX=30.)
865
866 COMMON/GG07/XB,Q2V,G,LQ
867
868 DATA ANUCL/40.,110.,197.,206./
869 DATA IMPAR
870 > /0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,6.5,7.,7.5,8.,
871 > 8.5,9.,9.5,10.,10.5,11.,11.5,12.,12.5,13.,13.5,14.,14.5,15./
872 DATA TA
873 > /0.0291662,0.0288636,0.0279354,0.0263044,0.0238341,0.0203565,
874 > 0.0158333,0.0107343,0.00617758,0.00307859,0.00139738,
875 > 0.000603839,0.000254843,0.000106329,4.40965e-05,1.82206e-05,
876 > 7.50949e-06,3.08881e-06,1.26837e-06,5.20082e-07,2.12978e-07,
877 > 8.71154e-08,3.55956e-08,1.45304e-08,5.92622e-09,2.41504e-09,
878 > 9.83429e-10,4.00184e-10,1.62741e-10,6.61407e-11,2.68657e-11,
879 > 0.0153487,0.0152775,0.0150614,0.0146926,0.0141561,0.0134268,
880 > 0.0124651,0.0112115,0.00959258,0.00756907,0.00527514,
881 > 0.00313482,0.00159852,0.000732479,0.000316628,0.000133116,
882 > 5.52508e-05,2.27904e-05,9.3692e-06,3.84349e-06,1.57423e-06,
883 > 6.43962e-07,2.63132e-07,1.07414e-07,4.38089e-08,1.78529e-08,
884 > 7.2699e-09,2.95832e-09,1.20304e-09,4.8894e-10,1.98603e-10,
885 > 0.0105299,0.010498,0.0104017,0.010239,0.010006,0.00969677,
886 > 0.00930205,0.00880781,0.00819264,0.00742455,0.00646058,
887 > 0.00526252,0.00385835,0.00243988,0.00131252,0.000621079,
888 > 0.000272322,0.000114993,4.77312e-05,1.96569e-05,8.06377e-06,
889 > 3.30066e-06,1.34904e-06,5.50751e-07,2.24633e-07,9.15434e-08,
890 > 3.72777e-08,1.51694e-08,6.16885e-09,2.50714e-09,1.01838e-09,
891 > 0.0102179,0.0101879,0.0100975,0.00994467,0.00972606,
892 > 0.00943628,0.0090671,0.00860605,0.00803416,0.00732289,
893 > 0.00643266,0.00532296,0.00400025,0.00261264,0.00144993,
894 > 0.00070098,0.000310867,0.000131949,5.48887e-05,2.26247e-05,
895 > 9.28457e-06,3.80091e-06,1.55358e-06,6.34273e-07,2.58701e-07,
896 > 1.05427e-07,4.29316e-08,1.74701e-08,7.10447e-09,2.88739e-09,
897 > 1.17283e-09/
898 DATA Q2V
899 > /4.000, 6.000, 9.000, 13.500, 20.250, 30.375, 45.562,
900 > 68.344, 102.516, 153.773, 230.660, 345.990, 518.985/
901 DATA XB
902 > /0.99999998E-05, 0.13000000E-04, 0.16899999E-04,
903 > 0.21970000E-04, 0.28561000E-04, 0.37129299E-04,
904 > 0.48268099E-04, 0.62748499E-04, 0.81573096E-04,
905 > 0.10604500E-03, 0.13785800E-03, 0.17921600E-03,
906 > 0.23298099E-03, 0.30287501E-03, 0.39373801E-03,
907 > 0.51185902E-03, 0.66541700E-03, 0.86504198E-03,
908 > 0.11245500E-02, 0.14619200E-02, 0.19005000E-02,
909 > 0.24706500E-02, 0.32118401E-02, 0.41753901E-02,
910 > 0.54280101E-02, 0.70564100E-02, 0.91733299E-02,
911 > 0.11925300E-01, 0.15502900E-01, 0.20153800E-01,
912 > 0.26200000E-01, 0.34059901E-01, 0.44277899E-01,
913 > 0.57561301E-01, 0.74829698E-01, 0.97278602E-01/
914 DATA G
915 > /2.542140,2.457280,2.372200,2.283920,2.195320,2.105540,
916 > 2.014040,1.921480,1.827600,1.733930,1.639990,1.546170,
917 > 1.453220,1.361040,1.268750,1.178500,1.088730,1.000740,
918 > 0.914606,0.830842,0.749315,0.671419,0.596883,0.526173,
919 > 0.460067,0.398489,0.341895,0.289998,0.242628,0.199405,
920 > 0.159625,0.122521,0.087252,0.053699,0.022804,0.000414,
921 > 1.825600,1.768440,1.712550,1.654660,1.595830,1.533980,
922 > 1.472860,1.409910,1.346060,1.282030,1.217520,1.153050,
923 > 1.089030,1.025930,0.962709,0.900094,0.837787,0.776733,
924 > 0.716072,0.656570,0.598295,0.541607,0.486661,0.433846,
925 > 0.383421,0.335585,0.290685,0.248542,0.209305,0.172614,
926 > 0.138166,0.105552,0.074246,0.044592,0.017935,0.000291,
927 > 1.445600,1.405920,1.364150,1.322050,1.276900,1.232020,
928 > 1.184670,1.136620,1.086950,1.037370,0.987225,0.937488,
929 > 0.888244,0.839061,0.790370,0.741798,0.693426,0.645511,
930 > 0.598298,0.551516,0.505500,0.460400,0.416326,0.373488,
931 > 0.332281,0.292614,0.254859,0.218944,0.184976,0.152759,
932 > 0.122077,0.092747,0.064477,0.037873,0.014527,0.000216,
933 > 1.209700,1.179300,1.146560,1.112700,1.077790,1.041400,
934 > 1.003150,0.964598,0.923825,0.882423,0.841304,0.800381,
935 > 0.759486,0.718947,0.678587,0.638579,0.598671,0.558978,
936 > 0.519524,0.480617,0.442116,0.404231,0.367057,0.330599,
937 > 0.295291,0.261100,0.228220,0.196605,0.166367,0.137331,
938 > 0.109496,0.082677,0.056835,0.032696,0.012013,0.000165,
939 > 1.049390,1.024810,0.998839,0.971705,0.942633,0.912223,
940 > 0.880504,0.847236,0.812716,0.777477,0.741670,0.706578,
941 > 0.671021,0.636183,0.601384,0.566923,0.532196,0.497956,
942 > 0.463884,0.430069,0.396612,0.363536,0.330915,0.298900,
943 > 0.267627,0.237155,0.207697,0.179169,0.151566,0.124961,
944 > 0.099266,0.074436,0.050544,0.028462,0.010013,0.000127,
945 > 0.943086,0.922388,0.900700,0.877793,0.853134,0.826808,
946 > 0.798937,0.769284,0.738709,0.706829,0.674818,0.643041,
947 > 0.611250,0.580168,0.548982,0.517908,0.486720,0.456021,
948 > 0.425312,0.394899,0.364616,0.334682,0.305173,0.276042,
949 > 0.247465,0.219584,0.192474,0.166033,0.140391,0.115481,
950 > 0.091347,0.068001,0.045636,0.025180,0.008510,0.000101,
951 > 0.858923,0.841883,0.823827,0.803952,0.782124,0.759447,
952 > 0.734659,0.708141,0.680281,0.651690,0.622372,0.593451,
953 > 0.564610,0.536045,0.507521,0.479248,0.450845,0.422731,
954 > 0.394687,0.366774,0.339095,0.311580,0.284350,0.257527,
955 > 0.231101,0.205195,0.179929,0.155178,0.131077,0.107555,
956 > 0.084712,0.062629,0.041557,0.022497,0.007326,0.000082,
957 > 0.790086,0.776082,0.760408,0.743282,0.724441,0.703812,
958 > 0.681673,0.657893,0.632354,0.605918,0.578976,0.552364,
959 > 0.525820,0.499665,0.473306,0.447125,0.421027,0.394926,
960 > 0.369017,0.343250,0.317596,0.292061,0.266805,0.241770,
961 > 0.217121,0.192830,0.169123,0.145809,0.122984,0.100685,
962 > 0.078967,0.057983,0.038065,0.020236,0.006364,0.000067,
963 > 0.732957,0.720881,0.707633,0.692613,0.676200,0.657914,
964 > 0.637626,0.615846,0.592373,0.567813,0.542877,0.518036,
965 > 0.493386,0.468899,0.444596,0.420291,0.395908,0.371641,
966 > 0.347423,0.323358,0.299396,0.275511,0.251830,0.228313,
967 > 0.205146,0.182274,0.159813,0.137693,0.115983,0.094703,
968 > 0.073974,0.053979,0.035068,0.018327,0.005580,0.000056,
969 > 0.684126,0.674389,0.662432,0.649473,0.634783,0.618458,
970 > 0.600173,0.580108,0.558464,0.535421,0.512051,0.488760,
971 > 0.465717,0.442837,0.419997,0.397309,0.374393,0.351598,
972 > 0.328943,0.306305,0.283685,0.261210,0.238875,0.216696,
973 > 0.194737,0.173035,0.151692,0.130607,0.109864,0.089483,
974 > 0.069628,0.050489,0.032492,0.016714,0.004940,0.000047,
975 > 0.642684,0.634012,0.624311,0.612687,0.599751,0.584885,
976 > 0.568165,0.549716,0.529352,0.507582,0.485673,0.463844,
977 > 0.442031,0.420480,0.398892,0.377446,0.355834,0.334383,
978 > 0.312903,0.291520,0.270128,0.248790,0.227632,0.206518,
979 > 0.185624,0.164939,0.144527,0.124348,0.104415,0.084843,
980 > 0.065754,0.047403,0.030217,0.015305,0.004396,0.000040,
981 > 0.606784,0.599501,0.591216,0.581053,0.569262,0.555644,
982 > 0.540298,0.522979,0.503888,0.483449,0.462725,0.441941,
983 > 0.421308,0.400965,0.380515,0.360166,0.339667,0.319306,
984 > 0.298896,0.278531,0.258224,0.237932,0.217694,0.197574,
985 > 0.177567,0.157731,0.138168,0.118766,0.099587,0.080708,
986 > 0.062310,0.044656,0.028210,0.014081,0.003935,0.000035,
987 > 0.575250,0.569055,0.561834,0.552756,0.542358,0.530173,
988 > 0.515880,0.499697,0.481592,0.462197,0.442485,0.422845,
989 > 0.403137,0.383704,0.364332,0.344943,0.325467,0.306020,
990 > 0.286496,0.267106,0.247644,0.228251,0.208881,0.189557,
991 > 0.170377,0.151316,0.132475,0.113748,0.095226,0.076992,
992 > 0.059214,0.042195,0.026424,0.013004,0.003541,0.000030,
993 > 2.510860,2.424890,2.338800,2.251430,2.161030,2.070080,
994 > 1.976500,1.884380,1.789520,1.695430,1.600900,1.507110,
995 > 1.413980,1.321430,1.229430,1.138950,1.049380,0.961914,
996 > 0.875863,0.792403,0.711728,0.634220,0.560088,0.490193,
997 > 0.424574,0.363647,0.307681,0.256400,0.209654,0.167105,
998 > 0.128273,0.092911,0.060892,0.033178,0.011804,0.000172,
999 > 1.797890,1.740890,1.683730,1.625330,1.564220,1.502940,
1000 > 1.439440,1.375560,1.311540,1.247060,1.182260,1.117820,
1001 > 1.053940,0.990598,0.927193,0.864982,0.802868,0.741596,
1002 > 0.681232,0.621931,0.563990,0.507817,0.452982,0.400714,
1003 > 0.350809,0.303514,0.259139,0.217505,0.178849,0.142987,
1004 > 0.109689,0.079022,0.051121,0.027212,0.009208,0.000121,
1005 > 1.421860,1.380980,1.338550,1.294150,1.249120,1.201850,
1006 > 1.153940,1.105280,1.054910,1.004680,0.954677,0.905086,
1007 > 0.855141,0.806467,0.757283,0.708817,0.660443,0.613233,
1008 > 0.565935,0.519307,0.473770,0.428741,0.385002,0.342676,
1009 > 0.301797,0.262591,0.225442,0.189974,0.156683,0.125290,
1010 > 0.095921,0.068630,0.043899,0.022874,0.007411,0.000090,
1011 > 1.187690,1.155530,1.122280,1.087030,1.051020,1.013270,
1012 > 0.973852,0.933953,0.892905,0.851713,0.810455,0.769192,
1013 > 0.728379,0.687798,0.647423,0.607517,0.567475,0.528048,
1014 > 0.488850,0.450119,0.411897,0.374353,0.337466,0.301396,
1015 > 0.266522,0.232713,0.200387,0.169339,0.139725,0.111661,
1016 > 0.085284,0.060612,0.038324,0.019572,0.006096,0.000069,
1017 > 1.028240,1.002610,0.975714,0.947100,0.916605,0.885287,
1018 > 0.852214,0.817988,0.782800,0.747390,0.711939,0.676325,
1019 > 0.641102,0.606200,0.571466,0.536994,0.502546,0.468438,
1020 > 0.434547,0.400992,0.367663,0.334910,0.302578,0.270890,
1021 > 0.240102,0.210041,0.181153,0.153175,0.126406,0.100826,
1022 > 0.076605,0.054074,0.033774,0.016896,0.005055,0.000053,
1023 > 0.923239,0.901357,0.878307,0.853928,0.827994,0.800448,
1024 > 0.770833,0.740864,0.709208,0.677436,0.645505,0.613815,
1025 > 0.582211,0.551121,0.519863,0.488847,0.457919,0.427291,
1026 > 0.396728,0.366593,0.336576,0.306991,0.277689,0.248974,
1027 > 0.220914,0.193485,0.166892,0.141077,0.116297,0.092513,
1028 > 0.069960,0.049024,0.030239,0.014840,0.004277,0.000042,
1029 > 0.840109,0.821572,0.802221,0.780737,0.758007,0.733570,
1030 > 0.707391,0.680060,0.651741,0.622598,0.593525,0.564790,
1031 > 0.535923,0.507607,0.479180,0.450770,0.422703,0.394808,
1032 > 0.366922,0.339201,0.311718,0.284555,0.257727,0.231240,
1033 > 0.205291,0.179821,0.155134,0.131145,0.107877,0.085551,
1034 > 0.064423,0.044815,0.027334,0.013173,0.003668,0.000034,
1035 > 0.771898,0.756403,0.739454,0.720917,0.700462,0.678667,
1036 > 0.655360,0.630295,0.604073,0.577412,0.550894,0.524237,
1037 > 0.497613,0.471552,0.445325,0.419311,0.393354,0.367582,
1038 > 0.341816,0.316292,0.290960,0.265746,0.240859,0.216190,
1039 > 0.191974,0.168255,0.145100,0.122534,0.100715,0.079595,
1040 > 0.059666,0.041209,0.024869,0.011778,0.003176,0.000028,
1041 > 0.715432,0.701697,0.687335,0.671043,0.652774,0.633196,
1042 > 0.611618,0.588703,0.564555,0.539941,0.515046,0.490273,
1043 > 0.465666,0.441417,0.417189,0.393083,0.368842,0.344918,
1044 > 0.320853,0.296990,0.273346,0.249853,0.226554,0.203403,
1045 > 0.180727,0.158303,0.136484,0.115160,0.094448,0.074509,
1046 > 0.055551,0.038122,0.022769,0.010610,0.002775,0.000023,
1047 > 0.666960,0.655627,0.643169,0.628497,0.612317,0.594375,
1048 > 0.574670,0.553372,0.531024,0.507805,0.484742,0.461630,
1049 > 0.438595,0.415762,0.393171,0.370479,0.347803,0.325380,
1050 > 0.302957,0.280541,0.258199,0.236074,0.214141,0.192375,
1051 > 0.170880,0.149710,0.129011,0.108739,0.089084,0.069999,
1052 > 0.052020,0.035461,0.020980,0.009629,0.002450,0.000020,
1053 > 0.626089,0.616427,0.605000,0.592191,0.577470,0.561221,
1054 > 0.543098,0.523316,0.502350,0.480560,0.458749,0.436948,
1055 > 0.415413,0.393933,0.372648,0.351193,0.329853,0.308661,
1056 > 0.287323,0.266214,0.245090,0.224191,0.203371,0.182736,
1057 > 0.162324,0.142202,0.122478,0.103078,0.084237,0.066064,
1058 > 0.048902,0.033114,0.019409,0.008778,0.002175,0.000017,
1059 > 0.590664,0.582220,0.572479,0.560768,0.547507,0.532490,
1060 > 0.515657,0.497052,0.477369,0.456893,0.436199,0.415686,
1061 > 0.395192,0.374819,0.354542,0.334339,0.314120,0.293919,
1062 > 0.273897,0.253708,0.233747,0.213772,0.193977,0.174220,
1063 > 0.154737,0.135524,0.116625,0.098102,0.080014,0.062548,
1064 > 0.046083,0.031041,0.018032,0.008042,0.001942,0.000014,
1065 > 0.559453,0.552473,0.543647,0.533124,0.521070,0.507229,
1066 > 0.491452,0.474143,0.455506,0.435945,0.416256,0.396977,
1067 > 0.377326,0.358052,0.338792,0.319598,0.300262,0.281122,
1068 > 0.261883,0.242706,0.223579,0.204566,0.185577,0.166710,
1069 > 0.148088,0.129570,0.111434,0.093593,0.076187,0.059415,
1070 > 0.043600,0.029214,0.016816,0.007398,0.001744,0.000012,
1071 > 2.494430,2.407750,2.320680,2.230570,2.139410,2.045820,
1072 > 1.954300,1.859160,1.764370,1.669990,1.574490,1.481120,
1073 > 1.386830,1.294720,1.202200,1.111670,1.021960,0.934479,
1074 > 0.848525,0.765125,0.684539,0.607178,0.533502,0.463759,
1075 > 0.398700,0.338188,0.282593,0.231728,0.185473,0.143600,
1076 > 0.105860,0.072262,0.043334,0.020543,0.005838,0.000062,
1077 > 1.781830,1.724140,1.666660,1.605560,1.544230,1.481280,
1078 > 1.418440,1.353130,1.288210,1.223120,1.158330,1.093950,
1079 > 1.029490,0.965787,0.902413,0.840091,0.777784,0.716620,
1080 > 0.656180,0.597156,0.539338,0.483075,0.428748,0.376710,
1081 > 0.327154,0.280125,0.236105,0.194961,0.156820,0.121629,
1082 > 0.089495,0.060685,0.035892,0.016616,0.004508,0.000044,
1083 > 1.407480,1.365110,1.320730,1.276090,1.228930,1.181910,
1084 > 1.132930,1.082570,1.032770,0.982143,0.932084,0.882066,
1085 > 0.832465,0.783250,0.734272,0.685821,0.637505,0.589736,
1086 > 0.542652,0.496248,0.450545,0.405906,0.362451,0.320187,
1087 > 0.279650,0.240757,0.203901,0.169026,0.136216,0.105650,
1088 > 0.077490,0.052176,0.030469,0.013808,0.003599,0.000032,
1089 > 1.173640,1.140320,1.105880,1.069540,1.032320,0.993897,
1090 > 0.953913,0.913021,0.871617,0.829958,0.788517,0.747238,
1091 > 0.706258,0.666094,0.625534,0.585422,0.545512,0.506014,
1092 > 0.466972,0.428374,0.390208,0.352773,0.316057,0.280282,
1093 > 0.245648,0.212195,0.180160,0.149604,0.120645,0.093445,
1094 > 0.068267,0.045631,0.026330,0.011698,0.002940,0.000025,
1095 > 1.014900,0.987942,0.959684,0.930270,0.899044,0.866678,
1096 > 0.832642,0.797620,0.762091,0.726232,0.690987,0.655077,
1097 > 0.620091,0.585191,0.550543,0.515986,0.481432,0.447342,
1098 > 0.413584,0.380034,0.346933,0.314295,0.282265,0.250834,
1099 > 0.220232,0.190535,0.162003,0.134530,0.108400,0.083726,
1100 > 0.060878,0.040604,0.022978,0.010002,0.002422,0.000019,
1101 > 0.909560,0.887017,0.862965,0.837572,0.810211,0.781875,
1102 > 0.751746,0.720666,0.688703,0.656776,0.624847,0.593122,
1103 > 0.561398,0.530391,0.499138,0.468280,0.437229,0.406831,
1104 > 0.376391,0.346350,0.316458,0.286984,0.257979,0.229524,
1105 > 0.201720,0.174606,0.148473,0.123248,0.099111,0.076308,
1106 > 0.055184,0.036277,0.020402,0.008711,0.002038,0.000015,
1107 > 0.827323,0.808240,0.787194,0.764785,0.741102,0.715710,
1108 > 0.688724,0.660472,0.631563,0.602307,0.573344,0.544275,
1109 > 0.515643,0.487203,0.458901,0.430660,0.402505,0.374606,
1110 > 0.346855,0.319408,0.292040,0.265141,0.238553,0.212326,
1111 > 0.186666,0.161641,0.137398,0.113940,0.091458,0.070180,
1112 > 0.050481,0.032923,0.018295,0.007672,0.001739,0.000012,
1113 > 0.759389,0.742705,0.724897,0.705170,0.684041,0.661178,
1114 > 0.636646,0.610897,0.584460,0.557556,0.530715,0.504308,
1115 > 0.477779,0.451618,0.425597,0.399640,0.373702,0.348041,
1116 > 0.322356,0.297020,0.271779,0.246872,0.222152,0.197859,
1117 > 0.173964,0.150642,0.127976,0.105994,0.084903,0.064934,
1118 > 0.046463,0.030071,0.016527,0.006814,0.001499,0.000010,
1119 > 0.703019,0.688814,0.673059,0.655722,0.636535,0.615667,
1120 > 0.593245,0.569605,0.545086,0.520194,0.495411,0.470863,
1121 > 0.446201,0.421921,0.397793,0.373652,0.349528,0.325620,
1122 > 0.301845,0.278188,0.254649,0.231349,0.208318,0.185538,
1123 > 0.163153,0.141239,0.119901,0.099172,0.079281,0.060422,
1124 > 0.043027,0.027647,0.015032,0.006099,0.001305,0.000008,
1125 > 0.655549,0.643053,0.628966,0.613794,0.596440,0.577313,
1126 > 0.556691,0.534653,0.511805,0.488542,0.465448,0.442353,
1127 > 0.419418,0.396807,0.374153,0.351526,0.328993,0.306561,
1128 > 0.284259,0.262079,0.240000,0.218086,0.196361,0.174927,
1129 > 0.153790,0.133102,0.112916,0.093252,0.074414,0.056530,
1130 > 0.040069,0.025572,0.013769,0.005504,0.001148,0.000007,
1131 > 0.614886,0.603957,0.591709,0.577655,0.561834,0.544400,
1132 > 0.525193,0.504891,0.483435,0.461525,0.439729,0.418166,
1133 > 0.396478,0.375064,0.353725,0.332565,0.311308,0.290188,
1134 > 0.269125,0.248171,0.227323,0.206600,0.186070,0.165736,
1135 > 0.145686,0.125993,0.106807,0.088112,0.070133,0.053112,
1136 > 0.037464,0.023751,0.012665,0.004991,0.001015,0.000006,
1137 > 0.579569,0.570010,0.559085,0.546465,0.532125,0.516026,
1138 > 0.498089,0.478853,0.458681,0.438173,0.417404,0.396958,
1139 > 0.376499,0.356301,0.336115,0.316014,0.295930,0.275886,
1140 > 0.255871,0.236010,0.216223,0.196524,0.176991,0.157626,
1141 > 0.138515,0.119742,0.101404,0.083545,0.066339,0.050078,
1142 > 0.035167,0.022152,0.011703,0.004550,0.000904,0.000005,
1143 > 0.548753,0.540580,0.530520,0.519276,0.506004,0.491061,
1144 > 0.474436,0.456132,0.436996,0.417488,0.397755,0.378421,
1145 > 0.359085,0.339851,0.320700,0.301518,0.282357,0.263296,
1146 > 0.244273,0.225311,0.206421,0.187624,0.168959,0.150424,
1147 > 0.132145,0.114172,0.096557,0.079446,0.062967,0.047378,
1148 > 0.033126,0.020736,0.010987,0.004166,0.000809,0.000004,
1149 > 2.486670,2.401390,2.313150,2.223100,2.132210,2.040620,
1150 > 1.946240,1.852960,1.758660,1.663790,1.569300,1.475400,
1151 > 1.381630,1.288980,1.197380,1.106900,1.017420,0.930068,
1152 > 0.844430,0.760955,0.680549,0.603582,0.530075,0.460570,
1153 > 0.395641,0.335338,0.279919,0.229212,0.183101,0.141385,
1154 > 0.103826,0.070754,0.041898,0.019580,0.005430,0.000056,
1155 > 1.776150,1.719900,1.661320,1.600440,1.539110,1.477010,
1156 > 1.412790,1.348550,1.283220,1.218900,1.153830,1.088770,
1157 > 1.025080,0.961812,0.898582,0.836002,0.773940,0.712963,
1158 > 0.652650,0.593610,0.536099,0.480001,0.425916,0.373940,
1159 > 0.324454,0.277574,0.233716,0.192654,0.154624,0.119622,
1160 > 0.087676,0.059112,0.034659,0.015816,0.004189,0.000039,
1161 > 1.403340,1.361000,1.316930,1.271610,1.225380,1.177160,
1162 > 1.128240,1.078580,1.028640,0.977947,0.928063,0.878521,
1163 > 0.828550,0.779706,0.730569,0.682199,0.633948,0.586522,
1164 > 0.539401,0.493171,0.447639,0.403145,0.359673,0.317729,
1165 > 0.277222,0.238454,0.201716,0.166913,0.134251,0.103840,
1166 > 0.075864,0.050772,0.029388,0.013129,0.003341,0.000029,
1167 > 1.169950,1.136380,1.102100,1.065820,1.028190,0.989956,
1168 > 0.950189,0.908974,0.867893,0.826368,0.784940,0.743815,
1169 > 0.702872,0.662550,0.622187,0.582284,0.542431,0.503183,
1170 > 0.464112,0.425583,0.387521,0.350185,0.313576,0.277943,
1171 > 0.243398,0.210036,0.178125,0.147636,0.118821,0.091774,
1172 > 0.066794,0.044372,0.025371,0.011112,0.002727,0.000022,
1173 > 1.011370,0.984724,0.956435,0.926814,0.895443,0.862964,
1174 > 0.829274,0.794479,0.758875,0.722999,0.687335,0.652160,
1175 > 0.616882,0.582106,0.547373,0.512881,0.478596,0.444627,
1176 > 0.410813,0.377446,0.344502,0.311951,0.279953,0.248594,
1177 > 0.218124,0.188529,0.160091,0.132732,0.106696,0.082174,
1178 > 0.059499,0.039204,0.022122,0.009492,0.002245,0.000017,
1179 > 0.906689,0.883943,0.860162,0.834079,0.807067,0.778589,
1180 > 0.748446,0.717460,0.685716,0.653595,0.621726,0.590095,
1181 > 0.558631,0.527321,0.496368,0.465511,0.434680,0.404249,
1182 > 0.373821,0.343908,0.314123,0.284707,0.255819,0.227414,
1183 > 0.199673,0.172714,0.146637,0.121510,0.097519,0.074849,
1184 > 0.053902,0.035215,0.019623,0.008259,0.001888,0.000014,
1185 > 0.824084,0.805081,0.784198,0.761980,0.737877,0.712401,
1186 > 0.685185,0.657419,0.628573,0.599438,0.570456,0.541554,
1187 > 0.512854,0.484478,0.456174,0.428019,0.399930,0.372215,
1188 > 0.344523,0.317092,0.289826,0.262998,0.236396,0.210320,
1189 > 0.184715,0.159801,0.135667,0.112306,0.089938,0.068791,
1190 > 0.049274,0.031936,0.017585,0.007270,0.001610,0.000011,
1191 > 0.756740,0.739909,0.722048,0.702184,0.681164,0.658293,
1192 > 0.633526,0.608140,0.581469,0.554572,0.527904,0.501517,
1193 > 0.475145,0.448970,0.422970,0.397121,0.371248,0.345564,
1194 > 0.320059,0.294788,0.269623,0.244755,0.220148,0.195923,
1195 > 0.172134,0.148871,0.126312,0.104417,0.083450,0.063624,
1196 > 0.045337,0.029156,0.015873,0.006450,0.001387,0.000009,
1197 > 0.700680,0.686067,0.670562,0.652740,0.633579,0.612848,
1198 > 0.590476,0.566706,0.542207,0.517556,0.492696,0.468105,
1199 > 0.443562,0.419350,0.395237,0.371204,0.347185,0.323339,
1200 > 0.299543,0.276011,0.252614,0.229355,0.206336,0.183677,
1201 > 0.161386,0.139527,0.118283,0.097686,0.077903,0.059183,
1202 > 0.041956,0.026789,0.014429,0.005771,0.001207,0.000008,
1203 > 0.652991,0.640812,0.626563,0.610980,0.593505,0.574586,
1204 > 0.553939,0.531965,0.509196,0.485938,0.462796,0.439793,
1205 > 0.416876,0.394269,0.371716,0.349126,0.326720,0.304355,
1206 > 0.282067,0.259966,0.237971,0.216140,0.194490,0.173138,
1207 > 0.152089,0.131483,0.111379,0.091839,0.073081,0.055347,
1208 > 0.039055,0.024764,0.013208,0.005204,0.001061,0.000006,
1209 > 0.612493,0.601527,0.589260,0.574987,0.559293,0.541793,
1210 > 0.522714,0.502071,0.480580,0.459021,0.437202,0.415483,
1211 > 0.394052,0.372703,0.351511,0.330271,0.309080,0.288012,
1212 > 0.267062,0.246156,0.225353,0.204703,0.184235,0.163954,
1213 > 0.144013,0.124421,0.105317,0.086709,0.068861,0.051973,
1214 > 0.036498,0.022991,0.012143,0.004716,0.000938,0.000005,
1215 > 0.577224,0.567718,0.556768,0.543992,0.529614,0.513264,
1216 > 0.495593,0.476332,0.456118,0.435643,0.414958,0.394475,
1217 > 0.374159,0.353948,0.333917,0.313860,0.293761,0.273842,
1218 > 0.253865,0.234066,0.214285,0.194711,0.175202,0.155932,
1219 > 0.136903,0.118198,0.099962,0.082192,0.065107,0.048988,
1220 > 0.034247,0.021430,0.011215,0.004297,0.000835,0.000005,
1221 > 0.546580,0.538342,0.528410,0.516932,0.503593,0.488594,
1222 > 0.471844,0.453700,0.434588,0.415081,0.395496,0.376054,
1223 > 0.356702,0.337597,0.318488,0.299383,0.280225,0.261271,
1224 > 0.242274,0.223404,0.204576,0.185839,0.167243,0.148787,
1225 > 0.130586,0.112676,0.095184,0.078140,0.061781,0.046329,
1226 > 0.032244,0.020050,0.010400,0.003932,0.000747,0.000004/
1227 DATA LQ
1228 > /0.514104,0.515832,0.517091,0.517536,0.517832,0.517221,
1229 > 0.516574,0.515224,0.513287,0.510858,0.508000,0.504917,
1230 > 0.501246,0.497403,0.492983,0.488152,0.482358,0.476122,
1231 > 0.468649,0.460339,0.450645,0.439812,0.427166,0.412836,
1232 > 0.396439,0.377753,0.356388,0.331566,0.302870,0.269780,
1233 > 0.231439,0.187668,0.138736,0.086850,0.036211,0.000604,
1234 > 0.579741,0.577476,0.575013,0.571386,0.567464,0.562828,
1235 > 0.557414,0.551410,0.544332,0.536891,0.528953,0.521049,
1236 > 0.512612,0.503967,0.494945,0.485619,0.475531,0.464958,
1237 > 0.453985,0.442176,0.429817,0.416392,0.402026,0.386267,
1238 > 0.369251,0.350413,0.329631,0.306106,0.279316,0.248539,
1239 > 0.213030,0.172492,0.127054,0.078865,0.032177,0.000513,
1240 > 0.617205,0.612889,0.607448,0.601280,0.593987,0.585725,
1241 > 0.577239,0.567645,0.557436,0.546468,0.535124,0.523690,
1242 > 0.511913,0.500192,0.488271,0.476132,0.463482,0.450653,
1243 > 0.437460,0.423920,0.409767,0.395178,0.379804,0.363620,
1244 > 0.346492,0.327958,0.307882,0.285606,0.260283,0.231441,
1245 > 0.198228,0.160224,0.117620,0.072457,0.029003,0.000445,
1246 > 0.635918,0.629589,0.622092,0.614001,0.604751,0.594717,
1247 > 0.583683,0.572041,0.559188,0.546115,0.532523,0.518899,
1248 > 0.505201,0.491782,0.477675,0.463752,0.449738,0.435623,
1249 > 0.421062,0.406484,0.391587,0.376279,0.360631,0.344349,
1250 > 0.327299,0.309280,0.289970,0.268596,0.244665,0.217365,
1251 > 0.185991,0.150104,0.109826,0.067166,0.026423,0.000392,
1252 > 0.642211,0.634584,0.625818,0.616705,0.605951,0.594718,
1253 > 0.582351,0.569037,0.554986,0.540276,0.525303,0.510521,
1254 > 0.495596,0.480637,0.465827,0.451058,0.435872,0.421066,
1255 > 0.405880,0.390796,0.375397,0.359903,0.344185,0.327987,
1256 > 0.311213,0.293614,0.274949,0.254454,0.231575,0.205516,
1257 > 0.175652,0.141477,0.103138,0.062618,0.024219,0.000348,
1258 > 0.639180,0.631092,0.621893,0.612137,0.600935,0.589006,
1259 > 0.576122,0.561754,0.547048,0.531226,0.515512,0.499878,
1260 > 0.484089,0.468662,0.453258,0.437735,0.422292,0.407026,
1261 > 0.391544,0.376258,0.360735,0.345208,0.329470,0.313579,
1262 > 0.297158,0.280066,0.261996,0.242200,0.220305,0.195335,
1263 > 0.166765,0.134023,0.097367,0.058717,0.022356,0.000312,
1264 > 0.631156,0.622788,0.613580,0.603527,0.592277,0.579779,
1265 > 0.566484,0.551813,0.536423,0.520508,0.504251,0.488021,
1266 > 0.472093,0.456329,0.440551,0.424908,0.409169,0.393745,
1267 > 0.378229,0.362803,0.347473,0.331992,0.316442,0.300771,
1268 > 0.284762,0.268174,0.250682,0.231645,0.210500,0.186498,
1269 > 0.158976,0.127569,0.092372,0.055348,0.020770,0.000282,
1270 > 0.619698,0.611519,0.602336,0.591943,0.580814,0.568514,
1271 > 0.554982,0.540509,0.524970,0.508623,0.492412,0.476046,
1272 > 0.459935,0.444036,0.428166,0.412462,0.396815,0.381305,
1273 > 0.365892,0.350546,0.335395,0.320121,0.304877,0.289437,
1274 > 0.273841,0.257707,0.240774,0.222305,0.201894,0.178726,
1275 > 0.152209,0.121898,0.087984,0.052399,0.019398,0.000257,
1276 > 0.606343,0.598499,0.589365,0.579821,0.568434,0.556327,
1277 > 0.542872,0.528454,0.512900,0.496651,0.480390,0.464149,
1278 > 0.447972,0.432310,0.416261,0.400748,0.385116,0.369727,
1279 > 0.354542,0.339434,0.324280,0.309328,0.294433,0.279362,
1280 > 0.264086,0.248363,0.231853,0.214025,0.194261,0.171833,
1281 > 0.146197,0.116876,0.084113,0.049805,0.018206,0.000236,
1282 > 0.591896,0.584256,0.575797,0.566248,0.555494,0.543797,
1283 > 0.530664,0.516396,0.501105,0.484951,0.468541,0.452368,
1284 > 0.436358,0.420764,0.405075,0.389617,0.374238,0.359020,
1285 > 0.343972,0.329043,0.314268,0.299607,0.284884,0.270198,
1286 > 0.255337,0.239949,0.223985,0.206644,0.187466,0.165701,
1287 > 0.140843,0.112412,0.080664,0.047512,0.017165,0.000218,
1288 > 0.576697,0.569585,0.561998,0.552972,0.542790,0.531174,
1289 > 0.518335,0.504287,0.489160,0.473285,0.457139,0.441230,
1290 > 0.425520,0.409903,0.394534,0.379177,0.364141,0.349046,
1291 > 0.334214,0.319602,0.305086,0.290690,0.276241,0.261890,
1292 > 0.247381,0.232411,0.216817,0.199969,0.181340,0.160182,
1293 > 0.135997,0.108375,0.077551,0.045446,0.016236,0.000202,
1294 > 0.561779,0.555491,0.548244,0.539719,0.529917,0.518731,
1295 > 0.506334,0.492791,0.477798,0.462049,0.446208,0.430534,
1296 > 0.414938,0.399619,0.384500,0.369572,0.354611,0.339912,
1297 > 0.325296,0.310898,0.296649,0.282465,0.268439,0.254375,
1298 > 0.240158,0.225578,0.210376,0.193910,0.175777,0.155185,
1299 > 0.131618,0.104727,0.074745,0.043586,0.015406,0.000188,
1300 > 0.547406,0.541522,0.534772,0.526637,0.517320,0.506665,
1301 > 0.494929,0.481541,0.466958,0.451581,0.435977,0.420382,
1302 > 0.405163,0.390039,0.375242,0.360562,0.345699,0.331352,
1303 > 0.316934,0.302875,0.288897,0.275075,0.261211,0.247494,
1304 > 0.233546,0.219335,0.204477,0.188465,0.170694,0.150613,
1305 > 0.127623,0.101395,0.072188,0.041898,0.014659,0.000176,
1306 > 0.510949,0.512461,0.513125,0.513775,0.513627,0.512953,
1307 > 0.511676,0.510153,0.508026,0.505586,0.502488,0.499514,
1308 > 0.495918,0.491777,0.487221,0.482318,0.476322,0.469865,
1309 > 0.462254,0.453539,0.443579,0.431957,0.418405,0.402714,
1310 > 0.384513,0.363211,0.338441,0.309525,0.276023,0.237418,
1311 > 0.194095,0.146991,0.098692,0.053921,0.018675,0.000251,
1312 > 0.573840,0.571539,0.568660,0.564636,0.560125,0.555429,
1313 > 0.549250,0.542654,0.535827,0.528339,0.520207,0.512007,
1314 > 0.503505,0.494854,0.485721,0.476289,0.466253,0.455649,
1315 > 0.444534,0.432471,0.419685,0.405709,0.390526,0.373656,
1316 > 0.355038,0.334208,0.310578,0.283443,0.252376,0.217005,
1317 > 0.177144,0.133976,0.089700,0.048645,0.016531,0.000213,
1318 > 0.610274,0.604837,0.598760,0.591967,0.584049,0.575825,
1319 > 0.566634,0.556589,0.545957,0.534847,0.523533,0.511998,
1320 > 0.500328,0.488702,0.476646,0.464532,0.451851,0.438900,
1321 > 0.425698,0.412022,0.397498,0.382650,0.366763,0.349574,
1322 > 0.330976,0.310680,0.288132,0.262614,0.233680,0.200673,
1323 > 0.163728,0.123628,0.082475,0.044455,0.014854,0.000185,
1324 > 0.627192,0.620322,0.611968,0.603362,0.593426,0.582683,
1325 > 0.571079,0.558798,0.545824,0.532347,0.519069,0.505318,
1326 > 0.491597,0.478201,0.464023,0.450314,0.436183,0.422097,
1327 > 0.407576,0.392963,0.377856,0.362363,0.346096,0.329071,
1328 > 0.310942,0.291304,0.269699,0.245618,0.218210,0.187279,
1329 > 0.152695,0.115069,0.076563,0.041003,0.013496,0.000163,
1330 > 0.631625,0.623836,0.614536,0.604467,0.593362,0.581189,
1331 > 0.568126,0.554473,0.540005,0.525342,0.510320,0.495556,
1332 > 0.480543,0.465856,0.450901,0.436047,0.421171,0.406166,
1333 > 0.391154,0.375993,0.360616,0.344868,0.328695,0.311826,
1334 > 0.294125,0.275057,0.254485,0.231491,0.205428,0.176164,
1335 > 0.143362,0.107930,0.071534,0.038049,0.012338,0.000145,
1336 > 0.628674,0.619289,0.609887,0.598892,0.587603,0.574461,
1337 > 0.560634,0.545957,0.530648,0.515019,0.499229,0.483499,
1338 > 0.467857,0.452587,0.437187,0.421817,0.406489,0.391076,
1339 > 0.375987,0.360550,0.345003,0.329399,0.313389,0.296812,
1340 > 0.279488,0.261191,0.241276,0.219209,0.194576,0.166551,
1341 > 0.135508,0.101696,0.067217,0.035540,0.011363,0.000130,
1342 > 0.619655,0.610842,0.600741,0.589506,0.577444,0.564098,
1343 > 0.550163,0.534861,0.519180,0.503427,0.486955,0.471011,
1344 > 0.455040,0.439305,0.423732,0.408132,0.392629,0.377169,
1345 > 0.361831,0.346498,0.331073,0.315724,0.299859,0.283638,
1346 > 0.266810,0.249063,0.229985,0.208716,0.185019,0.158309,
1347 > 0.128580,0.096355,0.063475,0.033376,0.010535,0.000117,
1348 > 0.607474,0.598863,0.589070,0.578071,0.565959,0.552353,
1349 > 0.538238,0.523141,0.507028,0.491211,0.474478,0.458296,
1350 > 0.442260,0.426619,0.410731,0.395250,0.379629,0.364356,
1351 > 0.348971,0.333841,0.318608,0.303265,0.287813,0.272025,
1352 > 0.255702,0.238406,0.219910,0.199549,0.176755,0.151092,
1353 > 0.122529,0.091666,0.060228,0.031496,0.009822,0.000107,
1354 > 0.593760,0.585285,0.576038,0.564716,0.553007,0.539720,
1355 > 0.525715,0.510544,0.494716,0.478059,0.461888,0.445882,
1356 > 0.429808,0.413922,0.398545,0.382937,0.367557,0.352300,
1357 > 0.337238,0.322197,0.307168,0.292182,0.277095,0.261674,
1358 > 0.245761,0.228939,0.211113,0.191442,0.169413,0.144749,
1359 > 0.117225,0.087581,0.057361,0.029844,0.009202,0.000098,
1360 > 0.579324,0.570883,0.561883,0.551297,0.539642,0.526840,
1361 > 0.512914,0.497903,0.481998,0.465838,0.449788,0.433693,
1362 > 0.417703,0.402286,0.386792,0.371437,0.356215,0.341283,
1363 > 0.326313,0.311637,0.296947,0.282140,0.267321,0.252270,
1364 > 0.236761,0.220580,0.203236,0.184230,0.162954,0.139137,
1365 > 0.112557,0.083990,0.054824,0.028391,0.008663,0.000090,
1366 > 0.564315,0.556562,0.547782,0.537905,0.526572,0.514011,
1367 > 0.500033,0.485378,0.470013,0.453915,0.438022,0.422153,
1368 > 0.406530,0.391112,0.375843,0.360792,0.345939,0.331089,
1369 > 0.316383,0.301913,0.287442,0.273075,0.258564,0.243883,
1370 > 0.228827,0.213013,0.196181,0.177712,0.157167,0.134036,
1371 > 0.108376,0.080670,0.052555,0.027085,0.008183,0.000084,
1372 > 0.549474,0.542540,0.534001,0.524222,0.513538,0.501429,
1373 > 0.487998,0.473665,0.458530,0.442760,0.426794,0.411187,
1374 > 0.395839,0.380807,0.365652,0.350917,0.336154,0.321565,
1375 > 0.307328,0.293072,0.278886,0.264801,0.250668,0.236339,
1376 > 0.221517,0.206147,0.189835,0.171902,0.151860,0.129450,
1377 > 0.104545,0.077707,0.050505,0.025913,0.007753,0.000078,
1378 > 0.535066,0.528539,0.520546,0.511350,0.500893,0.489256,
1379 > 0.476407,0.462092,0.447135,0.431648,0.416302,0.400932,
1380 > 0.385885,0.370871,0.356128,0.341562,0.327354,0.312933,
1381 > 0.298890,0.284864,0.270948,0.257177,0.243326,0.229321,
1382 > 0.214967,0.199986,0.183994,0.166574,0.147078,0.125313,
1383 > 0.101089,0.075065,0.048640,0.024855,0.007368,0.000073,
1384 > 0.509409,0.510722,0.511674,0.511738,0.511286,0.510766,
1385 > 0.509198,0.507564,0.505152,0.502704,0.499650,0.496452,
1386 > 0.492756,0.488747,0.483986,0.478817,0.472835,0.466164,
1387 > 0.458412,0.449362,0.438892,0.426796,0.412430,0.395391,
1388 > 0.375497,0.352110,0.324651,0.292275,0.254919,0.212461,
1389 > 0.165886,0.117521,0.071431,0.033526,0.009193,0.000091,
1390 > 0.570888,0.568118,0.564921,0.561074,0.555970,0.550478,
1391 > 0.544438,0.537826,0.530394,0.523005,0.514685,0.506342,
1392 > 0.497778,0.489102,0.479852,0.470348,0.460202,0.449558,
1393 > 0.438121,0.426046,0.412891,0.398444,0.382493,0.364765,
1394 > 0.344689,0.321925,0.295886,0.265915,0.231610,0.192817,
1395 > 0.150406,0.106352,0.064426,0.030028,0.008099,0.000077,
1396 > 0.605452,0.599834,0.593296,0.586627,0.578039,0.569200,
1397 > 0.559754,0.549624,0.538672,0.527535,0.516054,0.504396,
1398 > 0.492657,0.481104,0.468974,0.456682,0.443949,0.431100,
1399 > 0.417683,0.403750,0.389161,0.373852,0.357151,0.339250,
1400 > 0.319633,0.297684,0.273111,0.245042,0.213167,0.177272,
1401 > 0.138115,0.097498,0.058870,0.027271,0.007249,0.000067,
1402 > 0.621534,0.614157,0.605835,0.596272,0.585946,0.574335,
1403 > 0.562822,0.550233,0.536841,0.523598,0.509920,0.496293,
1404 > 0.482476,0.468781,0.455150,0.441275,0.427053,0.412946,
1405 > 0.398426,0.383702,0.368436,0.352438,0.335827,0.317985,
1406 > 0.298878,0.277758,0.254411,0.228040,0.198146,0.164659,
1407 > 0.128112,0.090261,0.054331,0.025017,0.006563,0.000059,
1408 > 0.625646,0.616575,0.607120,0.596395,0.584780,0.572447,
1409 > 0.558731,0.544391,0.529948,0.515041,0.499832,0.485344,
1410 > 0.470145,0.455555,0.440783,0.425851,0.410922,0.396061,
1411 > 0.381038,0.365825,0.350234,0.334195,0.317648,0.300202,
1412 > 0.281570,0.261329,0.238985,0.213898,0.185675,0.154079,
1413 > 0.119726,0.084159,0.050488,0.023108,0.005981,0.000052,
1414 > 0.621816,0.611958,0.601503,0.590515,0.577713,0.564667,
1415 > 0.550311,0.535181,0.519529,0.503933,0.488235,0.472641,
1416 > 0.456837,0.441572,0.426210,0.410935,0.395680,0.380387,
1417 > 0.365098,0.349710,0.334114,0.318243,0.301798,0.284744,
1418 > 0.266642,0.247149,0.225778,0.201848,0.175065,0.145114,
1419 > 0.112575,0.078976,0.047163,0.021485,0.005492,0.000047,
1420 > 0.612456,0.602615,0.592112,0.580629,0.567547,0.554000,
1421 > 0.539266,0.523505,0.507494,0.491317,0.475170,0.459335,
1422 > 0.443356,0.427588,0.412089,0.396613,0.381018,0.365836,
1423 > 0.350418,0.335168,0.319712,0.303995,0.287942,0.271269,
1424 > 0.253737,0.234846,0.214374,0.191514,0.165911,0.137420,
1425 > 0.106436,0.074522,0.044401,0.020096,0.005079,0.000042,
1426 > 0.599906,0.590658,0.579871,0.568198,0.555552,0.541752,
1427 > 0.526573,0.511088,0.494690,0.478321,0.462067,0.446064,
1428 > 0.430078,0.414174,0.398700,0.383180,0.367731,0.352392,
1429 > 0.337125,0.322012,0.306729,0.291406,0.275582,0.259369,
1430 > 0.242384,0.224158,0.204420,0.182500,0.157984,0.130669,
1431 > 0.101087,0.070603,0.041960,0.018893,0.004724,0.000039,
1432 > 0.586385,0.576832,0.566175,0.554843,0.542371,0.528570,
1433 > 0.513481,0.497995,0.481584,0.465454,0.449166,0.432955,
1434 > 0.417165,0.401569,0.386012,0.370517,0.355274,0.340119,
1435 > 0.324992,0.310107,0.295066,0.280058,0.264679,0.248876,
1436 > 0.232336,0.214712,0.195681,0.174542,0.150960,0.124776,
1437 > 0.096383,0.067222,0.039821,0.017842,0.004416,0.000035,
1438 > 0.571562,0.562230,0.553483,0.541106,0.528820,0.515154,
1439 > 0.500724,0.484959,0.468895,0.452796,0.436478,0.420645,
1440 > 0.404877,0.389321,0.373971,0.358762,0.343687,0.328774,
1441 > 0.313968,0.299216,0.284619,0.269828,0.254780,0.239398,
1442 > 0.223345,0.206321,0.187886,0.167537,0.144752,0.119523,
1443 > 0.092229,0.064212,0.037929,0.016920,0.004149,0.000033,
1444 > 0.556430,0.547716,0.538118,0.527507,0.515330,0.502167,
1445 > 0.487523,0.472378,0.456637,0.440513,0.424536,0.408900,
1446 > 0.393308,0.378059,0.362874,0.347899,0.333001,0.318366,
1447 > 0.303866,0.289409,0.275001,0.260570,0.245889,0.230899,
1448 > 0.215327,0.198784,0.180938,0.161191,0.139240,0.114867,
1449 > 0.088515,0.061527,0.036242,0.016094,0.003912,0.000030,
1450 > 0.541088,0.533114,0.524472,0.513843,0.502197,0.489243,
1451 > 0.475392,0.460194,0.444773,0.428968,0.413337,0.397800,
1452 > 0.382488,0.367387,0.352568,0.337799,0.323189,0.308771,
1453 > 0.294545,0.280383,0.266247,0.252107,0.237863,0.223231,
1454 > 0.208050,0.191998,0.174706,0.155516,0.134237,0.110630,
1455 > 0.085174,0.059101,0.034727,0.015356,0.003700,0.000028,
1456 > 0.526746,0.519306,0.510794,0.500680,0.489449,0.477026,
1457 > 0.463177,0.448622,0.433350,0.417867,0.402498,0.387171,
1458 > 0.372235,0.357466,0.342866,0.328463,0.314082,0.299986,
1459 > 0.285991,0.272106,0.258305,0.244422,0.230511,0.216246,
1460 > 0.201466,0.185822,0.169011,0.150402,0.129717,0.106819,
1461 > 0.082146,0.056897,0.033348,0.014690,0.003511,0.000026,
1462 > 0.508294,0.509537,0.510232,0.510364,0.510028,0.509066,
1463 > 0.507574,0.506096,0.503764,0.501269,0.498261,0.494940,
1464 > 0.491357,0.487139,0.482585,0.477518,0.471487,0.464685,
1465 > 0.456954,0.448138,0.437547,0.425359,0.410953,0.393893,
1466 > 0.373877,0.350282,0.322653,0.290145,0.252556,0.209889,
1467 > 0.163221,0.114885,0.069163,0.031966,0.008546,0.000081,
1468 > 0.569358,0.566803,0.562983,0.559001,0.554723,0.548955,
1469 > 0.542742,0.536150,0.528626,0.520981,0.513074,0.504898,
1470 > 0.496364,0.487653,0.478428,0.468815,0.458593,0.447996,
1471 > 0.436563,0.424615,0.411428,0.396895,0.381014,0.363131,
1472 > 0.342953,0.320173,0.294051,0.263852,0.229272,0.190367,
1473 > 0.147875,0.103902,0.062322,0.028615,0.007526,0.000069,
1474 > 0.603670,0.598438,0.591834,0.584727,0.576222,0.567680,
1475 > 0.558196,0.547591,0.536787,0.525654,0.514148,0.502844,
1476 > 0.491102,0.479283,0.467390,0.455046,0.442422,0.429358,
1477 > 0.416062,0.402212,0.387705,0.372222,0.355668,0.337678,
1478 > 0.318005,0.295922,0.271145,0.242997,0.210936,0.174968,
1479 > 0.135724,0.095201,0.057278,0.025971,0.006733,0.000060,
1480 > 0.619938,0.612660,0.603237,0.594399,0.584220,0.572945,
1481 > 0.560906,0.547974,0.535000,0.521640,0.508074,0.494503,
1482 > 0.480788,0.467052,0.453309,0.439430,0.425326,0.411154,
1483 > 0.396834,0.381957,0.366731,0.350903,0.334155,0.316435,
1484 > 0.297197,0.275976,0.252539,0.225982,0.196020,0.162415,
1485 > 0.125838,0.088092,0.052512,0.023816,0.006094,0.000053,
1486 > 0.623752,0.615162,0.604853,0.594199,0.582617,0.570319,
1487 > 0.556672,0.542512,0.527993,0.513108,0.498185,0.483234,
1488 > 0.468378,0.453629,0.438899,0.424139,0.409158,0.394331,
1489 > 0.379329,0.364110,0.348592,0.332662,0.316045,0.298521,
1490 > 0.279860,0.259548,0.237136,0.211942,0.183600,0.151947,
1491 > 0.117534,0.082104,0.048767,0.021982,0.005551,0.000047,
1492 > 0.619727,0.610281,0.600069,0.588100,0.575830,0.562394,
1493 > 0.548038,0.533038,0.517667,0.501730,0.486231,0.470564,
1494 > 0.454906,0.439606,0.424254,0.409150,0.393803,0.378587,
1495 > 0.363406,0.347957,0.332413,0.316594,0.300201,0.283135,
1496 > 0.264987,0.245411,0.223977,0.199987,0.173072,0.143020,
1497 > 0.110486,0.077018,0.045586,0.020431,0.005096,0.000042,
1498 > 0.610495,0.600588,0.589924,0.578181,0.565270,0.551795,
1499 > 0.537196,0.521624,0.505437,0.489150,0.472983,0.457321,
1500 > 0.441422,0.425796,0.410213,0.394710,0.379352,0.363971,
1501 > 0.348802,0.333422,0.317927,0.302300,0.286286,0.269619,
1502 > 0.252021,0.233168,0.212605,0.189631,0.163957,0.135368,
1503 > 0.104416,0.072632,0.042854,0.019101,0.004711,0.000038,
1504 > 0.597856,0.588522,0.577800,0.566289,0.553520,0.539436,
1505 > 0.524435,0.508812,0.492453,0.476231,0.459950,0.443865,
1506 > 0.427942,0.412263,0.396763,0.381179,0.365851,0.350734,
1507 > 0.335459,0.320164,0.305048,0.289712,0.274016,0.257736,
1508 > 0.240735,0.222511,0.202688,0.180671,0.156046,0.128701,
1509 > 0.099140,0.068812,0.040482,0.017948,0.004381,0.000035,
1510 > 0.584084,0.574752,0.564446,0.552705,0.540225,0.526223,
1511 > 0.511434,0.495863,0.479694,0.463424,0.447226,0.431058,
1512 > 0.415167,0.399560,0.384026,0.368648,0.353514,0.338375,
1513 > 0.323406,0.308441,0.293400,0.278400,0.263020,0.247257,
1514 > 0.230713,0.213051,0.193943,0.172756,0.149112,0.122832,
1515 > 0.094502,0.065475,0.038404,0.016946,0.004094,0.000032,
1516 > 0.569517,0.560385,0.550236,0.539113,0.526421,0.513222,
1517 > 0.498378,0.482790,0.466696,0.450744,0.434441,0.418603,
1518 > 0.402901,0.387328,0.372129,0.356950,0.341855,0.326972,
1519 > 0.312231,0.297539,0.282901,0.268148,0.253147,0.237801,
1520 > 0.221735,0.204645,0.186195,0.165769,0.142932,0.117672,
1521 > 0.090399,0.062530,0.036573,0.016065,0.003846,0.000029,
1522 > 0.554361,0.545724,0.536138,0.525127,0.513348,0.500011,
1523 > 0.485538,0.470383,0.454424,0.438520,0.422677,0.406888,
1524 > 0.391373,0.376154,0.360890,0.345988,0.331296,0.316666,
1525 > 0.302083,0.287727,0.273311,0.258904,0.244293,0.229265,
1526 > 0.213676,0.197142,0.179286,0.159481,0.137449,0.113024,
1527 > 0.086731,0.059884,0.034929,0.015279,0.003625,0.000027,
1528 > 0.539566,0.531375,0.522021,0.511690,0.499982,0.487119,
1529 > 0.474212,0.458150,0.442576,0.426822,0.411177,0.395600,
1530 > 0.380520,0.365481,0.350710,0.336033,0.321423,0.307074,
1531 > 0.292885,0.278746,0.264660,0.250555,0.236263,0.221649,
1532 > 0.206449,0.190350,0.173011,0.154255,0.132811,0.108865,
1533 > 0.083444,0.057506,0.033461,0.014573,0.003428,0.000025,
1534 > 0.524865,0.517319,0.508721,0.498692,0.487348,0.474698,
1535 > 0.461245,0.446504,0.431299,0.416072,0.400507,0.385278,
1536 > 0.370252,0.355605,0.341042,0.326661,0.312354,0.298301,
1537 > 0.284289,0.270473,0.256711,0.242903,0.228985,0.214684,
1538 > 0.199885,0.184233,0.167367,0.148740,0.128003,0.105093,
1539 > 0.080442,0.055358,0.032128,0.013936,0.003252,0.000024/
1540
1541C IF (X.GE.XMAXX.OR.Q2.LT.Q2MIN.OR.Q2.GT.Q2MAX) THEN
1542C WRITE(*,*) 'X or Q2 out of range!!'
1543C WRITE(*,*) 'Valid range is:'
1544C WRITE(*,*) '1E-05 < X < 0.95 and'
1545C WRITE(*,*) '4 < Q2 < 520 GeV**2'
1546C RETURN
1547C END IF
1548
1549 IF (X.GE.XMAX.OR.B.GT.BMAX) THEN
1550 RES(1) = 1.0
1551 RES(2) = 1.0
1552 RETURN
1553 END IF
1554 IF(Q2.LT.Q2MIN) THEN
1555 Q2 = Q2MIN
1556 ELSE IF(Q2.GT.Q2MAX) THEN
1557 Q2 = Q2MAX
1558 ENDIF
1559
1560 CALL GGINTER(INUCL,X,Q2,SHAD)
1561
1562 DO I=1,31
1563 TATMP(I) = TA(I,INUCL)
1564 ENDDO
1565 CALL SPLINE(TMAX,IMPAR,TATMP,C,D,E)
1566 TAF = SEVAL(TMAX,B,IMPAR,TATMP,C,D,E)
1567
1568 DO IK=1,2
1569 RES(IK) = 1.0/(1.0 + (ANUCL(INUCL)-1.0)*SHAD(IK)*TAF)
1570 ENDDO
1571
1572 RETURN
1573 END
1574
1575 SUBROUTINE GGINTER(INUCL,X,Q2,SHAD)
1576 IMPLICIT DOUBLE PRECISION(A-H,L-Z)
1577 IMPLICIT INTEGER(I-K)
1578
1579 DIMENSION SHAD(2)
1580 DIMENSION G(36,13,4),LQ(36,13,4)
1581 DIMENSION XB(36),Q2V(13)
1582 DIMENSION GQTMP(13),LQQTMP(13)
1583 DIMENSION GXTMP(36),LQXTMP(36)
1584 DIMENSION C(100),D(100),E(100)
1585
1586 PARAMETER(INMAX=13,IMMAX=36)
1587
1588 COMMON/GG07/ XB,Q2V,G,LQ
1589
1590 DO K=1,36
1591 DO J=1,13
1592 GQTMP(J) = G(K,J,INUCL)
1593 LQQTMP(J) = LQ(K,J,INUCL)
1594 ENDDO
1595 CALL SPLINE(INMAX,Q2V,GQTMP,C,D,E)
1596 GXTMP(K) = SEVAL(INMAX,Q2,Q2V,GQTMP,C,D,E)
1597 CALL SPLINE(INMAX,Q2V,LQQTMP,C,D,E)
1598 LQXTMP(K) = SEVAL(INMAX,Q2,Q2V,LQQTMP,C,D,E)
1599 ENDDO
1600
1601 CALL SPLINE(IMMAX,XB,GXTMP,C,D,E)
1602 SHAD(1) = SEVAL(IMMAX,X,XB,GXTMP,C,D,E)
1603 CALL SPLINE(IMMAX,XB,LQXTMP,C,D,E)
1604 SHAD(2) = SEVAL(IMMAX,X,XB,LQXTMP,C,D,E)
1605
1606 RETURN
1607 END
1608
1609C ---------------------------------------------------------------------
1610 SUBROUTINE SPLINE(N,X,Y,B,C,D)
1611C ---------------------------------------------------------------------
1612c***************************************************************************
1613C CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
1614C INTERPOLATION SUBROUTINES ARE TAKEN FROM
1615C G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
1616C COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
1617 IMPLICIT double precision (A-H,O-Z)
1618 DIMENSION X(100),Y(100),B(100),C(100),D(100)
1619 NM1=N-1
1620 IF(N.LT.2) RETURN
1621 IF(N.LT.3) GO TO 250
1622 D(1)=X(2)-X(1)
1623 C(2)=(Y(2)-Y(1))/D(1)
1624 DO 210 I=2,NM1
1625 D(I)=X(I+1)-X(I)
1626 B(I)=2.0D0*(D(I-1)+D(I))
1627 C(I+1)=(Y(I+1)-Y(I))/D(I)
1628 C(I)=C(I+1)-C(I)
1629 210 CONTINUE
1630 B(1)=-D(1)
1631 B(N)=-D(N-1)
1632 C(1)=0.0D0
1633 C(N)=0.0D0
1634 IF(N.EQ.3) GO TO 215
1635 C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
1636 C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
1637 C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
1638 C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
1639 215 CONTINUE
1640 DO 220 I=2,N
1641 T=D(I-1)/B(I-1)
1642 B(I)=B(I)-T*D(I-1)
1643 C(I)=C(I)-T*C(I-1)
1644 220 CONTINUE
1645 C(N)=C(N)/B(N)
1646 DO 230 IB=1,NM1
1647 I=N-IB
1648 C(I)=(C(I)-D(I)*C(I+1))/B(I)
1649 230 CONTINUE
1650 B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0D0*C(N))
1651 DO 240 I=1,NM1
1652 B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0D0*C(I))
1653 D(I)=(C(I+1)-C(I))/D(I)
1654 C(I)=3.0D0*C(I)
1655 240 CONTINUE
1656 C(N)=3.0D0*C(N)
1657 D(N)=D(N-1)
1658 RETURN
1659 250 CONTINUE
1660 B(1)=(Y(2)-Y(1))/(X(2)-X(1))
1661 C(1)=0.0D0
1662 D(1)=0.0D0
1663 B(2)=B(1)
1664 C(2)=0.0D0
1665 D(2)=0.0D0
1666 RETURN
1667 END
1668c
1669c***************************************************************************
1670C ---------------------------------------------------------------------
1671 double precision FUNCTION SEVAL(N,XX,X,Y,B,C,D)
1672C ---------------------------------------------------------------------
1673c***************************************************************************
1674C CALCULATE THE DISTRIBUTION AT XX BY CUBIC SPLINE INTERPOLATION.
1675 implicit double precision(A-H,O-Z)
1676 DIMENSION X(100),Y(100),B(100),C(100),D(100)
1677 DATA I/1/
1678 IF(I.GE.N) I=1
1679 IF(XX.LT.X(I)) GO TO 310
1680 IF(XX.LE.X(I+1)) GO TO 330
1681 310 CONTINUE
1682 I=1
1683 J=N+1
1684 320 CONTINUE
1685 K=(I+J)/2
1686 IF(XX.LT.X(K)) J=K
1687 IF(XX.GE.X(K)) I=K
1688 IF(J.GT.I+1) GO TO 320
1689 330 CONTINUE
1690 DX=XX-X(I)
1691 SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))
1692 RETURN
1693 END
1694c
1695c***************************************************************************
1696
1697
1698* function to generate gauss distribution
1699 double precision function gaus(x0,sig)
1700 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1701 41 u1=pyr(0)
1702 u2=pyr(0)
1703 v1=2.d0*u1-1.d0
1704 v2=2.d0*u2-1.d0
1705 s=v1**2+v2**2
1706 if(s.gt.1) go to 41
1707 gaus=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
1708 return
1709 end