program main_WUp
C Target in the Up position on the West side
implicit none
include 'const.inc'
integer loop, dncross, nlamx, nlamy
integer ngbx, ngbx0,ngbx1,ngbx2,ngbx3
integer ngby, ngby0,ngby1,ngby2,ngby3
integer nibx, nibx0,nibx1,nibx2,nibx3
integer niby, niby0,niby1,niby2,niby3
integer nobx, nobx0,nobx1,nobx2,nobx3
integer noby, noby0,noby1,noby2,noby3
integer nrfl1,nrfl2
integer scatflg1, scatflg2
double precision ep, wp, kp, k1, thp,php,qlast
double precision dcrosw,dcrose,dscrosw,dscrose
integer i, j, counter,lim,Event,cin,cot,trig !bec
double precision Nenscat1,Ntar1,Nenscat2,Ntar2,Ntar12,NEdet,Nxo
double precision encnt, Nsdet,Nthdet, Ntar1out, Ntar2out
double precision Nalscat1, Nalscat2, Nthscat1, Nthscat2
double precision Nalout1, Nalout2, Nalin1, Nalin2, Nairin1 !rcm
double precision Nairout1, Nairscat1
double precision NPSMout, Nincoilout, NOutcoilout, NTargetout ! bec
double precision numtar1flg, numtar2flg, numendif1, numbeamsc
double precision numendif2, numendifdet, numbounce
double precision numbounce2, numbounce3,Ndetin
double precision Nemp1in,Nemp2in,Nemp1out,Nemp2out, NtarWEin
double precision INPKFlux,PSMAVGFlux,PSMPKFlux
double precision INCOAVGFlux,INCOPKFlux,TARAVGFlux,TARPKFlux
double precision OUTCOAVGFlux,OUTCOPKFlux,DETAVGFlux,DETPKFlux
double precision xmax,ymax,capflpsm,capflinc,capflout,capfldet
double precision fsw,fse,ufsw,ufse,rotnsw,rotnse,urotnsw,urotnse
double precision rotsw,rotse,urotsw,urotse
double precision fxw,fxe,ufxw,ufxe
double precision rotxw,rotxe,urotxw,urotxe
c
real xlam, HRNDM1, vect(7), lm(3), const, xinit
real yinit, plr, xtest, qvec(2)
character*28 outfname
double precision th0, ph0, thindex, enstart, entar1, entar2
double precision en0, x0, y0, z0, thx0, thy0, xstart
double precision en1, x1, y1, z1, thx1, thy1
real xtmp, ytmp
double precision time1
logical tar1flg, tar2flg,enflg1,enflg2,thflg1,thflg2,alflg1,alflg2
logical airflg1, histflg, histflg2, empflg1, empflg2, bounceflg
logical bounceflg2, bounceflg3, beamscflg
integer ranno
integer seedc
double precision Nempty1, Nempty2
C parameter(ranno=1)
double precision r2d, lambda, sig, intlen, thscat
double precision phi, ztar1,scwgt
double precision pltot,thtot,pl_del,rotang,pltottar
double precision lsrot,lsrotsq,lsrotavg,lsrottsc,lsrottscsq
double precision lsrotgsc,lsrotgscsq,lsde,lsdesq
double precision lsrotxsc,lsrotxscsq
double precision rsrot,rsrotsq,rsrotavg,rsrottsc,rsrottscsq
double precision rsrotnsc,rsrotnscsq,lsrotnsc,lsrotnscsq
double precision rsrotgsc,rsrotgscsq,rsde,rsdesq
double precision rsrotxsc,rsrotxscsq
double precision thxinput, thyinput
double precision nr,nrts,nrgs,nl,nlts,nlgs,nlxs,nrxs,nder,ndel,nb
double precision nrsq,nrtssq,nrgssq,nlsq,nltssq,nlgssq
double precision nrxssq,nlxssq
double precision ndersq,ndelsq,nrtsraw,nltsraw,nrgsraw,nlgsraw
double precision nrxsraw,nlxsraw
double precision nrnsc,nlnsc,numE2W, numW2E,nrnscsq,nlnscsq
double precision lall,rall,lallsq,rallsq,lsc,rsc,lscsq,rscsq
double precision nrPA,nrPAsq,rPAavg,nlPA,nlPAsq,lPAavg !rcm PAsin
double precision uncrsrotavg,uncrPAavg,unclsrotavg,unclPAavg
double precision nrnscPA,nrnscPAsq,rsnscPAavg,rsrotnscavg
double precision nlnscPA,nlnscPasq,lsnscPAavg,lsrotnscavg
double precision uncrsnsc,uncrsnscPA,unclsnsc,unclsnscPA
double precision nrtsPA,nrtsPAsq,rtsPAavg,rsrottscavg
double precision nltsPA,nltsPAsq,ltsPAavg,lsrottscavg
double precision uncrstsc,uncrstscPA,unclstsc,unclstscPA
double precision nrxsPA,nrxsPAsq,rxsPAavg,rsrotxscavg
double precision nlxsPA,nlxsPAsq,lxsPAavg,lsrotxscavg
double precision uncrsxsc,uncrsxscPA,unclsxsc,unclsxscPA
double precision nrgsPA,nrgsPAsq,rgsPAavg,rsrotgscavg
double precision nlgsPA,nlgsPAsq,lgsPAavg,lsrotgscavg
double precision uncrsgsc,uncrsgscPA,unclsgsc,unclsgscPA
double precision nrallPA,nrallPAsq,nlallPA,nlallPAsq,rallavg
double precision rallPAavg,lallavg,lallPAavg,nrscPA,rscPAavg
double precision rscavg,nlscPA,lscPAavg,lscavg,nlscPAsq,nrscPAsq
double precision uncrallavg,uncrallPAavg,unclallavg,unclallPAavg
double precision argrsrot,argrsrotPA,argrsnsc,argrsnscPA
double precision uncrscavg,uncrscPAavg,unclscavg,unclscPAavg
double precision arglsrot,arglsrotPA,arglsnsc,arglsnscPA
double precision argltsc,argltscPA,arglgsc,arglgscPA
double precision argrtsc,argrtscPA,argrgsc,argrgscPA
double precision arglxsc,arglxscPA
double precision argrxsc,argrxscPA
double precision argrall,argrallPA,arglall,arglallPA
double precision zdet, rsqrt,ratan,Nguideout
double precision NguideoutE,NguideoutW,NincoiloutE,NincoiloutW
double precision NoutcoiloutW,NoutcoiloutE
double precision zlm
integer beamflg,nbx,nby,scflag,gsc,isc,osc
logical nanflg
real beamlength, tarlength
c----
double precision lwde,rwde,lwdesq,rwdesq,lpl,lplsq,rpl,rplsq
double precision lplns,rplns,lplnssq,rplnssq
integer nlall,nrall,nlsc,nrsc,ndelw,nderw
c----
double precision IntBz,temp,norm !for gradient runs -KG
double precision PA, xtar1,xtar2 !Polarization product -KG
CHARACTER day*2,month*2,hour*2,minute*2, filename*20
c-----------------------------------
integer ispec
double precision spec(1000,2), alambda
c-----------------------------------
integer qmin, qmax, hmax
integer hws1, hws2
double precision k0, qdummy, q, q1, wp1, ep1, kp1, th_test1
double precision q2, wp2, ep2, kp2, th_test2
double precision qtest, wtest
real wgt, wgt2
c-----------------------------------
c
c -------------------------
c fill the ntuple 'spin.rz' //// ntuple variables ....
integer max_paw, iquest
integer istat, icycle, datetime(8)
parameter(max_paw = 2000000) ! PAW common block size
c parameter(max_paw = 1000000) ! PAW common block size
real quest, p
character date*12, time*12, zone*12
c common/quest/iquest(100)
common/pawc/p(max_paw)
common/bfield/temp(1101),norm(1101)
if (bzgradnum.eq.1) then
open(unit=37,file=bzfile1,status='old')
else if (bzgradnum.eq.2) then
open(unit=37,file=bzfile2,status='old')
else if (bzgradnum.eq.3) then
open(unit=37,file=bzfile3,status='old')
else if (bzgradnum.eq.4) then
open(unit=37,file=bzfile4,status='old')
end if
if (bzgradnum.ne.0) then
read(37,*)
do nb = 1, 1101
read(37,*) temp(nb)
c write(6,*) temp(nb)
end do
read(37,*)
do nb = 1, 1101
read(37,*) norm(nb)
c write(6,*) norm(nb)
end do
end if
close(unit=37)
c time1=time()
c ranno=mod(time1,2000000)
c ranno=mod(time1,1089926518)
CALL Date_and_time(date,time,zone,datetime) !bec
day= CHAR(INT(datetime(3)/10)+48)//CHAR(mod(datetime(3),10)+48)
month= CHAR(int(datetime(2)/10)+48)//CHAR(MOD(datetime(2),10)+48)
hour= CHAR(INT(datetime(5)/10)+48)//CHAR(mod(datetime(5),10)+48)
minute= CHAR(INT(datetime(6)/10)+48)//CHAR(mod(datetime(6),10)+48)
filename ='xportWUp-'//month//day//'-'//hour//minute//'.o'
open(unit=3,FILE=filename,Status='new')
write(*,*) ' Nevents: ', nevts ! max_paw !bec
write(*,*) ' Start time: ', date,time !bec
write(3,*) ' Nevents: ', nevts ! max_paw !bec
write(3,*) ' Start time: ', date,time !bec
write(*,*) ' Upstream target on West (Right), x>0.'
write(3,*) ' Upstream target on West (Right), x>0.'
c print *,"asking for seed"
c read(*,*) ranno;
c ranno=seedc()
ranno=datetime(6)*datetime(7)*datetime(8)+datetime(7)+datetime(8) !bec
write(*,*) 'ranNo is ', ranno
write(3,*) 'ranNo is ', ranno
c stop
c stop/N/B/msarsour/spin/new_ih/
call hlimit(max_paw)
c iquest(10) = 6500000
c
c---------------------------
beamlength = talz1+talz2+2.0*talz3+2.0*talz4+2.0*talz5+
& 4.0*talz6+talz7+
& tairz1+tairz2+tairz3+tairz4+tairz5+tairz6+
& gl+2.0*tarz+iz0+oz+sz
tarlength = 2.0*tarz + 2.0*talz4 + talz3 + 2.0*talz5 + tairz3
write(6,*) ' Calculated beamlength = ',beamlength
write(6,*) ' Calculated length of target region = ',tarlength
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
call hrget(99,'hbooks/lam.hbook','') !New initial spectrum after PSM -KG
c call hrget(102,'hbooks/xyspot.hbook','')
call hrget(151,'hbooks/gradx.hbook','')
call hrget(152,'hbooks/gradx.hbook','')
call hrget(161,'hbooks/ng6-parab-xy.hbook','')
call hrget(162,'hbooks/ng6-parab-xy.hbook','')
c call hrget(161,'hbooks/ngc-parab8x8-xy.hbook','')
c call hrget(162,'hbooks/ngc-parab8x8-xy.hbook','')
c call hrget(161,'hbooks/ngc-parab11x11-xy.hbook','')
c call hrget(162,'hbooks/ngc-parab11x11-xy.hbook','')
call hrget(699,'hbooks/sqw.hbook','') ! New S(q,w)
cbec10 call hrget(700,'hbooks/sqw.hbook','') ! New S(q,w)
call hrget(200,'hbooks/psm.hbook','') !x-y out of PSM -bec
cc-----------------------------------------------------------------------
c read(*,*) ranno
call rluxgo(3,ranno,0,0) ! RANLUX luxury level 2
c-----------------------------------------------------------------------
c Print *, "setting up Hbooks"
call hbook1(101,'lamda_start',100,0,16.0,0.) !-KG
call hbook1(305,'West q sim(1/Ang)',100,0,0.1,0.) ! bec test q choices
call hbook1(306,'West w sim(1/psec)',100,-0.5,0.5,0.) ! bec test q choices
call hbook1(307,'West Eo sim(meV)',3000,0.,40.0,0.) ! bec test q choices
call hbook1(308,'West Ep sim(meV)',3000,0.,40.0,0.) ! bec test q choices
call hbook1(309,'West ko sim(1/Ang)',100,0,4.0,0.) ! bec test q choices
call hbook1(310,'West kp sim(1/Ang)',100,0,4.0,0.) ! bec test q choices
call hbook1(205,'East q sim(1/Ang)',100,0,0.1,0.) ! bec test q choices
call hbook1(206,'East w sim(1/psec)',100,-0.5,0.5,0.) ! bec test q choices
call hbook1(207,'East Eo sim(meV)',3000,0.,40.0,0.) ! bec test q choices
call hbook1(208,'East Ep sim(meV)',3000,0.,40.0,0.) ! bec test q choices
call hbook1(209,'East ko sim(1/Ang)',100,0,4.0,0.) ! bec test q choices
call hbook1(210,'East kp sim(1/Ang)',100,0,4.0,0.) ! bec test q choices
call hbook1(1001,'x (cm) Uptar ',100,-5.5,5.5,0.)
call hbook1(1,'x (cm) start',100,-5.5,5.5,0.)
call hbook1(10,'y (cm) start',100,-5.5,5.5,0.)
call hbook1(2,'x (cm) guide out',100,-5.5,5.5,0.)
call hbook1(3,'x (cm) incoil out',100,-5.5,5.5,0.)
call hbook1(4,'x (cm) target in',100,-5.5,5.5,0.)
call hbook1(5,'x (cm) outcoil out',100,-5.5,5.5,0.)
call hbook1(6,'x (cm) ASM out',100,-5.5,5.5,0.)
call hbook1(20,'E (meV) start Both',2000,0.,40.,0.)
call hbook1(31,' thx start (rad)',200,-0.1,0.1,0.)
call hbook1(32,' thy start (rad)',200,-0.1,0.1,0.)
call hbook2(2220,'Flux(xy) PSM out',300,-2.98,3.0,300,
& -2.98,3.0,0.) !-bec
call hbook2(2221,'Flux(xy) incoil out',300,-2.98,3.0,300,
& -2.98,3.0,0.) !-bec
call hbook2(2222,'Flux(xy) target out',350,-3.48,3.5,350,
& -3.48,3.5,0.) !-bec
call hbook2(3333,'Flux(xy) outcoil out',350,-3.48,3.5,350,
& -3.48,3.5,0.) !-bec
call hbook2(4444,'Flux(xy) Detector',350,-3.48,3.5,350,
& -3.48,3.5,0.) !-bec
call hbook1(441,'Flux(x) PSM (/s)',300,-3.0,3.0,0.)
call hbook1(442,'Flux(y) PSM (/s)',300,-3.0,3.0,0.)
call hbook1(443,'Flux(x) InGuide (/s)',300,-3.0,3.0,0.)
call hbook1(444,'Flux(y) InGuide (/s)',300,-3.0,3.0,0.)
call hbook1(445,'Flux(x) OutGuide (/s)',350,-3.5,3.5,0.)
call hbook1(446,'Flux(y) OutGuide (/s)',350,-3.5,3.5,0.)
call hbook1(1002,'y (cm) UpWtar',100,-6.,6.,0.)
c call hbook1(1003,'thx (mrad) Up',160,-1.6,1.6,0.)
c call hbook1(1004,'thy (mrad) Up',160,-1.6,1.6,0.)
call hbook1(1005,'E (meV)/ non-scat tar1',2000,0.,40.,0.)
c call hbook1(1006,'Path len (cm)/ non-scat UpB',100,35.0,100.0,0.)
c call hbook1(1007,'th_rot (mrad)/ non-scat UpB',100,0.,10000.,0.)
c call hbook1(1008,'Path len (cm)/ non-scat UpR',100,35.0,100.,0.)
c call hbook1(1009,'th_rot (mrad)/ non-scat UpR',100,0.,10000.,0.)
c call hbook1(1010,'Path len (cm)/ non-scat UL',100,35.0,100.0,0.)
c call hbook1(1011,'th_rot (mrad)/ non-scat UL',100,0.,10000.,0.)
call hbook1(1012,'E (meV)/ scat UpWE',2000,0.,40.,0.)
c call hbook1(1013,'Path len (cm)/ scat UpB',100,35.0,100.0,0.)
c call hbook1(1014,'th_rot (mrad)/ scat UpB',100,0.,10000.,0.)
c call hbook1(1015,'Path len (cm)/ scat UpR',100,35.0,100.0,0.)
c call hbook1(1016,'th_rot (mrad)/ scat UpR',100,0.,10000.,0.)
c call hbook1(1017,'Path len (cm)/ scat UpL',100,35.0,100.0,0.)
c call hbook1(1018,'th_rot (mrad)/ scat UpL',100,0.,10000.,0.)
*
call hbook1(2001,'x (cm) Dn',100,-5.5,5.5,0.)
call hbook1(2002,'y (cm) Dn',100,-6.,6.,0.)
c call hbook1(2003,'thx (mrad) Dn',160,-1.6,1.6,0.)
c call hbook1(2004,'thy (mrad) Dn',160,-1.6,1.6,0.)
call hbook1(2005,'E (meV)/ non-scat DnWE',2000,0.,40.,0.)
c call hbook1(2006,'Path len (cm)/ non-scat DnB',100,41.5,41.7,0.)
c call hbook1(2007,'th_rot (mrad)/ non-scat DnB',100,0.,1000.,0.)
c call hbook1(2008,'Path len (cm)/ non-scat DnR',100,41.5,41.7,0.)
c call hbook1(2009,'th_rot (mrad)/ non-scat DnR',100,0.,1000.,0.)
c call hbook1(2010,'Path len (cm)/ non-scat DnL',100,41.5,41.7,0.)
c call hbook1(2011,'th_rot (mrad)/ non-scat DnL',100,0.,1000.,0.)
call hbook1(2012,'E (meV)/ scat DnWE',2000,0.,40.,0.)
c call hbook1(2013,'Path len (cm)/ scat DnB',100,41.,45.,0.)
c call hbook1(2014,'th_rot (mrad)/ scat DnB',100,0.,1000.,0.)
c call hbook1(2015,'Path len (cm)/ scat DnR',100,41.,45.,0.)
c call hbook1(2016,'th_rot (mrad)/ scat DnR',100,0.,1000.,0.)
c call hbook1(2017,'Path len (cm)/ scat DnL',100,41.,45.,0.)
c call hbook1(2018,'th_rot (mrad)/ scat DnL',100,0.,1000.,0.)
*
call hbook1(3001,'x (cm) Det',100,-6.0,6.0,0.)
call hbook1(3002,'y (cm) Det',100,-6.,6.,0.)
call hbook1(3003,'thx (rad) Det',200,-0.15,0.15,0.)
call hbook1(3004,'thy (rad) Det',200,-0.15,0.15,0.)
call hbook1(3005,'E (meV) DetWE',2000,0.,40.,0.)
call hbook1(4004,'Lambda (Angs) Initial',1000,0.,16.,0.) !KG
call hbook1(4005,'Lambda (Angs) DetWE',1000,0.,16.,0.)
c call hbook1(3006,'Path length (cm) DeWE',200,480.58,480.70,0.) !KG
c call hbook1(3007,'th_rotation (mrad) DetWE',200,-400.,0.0,0.) !KG
c call hbook1(3008,'Path length (cm) DetEast',200,114.58,114.7,0.)
c call hbook1(3008,'Path length (cm) DetEast',200,480.50,480.81,0.)
c call hbook1(3009,'th_rotation (mrad) DetEast',200,-5.,5.0,0.)
c call hbook1(3010,'Path length (cm) DetWest',200,114.58,114.7,0.)
c call hbook1(3010,'Path length (cm) DetWest',200,480.50,480.81,0.)
c call hbook1(3011,'th_rotation (mrad) DetWest',200,-5.,5.0,0.)
call hbook1(3012,'E (meV) DetEast',2000,0.,40.,0.)
call hbook1(3013,'E (meV) DetWest',2000,0.,40.,0.)
c call hbook1(3014,'E (meV) scat DetWE',2000,0.,40.,0.)
c call hbook1(3020,'E_loss (meV) scat DetWE',2000,-0.4,0.4,0.)
c call hidopt(3020,'logy')
call hbook1(3021,'E_loss (meV) Allscat DetE',500,-0.5,2.0,0.)
call hbook1(3022,'E_loss (meV) Allscat DetW',500,-0.5,2.0,0.)
call hbook1(3023,'E_loss (meV) Tarscat DetE',1000,-0.5,0.5,0.)
call hbook1(3024,'E_loss (meV) Tarscat DetW',1000,-0.5,0.5,0.)
c
c adding more spectra below for diagnostics bec
call hbook2(8010,'sig tot (b/atom) vs. q (1/Ang)',50,0.0,0.2,
& 50,0.0,1.0,0.)
call hbook1(4010,'West sig tot (b/atom)',200,0.,1.,0.)
c call hbook1(4011,'West interaction len (cm)',200,0.,2000.,0.)
call hbook1(4020,'East sig tot (b/atom)',200,0.,1.,0.)
c call hbook1(4021,'East interaction len (cm)',200,0.,2000.,0.)
call hbook1(4012,'theta scattered (deg) U',200,0.,180.,0.)
c call hbook1(4015,'theta-p (deg) U',200,0.,180.,0.)
c call hbook1(4016,'phi-p (deg) U',200,0.,360.,0.)
call hbook1(4017,'th-0 (deg)',200,0.,10.,0.)
c call hbook1(4018,'phi-0 (deg) ',200,0.,360.,0.)
call hbook1(5012,'theta scattered (deg) D',200,0.,180.,0.)
c call hbook1(5015,'theta-p (deg) D',200,0.,180.,0.)
c call hbook1(5016,'phi-p (deg) D',200,0.,360.,0.)
call hbook1(5017,'theta-x (deg) U',200,-2.,2.,0.)
call hbook1(5018,'theta-y (deg) U',200,-2.,2.,0.)
call hbook1(5019,'theta-x (deg) D',200,-2.,2.,0.)
call hbook1(5020,'theta-y (deg) D',200,-2.,2.,0.)
c call hbook1(5021,'phi (deg) af-guide',200,0.,360.,0.)
c call hbook1(5022,'phi (deg) af-target 1',200,0.,360.,0.)
c call hbook1(5023,'phi (deg) af-target 2',200,0.,360.,0.)
call hbook1(6001,'Path length (cm) enter Up WE',200,213.0,
& 216.0,0.)
c call hbook1(6003,'Path length (cm) Det scat WE',300,480.50,
c & 480.90,0.)
c call hbook1(6004,'th_rotation (mrad) Det scat WE',400,-800.,
c & 0.0,0.)
call hbook1(6005,'PathLenTar (cm) Det scat West',300,
& tarlength-0.1,tarlength+0.5,0.)
call hbook1(6007,'PathLenTar (cm) Det scat East',300,
& tarlength-0.1,tarlength+0.5,0.)
c call hbook1(6009,'Path length (cm) Det no scat WE',300,480.50,
c & 480.90,0.)
c call hbook1(6010,'th_rotation (mrad) Det no scat WE',400,-800.,
c & 0.0,0.)
call hbook1(6011,'PathLenTar (cm) Det no scat East',300,
& tarlength-0.1,tarlength+0.5,0.)
call hbook1(6013,'PathLenTar (cm) Det no scat West',300,
& tarlength-0.1,tarlength+0.5,0.)
call hbook1(6015,'Crossover PathLenTar (cm) Det scat West',
& 300,tarlength-0.1,tarlength+0.5,0.)
call hbook1(6017,'Crossover PathLenTar (cm) Det scat East',
& 300,tarlength-0.1,tarlength+0.5,0.)
call hbook1(6035,'Crossover PathLenTar (cm) DetEast',
& 300,tarlength-0.1,tarlength+0.5,0.)
call hbook1(6037,'Crossover PathLenTar (cm) DetWest',
& 300,tarlength-0.1,tarlength+0.5,0.)
c
c Bz gradient and Pi-coil adjustments
c
if (bzgradnum.ne.0) then
c
if (picoilflg) then ! Bz gradient, Picoil ON
c call hbook1(3009,'th_rotation (mrad) DetEast',300,-800.,10.0,0.)
c call hbook1(3011,'th_rotation (mrad) DetWest',300,-800.,10.0,0.)
call hbook1(6002,'th_rotation (mrad) enter Up WE',400,-100.,
& 100.,0.)
call hbook1(6006,'th_rotation (mrad) Det scat West',400,-800.,
& 0.0,0.)
call hbook1(6008,'th_rotation (mrad) Det scat East',400,-800.,
& 0.0,0.)
call hbook1(6012,'th_rotation (mrad) Det no scat East',400,-800.,
& 0.0,0.)
call hbook1(6014,'th_rotation (mrad) Det no scat West',
& 400,-800.,0.,0.)
call hbook1(6016,'Crossover th_rotation (mrad) Det scat West',
& 400,-800.,0.0,0.)
call hbook1(6018,'Crossover th_rotation (mrad) Det scat East',
& 400,-800.,0.0,0.)
call hbook1(6036,'Crossover th_rotation (mrad) DetEast',
& 400,-800.,0.0,0.)
call hbook1(6038,'Crossover th_rotation (mrad) DetWest',
& 400,-800.,0.0,0.)
c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40.,
c & 300,-800.,0.,0.)
c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40.,
c & 300,-800.,0.,0.)
c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40.,
c & 300,-800.,0.,0.)
c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300,
c & 480.50,480.90,300,-800.,0.,0.)
c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40.,
c & 300,-800.,0.,0.)
c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300,
c & 480.50,480.90,300,-800.,0.,0.)
c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West',
c & 800,0.,40.,300,-800.,0.,0.)
c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West',
c & 300,480.50,480.90,300,-800.,0.,0.)
c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East',
c & 800,0.,40.,300,-800.,0.,0.)
c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East',
c & 300,480.50,480.90,300,-800.,0,0.)
else ! Bz gradient, Picoil OFF
c call hbook1(3009,'th_rotation (mrad) DetEast',300,-10.,80.0,0.)
c call hbook1(3011,'th_rotation (mrad) DetWest',300,-10.,80.0,0.)
call hbook1(6002,'th_rotation (mrad) enter Up WE',400,0.,
& 60.,0.)
call hbook1(6006,'th_rotation (mrad) Det scat West',400,0.,
& 60.0,0.)
call hbook1(6008,'th_rotation (mrad) Det scat East',400,0.,
& 60.0,0.)
call hbook1(6012,'th_rotation (mrad) Det no scat East',400,0.,
& 60.0,0.)
call hbook1(6014,'th_rotation (mrad) Det no scat West',
& 400,0.,60.,0.)
call hbook1(6016,'Crossover th_rotation (mrad) Det scat West',
& 400,0.,60.0,0.)
call hbook1(6018,'Crossover th_rotation (mrad) Det scat East',
& 400,0.,60.0,0.)
call hbook1(6036,'Crossover th_rotation (mrad) DetEast',
& 400,0.,60.0,0.)
call hbook1(6038,'Crossover th_rotation (mrad) DetWest',
& 400,0.,60.0,0.)
c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40.,
c & 300,-800.,0.,0.)
c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40.,
c & 300,-800.,0.,0.)
c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40.,
c & 300,0.,60.,0.)
c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300,
c & 480.50,480.90,300,0.,60.,0.)
c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40.,
c & 300,0.,60.,0.)
c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300,
c & 480.50,480.90,300,0.,60.,0.)
c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West',
c & 800,0.,40.,300,0.,60.,0.)
c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West',
c & 300,480.50,480.90,300,0.,60.,0.)
c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East',
c & 800,0.,40.,300,0.,60.,0.)
c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East',
c & 300,480.50,480.90,300,0.,60,0.)
end if
else ! No Bz gradient
If (picoilflg) then ! NO Bz gradient, Picoil ON
c
call hbook1(3009,'th_rotation (mrad) DetEast',300,-10.,10.0,0.)
call hbook1(3011,'th_rotation (mrad) DetWest',300,-10.,10.0,0.)
call hbook1(6002,'th_rotation (mrad) enter Up WE',300,-10.,
& 20.,0.)
call hbook1(6006,'th_rotation (mrad) Det scat West',300,-10.,
& 10.0,0.)
call hbook1(6008,'th_rotation (mrad) Det scat East',300,-10.,
& 10.0,0.)
call hbook1(6012,'th_rotation (mrad) Det no scat East',300,-10.,
& 10.0,0.)
call hbook1(6014,'th_rotation (mrad) Det no scat West',
& 300,-10.,10.,0.)
call hbook1(6016,'Crossover th_rotation (mrad) Det scat West',
& 300,-10.,10.0,0.)
call hbook1(6018,'Crossover th_rotation (mrad) Det scat East',
& 300,-10.,10.0,0.)
call hbook1(6036,'Crossover th_rotation (mrad) DetEast',
& 300,-10.,10.0,0.)
call hbook1(6038,'Crossover th_rotation (mrad) DetWest',
& 300,-10.,10.0,0.)
c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40.,
c & 300,-5.,50.,0.)
c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40.,
c & 300,-5.,50.,0.)
c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40.,
c & 300,-10.,10.,0.)
c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300,
c & 480.50,480.90,300,-10.,10.,0.)
c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40.,
c & 300,-10.,10.,0.)
c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300,
c & 480.50,480.90,300,-10.,10.,0.)
c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West',
c & 800,0.,40.,300,-10.,10.,0.)
c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West',
c & 300,480.50,480.90,300,-10.,10.,0.)
c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East',
c & 800,0.,40.,300,-10.,10.,0.)
c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East',
c & 300,480.50,480.90,300,-10.,10,0.)
else ! NO Bz gradient, Picoil OFF
call hbook1(3009,'th_rotation (mrad) DetEast',300,-5.,50.0,0.)
call hbook1(3011,'th_rotation (mrad) DetWest',300,-5.,50.0,0.)
call hbook1(6002,'th_rotation (mrad) enter Up WE',300,-5.,
& 50.,0.)
call hbook1(6006,'th_rotation (mrad) Det scat West',300,-5.,
& 50.0,0.)
call hbook1(6008,'th_rotation (mrad) Det scat East',300,-5.,
& 50.0,0.)
call hbook1(6012,'th_rotation (mrad) Det no scat East',300,-5.,
& 50.0,0.)
call hbook1(6014,'th_rotation (mrad) Det no scat West',
& 300,-5.,50.,0.)
call hbook1(6016,'Crossover th_rotation (mrad) Det scat West',
& 300,-5.,50.0,0.)
call hbook1(6018,'Crossover th_rotation (mrad) Det scat East',
& 300,-5.,50.0,0.)
call hbook1(6036,'Crossover th_rotation (mrad) DetEast',
& 300,-5.,50.0,0.)
call hbook1(6038,'Crossover th_rotation (mrad) DetWest',
& 300,-5.,50.0,0.)
c call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40.,
c & 300,-5.,50.,0.)
c call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40.,
c & 300,-5.,50.,0.)
c call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40.,
c & 300,-5.,50.,0.)
c call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300,
c & 480.50,480.90,300,-5.,50.,0.)
c call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40.,
c & 300,-5.,50.,0.)
c call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300,
c & 480.50,480.90,300,-5.,50.,0.)
c call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West',
c & 800,0.,40.,300,-5.,50.,0.)
c call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West',
c & 300,480.50,480.90,300,-5.,50.,0.)
c call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East',
c & 800,0.,40.,300,-5.,50.,0.)
c call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East',
c & 300,480.50,480.90,300,-5.,50,0.)
end if
end if
call hbook1(7001,'dE / 1st target',400,-1.,1.,0.)
call hbook1(7002,'dE / 2nd target',400,-1.,1.,0.)
call hbook1(7005,' Wgt Factor All',200,0.,2.,0.)
call hbook1(7006,' Wgt Factor Any Scat',200,0.,1.1,0.)
call hbook1(7007,' Wgt Factor Tar Scat',200,0.,1.1,0.)
c call hbook2(8011,'Thx vs. E (meV) scat DetB',100,0.,40.,
c & 200,-0.1,0.1,0.)
c call hbook2(8012,'Thx vs. E-loss (meV) scat DetB',100,
c & -0.4,0.4,200,-0.1,0.1,0.)
call hbook2(9999,'PA vs.WaveLength (Ang)',160,0.,16.,102,
& 0.,1.02,0.) !-KG
c call hbook2(9998,'omega (1/psec) vs. q (1/Ang)',100,0.,0.2,100,
c & -2.0,2.0,0.) !-KG
c
c------------Input file -------------------------
c open(39,file='inp.dat',status='old')
c read(39,*) gl, gin, T
c------------------------------------------------
r2d=180./pi
c------------
counter = 0
Event = 0
cin=0
cot=0
Nenscat1=0 !bec
Ntar1=0 !bec
Nthscat1=0
Nenscat2=0 !bec
Ntar2=0 !bec
Nthscat2=0
Nsdet=0 !bec
NEdet=0 !bec
nrtsraw=0;nltsraw=0;nrgsraw=0;nlgsraw=0
NTargetout=0; Ntar1out=0; Ntar2out=0 !bec
Nthdet=0
NPSMout=0 ; Nincoilout=0 ; NOutcoilout = 0 !bec
lsrot=0;lsrotsq=0;lsrottsc=0;lsrottscsq=0;lsrotgsc=0;lsrotgscsq=0
lsde=0;lsdesq=0;rsrot=0;rsrotsq=0;rsrottsc=0;rsrottscsq=0;
rsrotavg=0;lsrotavg=0
nrnsc=0;nlnsc=0;rsrotnsc=0;lsrotnsc=0;rsrotnscsq=0;lsrotnscsq=0
rsrotgsc=0
rsrotgscsq=0;rsde=0;rsdesq=0;nr=0;nrts=0;nrgs=0;nl=0;nlts=0;nlgs=0
nrPA=0;nrPAsq=0;rPAavg=0;nlPA=0;nlPAsq=0;lPAavg=0 !rcm PAsin
uncrsrotavg=0;uncrPAavg=0;unclsrotavg=0;unclPAavg=0
nrnscPA=0;nrnscPAsq=0;rsnscPAavg=0;rsrotnscavg=0
nlnscPA=0;nrnscPAsq=0;lsnscPAavg=0;lsrotnscavg=0
uncrsnsc=0;uncrsnscPA=0;unclsnsc=0;unclsnscPA=0
nrtsPA=0;nrtsPAsq=0;rtsPAavg=0;rsrottscavg=0;uncrstsc=0
uncrstscPA=0;nltsPA=0;nltsPAsq=0;ltsPAavg=0;lsrottscavg=0
uncrstsc=0;uncrstscPA=0;unclstsc=0;unclstscPA=0
nrgsPA=0;nrgsPAsq=0;rgsPAavg=0;rsrotgscavg=0;uncrsgsc=0;
uncrsgscPA=0;nlgsPA=0;nlgsPAsq=0;lgsPAavg=0;lsrotgscavg=0
uncrsxscPA=0; unclsxscPA=0
unclsgsc=0;unclsgscPA=0;nrallPA=0;nrallPAsq=0;nlallPA=0;rallavg=0
uncrscavg=0;uncrscPAavg=0;unclscavg=0;unclscPAavg=0
rallPAavg=0;lallavg=0;lallPAavg=0;nrscPA=0;rscPAavg=0;rscavg=0
nlscPA=0;lscPAavg=0;lscavg=0;nlallPAsq=0;uncrallavg=0;
uncrallPAavg=0;unclallavg=0;unclallPAavg=0;nlscPAsq=0;nrscPAsq=0
nder=0;ndel=0 ; rpl=0.0; lpl=0.0;rplns=0.0;lplns=0.0
rplsq=0.0;lplsq=0.0;rplnssq=0.0;lplnssq=0.0
Nalscat1=0;Nalscat2=0;Nalout1=0;Nalout2=0; Nalin1=0; Nalin2=0 !rcm
tar2flg = .false. ! bec
Nairin1 = 0; Nairout1 = 0; Nairscat1 = 0
histflg = .false. ; histflg2 = .false.
numtar1flg = 0; numendif1 = 0
numtar2flg = 0; numendif2 = 0 ; numbeamsc = 0.0
Nemp1in=0; Nemp2in=0; Nemp1out=0; Nemp2out=0; NtarWEin=0
Nxo=0 !-KG cross overs
numE2W=0 ; numW2E=0 ;numbounce=0.0
lwde=0.0;rwde=0.0;lwdesq=0.0;rwdesq=0.0
nlall=0;nrall=0;nlsc=0;nrsc=0;ndelw=0;nderw=0
ngbx=0;ngbx0=0;ngbx1=0;ngbx2=0;ngbx3=0
ngby=0;ngby0=0;ngby1=0;ngby2=0;ngby3=0
xmax = 0.d0; ymax=0.d0
INPKFlux=0.d0;PSMAVGFlux=0.d0;PSMPKFlux=0.d0;INCOAVGFlux=0.d0
INCOPKFlux=0.d0;TARAVGFlux=0.d0;TARPKFlux=0.d0;OUTCOAVGFlux=0.d0
OUTCOPKFlux=0.d0;DETAVGFlux=0.d0;DETPKFlux=0.d0
nanflg=.false. ; scwgt = 1.0d0
c
do 101 i = 1, nevts
Event = 1
pltot = 0.0
pltottar = 0.0
thtot = 0.0
enflg1=.false.
thflg1=.false.
enflg2=.false.
thflg2=.false.
wgt = 1.0 ! assume non-scattered until see scattering in a target
wgt2 = 1.0 ! assume non-scattered until see scattering in a target
tar1flg = .false. ; tar2flg = .false. ; bounceflg = .false.
bounceflg2 = .false. ; bounceflg3 = .false.
scatflg1 = 0 ; scatflg2 = 0 ; beamscflg = .false.
c
c Generate initial neutron wavelength
c
c
if (lamflg) then
if (i.eq.1) then
write(*,*) ' Lambda from Distribution '
write(3,*) ' Lambda from Distribution '
endif
xlam = HRNDM1(99) !New initial spectrum after PSM -KG
call hf1(101,real(xlam),1.) !-KG
if(xlam.lt.1.) goto 100
if(xlam.gt.20) goto 100
else
if (i.eq.1) then
write(*,*) ' All Lambda = 5Ang '
write(3,*) ' All Lambda = 5Ang '
endif
xlam = 5.0d0
end if
c
c call hf1(4004,real(xlam),1.) !Initial lambda histogram -KG
c
cbec
c
c The energy calculated from the neutron's wavelength
c
en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2
enstart = en0
if (en0.lt.0.d0) then
write(6,*) ' en<0, en=1e-6'
goto 100
end if
call hf1(20,real(en0),1.)
if (en0.le.1.0) encnt =encnt+1
c
c Generates the initial trajectory of the neutron.
c
321 call ranlux(vect,7)
c if (vect(5).eq.0.50) write(*,*) 'vect(5)=0.50' !for wq.f -KG
c vect(6) determines initial E or W, vect(7) gives the exact postition -KG
c
c thindex = ng6index*xlam ! NG6 about 1.7mrad/ang
c
c if (thindex.lt.1E-7) thindex = 1E-7
c if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) then
c write(6,*) ' Arccos(theta_start)>1 '
c goto 321
c end if
c
c
c th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec
cc write(99,*) vect(1), th0
c
c if (th0.lt.1.E-8) th0=1.E-8
c
c ph0 = 2.0d0*pi*vect(2) ! azimuthal angle (radians)
c call hf1(4017,real(th0*180./3.14159),1.)
cc call hf1(4018,real(ph0*180./3.14159),1.)
c
if (i.eq.1) then
if (hegastar) write(6,*) ' !! He gas target run !! '
if (hegastar) write(3,*) ' !! He gas target run !! '
if (airtar) write(6,*) ' !! Air target run !! '
if (airtar) write(3,*) ' !! Air target run !! '
if (bzgradnum.eq.0) then
write(6,*) ' ! NO Bz gradient in target region !'
write(3,*) ' ! NO Bz gradient in target region !'
else if (bzgradnum.eq.1) then
write(6,*) ' !!!!Running with Bz gradient: ',bzfile1,' !!!!'
write(3,*) ' !!!!Running with Bz gradient: ',bzfile1,' !!!!'
else if (bzgradnum.eq.2) then
write(6,*) ' !!!!Running with Bz gradient: ',bzfile2,' !!!!'
write(3,*) ' !!!!Running with Bz gradient: ',bzfile2,' !!!!'
else if (bzgradnum.eq.3) then
write(6,*) ' !!!!Running with Bz gradient: ',bzfile3,' !!!!'
write(3,*) ' !!!!Running with Bz gradient: ',bzfile3,' !!!!'
else if (bzgradnum.eq.4) then
write(6,*) ' !!!!Running with Bz gradient: ',bzfile4,' !!!!'
write(3,*) ' !!!!Running with Bz gradient: ',bzfile4,' !!!!'
end if
if (picoilflg) then
write(6,*) ' ! pi-coil ON ! '
write(3,*) ' ! pi-coil ON ! '
else
write(6,*) ' ! pi-coil OFF ! '
write(3,*) ' ! pi-coil OFF ! '
end if
if (sptmflg) then
write(6,*) ' ! Septum Upstream ! '
write(3,*) ' ! Septum Upstream ! '
else
write(6,*) ' ! No Septum Upstream ! '
write(3,*) ' ! No Septum Upstream ! '
end if
write(6,*)
write(3,*)
end if
If (xydist.eq.0) then
ccccccccccccccccccccccc Flat xy-distribution No x-y gradient
if (i.eq.1) then
write(*,*) ' ! Running with Flat xy-distribution !'
write(3,*) ' ! Running with Flat xy-distribution !'
end if
xtest = xin*(2.*vect(3)-1.) !-KG
x0 = xtest
xstart = x0 !-KG
y0 = yin*(2.*vect(4)-1.)
elseif (xydist.eq.1) then
ccccccccc gradient
if (i.eq.1) then
write(*,*) ' !!! Running with Position (x) gradient !'
write(3,*) ' !!! Running with Position (x) gradient !'
end if
cc x0 = 3.0d0 - 2.d0*HRNDM1(151) !bec R to L gradient
cc
cc x0 = 2.d0*HRNDM1(151) - 3.0d0 !bec L to R gradient
cc this gives double gradient -- gradient across each half...
xtest = xin*(2.*vect(3)-1.) !-KG
x0 = xtest
if(xtest.gt.0.) then
x0 = HRNDM1(151) ! gives gradient only acrros half
cc x0 = 3.d0 - HRNDM1(151)
elseif(xtest.lt.0) then
x0 = HRNDM1(151) - 3.
ccc x0 = -HRNDM1(151)
endif
xstart = x0
ccccccccc
elseif (xydist.eq.2) then
cccccccccccc x-y from PSM image
if (i.eq.1) then
write(*,*) ' ! xy-distribution from psm.hbook !'
write(3,*) ' ! xy-distribution from psm.hbook !'
end if
call HRNDM2(200,xtmp,ytmp)
x0= dble(xtmp)
y0 = dble(ytmp)
xtest = x0
xstart = x0
if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100
if ((x0**2+y0**2).le.0.625**2)
& INPKflux = INPKflux + 1/pi/0.625**2/nevts
cccccccccccc
elseif (xydist.eq.3) then
ccccccccc Inverted Parabola x - y distributions
if (i.eq.1) then
write(6,*) ' !!! Inverted Parabola x-y dist!!'
write(3,*) ' !!! Inverted Parabola x-y dist!!'
end if
xtest = HRNDM1(161)
y0 = HRNDM1(162)
x0 = xtest
xstart = xtest
if (x0.gt.xmax) xmax = x0
if (y0.gt.ymax) ymax = y0
if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100
if (x0**2+y0**2.le.0.625**2)
& INPKflux = INPKflux + 1/pi/0.625**2/nevts
ccccccccc
end if
z0= 0.0d0
cc
cc Calculate the incident angles in the xz-plane and yz-plane
cc note: these are not "direction cosines"= angle between momentum and x (or y)
cc
c thx0 = ratan(dtan(th0)*dcos(ph0))
c thy0 = ratan(dtan(th0)*dsin(ph0))
CC************************************************************************
c************for nSpinSim3.1.3 choosing thx0 and thy0 evening withing their
c respective critical angles..
if (angdist.eq.0) then
c isotropic
if (i.eq.1) then
write(6,*) ' ! Isotropic angles!',ng6indexx,ng6indexy
write(3,*) ' ! Isotropic angles!',ng6indexx,ng6indexy
end if
thindex = ng6indexx*xlam ! NG6 about 1.7mrad/ang
ph0 = 2.0d0*pi*vect(2) ! azimuthal angle (radians)
if (thindex.lt.1E-7) thindex = 1E-7
if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) then
write(6,*) ' Arccos(theta_start)>1 '
goto 321
end if
th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec
thx0 = ratan(dtan(th0)*dcos(ph0))
thy0 = ratan(dtan(th0)*dsin(ph0))
elseif (angdist.eq.1) then
c
c flat
c
if (i.eq.1) then
write(6,*) ' ! Flat angular dist.!',ng6indexx,ng6indexy
write(3,*) ' ! Flat angular dist.!',ng6indexx,ng6indexy
end if
thx0 = ng6indexx*xlam*(vect(1)-0.5d0)*2.0d0
thy0 = ng6indexy*xlam*(vect(2)-0.5d0)*2.0d0
elseif (angdist.eq.2) then
c
c cosine dist
if (i.eq.1) then
write(6,*); write(3,*)
write(6,*) ' ! Cosine Angular Dist.!',ng6indexx,ng6indexy
write(3,*) ' ! Cosine Angular Dist.!',ng6indexx,ng6indexy
write(6,*); write(3,*)
end if
!!!!set maximum divergence in const.inc ng6indexx,y
thx0 = 2*ng6indexx*xlam/pi*dasin(2.0d0*vect(1)-1.d0)
thy0 = 2*ng6indexy*xlam/pi*dasin(2.0d0*vect(2)-1.d0)
elseif (angdist.eq.3) then
C********for nSPinSim3.1.4 choosing thx0 and thy0 from hbooks, interpolated
C******distributions from J. Cook simulations of NG6 and NGC
if (i.eq.1) then
write(6,*) ' ! Angular Dist from J. Cook !'
write(3,*) ' ! Agular Dist from J. Cook !'
end if
if(xlam.lt.2.or.xlam.gt.12.) goto 101
nlamx = int(xlam*10+1)+980
nlamy = int(xlam*10+1)+1080
call hrget(nlamx,'hbooks/ngcx.hbook','')
thx0 = HRNDM1(nlamx)
call hdelet(nlamx)
call hrget(nlamy,'hbooks/ngcy.hbook','')
thy0 = HRNDM1(nlamy)
call hdelet(nlamy)
C**************************************************************
CC****************************************************************
c
end if
thxinput = thx0
thyinput = thy0
call hf1(31,real(thxinput),1.)
call hf1(32,real(thyinput),1.)
NPSMout = NPSMout + 1 ! bec count exiting PSM
PSMAVGflux = PSMAVGflux + 1/xin/yin/4.0d0/nevts
if (x0**2+y0**2.le.0.625**2)
& PSMPKflux = PSMPKflux + 1/pi/0.625**2/nevts
capflpsm = capflpsm + sqrt(25.0)/rsqrt(en0)
call hf1(1,real(x0),1.)
call hf1(10,real(y0),1.)
ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area
ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than
ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01
cccccccc =1.2869E8*3.1356E5*1.01=4.076E12
call hf2(2220,real(x0),real(y0),4.076E12/nevts)
if (dabs(y0).le.0.5) call hf1(441,real(x0),4.076E12/nevts*0.02)
if (dabs(x0).le.0.5) call hf1(442,real(y0),4.076E12/nevts*0.02)
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
2341 format(' EastDn ',1x,f9.2,1x,f9.2,1x,f9.2)
2342 format(' WestUp ',1x,f9.2,1x,f9.2,1x,f9.2)
********************************************************************
* Neutrons passing through the gap (3.0 cm) after the PSM before input guides
* this is being added in to old code, so x0, y0, z0, etc just redefined as
* result of this 3 cm transport
*
if((x0.ge.tairx1.and.x0.le.tairx2).and.
& (y0.ge.tairy1.and.y0.le.tairy2)) Nairin1 = Nairin1 + 1
c
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,tairz1, en1, x1, y1, z1,
& thx1,thy1,airflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) Nairscat1 = Nairscat1+1 ! count scatters with energy change bec
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg AIR 1',airflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
c write(6,*) x1,airflg1,beamscflg
if (en1.lt.0.d0) then
write(6,*) 'ag1 en<1, en=1e-6'
en1 = 0.0000001
end if
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot ! mrad
if (pltot.gt.1300) write(6,*) '0th air gap',pltot
Nairout1 = Nairout1 + 1
if (dabs(thtot).gt.1.0e9) then
write(6,*) 'Air thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through FIRST Al Window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
if((x0.ge.talx1.and.x0.le.talx2).and.
& (y0.ge.taly1.and.y0.le.taly2)) Nalin1 = Nalin1 + 1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) Nalscat1 = Nalscat1+1 ! count scatters with energy change bec
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/airflg AL 1',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
if (en1.lt.0.d0) then
write(6,*) 'Al1 en<1, en=1e-6'
en1 = 0.0000001
end if
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al 1 ',pltot
Nalout1 = Nalout1 + 1
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
c Neutrons passing through the guide
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
beamflg = 1
if (sptmflg) then
if ((x0.ge.gx1.and.x0.le.gx2)
& .and.(y0.ge.gy1.and.y0.lt.gy2)) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,gx1,gx2,gy1,gy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
if((x0*x1.lt.0d0).and.x1.gt.-1d5)
& write(6,*) ' Guide Crossover ',x0,x1 !check for crossovers
else if ((x0.ge.gx3.and.x0.le.gx4).and.
& (y0.ge.gy1.and.y0.lt.gy2)) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,gx3,gx4,gy1,gy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
if((x0*x1.lt.0d0).and.x1.gt.-1d5)
& write(6,*) ' Guide Crossover ',x0,x1 !check for crossovers
else
goto 100
end if
else
if((x0.ge.gx1.and.x0.le.gx4).and.(y0.ge.gy1.and.y0.le.gy2)) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,gx1,gx4,gy1,gy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
else
goto 100
end if
end if
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
call hf1(2,real(x1),1.)
if (nbx.eq.0) ngbx0=ngbx0+1
if (nbx.eq.1) ngbx1=ngbx1+1
if (nbx.eq.2) ngbx2=ngbx2+1
if (nbx.eq.3) ngbx3=ngbx3+1
if (nby.eq.0) ngby0=ngby0+1
if (nby.eq.1) ngby1=ngby1+1
if (nby.eq.2) ngby2=ngby2+1
if (nby.eq.3) ngby3=ngby3+1
if (nbx.gt.3) write(6,*) ' 4+ Guide xBounces ',nbx
if (nby.gt.3) write(6,*) ' 4+ Guide yBounces ',nby
gsc = gsc + scflag
if(scflag.gt.5) write(6,*)' 5+ Guide Scatters ',scflag
if (scflag.gt.0) beamscflg = .true.
if ( ((en0.eq.en1).and.(scflag.gt.0)).or.
& ((en0.ne.en1).and.(scflag.eq.0)) )
& write(6,*) ' Edif/scatflg GUIDE',scflag,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
c This path length is now OK for reflected neutrons -KG
c pl_del = gl/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2)))
c pltot = pl_del + pltot
c thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
pltot = pl_del + pltot
thtot = thtot + rotang
if (pltot.gt.1300) write(6,*) 'guide ',pltot
Nguideout = Nguideout+1
if(x1.gt.gx1.and.x1.lt.gx2) NguideoutE = NguideoutE+wgt
if(x1.gt.gx3.and.x1.lt.gx4) NguideoutW = NguideoutW+wgt
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' guide thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true.
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
********************************************************************
* Neutrons passing through SECOND Al Window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
if((x0.ge.talx1.and.x0.le.talx2).and.
& (y0.ge.taly1.and.y0.le.taly2)) Nalin2 = Nalin2 + 1 ! rcm
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz2, en1, x1, y1, z1,
& thx1,thy1,alflg2)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg2) Nalscat2 = Nalscat2+1 !rcm
if (alflg2) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg2)).or.
& ((en0.ne.en1).and.(.not.(alflg2))) )
& write(6,*) ' Edif/alflg Al 2',alflg2,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
if(en1.lt.0.d0) then
write(6,*) 'All en<1, en=1e-6'
en1 = 0.0000001
end if
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al 2 ',pltot
Nalout2 = Nalout2 + 1
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
*************************************************************************
* Neutrons passing through Al wire input coil entrance
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 3',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR 1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al input ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
************************************************************************
* Neutrons passing through the input coil
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
beamflg = 2
if (sptmflg) then
if ((x0.ge.ix1.and.x0.le.ix2).and.
& (y0.ge.iy1.and.y0.le.iy2)) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,ix1,ix2,iy1,iy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
if((x0*x1.lt.0d0).and.x1.gt.-1d5)
& write(6,*) ' Incoil Crossover ',x0,x1 !check for crossovers
else if ((x0.ge.ix3.and.x0.le.ix4)
& .and.(y0.ge.iy1.and.y0.le.iy2)) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,ix3,ix4,iy1,iy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
if((x0*x1.lt.0d0).and.x1.gt.-1d5)
& write(6,*) ' Incoil Crossover ',x0,x1 !check for crossovers
else
goto 100
end if
else
if((x0.ge.ix1.and.x0.le.ix4).and.(y0.ge.iy1.and.y0.le.iy2)) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,ix1,ix4,iy1,iy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
else
goto 100
end if
end if
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
Nincoilout = Nincoilout + 1 ! bec counttransmission through guide and Input coil
if(x1.gt.ix1.and.x1.lt.ix2) NincoiloutE = NincoiloutE+wgt
if(x1.gt.ix3.and.x1.lt.ix4) NincoiloutW = NincoiloutW+wgt
call hf1(3,real(x1),1.)
call hf2(1111,real(x1),real(y1),1.)
if (nbx.eq.0) nibx0=nibx0+1
if (nbx.eq.1) nibx1=nibx1+1
if (nbx.eq.2) nibx2=nibx2+1
if (nbx.eq.3) nibx3=nibx3+1
if (nby.eq.0) niby0=niby0+1
if (nby.eq.1) niby1=niby1+1
if (nby.eq.2) niby2=niby2+1
if (nby.eq.3) niby3=niby3+1
if (nbx.gt.3) write(6,*) ' 4+ Incoil xBounces ',nbx
if (nby.gt.3) write(6,*) ' 4+ Incoil yBounces ',nby
isc = isc + scflag
if(scflag.gt.5) write(6,*)' 5+ Incoil Scatters ',scflag
if (scflag.gt.0) beamscflg = .true.
if ( ((en0.eq.en1).and.(scflag.gt.0)).or.
& ((en0.ne.en1).and.(scflag.eq.0)) )
& write(6,*) ' Edif/scatflg INPUT',scflag,en0,en1
c if (en0.ne.en1) write(6,*) ' edif INCOIL ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
c This path length is now OK for reflected neutrons -KG
cc pl_del = iz0/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2))) bec2011 no longer needed with new routines!
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'incoil ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' incoil thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
*
if ((bounceflg).and.((thx1.eq.-1.0*thx0).or.
& (thy1.eq.-1.0*thy0))) bounceflg2 = .true.
if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true.
capflinc = capflinc + sqrt(25.0)/rsqrt(en1)
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through Al wire input coil exit
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 3',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al inexit ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al input ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
********************************************************************
* Neutrons passing through the gap (1st 1/2 of 10.0 cm) between the input
* coil and the target.
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,0.5*tairz2, en1, x1, y1, z1,
& thx1,thy1,airflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg Al 3',airflg1,en0,en1
c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers
c if (en0.ne.en1) write(6,*) ' edif in/tar ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
if(thx1.ge.0..and.thy1.ge.0.) then
phi = ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.ge.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.lt.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.ge.0..and.thy1.lt.0.) then
phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1))
endif
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'gap 1 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
*
********************************************************************
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
******************* INput Coil Exit Flux Measurement***************
if (sptmflg) then
INCOAVGflux = INCOAVGflux + 1/((ix2-ix1)+(ix4-ix3))/(iy2-iy1)/nevts
if ((x1+(ix2+ix1)/2.)**2+y1**2.le.0.625**2)
& INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts
else
INCOAVGflux = INCOAVGflux + 1/(ix4-ix1)/(iy2-iy1)/nevts
if (x1**2+y1**2.le.0.625**2)
& INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts
end if
ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area
ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than
ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01
cccccccc =1.2869E8*3.1356E5*1.01=4.076E12
call hf2(2221,real(x1),real(y1),4.076E12/nevts)
if (dabs(y1).le.0.5) call hf1(443,real(x1),4.076E12/nevts*0.02)
if (dabs(x1).le.0.5) call hf1(444,real(y1),4.076E12/nevts*0.02)
********************************************************************
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through the gap (2nd 1/2 of 10.0 cm) between the input
* coil and the target.
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,0.5*tairz2, en1, x1, y1, z1,
& thx1,thy1,airflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg AIR in/tar',airflg1,en0,en1
c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers
c if (en0.ne.en1) write(6,*) ' edif in/tar2 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call hf1(4,real(x1),1.)
call hf1(5017,real(180./3.14159*thx1),1.)
call hf1(5018,real(180./3.14159*thy1),1.)
if(thx1.ge.0..and.thy1.ge.0.) then
phi = ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.ge.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.lt.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.ge.0..and.thy1.lt.0.) then
phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1))
endif
c call hf1(5021,real(180./3.14159*phi),1.)
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
call hf1(6001,real(pltot),1.)
call hf1(6002,real(thtot),1.)
if (pltot.gt.1300) write(6,*) 'gap 1 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
*
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through THIRD Al window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz3, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 3',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al 3 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
if (en1.lt.0.d0) then
write(6,*) 'Al2 en<1, en=1e-6'
en1 = 0.0000001
end if
if (alflg1) beamscflg = .true.
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al 3 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
********************************************************************
NtarWEin = NtarWEin + 1
xtar1 = x1
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
c print *, " first half target"
* Neutrons passing through the 1st half of the target (41.6 cm).
C Note: below checks only where the neutron is in x direction.
C Once it enters target target.f checks the y-direction too, so might
C have neutrons that get into target.f that don't run target.f really,
C since the y-values are out of range, hence sig or thscat values get
C nothing put in them
*
* Defining x<0 to imply West which in main_WUp.exe has upstream target
C
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
ztar1=z0 ! begining of 1st target where Bz-gradient starts
c
enflg1 = .false.
thflg1 = .false.
empflg1 = 0
if( (x0.ge.tarx1.and.x0.le.tarx2).and.
& (y0.ge.tary1.and.y0.le.tary2) ) then ! East=neg. x, empty UPstream
Nemp1in=Nemp1in+wgt
call empty(en0, x0, y0, z0, thx0, thy0, tarx1, tarx2,
& en1, x1, y1, z1, thx1, thy1,empflg1,ztar1,pl_del,rotang)
if (en1.lt.0.d0) goto 100
if (empflg1.gt.0) Nempty1 = Nempty1+1.0 ! not trajectory dependent -- no wgt
if (empflg1.gt.0) beamscflg = .true.
if ((x1.ge.tarx1.and.x1.le.tarx2).and.
& (y1.ge.tary1.and.y1.le.tary2))
& Nemp1out = Nemp1out+wgt ! out of empty west side only
elseif( (x0.ge.tarx3.and.x0.le.tarx4).and.
& (y0.ge.tary1.and.y0.le.tary2) ) then ! West=pos. x, Target UP
Ntar1=Ntar1+wgt ! entering target, before wgt change
sig = 0.0 ; thscat = 0.0
thp = 0.0; php = 0.0
c Given current energy determine q from S(q) and delE from omega(q)
c Doing this only if in west side where target is located
c
k0 = rsqrt(2*m_n*en0)/h_par
loop = 0
qmin = 0
qmax = int(k0*1e3)
if(qmax.lt.100) then
hmax = qmax
else
hmax=100
endif
call hcopyr(699,302,"",qmin,hmax,0.,0.,"")
cbec10 call hcopyr(700,302,"",qmin,hmax,0.,0.,"")
501 q=HRNDM1(302)
if(q.eq.0.) goto 501
hws1 = q*1e3+700
call hrget(hws1,'hbooks/sqw.hbook','')
wp1=HRNDM1(hws1)
call hdelet(hws1)
ep1 = en0 - (h_par*wp1)
if (ep1.lt.0.d0) then
write(6,*) ' 1st tar, ep1<0 ',ep1
goto 501
endif
kp1 = rsqrt(2*m_n*ep1)/h_par
th_test1 = (k0**2 + kp1**2 - q**2)/(2.*k0*kp1)
loop = loop + 1
if(ep1.lt.0.or.abs(th_test1).gt.1.0) goto 501
call hdelet(302)
call hf1(305,real(q),1.)
call hf1(306,real(wp1),1.)
call hf1(307,real(en0),1.)
call hf1(308,real(ep1),1.)
call hf1(309,real(k0),1.)
call hf1(310,real(kp1),1.)
call target(q, en0, x0, y0, z0, thx0, thy0, tarx3, tarx4,
& en1, x1, y1, z1, thx1, thy1, sig, thscat,
& scatflg1,nrfl1,ztar1,pl_del,rotang,scwgt) ! bec
if (en1.lt.-1000.d0) goto 100
if(en1.ne.en0) enflg1=.true.
If ((.not.hegastar).and.(.not.airtar)) wgt2=wgtA/sig*scwgt
c If ((.not.hegastar).and.(.not.airtar)) wgt2=1.0d0
cbec fall2011 if((en0.gt.-100).and.(en1.gt.-100).and.(en0.ne.en1)) then ! scattered
if (scatflg1.gt.0) then ! scattered in target.f/scatt.f
tar1flg= .true.
c write(81,*) en0,en1,en0-en1
c determine the scattering weighting factor ratio of integral of S(q,w)
c over limited q range (q<0.1 1/Ang) to full
c integral (= energy dependent cross section from Sommers)
If ((.not.hegastar).and.(.not.airtar)) wgt=wgtA/sig*scwgt
c If ((.not.hegastar).and.(.not.airtar)) wgt=1.0d0
call hf1(7001,real(en0-en1),wgt)
endif
if ( (enflg1.and.(scatflg1.eq.0)).or.
& (.not.(enflg1).and.(scatflg1.gt.0)) )
& write(6,*) ' enflg/scatflg ',scatflg1,en0,en1
if (wgt.lt.0) write(6,*) ' 1st tar wgt<0 ', wgt, sig, en1
if (enflg1) Nenscat1 = Nenscat1+1.0 ! not trajectory dependent -- no wgt
c if (thflg1) Nthscat1 = Nthscat1+1.0 ! not trajectory dependent -- no wgt
c capture diagnostics for target interaction stuff bec
call hf1(4010,real(sig),1.0) !not trajectory dependent -- no wgt
call hf2(8010,real(q),real(sig),1.0) !temporary!not trajectory dependent -- no wgt
call hf1(4012,real(180./3.14159*thscat),wgt)
call hf1(4015,real(180./3.14159*thp),wgt)
c call hf1(4016,real(180./3.14159*php),wgt)
if ((x1.ge.tarx3.and.x1.le.tarx4).and.
& (y1.ge.tary1.and.y1.le.tary2)) Ntar1out = Ntar1out+wgt ! out of target west side only
else
x1 = -1.d6
y1 = -1.d6
z1 = -1.d6
thx1 = -1.d6
thy1 = -1.d6
en1 = -1.d6
endif
*
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
*
entar1 = en1
call hf1(1001,real(x1),wgt)
call hf1(1002,real(y1),wgt)
c call hf1(1003,real(thx1),wgt)
c call hf1(1004,real(thy1),wgt)
cccccccccccccc target rotation
c write(6,*) ' pl_del ', pl_del, ' rotang ',rotang
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
cbec fall2011 if(en0.eq.en1) then
if (scatflg1.eq.0) then
call hf1(1005,real(en1),wgt)
else
call hf1(1012,real(en1),wgt)
endif
cccccccccc end target rotation
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' tar1 thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
c write(*,*) 'tar 1 ends at z = ',z1 !-KG
*
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through FOURTH Al window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz4, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 4',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al 4 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
if (pltot.gt.1300) write(6,*) 'Al 4 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al tar1 thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
********************************************************************
ctemp if (wgt.ne.1.0) write(17,*) ' tar1 out wgt ',wgt,tar1flg,tar2flg
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through Al pi-coil wire support
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz5, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al pi1',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al pi1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
if (pltot.gt.1300) write(6,*) 'Al Pi ',pltot,x0,y0,z0,x1,y1,z1
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
***********************************************************************
c print *," Gap btw targets"
* Neutrons passing thruogh the gap (pi-coil = 7.3 cm) between the
* two tagets = 8.0cm.
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
cccc 1st half of air gap
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,0.5*tairz3, en1, x1, y1, z1,
& thx1,thy1,airflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg AIR tar1',airflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR tar1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is starting point of gradient field if being used
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
call hf1(5019,real(180./3.14159*thx1),wgt)
call hf1(5020,real(180./3.14159*thy1),wgt)
if(thx1.ge.0..and.thy1.ge.0.) then
phi = ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.ge.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.lt.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.ge.0..and.thy1.lt.0.) then
phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1))
endif
c call hf1(5022,real(180./3.14159*phi),wgt)
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
*********************************************************************
* Neutrons passing through Al Wire Pi Coil support
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz5, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al pi2',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al pi2 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
if (pltot.gt.1300) write(6,*) 'Al Pi ',pltot,x0,y0,z0,x1,y1,z1
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
***********************************************************************
c do NOT allow neutrons from one target cell into opposite side cell
if( ((x0.ge.tarx1).and.(x0.le.tarx2)).and.((x1.lt.tarx1).or.
& (x1.gt.tarx2)) )
& goto 100
if( ((x0.ge.tarx3).and.(x0.le.tarx4)).and.((x1.lt.tarx3).or.
& (x1.gt.tarx4)) )
& goto 100
c*************PI Coil reversal of rotation angle
c using 1/2 gap length -- rotate to halfway point, reverse, rotate the other half
c reverse thtot in center of gap and rotate the other half of gap length
if (picoilflg) then ! pi-coil ON
thtot = - thtot
end if
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,0.5*tairz3, en1, x1, y1, z1,
& thx1,thy1,airflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg AIR tar2',airflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR tar2 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pltot ',pltot
write(6,*) ' en0 ',en0,', en ',en1
write(6,*) ' thxy01',thx0,thy0, thx1, thy1
end if
ctemp
cc if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through FIFTH Al window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz4, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 5',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al 5 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
call rotation(en0,x0,y0,z0,x1,y1,z1,ztar1,rotang) ! ztar1 is tarting point of gradient field if being used
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
if (pltot.gt.1300) write(6,*) 'Al 5 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
********************************************************************
xtar2 = x1
cctemp if (wgt.ne.1.0) write(17,*) ' tar2 in wgt ',wgt,tar1flg,tar2flg
ctemp
cc if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
c print *, " second half target"
* Neutrons passing thruogh the 2nd part of the target (41.6 cm).
*
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
enflg2 = .false.
thflg2 = .false.
empflg2 = 0
if( (x0.gt.tarx1.and.x0.lt.tarx2).and.
& (y0.ge.tary1.and.y0.le.tary2) ) then ! East = neg. x, Target down
Ntar2=Ntar2+wgt ! entering target, before wgt change
if (tar1flg) then ! if scat in t1, now entered both targets
Ntar12=Ntar12+wgt
endif
sig = 0.0 ; thscat = 0.0
thp = 0.0; php = 0.0
c
c Given current energy determine q from S(q) and delE from omega(q)
c Doing this only if in east side where target is located
c
k0 = rsqrt(2*m_n*en0)/h_par
loop = 0
qmin = 0
qmax = int(k0*1e3)
if(qmax.lt.100) then
hmax = qmax
else
hmax=100
endif
call hcopyr(699,302,"",qmin,hmax,0.,0.,"")
cbec10 call hcopyr(700,302,"",qmin,hmax,0.,0.,"")
502 q=HRNDM1(302)
if(q.eq.0.) goto 502
hws2 = q*1e3+700
call hrget(hws2,'hbooks/sqw.hbook','')
wp2=HRNDM1(hws2)
call hdelet(hws2)
ep2 = en0 - (h_par*wp2)
if (ep2.lt.0.d0) then
write(6,*) ' 2nd tar, ep2<0 ',ep2
goto 502
endif
kp2 = rsqrt(2*m_n*ep2)/h_par
th_test2 = (k0**2 + kp2**2 - q**2)/(2.*k0*kp2)
loop = loop + 1
if(ep2.lt.0.or.abs(th_test2).gt.1.0) goto 502
call hdelet(302)
call hf1(205,real(q),wgt)
call hf1(206,real(wp2),wgt)
call hf1(207,real(en0),wgt)
call hf1(208,real(ep2),wgt)
call hf1(209,real(k0),wgt)
call hf1(210,real(kp2),wgt)
call target(q, en0, x0, y0, z0, thx0, thy0, tarx1, tarx2,
& en1, x1, y1, z1, thx1, thy1, sig, thscat,
& scatflg2, nrfl2, ztar1, pl_del, rotang,scwgt) ! bec
if (en1.lt.-1000.d0) goto 100
if (en1.ne.en0) enflg2 = .true.
If ((.not.hegastar).and.(.not.airtar)) wgt2=wgtA/sig*scwgt
c If ((.not.hegastar).and.(.not.airtar)) wgt2=1.0d0
c
cbec fall2011 if(en0.gt.-100.and.en1.gt.-100.and.en0.ne.en1) then ! scattered
if (scatflg2.gt.0) then
tar2flg = .true.
c write(81,*) en0,en1,en0-en1
c determine the scattering weighting factor ratio of integral of S(q,w)
c over limited q range (q<0.1 1/Ang) to full
c integral (= energy dependent cross section from Sommers)
If ((.not.hegastar).and.(.not.airtar)) wgt=wgtA/sig*scwgt
c If ((.not.hegastar).and.(.not.airtar)) wgt=1.0d0
call hf1(7002,real(en0-en1),wgt)
endif
if ( (enflg2.and.(scatflg2.eq.0)).or.
& (.not.(enflg2).and.(scatflg2.gt.0)) )
& write(6,*) ' enflg/scatflg ',scatflg2,en0,en1
if (wgt.lt.0) write(6,*) ' 2nd tar wgt<0 ',wgt, sig, en1
c
if (enflg2) Nenscat2 = Nenscat2+1.0 ! not trajectory dependent -- no wgt
c if (thflg2) Nthscat2 = Nthscat2+1.0 ! not trajectory dependent -- no wgt
c capture diagnostics for target interaction stuff bec
call hf1(4020,real(sig),1.0)
c call hf1(4021,real(intlen),wgt)
call hf1(5012,real(180./3.14159*thscat),wgt)
call hf1(5015,real(180./3.14159*thp),wgt)
c call hf1(5016,real(180./3.14159*php),wgt)
call hf2(9998,real(q),real(wp1),wgt) ! Dispersion relation -KG
if ((x1.ge.tarx1.and.x1.le.tarx2).and.
& (y1.ge.tary1.and.y1.le.tary2)) Ntar2out = Ntar2out+wgt ! out of target west side only
elseif( (x0.gt.tarx3.and.x0.lt.tarx4).and.
& (y0.ge.tary1.and.y0.le.tary2) ) then ! West = pos. x, empty DOWN
Nemp2in=Nemp2in+wgt
call empty(en1, x0, y0, z0, thx0, thy0, tarx3, tarx4,
& en1, x1, y1, z1, thx1, thy1, empflg2,ztar1,pl_del,rotang)
if (en1.lt.0.d0) goto 100
if (empflg2.gt.0) Nempty2 = Nempty2+1.0 ! not trajectory dependent -- no wgt
if (empflg2.gt.0) beamscflg = .true.
if ((x1.ge.tarx3.and.x1.le.tarx4).and.
& (y1.ge.tary1.and.y1.le.tary2))
& Nemp2out = Nemp2out+wgt ! out of empty east side only
else
x1 = -1.d6
y1 = -1.d6
z1 = -1.d6
thx1 = -1.d6
thy1 = -1.d6
en1 = -1.d6
endif
*
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
*
entar2 = en1
call hf1(2001,real(x1),wgt)
call hf1(2002,real(y1),wgt)
ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area
ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than
ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01
cccccccc =1.2869E8*3.1356E5*1.01=4.076E12
call hf2(2222,real(x1),real(y1),4.076E12*wgt/nevts)
c call hf1(2003,real(thx1),wgt)
c call hf1(2004,real(thy1),wgt)
cccccccc target 2 rotation
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = thtot + rotang
cbec fall2011 if(en1.eq.en0) then
if (scatflg2.eq.0) then
call hf1(2005,real(en1),wgt)
else
call hf1(2012,real(en1),wgt)
endif
cccccccccccc end target 2 rotaton
if (pltot.gt.1300) write(6,*) 'target 2 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' tar2 thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
c write(*,*) 'tar 2 ends at z = ', z1
*
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
cc if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
ctemp if (wgt.ne.1.0) write(17,*) ' tar2 out wgt ',wgt,tar1flg,tar2flg
NTargetout = NTargetout + wgt
TARAVGflux = TARAVGflux +
& wgt/((tarx2-tarx1)+(tarx4-tarx3))/(tary2-tary1)/nevts
if ((x1-0.5*(tarx2-tarx1))**2+y1**2.le.0.625**2)
& TARPKflux = TARPKflux + wgt/pi/0.625**2/nevts
********************************************************************
********************************************************************
* Neutrons passing through SIXTH Al window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz3, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 6',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al 6 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
pltottar = pl_del + pltottar
thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al 6 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al tar2 thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
cc if (x1.lt.0) write(6,2341) z1,pltot,thtot
cc if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
***********************************************************************
* Neutrons passing thruogh the gap (13.0 cm) between the taget and
* the output coil.
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,tairz4, en1, x1, y1, z1,
& thx1,thy1,airflg1)
c call collimator(x0,x1,y0,y1) !boroflex septum to prevent crossovers
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg AIR tar/out',airflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR tar/out ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
if(thx1.ge.0..and.thy1.ge.0.) then
phi = ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.ge.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.lt.0..and.thy1.lt.0.) then
phi = pi + ratan(dtan(thy1)/dtan(thx1))
elseif(thx1.ge.0..and.thy1.lt.0.) then
phi = 2.*pi + ratan(dtan(thy1)/dtan(thx1))
endif
c call hf1(5023,real(180./3.14159*phi),1.)
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'gap 3 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
********************************************************************
ctemp
cc if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through Al Wire output coil entrance
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al outcoil in',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al out1 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al output ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through the output coil.
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
cccccccccccccc If you want to put outcoil guide at an angle ccccc
c if (i.eq.1) write(6,*) ' !!! outcoil ang offset, thy->2mrad '
c thy0 = thy0+0.002
cccccccccccccccccccccc
beamflg = 3
if((x0.ge.ox1.and.x0.le.ox4).and.(y0.ge.oy1.and.y0.le.oy2)) then
if (x0.ge.ox1.and.x0.le.ox2) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,ox1,ox2,oy1,oy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
if((x0*x1.lt.0d0).and.x1.gt.-1d5)
& write(6,*) ' Outcoil Crossover ',x0,x1 !check for crossovers
else if (x0.ge.ox3.and.x0.le.ox4) then
call beampipe(beamflg,en0,x0,y0,z0,
& thx0,thy0,ox3,ox4,oy1,oy2,zlm,en1,x1,y1,z1,
& thx1,thy1,scflag,nbx,nby,pl_del,rotang)
if((x0*x1.lt.0d0).and.(x1.gt.-1d5))
& write(6,*) 'Outcoil Crossover ',x0,x1 !check for crossovers
else
goto 100
end if
else
goto 100
end if
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
call hf1(5,real(x1),wgt)
call hf2(3333,real(x1),real(y1),wgt)
if (nbx.eq.0) nobx0=nobx0+1
if (nbx.eq.1) nobx1=nobx1+1
if (nbx.eq.2) nobx2=nobx2+1
if (nbx.eq.3) nobx3=nobx3+1
if (nby.eq.0) noby0=noby0+1
if (nby.eq.1) noby1=noby1+1
if (nby.eq.2) noby2=noby2+1
if (nby.eq.3) noby3=noby3+1
if (nbx.gt.3) write(6,*) ' 4+ Outcoil xBounces ',nbx
if (nby.gt.3) write(6,*) ' 4+ Outcoil yBounces ',nby
osc = osc + scflag
if(scflag.gt.5) write(6,*)' 5+ Outcoil Scatters ',scflag
if (scflag.gt.0) beamscflg = .true.
if ( ((en0.eq.en1).and.(scflag.gt.0)).or.
& ((en0.ne.en1).and.(scflag.eq.0)) )
& write(6,*) ' Edif/scatflg INPUT',scflag,en0,en1
c if (en0.ne.en1) write(6,*) ' edif out ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
c This path length is not OK for reflected neutrons -KG bec2011 no longer needed with new routines
c pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
c pl_del = oz/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2)))
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'out coil ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
end if
c
if ((bounceflg).and.(bounceflg2).and.((thx1.eq.-1.0*thx0).or.
& (thy1.eq.-1.0*thy0))) bounceflg3 = .true.
if ((bounceflg).and.((thx1.eq.-1.0*thx0).or.
& (thy1.eq.-1.0*thy0))) bounceflg2 = .true.
if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true.
capflout = capflout + sqrt(25.0)/rsqrt(en1)
*
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
cc if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through Al Wire output coil exit
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz6, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al outcoil',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al out2 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al output ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
************************************************************************
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
************************************************************************
* Neutrons passing through the gap (10.0 cm) between the output coil
* and the analyzing super mirror (ASM).
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,tairz5, en1, x1, y1, z1,
& thx1,thy1,airflg1)
c call collimator(x0,x1,y0,y1) !boroflex septum to prevent crossovers
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg Al 3',airflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR out/asm ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
NOutcoilout = NOutcoilout + wgt
if(x1.gt.ox1.and.x1.lt.ox2) NoutcoiloutE = NoutcoiloutE+wgt
if(x1.gt.ox3.and.x1.lt.ox4) NoutcoiloutW = NoutcoiloutW+wgt
OUTCOAVGflux = OUTCOAVGflux +
& wgt/((ox2-ox1)+(ox4-ox3))/(oy2-oy1)/nevts
if ((x1+0.5*(ox1-ox2))**2+y1**2.le.0.625**2)
& OUTCOPKflux = OUTCOPKflux + wgt/pi/0.625**2/nevts
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'gap 5 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
*************Output of OuputCoil Flux measurement******************
ccccccc scale factor is fluxscale*totcnt from readimage.f *~1.028 since readimage active area
ccccccc from readimage (>10cunts) = 24.08cm2 turns out to be slightly smaller than
ccccccc active PSM area here 4*(2.75X2.25)....smaller effect than 24.75/24.08=1.028 --> 1.01
cccccccc =1.2869E8*3.1356E5*1.01=4.076E12
call hf2(3333,real(x1),real(y1),4.076E12*wgt/nevts)
if (dabs(y1).le.0.5)
& call hf1(445,real(x1),4.076E12*wgt/0.9976/nevts*0.02)
if (dabs(x1+0.9).le.0.5)
& call hf1(446,real(y1),4.076E12*wgt/0.9976/nevts*0.02)
********************************************************************
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
* Neutrons passing through the super mirror."
c print *," Super Mirror"
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call sm(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
cot=cot+wgt
call hf1(6,real(x1),1.)
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'ASM ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' ASM thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
*
if ((bounceflg).and.(bounceflg2).and.((thx1.eq.-1.0*thx0).or.
& (thy1.eq.-1.0*thy0))) bounceflg3 = .true.
if ((bounceflg).and.((thx1.eq.-1.0*thx0).or.
& (thy1.eq.-1.0*thy0))) bounceflg2 = .true.
if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true.
c
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
*****************************************************
* Air gap (10cm) before Detector *
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call airgap(en0, x0, y0, z0, thx0, thy0,
& tairx1,tairx2,tairy1,tairy2,tairz6, en1, x1, y1, z1,
& thx1,thy1,airflg1)
c call collimator(x0,x1,y0,y1) !boroflex septum to prevent crossovers
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (airflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(airflg1)).or.
& ((en0.ne.en1).and.(.not.(airflg1))) )
& write(6,*) ' Edif/airflg Al 3',airflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif AIR 6 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'gap 6 ',pltot
if (pltot.gt.1300) goto 100
if (x1**2+y1**2.gt.detrad**2) goto 100
Ndetin = Ndetin + wgt
DETAVGflux = DETAVGflux + wgt/pi/detrad**2/nevts
if (x1**2+y1**2.le.0.625**2)
& DETPKflux = DETPKflux + wgt/pi/0.625**2/nevts
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
if (((dabs(x1).lt.1d0).and.(dabs(y1).lt.1d0)).and.
& ((dabs(thx1).lt.1d-1).and.(dabs(thy1).lt.1d-1)))
& zdet = z1
capfldet = capfldet + sqrt(25.0)/rsqrt(en1)
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
********************************************************************
********************************************************************
* Neutrons passing through SEVENTH Al window
*
c Update x,y,z,thx,thy,en values to new values
x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
call alwindow(en0, x0, y0, z0, thx0, thy0,
& talx1,talx2,taly1,taly2,talz7, en1, x1, y1, z1,
& thx1,thy1,alflg1)
if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
if (alflg1) beamscflg = .true.
if ( ((en0.eq.en1).and.(alflg1)).or.
& ((en0.ne.en1).and.(.not.(alflg1))) )
& write(6,*) ' Edif/alflg Al 7',alflg1,en0,en1
c if (en0.ne.en1) write(6,*) ' edif Al 7 ',airflg1,alflg1,
c & beamscflg,tar1flg,tar2flg
pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
pltot = pl_del + pltot
thtot = gy*mag_B_beam*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot
if (pltot.gt.1300) write(6,*) 'Al 7 ',pltot
if (dabs(thtot).gt.1.0e9) then
write(6,*) ' Al out thtot ', thtot, ', pl_del ',pl_del
write(6,*) ' en ',en1, ', pltot ',pltot
write(6,*) ' z ', z1
end if
ctemp
c if (x1.lt.0) write(6,*) ' Al7 EastDn z1 ',z1,', pathlen ',pltot
c if (x1.gt.0) write(6,*) ' Al7 WestUp z1 ',z1,', pathlen ',pltot
ctemp
*******************************************************************
ctemp
c if (x1.lt.0) write(6,2341) z1,pltot,thtot
c if (x1.gt.0) write(6,2342) z1,pltot,thtot
ctemp
*******************************************************************
c Apply the polarization product -KG
lambda=(2.*pi*h_par)/rsqrt(2*m_n*en1)
c if((lambda.lt.2.).or.(lambda.gt.14.)) then
c write(6,*)'lambda = ', lambda, ' out of range (2,14)Ang'
c endif
c PA=0.27873+3.5461*(lambda/10.)-6.2591*(lambda/10.)**2
c &+4.5036*(lambda/10.)**3-1.1912*(lambda/10.)**4
c Central PA. See Elog Sim. post 30
cc-----------------------------------------------
c Divide into 3 regions: west, middle and east
if (PAflg) then
if(x1.lt.(-0.95)) then !East
PA=-2.8274+47.729*(lambda/10.)-243.93*(lambda/10.)**2
& +614.37*(lambda/10.)**3-845.19*(lambda/10.)**4
& +649.46*(lambda/10.)**5-262.23*(lambda/10.)**6
& +43.397*(lambda/10.)**7
endif
if((x1.ge.(-0.95)).and.(x1.le.0.95)) then !Middle
PA=0.27873+3.5461*(lambda/10.)-6.2591*(lambda/10.)**2
& +4.5036*(lambda/10.)**3-1.1912*(lambda/10.)**4
endif
if(x1.gt.0.95) then !West
PA=0.27757+3.3155*(lambda/10.)-5.3255*(lambda/10.)**2
& +3.6243*(lambda/10.)**3-0.94529*(lambda/10.)**4
endif
else
PA=1
endif
cc-----------------------------------------------
c if((PA.lt.0.).or.(PA.gt.1.02)) then
c write(6,*)'PA = ',PA,' out of range (0,1)'
c endif c
c
c if(vect(6).gt.PA) goto 100
c
call hf2(9999,real(lambda),real(PA),wgt)
********************************************************************
*************** End of Beamline **********************************
*
c write(6,*) gy*mag_B_beam*480.6*0.1/rsqrt(2.*en1/m_n) ,
c & gy*mag_B_beam*pltot*0.1/rsqrt(2.*en1/m_n) ! mrad
c & , pltot,thtot
if((xstart*x1).lt.0.) Nxo=Nxo+wgt !count cross-over neutrons -KG
c if((xstart*x1).lt.0.) Nxo=Nxo+1 !-rcm
c count bounced neutrons
if (bounceflg) numbounce = numbounce + wgt
if (bounceflg2) numbounce2 = numbounce2 + wgt
if (bounceflg3) numbounce3 = numbounce3 + wgt
********************************************************************
call hf1(3001,real(x1),wgt)
call hf1(3002,real(y1),wgt)
call hf2(4444,real(x1),real(y1),wgt)
call hf1(3003,real(thx1),wgt)
call hf1(3004,real(thy1),wgt)
call hf1(3005,real(en1),wgt)
c lambda=(2.*pi*h_par)/rsqrt(2*m_n*en1) !done in PA application -KG
call hf1(4005,real(lambda),wgt)
c call hf1(3006,real(pltot),wgt)
c call hf1(3007,real(thtot),wgt)
************Rotations and Histograms for ALL and Non-Scattered Neutrons****************
if(x1.ge.0.) then ! x>0 = West (=old right)
call hf1(3010,real(pltot),wgt)
call hf1(3011,real(thtot),wgt)
call hf1(3013,real(en1),wgt)
call hf2(6020,real(enstart),real(thtot),wgt)
call hf2(6021,real(pltot),real(thtot),wgt)
if ((xtest.lt.0).and.(x1.gt.0)) then
call hf1(6035, real(pltottar),wgt)
call hf1(6036, real(thtot),wgt)
call hf2(6040,real(enstart),real(thtot),wgt)
call hf2(6041,real(pltot),real(thtot),wgt)
end if
if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting
rall = rall + PA*dsin(thtot*1d-3) !rcm PAsin
rallsq = rallsq + (PA*dsin(thtot*1d-3))**2
nrall = nrall + 1
nrallPA = nrallPA + PA
nrallPasq = nrallPasq + PA**2
rallavg = rall/nrall
rallPAavg = nrallPA/nrall
end if
rsrot = rsrot + wgt*PA*(dsin(thtot*1d-3)) !rcm PAsin
rsrotsq = rsrotsq + wgt*(PA*dsin(thtot*1d-3))**2
nr = nr +wgt
nrsq = nrsq + wgt**2
nrPA = nrPA + wgt*PA
nrPAsq = nrPAsq + wgt*PA**2
rPAavg = nrPA/nr
rsrotavg = rsrot/nr
!rcm
if (.not.((beamscflg.or.tar1flg).or.tar2flg)) then ! nonscattered
rsrotnsc = rsrotnsc + wgt*PA*dsin(thtot*1d-3)
rsrotnscsq = rsrotnscsq + wgt*(PA*dsin(thtot*1d-3))**2
nrnsc = nrnsc + wgt
nrnscsq = nrnscsq + wgt**2
nrnscPA = nrnscPA + wgt*PA
nrnscPAsq = nrnscPAsq + wgt*PA**2
rsnscPAavg = nrnscPA/nrnsc
rsrotnscavg = rsrotnsc/nrnsc
call hf1(6013, real(pltottar),wgt)
call hf1(6014, real(thtot),wgt)
call hf1(3017,real(en1),wgt)
c call hf2(6022,real(entar2),real(thtot),wgt)
c call hf2(6023,real(pltot),real(thtot),wgt)
c call hf2(6026,real(pltot),real(en1),wgt)
rplns = rplns + pltottar
rplnssq = rplnssq + pltottar**2
if (en1.ne.enstart) write(6,*) 'W nonscat delE ',en1-enstart,
& beamscflg,tar1flg,tar2flg
endif
else ! x<0 = East (=old left)
call hf1(3008,real(pltot),wgt)
call hf1(3009,real(thtot),wgt)
call hf1(3012,real(en1),wgt)
call hf2(6022,real(enstart),real(thtot),1.0) ! not trajectory dependent -- no wgt
call hf2(6023,real(pltot),real(thtot),wgt)
if ((xtest.gt.0).and.(x1.lt.0)) then
call hf1(6037, real(pltottar),wgt)
call hf1(6038, real(thtot),wgt)
call hf2(6042,real(enstart),real(thtot),1.0) ! not trajectory dependent -- no wgt
call hf2(6043,real(pltot),real(thtot),wgt)
end if
if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting
lall = lall + PA*dsin(thtot*1d-3)
lallsq = lallsq + (PA*dsin(thtot*1d-3))**2
nlall = nlall + 1
nlallPA = nlallPA + PA
nlallPAsq =nlallPAsq + PA**2
lallavg = lall/nlall
lallPAavg = nlallPA/nlall
end if
lsrot = lsrot + wgt*PA*dsin(thtot*1d-3) !rcm PAsin
lsrotsq = lsrotsq + wgt*(PA*dsin(thtot*1d-3))**2
nl = nl + wgt
nlsq = nlsq + wgt**2
nlPA = nlPA + wgt*PA
nlPAsq = nlPAsq + wgt*PA**2
lPAavg = nlPA/nl
lsrotavg = lsrot/nl
if (.not.((beamscflg.or.tar1flg).or.tar2flg)) then !rcm PAsin nonscattered
lsrotnsc = lsrotnsc + wgt*PA*dsin(thtot*1d-3)
lsrotnscsq = lsrotnscsq + wgt*(PA*dsin(thtot*1d-3))**2
nlnsc = nlnsc + wgt
nlnscsq = nlnscsq + wgt**2
nlnscPA = nlnscPA + wgt*PA
nlnscPAsq = nlnscPAsq + wgt*PA**2
lsnscPAavg = nlnscPA/nlnsc
lsrotnscavg = lsrotnsc/nlnsc
c call hf1(6009, real(pltot),wgt)
c call hf1(6010, real(thtot),wgt)
call hf1(6011, real(pltottar),wgt)
call hf1(6012, real(thtot),wgt)
call hf1(3018,real(en1),wgt)
c call hf2(6024,real(entar2),real(thtot),wgt)
c call hf2(6025,real(pltot),real(thtot),wgt)
c call hf2(6027,real(pltot),real(en1),wgt)
lplns = lplns + wgt*pltottar
lplnssq = lplnssq + wgt*pltottar**2
if (en1.ne.enstart) write(6,*) 'E nonscat delE ',en1-enstart,
& beamscflg,tar1flg,tar2flg
if (tar1flg) write(6,*) ' prob with tar1flg'
if (tar2flg) write(6,*) ' prob with tar2flg'
if (beamscflg) write(6,*) ' prob with beamscflg'
endif
endif
call hf1(7005, real(wgt),1.0)
if (wgt.lt.0) write(6,*) ' end loop wgt<0 ',wgt, sig, en1
if (enflg1) numendif1 = numendif1 + wgt
if (tar1flg) numtar1flg = numtar1flg + wgt
if (enflg2) numendif2 = numendif2 + wgt
if (tar2flg) numtar2flg = numtar2flg + wgt
if (en1.ne.enstart) numendifdet = numendifdet +wgt
if (beamscflg) numbeamsc = numbeamsc + wgt
if ( (en1.ne.enstart).and.((.not.(tar1flg)).and.(.not.(tar2flg))
&.and.(.not.(beamscflg))) ) write(6,*) ' Num Edif Prob ',en1,
& enstart,tar1flg,tar2flg,beamscflg
if ( (en1.eq.enstart).and.((tar1flg).or.(tar2flg)
&.or.(beamscflg)) ) write(6,*) ' Num Flg Prob ',en1,
& enstart,tar1flg,tar2flg,beamscflg
c if ( ((tar1flg).and.(tar2flg)).or.
c & ((tar1flg).and.(beamscflg)).or.
c & ((tar2flg).and.(beamscflg)) ) write(6,*) ' Dble Flg Prob ',en1,
c & enstart,tar1flg,tar2flg,beamscflg
****************Rotations and Histograms for Scattered neutrons****************
if ((beamscflg.or.tar1flg).or.tar2flg) then ! scattered
NEdet=NEdet+wgt
call hf1(7006, real(wgt),1.0)
c call hf1(6003, real(pltot),wgt)
c call hf1(6004, real(thtot),wgt)
if (x1.lt.0.) then ! x<0 = East (=old left)
call hf2(6031,real(enstart),real(thtot),wgt)
call hf1(6007, real(pltottar),wgt)
call hf1(6008, real(thtot),wgt)
if ((xtest.gt.0).and.(x1.lt.0)) then
call hf1(6017, real(pltottar),wgt)
call hf1(6018, real(thtot),wgt)
end if
call hf1(3015,real(en1),wgt)
c call hf2(6024,real(enstart),real(thtot),wgt)
c call hf2(6025,real(pltot),real(thtot),wgt)
c call hf2(6027,real(pltot),real(en1),wgt)
if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting
lsc = lsc + PA*dsin(thtot*1d-3)
lscsq = lscsq + (PA*dsin(thtot*1d-3))**2
nlsc = nlsc + 1
nlscPA = nlscPA + PA
nlscPAsq = nlscPAsq + PA**2
lscPAavg = nlscPA/nlsc
lscavg = lsc/nlsc
end if
if (tar2flg.and.(.not.(tar1flg))) then ! scat in East Dn only and on E side still
call hf1(7007, real(wgt),1.0)
lsrottsc = lsrottsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin
lsrottscsq = lsrottscsq + wgt*(PA*dsin(thtot*1d-3))**2
nltsraw = nltsraw + 1
nlts = nlts + wgt
nltssq = nltssq + wgt**2
nltsPA = nltsPA + wgt*PA
nltsPAsq = nltsPAsq + wgt*PA**2
ltsPAavg = nltsPA/nlts
lsrottscavg = lsrottsc/nlts
lwde = lwde + wgt*(enstart-en1)
lwdesq = lwdesq + wgt*(enstart-en1)**2
lpl = lpl + wgt*pltottar
lplsq = lplsq + wgt*pltottar**2
call hf1(3023,real(enstart-en1),wgt)
elseif (tar1flg) then ! scat in West Up i.e. a XTS (cross over tar scat)
lsrotxsc = lsrotxsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin
lsrotxscsq = lsrotxscsq + wgt*(PA*dsin(thtot*1d-3))**2
nlxsraw = nlxsraw + 1
nlxs = nlxs + wgt
nlxssq = nlxssq + wgt**2
nlxsPA = nlxsPA + wgt*PA
nlxsPAsq = nlxsPAsq + wgt*PA**2
lxsPAavg = nlxsPA/nlxs
lsrotxscavg = lsrotxsc/nlxs
ctemp if (wgt.eq.1.0) write(6,*) ' tar2 wgt=1 '
else ! no target scat, just Al or Gas scat
if (.not.(beamscflg)) write(6,*) ' prob with beamscflg ls'
lsrotgsc = lsrotgsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin
lsrotgscsq = lsrotgscsq + wgt*(PA*dsin(thtot*1d-3))**2
nlgsraw = nlgsraw + 1
nlgs = nlgs + wgt
nlgssq = nlgssq + wgt**2
nlgsPA = nlgsPA + wgt*PA
nlgsPAsq = nlgsPAsq + wgt*PA**2
lgsPAavg = nlgsPA/nlgs
lsrotgscavg = lsrotgsc/nlgs
ctemp if (wgt.ne.1.0) write(6,*) ' tar2 wgt',wgt,tar1flg,tar2flg,xtar1,
c & xtar2,x1
end if
c call hf1(3020,real(enstart-en1),wgt)
call hf1(3021,real(enstart-en1),wgt)
c call hf2(8011,real(en1),real(thx1),wgt)
c call hf2(8012,real(enstart-en1),real(thx1),wgt)
lsde = lsde + wgt*(enstart-en1)
lsdesq = lsdesq + wgt*(enstart-en1)**2
ndel = ndel + wgt
ndelsq = ndelsq + wgt**2
c write(6,*) ' 1W eloss tar1 tar2 ', entar1-entar2
else ! x>0 = West (=old right)
call hf2(6030,real(enstart),real(thtot),wgt)
call hf1(6005, real(pltottar),wgt)
call hf1(6006, real(thtot),wgt)
if ((xtest.lt.0).and.(x1.gt.0)) then
call hf1(6015, real(pltottar),wgt)
call hf1(6016, real(thtot),wgt)
end if
call hf1(3016,real(en1),wgt)
c call hf2(6022,real(enstart),real(thtot),wgt)
c call hf2(6023,real(pltot),real(thtot),wgt)
c call hf2(6026,real(pltot),real(en1),wgt)
if (.not.((tar1flg.or.tar2flg).and.(vect(5).gt.wgt2))) then ! alternate weighting
rsc = rsc + PA*dsin(thtot*1d-3)
rscsq = rscsq + (PA*dsin(thtot*1d-3))**2
nrsc = nrsc + 1
nrscPA = nrscPA + PA
nrscPAsq = nrscPAsq + PA**2
rscPAavg = nrscPA/nrsc
rscavg = rsc/nrsc
end if
if (tar1flg.and.(.not.(tar2flg))) then ! scat West Up only and still in West
call hf1(7007, real(wgt),1.0) !rcm PAsin
rsrottsc = rsrottsc + wgt*PA*dsin(thtot*1d-3)
rsrottscsq = rsrottscsq + wgt*(PA*dsin(thtot*1d-3))**2
nrtsraw = nrtsraw + 1
nrts = nrts + wgt
nrtssq = nrtssq + wgt**2
nrtsPA = nrtsPA + wgt*PA
nrtsPAsq = nrtsPAsq + wgt*PA**2
rtsPAavg = nrtsPA/nrts
rsrottscavg = rsrottsc/nrts
rwde = rwde + wgt*(enstart-en1)
rwdesq = rwdesq + wgt*(enstart-en1)**2
rpl = rpl + wgt*pltottar
rplsq = rplsq + wgt*pltottar**2
call hf1(3024,real(enstart-en1),wgt)
else if (tar2flg) then ! scat East Dn so an XTS whether also scat in up or not
rsrotxsc = rsrotxsc + wgt*PA*dsin(thtot*1d-3)
rsrotxscsq = rsrotxscsq + wgt*(PA*dsin(thtot*1d-3))**2
nrxsraw = nrxsraw + 1
nrxs = nrxs + wgt
nrxssq = nrxssq + wgt**2
nrxsPA = nrxsPA + wgt*PA
nrxsPAsq = nrxsPAsq + wgt*PA**2
rxsPAavg = nrxsPA/nrxs
rsrotxscavg = rsrotxsc/nrxs
ctemp if (wgt.eq.1.0) write(6,*) ' tar1 wgt=1 '
else
if (.not.(beamscflg)) write(6,*) ' prob with beamscflg rs'
rsrotgsc = rsrotgsc + wgt*PA*dsin(thtot*1d-3) !rcm PAsin
rsrotgscsq = rsrotgscsq + wgt*(PA*dsin(thtot*1d-3))**2
nrgsraw = nrgsraw + 1
nrgs = nrgs + wgt
nrgssq = nrgssq + wgt**2
nrgsPA = nrgsPA + wgt*PA
nrgsPAsq = nrgsPAsq + wgt*PA**2
rgsPAavg = nrgsPA/nrgs
rsrotgscavg = rsrotgsc/nrgs
ctemp if (wgt.ne.1.0) write(6,*) ' tar1 wgt',wgt,tar1flg,tar2flg,xtar1,
c & xtar2,x1
end if
c call hf1(3020,real(enstart-en1),wgt)
call hf1(3022,real(enstart-en1),wgt)
c call hf2(8011,real(en1),real(thx1),wgt)
c call hf2(8012,real(enstart-en1),real(thx1),wgt)
rsde = rsde + wgt*(enstart-en1)
rsdesq = rsdesq + wgt*(enstart-en1)**2
nder = nder + wgt
ndersq = ndersq + wgt**2
c write(6,*) ' 2E eloss tar1 tar2 ', entar1-entar2
endif
endif
********************************************************************
*
alambda = 2*pi*h_par/rsqrt(2*m_n*en1)
c ispec = int(alambda/0.016)
ispec = int(alambda/0.030)
if (ispec.gt.1000) then
write(6,*) ' ispec:',ispec,' set to 1000, alambda:',
& alambda
ispec = 1000
end if
spec(ispec,1) = alambda
spec(ispec,2) = spec(ispec,2) + 1
********************************************************************
if ((xtar1*xtar2.gt.0.0).and. ! xtar1 = x-pos before UPtar
& ((x1*xtar1.lt.0.0).or.(x1*xtar2.lt.0.0)) ) then
c write(6,*) xtar1, xtar2, x1
dncross = dncross + wgt
end if
c if (abs(thtot).gt.300.0) then
c write(6,941) xtest,x1,enstart,en1,pltot,thtot
c endif
941 format(1x,1pe9.2,2x,1pe9.2,2x,1pe9.2,2x,1pe9.2,2x,
& 1pe12.5,2x,1pe9.2)
if ((xtest.lt.0).and.(x1.gt.0)) numE2W = numE2W + wgt
if ((xtest.gt.0).and.(x1.lt.0)) numW2E = numW2E + wgt
34 format(i12,2x,a6)
ctemp
c if (x1.lt.0) write(6,*) 'End EastDn ',z1,pltot ,thtot
c if (x1.gt.0) write(6,*) 'End WestUp ',z1,pltot,thtot
ctemp
if (pltottar.gt.(tarlength+2.0)) then
write(6,*) ' LongPlL ',z1,pltottar,thx1,thy1
end if
100 do j=1,10
lim = float(j)*float(nevts)/10.
if(float(i).eq.lim) then
write(*,34) i, 'events'
c write(*,*) 'rsrot=',rsrot,' rsrotsq=',rsrotsq,' nrPA=',nrPA !-KG !rcm
c write(*,*) 'lsrot=',lsrot,' lsrotsq=',lsrotsq,' nlPa=',nlPA !-KG !rcm
endif
enddo
101 continue
do j=1, 1000
cbec write(71,*) j, spec(j,1), spec(j,2)
enddo
c write(*,*) 'calling nhe3'
call nhe3(spec)
print *
print *, " through loop, writing to ", filename
write(6,*) ' num of tar1flg ',numtar1flg
write(6,*) ' num of endif 1',numendif1
write(3,*) ' num of tar1flg ',numtar1flg
write(3,*) ' num of endif 1',numendif1
write(6,*) ' num of tar2flg ',numtar2flg
write(6,*) ' num of endif 2',numendif2
write(3,*) ' num of tar2flg ',numtar2flg
write(3,*) ' num of endif 2',numendif2
write(6,*) ' num of beamscflg ',numbeamsc
write(3,*) ' num of beamscflg ',numbeamsc
write(6,*) ' num Tar scat West ',nrts
write(6,*) ' num Tar scat East ',nlts
write(3,*) ' num Tar scat West ',nrts
write(3,*) ' num Tar scat East ',nlts
write(6,*) ' num Gas scat West ',nrgs
write(6,*) ' num Gas scat East ',nlgs
write(3,*) ' num Gas scat West ',nrgs
write(3,*) ' num Gas scat East ',nlgs
write(6,*) ' num of endet.ne.enstart ',numendifdet
write(3,*) ' num of endet.ne.enstart ',numendifdet
write(6,*) ' num E to W at det ',numE2W
write(6,*) ' num W to E at det ',numW2E
write(3,*) ' num E to W at det ',numE2W
write(3,*) ' num W to E at det ',numW2E
print *
if (hegastar) then
write(6,*) ' 4He gas target ', ndenhegas
write(3,*) ' 4He gas target ', ndenhegas
elseif (airtar) then
write(6,*) ' Air target ', ndenair
write(3,*) ' Air target ', ndenair
else
write(6,*) ' Liquid He target; Temp= ',T,' K '
write(3,*) ' Liquid He target; Temp= ',T,' K '
endif
print *
write(*,*) 'input ',nevts
write(*,*) 'PSM output ',NPSMout,NPSMout/NPSMout
write(*,*) 'Air gap 1 out ',Nairout1,Nairout1/NPSMout
write(*,*) 'Al window 1 out ',Nalout1,Nalout1/NPSMout
write(*,*) 'Inguide Output ',Nguideout,Nguideout/NPSMout ! rcm
write(*,*) 'Al window 2 out ',Nalout2,Nalout2/NPSMout ! rcm
write(*,*) 'Input Coil output ',Nincoilout,Nincoilout/NPSMout
write(*,*) 'WE target input ',NtarWEin,NtarWEin/NPSMout
write(*,*) 'Up West target IN ',Ntar1,Ntar1/NPSMout
write(*,*) 'Up West target OUT ',Ntar1out,Ntar1out/NPSMout
write(*,*) 'Up East empty IN ',Nemp1in,Nemp1in/NPSMout
write(*,*) 'Up East empty Out ',Nemp1out,Nemp1out/NPSMout
write(*,*) 'Dn West empty IN ',Nemp2in,Nemp2in/NPSMout
write(*,*) 'DN West empty Out ',Nemp2out,Nemp2out/NPSMout
write(*,*) 'Dn East target IN ',Ntar2,Ntar2/NPSMout
write(*,*) 'Dn East target OUT ',Ntar2out,Ntar2out/NPSMout
write(*,*) 'E-W Target Outputs ',NTargetout,NTargetout/NPSMout
write(*,*) 'Output Coil output ',NOutcoilout,NOutcoilout/NPSMout
write(*,*) 'ASM output (geom.) ', dble(cot), cot/NPSMout
write(*,*) 'ASM out scattered ',NEdet,NEdet/NPSMout
write(*,*) 'Detector input ',Ndetin,Ndetin/NPSMout
write(*,*)
write(*,5551) 0.25/xmax/ymax,INPKflux,0.25/xmax/ymax/INPKflux
write(*,5552) PSMAVGflux,PSMPKflux,PSMAVGflux/PSMPKflux
write(*,5553) INCOAVGflux,INCOPKflux,INCOAVGflux/INCOPKflux
write(*,5554) TARAVGflux,TARPKflux,TARAVGFlux/TARPKFlux
write(*,5555) OUTCOAVGflux,OUTCOPKflux,OUTCOAVGflux/OUTCOPKflux
write(*,5556) DETAVGflux,DETPKflux,DETAVGFlux/DETPKFlux
write(*,*)
write(*,5560) INCOAVGFlux/PSMAVGflux,INCOPKFlux/PSMPKflux,
& INCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflinc/nincoilout
write(*,5561) OUTCOAVGFlux/PSMAVGflux,OUTCOPKFlux/PSMPKflux,
& OUTCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout
write(*,5562) DETAVGFlux/PSMAVGflux,DETPKFlux/PSMPKflux,
& DETPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout
write(*,*)
write(3,5560) INCOAVGFlux/PSMAVGflux,INCOPKFlux/PSMPKflux,
& INCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflinc/nincoilout
write(3,5561) OUTCOAVGFlux/PSMAVGflux,OUTCOPKFlux/PSMPKflux,
& OUTCOPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout
write(3,5562) DETAVGFlux/PSMAVGflux,DETPKFlux/PSMPKflux,
& DETPKFlux/PSMPKflux/capflpsm*nPSMout*capflout/NOutcoilout
write(3,*)
5560 format(1x,' Incoil/PSM av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', CAPpk',1x,1pe10.3)
5561 format(1x,' Outcoil/PSM av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', CAPpk',1x,1pe10.3)
5562 format(1x,' Detector/PSM av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', CAPpk',1x,1pe10.3)
write(*,*) ' Cap Flux ratio PSM ', capflpsm/nPSMout
write(*,*) ' Cap Flux ratio Incoil ', capflinc/nincoilout
write(*,*) ' Cap Flux ratio Outcoil ', capflout/NOutcoilout
write(*,*) ' Cap Flux ratio Detector ', capfldet/ndetin
write(*,*)
write(*,5651) xmax,ymax
write(*,5652) xin,yin
if (sptmflg) write(*,5659) ix1,ix2,ix3,ix4,iy1,iy2
if (.not.(sptmflg)) write(*,5653) ix1,ix4,iy1,iy2
write(*,5654) tarx1,tarx2,tarx3,tarx4
write(*,5655) tary1,tary2
write(*,5656) ox1,ox2,ox3,ox4
write(*,5657) oy1,oy2
write(*,5658) detrad
write(*,*)
5551 format(1x,' Input Flux av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', av/pk',1x,1pe10.3)
5552 format(1x,' PSM Flux av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', av/pk',1x,1pe10.3)
5553 format(1x,' InCoil Flux av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', av/pk',1x,1pe10.3)
5554 format(1x,' Target Flux av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', av/pk',1x,1pe10.3)
5555 format(1x,' Outcoil Flux av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', av/pk',1x,1pe10.3)
5556 format(1x,' Detector Flux av',1x,1pe10.3,', pk',1x,1pe10.3,
& ', av/pk',1x,1pe10.3)
5651 format(1x,' Input Dim (cm)',1x,1f5.2,1x,1f5.2)
5652 format(1x,' PSM Dim (cm)',1x,1f5.2,1x,1f5.2)
5653 format(1x,' Incoil Dim (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2,
& 1x,1f5.2)
5659 format(1x,' Incoil Dim (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2,
& 1x,1f5.2,1x,1f5.2,1x,1f5.2)
5654 format(1x,' Target Dim X (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2,
& 1x,1f5.2)
5655 format(1x,' Target Dim Y (cm)',1x,1f5.2,1x,1f5.2)
5656 format(1x,' Outcol Dim X (cm)',1x,1f5.2,1x,1f5.2,1x,1f5.2,
& 1x,1f5.2)
5657 format(1x,' Outcoil Dim Y (cm)',1x,1f5.2,1x,1f5.2)
5658 format(1x,' Detector radius (cm)',1x,1f5.2)
write(3,*) 'input ',nevts
write(3,*) 'PSM output ',NPSMout,NPSMout/NPSMout
write(3,*) 'Air gap 1 out ',Nairout1,Nairout1/NPSMout
write(3,*) 'Al window 1 out ',Nalout1,Nalout1/NPSMout
write(3,*) 'Inguide Output ',Nguideout,Nguideout/NPSMout ! rcm
write(3,*) 'Al window 2 out ',Nalout2,Nalout2/NPSMout ! rcm
write(3,*) 'Input Coil output ',Nincoilout,Nincoilout/NPSMout
write(3,*) 'WE target input ',NtarWEin,NtarWEin/NPSMout
write(3,*) 'Up West target IN ',Ntar1,Ntar1/NPSMout
write(3,*) 'Up West target OUT ',Ntar1out,Ntar1out/NPSMout
write(3,*) 'Up East empty IN ',Nemp1in,Nemp1in/NPSMout
write(3,*) 'Up East empty Out ',Nemp1out,Nemp1out/NPSMout
write(3,*) 'Dn West empty IN ',Nemp2in,Nemp2in/NPSMout
write(3,*) 'DN West empty Out ',Nemp2out,Nemp2out/NPSMout
write(3,*) 'Dn East target IN ',Ntar2,Ntar2/NPSMout
write(3,*) 'Dn East target OUT ',Ntar2out,Ntar2out/NPSMout
write(3,*) 'E-W Target Outputs ',NTargetout,NTargetout/NPSMout
write(3,*) 'Output Coil output ',NOutcoilout,NOutcoilout/NPSMout
write(3,*) 'ASM output (geom.) ', dble(cot), cot/NPSMout
write(3,*) 'ASM out scattered ',NEdet,NEdet/NPSMout
write(3,*) 'Detector input ',Ndetin,Ndetin/NPSMout
write(3,*)
write(3,5551) 0.25/xmax/ymax,INPKflux,0.25/xmax/ymax/INPKflux
write(3,5552) PSMAVGflux,PSMPKflux,PSMAVGflux/PSMPKflux
write(3,5553) INCOAVGflux,INCOPKflux,INCOAVGflux/INCOPKflux
write(3,5554) TARAVGflux,TARPKflux,TARAVGFlux/TARPKFlux
write(3,5555) OUTCOAVGflux,OUTCOPKflux,OUTCOAVGflux/OUTCOPKflux
write(3,5556) DETAVGflux,DETPKflux,DETAVGFlux/DETPKFlux
write(3,*)
write(3,*) ' Incoil/PSM av',INCOAVGFlux/PSMAVGflux,
& ', pk',INCOPKFlux/PSMPKflux
write(3,*) ' Outcoil/PSM av',OUTCOAVGFlux/PSMAVGflux,
& ', pk',OUTCOPKFlux/PSMPKflux
write(3,*) ' Detector/PSM av',DETAVGFlux/PSMAVGflux,
& ', pk',DETPKFlux/PSMPKflux
write(3,*)
write(3,*) ' Cap Flux ratio PSM ', capflpsm/nPSMout
write(3,*) ' Cap Flux ratio Incoil ', capflinc/nincoilout
write(3,*) ' Cap Flux ratio Outcoil ', capflout/NOutcoilout
write(3,*) ' Cap Flux ratio Detector ', capfldet/ndetin
write(3,*)
write(3,5651) xmax,ymax
write(3,5652) xin,yin
if (sptmflg) write(3,5659) ix1,ix2,ix3,ix4,iy1,iy2
if (.not.(sptmflg)) write(3,5653) ix1,ix4,iy1,iy2
write(3,5654) tarx1,tarx2,tarx3,tarx4
write(3,5655) tary1,tary2
write(3,5656) ox1,ox2,ox3,ox4
write(3,5657) oy1,oy2
write(3,5658) detrad
write(3,*)
write(*,*) 'Number scattered in Al window 1 ',Nalscat1
write(3,*) 'Number scattered in Al window 1 ',Nalscat1
write(*,*) 'Number scattered in Al window 2 ',Nalscat2
write(3,*) 'Number scattered in Al window 2 ',Nalscat2 !rcm
write(*,*) 'Number scattered in Air gap 1 ',Nairscat1
write(3,*) 'Number scattered in Air gap 1 ',Nairscat1
write(*,*) 'Number scattered in UpTar ',Nenscat1
write(3,*) 'Number scattered in UpTar ',Nenscat1
write(*,*) 'Number scattered in DnTar ',Nenscat2
write(3,*) 'Number scattered in DnTar ',Nenscat2
write(*,*) 'Number scattered in EmptyTar Up ',Nempty1
write(*,*) 'Number scattered in EmptyTar Dn ',Nempty2
write(3,*) 'Number scattered in EmptyTar Up ',Nempty1
write(3,*) 'Number scattered in EmptyTar Dn ',Nempty2
c write(*,*) 'Number Theta scattered in Target 1 ',Nthscat1
c write(3,*) 'Number Theta scattered in Target 1 ',Nthscat1
c write(*,*) 'Number Theta scattered in Target 2 ',Nthscat2
c write(3,*) 'Number Theta scattered in Target 2 ',Nthscat2
c
c -------------------------
if(ranno<1e1) then
Write(outfname,'(1Hh,I1,6H.hbook)')ranno
elseif(ranno<1e2) then
Write(outfname,'(1Hh,I2,6H.hbook)')ranno
elseif(ranno<1e3) then
Write(outfname,'(1Hh,I3,6H.hbook)')ranno
elseif(ranno<1e4) then
Write(outfname,'(1Hh,I4,6H.hbook)')ranno
elseif(ranno<1e5) then
Write(outfname,'(1Hh,I5,6H.hbook)')ranno
elseif(ranno<1e6) then
Write(outfname,'(1Hh,I6,6H.hbook)')ranno
elseif(ranno<1e7) then
Write(outfname,'(1Hh,I7,6H.hbook)')ranno
elseif(ranno<1e8) then
Write(outfname,'(1Hh,I8,6H.hbook)')ranno
endif
call hrput(0,outfname,'N')
c -------------------------
c Adding some printouts at end of run -- bec
print *
Call Date_and_time(date,time,zone,datetime)
write(6,*) ' End time: ', date, time
write(6,*) ' HBOOK filename: ',outfname
write(6,*)
c*********averages -- Note: Now with weights, nr, nl, etc are all sums of
c******* weights, see logbook (bec)
fsw = nrts/nr
ufsw = fsw*rsqrt(1.0/dble(nrtsraw)+1.0/dble(nr))
fse = nlts/nl
ufse = fse*rsqrt(1.0/dble(nltsraw)+1.0/dble(nl))
fxw = nrxs/nr
ufxw = fxw*rsqrt(1.0/dble(nrxsraw)+1.0/dble(nr))
fxe = nlxs/nl
ufxe = fxe*rsqrt(1.0/dble(nlxsraw)+1.0/dble(nl))
write(6,755) Nenscat1/Ntar1
write(6,756) Nenscat2/Ntar2
write(6,757) Ntar12/Ntar1
write(6,758) cot/float(Nevts)
write(6,759) fsw,ufsw
write(6,760) fse,ufse
write(6,761) nrgs/nr,nrgs/nr*dsqrt(1/nrgsraw+1/nr)
write(6,762) nlgs/nl,nlgs/nl*dsqrt(1/nlgsraw+1/nl)
write(6,765) fxw,ufxw
write(6,766) fxe,ufxe
write(6,767) nrnsc/nr,nrnsc/nr*dsqrt(1/nrnsc+1/nr)
write(6,768) nlnsc/nl,nlnsc/nl*dsqrt(1/nlnsc+1/nl)
write(6,*)
write(6,*) ' StDev All West ',
& rsqrt((rsrotsq-rsrot**2/nr)/nr)*1.0d3
write(6,*) ' St Err All West ',
& rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr*1.0d3 !rcm
c write(6,*) ' nrPA ',nrPA , ' nrPAsq ',nrPAsq
c write(6,*) 'rsqrt(nrPAsq)/nrPA ', rSqrt(nrPAsq)/nrPA
c
write(6,*) ' StDev NoScat West ',
& rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*1.0d3
write(6,*) ' St Err NoScat West ',
& rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*
& rSqrt(nrnscsq)/nrnsc*1.0d3
c write(3,*)
write(3,*) ' StDev XOvTar Scat West ',
& rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1.0d3
write(3,*) ' St Err XOvTar Scat West ',
& rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1.0d3
& *rSqrt(nrxssq)/nrxs
c write(3,*) ' nrxs ',nrxs , ' nrxssq ',nrxssq
c write(3,*) 'rsqrt(nrxssq)/nrxs ', rSqrt(nrxssq)/nrxs
if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then
c write(6,*)
write(6,*) ' StDev Tar Scat West ',
& rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1.0d3
write(6,*) ' St Err Tar West ',
& rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1.0d3
& *rSqrt(nrtssq)/nrts
c write(6,*) ' nrts ',nrts , ' nrtssq ',nrtssq
c write(6,*) 'rsqrt(nrtssq)/nrts ', rSqrt(nrtssq)/nrts
c write(6,*)
end if
if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then
write(6,*) ' StDev Gas Scat West ',
& rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1.0d3
write(6,*) ' St Err Gas West ',
& rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1.0d3
& *rSqrt(nrgssq)/nrgs
c write(6,*) ' nrgs ',nrgs , ' nrgssq ',nrgssq
c write(6,*) 'rsqrt(nrgssq)/nrgs ', rSqrt(nrgssq)/nrgs
end if
write(6,*)
write(6,*)
write(6,713) nr, nl
write(6,714) nr/nl
write(6,717) nr-nl,rsqrt(nlsq+nrsq)
write(6,*)
write(6,735) NguideoutW,NguideoutE
write(6,736) NincoiloutW,NincoiloutE
write(6,737) NoutcoiloutW,NoutcoiloutE
write(6,*)
write(6,770) nrnsc , nlnsc
write(6,*)
write(6,700) zdet, mag_B,mag_B_beam, bzgradnum
write(6,*)
write(6,*) ' ********************** Average Rotation Values
&***********************'
write(6,*)
argrsrot = (rsrotsq-rsrot**2/nr)/nr
argrsrotPA = (nrPAsq-nrPA**2/nr)/nr
c if (argrsrot.lt.0) then
c write(6,*) 'argrsrot',argrsrot
c argrsrot =0
c endif
c if (argrsrotPA.lt.0) then
c write(6,*) 'argrsrotPA',argrsrotPA
c argrsrotPA =0
c endif
uncrsrotavg = rsqrt(argrsrot)*rsqrt(nrsq)/nr
uncrPAavg = rsqrt(argrsrotPA)*rsqrt(nrsq)/nr
arglsrot = (lsrotsq-lsrot**2/nl)/nl
arglsrotPA = (nlPAsq-nlPA**2/nl)/nl
c if (arglsrot.lt.0) then
c write(6,*) 'arglsrot',arglsrot
c arglsrot =0
c endif
c if (arglsrotPA.lt.0) then
c write(6,*) 'arglsrotPA',arglsrotPA
c arglsrotPA =0
c endif
unclsrotavg = rsqrt(arglsrot)*rsqrt(nlsq)/nl
unclPAavg = rsqrt(arglsrotPA)*rsqrt(nlsq)/nl
write(6,704) real(rPAavg*100),dasin(rsrotavg/rPAavg)*1d3, !rcm
& rsqrt((uncrsrotavg/rPAavg)**2+
& (uncrPAavg*rsrotavg/rPAavg**2)**2)/
& rsqrt(1-(rsrotavg/rPAavg)**2)*1d3
write(6,703) real(lPAavg*100),dasin(lsrotavg/lPAavg)*1d3,
& rsqrt((unclsrotavg/lPAavg)**2+
& (unclPAavg*lsrotavg/lPAavg**2)**2)/
& rsqrt(1-(lsrotavg/lPAavg)**2)*1d3
write(6,711) dasin(rsrotavg/rPAavg)*1d3-
& dasin(lsrotavg/lPAavg)*1d3,
& rsqrt( ((uncrsrotavg/rPAavg)**2+
& (uncrPAavg*rsrotavg/rPAavg**2)**2)/
& (1-(rsrotavg/rPAavg)**2) + ((unclsrotavg/lPAavg)**2+
& (unclPAavg*lsrotavg/lPAavg**2)**2)/
& (1-(lsrotavg/lPAavg)**2))*1d3
c write(6,704) rsrot/nr,
c & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr
c write(6,703) lsrot/nl,
c & rsqrt((lsrotsq-lsrot**2/nl)/nl)*rSqrt(nlsq)/nl
c write(6,711) rsrot/nr-lsrot/nl,
c & rsqrt( (lsrotsq-lsrot**2/nl)/nl*nlsq/nl**2
c & + (rsrotsq-rsrot**2/nr)/nr*nrsq/nr**2 )
c if ((nrall.gt.1).and.(nlall.gt.1)) then
c write(6,*)
c write(6,804) rall/float(nrall), ! alternate weighting method
c & rsqrt((rallsq-rall**2/nrall)/nrall/(nrall-1))
c write(6,803) lall/nlall,
c & rsqrt((lallsq-lall**2/nlall)/nlall/(nlall-1))
c write(6,811) rall/nrall-lall/nlall,
c & rsqrt( (lallsq-lall**2/nlall)/nlall/(nlall-1)
c & + (rallsq-rall**2/nrall)/nrall/(nlall-1) )
c else
c write(6,*) 'No alternate weighting data '
c end if
uncrallavg = rsqrt((rallsq-rall**2/nrall)/ !rcm PAsin
& (nrall-1)/nrall)
uncrallPAavg = rsqrt((nrallPAsq - nrallPA**2/nrall)/
& (nrall-1)/nrall)
unclallavg = rsqrt((lallsq-lall**2/nlall)/
& (nlall-1)/nlall)
unclallPAavg = rsqrt((nlallPAsq - nlallPA**2/nlall)/
& (nlall-1)/nlall)
write(6,*)
if ((nrall.gt.1).and.(nlall.gt.1)) then
write(6,804) real(rallPAavg*100),dasin(rallavg/rallPAavg)*1d3, ! alternate weighting method
& rsqrt( (uncrallavg/rallPAavg)**2+
& (uncrallPAavg*rallavg/rallPAavg**2)**2 )/
& rsqrt(1-(rallavg/rallPAavg)**2)*1d3
write(6,803) real(lallPAavg*100),dasin(lallavg/lallPAavg)*1d3,
& rsqrt( (unclallavg/lallPAavg)**2+
& (unclallPAavg*lallavg/lallPAavg**2)**2 )/
& rsqrt(1-(lallavg/lallPAavg)**2)*1d3
write(6,811) dasin(rallavg/rallPAavg)*1d3-
& dasin(lallavg/lallPAavg)*1d3,
& rsqrt( ((uncrallavg/rallPAavg)**2+
& (uncrallPAavg*rallavg/rallPAavg**2)**2)/
& (1-(rallavg/rallPAavg)**2) + ((unclallavg/lallPAavg)**2+
& (unclallPAavg*lallavg/lallPAavg**2)**2)/
& (1-(lallavg/lallPAavg)**2))*1d3
else
write(6,*) 'No alternate weighting data '
end if
write(6,*)
argrsnsc = (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc
argrsnscPA = (nrnscPAsq-nrnscPA**2/nrnsc)/nrnsc
c if (argrsnsc.lt.0) then
c write(6,*) 'argrsnsc',argrsnsc
c argrsnsc =0
c endif
c if (argrsnscPA.lt.0) then
c write(6,*) 'argrsnscPA',argrsnscPA
c argrsnscPA =0
c endif
uncrsnsc = rsqrt(argrsnsc)*rsqrt(nrnscsq)/nrnsc
uncrsnscPA = rsqrt(argrsnscPA)*rsqrt(nrnscsq)/nrnsc
arglsnsc = (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc
arglsnscPA = (nlnscPAsq-nlnscPA**2/nlnsc)/nlnsc
c if (arglsnsc.lt.0) then
c write(6,*) 'arglsnsc',arglsnsc
c arglsnsc =0
c endif
c if (arglsnscPA.lt.0) then
c write(6,*) 'arglsnscPA',arglsnscPA
c arglsnscPA =0
c endif
unclsnsc = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc
unclsnscPA = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc
write(6,726) real(rsnscPAavg*100),
& dasin(rsrotnscavg/rsnscPAavg)*1d3,
& rsqrt((uncrsnsc/rsnscPAavg)**2+
& (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/
& rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3
write(6,725) real(lsnscPAavg*100),
& dasin(lsrotnscavg/lsnscPAavg)*1d3,
& rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm
& (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/
& rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3
write(6,732)dasin(rsrotnscavg/rsnscPAavg)*1d3-
& dasin(lsrotnscavg/lsnscPAavg)*1d3,
& rsqrt( ((uncrsnsc/rsnscPAavg)**2+
& (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/
& (1-(rsrotnscavg/rsnscPAavg)**2) + ((unclsnsc/lsnscPAavg)**2+
& (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/
& (1-(lsrotnscavg/lsnscPAavg)**2) )*1d3
c write(6,726) rsrotnsc/nrnsc, !rcm
c & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)
c & *rSqrt(nrnscsq)/nrnsc
c write(6,725) lsrotnsc/nlnsc,
c & rsqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc)
c & *rSqrt(nlnscsq)/nlnsc
c write(6,732) rsrotnsc/nrnsc - lsrotnsc/nlnsc,
c & rsqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc
c & *nlnscsq/nlnsc**2
c & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc
c & *nrnscsq/nrnsc**2)
write(6,*)
if ((nrsc.gt.1).and.(nlsc.gt.1)) then !rcm
uncrscavg = rsqrt((rscsq-rsc**2/nrsc)/
& (nrsc-1)/nrsc)
uncrscPAavg = rsqrt((nrscPAsq - nrscPA**2/nrsc)/
& (nrsc-1)/nrsc)
unclscavg = rsqrt((lscsq-lsc**2/nlsc)/
& (nlsc-1)/nlsc)
unclscPAavg = rsqrt((nlscPAsq - nlscPA**2/nlsc)/
& (nlsc-1)/nlsc)
write(6,806) real(rscPAavg*100),dasin(rscavg/rscPAavg)*1d3, ! alternate weighting method
& rsqrt( (uncrscavg/rscPAavg)**2+
& (uncrscPAavg*rscavg/rscPAavg**2)**2 )/
& rsqrt(1-(rscavg/rscPAavg)**2)*1d3
write(6,805) real(lscPAavg*100),dasin(lscavg/lscPAavg)*1d3,
& rsqrt( (unclscavg/lscPAavg)**2+
& (unclscPAavg*lscavg/lscPAavg**2)**2 )/
& rsqrt(1-(lscavg/lscPAavg)**2)*1d3
write(6,812) dasin(rscavg/rscPAavg)*1d3-dasin(lscavg/lscPAavg)*1d3,
& rsqrt( ((uncrscavg/rscPAavg)**2+(uncrscPAavg*rscavg/
& rscPAavg**2)**2)/(1-(rscavg/rscPAavg)**2) +
& ((unclscavg/lscPAavg)**2+(unclscPAavg*lscavg/
& lscPAavg**2)**2)/(1-(lscavg/lscPAavg)**2) )*1d3
else
write(6,*) ' No alternate weighting scatters'
end if
write(6,*)
if ((nlts.gt.0).and.(nrts.gt.0)) then
argrtsc = (rsrottscsq-rsrottsc**2/nrts)/nrts
argrtscPA = (nrtsPAsq-nrtsPA**2/nrts)/nrts
uncrstsc = rsqrt(argrtsc)*rsqrt(nrtssq)/nrts
uncrstscPA = rsqrt(argrtscPA)*rsqrt(nrtssq)/nrts
argltsc = (lsrottscsq-lsrottsc**2/nlts)/nlts
argltscPA = (nltsPAsq-nltsPA**2/nlts)/nlts
unclstsc = rsqrt(argltsc)*rsqrt(nltssq)/nlts
unclstscPA = rsqrt(argltscPA)*rsqrt(nltssq)/nlts
write(6,706) real(rtsPAavg*100),dasin(rsrottscavg/rtsPAavg)*1d3,
& rsqrt((uncrstsc/rtsPAavg)**2+
& (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/
& rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3
write(6,705) real(ltsPAavg*100),dasin(lsrottscavg/ltsPAavg)*1d3,
& rsqrt((unclstsc/ltsPAavg)**2+
& (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/
& rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3
write(6,712) dasin(rsrottscavg/rtsPAavg)*1d3-
& dasin(lsrottscavg/ltsPAavg)*1d3,
& rsqrt( ((uncrstsc/rtsPAavg)**2+(uncrstscPA*rsrottscavg/
& rtsPAavg**2)**2)/(1-(rsrottscavg/rtsPAavg)**2) +
& ((unclstsc/ltsPAavg)**2+(unclstscPA*lsrottscavg/
& ltsPAavg**2)**2)/(1-(lsrottscavg/ltsPAavg)**2) )*1d3
end if
c write(6,706) rsrottsc/nrts,
c & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)
c & *rSqrt(nrtssq)/nrts
c write(6,705) lsrottsc/nlts,
c & rsqrt((lsrottscsq-lsrottsc**2/nlts)/nlts)
c & *rSqrt(nltssq)/nrts
c write(6,712) rsrottsc/nrts - lsrottsc/nlts ,
c & rsqrt( (lsrottscsq-lsrottsc**2/nlts)/nlts
c & *nltssq/nlts**2
c & + (rsrottscsq-rsrottsc**2/nrts)/nrts
c & *nrtssq/nrts**2 )
write(6,*)
if ((nlgs.gt.0).and.(nrgs.gt.0)) then
argrgsc = (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs
argrgscPA = (nrgsPAsq-nrgsPA**2/nrgs)/nrgs
uncrsgsc = rsqrt(argrgsc)* rsqrt(nrgssq)/nrgs
uncrsgscPA = rsqrt(argrgscPA)*rsqrt(nrgssq)/nrgs
arglgsc = (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs
arglgscPA = (nlgsPAsq-nlgsPA**2/nlgs)/nlgs
unclsgsc = rsqrt(arglgsc)* rsqrt(nlgssq)/nlgs
unclsgscPA = rsqrt(arglgscPA)*rsqrt(nlgssq)/nlgs
write(6,7006) real(rgsPAavg*100),dasin(rsrotgscavg/rgsPAavg)*1d3,
& rsqrt((uncrsgsc/rgsPAavg)**2+
& (uncrsgscPA*rsrotgscavg/rgsPAavg**2)**2)/
& rsqrt(1-(rsrotgscavg/rgsPAavg)**2)*1d3
write(6,7005) real(lgsPAavg*100),dasin(lsrotgscavg/lgsPAavg)*1d3,
& rsqrt((unclsgsc/lgsPAavg)**2+
& (unclsgscPA*lsrotgscavg/lgsPAavg**2)**2)/
& rsqrt(1-(lsrotgscavg/lgsPAavg)**2)*1d3
write(6,7012) dasin(rsrotgscavg/rgsPAavg)*1d3-
& dasin(lsrotgscavg/lgsPAavg)*1d3,
& rsqrt( ((uncrsgsc/rgsPAavg)**2+(uncrsgscPA*rsrotgscavg/
& rgsPAavg**2)**2)/(1-(rsrotgscavg/rgsPAavg)**2) +
& ((unclsgsc/lgsPAavg)**2+(unclsgscPA*lsrotgscavg/
& lgsPAavg**2)**2)/(1-(lsrotgscavg/lgsPAavg)**2) )*1d3
end if
write(6,*)
if ((nlxs.gt.0).and.(nrxs.gt.0)) then
argrxsc = (rsrotxscsq-rsrotxsc**2/nrxs)/nrxs
argrxscPA = (nrxsPAsq-nrxsPA**2/nrxs)/nrxs
uncrsxsc = rsqrt(argrxsc)* rsqrt(nrxssq)/nrxs
uncrsxscPA = rsqrt(argrxscPA)*rsqrt(nrxssq)/nrxs
arglxsc = (lsrotxscsq-lsrotxsc**2/nlxs)/nlxs
arglxscPA = (nlxsPAsq-nlxsPA**2/nlxs)/nlxs
unclsxsc = rsqrt(arglxsc)* rsqrt(nlxssq)/nlxs
unclsxscPA = rsqrt(arglxscPA)*rsqrt(nlxssq)/nlxs
write(6,7008) real(rxsPAavg*100),dasin(rsrotxscavg/rxsPAavg)*1d3,
& rsqrt((uncrsxsc/rxsPAavg)**2+
& (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/
& rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3
write(6,7007) real(lxsPAavg*100),dasin(lsrotxscavg/lxsPAavg)*1d3,
& rsqrt((unclsxsc/lxsPAavg)**2+
& (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/
& rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3
write(6,7013) dasin(rsrotxscavg/rxsPAavg)*1d3-
& dasin(lsrotxscavg/lxsPAavg)*1d3,
& rsqrt( ((uncrsxsc/rxsPAavg)**2+(uncrsxscPA*rsrotxscavg/
& rxsPAavg**2)**2)/(1-(rsrotxscavg/rxsPAavg)**2) +
& ((unclsxsc/lxsPAavg)**2+(unclsxscPA*lsrotxscavg/
& lxsPAavg**2)**2)/(1-(lsrotxscavg/lxsPAavg)**2) )*1d3
end if
c write(6,7006) rsrotgsc/nrgs,
c & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)
c & *rSqrt(nrgssq)/nrgs
c write(6,7005) lsrotgsc/nlgs
c & rsqrt((lsrotgscsq-lsrotgsc**2/nlgs)/nlgs)
c & *rSqrt(nlgssq)/nrgs
c write(6,7012) rsrotgsc/nrgs - lsrotgsc/nlgs ,
c & rsqrt( (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs
c & *nlgssq/nlgs**2
c & + (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs
c & *nrgssq/nrgs**2 )
cccccccc estimated PNC angle scaling TScat and Xscat into average of ALL
if (nlxs.gt.0) then
rotnsw = dasin(rsrotnscavg/rsnscPAavg)*1d3
urotnsw = rsqrt((uncrsnsc/rsnscPAavg)**2+
& (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/
& rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3
rotnse =dasin(lsrotnscavg/lsnscPAavg)*1d3
urotnse = rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm
& (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/
& rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3
rotsw = dasin(rsrottscavg/rtsPAavg)*1d3
urotsw = rsqrt((uncrstsc/rtsPAavg)**2+
& (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/
& rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3
rotse = dasin(lsrottscavg/ltsPAavg)*1d3
urotse = rsqrt((unclstsc/ltsPAavg)**2+
& (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/
& rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3
rotxw = dasin(rsrotxscavg/rxsPAavg)*1d3
urotxw = rsqrt((uncrsxsc/rxsPAavg)**2+
& (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/
& rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3
rotxe = dasin(lsrotxscavg/lxsPAavg)*1d3
urotxe = rsqrt((unclsxsc/lxsPAavg)**2+
& (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/
& rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3
write(6,*)
write(6,*)
write(6,889) 0.5*( (fse-fsw)*0.5*(rotnsw+rotnse) +
& fsw*rotsw-fse*rotse ),
& 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 +
& (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 +
& (0.5*(fse-fsw)*urotnse)**2 + (0.5*(fse-fsw)*urotnsw)**2 +
& (fsw*urotsw)**2 + (fse*urotse)**2 )
write(6,*)
write(6,888) 0.5*( (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse) +
& fsw*rotsw-fse*rotse + fxw*rotxw-fxe*rotxe ),
& 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 +
& (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 +
& (0.5*(fse-fsw+fxe-fxw)*urotnse)**2 +
& (0.5*(fse-fsw+fxe-fxw)*urotnsw)**2 +
& (fsw*urotsw)**2 +(fse*urotse)**2 +
& (fxw*urotxw)**2 +(fxe*urotxe)**2 +
& (rotxw*ufxw)**2 +(rotxe*ufxe)**2 )
write(6,*)
write(6,890) (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse),
& rsqrt( (ufsw**2+ufse**2+ufxw**2+ufxe**2)*
& (0.5*(rotnsw+rotnse))**2
& + ((fse-fsw+fxe-fxw)*0.5)**2*(urotnsw**2+urotnse**2) )
write(6,893) (fse-fsw)*0.5*(rotnsw+rotnse),
& rsqrt( (ufsw**2+ufse**2)*(0.5*(rotnsw+rotnse))**2
& + ((fse-fsw)*0.5)**2*(urotnsw**2+urotnse**2)**2 )
write(6,891) fsw*rotsw-fse*rotse ,
& rsqrt( (ufsw*rotsw)**2 + (fsw*urotsw)**2
& + (ufse*rotse)**2 + (fse*urotse)**2 )
write(6,892) fxw*rotxw-fxe*rotxe ,
& rsqrt( (ufxw*rotxw)**2 + (fxw*urotxw)**2
& + (ufxe*rotxe)**2 + (fxe*urotxe)**2 )
write(6,*)
end if
write(6,*) ' ***************************************************',
&'******************'
write(6,*)
if (nder.gt.0) then
write(6,702) rsde/nder,
& rsqrt((rsdesq-rsde**2/nder)/nder)*rSqrt(ndersq)/nder
write(6,701) lsde/ndel,
& rsqrt((lsdesq-lsde**2/ndel)/ndel)*rSqrt(ndelsq)/ndel
end if
if ((nrts.gt.1).and.(nlts.gt.1)) then
write(6,802) rwde/nrts,
& rsqrt((rwdesq-rwde**2/nrts)/nrts)*rsqrt(nrtssq)/nrts
write(6,801) lwde/nlts,
& rsqrt((lwdesq-lwde**2/nlts)/nlts)*rsqrt(nltssq)/nlts
write(6,824) rwde/nrts-lwde/nlts,
& rsqrt( (lwdesq-lwde**2/nlts)/nlts*nltssq/nlts**2 +
& (rwdesq-rwde**2/nrts)/nrts*nrtssq/nrts**2 )
write(6,*)
write(6,822) rpl/nrts,
& rsqrt((rplsq-rpl**2/nrts)/nrts)*rsqrt(nrtssq)/nrts
write(6,821) lpl/nlts,
& rsqrt((lplsq-lpl**2/nlts)/nlts)*rsqrt(nltssq)/nlts
write(6,823) rpl/nrts-lpl/nlts,
& rsqrt( (lplsq-lpl**2/nlts)/nlts*nltssq/nlts**2 +
& (rplsq-rpl**2/nrts)/nrts*nrtssq/nrts**2 )
end if
write(6,*)
write(6,825) rplns/nrnsc,
& rsqrt((rplnssq-rplns**2/nrnsc)/nrnsc)*rSqrt(nrnscsq)/nrnsc
write(6,826) lplns/nlnsc,
& rsqrt((lplnssq-lplns**2/nlnsc)/nlnsc)*rSqrt(nlnscsq)/nlnsc
write(6,827) rplns/nrnsc-lplns/nlnsc,
& rsqrt( (rplnssq-rplns**2/nrnsc)/nrnsc*nrnscsq/nrnsc**2 +
& (lplnssq-lplns**2/nlnsc)/nlnsc*nlnscsq/nlnsc**2 )
write(6,*)
write(6,715) nrts, nlts
write(6,716) nrts/nlts
write(6,718) nrts-nlts,rSqrt(nltssq+nrtssq)
write(6,782) nrtsraw,nltsraw
write(6,*)
write(6,7015) nrgs, nlgs
write(6,7016) nrgs/nlgs
write(6,7018) nrgs-nlgs,rSqrt(nlgssq+nrgssq)
write(6,7082) nrgsraw,nlgsraw
write(6,*)
write(6,7020) nrxs, nlxs
write(6,7021) nrxs/nlxs
write(6,7022) nrxs-nlxs,rsqrt(nlxssq+nrxssq)
write(6,7023) nrxsraw,nlxsraw
write(3,*)
if ((nrsc.gt.0).and.(nlsc.gt.0)) then
write(6,815) nrsc, nlsc
write(6,816) nrsc/nlsc
write(6,818) nrsc-nlsc,sqrt(float(nlsc+nrsc))
else
write(6,*) ' NO scattering so limited output written '
endif
c 700 format(3x,' Rotations for Incoil Ouput to Outcoil input 114.6~cm with
c & B(mGauss)=',1f9.4)
700 format(3x,' PSMout to Det ',f9.4,' ~cm ',
& ' Bt(mG)=',1f6.2,'Bb(mG)=',1f6.2,' BzGradNum ',i2)
701 format(3x,' Avg E loss Det East Dn All Scat (meV)',1pe11.4,' +- '
& ,1pe11.4)
702 format(3x,' Avg E loss Det West Up All Scat (meV)',1pe11.4,' +- '
& ,1pe11.4)
801 format(3x,' Avg E loss Det EDn Target Scat (meV)',1pe11.4,' +- '
& ,1pe11.4)
802 format(3x,' Avg E loss Det WUp Target Scat (meV)',1pe11.4,' +- '
& ,1pe11.4)
821 format(3x,' Avg PathLenTar Det EDn Tar Scat (cm)',1pe15.8,' +- '
& ,1pe11.4)
822 format(3x,' Avg PathLenTar Det WUp Tar Scat (cm)',1pe15.8,' +- '
& ,1pe11.4)
823 format(3x,' Avg PathLenTar Det TS [WUp-EDn] (cm)',1pe15.8,' +- '
& ,1pe11.4)
824 format(3x,' Avg E Loss Det TS [WUp-EDn] (meV)',1pe15.8,' +- '
& ,1pe11.4)
825 format(3x,' Avg PathLenTar Det EDn NoScat (cm)',1pe15.8,' +- '
& ,1pe11.4)
826 format(3x,' Avg PathLenTar Det WUp NoScat (cm)',1pe15.8,' +- '
& ,1pe11.4)
827 format(3x,' Avg PathLenTar Det NS [WUp-EDn] (cm)',1pe15.8,' +- '
& ,1pe11.4)
703 format(3x,' PA ',f5.1,' East Dn (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
704 format(3x,' PA ',f5.1,' West Up (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
803 format(3x,' PA ',f5.1,' Rot East Dn (alt wgt) (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
804 format(3x,' PA ',f5.1,' Rot West Up (alt wgt) (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
705 format(3x,' PA ',f5.1,' Target Scat East Dn (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
706 format(3x,' PA ',f5.1,' Target Scat West Up (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
805 format(3x,' PA ',f5.1,' Scat East Dn (alt wgt) (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
806 format(3x,' PA ',f5.1,' Scat West Up (alt wgt) (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
7005 format(3x,' PA ',f5.1,' Gas Scat East Dn (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
7006 format(3x,' PA ',f5.1,' Gas Scat West Up (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
7007 format(3x,' PA ',f5.1,' Xover Tar Scat East Dn (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
7008 format(3x,' PA ',f5.1,' Xover Tar Scat West Up (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
725 format(3x,' PA ',f5.1,' NOScat East Dn (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
726 format(3x,' PA ',f5.1,' NOScat West Up (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
c 707 format(3x,' Avg. Rot Scat East Dn E<0.35meV (mrad) ',1pe11.4,
c & ' +- ',1pe11.4)
c 708 format(3x,' Avg. Rot Scat West Up E<0.35meV (mrad) ',1pe11.4,
c & ' +- ',1pe11.4)
c 709 format(3x,' Frac of rotation from scat East Dn ',1pe11.4,
c & ', oflow En ',1pe11.4)
c 710 format(3x,' Frac of rotation from scat West Up ',1pe11.4,
c & ', oflow En ',1pe11.4)
711 format(3x,' UpW-DnE (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
811 format(3x,' UpW-DnE (alt wgt) (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
712 format(3x,' UpW-DnE Tar scat (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
812 format(3x,' UpW-DnE Scat (alt wgt) (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
7012 format(3x,' UpW-DnE Gas scat (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
7013 format(3x,' UpW-DnE Xov Tar Scat (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
888 format(3x,' Estimated PNC angle from TS and XTS (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
889 format(3x,' Estimated PNC angle from TS only (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
890 format(3x,' (fse-fsw+fxe-fxw)*Theta-ns-avg (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
893 format(3x,' (fse-fsw)*Theta-ns-avg (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
891 format(3x,' fsw*thsw-fse*thse (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
892 format(3x,' fxw*thxw-fxe*thxe (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
732 format(3x,' UpW-DnE NOScat (mrad) ',1pe13.6,
& ' +- ',1pe9.2)
713 format(3x,' # in DetW Up',1pe13.6,2x,', # in DetE Dn ',1pe13.6)
714 format(3x,' Ratio # DetWUp/DetEDn ,'1pe13.6)
717 format(3x,' Det Counts UpW-DnE ',1pe11.4,' +- ',1pe11.4)
715 format(3x,' # Tarscat in DetW Up ',1pe11.4,2x,
& ', # Tarscat in Det E Dn ',1pe11.4)
716 format(3x,' Ratio # Tarscat DetWUp/DetEDn ,'1pe11.4)
718 format(3x,' TarScat Det Counts UpW-DnE ',1pe11.4,' +- ',1pe11.4)
815 format(3x,' AltWgt # Tarscat in DetW Up ',1pe11.4,2x,
& ', # Tarscat in Det E Dn ',1pe11.4)
816 format(3x,' AltWgt Ratio # Tarscat DetWUp/DetEDn ,'1pe11.4)
818 format(3x,' AltWgt TarScat Det Counts UpW-DnE ',1pe11.4,' +- '
& ,1pe11.4)
7015 format(2x,'# Gas/Al scat DetW Up ',1pe11.4,2x,
& ',# Gas/Al scat DetE Dn ',1pe11.4)
7016 format(3x,' Ratio # Gas/Al Scat DetWUp/DetEDn ,'1pe11.4)
7018 format(3x,' G/Al Scat Det Counts UpW-DnE ',1pe11.4,' +- ',
& 1pe11.4)
7020 format(2x,'# XovTarScat scat DetW Up ',1pe11.4,2x,
& ',# XovTarScat scat DetE Dn ',1pe11.4)
7021 format(3x,' Ratio # XovTarScat DetWUp/DetEDn ,'1pe11.4)
7022 format(3x,' XovTarScat Det Counts UpW-DnE ',1pe11.4,' +- ',
& 1pe11.4)
719 format(3x,'# scat w/ E<0.35meV at det West: ',1pe11.4,
& ' East: ',1pe11.4)
755 format(1x,' Frac # of scatters with Ef.ne.Ei tar1 ',1pe11.4)
756 format(1x,' Frac # of scatters with Ef.ne.Ei tar2 ',1pe11.4)
757 format(1x,' Frac num of tar1 entering tar2 ',1pe11.4)
758 format(1x,' Frac # of total reaching det ',1pe11.4)
759 format(1x,' Frac det hits TarScat W ',1pe11.4,' +- ',1pe11.4)
760 format(1x,' Frac det hits TarScat E ',1pe11.4,' +- ',1pe11.4)
761 format(1x,' Frac det hits G/Al Scat W ',1pe11.4,' +- ',1pe11.4)
762 format(1x,' Frac det hits G/Al Scat E ',1pe11.4,' +- ',1pe11.4)
763 format(1x,' Frac det hits TXall W ',1pe11.4,' +- ',1pe11.4)
764 format(1x,' Frac det hits TXall E ',1pe11.4,' +- ',1pe11.4)
765 format(1x,' Frac det hits TXScat W ',1pe11.4,' +- ',1pe11.4)
766 format(1x,' Frac det hits TXScat E ',1pe11.4,' +- ',1pe11.4)
767 format(1x,' Frac det hits NO Scat W ',1pe11.4,' +- ',1pe11.4)
768 format(1x,' Frac det hits NO Scat E ',1pe11.4,' +- ',1pe11.4)
770 format(3x,' num NO scat -- det W: ',1pe13.6,
& ' det E: ',1pe13.6)
780 format(3x,'Total No. of cross-over neutrons: ',1pe11.4)
782 format(3x,'Unwt # Tar scat DetW Up ',1pe11.4,2x,
& ', Unwt # Tar scat Det E Dn ',1pe11.4)
7082 format(3x,'UnWt # Gas scat DetW Up ',1pe11.4,2x,
& ', Unwt # Gas scat Det E Dn ',1pe11.4)
7023 format(3x,'UnWt # XovTarScat DetW Up ',1pe11.4,2x,
& ', Unwt # XovTarScat Det E Dn ',1pe11.4)
735 format(3x,' Guide output W ',1pe11.4,' Guide output E ',1pe11.4)
736 format(3x,' Incoil output W ',1pe11.4,' Incoil output E ',1pe11.4)
737 format(3x,' Outcoil output W ',1pe11.4,
& ' Outcoil output E ',1pe11.4)
c
write(3,*) ' End time: ', date, time
write(3,*) ' HBOOK filename: ',outfname
write(3,*)
c*********averages -- Note: Now with weights, nr, nl, etc are all sums of
c******* weights, see logbook (bec)
write(3,755) Nenscat1/Ntar1
write(3,756) Nenscat2/Ntar2
write(3,757) Ntar12/Ntar1
write(3,758) cot/float(Nevts)
write(3,759) fsw,ufsw
write(3,760) fse,ufse
write(3,761) nrgs/nr,nrgs/nr*dsqrt(1/nrgsraw+1/nr)
write(3,762) nlgs/nl,nlgs/nl*dsqrt(1/nlgsraw+1/nl)
write(3,765) fxw,ufxw
write(3,766) fxe,ufxe
write(3,767) nrnsc/nr,nrnsc/nr*dsqrt(1/nrnsc+1/nr)
write(3,768) nlnsc/nl,nlnsc/nl*dsqrt(1/nlnsc+1/nl)
write(3,*)
write(3,*) ' StDev All West ',
& rsqrt((rsrotsq-rsrot**2/nr)/nr)*1d3
write(3,*) ' St Err All West ',
& rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr*1d3
c write(3,*) ' nr ',nr , ' nrsq ',nrsq
c write(3,*) 'rsqrt(nrsq)/nr ', rSqrt(nrsq)/nr
c
write(3,*) ' StDev NoScat West ',
& rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*1d3
write(3,*) ' St Err NoScat West ',
& rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*
& rSqrt(nrnscsq)/nrnsc*1d3
c write(3,*) ' nr ',nr , ' nrsq ',nrsq
c write(3,*) 'rsqrt(nrsq)/nr ', rSqrt(nrsq)/nr
c
if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then
c write(3,*)
write(3,*) ' StDev Tar Scat West ',
& rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1d3
write(3,*) ' St Err Tar Scat West ',
& rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*1d3
& *rSqrt(nrtssq)/nrts
c write(3,*) ' nrts ',nrts , ' nrtssq ',nrtssq
c write(3,*) 'rsqrt(nrtssq)/nrts ', rSqrt(nrtssq)/nrts
end if
c write(3,*)
write(3,*) ' StDev XOvTar Scat West ',
& rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1d3
write(3,*) ' St Err XOvTar Scat West ',
& rsqrt((rsrotxscsq-rsrotxsc**2/nrxs)/nrxs)*1d3
& *rSqrt(nrxssq)/nrxs
c write(3,*) ' nrxs ',nrxs , ' nrxssq ',nrxssq
c write(3,*) 'rsqrt(nrxssq)/nrxs ', rSqrt(nrxssq)/nrxs
if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then
c write(3,*)
write(3,*) ' StDev Gas Scat West ',
& rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1d3
write(3,*) ' St Err Gas Scat West ',
& rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*1d3
& *rSqrt(nrgssq)/nrgs
c write(3,*) ' nrgs ',nrgs , ' nrgssq ',nrgssq
c write(3,*) 'rsqrt(nrgssq)/nrgs ', rSqrt(nrgssq)/nrgs
end if
write(3,*)
write(3,*)
write(3,713) nr, nl
write(3,714) nr/nl
write(3,717) nr-nl,rSqrt(nlsq+nrsq)
write(3,*)
write(3,735) NguideoutW,NguideoutE
write(3,736) NincoiloutW,NincoiloutE
write(3,737) NoutcoiloutW,NoutcoiloutE
write(3,*)
write(3,770) nrnsc , nlnsc
write(3,*)
write(3,700) zdet, mag_B,mag_B_beam, bzgradnum
write(3,*)
write(3,*) ' ********************** Average Rotation Values
&***********************'
write(3,*)
argrsrot = (rsrotsq-rsrot**2/nr)/nr
argrsrotPA = (nrPAsq-nrPA**2/nr)/nr
c if (argrsrot.lt.0) then
c write(3,*) 'argrsrot',argrsrot
c argrsrot =0
c endif
c if (argrsrotPA.lt.0) then
c write(3,*) 'argrsrotPA',argrsrotPA
c argrsrotPA =0
c endif
uncrsrotavg = rsqrt(argrsrot)*rsqrt(nrsq)/nr
uncrPAavg = rsqrt(argrsrotPA)*rsqrt(nrsq)/nr
arglsrot = (lsrotsq-lsrot**2/nl)/nl
arglsrotPA = (nlPAsq-nlPA**2/nl)/nl
c if (arglsrot.lt.0) then
c write(3,*) 'arglsrot',arglsrot
c arglsrot =0
c endif
c if (arglsrotPA.lt.0) then
c write(3,*) 'arglsrotPA',arglsrotPA
c arglsrotPA =0
c endif
unclsrotavg = rsqrt(arglsrot)*rsqrt(nlsq)/nl
unclPAavg = rsqrt(arglsrotPA)*rsqrt(nlsq)/nl
write(3,704) real(rPAavg*100),dasin(rsrotavg/rPAavg)*1d3, !rcm
& rsqrt((uncrsrotavg/rPAavg)**2+
& (uncrPAavg*rsrotavg/rPAavg**2)**2)/
& rsqrt(1-(rsrotavg/rPAavg)**2)*1d3
write(3,703) real(lPAavg*100),dasin(lsrotavg/lPAavg)*1d3,
& rsqrt((unclsrotavg/lPAavg)**2+
& (unclPAavg*lsrotavg/lPAavg**2)**2)/
& rsqrt(1-(lsrotavg/lPAavg)**2)*1d3
write(3,711) dasin(rsrotavg/rPAavg)*1d3-
& dasin(lsrotavg/lPAavg)*1d3,
& rsqrt( ((uncrsrotavg/rPAavg)**2+
& (uncrPAavg*rsrotavg/rPAavg**2)**2)/
& (1-(rsrotavg/rPAavg)**2) + ((unclsrotavg/lPAavg)**2+
& (unclPAavg*lsrotavg/lPAavg**2)**2)/
& (1-(lsrotavg/lPAavg)**2))*1d3
c write(3,704) rsrot/nr,
c & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr
c write(3,703) lsrot/nl,
c & rsqrt((lsrotsq-lsrot**2/nl)/nl)*rSqrt(nlsq)/nl
c write(3,711) rsrot/nr-lsrot/nl,
c & rsqrt( (lsrotsq-lsrot**2/nl)/nl*nlsq/nl**2
c & + (rsrotsq-rsrot**2/nr)/nr*nrsq/nr**2 )
c if ((nrall.gt.1).and.(nlall.gt.1)) then
c write(3,*)
c write(3,804) rall/float(nrall), ! alternate weighting method
c & rsqrt((rallsq-rall**2/nrall)/nrall/(nrall-1))
c write(3,803) lall/nlall,
c & rsqrt((lallsq-lall**2/nlall)/nlall/(nlall-1))
c write(3,811) rall/nrall-lall/nlall,
c & rsqrt( (lallsq-lall**2/nlall)/nlall/(nlall-1)
c & + (rallsq-rall**2/nrall)/nrall/(nlall-1) )
c else
c write(3,*) 'No alternate weighting data '
c end if
uncrallavg = rsqrt((rallsq-rall**2/nrall)/ !rcm PAsin
& (nrall-1)/nrall)
uncrallPAavg = rsqrt((nrallPAsq - nrallPA**2/nrall)/
& (nrall-1)/nrall)
unclallavg = rsqrt((lallsq-lall**2/nlall)/
& (nlall-1)/nlall)
unclallPAavg = rsqrt((nlallPAsq - nlallPA**2/nlall)/
& (nlall-1)/nlall)
write(3,*)
if ((nrall.gt.1).and.(nlall.gt.1)) then
write(3,804) real(rallPAavg*100),dasin(rallavg/rallPAavg)*1d3, ! alternate weighting method
& rsqrt( (uncrallavg/rallPAavg)**2+
& (uncrallPAavg*rallavg/rallPAavg**2)**2 )/
& rsqrt(1-(rallavg/rallPAavg)**2)*1d3
write(3,803) real(lallPAavg*100),dasin(lallavg/lallPAavg)*1d3,
& rsqrt( (unclallavg/lallPAavg)**2+
& (unclallPAavg*lallavg/lallPAavg**2)**2 )/
& rsqrt(1-(lallavg/lallPAavg)**2)*1d3
write(3,811) dasin(rallavg/rallPAavg)*1d3-
& dasin(lallavg/lallPAavg)*1d3,
& rsqrt( ((uncrallavg/rallPAavg)**2+
& (uncrallPAavg*rallavg/rallPAavg**2)**2)/
& (1-(rallavg/rallPAavg)**2) + ((unclallavg/lallPAavg)**2+
& (unclallPAavg*lallavg/lallPAavg**2)**2)/
& (1-(lallavg/lallPAavg)**2))*1d3
else
write(3,*) 'No alternate weighting data '
end if
write(3,*)
argrsnsc = (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc
argrsnscPA = (nrnscPAsq-nrnscPA**2/nrnsc)/nrnsc
c if (argrsnsc.lt.0) then
c write(3,*) 'argrsnsc',argrsnsc
c argrsnsc =0
c endif
c if (argrsnscPA.lt.0) then
c write(3,*) 'argrsnscPA',argrsnscPA
c argrsnscPA =0
c endif
uncrsnsc = rsqrt(argrsnsc)*rsqrt(nrnscsq)/nrnsc
uncrsnscPA = rsqrt(argrsnscPA)*rsqrt(nrnscsq)/nrnsc
arglsnsc = (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc
arglsnscPA = (nlnscPAsq-nlnscPA**2/nlnsc)/nlnsc
c if (arglsnsc.lt.0) then
c write(3,*) 'arglsnsc',arglsnsc
c arglsnsc =0
c endif
c if (arglsnscPA.lt.0) then
c write(3,*) 'arglsnscPA',arglsnscPA
c arglsnscPA =0
c endif
unclsnsc = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc
unclsnscPA = rsqrt(arglsnsc)*rsqrt(nlnscsq)/nlnsc
write(3,726) real(rsnscPAavg*100),
& dasin(rsrotnscavg/rsnscPAavg)*1d3,
& rsqrt((uncrsnsc/rsnscPAavg)**2+
& (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/
& rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3
write(3,725) real(lsnscPAavg*100),
& dasin(lsrotnscavg/lsnscPAavg)*1d3,
& rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm
& (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/
& rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3
write(3,732)dasin(rsrotnscavg/rsnscPAavg)*1d3-
& dasin(lsrotnscavg/lsnscPAavg)*1d3,
& rsqrt( ((uncrsnsc/rsnscPAavg)**2+
& (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/
& (1-(rsrotnscavg/rsnscPAavg)**2) + ((unclsnsc/lsnscPAavg)**2+
& (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/
& (1-(lsrotnscavg/lsnscPAavg)**2) )*1d3
c write(3,726) rsrotnsc/nrnsc, !rcm
c & rsqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)
c & *rSqrt(nrnscsq)/nrnsc
c write(3,725) lsrotnsc/nlnsc,
c & rsqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc)
c & *rSqrt(nlnscsq)/nlnsc
c write(3,732) rsrotnsc/nrnsc - lsrotnsc/nlnsc,
c & rsqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc
c & *nlnscsq/nlnsc**2
c & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc
c & *nrnscsq/nrnsc**2)
write(3,*)
if ((nrsc.gt.1).and.(nlsc.gt.1)) then !rcm
uncrscavg = rsqrt((rscsq-rsc**2/nrsc)/
& (nrsc-1)/nrsc)
uncrscPAavg = rsqrt((nrscPAsq - nrscPA**2/nrsc)/
& (nrsc-1)/nrsc)
unclscavg = rsqrt((lscsq-lsc**2/nlsc)/
& (nlsc-1)/nlsc)
unclscPAavg = rsqrt((nlscPAsq - nlscPA**2/nlsc)/
& (nlsc-1)/nlsc)
write(3,806) real(rscPAavg*100),dasin(rscavg/rscPAavg)*1d3, ! alternate weighting method
& rsqrt( (uncrscavg/rscPAavg)**2+
& (uncrscPAavg*rscavg/rscPAavg**2)**2 )/
& rsqrt(1-(rscavg/rscPAavg)**2)*1d3
write(3,805) real(lscPAavg*100),dasin(lscavg/lscPAavg)*1d3,
& rsqrt( (unclscavg/lscPAavg)**2+
& (unclscPAavg*lscavg/lscPAavg**2)**2 )/
& rsqrt(1-(lscavg/lscPAavg)**2)*1d3
write(3,812) dasin(rscavg/rscPAavg)*1d3-dasin(lscavg/lscPAavg)*1d3,
& rsqrt( ((uncrscavg/rscPAavg)**2+(uncrscPAavg*rscavg/
& rscPAavg**2)**2)/(1-(rscavg/rscPAavg)**2) +
& ((unclscavg/lscPAavg)**2+(unclscPAavg*lscavg/
& lscPAavg**2)**2)/(1-(lscavg/lscPAavg)**2) )*1d3
else
write(3,*) ' No alternate weighting scatters'
end if
write(3,*)
if ((nlts.gt.0).and.(nrts.gt.0)) then
argrtsc = (rsrottscsq-rsrottsc**2/nrts)/nrts
argrtscPA = (nrtsPAsq-nrtsPA**2/nrts)/nrts
uncrstsc = rsqrt(argrtsc)*rsqrt(nrtssq)/nrts
uncrstscPA = rsqrt(argrtscPA)*rsqrt(nrtssq)/nrts
argltsc = (lsrottscsq-lsrottsc**2/nlts)/nlts
argltscPA = (nltsPAsq-nltsPA**2/nlts)/nlts
unclstsc = rsqrt(argltsc)*rsqrt(nltssq)/nlts
unclstscPA = rsqrt(argltscPA)*rsqrt(nltssq)/nlts
write(3,706) real(rtsPAavg*100),dasin(rsrottscavg/rtsPAavg)*1d3,
& rsqrt((uncrstsc/rtsPAavg)**2+
& (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/
& rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3
write(3,705) real(ltsPAavg*100),dasin(lsrottscavg/ltsPAavg)*1d3,
& rsqrt((unclstsc/ltsPAavg)**2+
& (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/
& rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3
write(3,712) dasin(rsrottscavg/rtsPAavg)*1d3-
& dasin(lsrottscavg/ltsPAavg)*1d3,
& rsqrt( ((uncrstsc/rtsPAavg)**2+(uncrstscPA*rsrottscavg/
& rtsPAavg**2)**2)/(1-(rsrottscavg/rtsPAavg)**2) +
& ((unclstsc/ltsPAavg)**2+(unclstscPA*lsrottscavg/
& ltsPAavg**2)**2)/(1-(lsrottscavg/ltsPAavg)**2) )*1d3
end if
c write(3,706) rsrottsc/nrts,
c & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)
c & *rSqrt(nrtssq)/nrts
c write(3,705) lsrottsc/nlts,
c & rsqrt((lsrottscsq-lsrottsc**2/nlts)/nlts)
c & *rSqrt(nltssq)/nrts
c write(3,712) rsrottsc/nrts - lsrottsc/nlts ,
c & rsqrt( (lsrottscsq-lsrottsc**2/nlts)/nlts
c & *nltssq/nlts**2
c & + (rsrottscsq-rsrottsc**2/nrts)/nrts
c & *nrtssq/nrts**2 )
write(3,*)
if ((nlgs.gt.0).and.(nrgs.gt.0)) then
argrgsc = (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs
argrgscPA = (nrgsPAsq-nrgsPA**2/nrgs)/nrgs
uncrsgsc = rsqrt(argrgsc)* rsqrt(nrgssq)/nrgs
uncrsgscPA = rsqrt(argrgscPA)*rsqrt(nrgssq)/nrgs
arglgsc = (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs
arglgscPA = (nlgsPAsq-nlgsPA**2/nlgs)/nlgs
unclsgsc = rsqrt(arglgsc)* rsqrt(nlgssq)/nlgs
unclsgscPA = rsqrt(arglgscPA)*rsqrt(nlgssq)/nlgs
write(3,7006) real(rgsPAavg*100),dasin(rsrotgscavg/rgsPAavg)*1d3,
& rsqrt((uncrsgsc/rgsPAavg)**2+
& (uncrsgscPA*rsrotgscavg/rgsPAavg**2)**2)/
& rsqrt(1-(rsrotgscavg/rgsPAavg)**2)*1d3
write(3,7005) real(lgsPAavg*100),dasin(lsrotgscavg/lgsPAavg)*1d3,
& rsqrt((unclsgsc/lgsPAavg)**2+
& (unclsgscPA*lsrotgscavg/lgsPAavg**2)**2)/
& rsqrt(1-(lsrotgscavg/lgsPAavg)**2)*1d3
write(3,7012) dasin(rsrotgscavg/rgsPAavg)*1d3-
& dasin(lsrotgscavg/lgsPAavg)*1d3,
& rsqrt( ((uncrsgsc/rgsPAavg)**2+(uncrsgscPA*rsrotgscavg/
& rgsPAavg**2)**2)/(1-(rsrotgscavg/rgsPAavg)**2) +
& ((unclsgsc/lgsPAavg)**2+(unclsgscPA*lsrotgscavg/
& lgsPAavg**2)**2)/(1-(lsrotgscavg/lgsPAavg)**2) )*1d3
end if
write(3,*)
if ((nlxs.gt.0).and.(nrxs.gt.0)) then
argrxsc = (rsrotxscsq-rsrotxsc**2/nrxs)/nrxs
argrxscPA = (nrxsPAsq-nrxsPA**2/nrxs)/nrxs
uncrsxsc = rsqrt(argrxsc)* rsqrt(nrxssq)/nrxs
uncrsxscPA = rsqrt(argrxscPA)*rsqrt(nrxssq)/nrxs
arglxsc = (lsrotxscsq-lsrotxsc**2/nlxs)/nlxs
arglxscPA = (nlxsPAsq-nlxsPA**2/nlxs)/nlxs
unclsxsc = rsqrt(arglxsc)* rsqrt(nlxssq)/nlxs
unclsxscPA = rsqrt(arglxscPA)*rsqrt(nlxssq)/nlxs
write(3,7008) real(rxsPAavg*100),dasin(rsrotxscavg/rxsPAavg)*1d3,
& rsqrt((uncrsxsc/rxsPAavg)**2+
& (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/
& rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3
write(3,7007) real(lxsPAavg*100),dasin(lsrotxscavg/lxsPAavg)*1d3,
& rsqrt((unclsxsc/lxsPAavg)**2+
& (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/
& rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3
write(3,7013) dasin(rsrotxscavg/rxsPAavg)*1d3-
& dasin(lsrotxscavg/lxsPAavg)*1d3,
& rsqrt( ((uncrsxsc/rxsPAavg)**2+(uncrsxscPA*rsrotxscavg/
& rxsPAavg**2)**2)/(1-(rsrotxscavg/rxsPAavg)**2) +
& ((unclsxsc/lxsPAavg)**2+(unclsxscPA*lsrotxscavg/
& lxsPAavg**2)**2)/(1-(lsrotxscavg/lxsPAavg)**2) )*1d3
end if
c write(3,7006) rsrotgsc/nrgs,
c & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)
c & *rSqrt(nrgssq)/nrgs
c write(3,7005) lsrotgsc/nlgs
c & rsqrt((lsrotgscsq-lsrotgsc**2/nlgs)/nlgs)
c & *rSqrt(nlgssq)/nrgs
c write(3,7012) rsrotgsc/nrgs - lsrotgsc/nlgs ,
c & rsqrt( (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs
c & *nlgssq/nlgs**2
c & + (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs
c & *nrgssq/nrgs**2 )
cccccccc estimated PNC angle scaling TScat and Xscat into average of ALL
if (nlxs.gt.0) then
rotnsw = dasin(rsrotnscavg/rsnscPAavg)*1d3
urotnsw = rsqrt((uncrsnsc/rsnscPAavg)**2+
& (uncrsnscPA*rsrotnscavg/rsnscPAavg**2)**2)/
& rsqrt(1-(rsrotnscavg/rsnscPAavg)**2)*1d3
rotnse =dasin(lsrotnscavg/lsnscPAavg)*1d3
urotnse = rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm
& (unclsnscPA*lsrotnscavg/lsnscPAavg**2)**2)/
& rsqrt(1-(lsrotnscavg/lsnscPAavg)**2)*1d3
rotsw = dasin(rsrottscavg/rtsPAavg)*1d3
urotsw = rsqrt((uncrstsc/rtsPAavg)**2+
& (uncrstscPA*rsrottscavg/rtsPAavg**2)**2)/
& rsqrt(1-(rsrottscavg/rtsPAavg)**2)*1d3
rotse = dasin(lsrottscavg/ltsPAavg)*1d3
urotse = rsqrt((unclstsc/ltsPAavg)**2+
& (unclstscPA*lsrottscavg/ltsPAavg**2)**2)/
& rsqrt(1-(lsrottscavg/ltsPAavg)**2)*1d3
rotxw = dasin(rsrotxscavg/rxsPAavg)*1d3
urotxw = rsqrt((uncrsxsc/rxsPAavg)**2+
& (uncrsxscPA*rsrotxscavg/rxsPAavg**2)**2)/
& rsqrt(1-(rsrotxscavg/rxsPAavg)**2)*1d3
rotxe = dasin(lsrotxscavg/lxsPAavg)*1d3
urotxe = rsqrt((unclsxsc/lxsPAavg)**2+
& (unclsxscPA*lsrotxscavg/lxsPAavg**2)**2)/
& rsqrt(1-(lsrotxscavg/lxsPAavg)**2)*1d3
write(3,*)
write(3,*)
write(3,889) 0.5*( (fse-fsw)*0.5*(rotnsw+rotnse) +
& fsw*rotsw-fse*rotse ),
& 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 +
& (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 +
& (0.5*(fse-fsw)*urotnse)**2 + (0.5*(fse-fsw)*urotnsw)**2 +
& (fsw*urotsw)**2 + (fse*urotse)**2 )
write(3,*)
write(3,888) 0.5*( (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse) +
& fsw*rotsw-fse*rotse + fxw*rotxw-fxe*rotxe ),
& 0.5*rsqrt( ((0.5*(rotnsw+rotnse)-rotse)*ufse)**2 +
& (-(0.5*(rotnsw+rotnse)+rotsw)*ufsw)**2 +
& (0.5*(fse-fsw+fxe-fxw)*urotnse)**2 +
& (0.5*(fse-fsw+fxe-fxw)*urotnsw)**2 +
& (fsw*urotsw)**2 +(fse*urotse)**2 +
& (fxw*urotxw)**2 +(fxe*urotxe)**2 +
& (rotxw*ufxw)**2 +(rotxe*ufxe)**2 )
write(3,*)
write(3,890) (fse-fsw+fxe-fxw)*0.5*(rotnsw+rotnse),
& rsqrt( (ufsw**2+ufse**2+ufxw**2+ufxe**2)*
& (0.5*(rotnsw+rotnse))**2
& + ((fse-fsw+fxe-fxw)*0.5)**2*(urotnsw**2+urotnse**2) )
write(3,893) (fse-fsw)*0.5*(rotnsw+rotnse),
& rsqrt( (ufsw**2+ufse**2)*(0.5*(rotnsw+rotnse))**2
& + ((fse-fsw)*0.5)**2*(urotnsw**2+urotnse**2)**2 )
write(3,891) fsw*rotsw-fse*rotse ,
& rsqrt( (ufsw*rotsw)**2 + (fsw*urotsw)**2
& + (ufse*rotse)**2 + (fse*urotse)**2 )
write(3,892) fxw*rotxw-fxe*rotxe ,
& rsqrt( (ufxw*rotxw)**2 + (fxw*urotxw)**2
& + (ufxe*rotxe)**2 + (fxe*urotxe)**2 )
write(3,*)
end if
write(3,*) ' ***************************************************',
&'******************'
write(3,*)
if (nder.gt.0) then
write(3,702) rsde/nder,
& rsqrt((rsdesq-rsde**2/nder)/nder)*rSqrt(ndersq)/nder
write(3,701) lsde/ndel,
& rsqrt((lsdesq-lsde**2/ndel)/ndel)*rSqrt(ndelsq)/ndel
end if
if ((nrts.gt.1).and.(nlts.gt.1)) then
write(3,802) rwde/nrts,
& rsqrt((rwdesq-rwde**2/nrts)/nrts)*rsqrt(nrtssq)/nrts
write(3,801) lwde/nlts,
& rsqrt((lwdesq-lwde**2/nlts)/nlts)*rsqrt(nltssq)/nlts
write(3,824) rwde/nrts-lwde/nlts,
& rsqrt( (lwdesq-lwde**2/nlts)/nlts*nltssq/nlts**2 +
& (rwdesq-rwde**2/nrts)/nrts*nrtssq/nrts**2 )
write(3,*)
write(3,822) rpl/nrts,
& rsqrt((rplsq-rpl**2/nrts)/nrts)*rsqrt(nrtssq)/nrts
write(3,821) lpl/nlts,
& rsqrt((lplsq-lpl**2/nlts)/nlts)*rsqrt(nltssq)/nlts
write(3,823) rpl/nrts-lpl/nlts,
& rsqrt( (lplsq-lpl**2/nlts)/nlts*nltssq/nlts**2 +
& (rplsq-rpl**2/nrts)/nrts*nrtssq/nrts**2 )
end if
write(3,*)
write(3,*)
write(3,825) rplns/nrnsc,
& rsqrt((rplnssq-rplns**2/nrnsc)/nrnsc)*rSqrt(nrnscsq)/nrnsc
write(3,826) lplns/nlnsc,
& rsqrt((lplnssq-lplns**2/nlnsc)/nlnsc)*rSqrt(nlnscsq)/nlnsc
write(3,827) rplns/nrnsc-lplns/nlnsc,
& rsqrt( (rplnssq-rplns**2/nrnsc)/nrnsc*nrnscsq/nrnsc**2 +
& (lplnssq-lplns**2/nlnsc)/nlnsc*nlnscsq/nlnsc**2 )
write(3,*)
write(3,715) nrts, nlts
write(3,716) nrts/nlts
write(3,718) nrts-nlts,rSqrt(nltssq+nrtssq)
write(3,782) nrtsraw,nltsraw
write(3,*)
write(3,7015) nrgs, nlgs
write(3,7016) nrgs/nlgs
write(3,7018) nrgs-nlgs,rSqrt(nlgssq+nrgssq)
write(3,7082) nrgsraw,nlgsraw
write(3,*)
write(3,7020) nrxs, nlxs
write(3,7021) nrxs/nlxs
write(3,7022) nrxs-nlxs,rsqrt(nlxssq+nrxssq)
write(3,7023) nrxsraw,nlxsraw
write(3,*)
if ((nrsc.gt.0).and.(nlsc.gt.0)) then
write(3,815) nrsc, nlsc
write(3,816) nrsc/nlsc
write(3,818) nrsc-nlsc,sqrt(float(nlsc+nrsc))
else
write(3,*) ' NO scattering so limited output written '
endif
c---------
write(6,*);write(3,*)
write(6,*) ' Number that bounced ',numbounce
write(3,*) ' Number that bounced ',numbounce
write(6,*) ' Number that bounced in 2 guides ',numbounce2
write(3,*) ' Number that bounced in 2 guides ',numbounce2
write(6,*) ' Number that bounced in 3 guides ',numbounce3
write(3,*) ' Number that bounced in 3 guides ',numbounce3
write(6,*) ' Number crossed over after target area WGT: ',
& dncross
write(3,*) ' Number crossed over after target area WGT: ',
& dncross
write(3,780) Nxo
write(6,780) Nxo
write(6,*)
write(6,*) ' num guide1 bounced x,y 0x: ',ngbx0,ngby0
write(6,*) ' num guide1 bounced x,y 1x: ',ngbx1,ngby1
write(6,*) ' num guide1 bounced x,y 2x: ',ngbx2,ngby2
write(6,*) ' num guide1 bounced x,y 3x: ',ngbx3,ngby3
write(6,*) ' num guide1 scatters ',gsc
write(3,*)
write(3,*) ' num guide1 bounced x,y 0x: ',ngbx0,ngby0
write(3,*) ' num guide1 bounced x,y 1x: ',ngbx1,ngby1
write(3,*) ' num guide1 bounced x,y 2x: ',ngbx2,ngby2
write(3,*) ' num guide1 bounced x,y 3x: ',ngbx3,ngby3
write(3,*) ' num guide1 scatters ',gsc
write(6,*)
write(6,*) ' num incoil bounced x,y 0x: ',nibx0,niby0
write(6,*) ' num incoil bounced x,y 1x: ',nibx1,niby1
write(6,*) ' num incoil bounced x,y 2x: ',nibx2,niby2
write(6,*) ' num incoil bounced x,y 3x: ',nibx3,niby3
write(6,*) ' num incoil scatters ',isc
write(3,*)
write(3,*) ' num incoil bounced x,y 0x: ',nibx0,niby0
write(3,*) ' num incoil bounced x,y 1x: ',nibx1,niby1
write(3,*) ' num incoil bounced x,y 2x: ',nibx2,niby2
write(3,*) ' num incoil bounced x,y 3x: ',nibx3,niby3
write(3,*) ' num incoil scatters ',isc
write(6,*)
write(6,*) ' num outcoil bounced x,y 0x: ',nobx0,noby0
write(6,*) ' num outcoil bounced x,y 1x: ',nobx1,noby1
write(6,*) ' num outcoil bounced x,y 2x: ',nobx2,noby2
write(6,*) ' num outcoil bounced x,y 3x: ',nobx3,noby3
write(6,*) ' num outcoil scatters ',osc
write(3,*)
write(3,*) ' num outcoil bounced x,y 0x: ',nobx0,noby0
write(3,*) ' num outcoil bounced x,y 1x: ',nobx1,noby1
write(3,*) ' num outcoil bounced x,y 2x: ',nobx2,noby2
write(3,*) ' num outcoil bounced x,y 3x: ',nobx3,noby3
write(3,*) ' num outcoil scatters ',osc
write(6,*) ; write(3,*)
write(6,*) ' Fraction of det hits that bounced ',
& numbounce/cot
write(3,*) ' Fraction of det hits that bounced ',
& numbounce/cot
write(6,*) ' Fraction of det hits that bounced in 2 guides',
& numbounce2/cot
write(3,*) ' Fraction of det hits that bounced in 2 guides ',
& numbounce2/cot
write(6,*) ' Fraction of det hits that bounced in 3 guides ',
& numbounce3/cot
write(3,*) ' Fraction of det hits that bounced in 3 guides ',
& numbounce3/cot
write(6,*) ' Frac. of det hits that crossed after targets WGT ',
& dble(dncross)/cot
write(3,*) ' Frac. of det hits that crossed after targets WGT',
& dble(dncross)/cot
write(6,*) ' Frac. of det hits that crossed over anywhere ',
& Nxo/cot
write(3,*) ' Frac. of det hits that crossed over anywhere ',
& Nxo/cot
close(unit=3)
c bec return
end
c Functions to avoid domain errors -rcm
function rsqrt(B)
double precision rsqrt,B
nanflg = isnan(B)
if(B.lt.0.d0) then
write(17,*) ' argsqrt<0 ',B
write(17,*) ' argsqrt<0 ',B
rsqrt = 0
else if (nanflg.eq..true.) then
B = 0
x1 = -1d6; y1 = -1d6
else
rsqrt = sqrt(B)
end if
end
function ratan(A)
double precision ratan,A,Pi
Pi = 3.14159
nanflg = isnan(A)
if(A.gt.1.d10) then
ratan = Pi/2
else if(A.lt.-1.d10) then
ratan = 3/2*Pi
else if(nanflg.eq..true.)then
ratan = 0
x1 = -1d6; y1 = -1d6
else
ratan = atan(A)
end if
end
c function rtan(A)
c double precision rtan,A
c nanflg = isnan(A)
c if(nanflg.eq..true.) then
c rtan = 0
c x1 = -1d6; y1 = -1d6
c else
c rtan = tan(A)
c end if
c end
c function rsin(A)
c double precision rsin,A
c nanflg = isnan(A)
c if(nanflg.eq..true.) then
c rsin = 0
c x1 = -1d6; y1 = -1d6
c else
c rsin = sin(A)
c end if
c end
c
c function rasin(A)
c double precision rasin,A
c nanflg = isnan(A)
c if(nanflg.eq..true.) then
c rasin = 0
c x1 = -1d6; y1 = -1d6
c else
c rasin = tan(A)
c end if
c end
c
c function rcos(A)
c double precision rcos,A
c nanflg = isnan(A)
c if(nanflg.eq..true.) then
c rcos = 0
c x1 = -1d6; y1 = -1d6
c else
c rcos = cos(A)
c end if
c end
c
c function racos(A)
c double precision racos,A
c nanflg = isnan(A)
c if(nanflg.eq..true.) then
c racos = 0
c x1 = -1d6; y1 = -1d6
c else
c racos = tan(A)
c end if
c end
c
c
c