1 #include "BlockData/cblkGene.h" 7 real*8 xin(maxbin), yin(maxbin)
9 integer nbin0, icon, code, tcode
10 integer count, status, cond
11 integer,
external:: NARGS
14 if(count .eq. 3) cond=1
15 if(count .eq. 2) cond=0
16 if( count .ne. 2 .and. count .ne. 3 )
then 17 write(0,*)
'Usage: latFit.. code [c] < inputfile > outfile' 18 write(0,*)
' code = particle code,',
19 *
' inputfile is output from splitLat.csh ' 20 write(0,*)
' c: if exists, write only coeff. on stdout ' 21 write(0,*)
' if omitted, coeff. is output on stderr,',
22 *
' fitted (x,y) on stdout' 25 call getarg(1, buf, status)
38 read(*,*,end=100) xin(nbin0+1), yin(nbin0+1)
42 if(nbin0 .gt. maxbin)
then 43 write(0,*)
' too many lat data> ', maxbin
48 call copenfw2(minout,
"/dev/null", 1, icon)
52 call fitlat0(cond, tcode, code, xin, yin, nbin0)
55 subroutine fitlat0(cond, tcode, code, xin, yin, nbin)
59 integer nbin, code, tcode
60 real*8 xin(nbin), yin(nbin)
61 real*8 xuse(nbin), yuse(nbin)
62 real*8 prmout(nparam, 4, nregion)
65 integer n1, n2, region
79 if(xin(i) .gt.
x2(region) )
exit 80 if(xin(i) .lt.
x1(region) ) cycle
92 call fitlat1(region, code, xuse, yuse, n1,
93 *
param(1, code, region), prmout(1, code,region))
97 write(*,
'(5g12.4)') prmout(1, code,region),
98 * prmout(2,code, region),
99 * prmout(3, code, region), prmout(4,code, region),
108 write(0,
'(5g12.4)') prmout(1,code, region),
109 * prmout(2, code, region),
110 * prmout(3, code, region), prmout(4, code,region),
116 pw = prmout(4,code, region)
117 do while ( xx .le.
drx2(region) )
118 f=prmout(1,code,region)/
119 * xx**(prmout(2,code, region) +
121 * prmout(3,code, region)* xx**
pw )
129 subroutine fitlat1(region, code, xin, yin, n, prmin, prmout )
133 real*8 prmin(nparam), prmout(nparam)
135 real*8 xin(n), yin(n)
137 integer nlabel(nparam)
138 character*10 pname(nparam)
139 real*8 initval(nparam)
142 data nlabel/ 1, 2, 3, 4/
143 data pname/
'p',
'q',
'r',
'pw'/
144 data step/ 1., 0.001
d0, 0.0001
d0, 0.1
d0/
145 real*8 zero, one, three, four, five
146 data zero,one,three,four, five / 0., 1., 3.,4., 5. /
165 initval(i) = prmin(i)
168 call mninit( 5, minout, minsave)
178 call mnparm(nlabel(i), pname(i), initval(i), step(i),
179 *
low(i, code,region),
up(i, code,region), ierflg)
181 if (ierflg .ne. 0)
then 182 write (0,
'(a,i3)')
' unable to define parameter no.',i
187 call mnseti(
'lat as a function of core distance')
209 if(region .eq. 4 )
then 212 if(maxdiff .gt. 0.01)
then 213 badindex(maxindex)= -maxindex
218 if(maxdiff .gt. 0.033)
then 221 badindex(maxindex)= -maxindex
block data include Zlatfit h c fitting region data x1(1)/0.03/
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data h g *is for param c g data up(2, 1)/7.0d0/
subroutine mnexcm(FCN, COMAND, PLIST, LLIST, IERFLG, FUTIL)
subroutine fitlat1(region, code, xin, yin, n, prmin, prmout)
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
subroutine copenfw2(io, fnin, form, icon)
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer op real x1h common Zfitc * param
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer op real x1h common Zfitc drx1
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data h g *is for param c g data low(1, 1)/1.d-5/
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer npoint
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
subroutine mnparm(K, CNAMJ, UK, WK, A, B, IERFLG)
block data include Zlatfit h c fitting region data x2(1)/0.5/data x1(2)/0.3/
subroutine mninit(I1, I2, I3)
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer op real x1h common Zfitc * oparam
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer op real x1h common Zfitc drx2
subroutine latfnc(npar, gin, f, paramx, iflag)
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
subroutine fitlat0(cond, tcode, code, xin, yin, nbin)