#!/bin/sh # modify input - via rcp - see tmzc_old_ifk for older input 6-3-07 # or use sund.amtp.liv.ac.uk:/export/home/sund in place of /ukqcd6 # 27-4-06 bugfix to F,L order # 28-4-06 bugfix - choice of observables # based tmz - now 20 operators - up to 6x6 fits # set -x # uncomment for more in-line diagnostics # as is needs fitex.f as well as this "here document - script" # needs -l nag18 also # and input files in ${DATAPATH}.??.???? # full output in out.fitz # summary in summary.fitz (if IFBOOT=1) # tmpr1.dat is axis fit file # fort.30 fort.31 has some PCAC info - for plots ## fort.69 is pion average correlators, 70 is rho USERSDIR=`pwd` ## parameters to set #---------- choose gauge set-up labelled ifk 0,1,2,.. ## 0 INFO="24^3 48 twist 3.9 .1608 mu=0.004 nlong=4" ## 1 INFO=" 32^3 64 twist 4.05 k=.156985 mu=0.003 every 5 traj nlong=6" # 4 INFO=" 24^3 48 twist 3.9 .160856 mu=0.004 every 5 (*2 traj) nlong=6" # 41 INFO="24^3 48 3.9 dd 004 k=1 LL LF nlong=6" # 432 INFO="32^3 64 3.9 mu=.004 random t_s every 4(*2 traj) nlong=6" # 3 INFO=" 32^3 64 twist 4.05 k=0.157025 mu=0.003 every 5 traj nlong=6" # 323 INFO=" 32^3 64 twist 4.05 k=0.15701 mu=0.003 every 5 traj nlong=7" # 3231 INFO=" 32^3 64 k=1 4.05 k=0.15701 mu=0.003 every 5 traj nlong=7" # 326 INFO=" 32^3 64 twist 4.05 k=0.15701 mu=0.006 every 5 traj nlong=7" # 246 INFO=" 24^3 48 twist 4.05 k=0.15701 mu=0.006 every 5 traj nlong=7" # 206 INFO=" 20^3 48 twist 4.05 k=0.15701 mu=0.006 every 40 traj nlong=6" # 328 INFO=" 32^3 64 twist 4.05 k=0.15701 mu=0.008 every 5 traj nlong=7" # 3212 INFO="32^3 64 twist 4.05 k=0.15701 mu=0.012 every 5 traj nlong=7" # 85 INFO=" 24^3 48 twist 3.75 k=0.1660 mu=0.0085 every 5 traj nlong=5" # 851 INFO=" 24^3 48 twist 3.75 k=0.1660 mu=0.0085 every 5 traj nlong=5" # 852 INFO=" 24^3 48 twist 3.75 k=0.1660 mu=0.0085 every 5 traj nlong=5" # 10 INFO=" 24^3 48 twist 3.9 .160856 mu=0.01 every 5 (*2 traj) nlong=6" # 8 INFO="24^3 48 twist 3.9 .160856 mu=0.0085 every 5 (*2 traj) nlong=6" # 81 INFO="24^3 48 twist 3.9 .160856 mu=0.0085 every 5 (*2 traj) nlong=6" # 15 INFO=" 24^3 48 twist 3.9 .160856 mu=0.015 every 5 (*2 traj) nlong=6" # 64 INFO=" 24^3 48 twist 3.9 .160856 mu=0.0064 every 8 (*2 traj) nlong=6" # 99 INFO=" TEST " # 6''' INFO="20^3 48 twist 3.8 .16428 mu=0.006 every 5 (*2 traj)nlong=6" # 6' INFO="20^3 48 twist 3.8 .164084 mu=0.006 every 5 (*2 traj) nlong=6" # 60 INFO="20^3 48 twist 3.8 .16415 mu=0.006 every 5 (*2 traj) nlong=6" # 61 INFO="20^3 48 twist 3.8 .164099 mu=0.006 every 5 (*2 traj) nlong=6" # 6 INFO="20^3 48 twist 3.8 .164111 mu=0.006 every 5 (*2 traj) nlong=6" # 45 INFO="24^3 48 twist 3.8 .164111 mu=0.0045 every 5 (*2 traj)nlong=6" # 624 INFO="24^3 48 twist 3.8 .164111 mu=0.006 every 2 (*2 traj)nlong=6" # 11 INFO="24^3 48 twist 3.8 .164111 mu=0.0045 every 5 (*2 traj)nlong=6" # 9 INFO="20^3 48 twist 3.8 .164099 mu=0.009 every 5 (*2 traj) nlong=6" # 12 INFO="20^3 48 twist 3.8 .164099 mu=0.012 every 5 (*2 traj)nlong=6" # 15 INFO="20^3 48 twist 3.8 .164099 mu=0.012 every 5 (*2 traj) nlong=6" # 26 INFO="24^3 24 twist 3.9 sea 004 valence 026 nlong=6" # 6661 INFO="stout test" # 666 INFO=" stout test dec 2010 mu=0.006 K=13563" # 444 INFO="2+1+1 test" run1 and run2 # 446 INFO="2+1+1 test" # 448 2+1+1 L24T48_b1.90_k0.16326_mu0.008_musigma0.15_mudelta0.19 # 456 2+1+1 STOUT test 004 # 66 INFO="stout test" # 2112 2+1+1 2.1 48^3 002 # 455 2+1+1 1.95 32^3 # 435 2+1+1 1.95 32^3 0.035 # 3.75 STOUT 77 66 3.8 1476 006 776 1475 775 # 3.75 85 851(BGL) 852(QCDOC) # 3.8 20^3 6 (60, 61)(several k_v) 9 12 153 38(test) # 24^3 45 624 824 11 165 # 3.9 4 64 8 10 15(2 series) 99=TEST # 41 mom=1 81 mom=1 [L only choose FCASE=1 ] # 332 432 (32^3) # 27 charm # 4.05 323(rd) 3232(2 series-old ) 326 328 3212 # 3231 mom=1 # 246 206 # 4.2 482 IFK=435 echo " ifk="$IFK NOTAR=0 # NOTAR=1 # usable 4 64 10 15O 38 6 666 # NOTAR=1 assumes /tmp/cmi has out files already - does not reread #---- choose particle type--------------------- # iobval=1 is PION 2 is rho-a1 3 is a0 4 is b1 IOBVAL=1 if [ $IOBVAL -eq 1 ] then PARTICLE="pion" NPV=21 NFACP=6 fi if [ $IOBVAL -eq 2 ] then PARTICLE="rho-a1" NPV=21 NFACP=6 fi if [ $IOBVAL -eq 3 ] then PARTICLE="a0" NPV=3 NFACP=2 fi if [ $IOBVAL -eq 4 ] then PARTICLE="b1" NPV=3 NFACP=2 fi #----- t range to fit---------------------------------- # default: for rho etc (rho 4-12 5-15 5-16 6-18 2 exp IT1=4 IT1=6 #IT1=7 #IT1=8 IT1=10 #IT2=8 #IT2=12 #IT2=15 ## use 18 3.9 20 4.05 IT2=18 #IT2=20 # a0 b1 (b1 4-8 2 exp 3-8 4.05: 6-12 5-14 1 exp 4-12 2 exp #IT1=3 #IT1=4 IT1=5 #IT1=6 #IT2=8 #IT2=10 IT2=12 IT2=14 ## for pion (2 exp pi fits need careful intialisation) # 3.9 7-23 8-18; 4.05 8-31 10-23 if [ $IOBVAL -eq 1 ] then #IT1=5 #IT1=6 #IT1=7 #IT1=9 IT1=10 #IT1=12 #### IT1=3 #IT1=12 #IT2=15 #IT2=18 IT2=23 #IT2=20 IT2=31 #IT2=47 ### it2 gets set to LX4-1 if above this # to mimic pi_0 #IT1=7 #IT2=15 fi ### ------------------- fit decisions---------------------------- # choose paths to fit (and factorising coefficients): FCASE # default FCASE=0 # order used in fits: # ip1 /1,1,2,1,2,3,1,2,3,4,1,2,3,4,5,1,2,3,4,5,6 / # ip2 /1,2,2,3,3,3,4,4,4,4,5,5,5,5,5,6,6,6,6,6,6 / # default full fit ($NPV first values from list with $NFACP coefficients): IPTH='/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21/' # IPMAP maps chosen operators to 1,2,3,4,.. ie ipmap(ip1(ipth(ip))) IPMAP='/1,2,3,4,5,6,15*0/' # pi in 12 56 and tuned to zero in 34 pi: g5 g5g4 g4 # mpcac use 1234 ZV use 1256 # rho in 12 56 (and a1 in 34) rho gig4 gi gig5 FCASE=0 #FCASE=1234 #FCASE=12 # for IFK=41 81 3231 27 (no FF) so use L only # (care LF may be norm wrongly) Care J=1 # or for LL LF use 19 59 #FCASE=1 #FCASE=19 #FCASE=59 #FCASE=13 #FCASE=1256 #FCASE=56 #FCASE=34 #FCASE=2 #FCASE=26 case $FCASE in 1) ## first one only NPV=1 NFACP=1;; 12) ## first two only (iobval = 1 or 2) NPV=3 NFACP=2;; 2) ## F (no 2) only (iobval = 1 or 2) NPV=1 NFACP=1 IPTH='/3,2,1,11,12,15,16,17,20,21,11*1/' IPMAP='/2,1,3,4,5,6,15*0/';; 1234) ## first 4 only (iobval = 1 or 2) NPV=10 NFACP=4;; 1256) ## first 2 and last 2 (iobval = 1 or 2) eg pi I g4 or rho gig4 gig5 NPV=10 NFACP=4 IPTH='/1,2,3,11,12,15,16,17,20,21,11*1/' IPMAP='/1,2,5,6,3,4,15*0/';; 26) ## F.F ops 1 & 3 NPV=3 NFACP=2 IPTH='/3,17,21,18*1/' IPMAP='/3,1,5,6,4,2,15*0/' ;; 56) ## last 2 only (5 and 6)(iobval = 1 or 2) eg rho gig5 NPV=3 NFACP=2 IPTH='/15,20,21,18*1/' IPMAP='/3,4,5,6,1,2,15*0/' ;; 6) ## last 1 only (6)(iobval = 1 or 2) eg rho gig5 NPV=1 NFACP=1 IPTH='/21,20*1/' IPMAP='/3,4,5,6,2,1,15*0/' ;; 34) ## middle 2 only (3 and 4)(iobval = 1 or 2) eg a1 NPV=3 NFACP=2 IPTH='/6,9,10,11,12,15,16,17,20,21,11*1/' IPMAP='/3,4,1,2,5,6,15*0/';; 13) ## L only (iobval=1 or 2) NPV=3 NFACP=2 IPTH='/1,4,6,11,13,15,15*1/' IPMAP='/1,4,2,5,3,6,15*0/';; 19) ## LL LF only (care may average in 0 for FL) NPV=2 NFACP=2 IPTH='/1,2,6,11,13,15,15*1/' IPMAP='/1,2,3,4,5,6,15*0/';; 59) ## LL LF only (iobval=1 or 2 last spin, care may average in 0 for FL) NPV=2 NFACP=2 IPTH='/15,20,6,11,13,1,15*1/' IPMAP='/5,6,3,4,1,2,15*0/';; 135) ## L only (iobval=1 or 2) NPV=6 NFACP=3 IPTH='/1,4,6,11,13,15,15*1/' IPMAP='/1,4,2,5,3,6,15*0/';; 246) ## F only (iobval=1 or 2) NPV=6 NFACP=3 IPTH='/3,8,10,17,19,21,15*1/' IPMAP='/4,1,5,2,6,3,15*0/';; esac #--- choose if full jacknife/bootstrap errors (1; slower) or not (0)---- # 11 set boot list and write x; 99 use bootlist and write x IFBOOT=1 #---choose correlation (in t and path) model in fit------------ # see prologue to subroutine corfit (in fitex.f) for further information ## uncorrelated CORRELATION='icorrt=0,icorrp=0,idelt=4,ievals=10,xlmin=0.001' ## correlated # #CORRELATION='icorrt=3,icorrp=2,idelt=4,ievals=6,xlmin=0.001' #-- choose number of states to fit (typically 1-3)--------------- # pion NEXP=1 IT1=6 IT2=23 is OK - NFACP=6 # pion NEXP=2 IT1=4 IT2=23 is OK - NFACP=2 # rho NEXP=2 IT1=4 IT2=15 is OK - NFACP=2 # rho NEXP=2 IT1=5 IT2=15 is OK - NFACP=6 + KLudge # rho NEXP=3 IT1=4 IT2=15 is OK - NFACP=6 + Kludge NEXP=1 if [ $IOBVAL -eq 4 ] then NEXP=2 fi # Kludge - to make a1 coefficient in gig4 zero #### this is not exactly correct - except in continuum limit A1COEFF0=0 if [ $IOBVAL -eq 2 -a $NFACP -eq 6 -a $NEXP -gt 1111 ] then # kills 7 and 8 A1COEFF0=123 fi if [ $IOBVAL -eq 2 -a $FCASE -eq 1234 -a $NEXP -gt 1111 ] then # kills 5 and 6 - use 1234 but 1256 only if a1 is significant A1COEFF0=246 fi # note mixings for mom=/= 0 (e.p=0 e.e=-1 p.p=m^2) # if p_3 p=(0 0 p3 E) e=(0 0 E p3)/m # # pi g_5 g_5 p.g a1 e.g g_5 rho e.g e.g p.g # b1 e.{epgg} s 1 exotic p.g (e.{epgg} == e.g p.g g_5) # # mom in 3-direction: g_4 g_5, g_3 g_5 are pi and a1 # mom in 3-direction: g_1 g_4, g_1 g_3 are rho' and b1 # mom in 3-direction: g_2 g_4, g_2 g_3 are rho' and b1 # mom in 3-direction: g_3, g_4 are rho and exotic # others invariant (inc g_3 g_4 and g_1 g_2) #--- data prepared for reading in here --------- NBLOCK=1 TEMPPATH=/tmp/cmi if [ $IFK -eq 4 ] then INFO=" 24^3 48 rand .160856 mu=0.004 every 5 (*2 traj) nlong=6" ### see tmzc_old_ifk for other data sets (cyclic t source) LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=uurandom004 if [ $NOTAR -eq 0 ] then rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm fi cd $USERSDIR #DATAPATH=$TEMPPATH/outprcn000000uu DATAPATH=$TEMPPATH/outprc0404uu ## care VWI problem gauge 0990 - so remove rm ${DATAPATH}*0990 # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=900 # IGAUGE1=3400 IGSTEP=1 NGAUGE=4500 # NGAUGE=2300 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=8 # NBLOCK=32 fi if [ $IFK -eq 3 ] then INFO=" 32^3 64 twist 4.05 k=0.157025 mu=0.003 every 5 traj nlong=6" LX4VAL=64 DATAPATH=${USERSDIR}/carsten_3/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=100 IGSTEP=1 NGAUGE=95 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=64 # 109-112 ITSTEP=32 -- KLUDGE ITS0=0 fi if [ $IFK -eq 323 ] then INFO="32^3 64 twist 4.05 k=0.15701 mu=0.003 r2*5 traj nlong=7" LX4VAL=64 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=data4.05mu0.003rd cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR FNAME=data4.05mu0.003-2rd cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv NGAUGE=2000 NBLOCK=8 # time sources ITS0 to LX4VAL in steps of ITSTEP IGAUGE1=0 IGSTEP=1 ITSTEP=1 ITS0=0 fi if [ $IFK -eq 3231 ] then INFO="32^3 64 k=1 rand 4.05 k=0.15701 mu=0.003 every 10*5 traj nlong=7" LX4VAL=64 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=data4.05mu0.003conn-k cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprck0303dd NGAUGE=2000 NBLOCK=1 # time sources ITS0 to LX4VAL in steps of ITSTEP IGAUGE1=1 IGSTEP=1 ITSTEP=1 ITS0=0 fi if [ $IFK -eq 3232 ] then INFO="32^3 64 OLD 4.05 k=0.15701 mu=0.003 every 5 traj nlong=7" LX4VAL=64 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=data4.05m0.003 cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) # IGAUGE1=55 #### from 55-92 had nlong=6 so skip them IGAUGE1=93 # to check effect of equilibration # IGAUGE1=160 IGSTEP=1 NGAUGE=2000 NADD=1000 ### from shiftg.j ### NG=$IGAUGE1 while [ $NG -le $NADD ] do for IT in 00 32 do FFR=outprcv.${IT} if [ $NG -lt 100 ] then NGF=00$NG else # if 100-999 NGF=0$NG fi FF=${FFR}.${NGF} # echo $FF if [ -s $FF ] then NGP=`expr $NG + $NADD ` if [ $NGP -lt 1000 ] then NGPF=0$NGP else NGPF=$NGP fi ##### rename if used##### mv $FF $FFR.${NGPF} fi done NG=`expr $NG + $IGSTEP ` done ### ### FNAME=data4.05mu0.003-2 cp $DPATH/$FNAME.tgz $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR #DATAPATH=$TEMPPATH/outprcv DATAPATH=$TEMPPATH/outprcv NBLOCK=32 # time sources ITS0 to LX4VAL in steps of ITSTEP IGAUGE1=0 ITSTEP=32 ITS0=0 fi if [ $IFK -eq 206 ] then INFO="20^3 48 twist 4.05 k=0.15701 mu=0.006 every 20 traj nlong=6" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=data20at006 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv IGAUGE1=0 IGSTEP=1 NGAUGE=5000 NBLOCK=4 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 fi if [ $IFK -eq 246 ] then INFO="24^3 48 twist 4.05 k=0.15701 mu=0.006 every 5/10 traj nlong=7" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=data4.05mu0.006-L24 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar FNAME=data4.05mu0.006-2-L24 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv IGAUGE1=0 IGSTEP=1 NGAUGE=5000 NBLOCK=4 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 fi if [ $IFK -eq 326 ] then INFO="32^3 64 tw 4.05 k=0.15701 mu=0.006 r2*5/10*1 traj nlong=7" LX4VAL=64 DPATH=${USERSDIR} #FNAME=data4.05m0.006 # care 2 sets - measured differently # .006-10 .006-2-10 measured random t every 10 # .006 older fixed/cyclic t 0/32 0..400 every 5 # here 1-400 1710-5090 # disc seems to have 14-394 and 1710-4960 TEMPPATH=/tmp/cmi FNAME=data4.05mu0.006-10 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR FNAME=data4.05mu0.006-2-10 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm ## kill NaN rm o*4460 cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=0 IGSTEP=1 NGAUGE=5100 NBLOCK=8 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 fi if [ $IFK -eq 328 ] then INFO="32^3 64 twist 4.05 k=0.15701 mu=0.008 r10*1 traj nlong=7" LX4VAL=64 DPATH=${USERSDIR} FNAME=data4.05mu0.008-10 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv #DATAPATH=${USERSDIR}/carsten_323/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) # IGAUGE1=725 IGAUGE1=500 IGSTEP=10 NGAUGE=2000 NBLOCK=8 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 fi if [ $IFK -eq 3212 ] then INFO="32^3 64 tw 4.05 k=0.15701 mu=0.012 r5/10*1 traj nlong=7" LX4VAL=64 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi cd $TEMPPATH FNAME=data4.05mu0.012-10 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar FNAME=data4.05mu0.012-5 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar FNAME=data4.05mu0.012-2-10 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv IGAUGE1=500 IGSTEP=5 NGAUGE=2000 NBLOCK=16 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 fi if [ $IFK -eq 85 ] then INFO=" 24^3 48 Combined 3.75 k=0.1660 mu=0.0085 every 5 traj nlong=5" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi FNAME=data3.75m0.0085 cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar FNAME=data3.75m0.0085-2 cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv IGAUGE1=3005 IGSTEP=1 NGAUGE=3000 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=6 fi if [ $IFK -eq 851 ] then INFO=" 24^3 48 BGL 3.75 k=0.1660 mu=0.0085 every 5 traj nlong=5" LX4VAL=48 DPATH=${USERSDIR} FNAME=data3.75m0.0085-2 TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=4006 IGSTEP=1 NGAUGE=1000 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=6 fi if [ $IFK -eq 852 ] then INFO=" 24^3 48 QCDOC 3.75 k=0.1660 mu=0.0085 every 5 traj nlong=5" LX4VAL=48 DPATH=${USERSDIR} FNAME=data3.75m0.0085 TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=3005 IGSTEP=5 NGAUGE=201 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=24 ITS0=0 fi if [ $IFK -eq 80001 ] then # was used for spain alone - fix FNAME INFO="24^3 48 twist .160856 mu=0.0085 every 4 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} FNAME=corr.. TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tar.bz2 $TEMPPATH cd $TEMPPATH bunzip2 $FNAME.tar.bz2 tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=700 IGSTEP=5 NGAUGE=1000 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 fi if [ $IFK -eq 81 ] then ## LL LF only INFO="24^3 48 k=1 LL LF .160856 mu=0.0085 5 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} FNAME=data3.9mu0.0085nk TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar # wrong length cd $USERSDIR DATAPATH=$TEMPPATH/outprck085085dd ## also could use DATAPATH=$TEMPPATH/outprc085085uu mom=0 # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=750 IGSTEP=10 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=1 fi if [ $IFK -eq 8 ] then INFO="24^3 48 twist.160856 mu=0.0085 every 4(*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi cd $TEMPPATH FNAME=corr-24to3x48-b3.9-mu0.0085 # 700-1300 #rcp sunc.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.bz2 $TEMPPATH rcp sund.amtp.liv.ac.uk:/export/home/sund/cmi/connected/$FNAME.tar.bz2 $TEMPPATH #scp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.bz2 $TEMPPATH bunzip2 $FNAME.tar.bz2 tar -xf $FNAME.tar IGAUGE1=700 IGSTEP=1 NGAUGE=2000 NADD=2000 ### from shiftg.j ### NG=$IGAUGE1 while [ $NG -le $NADD ] do for IT in 00 12 24 36 do FFR=outprcv.${IT} if [ $NG -lt 1000 ] then NGF=0$NG else NGF=$NG fi FF=${FFR}.${NGF} # echo $FF if [ -s $FF ] then NGP=`expr $NG + $NADD ` if [ $NGP -lt 1000 ] then NGPF=0$NGP else NGPF=$NGP fi mv $FF $FFR.${NGPF} fi done NG=`expr $NG + $IGSTEP ` done ### ### echo "first part read in" FNAME=data3.9mu0.0085 #500-2600 #cp $DPATH/$FNAME.tgz $TEMPPATH #rcp sunc.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH rcp sund.amtp.liv.ac.uk:/export/home/sund/cmi/connected/$FNAME.tgz $TEMPPATH #scp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=500 IGSTEP=1 NGAUGE=4000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=8 NBLOCK=32 echo " all read in to /tmp/cmi" fi if [ $IFK -eq 10 ] then INFO="24^3 48 twist .160856 mu=0.010 every 5 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi if [ $NOTAR -eq 0 ] then FNAME=data3.9m0.01 cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar # wrong length fi cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) # 340 (1) 2744 exists on grid IGAUGE1=340 IGSTEP=1 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=4 NBLOCK=32 fi if [ $IFK -eq 15 ] then INFO=" 24^3 48 twist .160856 mu=0.015 every 5 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) # each 1000 - 2100 IGAUGE1=1000 IGSTEP=1 NGAUGE=2200 NADD=2200 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 if [ $NOTAR -eq 0 ] then ## problem - 2 sequeneces with same filenames - so add $NADD to one FNAME=data3.9mu0.015-s1 cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ### from shiftg.j ### ASSUME NG > 999 < 10000 NG=$IGAUGE1 while [ $NG -le $NGAUGE ] do for IT in 00 12 24 36 do FFR=outprcv.${IT} FF=${FFR}.${NG} # echo $FF if [ -s $FF ] then NGP=`expr $NG + $NADD ` mv $FF $FFR.${NGP} fi done NG=`expr $NG + $IGSTEP ` done ### ### cd $USERSDIR FNAME=data3.9mu0.015-s2 cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm fi cd $USERSDIR DATAPATH=$TEMPPATH/outprcv NGAUGE=`expr $NGAUGE + $NADD ` # NBLOCK=4 NBLOCK=32 # NBLOCK=8 fi if [ $IFK -eq 64 ] then INFO=" 24^3 48 twist .160856 mu=0.0064 every 8 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} FNAME=data3.9mu0.0064 TEMPPATH=/tmp/cmi # also 242 every 10 random t-slice done 11-08 #FNAME=data0064uu if [ $NOTAR -eq 0 ] then rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar # wrong length fi cd $USERSDIR DATAPATH=$TEMPPATH/outprcv #DATAPATH=$TEMPPATH/outprc00640064uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=1000 IGSTEP=1 NGAUGE=4000 # each t skips one forward # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP # ITSTEP=6 ITSTEP=1 ITS0=0 # NBLOCK=8 NBLOCK=32 fi if [ $IFK -eq 99 ] then INFO=" 16^3 32 test" LX4VAL=32 DPATH=${USERSDIR} FNAME=data3.9_1632 TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=600 IGSTEP=1 NGAUGE=250 # each t skips one forward # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=4 ITS0=0 NBLOCK=8 fi if [ $IFK -eq 77 ] then INFO=" 16^3 32 test STOUT" LX4VAL=32 DPATH=${USERSDIR} #FNAME=data.stout FNAME=stout.piplus TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar $TEMPPATH cd $TEMPPATH tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=500 IGSTEP=2 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 66 ] then INFO=" 16^3 32 3.75 .1512 006 " LX4VAL=32 DPATH=${USERSDIR} FNAME=data3.9mu0.0085nk #Q# diff on diff ux? TEMPPATH=/tmp/cmi FNAME=outprc0606uu rcp sune.amtp.liv.ac.uk:/home1/users/mcneile/TWISTED/beta3.75_H0.006_16x32_kappa0.1512/pi0_conn_beta3.75_H0.006_16x32_kappa0.1512/corr/${FNAME}* $TEMPPATH cd $USERSDIR DATAPATH=$TEMPPATH/$FNAME # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=700 IGSTEP=5 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 # choose nblock=1 to dump a4 info NBLOCK=2 fi if [ $IFK -eq 776 ] then INFO=" 24^3 48 test STOUT" LX4VAL=48 DPATH=${USERSDIR} FNAME=stout38_1476_006 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar $TEMPPATH cd $TEMPPATH tar -xf $FNAME.tar find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=600 IGSTEP=5 NGAUGE=1000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=1 fi if [ $IFK -eq 775 ] then INFO=" 24^3 48 test STOUT" LX4VAL=48 DPATH=${USERSDIR} FNAME=stout38_1475_006 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar $TEMPPATH cd $TEMPPATH tar -xf $FNAME.tar find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=600 IGSTEP=5 NGAUGE=1000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=1 fi if [ $IFK -eq 27 ] then INFO="3.9 24^348 mu=.27 charm NL=2" LX4VAL=48 DPATH=${USERSDIR} FNAME=charm TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar $TEMPPATH cd $TEMPPATH tar -xf $FNAME.tar rm out*1386 cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=1000 IGSTEP=2 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 60 ] then #INFO="20^3 48 twist 3.8 .16428 mu=0.006 every 5 (*2 traj) nlong=6" INFO="20^3 48 twist 3.8 .16415 mu=0.006 every 5 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} #FNAME=data3.8mu0.006 FNAME=data3.8k0.16415mu0.006 TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar # wrong length cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=0 IGSTEP=1 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 624 ] then INFO="24^3 48 twist 3.8 .164111 mu=0.006 every 2 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} FNAME=data38_24_006 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar # wrong length cd $USERSDIR DATAPATH=$TEMPPATH/outprcn006uu # to match gauges with v=/=s # v=/=s diff random t-vals #DATAPATH=$TEMPPATH/outprc164099006uu # v=/=s same random t-vals as s=v #DATAPATH=$TEMPPATH/outprc164120006uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=300 IGSTEP=4 NGAUGE=2500 # NGAUGE=69 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 # NBLOCK=1 fi if [ $IFK -eq 824 ] then INFO="24^3 48 twist 3.8 .164111 mu=0.008 every 2 (*2 traj) nlong=6" LX4VAL=48 DPATH=${USERSDIR} ## two series 1098-1846 (4) // FNAME=data3.8k0.164111mu0.008L24 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar NADD=1000 for FN in out* do NG=`echo $FN | cut -f 3 -d "."` TV=`echo $FN | cut -f 2 -d "."` ROOT=`echo $FN | cut -f 1 -d "."` NGP=`expr $NG + $NADD ` if [ $NGP -lt 1000 ] then NGPF=0$NGP else NGPF=$NGP fi #echo $ROOT.$TV.$NG $ROOT.$TV.$NGP mv $ROOT.$TV.$NG $ROOT.$TV.$NGP # done FNAME=data3.8k0.164111mu0.008L24-2 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar # wrong length cd $USERSDIR DATAPATH=$TEMPPATH/outprc0808uu # 1002-1802 (4) or so # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=1002 IGSTEP=4 NGAUGE=2500 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 6 ] then INFO="20^3 48 twist 3.8 0.164111 mu=0.006 every 5 (*2 traj) nlong=6" TEMPPATH=/tmp/cmi LX4VAL=48 DPATH=${USERSDIR} if [ $NOTAR -eq 0 ] then #FNAME=MN20to3x48-Kc0.164111-0.0060_29mar07 FNAME=mn.20to3x48-Kc0.164111-0.0060-0715-1757 # 400-2300 (+2000) rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar # care overlapping gauge numbers NADD=1000 for FN in out* do NG=`echo $FN | cut -f 3 -d "."` TV=`echo $FN | cut -f 2 -d "."` ROOT=`echo $FN | cut -f 1 -d "."` NGP=`expr $NG + $NADD ` if [ $NGP -lt 1000 ] then NGPF=0$NGP else NGPF=$NGP fi #echo $ROOT.$TV.$NG $ROOT.$TV.$NGP mv $ROOT.$TV.$NG $ROOT.$TV.$NGP # done ####FNAME=BGL20to3x48-Kc0.164111-0.0060_29mar07 FNAME=bgl.20to3x48-Kc0.164111-0.0060-0715-1757 ## 700-1750 ## 2270-2720 (5) BGL started from end of MN run # outprcn006uu #FNAME=data38_20_006 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar fi cd $USERSDIR ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv #DATAPATH=$TEMPPATH/outprcn006uu IGAUGE1=600 IGSTEP=1 NGAUGE=5000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=8 fi if [ $IFK -eq 61 ] then #INFO="20^3 48 twist 3.8 .16415 mu=0.006 every 5 (*2 traj) nlong=6" #FNAME=outprcvb38mu006k16415 #INFO="20^3 48 twist 3.8 .164084 mu=0.006 every 5 (*2 traj) nlong=6" #FNAME=b38k164084 INFO="20^3 48 twist 3.8 0.164099 mu=0.006 every 5 (*2 traj) nlong=6" TEMPPATH=/tmp/cmi LX4VAL=48 DPATH=${USERSDIR} # 300-1300 0 12 24 36 cyclic FNAME=data3.8mu0.006-2 cp $DPATH/$FNAME.tgz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH #gunzip $FNAME.tar.gz #tar -xf $FNAME.tar cd $USERSDIR # random time 500-4100 FNAME=data3.8mu0.006rd cp $DPATH/$FNAME.tgz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=80 # safer to start at 500 IGAUGE1=500 IGSTEP=1 NGAUGE=5000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=8 fi if [ $IFK -eq 9 ] then INFO="20^3 48 twist 3.8 0.164099 mu=0.009 every 5 (*2 traj) nlong=6" # 186-1814 0 12 24 36 cyclic FNAME=data3.8mu0.009 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=100 IGSTEP=1 NGAUGE=4000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=8 fi if [ $IFK -eq 11 ] then INFO="24^3 48 twist 3.8 0.164111 mu=0.011 every 4 (*2 traj) nlong=6" # random every 4 46- FNAME=data38_24_011 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar # corrupt (now removed from tar file) #rm outprcn011uu.17.1518 #rm outprcn011uu.??.1526 #rm outprcn011uu.??.1634 #rm outprcn011uu.??.2110 ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcn011uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) # series 1 1046-2998 series2 1046 -1906 (added 2000) IGAUGE1=1046 IGSTEP=4 NGAUGE=3000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 165 ] then INFO="24^3 48 twist 3.8 0.164111 mu=0.0165 every 4 (*2 traj) nlong=6" # random every 4 46- FNAME=data3.8k0.164111mu0.0165L24 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) # series 1 1046-2998 series2 1046 -1906 (added 2000) IGAUGE1=1046 IGSTEP=4 NGAUGE=3000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 45 ] then INFO="24^3 48 twist 3.8 0.164111 mu=0.0045 every 5 (*2 traj) nlong=6" # random every 4 46-542 FNAME=data38_24_0045 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcn0045uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=46 IGSTEP=4 NGAUGE=1000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 38 ] then INFO="20^3 48 twist 3.8 0.164111 mu=0.006 every ? nlong=6" FNAME=corr_20to3x48_beta3.8_kappa0.164111_mu0.0060 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi cd $TEMPPATH if [ $NOTAR -eq 0 ] then cp $DPATH/$FNAME.tar.bz2 $TEMPPATH/$FNAME.tar.bz2 bunzip2 $FNAME.tar.bz2 tar -xf $FNAME.tar fi ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=150 IGSTEP=1 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 12 ] then INFO="20^3 48 twist 3.8 0.164099 mu=0.012 every 5 (*2 traj) nlong=6" # 100-2588 0 12 24 36 cyclic FNAME=data3.8mu0.012 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=80 IGSTEP=1 NGAUGE=4000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=12 ITS0=0 NBLOCK=8 fi if [ $IFK -eq 153 ] then INFO="20^3 48 twist 3.8 0.164099 mu=0.015 every 5 (*2 traj) nlong=6" # 100-580 random FNAME=data3.8mu0.015 LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH/$FNAME.tar.gz cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=80 IGSTEP=1 NGAUGE=4000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=8 fi if [ $IFK -eq 26 ] then LX4VAL=48 DPATH=${USERSDIR} FNAME=uu026 INFO="24^3 24 twist 3.9 sea 004 valence 026 nlong=6" TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar # wrong length cd $USERSDIR DATAPATH=$TEMPPATH/outprcn000026uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=900 IGSTEP=1 NGAUGE=3000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=24 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 6661 ] then LX4VAL=32 DPATH=${USERSDIR} FNAME=stout INFO="stout test 16^3 32 k=0.15 mu=0.08" TEMPPATH=/tmp/cmi if [ $NOTAR -eq 0 ] then cp $DPATH/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar # wrong length fi cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=500 IGSTEP=1 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 666 ] then INFO=" 24^3 48 STOUT mul=0.006 dec 2010" LX4VAL=48 DPATH=${USERSDIR} TEMPPATH=/tmp/cmi # kill charged.00.0718 - crazy length if [ $NOTAR -eq 0 ] then FNAME=stout006c.tar rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.gz tar -xf $FNAME #rm charged.00.0718 fi cd $USERSDIR DATAPATH=$TEMPPATH/charged # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=100 IGAUGE1=400 IGSTEP=2 NGAUGE=800 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 # choose nblock=1 to dump a4 info NBLOCK=10 fi if [ $IFK -eq 44444 ] then LX4VAL=48 DPATH=${USERSDIR} #FNAME=datab1.90_k0.163335_mu0.004_musigma0.15_mudelta0.19 #FNAME=data_L24T48_b1.90_k0.163335_mu0.004_musigma0.185_mudelta0.236 # /home1/users/mcneile/TWISTED/pi0_conn_iwasaki_b1.90_L24T48_k0.16327_mu0.004_musigma0.15_mudelta0.19/pi0 FNAME=outprc0404uu INFO="2+1+1 test 004 .16327" TEMPPATH=/tmp/cmi #pi0_conn_iwasaki_b1.90_L24T48_k0.16327_mu0.004_musigma0.15_mudelta0.19 rcp sune.amtp.liv.ac.uk:/home1/users/mcneile/TWISTED/pi0_conn_iwasaki_b1.90_L24T48_k0.16327_mu0.004_musigma0.15_mudelta0.19/pi0_old/${FNAME}* $TEMPPATH cd $TEMPPATH # get rid of NaN rm $TEMPPATH/${FNAME}*1400 #gunzip $FNAME.tar.gz #tar -xf $FNAME.tar # wrong length cd $USERSDIR DATAPATH=$TEMPPATH/$FNAME # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=800 IGSTEP=10 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=1 fi if [ $IFK -eq 444 ] then LX4VAL=48 DPATH=${USERSDIR} INFO=" 24^3 48 2+1+1 every 10 (*2) 004" TEMPPATH=/tmp/cmi DATAPATH=$TEMPPATH/$FNAME # run1 FNAME=data211004uu.tar # run2 FNAME=data211004run2uu.tar rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.gz tar -xf $FNAME cd $USERSDIR DATAPATH=$TEMPPATH/outprc004004uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=10 IGSTEP=10 NGAUGE=400 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 446 ] then LX4VAL=48 DPATH=${USERSDIR} INFO=" 24^3 48 2+1+1 every 5 (*2) 006" TEMPPATH=/tmp/cmi DATAPATH=$TEMPPATH/$FNAME FNAME=data211006uu.tar rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.gz tar -xf $FNAME cd $USERSDIR DATAPATH=$TEMPPATH/outprc006006uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=600 IGSTEP=5 NGAUGE=600 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 448 ] then LX4VAL=48 DPATH=${USERSDIR} INFO=" 24^3 48 2+1+1 every 5 (*2) 008" TEMPPATH=/tmp/cmi #FNAME=outprc0808uu #rcp sune.amtp.liv.ac.uk:/home1/users/mcneile/TWISTED/pi0_conn_L24T48_b1.90_k0.16326_mu0.008_musigma0.15_mudelta0.19/corr/${FNAME}* $TEMPPATH DATAPATH=$TEMPPATH/$FNAME FNAME=data211008uu.tar rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.gz tar -xf $FNAME cd $USERSDIR DATAPATH=$TEMPPATH/outprc008008uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=700 IGSTEP=5 NGAUGE=3300 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 455 ] then LX4VAL=64 DPATH=${USERSDIR} INFO=" 32^3 64 2+1+1 every 5 (*2) 0055 b=1.95" TEMPPATH=/tmp/cmi DATAPATH=$TEMPPATH/$FNAME FNAME=data1.95mu0.0055L32 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/charged # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=600 IGSTEP=10 NGAUGE=330 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 435 ] then LX4VAL=64 DPATH=${USERSDIR} INFO=" 32^3 64 2+1+1 every 10 (*2) 0055 b=1.95 KO" TEMPPATH=/tmp/cmi DATAPATH=$TEMPPATH/$FNAME FNAME=B35.32_L32_T64_beta195_mul0035_musig135_nudel170_kappa1612400.light.conn rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/charged # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=1000 IGSTEP=10 NGAUGE=330 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=2 fi if [ $IFK -eq 456 ] then LX4VAL=48 DPATH=${USERSDIR} INFO=" 24^3 48 2+1+1 STOUT every 10 (*2) 006" TEMPPATH=/tmp/cmi DATAPATH=$TEMPPATH/$FNAME FNAME=data211s006uu.tar rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.gz tar -xf $FNAME cd $USERSDIR DATAPATH=$TEMPPATH/outprc006006uu # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=1000 IGSTEP=10 NGAUGE=400 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=1 fi if [ $IFK -eq 41 ] then ## note dd004k1_ave.tar.gz has 30 gauges ???6 with only mom ave for J=1 ## also 1775 has mom_ave only - so remove #see ifk=4 for more (dd) ones LX4VAL=48 DPATH=${USERSDIR} #FNAME=dd004k1 #INFO="24^3 48 3.9 dd 004 k=1 LL LF nlong=6" INFO="24^3 48 3.9 dd 004 random t k=1 LL LF nlong=6" FNAME=dd1random004 TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm # see above - for use with rho #rm outprck000000dd.??.1775 #rm outprck000000dd.??.???6 cd $USERSDIR #DATAPATH=$TEMPPATH/outprck000000dd DATAPATH=$TEMPPATH/outprck0404dd # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=900 IGSTEP=1 # test NGAUGE=4500 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 432 ] then LX4VAL=64 DPATH=${USERSDIR} FNAME=data3.9mu0.004L32 ## 2 series - every 2 saved - random t0: now all combined # care some 50 at end were dud (all same) last disc is 2931 # was to 3299 diff t-slices so hard to check. STOP at 3241 INFO="32^3 64 3.9 mu=.004 random t_s every 4 confs nlong=6" TEMPPATH=/tmp/cmi cp $DPATH/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar #cd $USERSDIR #FNAME=data3.9mu0.004L32-2 #cp $DPATH/$FNAME.tgz $TEMPPATH #cd $TEMPPATH #gunzip $FNAME.tgz #tar -xf $FNAME.tar # corrupt #rm outprcv.??.2305 #rm outprcv.??.2309 #rm outprcv.??.23[1-4]? rm outprcv.??.2353 # wrong update - so delete measurements (last 100 trajs) rm outprcv.??.33?? cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=600 # to avoid initialisation - could start at 800 not 600 # IGAUGE1=800 # IGAUGE1=1000 # IGAUGE1=1250 IGSTEP=1 NGAUGE=4000 ## max 3241 so 3241-600 change 3-7-07 NGAUGE=2642 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=16 fi if [ $IFK -eq 332 ] then LX4VAL=64 DPATH=${USERSDIR} #FNAME=data3.9mu0.003L32 INFO="32^3 64 3.9 mu=.003 random t_s every 4 confs nlong=6" TEMPPATH=/tmp/cmi FNAME=mu0.003-32 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar FNAME=mu0.003-2-32 rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tar.gz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tar.gz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/outprcv # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=0 IGSTEP=1 NGAUGE=2000 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 NBLOCK=4 fi if [ $IFK -eq 482 ] then INFO="48^3 96 twist 4.2 k=0.15 mu=0.002 r10*1 traj nlong=8" LX4VAL=96 DPATH=${USERSDIR} FNAME=data4.2mu0.002connected TEMPPATH=/tmp/cmi rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/$FNAME.tgz $TEMPPATH cd $TEMPPATH gunzip $FNAME.tgz tar -xf $FNAME.tar ## kill length 85 cases find * -size 85c | xargs rm cd $USERSDIR DATAPATH=$TEMPPATH/charged # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=0 IGSTEP=5 NGAUGE=200 NBLOCK=1 # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 fi if [ $IFK -eq 2112 ] then INFO=" 48^3 96 2+1+1 IWA mu 002 b2.1 " LX4VAL=96 DPATH=${USERSDIR} #Q# diff on diff ux? TEMPPATH=/tmp/cmi FNAME=datab2.1k0.156357mu0.002conn-charged rcp sune.amtp.liv.ac.uk:/ukqcd6/cmi/connected/${FNAME}.tgz $TEMPPATH/${FNAME}.tar cd $TEMPPATH #gunzip $FNAME.tar.gz tar -xf $FNAME.tar cd $USERSDIR DATAPATH=$TEMPPATH/charged # gauges IGAUGE1 in steps of IGSTEP for NGAUGE cases (see below to skip) IGAUGE1=500 IGSTEP=10 NGAUGE=200 # NOW USES "INQUIRE" --- run copyr in directory to clean up # time sources ITS0 to LX4VAL in steps of ITSTEP ITSTEP=1 ITS0=0 # choose nblock=1 to dump a4 info NBLOCK=1 fi LX4H=`expr $LX4VAL / 2 - 1 ` if [ $IT2 -gt $LX4H ] then IT2=$LX4H fi echo "particle="$PARTICLE #------ end of parameters ------------------------- ## as is runs in temp directory TDIR=cmi TEMPDIR=/tmp/$TDIR cd /tmp if [ -d $TDIR ] then echo " " else mkdir $TDIR fi cd $TEMPDIR if [ -f tmprf.f ] then rm tmprf.f fi #pn=$$ pn="zfit" echo `date` # (no " quotes on here document target causes shell substitution) cat > $TEMPDIR/fitex.h << EOFF c SET DATA STORAGE c npst itst nlst used overall - match length of data statements in TFIT c and MAIN to npst and in TMAT to npstph c PARAMETER (nlst=3450,npsth=21,itsth=33) c for beta=4.2 itmax=48 so increase PARAMETER (nlst=3450,npsth=21,itsth=49) PARAMETER (NPARSv=28, npstph=6, nrsth=2) c SET TIME EXTENT (used cfs in cosh fits ) PARAMETER (lx4val=$LX4VAL) c to activate kludges if required PARAMETER(ikludge=$A1COEFF0) EOFF cat > $TEMPDIR/tmprf.f << EOFH c****************************************************************** c c Program to fit energies to lattice correlations versus t c several different operators can be used (called paths ip) c As written it does analysis of meson data c c Structure: c MAIN determines what data is to be analysed c BOOTV, JACKV is "old code" to look at c effective mases etc. c c TMAT EIGENT BOOTV2 are old code to do Variational c analyses (assuming all*all observables exist) c c TFIT is a bridge between this code and the general c correlated fitting program FITEX c set some fit parameters here (nexp, idv, initial values) c (see GRAPH for vertical range and t-offset by hand) c c FITEX, CORFIT are documented in prologues to code c c****************************************************************** c MAIN PROGRAM (FORTRAN) c ************************************ double precision PLAV IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' PARAMETER (ntypevs=21,nfz=1,nmeson=40) parameter (lx4v=lx4val,lx4h=lx4v/2) parameter (nobss=(1+nfz+nfz*nfz)*4) parameter (npstp=npstph,itst=itsth,nrst=nrsth,npst=npsth,iprin=0) COMMON/EPP/PLAV(NPSTp,NPSTp,ITST,NLST,NRST),nconf common/fact/npv,ip1s(npst),ip2s(npst),ifs(npst), nfacp c COMMON/measur/c(nlst,0:lx4h,ntypevs,nobss), & ec(nlst,0:lx4h,ntypevs,nobss),ct(nlst,0:lx4h) Dimension cll(nlst,0:lx4h,ntypevs,nobss) CHARACTER*120 CONF LOGICAL there c info to TMAT here (not all used) COMMON/INFO/NLAT,NR,NRMN,NPVV, ITMAX,ITMAXP,ITVMAX,ITBIAS,IPRINv Dimension pcac(nlst) Integer ip1v(npst), ip2v(npst),imeson1(nmeson), & imeson2(nmeson),ifimag(nmeson), igaugest(nlst) c PA more accurate than AP c to average t>0 and t<0: c order PP PA AP AA 44 P4 4P A4 4A for pion like P=g5 A=g4g5 4=g4 c order 44 VV AA 4V V4 4A A4 VA AV for rho-a1 like 4=gig4 V=gi A=gig5 c order BB SS - total 20 B=gig4g5 S=I c correction 8-1-08 order of last two is SS BB c best choice is weaker coupling at sink - ie PA rather than AP c order of magnitudes P > 4 > A (4 mixes A) c order of magnitudes 4 ~ A > V (A mixes V) dimension corl(nlst, 0:lx4h, 2,3,2,3), dr(nlst),drr(nlst) c antisymmetric in time dimension sf(ntypevs), c00(nlst) c DATA sf /1.,-1.,-1.,1.,1.,-1.,-1.,-1.,-1., c & 1.,1.,1.,-1.,-1.,-1.,-1.,-1.,-1., 1. , 1. ,1. / c ifk=60 seems sf is different (ie commute g_5 each end so g5.g5g4;g5.g4 c etc ) c 8-11-06 DATA sf /1.,-1.,-1.,1.,1.,-1.,-1.,1.,1., & 1.,1.,1.,-1.,-1.,-1.,-1.,1.,1., 1. , 1. ,1. / c negative diagonal cases to fix: g4.g4 a1.a1 b1.b1 s.s - commute g_4 c checked by observation c 13-3-07 004 seems AA is negative - but 015 and 432 and ifk=6 favour c original sign - small effect on chi2 anyway c dimension sdv(ntypevs) c original DATA sdv /4*1.0,-1.0,6*1.0,-1.0,6*1.0,-1.0,-1.0,1.0/ c DATA sdv /3*1.0,2*-1.0,6*1.0,-1.0,6*1.0,-1.0,-1.0,1.0/ c order of itype for pion-like integer ipion1(9),ipion2(9) Data ipion1 /1,1,2,2,3,1,3,2,3/ Data ipion2 /1,2,1,2,3,3,1,3,2/ c order of itype for rho-like integer irho1(9),irho2(9) Data irho1 /1,2,3,1,2,1,3,2,3/ Data irho2 /1,2,3,2,1,3,1,3,2/ c order used in fits Data ip1v /1,1,2,1,2,3,1,2,3,4,1,2,3,4,5,1,2,3,4,5,6 / Data ip2v /1,2,2,3,3,3,4,4,4,4,5,5,5,5,5,6,6,6,6,6,6 / c files from mallprc - for info only here - numbers as cmi gamma's c but later 3 components averaged for vectors DATA imeson1 /1,1,5,5,13,1,13,5,13, 6,7,8,2,3,4,10,11,12, &6,7,8,2,3,4, 6,7,8,10,11,12, 2,3,4,10,11,12, &9,14,15,16/ DATA imeson2 /1,5,1,5,13,13,1,13,5, 6,7,8,2,3,4,10,11,12, &2,3,4,6,7,8, 10,11,12,6,7,8, 10,11,12,2,3,4, &9,14,15,16/ c those allowed by twisting (proportional to i mu ) DATA ifimag /0,0,0,0,0,1,1,1,1, 0,0,0,0,0,0,0,0,0, &0,0,0,0,0,0, 1,1,1,1,1,1, 1,1,1,1,1,1, 0,0,0,0 / c pi 1 5 13 as 1 3 5 (ie L) +1 for F c see later problem about off-diag signs for pion iout=77 ifk=$IFK write(*,*) & '${INFO}' write(*,*)' particle=${PARTICLE}' c NORM still -/+1 -/+1 in source, no (1/2K), CU inverter. so needs xnorm=1/2 itstep=$ITSTEP its0=$ITS0 igauge1=$IGAUGE1 ngauge=$NGAUGE igstep=$IGSTEP ntype=20 c iobval=1 is PION 2 is rho-a1 3 is a0 4 is b1 c pi fit 6-23 nexp=1 OK c rho-a1 fit 4-15 nexp=3 OK Care KLUDGE in funct1 in fitex.f to force c search "ne funct1" x7 x8 ==0 (a1) iobval=$IOBVAL c T RANGE it1=$IT1 it2=$IT2 ipr=0 igg=0 iga=0 igaugef=igauge1+(ngauge-1)*igstep Do igaugev=igauge1,igaugef,igstep igauge=igaugev lx4m=lx4v-1 c*** treat each its as different gauge - so far itstepv=itstep c KLUDGES to read in "erratic" data sets c TWIST - options to skip missing gauges , vary T0 etc c if(igauge .ne. 340) itstepv=itstep1 c if(igaugev .ge. 71) igauge=igaugev+2 c if(ifk .eq. 2 .and. igaugev .gt. 1335) igauge=igaugev+10 if ( ifk .eq. 1 ) then if ( igauge .le. 134) then itstepv=24 lx4m=25 its0=0 else itstepv=64 lx4m=63 its0=0 endif else if ( ifk .eq. 3 ) then if( igaugev .gt. 135) igauge=igaugev+5 if ( igauge .le. 121) then itstepv=24 lx4m=25 its0=0 else if( igauge .le. 124) then itstepv=64 lx4m=63 its0=24 else if( igauge .le. 128) then itstepv=24 lx4m=25 its0=0 else if( igauge .le. 135) then itstepv=64 lx4m=63 its0=24 else itstepv=64 lx4m=63 its0=0 endif endif c write(*,*) igauge,its0,itstepv,igaugev,lx4m, lx4val Do its=its0,lx4m,itstepv igaugex=igauge write(conf,906)its,igaugex 906 format('${DATAPATH}.',i2.2,'.',i4.4) if (ifk .eq .77 .or. ifk .eq. 27 .or .ifk .eq. 776 & .or. ifk .eq. 775) then c to cope with Craig's numbering of files if( its .le. 9 .and. igaugex .le. 999) write(conf,9061)its,igaugex if( its .le. 9 .and. igaugex .gt. 999) write(conf,9062)its,igaugex if( its .gt. 9 .and. igaugex .le. 999) write(conf,9063)its,igaugex if( its .gt. 9 .and. igaugex .gt. 999) write(conf,9064)its,igaugex 9061 format('${DATAPATH}.',i1,'.',i3) 9062 format('${DATAPATH}.',i1,'.',i4) 9063 format('${DATAPATH}.',i2,'.',i3) 9064 format('${DATAPATH}.',i2,'.',i4) endif c write(*,987) CONF inquire (file=conf, exist=there) if(there) then igg=igg+1 write(*,987) CONF 987 format(' input file is ',A) open(iout,file=conf) c read(iout,901)igauge0,it,ixstep,nobsv,lx4hv,nmeasv,lev,fp,nlong, c lx4hv corrupt in some early file outputs c add g5mu1 for non-deg cases (some only) if (ifk .eq. 482 ) then c no header - set by hand coupK=0.154073 lx4=lx4val lx1=lx4/2 nlong=8 nobsv=8 igauge0=igaugex else if (ifk .eq. 2112 ) then c no header - set by hand coupK=0.156357 lx4=lx4val lx1=lx4/2 nlong=7 nobsv=8 igauge0=igaugex else if (ifk .eq. 455 ) then c no header - set by hand coupK=0.1612360 lx4=lx4val lx1=lx4/2 nlong=7 nobsv=8 igauge0=igaugex else if (ifk .eq. 435 ) then c no header - set by hand coupK=0.16124 lx4=lx4val lx1=lx4/2 nlong=7 nobsv=8 igauge0=igaugex else if (ifk .eq. 666 ) then c no header - c coupK=0.13570 coupK=0.13563 g5mu=0.006 lx4=lx4val lx1=lx4/2 nlong=5 nobsv=8 igauge0=igaugex else if ( ifk .ne. 44400 ) then read(iout,9017) igauge0,it,ixstep,nobsv,nmeasv,lev,fp,nlong, & lx1,lx2,lx3,lx4,coupK,g5mu 901 format(7i5,f8.4,5i5,2f8.5) 9017 format(4i5,5x,2i5,f8.4,5i5,2f8.5) else read(iout,*) igauge0,it,lx1,lx2,lx3,lx4,coupK c gauge varied length? nobsv=8 lx4hv=lx4/2 9617 format(i3,5i3,f9.5) endif endif igaugest(igg)=igauge0 g5mu=g5mu/(2.0*coupk) c BUGFIX (and get precision correct in f_pi) if( ifk .eq. 4) g5mu=0.004 if( ifk .eq. 432) g5mu=0.004 if( ifk .eq. 332) g5mu=0.003 if( ifk .eq. 41) g5mu=0.004 if( ifk .eq. 64) g5mu=0.0064 if( ifk .eq. 8) g5mu=0.0085 if( ifk .eq. 10) g5mu=0.01 if( ifk .eq. 15) g5mu=0.015 if( ifk .eq. 99) g5mu=0.0075 if( ifk .eq. 323) g5mu=0.003 if( ifk .eq. 3231) g5mu=0.003 if( ifk .eq. 3232) g5mu=0.003 if( ifk .eq. 326) g5mu=0.006 if( ifk .eq. 246) g5mu=0.006 if( ifk .eq. 206) g5mu=0.006 if( ifk .eq. 328) g5mu=0.008 if( ifk .eq. 3212) g5mu=0.012 if( ifk .eq. 26) g5mu=0.004 if( ifk .eq. 26) g5mu1=0.026 if( ifk .eq. 6) g5mu=0.006 if( ifk .eq. 45) g5mu=0.0045 if( ifk .eq. 624) g5mu=0.006 if( ifk .eq. 824) g5mu=0.008 if( ifk .eq. 165) g5mu=0.0165 if( ifk .eq. 11) g5mu=0.011 if( ifk .eq. 38) g5mu=0.006 if( ifk .eq. 60) g5mu=0.006 if( ifk .eq. 61) g5mu=0.006 if( ifk .eq. 9) g5mu=0.009 if( ifk .eq. 12) g5mu=0.012 if( ifk .eq. 15) g5mu=0.015 if( ifk .eq. 153) g5mu=0.015 if( ifk .eq. 6661) g5mu=0.08 if( ifk .eq. 666) g5mu=0.006 if( ifk .eq. 444) g5mu=0.004 if( ifk .eq. 446) g5mu=0.006 if( ifk .eq. 456) g5mu=0.006 if( ifk .eq. 448) g5mu=0.008 if( ifk .eq. 77) g5mu=0.006 if( ifk .eq. 66) g5mu=0.006 if( ifk .eq. 776) g5mu=0.006 if( ifk .eq. 775) g5mu=0.006 if( ifk .eq. 27) g5mu=0.27 if( ifk .eq. 482) g5mu=0.002 if( ifk .eq. 2112) g5mu=0.002 if( ifk .eq. 455) g5mu=0.0055 if( ifk .eq. 435) g5mu=0.0035 if(ifk .ne. 15 .and. ifk .ne. 323 .and. ifk .ne. 8 .and. & ifk .ne. 6 .and. ifk .ne. 824 .and. & igaugex .ne. igauge0) & STOP 'wrong igauge' if(nobsv .gt. nobss) STOP ' nobss too small' if( ipr .eq. 1) then write(6,901) igauge0,it,ixstep,nobsv,lx4h,nmeasv,lev,fp,nlong, & lx1,lx2,lx3,lx4,coupK,g5mu write(6,*) ' K=',coupK,' g5mu=',g5mu write(6,*) ' nlong=', nlong, ' lev=',lev,' fp=',fp ipr=0 endif xnorm=0.5 c *** read in here one its one imeas nobsm=nobsv-1 ntypep=ntype if( ifk .ge. 1 .and. ifk .ne. 482 .and. ifk .ne. 455 & .and. ifk .ne. 435 & .and. ifk .ne. 2112 .and. ifk .ne. 666 ) ntypep=ntype+1 Do itype=1,ntypep Do iobsi=1,nobsm,2 iobs=(iobsi+1)/2 Do idt=0,lx4h if( $IFK .eq. 41 .and. $IOBVAL .eq. 2) then c total then parallel for rho b1 read(iout,9691) itypev,iobsv,idtv,cav1,cav2,cav3,cav4 c t> 0 then t<0 9691 format(3i3,4f20.12) ifparallel=0 if( ifparallel .eq. 1) then c parallel c write(*,*) ' parallel' cav1=cav3 cav2=cav4 else c perpendicular c write(*,*) ' perpendicular' cav1=(cav1-cav3)*0.5 cav2=(cav2-cav4)*0.5 endif else c ------------------regular case ---------------------- if(ifk .ne. 44400 .and. ifk .ne. 44800) then c t> 0 then t<0 if (ifk .eq. 482 .or. ifk .eq. 455 .or. ifk .eq. 435 .or. & ifk .eq. 2112 .or. ifk .eq. 666) then c beta=4.2 1.95 special cases read(iout,9096) itypev,iobsv,idtv,cav1,cav2 9096 format(2i3,i4,2e25.16) c write(6,9096) itypev,iobsv,idtv,cav1,cav2 else read(iout,9091) itypev,iobsv,idtv,cav1,cav2 9091 format(3i3,2f20.12) endif else read(iout,*) itypev,iobsv,idtv,cav1,cav2 c read(iout,9861) itypev,iobsv,idtv,cav1,cav2 9891 format(2i3,i4,2e16.7) endif endif if(itypev .ne. itype) STOP ' wrong itype' if(iobsv .ne. iobsi) STOP ' wrong iobs' if(idtv .ne. idt) STOP ' wrong idt' cll(igg,idt,itype,iobsi)=cav1*xnorm cll(igg,idt,itype,iobsi+1)=cav2*xnorm cav=cav1 if( idt .ne. 0 .and. idt .ne. lx4h) then cav=(cav1+cav2*sf(itype))*0.5 endif c(igg,idt,itype,iobs)=cav*xnorm ec(igg,idt,itype,iobs)=0.0 c save for pcac output versus gauge at its=0 c if(its .eq. 0 .and. idt .eq. 10 .and. itype .eq. 2 if(its .ge. 0 .and. idt .eq. 10 .and. itype .eq. 2 & .and. iobs .eq. 1) then iga=iga+1 c00(iga)=cav*xnorm endif Enddo Enddo Enddo Endif Enddo c end igaugev Enddo c dump to check signs LL iobs=1 do idt=6,12,6 write(6,*)' dump LL idt=',idt do itype=1,20 call avn(cll(1,idt,itype,iobs),cav1,sav1,ngauge) call avn(cll(1,idt,itype,iobs+1),cav2,sav2,ngauge) sav1=sav1/sqrt(float(ngauge-1)) sav2=sav2/sqrt(float(ngauge-1)) write(6,9391)itype,cav1,sav1,cav2,sav2 9391 format(i3,4f15.5) enddo enddo nga=iga do iga=1,nga write(32,9494)iga,c00(iga) enddo c STOP ngauge=igg if( ngauge .gt. nlst) STOP ' nlst too small' c*** conserved vector current c check VWI if( ifk .gt. 0) then itdr=10 xmu=g5mu if(ifk .eq. 26) xmu=(g5mu+g5mu1)*0.5 Do igg=1,ngauge Do idt=1,lx4h idtm=idt-1 iob1=1 c iob1=5 check iob1=5 gives same Z_V but with bigger error c t > t0 tm....t....tp dv=cll(igg,idt,21,iob1)- & cll(igg,idtm,21,iob1) cp=cll(igg,idt,1,iob1) cpmu=cp*2.0*xmu r=-cll(igg,idt,21,iob1)/cll(igg,idt,6,iob1) if(idt .eq. itdr) dr(igg)=dv/cpmu if(idt .eq. itdr) drr(igg)=r c t < t0 tmm....tm....t (now averaged above) iob1=2 dvm=cll(igg,idt,21,iob1)- & cll(igg,idtm,21,iob1) cpm=cll(igg,idt,1,iob1) cpmum=cpm*2.0*xmu if(idt .eq. itdr .and. abs(dr(igg)-1.0) .gt. 0.0001 ) & write(6,969) itdr, igg,igaugest(igg),dv,cpmu, dr(igg) 969 format(' VWI check failure',3i5,6f12.5) Enddo Enddo c STOP c VWI check call avn(dr,a,s,ngauge) write(*,*)' VWI check ', a, s, itdr c Z_V eval c here V is evaluated it to it+1 - so factors exp(-m_pi/2) needed c for comarison with g_4 c call avn(drr,a,s,ngauge) nblock=$NBLOCK call avratb(cll(1,itdr,21,iob1),cll(1,itdr,6,iob1) & ,ngauge,a,s,nblock) a=-a write(*,*)' Z_V, err, t_val ', a, s, itdr write(25,925) itdr,a,s 925 format('Z_V*exp(-m_pi/2) evaluated at t=',i2, & ' is ',f8.5,'+/-',f8.5) c another way to define Z_V is -2 mu /d with d as symmetric diff c #1 / d #6 c ### note c is anti/symmetrized, cll is not idt=11 iob=1 itp2=6 c try 7 instead - error 10* bigger itp2=7 idtp=idt+1 idtm=idt-1 Do igg=1,ngauge pcac(igg)=(-c(igg,idtp,itp2,iob)+c(igg,idtm,itp2,iob))*0.5 drr(igg)=c(igg,idt,1,iob)/pcac(igg) Enddo c call avratb(cll(1,idt,1,iob),pcac(1) c & ,ngauge,a,s,nblock) call avratb(c(1,idt,1,iob),pcac(1) & ,ngauge,a,s,nblock) a=a*g5mu*2.0 s=s*g5mu*2.0 write(*,*)' Z_V, err, t_val ', a, s, idt write(25,9251) idt,a,s 9251 format('Z_V evaluated at t=',i2, & ' is ',f8.5,'+/-',f8.5) call avn(drr(1),a,s,ngauge) a=a*g5mu*2.0 s=s*g5mu*2.0/sqrt(float(ngauge-1)) write(*,*)' Z_V, err, t_val ', a, s, idt endif c**** some PCAC TWIST ratios without fitting **** c do ratios 6 is Pg4 (7 is g4P) 2 is PA 3 is AP: 2 smaller errors itp1=6 itp2=3 itp2=2 c itp2=2 seems to have smaller variance in fort.31 ifk=1 (also 624) if(iprin .gt. 1) write(*,*) itp1,itp2 if(iprin .gt. 1)write(6,*) (c(igg,3,itp1,1),igg=1,20) Do iob=1,4 Do idt=1,23 call avrat(c(1,idt,itp2,iob),c(1,idt,itp1,iob),nlst,ngauge,a,s) if(iprin .gt. 1)write(6,9395) iob, idt, a,s 9395 format(2i5,4f15.5) Enddo Enddo c pcac mass - finite difference PA /PP c P source A sink: L sink L/F source 1,3 c A source P sink: L source L/F sink 1,2 xmpi=0.0 if(ifk .eq. 1) xmpi=0.111 if(ifk .eq. 3) xmpi=0.106 if(ifk .eq. 323) xmpi=0.111 if(ifk .eq. 3231) xmpi=0.111 if(ifk .eq. 3232) xmpi=0.111 if(ifk .eq. 326) xmpi=0.144 if(ifk .eq. 246) xmpi=0.144 if(ifk .eq. 206) xmpi=0.144 if(ifk .eq. 328) xmpi=0.164 if(ifk .eq. 3212) xmpi=0.201 if(ifk .eq. 85) xmpi=0.226 if(ifk .eq. 851) xmpi=0.226 if(ifk .eq. 852) xmpi=0.226 if(ifk .eq. 4) xmpi=0.1358 if(ifk .eq. 432) xmpi=0.136 if(ifk .eq. 332) xmpi=0.115 if(ifk .eq. 41) xmpi=0.295 if(ifk .eq. 64) xmpi=0.1712 if(ifk .eq. 8) xmpi=0.194 if(ifk .eq. 81) xmpi=0.194 if(ifk .eq. 10) xmpi=0.211 if(ifk .eq. 15) xmpi=0.259 if(ifk .eq. 45) xmpi=0.162 if(ifk .eq. 60) xmpi=0.222 if(ifk .eq. 61) xmpi=0.184 if(ifk .eq. 6) xmpi=0.181 if(ifk .eq. 624) xmpi=0.181 if(ifk .eq. 824) xmpi=0.21 if(ifk .eq. 11) xmpi=0.242 if(ifk .eq. 165) xmpi=0.296 if(ifk .eq. 38) xmpi=0.184 if(ifk .eq. 6661) xmpi=1.12 if(ifk .eq. 666) xmpi=0.20 if(ifk .eq. 444) xmpi=0.11 if(ifk .eq. 446) xmpi=0.17 if(ifk .eq. 456) xmpi=0.17 if(ifk .eq. 448) xmpi=0.27 if(ifk .eq. 9) xmpi=0.222 if(ifk .eq. 12) xmpi=0.254 if(ifk .eq. 153) xmpi=0.283 if(ifk .eq. 26) xmpi=0.260 if(ifk .eq. 77) xmpi=0.2445 if(ifk .eq. 66) xmpi=0.18 if(ifk .eq. 776) xmpi=0.174 if(ifk .eq. 775) xmpi=0.174 if(ifk .eq. 27) xmpi=1.3 if(ifk .eq. 482) xmpi=0.073 if(ifk .eq. 2112) xmpi=0.073 if(ifk .eq. 455) xmpi=0.155 if(ifk .eq. 435) xmpi=0.124 if( xmpi .lt. 0.0001 ) STOP "xmpi not defined" c iob=1 LL 2 LF 3 FL 4 FF Do iob=1,3 Do idt=2,23 idtp=idt+1 idtm=idt-1 Do igg=1,ngauge pcac(igg)=(-c(igg,idtp,itp2,iob)+c(igg,idtm,itp2,iob))*0.5*0.5 Enddo call avratb(pcac(1),c(1,idt,1,iob),ngauge,a,s,nblock) c m_pcac assuming exponentials as pion mass PA /PP (Less noisy) call avratb(c(1,idt,itp2,iob),c(1,idt,1,iob),ngauge,an,sn,nblock) a1=an*xmpi*0.5 s1=sn*xmpi*0.5 c 22-6-06 correct for sinh/cosh?? - but d/dt of sinh is cosh - so no effect c exprat=exp(-xmpi*(LX4-2*idt)) c tanch=(1.0-exprat)/(1.0+exprat) c a1=a1/tanch c s1=s1/tanch t=float(idt)+0.1*(iob-1) write(30,9397)t,a1,s1 write(31,9397)t,a,s 9397 format(3f15.5) if(iprin .gt. 1)write(6,9395) iob, idt, a,s,a1,s1 c9395 format(2i5,2f15.5) Enddo Enddo c m PCAC versus gauge number (pick t or could average here) idtv=10 c idtv=21 iob=1 c *** change to PA/ ave_of_PP call avn(c(1,idtv,1,iob),ppav,eppav,ngauge) Do igg=1,ngauge c xmpcac=xmpi*0.5*c(igg,idtv,itp2,iob)/c(igg,idtv,1,iob) xmpcac=xmpi*0.5*c(igg,idtv,itp2,iob)/ppav write(31,9494) igaugest(igg), xmpcac 9494 format(i5,f15.5) Enddo c STOP c other mixing ratios - c rho AA is 9+2 VV is 9+3 AV is 9+8 VA is 9+9 itp1=2+9 itp2=8+9 if(iprin .gt. 1) write(*,*) itp1,itp2 Do iob=1,4 Do idt=1,23 call avrat(c(1,idt,itp2,iob),c(1,idt,itp1,iob),nlst,ngauge,a,s) if(iprin .gt. 1)write(6,9395) iob, idt, a,s c9395 format(2i5,4f15.5) Enddo Enddo c STOP c sum over pion to test GMO type identities (check ratio to exp(-m_pi t) 0 c compare with gd state pion piece only,, and GMO-type sum if (ifk .eq. 4) then itype=1 iobs=1 idt=lx4h call avn(cll(1,idt,itype,iobs),cav1,sav1,ngauge) sum=cav1 lx4hm=lx4h-1 c *** fix here by hand from fit **** fpi=0.0653 xnn=(0.5/xmpi)*(fpi*xmpi**2/(4.0*coupK*g5mu))**2 rat=xnn*( exp(-xmpi*idt)+exp(-xmpi*(lx4-idt)) ) write(6,*)' pion sum to t=',idt, sum, rat do idt=lx4hm,1,-1 ratc=2.0*xnn*( exp(-xmpi*idt)+exp(-xmpi*(lx4-idt)) ) rat=rat+ratc call avn(cll(1,idt,itype,iobs),cav1,sav1,ngauge) call avn(cll(1,idt,itype,iobs+1),cav2,sav2,ngauge) sum=sum+cav1+cav2 write(6,*)' pion sum to t=',idt, sum, rat enddo idt=0 call avn(cll(1,idt,itype,iobs),cav1,sav1,ngauge) sum=sum+cav1 ratc=xnn* exp(-xmpi*idt) rat=rat+ratc write(6,*)' pion sum to t=',idt, sum, rat ratg=xnn*2.0/xmpi write(6,*)' GMO method:',ratg endif c STOP c BLOCKING OPTION c combine every istep to every istep*nblock if needed c 004 24^3 could block to same number of trajs (64) ie 64 / 6 c ie pre-block the 920-3398 every 4 c NOW ALL random t and consistent - so not used so 4 -> 4111 nblock=$NBLOCK nblock2=0 if( ifk .eq. 4111) nblock2=6 write(*,*)' nblock=',nblock,nblock2 if(nblock .gt. 1) then nobs=nobsv/2 ngg2=ngauge/nblock if(nblock2 .gt. 0) ngg2=ngauge/nblock2 iblock=0 igg1=1-nblock c seems to be bug here pre 31-1-07 - but irrel since wrote beyond ngg c do igg=1,ngauge do igg=1,ngg2 if(ifk. eq. 4111 .and. igaugest(igg1+nblock) .gt. 3398-nblock) & nblock=nblock2 igg1=igg1+nblock xnb=1.0D0/float(nblock) if( igg1+nblock .le. ngauge) then iblock=iblock+1 write(*,*) igg,igg1,iblock,igaugest(igg1) do iobs=1,nobs do itype=1,ntype do idt=0,lx4h c copies to smaller igg value so overwrites OK csum=0.0 do ibl=1,nblock csum=csum+c(igg1+ibl,idt,itype,iobs) enddo c(igg,idt,itype,iobs)=csum*xnb Enddo Enddo Enddo Endif Enddo c ngauge=ngg ngauge=iblock endif c set up PLAV and do simple effective masses (BOOTV) ir=1 nconf=ngauge nconfa=nconf itmax=lx4h CHOOSE Quantum numbers here - and path choices c now just gamma x gamma (LL LF FL FF elsewhere) c create symmetric iob.itype Do igg=1,nconf Do idt=0,lx4h Do iob1=1,2 Do iob2=1,2 c order written (by mallprc.f) is LL LF FL FF (source-sink) c exchange order 27-4-06 so 1 2 3 4 is LL LF FL FF c iob=(iob2-1)*2+iob1 iob=(iob1-1)*2+iob2 npv2=9 if( iobval .ge. 3) npv2=3 Do iop=1,npv2 iop1=1 iop2=1 nfacp=2 npv=3 if(iobval .eq. 1) then itype=iop iop1=ipion1(iop) iop2=ipion2(iop) nfacp=6 npv=21 elseif(iobval .eq. 2) then itype=iop+9 iop1=irho1(iop) iop2=irho2(iop) nfacp=6 npv=21 elseif(iobval .eq. 3) then c ****before 8-1-08 this was wrong way round****** 3 is a0 4 is b1 c itype=20 c so swap #19 --- a0 --- 3 #20 ---b1 ---4 itype=19 elseif(iobval .eq. 4) then c itype=19 itype=20 else STOP Endif corl(igg,idt,iob1,iop1,iob2,iop2)= c(igg,idt,itype,iob) & *sdv(itype) c carsten's output (v1) has 2 22 19 20 opp sign to CM if(ifk .eq. 44400)then if(iobval .eq. 1 .and. iop2 .eq. 2) & corl(igg,idt,iob1,iop1,iob2,iop2)= & - corl(igg,idt,iob1,iop1,iob2,iop2) endif Enddo Enddo Enddo Enddo Enddo c diagnostic test of pion factorisation: idt=9 write(6,*) 'pion fact test idt=',idt do iop1=1,3 do iop2=1,3 do iob1=1,2 do iob2=1,2 call avn(corl(1,idt,iob1,iop1,iob1,iop1),ca11,s,ngauge) call avn(corl(1,idt,iob2,iop2,iob2,iop2),ca22,s,ngauge) call avn(corl(1,idt,iob1,iop1,iob2,iop2),ca12,s,ngauge) r=ca12/sqrt(abs(ca11*ca22)) write(6,9393) iob1,iob2,iop1,iop2, r,ca11,ca22,ca12 9393 format('LF',2i2,'PA4',2i2,4f15.5) Enddo Enddo Enddo Enddo c symmetrize / pick which is less noisy & set up plav do 400 ip=1,npv ip1=ip1v(ip) ip2=ip2v(ip) ip1s(ip)=ip1 ip2s(ip)=ip2 nop=nfacp/2 iop1=(ip1-1)/2 +1 iop2=(ip2-1)/2 +1 iob1=ip1-(iop1-1)*2 iob2=ip2-(iop2-1)*2 c COSH fits, some SINH ifs(ip)=0 c need sinh when sf=-1 12 13 not 23 CHANGE from all off-diag is sinh 11-06 if(iop1 .ne. iop2 .and. iop1 .eq. 1) ifs(ip)=1 write(6,905) ip,iob1,iop1,iob2,iop2,nfacp 905 format(' path number ',i2,' obs1=',2i3,' obs2=',2i3, & ' nfacp=',i5) nconf=nconfa do icv=1,nconfa ic=icv do jt=0,itmax if ( ifk .eq. 27) then c ( favour LF rel FL etc here - but symmetrize on ip ?) ct(ic,jt)=corl(ic,jt,iob1,iop1,iob2,iop2) if(iob2 .gt. iob1) & ct(ic,jt)=corl(ic,jt,iob2,iop1,iob1,iop2) else c symmetrize ct(ic,jt)=(corl(ic,jt,iob1,iop1,iob2,iop2)+ & corl(ic,jt,iob2,iop2,iob1,iop1))*0.5 endif c best to choose smallest at sink for off-diag (order source-sink) c swap 28-4-06 to this below for pion c for pion order P g4 a 1 3 2 keep order for sign c CARE 3,4 . 5,6 can be -ve - so option to swap by hand xswap=1.0 c seems swap preferred ifk=0,1,2,3 - not obvious why c ifk=0 6-18 pi 6x6 xswap=-1 chi2=257 xswap=1 chi2=669 (dof 273-7) c ifk=4 6-18 pi 6x6 xswap=-1 chi2=411 xswap=1 chi2=400 (dof 273-7) nbl=8 c ifk=4 variational sick unless xswap=-1 c 11-06 fix 44 t symmetry - then xswap=1 is OK (also checked rho) c set iobval to 1 to swap, else 11 (irrel for nfacp < 6 fits) if(iobval .eq. 11 .and. iop1 .eq. 2 .and. iop2 .eq. 3) & xswap=-1.0 if(iobval .eq. 1 .and. iop1 .eq. 1 .and. iop2 .eq. 3) & ct(ic,jt)=corl(ic,jt,iob1,iop1,iob2,iop2) c A4 vs 4A not much different in error ifk=60 if(iobval .eq. 1 .and. iop1 .eq. 2 .and. iop2 .eq. 3) & ct(ic,jt)=corl(ic,jt,iob2,iop2,iob1,iop1)*xswap if(iobval .eq. 1 .and. iop1 .eq. 1 .and. iop2 .eq. 2) & ct(ic,jt)=corl(ic,jt,iob1,iop1,iob2,iop2) c test PA -> AP swap c & ct(ic,jt)=corl(ic,jt,iob2,iop2,iob1,iop1) c for rho - a1 order rho_4 a_1 rho is 1 3 2 also c find error slightly smaller 14 16 17 so could use them c but differences are not so big - so could just average (ie set 2 to 92 if(iobval .eq. 92 .and. iop1 .eq. 1 .and. iop2 .eq. 3) & ct(ic,jt)=corl(ic,jt,iob1,iop1,iob2,iop2) if(iobval .eq. 92 .and. iop1 .eq. 2 .and. iop2 .eq. 3) & ct(ic,jt)=corl(ic,jt,iob2,iop2,iob1,iop1) if(iobval .eq. 92 .and. iop1 .eq. 1 .and. iop2 .eq. 2) & ct(ic,jt)=corl(ic,jt,iob1,iop1,iob2,iop2) plav(ip1,ip2,jt+1,ic,ir)=ct(ic,jt) c Bugfix 27-11-97 plav 12 not 21 defined - put symmetric values here (TMAT) if(ip1 .ne. ip2) & plav(ip2,ip1,jt+1,ic,ir)=ct(ic,jt) c test output option 4 (pion) is PA LL if(ip.eq.4 .and. jt .eq. 10) write(6,9209) ic,ct(ic,jt),0.0 9209 format(i5,2f15.7) enddo enddo c c For small nconf JackKnife is more stable than Bootstrap (used for graph) c c call bootv(ct(1,0),nlst,nconf,itmax,ip) call jackv(ct(1,0),nlst,nconf,itmax,ip) Do jt=0,lx4h call avn(ct(1,jt),cav,sav,nconf) c convert to stoch convention c for rho PVSAB is sum of 3 polarisations - so /3 c assume TMQCD has M=1+K.D +2K mu ... used to eval corr c need M=(1/2k) +D/2 + mu ... so corrc= (2k)^2 corr if(iobval .eq. 2) xkf=0.3333333333/(coupk*2.0)**2 if(iobval .eq. 4) xkf=0.3333333333/(coupk*2.0)**2 if(iobval .eq. 1) xkf=1.0/(coupk*2.0)**2 if(iobval .eq. 3) xkf=1.0/(coupk*2.0)**2 Care - choose convention needed here xkf=1.0 c choose rho as tmzc0 if(iobval .eq. 2) xkf=0.333333333 cav=cav*xkf sav=sav*xkf/sqrt(float(nconf-1)) write(6,939) jt,ip,cav,sav c for g5 g5 pion if( ip .le. 3) write(69,939) jt,ip,cav,sav c will also write out a0 if iobval=3 to fort.69 c for gig5 gig5 rho (current-current) if( ip1 .ge. 5 .and. ip2 .ge. 5) write(70,939) jt,ip,cav,sav c for g4.g4 pion LL, LF, FF if( ip .eq. 15) write(68,939) jt,1,cav,sav if( ip .eq. 20) write(68,939) jt,2,cav,sav if( ip .eq. 21) write(68,939) jt,3,cav,sav Enddo 939 format(2i6,2f17.6) c dump test if(ip.eq.3 .and. iprin .gt. 20)then write(6,909)(ct(ic,0),ic=1,nconf) 909 format(5f10.5) endif 400 continue c evaluate ratios of parity violating amplitudes if ( iobval .eq. 1 ) then c pi-X 4A/44 order P A 4 ip2=3 ip1=5 elseif( iobval .eq. 2 ) then c rho-a1 AV/AA order 4 V A ip2=3 ip1=5 endif if(iprin .gt. 1) write(*,*) ' mixing of parity odd', ip1,ip2 Do idt=1,lx4h Do ic=1,ngauge ct(ic,1) = plav(ip1,ip1,idt+1,ic,ir) ct(ic,2) = plav(ip1,ip2,idt+1,ic,ir) Enddo tidt=idt+0.1 call avrat(ct(1,2),ct(1,1),nlst,ngauge,a,s) if(iprin .gt. 1)write(6,9315) tidt, a,s 9315 format(f5.1,4f15.5) Enddo c STOP c Variational analysis here c stuff to pass to TMAT via /INFO/ NLAT=ngauge NLAT=nconf C NLAT BLOCKS IN ITBIAS-1/ITBIAS BASIS itmax=lx4h ITMAXP=ITMAX+1 c itvmax NRMN NR not used now NRMN=1 NR=1 ITVMAX=3 ITBIAS=2 ITBIAS=3 ITBIAS=4 IPRINv=3 NPVV=nfacp call TMAT(ir) call BOOTV2(ir) c CM's fancy fitting code now npv paths from nfacp.nfacp matrix ir=1 npv=$NPV nfacp=$NFACP c call with average vaue if g5mu=/= g5mu1 if(ifk .eq. 26) g5mu=(g5mu+g5mu1)*0.5 call tfit(ir,it1,it2,g5mu,coupK) STOP END C--------------------------------------------------------- C from tmpr 14-12-95 **here exp assumed*** - Now cosh - not double checked c 8-6-00 now cosh/sinh C CALCULATES EFF POTENTIAL = mass C METHOD IS NBS BOOTSTRAP CALLS c (set ir for xeff, ip not used) C--------------------------------------------------------- SUBROUTINE BOOTV(ad,nconfd,nconf,itmax,ip) IMPLICIT double precision (A-H,O-Z) include 'fitex.h' PARAMETER (itst=itsth,npstx=npsth,nrst=nrsth,LTHP=lx4val/2) PARAMETER (NBS=99,lx4=lx4val,lx4h=lx4/2,lx4hp=lx4h+1) COMMON/mass/xeff(npstx,itst,nrst),xerr(npstx,itst,nrst) common/fact/npv,ip1(npsth),ip2(npsth),ifs(npsth), nfacp DIMENSION ad(nconfd,0:itmax),XMAV(0:ITST),A(0:ITST) & ,XM(NBS,0:ITST),XMT(NBS,0:ITST) C----CALCULATE AVERAGES BEFORE BOOTSTRAP LOOP ir=1 iprin=12 c iprin=2 nlat=nconf itmin=0 itmaxm=itmax-1 DO 10 IT=ITMIN,ITMAXm AA=0.0 BB=0.0 DO 40 NL=1,NLAT AA=AA+ad(nl,it) BB=BB+ad(nl,it+1) 40 CONTINUE c xmav(it)=-log(bb/aa) cosh xmav(it)=cfs(ifs(ip),LTHP-it,bb/aa) xeff(ip,it+1,ir)=xmav(it) 10 CONTINUE c write(6,909) (xmav(it),it=0,itmaxm) 909 format(6g15.5) c ---- START BOOT STRAP LOOP ------- AN=FLOAT(NLAT) XNBS=1.0/FLOAT(NBS) C C PICK NBS SETS OF NLAT RANDOMLY FROM NLAT ENSEMBLES C (WITH REPLACEMENT) C-BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB DO 80 KSE =1,NBS DO 21 IT=itmin,ITMAX A(IT)=0.0 21 CONTINUE DO 410 INL=1,NLAT NL = INT(RNDM2(DUMMY) * AN ) + 1 c write(6,959)nl 959 format(i20) DO 210 IT=itmin,ITMAX A(IT)=A(IT)+ad(nl,it) 210 CONTINUE 410 CONTINUE DO 221 IT=itmin,ITMAXm c XM(KSE,IT)=-log(A(it+1)/A(it)) xm(kse,it)=cfs(ifs(ip),LTHP-it,A(it+1)/A(it)) 221 CONTINUE C----------------------------------------------------------- DO 222 IT=itmin+1,ITMAXm XMT(KSE,IT)=-XM(KSE,IT)+XM(KSE,IT-1) 222 CONTINUE C 22 CONTINUE 80 CONTINUE C-BBBBBBBBBBBBB END BOOTSTRAP LOOP BBBBBBBBBBBBBBBBBBBBBB- 988 FORMAT(6G12.4) C C --SORT OUT POTENTIALS , ASYMPTOTIC CORRECTIONS, ETC C CALL AVN(XM(1,itmin),AV,S,nbs) xerr(ip,itmin+1,ir)=s IF(IPRIN.GT.1) & WRITE(6,97) itmin+1,AV,XMAV(itmin),S DO 24 IT=itmin+1,ITMAXm CALL AVN(XM(1,IT),AV,S,nbs) xerr(ip,it+1,ir)=s CALL AVN(XMT(1,IT),AVT,ST,nbs) XMAVT=-XMAV(IT)+XMAV(IT-1) IF(IPRIN.GT.1) & WRITE(6,97) IT+1,AV,XMAV(IT),S,XMAVT,ST tv=float(it)+0.5 c for m_eff plots rho gi g4 g5 set iprin=12 above IF(IPRIN.GT.11 .and.(ip.eq.15.or.ip.eq.20.or.ip.eq.21)) & WRITE(11,9777) tv,XMAV(IT),S 9777 format(3f15.5) 24 CONTINUE 241 CONTINUE 97 FORMAT(' TMAX=',I2,2F9.5,'+/-',5F9.5) 2001 continue RETURN END C--------------------------------------------------------- c **here exp as is *** C CALCULATES EFF POTENTIAL = mass c mod BOOTV 29-2-96 to do jackknife c (set ir for xeff, ip not used) C--------------------------------------------------------- SUBROUTINE JACKV(ad,nconfd,nconf,itmax,ip) IMPLICIT DOUBLE PRECISION (A-H,O-Z) include 'fitex.h' PARAMETER (itst=itsth,npstx=npsth,nrst=nrsth,LTHP=lx4val/2) PARAMETER (NBS=nlst,lx4=lx4val,lx4h=lx4/2,lx4hp=lx4h+1) COMMON/mass/xeff(npstx,itst,nrst),xerr(npstx,itst,nrst) common/fact/npv,ip1(npsth),ip2(npsth),ifs(npsth), nfacp DIMENSION ad(nconfd,0:itmax),XMAV(0:ITST),A(0:ITST) & ,XM(NBS,0:ITST),XMT(NBS,0:ITST) C----CALCULATE AVERAGES BEFORE JACKNIFE LOOP nlat=nconf xnlat=nlat xnlatm=nlat-1 xnlatr=xnlatm/xnlat ir=1 cc set > 11 to have for.11 output for meff graph iprin=2 c iprin=12 itmin=0 itmaxm=itmax-1 DO 10 IT=ITMIN,ITMAXm AA=0.0 BB=0.0 DO 40 NL=1,NLAT AA=AA+ad(nl,it) BB=BB+ad(nl,it+1) 40 CONTINUE c xmav(it)=-log(bb/aa) cosh xmav(it)=cfs(ifs(ip),LTHP-it,bb/aa) xeff(ip,it+1,ir)=xmav(it) 10 CONTINUE c write(6,909) (xmav(it),it=0,itmaxm) 909 format(6g15.5) AN=FLOAT(NLAT) XNBS=1.0/FLOAT(nlat) C C-Single elimination Jack-Knife C C----------------------------------------------------------- DO 80 KSE =1,nlat DO 21 IT=itmin,ITMAX A(IT)=0.0 21 CONTINUE DO 410 NL=1,NLAT if(nl.ne.kse)then DO 210 IT=itmin,ITMAX A(IT)=A(IT)+ad(nl,it) 210 CONTINUE endif 410 CONTINUE DO 221 IT=itmin,ITMAXm c XM(KSE,IT)=-log(A(it+1)/A(it)) xm(kse,it)=cfs(ifs(ip),LTHP-it,A(it+1)/A(it)) 221 CONTINUE DO 222 IT=itmin+1,ITMAXm XMT(KSE,IT)=xnlat*(-XMav(IT)+XMav(IT-1)) & -xnlatm*(-XM(KSE,IT)+XM(KSE,IT-1)) 222 CONTINUE DO 224 IT=itmin,ITMAXm XM(KSE,IT)=xnlat*xmav(it)-xnlatm*xm(kse,it) 224 CONTINUE C 22 CONTINUE 80 CONTINUE C-------------------------------------------------------- 988 FORMAT(6G12.4) C sff=1.0/sqrt(float(nlat-1)) CALL AVN(XM(1,itmin),AV,S,nlat) s=s*sff xerr(ip,itmin+1,ir)=s IF(IPRIN.GT.1) & WRITE(6,97) itmin+1,AV,XMAV(itmin),S DO 24 IT=itmin+1,ITMAXm CALL AVN(XM(1,IT),AV,S,nlat) s=s*sff xerr(ip,it+1,ir)=s CALL AVN(XMT(1,IT),AVT,ST,nlat) st=st*sff XMAVT=-XMAV(IT)+XMAV(IT-1) IF(IPRIN.GT.1) & WRITE(6,97) IT+1,AV,XMAV(IT),S,XMAVT,ST tv=float(it)+0.5 c for m_eff plots rho gi g4 g5 set iprin=12 above IF(IPRIN.GT.11 .and.(ip.eq.15.or.ip.eq.20.or.ip.eq.21)) & WRITE(11,9777) tv,XMAV(IT),S 9777 format(3f15.5) 24 CONTINUE 241 CONTINUE 97 FORMAT(' TMAX=',I2,2F9.5,'+/-',5F9.5) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TAKE PLAV(IP,IP,IT,NL,IR) TO W(NLAT,ITMAX,NR) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TMAT(IC) IMPLICIT REAL*8 (A-H,O-Z) include 'fitex.h' PARAMETER (itst=itsth,npst=npstph,nrst=nrsth) COMMON/INFO/NLAT,NR,NRMN,NPVI, ITMAX,ITMAXP,ITVMAX,ITBIAS,IPRIN COMMON/WLOOP/W(NLST,ITST,NRST) COMMON/EPP/PLAV(NPST,NPST,ITST,NLST,NRST), NCONFV COMMON/TJ/TM(NPST,NPST,ITST,NLST), 1 TOTAL(NPST,NPST,ITST),EVALS(NPST),EVECTS(NPST,NPST),EEVAL(NPST) REAL*8 WKSPCE(NPST),AAD(NPST,NPST),BBD(NPST,NPST),ADI(NPST,NPST) REAL*8 BDI(NPST,NPST),CDI(NPST,NPST) Integer ipu(npst) c could reorder here *** match size to npstph *****npv can be smaller Data ipu /1,2,3,4,5,6/ npv=npvi c npv=2 write(6,90) (ipu(j),j=1,npv) 90 format(' paths used=',10i4) write(*,*)' basis used ',itbias,'/',itbias-1 C SETS EIGENVECTORS TO IDENTITY: USE IF DIAGONAL ONLY WANTED DO 30 I=1,NPV DO 30 J=1,NPV EVECTS(I,J)=0.0 IF(I.EQ.J)EVECTS(I,J)=1.0 DO 30 K=1,ITMAXP 30 TOTAL(I,J,K)=0.0 *READ IN, SYMMETRISE AND TOTALISE TM'S XNV=1.0/NLAT XNV=xnv*1000.0 DO 40 NL=1,NLAT DO 41 I=1,NPV DO 41 J=1,NPV Iu=ipu(i) Ju=ipu(j) DO 41 K=1,ITMAXP TM(I,J,K,NL)=(PLAV(Iu,Ju,K,NL,IC)+PLAV(Ju,Iu,K,NL,IC))*XNV*0.5 TOTAL(I,J,K)=TOTAL(I,J,K)+TM(I,J,K,NL) 41 CONTINUE 40 CONTINUE iprin=3 if(iprin .gt. 3) then do 173 k=1,itmaxp write(6,976)k 976 format(' it=',i2) do 174 i=1,npv write(6,977) (total(i,j,k),j=1,npv) 977 format(7g11.4) 174 continue 173 continue endif EV=0.0 EVL=0.0 EVLERR=0.0 C <------ FIX BASIS HERE KK=ITBIAS *CALCULATE INVERSE OF MATRIX TOTALISED DO 66 JD=1,NPV DO 67 ID=1,NPV BDI(ID,JD)=0.D0 IF(JD.EQ.ID) BDI(ID,JD)=1.D0 ADI(ID,JD)=TOTAL(ID,JD,KK) 67 CONTINUE 66 CONTINUE IFAIL=0 CALL F04AEF(ADI,NPST,BDI,NPST,NPV,NPV,CDI,NPST,WKSPCE, 1 AAD,NPST,BBD,NPST,IFAIL) DO 68 ID=1,NPV DO 69 JD=1,NPV TOTAL(ID,JD,KK)=CDI(ID,JD) 69 CONTINUE 68 CONTINUE *CALCULATE EIGENVECTORS OF INVTOTAL*TOTAL CALL EIGENT(KK,KK+1,NPV,IC) RETURN END C*********************************************************************** CALCULATE EIGENVECTORS OF INVTOTAL*TOTAL C*********************************************************************** SUBROUTINE EIGENT(KK,KL,NPV,IC) IMPLICIT REAL*8 (A-H,O-Z) include 'fitex.h' PARAMETER (itst=itsth,npst=npstph,nrst=nrsth,lwork=npst*npst) PARAMETER (nparms=npst*(npst+1)) COMMON/INFO/NLAT,NR,NRMN,NP, ITMAX,ITMAXP,ITVMAX,ITBIAS,IPRIN COMMON/WLOOP/W(NLST,ITST,NRST) Common/start/x(nparms) COMMON/EXCIT/XL(NRST), V0(NRST) COMMON/TJ/TM(NPST,NPST,ITST,NLST), 1 TOTAL(NPST,NPST,ITST),EVALS(NPST),EVECTS(NPST,NPST),EEVAL(NPST) DIMENSION EVECTSI(NPST,NPST) INTEGER INTGER(NPST) DIMENSION TEMP(NPST,NPST),EVALSI(NPST) Dimension IPIV(npst), work(lwork) DO 10 I=1,NPV DO 10 J=1,NPV TEMP(I,J)=0.0 DO 10 K=1,NPV 10 TEMP(I,J)=TEMP(I,J)+TOTAL(I,K,KK)*TOTAL(K,J,KL) IFAIL=0 CALL F02AGF(TEMP,NPST,NPV,EVALS,EVALSI,EVECTS,NPST,EVECTSI, 1 NPST, INTGER,IFAIL) *...AND ORDER IN DECREASING SIZE OF EIGENVALUE 20 NCH=0 IF(NPV.GT.1) THEN DO 40 I=1,NPV-1 IF(EVALS(I).GE.EVALS(I+1))GOTO 40 NCH=1 STORE=EVALS(I) EVALS(I)=EVALS(I+1) EVALS(I+1)=STORE DO 30 K=1,NPV STORE=EVECTS(K,I) EVECTS(K,I)=EVECTS(K,I+1) 30 EVECTS(K,I+1)=STORE 40 CONTINUE IF(NCH.EQ.1) GOTO 20 ENDIF C <---BIGGEST E VALUE EVECTS(IP,NE) XM1=-LOG(EVALS(1)) IF(NPV.GT.1) XM2=-LOG(EVALS(2)) XL(IC)=EVALS(2)/EVALS(1) V0(IC)=XM1 IF(IPRIN.GT.1)WRITE(6,98)IC,XM1,XM2,XL(IC) do iev=1,npv xm=-log(evals(iev)) if(iprin.gt.2) write(6,98) iev, xm IF(IPRIN.GT.1)WRITE(6,988)iev,(EVECTS(IP,iev),IP=1,NPV) enddo 98 FORMAT(I6,3F12.6) c 988 FORMAT(6F12.6) 988 FORMAT(I5,6F12.6) c code here from tmball in ../bm if(iprin .gt. 1) write(*,*) ' inverse' c Au = L B u thus B^{-1} A u = L u and uAu = L uBu c now for fitting B=vv A=vLv (both summed on evect) c which implies vu=1 so solve using v=u^{-1} c norm: take u-> u/c v-> vc so uBu=c^2 B=c^2 vv c invert evects for use in fitting c see www.nag.co.uk do i=1,npv do j=1,npv temp(i,j)=evects(i,j) enddo enddo call f07adf(npv,npv,temp,npst,ipiv,info) if(info .ne. 0) write(*,*) 'f07adf info=',info call f07ajf(npv,temp,npst,ipiv,work,lwork,info) c if(info .ne. 0) write(*,*) 'f07ajf info=',info write(*,*) 'v coeff eig_v v ; ip->' Do iv=1,npv IF(IPRIN.GT.1)WRITE(6,988) iv,(TEMP(iv,ip),IP=1,NPV) Enddo c check Do iv1=1,npv Do iv2=1,npv d=0 Do i=1,npv d=d+temp(iv1,i)*evects(i,iv2) Enddo if( abs(d) .gt. 0.00001 .and. iv1 .ne. iv2) & write(*,*)iv1,iv2,d Enddo Enddo c normalise v for fitting to 11 at t=KK b=0.0 Do nl=1,nlat b=b+tm(1,1,kk,nl) Enddo xn=b/float(nlat) Do iv=1,npv b=0.0 Do nl=1,nlat Do i=1,npv Do j=1,npv b=b+evects(i,iv)*tm(i,j,kk,nl)*evects(j,iv) Enddo Enddo Enddo b=b/float(nlat) if(b.gt.0) then xc=sqrt(b/xn) else write(*,*)' negative norm for iv=',iv xc=0.0 c stop endif Do i=1,npv temp(iv,i)=xc*temp(iv,i) Enddo Enddo c set up starting values for fit (temp(iv,i) evals(i) Do iv=1,npv Do i=1,npv ivv=(iv-1)*npv+i x(ivv)=temp(iv,i) Enddo ivv=npv*npv+iv x(ivv)=-log(evals(iv)) Enddo write(*,*) ' x coeffs for fit ip -> ' Do iv=1,npv IF(IPRIN.GT.1)WRITE(6,988) iv,(TEMP(iv,ip),IP=1,NPV) Enddo C -- FOLLOWING FROM ERROR -------- write(*,*)' u coeff ip ->' Do NE=1,2 c NE=1 IF(IPRIN.GT.1)WRITE(6,988)NE,(EVECTS(IP,NE),IP=1,NPV) DO 33 KN=1,ITMAXP C: <----- LOOP OVER LATTICES USED DO 33 NL=1,NLAT *CALCULATE A= (EVECT)(EVECT) A=0.0 DO 23 I=1,NPV STORE=0.0 DO 13 K=1,NPV 13 STORE=STORE+TM(I,K,KN,NL)*EVECTS(K,NE) A=A+EVECTS(I,NE)*STORE 23 CONTINUE W(NL,KN,NE)=A 33 CONTINUE Enddo RETURN END C--------------------------------------------------------- C Assumes estimate of e0,e1 available (XLAM from XL from TMAT) C 21-1-92 FIX CORRECTION NEGATIVE ==> ZERO C 23-7-92 SUMMARY: C INPUT WLOOP (FROM ITBIAS VARIATIONAL BEST EIGENVECTOR) C XL LAM1/LAM0 FROM VARIATIONAL C CALCULATES EFF POTENTIAL = mass C METHOD IS NBS BOOTSTRAP CALLS -- here t 1,...T/2+1 C--------------------------------------------------------- SUBROUTINE BOOTV2(IR) IMPLICIT REAL*8 (A-H,O-Z) include 'fitex.h' PARAMETER (itst=itsth,nrst=nrsth,LTHP=lx4val/2) PARAMETER (NBS=99) COMMON/INFO/NLAT,NR,NRMN,NP, ITMAX,ITMAXP,ITVMAX,ITBIAS,IPRIN COMMON/WLOOP/W(NLST,ITST,NRST) COMMON/EXCIT/XL(NRST),V0(NRST) DIMENSION XMAV(ITST,NRST),A(ITST,NRST) & ,XM(NBS,ITST,NRST),XMT(NBS,ITST,NRST),XMTC(NBS,ITST,NRST) C----CALCULATE AVERAGES BEFORE BOOTSTRAP LOOP c 19-12-05 fix up to do irmax eigenvectors c could pass iobval in here irmax=1 if( $IOBVAL .eq. 2) irmax=2 if( $IOBVAL .eq. 1) irmax=2 Do ir=1,irmax DO 10 IT=1,ITMAX AA=0.0 BB=0.0 DO 40 NL=1,NLAT AA=AA+W(NL,IT,IR) BB=BB+W(NL,IT+1,IR) 40 CONTINUE c XMAV(IT,IR)=-LOG(BB/AA) cosh xmav(it,ir)=cfs(0,LTHP-it+1,bb/aa) 10 CONTINUE 20 CONTINUE C ---- START BOOT STRAP LOOP ------- AN=FLOAT(NLAT) XNBS=1.0/FLOAT(NBS) C C PICK NBS SETS OF NLAT RANDOMLY FROM NLAT ENSEMBLES C (WITH REPLACEMENT) C-BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB DO 80 KSE =1,NBS DO 21 IT=1,ITMAXP A(IT,IR)=0.0 21 CONTINUE DO 410 INL=1,NLAT NL = INT(RNDM2(DUMMY) * AN ) + 1 DO 210 IT=1,ITMAXP A(IT,IR)=A(IT,IR)+W(NL,IT,IR) 210 CONTINUE 410 CONTINUE DO 221 IT=1,ITMAX c XM(KSE,IT,IR)=-LOG(A(IT+1,IR)/A(IT,IR)) xm(kse,it,ir)=cfs(0,LTHP-it+1,A(it+1,ir)/A(it,ir)) 221 CONTINUE XLAM=XL(IR) C----------------------------------------------------------- C C ATTEMPT TO EXTRAPOLATE EFF-ENERGY IN T: C LOG(W(T))=-T*M+LOG(1.0 + A*XL**T) XL=EVAL(2)/EVAL(1) C-LOG(W(T+1)/W(T))=M-LOG(1.0 + A*XL*XL**T)+LOG(1.0*A*XL**T) C IE XEFF(T) =M + A*XL**T *(1-XL) C SO XEFF(T-1)-XEFF(T)=A*XL**(T-1)*(1-XL)*(1-XL) C AND M=XEFF(T) - XL*( XEFF(T-1)-XEFF(T) )/ (1-XL) C MOD 21-1-91 REQUIRE MONOTONIC DECREASE OF V(T) SO CORRECTION C ZERO OTHERWISE C----------------------------------------------------------- DO 222 IT=2,ITMAX XMT(KSE,IT,IR)=-XM(KSE,IT,IR)+XM(KSE,IT-1,IR) XCOR=XLAM*XMT(KSE,IT,IR)/(1.-XLAM) IF(XCOR.LT.0.0)XCOR=0.0 XMTC(KSE,IT,IR)= XM(KSE,IT,IR)-XCOR 222 CONTINUE C 22 CONTINUE 80 CONTINUE C-BBBBBBBBBBBBB END BOOTSTRAP LOOP BBBBBBBBBBBBBBBBBBBBBB- 988 FORMAT(6G12.4) C C --SORT OUT POTENTIALS , ASYMPTOTIC CORRECTIONS, ETC C IF(IPRIN.GT.1) WRITE(6,93)IR,XL(IR) 93 FORMAT(' IR=',I2,' XLAM=',F9.5) CALL AVN(XM(1,1,IR),AV,S,NBS) IF(IPRIN.GT.1) & WRITE(6,97) 1,AV,XMAV(1,IR),S DO 24 IT=2,ITMAX CALL AVN(XM(1,IT,IR),AV,S,NBS) CALL AVN(XMT(1,IT,IR),AVT,ST,NBS) XMAVT=-XMAV(IT,IR)+XMAV(IT-1,IR) CALL AVN(XMTC(1,IT,IR),AVTC,STC,NBS) XLAM=XL(IR) XCOR=XLAM*XMAVT/(1.-XLAM) IF(XCOR.LT.0.0)XCOR=0.0 XMAVTC= XMAV(IT,IR)-XCOR IF(IPRIN.GT.1) & WRITE(6,97) IT,AV,XMAV(IT,IR),S,XMAVT,ST,XMAVTC,STC 24 CONTINUE 241 CONTINUE 97 FORMAT(' TMAX=',I2,2F9.5,'+/-',5F9.5) enddo RETURN END C-------------------------------------------- COPIED FROM TMPOT CM WROTE 1990 C DATA IN D: A IS AVE S IS S.D. F IS FACTOR TO MULT BY C-------------------------------------------- SUBROUTINE AUTO(D,A,S,F,N) REAL*8 D(N),A,S,AA,F,AN ,AAN A=D(1) AN=0.0 S=D(1)**2 DO 1 I=2,N A=A+D(I) AN=AN+D(I-1)*D(I) S=S+D(I)**2 1 CONTINUE A=A/N AA=(-A*A+S/N) C**<(DI-A)(DI+1-A)>=+A*A-A*-*A *** AAN=(AN-A*A*(N+1)+(D(1)+D(N))*A)/(N-1) F=SQRT( (AA+AAN)/(AA-AAN) ) S=0.0 IF(AA.LE.0.0)RETURN S=SQRT(AA) RETURN END c***************************************************** C c Returns average A and std dev S of D(i) i=1,N c This is sample s.d. (not parent s.d.) c c***************************************************** SUBROUTINE AVN(D,A,S,N) IMPLICIT double precision (A-H,O-Z) Dimension D(N) A=0.0 S=0.0 if(n.eq.1) then a=D(1) s=0.0 return endif DO 1 I=1,N A=A+D(I) S=S+D(I)**2 1 CONTINUE A=A/N AA=(-A*A+S/N) S=0.0 IF(AA.LE.0.0)RETURN S=SQRT(AA) RETURN END C***************************************************** C c Returns median value DM of data D(i), I=1,..N c and upper error DU and lower error DL values that contain c fraction f of data points. For 1 sigma errors use f=0.683 c c written cm 28-6-94 c c***************************************************** SUBROUTINE MEDAV(D,DM,DU,DL,F,N) IMPLICIT double precision (A-H,O-Z) Dimension D(N) if(f.lt.0.0.or.f.ge.1.0) stop 'medav: choice of F wrong' fw=(1.0-f)*0.5 xn=fw*n nxn=int(xn) N1rd=nxn+1 N1ru=nxn+1 c if xn is integer, average 2 values above/below if( abs(xn-float(nxn)).le.1.0E-8)n1rd=nxn xn=0.5*n nxn=int(xn) N2rd=nxn+1 N2ru=nxn+1 c if xn is integer, average 2 values above/below if( abs(xn-float(nxn)).le.1.0E-8)n2rd=nxn xn=(1.0-fw)*n nxn=int(xn) N3rd=nxn+1 N3ru=nxn+1 c if xn is integer, average 2 values above/below if( abs(xn-float(nxn)).le.1.0E-8)n3rd=nxn do 10 i=1,n Count numbers less than or equal each datum nc=0 do 1 nn=1,n if(d(I).ge.d(nn)) nc=nc+1 1 continue if(nc.eq.n1rd) i1rd=i if(nc.eq.n1ru) i1ru=i if(nc.eq.n2rd) i2rd=i if(nc.eq.n2ru) i2ru=i if(nc.eq.n3rd) i3rd=i if(nc.eq.n3ru) i3ru=i 10 continue dm=0.5*(d(i2rd)+d(i2ru)) dl=dm-0.5*(d(i1rd)+d(i1ru)) du=0.5*(d(i3rd)+d(i3ru))-dm return end C********************************************************** C*** FIT R = COSH(MT)/COSH(MT+1) C*** ITV IS DIFF LT/2 AND T ** (IT/IT+1 RATIO USED) c if itv=1 then lowest possible c 10-3-94 add itv=1 direct exact calc c 8-9-94 ip=0 cosh, 1 sinh, -1 LOG C********************************************************** FUNCTION CFS(ip,ITV,R) parameter (iprin=0) IMPLICIT double precision (A-H,O-Z) CFS=0.0 NI=0 c cosh if(ip.eq.0)then IF(R.LE.0.0.OR.R.GT.1.00) RETURN IF(ITV.eq.1) cfs=DLOG((1.0/R)+sqrt((1.0/r)**2-1.0)) if(itv.eq.1) RETURN XT=-DLOG(R) 1 CONTINUE NI=NI+1 R1=R*(1.0+EXP(-2.0*XT*ITV ) )/ (1.0+EXP(-2.0*XT*(ITV-1) ) ) X2=-DLOG(R1) IF (ABS(X2-XT).LT.0.00001) GOTO 2 XT=X2 IF(NI.lt.40) GOTO 1 RT=DEXP(-x2)/((1.0+EXP(-2.0*X2*ITV ) )/(1.0+EXP(-2.0*X2*(ITV-1)))) if(iprin .gt. 1) write(*,*)' cfs',RT,R 2 CONTINUE CFS=X2 RETURN c sinh elseif(ip.eq.1)then rat=(itv-1.0)/float(itv) if(R.gt.rat.or.R.le.0)return c cfs=0 if itv=0 ?? if(itv.eq.1)return XT=-DLOG(R) 10 CONTINUE NI=NI+1 R1=R*(1.0-EXP(-2.0*XT*ITV ) )/ (1.0-EXP(-2.0*XT*(ITV-1) ) ) X2=-DLOG(R1) IF (ABS(X2-XT).LT.0.00001) GOTO 20 XT=X2 IF(NI.lt.40) GOTO 10 RT=DEXP(-x2)/((1.0-EXP(-2.0*X2*ITV ) )/(1.0-EXP(-2.0*X2*(ITV-1)))) if(iprin .gt. 1) write(*,*)' cfs',RT,R 20 CONTINUE CFS=X2 RETURN c LOG elseif(ip.eq.-1)then IF(R.LE.0.0.OR.R.GT.1.00) RETURN CFS=-DLOG(R) else STOP ' wrong CFS option' endif END c-------------------------------------------------------------------- C CM's patent sloppy random number generator c upgraded a la Marsaglia Period > 2**94 c Combine 69069 32 bit congruent generator with c x(n)=x(n-3)-x(n-1) mod 2**31-69 ( assumes integers to 2**31 OK) c Default seeds included c-------------------------------------------------------------------- FUNCTION RNDM2() IMPLICIT double precision (A-H,O-Z) SAVE i,j,k,n c multiplier is 1/largest int-least int: to 0,..1 returned DATA I,J,K,N /521288629,362436069,16163801,1131199299/ m=i-j if(m.lt.0) m=m+2147483579 i=j j=k k=m n=n*69069 + 1013904243 m=m+n RNDM2=m*2.3283064E-10 + 0.5 RETURN c ENTRY ranset(is,js,ks,ns) i=1+iabs(is) j=1+iabs(js) k=1+iabs(ks) n=ns ranset=n RETURN END c*********************************************************** c Example bootstrap fit program c Calls FITEX - generic correlated fit program c FEVAL - fit function c DISPX - organise fit parameters X c GRAPH - plot Axis graphs of fit c AVN, MEDAV - averages, medians c RNDM2 - for bootstrap choices c*********************************************************** c C FIT T CORRELATION WITH NEXP EXPONENTIALS C FIT IT1 TO IT2 (IT CONVENTION 0,..,ITMAX) C INPUT: PLAV(ip1,ip2,it+1,nlat,ir) c XEFF(ip,it+1,ir) used for initial mass values in X c nlat c ip1(ip), ip2(ip) gives data factorisation, npv paths c ie Data =coeff(ip1)*coeff(ip2)*.. c ifs(ip) shows if Cosh (0), Sinh (1) or Log (-1) fit c (Cosh and Sinh need LT - lattice time extent) c c CHOICES: c path values ipth(1,..npth) for number of paths npath(ir) c related ipmap to make x parametrs consequtive c nfac1 nfac2 at t=0, t=T respectively c if nfac1=1 then non-factorising fit c if m(nexp) fixed at xdv0 (idv=0) or free (idv=1) c initial values m0 from xeff; m2 = xdv1 c correlation model: c icorrt, icorrp, idelt, ievals, xlmin (see CORFIT) c ifboot=1 is bootstrap of nbs calls c XD normalises input data (to itmin, for each ip) c iprin (print out amount) (none if iprin=0) c c OUTPUT: X() initial values, Data to fit DAT, common PARMS c care it 0,..convention except for PLAV XEFF C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TFIT(IR,it1,it2,g5mu,coupK) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' C--------------PARAMETERS----------------------------------- c inputs PARAMETER (nrst=nrsth,npstp=npstph,npstx=npsth,itst=itsth) PARAMETER (LT=lx4val) c rest c PARAMETER (icorrt=3,icorrp=2,idelt=4,ievals=6,xlmin=0.001) c PARAMETER (icorrt=0,icorrp=0,idelt=4,ievals=6,xlmin=0.001) PARAMETER (${CORRELATION}) c idv=0 fixes at xdv0 c Parameter (nexp=3,xdv0=1.25,idv=1,iprin=2,npst=npsth,xdv1=1.6) c Parameter (nexp=2,xdv0=1.25,idv=1,iprin=2,npst=npsth,xdv1=1.6) c Parameter (nexp=1,xdv0=1.15,idv=1,iprin=2,npst=npsth,xdv1=1.6) Parameter (nexp=$NEXP,xdv0=1.15,idv=1,iprin=2,npst=npsth,xdv1=1.6) c Parameter (nexp=$NEXP,xdv0=1.15,idv=0,iprin=2,npst=npsth,xdv1=1.6) c Parameter (nexp=$NEXP,xdv0=0.67,idv=0,iprin=2,npst=npsth,xdv1=1.6) c Parameter (nexp=$NEXP,xdv0=0.27,idv=0,iprin=2,npst=npsth,xdv1=1.6) c Parameter (nexp=3,xdv0=0.12,idv=0,iprin=2,npst=npsth,xdv1=1.6) c see ipth DATA statements too c see initial values too c ifboot=99 READIN list & Write X values; ifboot=11 WRITEOUT list c as 1x,25I3 in ibtab(nlst,nbs): boot.lst: 9-6-00 written (see Kieran) PARAMETER (nbs=99,ifboot=$IFBOOT) c PARAMETER (nbs=999,ifboot=$IFBOOT) PARAMETER (NPARS=nexp*(npst+1) ,nparv=npars+3) COMMON/EPP/PLAV(NPSTp,NPSTp,ITST,NLST,NRST),nconf C CORREL IN PLAV( P1,P2, T, CONFIG,R-VAL) C MAX NPV NPV ITMAXP NLAT NR COMMON/mass/xeff(npstx,itst,nrst),xerr(npstx,itst,nrst) common/fact/npv,ip1v(npst),ip2v(npst),ifs(npst),nfacp Common/start/xstart(nparv) c***output to FITEX*****(ipth needed by GRAPH) c X and DAT, COMMON/PARMS/ xdv,xlminv,npar,nexpv,npth,ltv, & icorrtv,icorrpv,ideltv,ievalsv,itmin,itmax,nlat,iprinv common/factf/nfac1c,nfac2c,ix1(npst),ix2(npst),ifsc(npst) & ,ipthc(npst) c***LOCAL***** DIMENSION res(nbs,nparv),xd(npst),dat(0:itst,npst,nlst), & ipth(npst) , x(npars),fita(0:itst,npst),xr(npst),xx(npars) & ,axx(nparv),xorig(npars),ipmap(npst) & ,ibtab(nlst,nbs) c **arrays need length npsth ** c choose paths (as listed in main program arrays ip1v, ip2v) to fit c all c DATA ipth/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15, c & 16,17,18,19,20,21/ c DATA ipmap /1,2,3,4,5,6,15*0/ DATA ipth$IPTH DATA ipmap$IPMAP c THIS now in prologue c mapping of coefficients (as listed in arrays ip1v ip2v) to 1,2,3,.. c if not all are needed c L only (needs npv=6, nfacp=3 in main program) c DATA ipth/1,4,6,11,13,15,15*1/ c map 1->1 3->2 5->3 if L only c DATA ipmap /1,4,2,5,3,6,15*0/ c F only (needs npv=6, nfacp=3 in main program) c DATA ipth/3,8,10,17,19,21,15*1/ c map 2->1 4->2 6->3 if F only c DATA ipmap /4,1,5,2,6,3,15*0/ c obs1 and obs3 (L & F) [ie g5; g4 for pion only] npv=10 nfacp=4 c Data ipth/1,2,3,11,12,15,16,17,20,21,11*1/ c map 5->3 6->4 for obs1 obs3 fits (pi or rho) c Data ipmap /1,2,5,6,3,4,15*0/ c a1 only choose ipponly=1 and fix below ipmap c Data ipth/ 6,9,10,11,12,15,16,17,20,21,11*1/ c map 3->1 4->2 for a1 fit c Data ipmap /3,4,1,2,5,6,15*0/ c iobs3 only c Data ipth/15,20,21,18*1/ c map 5->1 6->2 for obs3 fit only (eg rho) c Data ipmap /3,4,5,6,1,2,15*0/ copy to common (for FITEX) ITmin=IT1 ITmax=IT2 nfac1=nfacp nfac2=nfacp npth=npv icorrtv=icorrt icorrpv=icorrp ideltv=idelt ievalsv=ievals ltv=LT xdv=xdv0 xlminv=xlmin nexpv=nexp iprinv=iprin nlat=nconf nxmax=0 nfac1c=nfac1 nfac2c=nfac2 do 231 ip=1,npth ipv=ipth(ip) c set up mapping to X coefficients in fit (consequtive) ix1(ip)=ipmap(ip1v(ipv)) ix2(ip)=ipmap(ip2v(ipv)) ifsc(ip)=ifs(ipv) ipthc(ip)=ipth(ip) write(6,903)ip,ip1v(ipv),ip2v(ipv),ifs(ipv) & ,ipv,ix1(ip),ix2(ip) 903 format(' data=',i2,' p0=',i2,' pt=',i2,' ifs=',i2, & ' fit=',i2 ,' p0=',i2,' pt=',i2) 231 continue if(iprin.gt.0) write(6,936)ir,nconf 936 FORMAT(' Fits for observable IR=',i1,' number configs=',i4) if(iprin.gt.0) write(6,933)icorrt,icorrp,idelt,ievals, & nexp, idv,(ipth(i),i=1,npth) 933 FORMAT(' Correlations with icorrt=',i2,',icorrp=',I2, &' idelt=',i2,' number exact evals=',i3/ & ' number exponentials=',i2,' idv=',I1,' paths=',4(i3,',')) CONVENIENT NORMALISATION - sum over all latts c it is true t-value here 0,..except PLAV do 1311 ip=1,npth XD(ip)=0.0 do 131 nl=1,nlat XD(ip)=xd(ip)+ & PLAV(ip1v(ipth(ip)),ip2v(ipth(ip)),ITmin+1,nl,IR) 938 format(i5,f20.8) 131 continue xd(ip)=nlat/xd(ip) 1311 continue c xd=1.0 write(6,929)xd(1) 929 format(' Data scaled by XD(1):',g15.6) do 100 nl=1,nlat DO 70 Ip=1,npth do 71 it=itmin,itmax dat(It,ip,nl)=xd(1)* & PLAV(ip1v(ipth(ip)),ip2v(ipth(ip)),IT+1,nl,IR) 71 continue 70 CONTINUE 100 CONTINUE ir=1 c choose parameters C CHOOSE parameters here - initial values np=max(nfac1,nfac2) nm=nexp*np NPAR=nm+nexp If(Npar .gt. nparsv) STOP ' nexp, npath too big for NPARSV' if(idv.eq.0.and.nexp.ge.2) npar=npar-1 C PARAMETERS NP*nexp COEFFICIENTS + V0 + V1 + v.. C INITIAL VALUES ip1=1 itav=(3*itmin+itmax)/4 X (nm+1)= xeff(ipth(ip1),itav,ir) if(nexp.gt.1)then do 144 im=2,nexp x(nm+im)=xdv1 144 continue x(nm+2)=xdv0 endif c order LP LA SP SA for m0, then m1 etc do 101 ip=1,np x(ip)=1.0 if(nexp.gt.1) then do 110 im=2,nexp x( (im-1)*np+ip )=0.0 110 continue endif 101 continue if(np.eq.4.and.nexp.eq.20) then c 3-15 2 exp fit to all pion 10 gauges .040 twist .1608 c c 0.13804 0.86379 1.82660 0.02903 0.12795 c 1.05000 -0.49258 -0.62827 0.01183 0.06234 x(1)=0.86 x(2)= 1.80 x(3)= 0.029 x(4)= 0.128 x(5)= .49 x(6)= 0.62 x(7)= -0.012 x(8)= -0.062 x(9)= 0.138 x(10)= 1.05 endif if(np.eq.6.and.nexp.eq.20) then c pion 1608 .004 3-23 1 exp c 0.15787 0.91610 1.84712 0.04608 0.11080 0.08709 0.22647 c pion 1608 .004 3-15 2 exp PP only c 0.15647 0.89781 1.81594 c 1.14826 0.43553 0.34935 c pion 6x6 m' fixed 7-23 mu=0.004 IN USE c 0.13506 0.98507 1.06144 -0.00172 0.02520 0.09576 0.19741 c 0.45000 0.12982 0.13270 0.02619 0.03957 -0.00478 0.00170 c pion 6x6 2 exp 4-23 mu=0.004 c 0.13545 0.94461 1.01949 -0.00083 0.02484 0.09134 0.18887 c 0.80722 0.31674 0.11995 -0.01699 0.01396 0.00334 -0.01322 x(1)=0.916 x(2)= 1.06 x(3)= -0.002 x(4)= 0.025 x(5)= .096 x(6)= 0.20 x(7)= 0.13 x(8)= 0.13 x(9)= 0.026 x(10)= 0.04 x(11)= -0.005 x(12)= 0.002 x(13)= 0.135 x(14)= 0.45 endif if(np.eq.4.and.nexp.eq.20) then c pion 6x6 2 exp 4-23 mu=0.004 - try as 4x4 c 0.13545 0.94461 1.01949 -0.00083 0.02484 0.09134 0.18887 c 0.80722 0.31674 0.11995 -0.01699 0.01396 0.00334 -0.01322 x(1)=0.916 x(2)= 1.06 x(3)= .096 x(4)= 0.20 x(5)= 0.31 x(6)= 0.13 x(7)= 0.003 x(8)= -0.013 x(9)= 0.135 x(10)= 0.81 endif if(np.eq.6.and.nexp.eq.30) then c rho a1 4-15 1608 c 0.42843 -0.84647 -2.80651 0.26530 0.84600 0.81524 2.86567 c 0.53104 -0.12638 -0.11734 1.41090 3.40250 -0.30098 -1.01887 c 1.11203 0.33667 0.75258 0.36787 0.46559 0.76565 1.76925 x(1)=-0.846 x(2)= -2.806 x(3)= 0.265 x(4)= 0.846 x(5)= .815 x(6)= 2.866 c fix these in funct1 x(7)= 0.0 x(8)= 0.0 x(9)= 1.41 x(10)= 3.40 x(11)= -.3 x(12)= -1.0 x(13)= 0.33 x(14)= 0.75 x(15)= 0.36 x(16)= 0.45 x(17)= 0.77 x(18)= 1.78 x(19)= 0.428 x(20)= 0.53 x(21)= 1.11 endif if(np.eq.6.and.nexp.eq.30) then c rho a1 6-15 1608 mu=0.015 c 0.48405 -0.94858 -2.28119 0.14591 0.13268 1.17083 2.68048 c 0.66050 -0.00200 0.00194 1.10987 1.33869 -0.01763 -0.06288 c 0.97909 0.16659 0.09017 -0.06275 -0.05217 0.49628 0.41527 x(1)=-0.95 x(2)= -2.28 x(3)= 0.14 x(4)= 0.13 x(5)= 1.17 x(6)= 2.68 c fix these in funct1 x(7)= 0.0 x(8)= 0.0 x(9)= 1.11 x(10)= 1.34 x(11)= -.02 x(12)= -0.06 x(13)= 0.17 x(14)= 0.09 x(15)= -0.06 x(16)= -0.05 x(17)= 0.50 x(18)= 0.41 x(19)= 0.48 x(20)= 0.66 x(21)= 0.98 endif if(np.eq.6.and.nexp.eq.2) then c rho a1 3-15 160856 c 0.51363 -0.83034 -1.84153 0.15449 0.33951 0.96762 2.16132 c 0.66981 -0.01395 0.00864 1.39710 1.61016 -0.08939 -0.14778 x(1)=-0.846 x(2)= -1.806 x(3)= 0.15 x(4)= 0.3 x(5)= .96 x(6)= 2.166 c fix these in funct1 x(7)= 0.0 x(8)= 0.0 x(9)= 1.41 x(10)= 1.60 x(11)= -.09 x(12)= -.14 x(13)=0.51 x(14)=0.67 endif *** only use if tmat called with same number of paths as here Do ip=1,npars c x(ip)=xstart(ip) Enddo call fitex(dat,itmin,itmax,npth,x,npars,chi2) call dispx(x,axx,npars) Calculate functions of X of interest: eq f_P x_1 is P x_3 is A_4 c from t3d code: norm careful - assume 1 is LLP, 2 is LL A c 5 is L g_4 - need full twist axx(npar+1)=aXX(5)*exp(itmin*0.5*axx(nm+1))/ & (aXX(nm+1)*sqrt( XD(1)/(2.0*axx(nm+1) ) ) ) c *** use vector ward identity with conserved vector current: c *** needs *2 *2 kappa c #5 which is L_C g_4 is g5mu . L g_5 (#1) / m_pi c norm to fpi=130.7 MeV - charged Axial current. xfpi=2.0*coupK axx(npar+1)=xfpi*aXX(1)*2.0*g5mu*exp(itmin*0.5*axx(nm+1))/ & (aXX(nm+1)*axx(nm+1)*sqrt( XD(1)/(2.0*axx(nm+1) ) ) ) c GMOR fpi^2 mpi^2/4 mu -- care convention for f_pi c norm for ZP ZA- xnff=exp(itmin*0.5*axx(nm+1))/ & sqrt( XD(1)/(2.0*axx(nm+1) ) ) c modify 9-6-08 Do iexp=1,nexp xnff=exp(itmin*0.5*axx(nm+iexp))/ & sqrt( XD(1)/(2.0*axx(nm+iexp) ) ) ioff=(iexp-1)*np do ix=1,np axx(ix+ioff)=xnff*axx(ix+ioff) enddo enddo c mpcac axx(npar+2)=axx(nm+1)*0.5*axx(3)/axx(1) c ZV if ( $FCASE .eq. 1256) then axx(npar+2)=g5mu*2.0*axx(1)/(axx(nm+1)*axx(3)) c alt roman ZV: equiv if exp(m)-exp(-m)=2m (1+m^2/6) c axx(npar+2)=g5mu*4.0*axx(1)/( axx(3)* c & ( exp(axx(nm+1))-exp(-axx(nm+1)) ) ) c test 6x6 case for zv c if ( $FCASE .eq. 0) then c axx(npar+2)=g5mu*2.0*axx(1)/(axx(nm+1)*axx(5)) endif if( $NFACP .eq. 2 .and. $IOBVAL .eq. 1) then c GMOR axx(npar+2)=(axx(npar+1)*axx(nm+1))**2*0.25/g5mu c L only mpcac if($FCASE .eq. 13) & axx(npar+2)=axx(nm+1)*0.5*axx(2)/axx(1) endif c if rho=a1 case - mixing: tan**2 if( $IOBVAL .eq. 2 .and. $NFACP .eq. 6) then axx(npar+1)=-axx(3)*axx(5+6)/(axx(5)*axx(3+6)) endif axx(npar+3)=chi2 if(iprin.gt.0)then write(6,909)axx(npar+1),axx(npar+2) 909 format(' fpi=',f12.6,' m_pcac=',f12.6) call feval(fita,itst,Npst,x,npar) care lots of output : default 4 per line here (nobs) nobs=4 nmax=float(npth)/float(nobs)+0.99 do 5331 nop=1,nmax ngn=(nop-1)*nobs ngp=ngn+nobs if(ngp.gt.npth)ngp=npth write(6,946) (ip1v(ip),ip2v(ip),ip=ngn+1,ngp) 946 format(/' tm E_fit E_data+/-err E_fit E_data+/-err', & ' E_fit E_data+/-err' / 4(' p0=',i4,' pt=',i4) ) do 53 itp=itmin+1,itmax it=itp-1 do 531 ip=1,npth xr(ip)=cfs(ifsc(ip), (lt/2)-it,fita(itp,ip)/fita(it,ip)) 531 continue write(6,935) itp,(xr(ip), & xeff(ipth(ip),itp,ir),xerr(ipth(ip),itp,ir),ip=ngn+1,ngp) 935 format(i3,4(1x,f5.3,1x,f5.3,'(',f4.3,')' ) ) 53 continue 5331 continue call GRAPH(x,npars,ir) endif C**** if(ifboot.ge.1) then c BOOTSTRAP HERE c shuffle plav ===> DAT, use above X as start values,save res c in THIS case shuffle configs not samples nver=nlat/nconf do 1231 i=1,npar xorig(i)=x(i) 1231 continue AN=FLOAT(NCONF) XNBS=1.0/FLOAT(NBS) C READ ibout=77 open(ibout,FILE='${USERSDIR}/boot.lst') c 1-2-01 copied from Kieran's savepi list if(ifboot.eq.99)then read(ibout,944)nbsv if(nbsv.lt.nbs) stop 'nbs mismatch in reading bootstrap list' do 518 kse=1,nbs read(ibout,944)(ibtab(nl,kse),nl=1,nconf) 518 continue 944 format(1x,25I3) endif C-BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB DO 750 KSE =1,NBS DO 700 icf=1,Nconf c SHUFFLE iconf= INT(RNDM2(DUMMY) * AN ) + 1 do 701 inc=1,nver nli=(iconf-1)*nver+inc nl=(icf-1)*nver+inc c WRITE if(ifboot.eq.11) ibtab(nl,kse)=nli c READ if(ifboot.eq.99) nli=ibtab(nl,kse) c use xd(ip) from overall average to normalise DO 470 Ip=1,npth ir=1 do 471 it=itmin,itmax dat(It,ip,nl)=xd(1)* & PLAV(ip1v(ipth(ip)),ip2v(ipth(ip)),IT+1,nli,IR) 471 continue 470 CONTINUE 701 CONTINUE 700 CONTINUE ir=1 c OPTION TO USE ORIGINAL VALUES do 123 i=1,npar x(i)=xorig(i) 123 continue c assume initial X values exist from previous run iprinv=0 call fitex(dat,itmin,itmax,npth,x,npars,chi2) write(*,*)kse call dispx(x,xx,npars) do 111 i=1,npar res(kse,i)=xx(i) c can accumulate averages of functions of X here 111 continue c raw pcac mass res(kse,npar+2)=xx(nm+1)*0.5*xx(3)/xx(1) c ZV if ( $FCASE .eq. 1256) then res(kse,npar+2)=g5mu*2.0*xx(1)/(xx(nm+1)*xx(3)) c if ( $FCASE .eq. 1256 .or. $FCASE .eq. 0) then c res(kse,npar+2)=g5mu*2.0*xx(1)/(xx(nm+1)*xx(5)) endif res(kse,npar+3)=chi2 res(kse,npar+1)=XX(5)*exp(itmin*0.5*xx(nm+1))/ & (XX(nm+1)*sqrt( XD(1)/(2.0*xx(nm+1) ) ) ) c bugfix 12-7-06 from 0.5 to 2.0 res(kse,npar+1)=xfpi*2.0*xx(1)*g5mu*exp(itmin*0.5*XX(nm+1))/ & (xx(nm+1)*XX(nm+1)*sqrt( XD(1) & /(2.0*XX(nm+1) ) ) ) if( $NFACP .eq. 2 .and. $IOBVAL .eq. 1) then c GMOR res(kse,npar+2)=(res(kse,npar+1)*xx(nm+1))**2*0.25/g5mu c L only MPCAC if($FCASE .eq. 13) & res(kse,npar+2)=xx(nm+1)*0.5*xx(2)/xx(1) endif c if rho=a1 case - mixing: tan**2 if( $IOBVAL .eq. 2 .and. $NFACP .eq. 6) then res(kse,npar+1)=-xx(3)*xx(5+6)/(xx(5)*xx(3+6)) endif c norm for ZP ZA xnff=exp(itmin*0.5*res(kse,nm+1))/ & sqrt( XD(1)/(2.0*res(kse,nm+1) ) ) c modify 9-6-08 Do iexp=1,nexp xnff=exp(itmin*0.5*xx(nm+iexp))/ & sqrt( XD(1)/(2.0*xx(nm+iexp) ) ) ioff=(iexp-1)*np do ix=1,np res(kse,ix+ioff)=xnff*xx(ix+ioff) enddo enddo c call flush(6) 750 continue C-BBBBBBBBBBBBB END BOOTSTRAP LOOP BBBBBBBBBBBBBBBBBBBBBB- C WRITE if(ifboot.eq.11)then write(ibout,944)nbs do 517 kse=1,nbs write(ibout,944)(ibtab(nl,kse),nl=1,nconf) 517 continue endif close(ibout) c WRITE X values - for subsequent bootstrap analysis if(ifboot.eq.99.or.ifboot.eq.11)then ixout=88 open(ixout,FILE='${USERSDIR}/bootout${IFK}_${IOBVAL}.dat') write(ixout,944)nbs write(ixout,9441)(axx(iv),iv=1,npar+2) do 5171 kse=1,nbs write(ixout,9441)(res(kse,iv),iv=1,npar+2) 9441 format(5g15.7) 5171 continue close(ixout) endif c copy of salient output to fort.25 write(25,*)'Fits for meson ${PARTICLE}; iobval=${IOBVAL}' write(25,*)'ifk=${IFK};', & ' ${INFO}' write(25,*)'nexp= ',nexp, ' CASE=', ${FCASE} write(25,*)'${CORRELATION}' write(25,*)'itmin= ',itmin,' itmax= ',itmax write(25,*)'n_gauges = ',nlat,' in blocks of $NBLOCK' write(25,978) 978 format(' fit coefficients; masses; derived quantities') write(25,*)' nbs x val +err -err ; s.d.' write(6,978) c ideally use nbs/(nbs-1) in estimate of BS variance .5% difference! do 235 il=1,npar call avn(res(1,il),a0,s0,nbs) f=0.683 call medav(res(1,il),am,au,al,f,nbs) if(il.eq.nexp*np+1) write(6,959) if(il.eq.nexp*np+1) write(25,959) 959 format(' masses') write(6,977) nbs,il,a0,axx(il),s0 write(6,977) nbs,il,am,axx(il),au,al write(25,9770) nbs,il,axx(il),au,al,s0 9770 format(2i6,f12.5,' (+',f10.5,'-',f10.5,') ;',f10.5) 977 format(2i6,2f12.5,'(',f10.5,') ',f10.5) 235 continue if(idv .eq. 0) write(25,*) 'fixed excited mass=',xdv0 nex=0 if($IOBVAL .eq. 1) then if( $NFACP .eq. 6) write(25,*)' f_pi, m_PCAC' if( $NFACP .eq. 6) nex=2 if( $NFACP .eq. 2) write(25,*)' f_pi, GMOR' if( $NFACP .eq. 2) nex=2 if( $NFACP .eq. 1) write(25,*)' f_pi, GMOR' if( $NFACP .eq. 1) nex=2 if( $NFACP .eq. 4 .and. $FCASE .eq. 1234) & write(25,*)' f_pi, m_PCAC' if( $NFACP .eq. 4 .and. $FCASE .eq. 1256) & write(25,*)' f_pi, Z_V' if( $NFACP .eq. 4) nex=2 endif if($IOBVAL .eq. 2) then if($NFACP .eq. 6) write(25,*)' tan^2 alpha' if($NFACP .eq. 6) nex=1 endif do il=npar+1, npar+nex call avn(res(1,il),a0,s0,nbs) f=0.683 call medav(res(1,il),am,au,al,f,nbs) write(6,977) nbs,il,a0,axx(il),s0 write(6,977) nbs,il,am,axx(il),au,al write(25,9770) nbs,il,axx(il),au,al,s0 enddo ndof=npth*(itmax-itmin+1) - npar write(25,9771) axx(npar+3),ndof write(6,9771) axx(npar+3),ndof 9771 format(' chi^2=',f12.5,'/dof: ',i5) endif return end c**************************************************************** c CALLED TFIT c display effective masses, errors and fits c use AXIS conventions (these are GNUPLOT plus extras) called by TFIT Plot Data wherever available, fit only where fitted c offset different paths by t steps c care if called several times - different ir to fort.7, 8.. c CALLS FEVAL c**************************************************************** subroutine graph(x,npars,ir) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' parameter (itst=itsth,npst=npsth,nrst=nrsth,npstx=npsth) dimension x(npars),fita(0:itst,npst) COMMON/PARMS/ xdv,xlmin,npar,nexp,npth,lt, & icorrt,icorrp,idelt,ievals,itmin,itmax,nlat,iprin COMMON/mass/xeff(npstx,itst,nrst),xerr(npstx,itst,nrst) common/factf/nfac1c,nfac2c,ix1(npst),ix2(npst),ifs(npst) & ,ipth(npst) call feval(fita,itst,Npst,x,npar) iuf=7+ir-1 c *** set increment for nice plots here tstep=0.02 titmax=itmax+1 titmin=itmin-1 write(iuf,911) titmin,titmax 911 FORMAT('# w 0.9 '/'# h 0.9'/'# r 0.0'/'# e '/'# x ',2f12.5) CCC **** set by hand ****** less for pi than rho etc ymin=0.0 ymax=0.25 iobval=$IOBVAL if( iobval .eq. 2) ymax=0.8 if( iobval .eq. 2) tstep=0.03 write(iuf,9110) ymin, ymax 9110 Format('# y ',2f12.5) c PLOT DATA WITH ERRORS xav=-10.0 do 225 ip=1,npth ipv=ipth(ip) do 125 i=itmin,itmax tplot=i+tstep*(ip-1) if(xeff(ipv,i,ir).gt.xav) & WRITE(iuf,988) tplot,xeff(ipv,i,ir),xerr(ipv,i,ir) 988 FORMAT(3G12.4) 125 continue 225 continue write(iuf,912) 912 format('# e0 '/'# s'/'# cs 1.0') c PLOT DATA CENTRAL POINTS c care '...' cant cope with \ NEED \\ to quote \ c care here document translation messes this up - need \\ twice? do 235 ip=1,npth ipv=ipth(ip) 901 format('# c "\\\\pl"') 902 format('# c "\\\\cr"') 903 format('# c "\\\\di"') 904 format('# c "\\\\oc"') 905 format('# c "\\\\sq"') 906 format('# c "\\\\bu"') if(ip.eq.1)write(iuf,901) if(ip.eq.2)write(iuf,902) if(ip.eq.3)write(iuf,903) if(ip.eq.4)write(iuf,904) if(ip.eq.5)write(iuf,905) if(ip.ge.6)write(iuf,906) do 135 i=itmin,itmax tplot=i+tstep*(ip-1) if(xeff(ipv,i,ir).gt.xav) & WRITE(iuf,938) tplot,xeff(ipv,i,ir) 938 format(2g12.4,/'#') 135 continue 235 continue write(iuf,922) 922 format('# e0 '/'# s'/'# c0') c**** add titles by hand ****** c 1 full 2 dotted 3 dash 931 format('# m 1') 932 format('# m 2') 933 format('# m 3') 934 format('# m 4') 935 format('# m 5') do 531 ip=1,npth if(ip.eq.1)write(iuf,931) if(ip.eq.2)write(iuf,932) if(ip.eq.3)write(iuf,933) if(ip.eq.4)write(iuf,934) if(ip.ge.5)write(iuf,935) do 53 itp=itmin+1,itmax it=itp-1 tplot=itp xr=cfs(ifs(ip), (lt/2)-it,fita(itp,ip)/fita(it,ip)) c xr=-log(fita(itp,ip)/fita(it,ip)) if(xr.gt.xav) & write(iuf,988) tplot,xr 53 continue write(iuf,977) 977 format('#') 531 continue return end c**************************************************************** c c CALLED TFIT c display parameters X, order them in increasing mass (as XX) c match conventions of FEVAL c c**************************************************************** subroutine dispx(x,xx,npars) IMPLICIT double precision (A-H,O-Z) Include 'fitex.h' PARAMETER (npst=npsth) dimension x(npars),xx(npars),xm(nparsv),imv(nparsv) COMMON/PARMS/ xdv,xlmin,npar,nexp,npth,lt, & icorrt,icorrp,idelt,ievals,itmin,itmax,nlat,iprin common/factf/nfac1,nfac2,ip1v(npst),ip2v(npst),ifs(npst) & , ipth(npst) np=max(nfac1,nfac2) nm=np*nexp xm(1)=x(nm+1) if(iprin.gt.0) write(6,900) 900 format(' FIT PARAMETERS mass, coeff(path ==>) ') if(nexp.eq.1) then do 1 ip=1,npar xx(ip)=x(ip) 1 continue if(iprin.gt.0) write(6,99) xm(1),(xx(ip),ip=1,np) else xm(1)=x(nm+1) do 2 im=2,nexp xm(im)=x(nm+im) 2 continue if(npar.eq.nm+nexp-1) xm(nexp)=xdv do 4 i=1,nexp imv(i)=i 4 continue c order xm values 20 continue nch=0 do 30 i=2,nexp im=i-1 c mod 3-08 if( xm(im).lt.xm(i) ) then nch=1 c swap i and im st=xm(im) xm(im)=xm(i) xm(i)=st ist=imv(im) imv(im)=imv(i) imv(i)=ist endif 30 continue if(nch.eq.1) goto 20 c store in xx and display do 10 i=1,nexp xx(nm+i)=xm(i) ivv=np*(i-1) do 11 ip=1,np xx(ivv+ip)=x(np*(imv(i)-1)+ip) 11 continue if(iprin.gt.0) write(6,99) xm(i),(xx(ivv+ip),ip=1,np) 99 format(f10.5,1x,6f10.5) 10 continue endif return end c***************************************************** C c Returns jacknife error on DC/C1 c c c***************************************************** SUBROUTINE AVRAT(DC,C1,NS,N,A,S) IMPLICIT double precision (A-H,O-Z) Dimension C1(NS), DC(NS) A=0.0 S=0.0 if(n.eq.1) then a=DC(1)/C1(1) return endif xn=1.0D0/float(N) xm=1.0D0/float(N-1) ds=0.0 cs1=0.0 dcs=0.0 DO K=1,N cs1=cs1+c1(k) dcs=dcs+dc(k) Enddo A=(dcs/cs1) c jacknife - single elimination ss=0.0 sa=0.0 Do K=1,n c prefer xm from N-1 here v=( (dcs-dc(k)) / (cs1-c1(k)) ) c based on avrat c write(*,*)a,v dv=a*n-v*(n-1) sa=sa+dv ss=ss+dv*dv Enddo ss=ss-sa*sa/float(N) If( n .gt. 1) S=SQRT(ss/float(n*(n-1))) RETURN END c***************************************************** C c Returns jacknife error on DC/C1 blocked c c c***************************************************** SUBROUTINE AVRATBO(DC0,C10,N0,A,S,nb) IMPLICIT double precision (A-H,O-Z) Parameter (Nv=4000) Dimension C10(N0), DC0(N0) Dimension C1(Nv), DC(Nv) n=N0/NB xn=1.0D0/float(nb) if( n .gt. nv) STOP ' nv too small in avratb' k=0 do i=1,n c1(i)=0.0 dc(i)=0.0 do j=1,nb k=k+1 c1(i)=c1(i)+c10(k)*xn dc(i)=dc(i)+dc0(k)*xn Enddo Enddo A=0.0 S=0.0 if(n.eq.1) then a=DC(1)/C1(1) return endif xn=1.0D0/float(N) xm=1.0D0/float(N-1) ds=0.0 cs1=0.0 dcs=0.0 DO K=1,N cs1=cs1+c1(k) dcs=dcs+dc(k) Enddo A=(dcs/cs1) c jacknife - single elimination ss=0.0 sa=0.0 Do K=1,n c prefer xm from N-1 here v=( (dcs-dc(k)) / (cs1-c1(k)) ) c based on avrat c write(*,*)a,v dv=a*n-v*float(n-1) sa=sa+dv ss=ss+dv*dv Enddo ss=ss-sa*sa/float(N) If( n .gt. 1) S=SQRT(ss/float(n*(n-1))) RETURN END c***************************************************** C c Returns jacknife error on DC/C1 blocked c NEW - seems to be same for Z_V c c***************************************************** SUBROUTINE AVRATB(DC0,C10,N0,A,S,nb) IMPLICIT double precision (A-H,O-Z) Parameter (Nv=4000) Dimension C10(N0), DC0(N0) Dimension C1(Nv), DC(Nv) n=N0/NB xn=1.0D0/float(nb) if( n .gt. nv) STOP ' nv too small in avratb' k=0 do i=1,n c1(i)=0.0 dc(i)=0.0 do j=1,nb k=k+1 c1(i)=c1(i)+c10(k)*xn dc(i)=dc(i)+dc0(k)*xn Enddo Enddo A=0.0 S=0.0 if(n.eq.1) then a=DC(1)/C1(1) return endif xn=1.0D0/float(N) xm=1.0D0/float(N-1) ds=0.0 cs1=0.0 dcs=0.0 DO K=1,N cs1=cs1+c1(k) dcs=dcs+dc(k) Enddo A=(dcs/cs1) c jacknife - single elimination ss=0.0 sa=0.0 Do K=1,n c prefer xm from N-1 here v=( (dcs-dc(k)) / (cs1-c1(k)) ) c based on avrat c write(*,*)a,v dv=a-v c sa=sa+dv ss=ss+dv*dv Enddo ss=ss/float(n) c ss=ss-sa*sa/float(N) c If( n .gt. 1) S=SQRT(ss/float(n*(n-1))) If( n .gt. 1) S=SQRT(float(n-1)*ss) RETURN END EOFH # # # don't assume fitex.o exists on $USERSDIR cp $USERSDIR/fitex.f fitex.f f77 $TEMPDIR/fitex.f -c > out.$pn 2> err.$pn f77 $TEMPDIR/tmprf.f -c > out.$pn 2> err.$pn #f77 tmprf.o fitex.o -o my$pn -lnag18 # >> out.$pn 2 >> err.$pn f77 tmprf.o fitex.o -o my$pn -L /apps/nag/lib -lnag18 # >> out.$pn 2 >> err.$pn my$pn >> $TEMPDIR/out.$pn echo `date` >> $TEMPDIR/out.$pn cp $TEMPDIR/out.$pn $USERSDIR/out.$pn$IFK cp $TEMPDIR/err.$pn $USERSDIR/err.$pn$IFK cp $TEMPDIR/fort.25 $USERSDIR/summary.$pn$IFK # send fort.7 , 8 etc AXIS file to SUNC # only more than 1 if nr > 1 GC=1 while [ $GC -le 1 ] do GF=`expr $GC + 6 ` if [ -f fort.$GF ] then cp $TEMPDIR/fort.$GF $USERSDIR/tmpr$GC.dat # send axis plot to local machine if [ $TDIR -eq "cmi" ] then rcp fort.$GF sune.amtp.liv.uk:plot/tmpr$GC.dat fi rm fort.$GF fi GC=`expr $GC + 1 ` GF=`expr $GC + 6 ` done cp fort.11 $USERSDIR cp fort.30 $USERSDIR cp fort.31 $USERSDIR cp fort.32 $USERSDIR cp fort.68 $USERSDIR cp fort.69 $USERSDIR cp fort.70 $USERSDIR cp fort.31 $USERSDIR/mpcac_v_g$IFK rcp fort.31 sune.amtp.liv.uk:tmp/mpcac_v_g$IFK cd $USERSDIR #rm my$pn out.$pn err.$pn tmprf.f tmprf.o core fitex.* NOTAR=0 NOTAR=1 if [ $NOTAR -eq 0 ] then rm $TEMPPATH/* fi