program main_EUp 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 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 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 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 integer ranno integer seedc double precision Nempty1, Nempty2 C parameter(ranno=1) double precision r2d, lambda, sig, intlen, thscat double precision phi, ztar1 double precision pltot,thtot,pl_del,rotang double precision lsrot,lsrotsq,lsrotavg,lsrottsc,lsrottscsq double precision lsrotgsc,lsrotgscsq,lsde,lsdesq double precision rsrot,rsrotsq,rsrotavg,rsrottsc,rsrottscsq double precision rsrotnsc,rsrotnscsq,lsrotnsc,lsrotnscsq double precision rsrotgsc,rsrotgscsq,rsde,rsdesq double precision thxinput, thyinput double precision nr,nrts,nrgs,nl,nlts,nlgs,nder,ndel, nb double precision nrsq,nrtssq,nrgssq,nlsq,nltssq,nlgssq double precision ndersq,ndelsq,nrtsraw,nltsraw,nrgsraw,nlgsraw 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 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 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 c---- double precision lwde,rwde,lwdesq,rwdesq 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----------------------------------------------------------------------- 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 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,-3.5,3.5,0.) call hbook1(3002,'y (cm) Det',100,-6.,6.,0.) call hbook1(3003,'thx (rad) Det',200,-0.1,0.1,0.) call hbook1(3004,'thy (rad) Det',200,-0.1,0.1,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) scat DetE',2000,-3.,3.,0.) call hbook1(3022,'E_loss (meV) scat DetW',2000,-3.,3.,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.) 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.) 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,212.0, & 215.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,'Path length (cm) Det scat West',300,492.20, & 492.60,0.) call hbook1(6007,'Path length (cm) Det scat East',300,492.20, & 492.60,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,'Path length (cm) Det no scat East',300,492.20, & 492.60,0.) call hbook1(6013,'Path length (cm) Det no scat West',300,492.20, & 492.60,0.) call hbook1(6015,'Crossover Path length (cm) Det scat West', & 300,480.50,480.90,0.) call hbook1(6017,'Crossover Path length (cm) Det scat East', & 300,480.50,480.90,0.) call hbook1(6035,'Crossover Path length (cm) DetEast', & 300,480.50,480.90,0.) call hbook1(6037,'Crossover Path length (cm) DetWest', & 300,480.50,480.90,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 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 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 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. c do 101 i = 1, nevts Event = 1 pltot = 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. c c Generate initial neutron wavelength 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 cbec2010temp c if (i.eq.1) then c write(*,*) ' All Lambda = 5Ang ' c write(3,*) ' All Lambda = 5Ang ' c endif c xlam = 5.0d0 cebec2010temp call hf1(4004,real(xlam),1.) !Initial lambda histogram -KG 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 ccccccccccccccccccccccc No x-y gradient c if (i.eq.1) then c write(*,*) ' ! Running with No Position (x) gradient !' c write(3,*) ' ! Running with No Position (x) gradient !' c end if c xtest = xin*(2.*vect(3)-1.) !-KG c x0 = xtest c xstart = x0 !-KG c y0 = yin*(2.*vect(4)-1.) c ccccccccc gradient cc 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 cc this gives double gradient -- gradient across each half... cc if(xtest.gt.0.) then cc x0 = HRNDM1(151) ! gives gradient only acrros half ccc x0 = 3.d0 - HRNDM1(151) cc elseif(xtest.lt.0) then cc x0 = HRNDM1(151) - 3. ccc x0 = -HRNDM1(151) cc endif cc cc xstart = x0 ccccccccc cccccccccccc x-y from PSM image c 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 ccccccccc Inverted Parabola x - y distributions c c if (i.eq.1) then c write(6,*) ' !!! Inverted Parabola x-y dist!!' c write(3,*) ' !!! Inverted Parabola x-y dist!!' c end if c c xtest = HRNDM1(161) c y0 = HRNDM1(162) c x0 = xtest c xstart = xtest c if (x0.gt.xmax) xmax = x0 c if (y0.gt.ymax) ymax = y0 c c if ((abs(x0).gt.xin).or.(abs(y0).gt.yin)) goto 100 c c if (x0**2+y0**2.le.0.625**2) c & INPKflux = INPKflux + 1/pi/0.625**2/nevts c c ccccccccc 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.. c isotropic 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)) c c flat c thx0 = ng6indexx*xlam*(vect(1)-0.5d0)*2.0d0 c thy0 = ng6indexy*xlam*(vect(2)-0.5d0)*2.0d0 c cosine dist c !!!!set maximum divergence in const.inc ng6indexx,y c thx0 = 2*ng6indexx*xlam/pi*dasin(2.0d0*vect(1)-1.d0) c thy0 = 2*ng6indexy*xlam/pi*dasin(2.0d0*vect(2)-1.d0) C C********for nSPinSim3.1.4 choosing thx0 and thy0 from hbooks, interpolated C******distributions from J. Cook simulations of NG6 and NGC c if(xlam.lt.2.or.xlam.gt.12.) goto 101 c nlamx = int(xlam*10+1)+980 c nlamy = int(xlam*10+1)+1080 c call hrget(nlamx,'hbooks/ngcx.hbook','') c thx0 = HRNDM1(nlamx) c call hdelet(nlamx) c call hrget(nlamy,'hbooks/ngcy.hbook','') c thy0 = HRNDM1(nlamy) c call hdelet(nlamy) C************************************************************** CC**************************************************************** c thxinput = thx0 thyinput = thy0 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) ******************************************************************** * 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 (airflg1) Nairscat1 = Nairscat1+1 ! count scatters with energy change bec if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot ! mrad if (pltot.gt.500) 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 ******************************************************************** * 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 (alflg1) Nalscat1 = Nalscat1+1 ! count scatters with energy change bec if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if (en1.lt.0.d0) then write(6,*) 'Al1 en<1, en=1e-6' en1 = 0.0000001 end if pl_del = 0.d0 pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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 ******************************************************************** 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 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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot pltot = pl_del + pltot thtot = thtot + rotang if (pltot.gt.500) 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. ******************************************************************** ******************************************************************** * 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 (alflg2) Nalscat2 = Nalscat2+1 !rcm if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if(en1.lt.0.d0) then write(6,*) 'All en<1, en=1e-6' en1 = 0.0000001 end if pl_del = 0.d0 pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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 ************************************************************************* * 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 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.500) 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 ************************************************************************ * 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 INCOAVGflux = INCOAVGflux + 1/(ix2-ix1)/(iy2-iy1)/nevts if ((2.0+x1)**2+y1**2.le.0.625**2) & INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts 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 c This path length is now OK for reflected neutrons -KG pl_del = iz0/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2))) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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) ******************************************************************** * 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 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.500) 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 ******************************************************************** ******************************************************************** * 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) c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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 * ******************************************************************** ******************* INput Coil Exit Flux Measurement*************** INCOAVGflux = INCOAVGflux + 1/(ix2-ix1)/(iy2-iy1)/nevts if (x1**2+y1**2.le.0.625**2) & INCOPKflux = INCOPKflux + 1/pi/0.625**2/nevts 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) ******************************************************************** ******************************************************************** * 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) c call collimator(x0,x1,y0,y1) !Boroflex septum to prevent crossovers if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 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*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.500) 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 * ******************************************************************** * 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 (en1.lt.0.d0) then write(6,*) 'Al2 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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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 ******************************************************************** 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 c Given current energy determine q from S(q) and delE from omega(q) k0 = rsqrt(2*m_n*en0)/h_par c c enflg1 = .false. thflg1 = .false. empflg1 = .false. if( (x0.ge.tarx1.and.x0.le.tarx2) ) then ! East=neg. x, Target UPstream sig = 0.0 ; thscat = 0.0 thp = 0.0; php = 0.0 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, tarx1, tarx2, & en1, x1, y1, z1, thx1, thy1, sig, thscat, & scatflg1,nrfl1,ztar1,pl_del,rotang) ! bec if(en1.ne.en0) enflg1=.true. If ((.not.hegastar).and.(.not.airtar)) wgt2 = wgtA/sig if ((en1.gt.-1d5).and.(abs(wgt2).gt.1)) & write(*,*) ' wgt2 ',wgt2, en0, en1, sig ! bec 2010 c make sure really went in target -- we know x-position is ok, check y-position if (y0.ge.tary1.and.y0.le.tary2) then Ntar1=Ntar1+wgt ! entering target, before wgt change if((en0.gt.-100).and.(en1.gt.-100).and.(en0.ne.en1)) then ! scattered 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 call hf1(7001,real(en0-en1),wgt) endif 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) end if if ((x1.ge.tarx1.and.x1.le.tarx2).and. & (y1.ge.tary1.and.y1.le.tary2)) & Ntar1out = Ntar1out+wgt ! out of target west side only elseif( (x0.ge.tarx3.and.x0.le.tarx4) ) then ! West=pos. x, Empty UP if (y0.ge.tary1.and.y0.le.tary2) Nemp1in=Nemp1in+wgt call empty(en0, x0, y0, z0, thx0, thy0, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1,empflg1, ztar1,pl_del,rotang) if (en1.lt.0.d0) goto 100 c write(6,*) 'tar1 en<1, en=1e-6' c en1 = 0.0000001 c end if if (empflg1) Nempty1 = Nempty1+1.0 ! not trajectory dependent -- no wgt if ((x1.ge.tarx3.and.x1.le.tarx4).and. & (y1.ge.tary1.and.y1.le.tary2)) & Nemp1out = Nemp1out+wgt ! out of empty 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 thtot = thtot + rotang if(en0.eq.en1) then call hf1(1005,real(en1),wgt) else enflg1=.true. 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 * ******************************************************************** * 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 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.500) 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 ******************************************************************** * 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 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.500) write(6,*) 'Al Pi ',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 *********************************************************************** 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 call airgap(en0, x0, y0, z0, thx0, thy0, & tairx1,tairx2,tairy1,tairy2,tairz3, en1, x1, y1, z1, & thx1,thy1,airflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 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 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) ********************************************************************* * 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 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.500) write(6,*) 'Al Pi ',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 c*************PI Coil reversal of rotation angle c using 1/2 gap length -- rotate to halfway point, reverse, rotate the other half pl_del = 0.5*rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot call rotation(en0,x0,y0,z0,0.5*(x1+x0), & 0.5*(y1+y0),0.5*(z1+z0),ztar1,rotang) 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 c reverse thtot in center of gap and rotate the other half of gap length if (picoilflg) then ! pi-coil ON thtot = - thtot pltot = pl_del + pltot else ! pi-coil OFF pltot = pl_del + pltot end if call rotation(en0,0.5*(x1+x0),0.5*(y1+y0),0.5*(z1+z0), & x1,y1,z1,ztar1,rotang) thtot = thtot + rotang if (pltot.gt.500) write(6,*) 'gap 2 ',pltot if (dabs(thtot).gt.1.0e9) then write(6,*) ' targap thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) ' z ', z1 end if * ******************************************************************** * 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 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.500) 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 ******************************************************************** 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 write(*,*) 'tar 2 starts at z = ',z0 c c call hf1(203,real(x0),1.) enflg2 = .false. thflg2 = .false. empflg2 = .false. if( (x0.gt.tarx1.and.x0.lt.tarx2) ) then ! East = neg. x, Empty down call empty(en1, x0, y0, z0, thx0, thy0, tarx1, tarx2, & en1, x1, y1, z1, thx1, thy1, empflg2, ztar1, pl_del, rotang) if (empflg2) Nempty2 = Nempty2+1.0 !not trajectory dependent -- no wgt if (y0.ge.tary1.and.y0.le.tary2) Nemp2in=Nemp2in+wgt if ((x1.ge.tarx1.and.x1.le.tarx2).and. & (y1.ge.tary1.and.y1.le.tary2)) & Nemp2out = Nemp2out+wgt ! out of empty east side only elseif( (x0.gt.tarx3.and.x0.lt.tarx4) ) then ! West = pos. x, empty DOWN if (enflg1) 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) hws2 = q*1e3+700 if(hws2.eq.700) goto 502 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 k0 = rsqrt(2*m_n*en0)/h_par 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, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1, sig, thscat, & scatflg2, nrfl2, ztar1, pl_del, rotang) ! bec if (en1.ne.en0) enflg2 = .true. If ((.not.hegastar).and.(.not.airtar)) wgt2 = wgtA/sig c make sure really went in target -- we know x-position is ok, check y-position if (y0.ge.tary1.and.y0.le.tary2) then Ntar2=Ntar2+wgt ! entering target, before wgt change c if(en0.gt.-100.and.en1.gt.-100.and.en0.ne.en1) then ! scattered 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 call hf1(7002,real(en0-en1),wgt) endif 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 end if if ((x1.ge.tarx3.and.x1.le.tarx4).and. & (y1.ge.tary1.and.y1.le.tary2)) Ntar2out = Ntar2out+wgt ! out of target west side only else x1 = -1.d6 y1 = -1.d6 z1 = -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 thtot = thtot + rotang if(en1.eq.en0) then call hf1(2005,real(en1),wgt) else enflg2=.true. call hf1(2012,real(en1),wgt) endif cccccccccccc end target 2 rotaton if (pltot.gt.500) 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 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 if (en1.ne.en0) then call hf1(5024,real(180./3.14159*thx1),wgt) call hf1(5025,real(180./3.14159*thy1),wgt) end if call hf1(5026,real(180./3.14159*thx1),wgt) call hf1(5027,real(180./3.14159*thy1),wgt) ******************************************************************** ******************************************************************** * 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 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.500) 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 *********************************************************************** * 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(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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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 ******************************************************************** ******************************************************************** * 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 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.500) 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 ******************************************************************** * 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 c 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 c This path length is not OK for reflected neutrons -KG c pl_del = rsqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pl_del = oz/dcos(ratan(rsqrt((dtan(thx0))**2+(dtan(thy0))**2))) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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) * ******************************************************************** * 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 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.500) 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 ************************************************************************ ************************************************************************ * 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 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**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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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) ******************************************************************** ******************************************************************** * 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*pl_del*0.1/rsqrt(2.*en1/m_n) + thtot if (pltot.gt.500) 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 ******************************************************************** * 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 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.500) 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 ******************************************************************* * 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 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.600) write(6,*) 'gap 6 ',pltot if (pltot.gt.600) 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) ******************************************************************** ******************************************************************* 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*480.6*0.1/rsqrt(2.*en1/m_n) , c & gy*mag_B*pltot*0.1/rsqrt(2.*en1/m_n) ! mrad c & , pltot,thtot c if((xstart*x1).lt.0.) Nxo=Nxo+wgt !count cross-over neutrons -KG 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 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(pltot),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.((enflg1.or.enflg2).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 (en1.eq.enstart) then 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 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(pltot),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.((enflg1.or.enflg2).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 (en1.eq.enstart) then !rcm PAsin 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 endif endif call hf1(7005, real(wgt),1.0) if (wgt.lt.0) write(6,*) ' end loop wgt<0 ',wgt, sig, en1 if (entar1.ne.enstart) numendif1 = numendif1 + wgt if (enflg1) numtar1flg = numtar1flg + wgt if (entar2.ne.entar1) numendif2 = numendif2 + wgt if (enflg2) numtar2flg = numtar2flg + wgt if (en1.ne.enstart) numendifdet = numendifdet +wgt ****************Rotations and Histograms for Scattered (and NonScat) neutrons**************** if(en1.ne.enstart) then ! scattered and changed energy anywhere 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(pltot),wgt) call hf1(6008, real(thtot),wgt) if ((xtest.gt.0).and.(x1.lt.0)) then call hf1(6017, real(pltot),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.((enflg1.or.enflg2).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 lwde = lwde + (enstart-en1) lwdesq = lwdesq + (enstart-en1)**2 ndelw = ndelw + 1 end if if (tar2flg) then ! East Dn 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 ctemp if (wgt.eq.1.0) write(6,*) ' tar2 wgt=1 ' else 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(pltot),wgt) call hf1(6006, real(thtot),wgt) if ((xtest.lt.0).and.(x1.gt.0)) then call hf1(6015, real(pltot),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.((enflg1.or.enflg2).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 rwde = rwde + (enstart-en1) rwdesq = rwdesq + (enstart-en1)**2 nderw = nderw + 1 end if if (tar1flg) then ! West Up 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 ctemp if (wgt.eq.1.0) write(6,*) ' tar1 wgt=1 ' else 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 else ! non-scattered neutrons c call hf1(6009, real(pltot),wgt) c call hf1(6010, real(thtot),wgt) if (x1.lt.0.) then call hf1(6011, real(pltot),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) else call hf1(6013, real(pltot),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) 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. & ((x1*xtar1.lt.0.0).or.(x1*xtar2.lt.0.0)) ) then c write(6,*) xtar1, xtar2, x1 dncross = dncross + 1 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) 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 enflg1 ',numtar1flg write(6,*) ' num of endif 1',numendif1 write(3,*) ' num of enflg1 ',numtar1flg write(3,*) ' num of endif 1',numendif1 write(6,*) ' num of enflg2 ',numtar2flg write(6,*) ' num of endif 2',numendif2 write(3,*) ' num of enflg1 ',numtar2flg write(3,*) ' num of endif 2',numendif2 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 East target IN ',Ntar1,Ntar1/NPSMout write(*,*) 'Up East target OUT ',Ntar1out,Ntar1out/NPSMout write(*,*) 'Up West empty IN ',Nemp1in,Nemp1in/NPSMout write(*,*) 'Up West empty Out ',Nemp1out,Nemp1out/NPSMout write(*,*) 'Dn East empty IN ',Nemp2in,Nemp2in/NPSMout write(*,*) 'DN East empty Out ',Nemp2out,Nemp2out/NPSMout write(*,*) 'Dn West target IN ',Ntar2,Ntar2/NPSMout write(*,*) 'Dn West 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 write(*,5653) ix1,ix2,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,1pe10.3,1x,1pe10.3) 5652 format(1x,' PSM Dim (cm)',1x,1pe10.3,1x,1pe10.3) 5653 format(1x,' Incoil Dim (cm)',1x,1pe10.3,1x,1pe10.3,1x,1pe10.3, & 1x,1pe10.3) 5654 format(1x,' Target Dim X (cm)',1x,1pe10.3,1x,1pe10.3,1x,1pe10.3, & 1x,1pe10.3) 5655 format(1x,' Target Dim Y (cm)',1x,1pe10.3,1x,1pe10.3) 5656 format(1x,' Outcol Dim X (cm)',1x,1pe10.3,1x,1pe10.3,1x,1pe10.3, & 1x,1pe10.3) 5657 format(1x,' Outcoil Dim Y (cm)',1x,1pe10.3,1x,1pe10.3) 5658 format(1x,' Detector radius (cm)',1x,1pe10.3) 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 East target IN ',Ntar1,Ntar1/NPSMout write(3,*) 'Up East target OUT ',Ntar1out,Ntar1out/NPSMout write(3,*) 'Up West empty IN ',Nemp1in,Nemp1in/NPSMout write(3,*) 'Up West empty Out ',Nemp1out,Nemp1out/NPSMout write(3,*) 'Dn East empty IN ',Nemp2in,Nemp2in/NPSMout write(3,*) 'DN East empty Out ',Nemp2out,Nemp2out/NPSMout write(3,*) 'Dn West target IN ',Ntar2,Ntar2/NPSMout write(3,*) 'Dn West 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 write(3,5653) ix1,ix2,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) write(6,755) Nenscat1/Ntar1 write(6,756) Nenscat2/Ntar2 write(6,757) Ntar12/Ntar1 write(6,758) cot/float(Nevts) write(6,759) nrts/nr write(6,760) nlts/nl write(6,761) nrgs/nr write(6,762) nlgs/nl write(6,*) write(6,*) ' StDev All West ', & rsqrt((rsrotsq-rsrot**2/nr)/nr) write(6,*) ' St Err All West ', & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr !rcm c write(6,*) ' nrPA ',nrPA , ' nrPAsq ',nrPAsq c write(6,*) 'rsqrt(nrPAsq)/nrPA ', rSqrt(nrPAsq)/nrPA if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then c write(6,*) write(6,*) ' StDev Tar Scat West ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) write(6,*) ' St Err All West ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) & *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) write(6,*) ' St Err All West ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) & *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, 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*rsrotavg/rsnscPAavg**2)**2)/ & rsqrt(1-(rsrotavg/rPAavg)**2)*1d3 write(6,725) real(lsnscPAavg*100), & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt((unclsnsc/lsnscPAavg)**2+ !rcm & (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/ & rsqrt(1-(lsrotavg/lPAavg)**2)*1d3 write(6,732)dasin(rsrotnscavg/rsnscPAavg)*1d3- & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt( ((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotavg/rsnscPAavg**2)**2)/ & (1-(rsrotavg/rsnscPAavg)**2) + ((unclsnsc/lsnscPAavg)**2+ & (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/ & (1-(lsrotavg/lPAavg)**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,*) 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) if ((nrsc.gt.1).and.(nlsc.gt.1)) then !rcm 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,*) argrtsc = (rsrottscsq-rsrottsc**2/nrts)/nrts argrtscPA = (nrtsPAsq-nrtsPA**2/nrts)/nrts c if (argrtsc.lt.0) then c write(6,*) 'argrtsc',argrtsc c argrtsc =0 c endif c if (argrtscPA.lt.0) then c write(6,*) 'argrtscPA',argrtscPA c argrtscPA =0 c endif uncrstsc = rsqrt(argrtsc)*rsqrt(nrtssq)/nrts uncrstscPA = rsqrt(argrtscPA)*rsqrt(nrtssq)/nrts argltsc = (lsrottscsq-lsrottsc**2/nlts)/nlts argltscPA = (nltsPAsq-nltsPA**2/nlts)/nlts c if (argltsc.lt.0) then c write(6,*) 'argltsc',argltsc c argltsc =0 c endif c if (argltscPA.lt.0) then c write(6,*) 'argltscPA',argltscPA c argltscPA =0 c endif unclstsc = rsqrt(argltsc)*rsqrt(nltssq)/nlts unclstscPA = rsqrt(argltscPA)*rsqrt(nltssq)/nlts if (nlts.gt.0) then 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,*) argrgsc = (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs argrgscPA = (nrgsPAsq-nrgsPA**2/nrgs)/nrgs c if (argrgsc.lt.0) then c write(6,*) 'argrgsc',argrgsc c argrgsc =0 c endif c if (argrgscPA.lt.0) then c write(6,*) 'argrgscPA',argrgscPA c argrgscPA =0 c endif uncrsgsc = rsqrt(argrgsc)* rsqrt(nrgssq)/nrgs uncrsgscPA = rsqrt(argrgscPA)*rsqrt(nrgssq)/nrgs arglgsc = (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs arglgscPA = (nlgsPAsq-nlgsPA**2/nlgs)/nlgs c if (arglgsc.lt.0) then c write(6,*) 'arglgsc',arglgsc c arglgsc =0 c endif c if (arglgscPA.lt.0) then c write(6,*) 'arglgscPA',arglgscPA c arglgscPA =0 c endif unclsgsc = rsqrt(arglgsc)* rsqrt(nlgssq)/nlgs unclsgscPA = rsqrt(arglgscPA)*rsqrt(nlgssq)/nlgs if (nlgs.gt.0) then 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 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 ) write(6,*) write(6,*) ' ***************************************************', &'******************' write(6,*) 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 if ((nderw.gt.1).and.(ndelw.gt.1)) then write(6,802) rwde/nderw, & rsqrt((rwdesq-rwde**2/nderw)/nderw/(nderw-1)) write(6,801) lwde/ndelw, & rsqrt((lwdesq-lwde**2/ndelw)/ndelw/(ndelw-1)) end if 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,*) 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,' Rotations: PSMout to ASMout ',f9.4,' ~cm ', & ' B(mG)=',1f9.4,' BzGradNum ',i2) 701 format(3x,' Avg. Energy loss at Det West Dn (meV)',1pe11.4,' +- ' & ,1pe11.4) 702 format(3x,' Avg. Energy loss at Det East Up (meV)',1pe11.4,' +- ' & ,1pe11.4) 801 format(3x,' Avg. E loss at Det WDn (alt. wgt) (meV)',1pe11.4,' +- ' & ,1pe11.4) 802 format(3x,' Avg. E loss at Det EUp (alt. wgt) (meV)',1pe11.4,' +- ' & ,1pe11.4) 703 format(3x,' PA ',f5.1,' West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 704 format(3x,' PA ',f5.1,' East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 803 format(3x,' PA ',f5.1,' Rot West Dn (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 804 format(3x,' PA ',f5.1,' Rot East Up (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 705 format(3x,' PA ',f5.1,' Target Scat West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 706 format(3x,' PA ',f5.1,' Target Scat East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 805 format(3x,' PA ',f5.1' Scat West Dn (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 806 format(3x,' PA ',f5.1,' Scat East Up (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 7005 format(3x,' PA ',f5.1,' Gas Scat West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 7006 format(3x,' PA ',f5.1,' Gas Scat East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 725 format(3x,' PA ',f5.1,' NOScat West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 726 format(3x,' PA ',f5.1,' NOScat East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) c 707 format(3x,' Avg. Rot Scat West Dn E<0.35meV (mrad) ',1pe11.4, c & ' +- ',1pe11.4) c 708 format(3x,' Avg. Rot Scat East Up E<0.35meV (mrad) ',1pe11.4, c & ' +- ',1pe11.4) c 709 format(3x,' Frac of rotation from scat West Dn ',1pe11.4, c & ', oflow En ',1pe11.4) c 710 format(3x,' Frac of rotation from scat East Up ',1pe11.4, c & ', oflow En ',1pe11.4) 711 format(3x,' UpE-DnW (mrad) ',1pe11.4, & ' +- ',1pe11.4) 811 format(3x,' UpE-DnW (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 712 format(3x,' UpE-DnW Tar scat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 812 format(3x,' UpE-DnW Scat (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 7012 format(3x,' UpE-DnW Gas scat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 732 format(3x,' UpE-DnW NOScat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 713 format(3x,' # in DetE Up',1pe13.6,2x,', # in DetW Dn ',1pe13.6) 714 format(3x,' Ratio # DetEUp/DetWDn ,'1pe13.6) 717 format(3x,' Det Counts UpE-DnW ',1pe11.4,' +- ',1pe11.4) 715 format(3x,' # Tarscat in DetE Up ',1pe11.4,2x, & ', # Tarscat in Det W Dn ',1pe11.4) 716 format(3x,' Ratio # Tarscat DetEUp/DetWDn ,'1pe11.4) 718 format(3x,' TarScat Det Counts UpE-DnW ',1pe11.4,' +- ',1pe11.4) 815 format(3x,' AltWgt # Tarscat in DetE Up ',1pe11.4,2x, & ', # Tarscat in Det W Dn ',1pe11.4) 816 format(3x,' AltWgt Ratio # Tarscat DetEUp/DetWDn ,'1pe11.4) 818 format(3x,' AltWgt TarScat Det Counts UpE-DnW ',1pe11.4,' +- ' & ,1pe11.4) 7015 format(2x,'# Gas/Al/Xdn scat DetE Up ',1pe11.4,2x, & ',# Gas/Al/Xdn scat DetW Dn ',1pe11.4) 7016 format(3x,' Ratio # Gas/Al/tarcross Scat DetEUp/DetWDn ,'1pe11.4) 7018 format(3x,' G/Al/TX Scat Det Counts UpE-DnW ',1pe11.4,' +- ', & 1pe11.4) 719 format(3x,'# scat w/ E<0.35meV at det East: ',1pe11.4, & ' West: ',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 of det hits that had TarScat E ',1pe11.4) 760 format(1x,' Frac of det hits that had TarScat W ',1pe11.4) 761 format(1x,' Frac of det hits that had G/Al/TXScat E ',1pe11.4) 762 format(1x,' Frac of det hits that had G/Al/TXScat W',1pe11.4) 770 format(3x,' num NO scat -- det E: ',1pe11.4, & ' det W: ',1pe11.4) 780 format(3x,'Total No. of cross-over neutrons: ',1pe11.4) 782 format(3x,'Unwt # Tar scat DetE Up ',1pe11.4,2x, & ', Unwt # Tar scat Det W Dn ',1pe11.4) 7082 format(3x,'UnWt # Gas scat DetE Up ',1pe11.4,2x, & ', Unwt # Gas scat Det W Dn ',1pe11.4) 735 format(3x,' Guide output E ',1pe11.4,' Guide output W ',1pe11.4) 736 format(3x,' Incoil output E ',1pe11.4,' Incoil output W ',1pe11.4) 737 format(3x,' Outcoil output E ',1pe11.4, & ' Outcoil output W ',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) nrts/nr !rcm write(3,760) nlts/nl write(3,761) nrgs/nr write(3,762) nlgs/nl write(3,*) write(3,*) ' StDev All East ', & rsqrt((rsrotsq-rsrot**2/nr)/nr) write(3,*) ' St Err All East ', & rsqrt((rsrotsq-rsrot**2/nr)/nr)*rSqrt(nrsq)/nr c write(3,*) ' nr ',nr , ' nrsq ',nrsq c write(3,*) 'rsqrt(nrsq)/nr ', rSqrt(nrsq)/nr if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then c write(3,*) write(3,*) ' StDev Tar Scat East ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) write(3,*) ' St Err All East ', & rsqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) & *rSqrt(nrtssq)/nrts c write(3,*) ' nrts ',nrts , ' nrtssq ',nrtssq c write(3,*) 'rsqrt(nrtssq)/nrts ', rSqrt(nrtssq)/nrts end if if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then c write(3,*) write(3,*) ' StDev Gas Scat East ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) write(3,*) ' St Err All East ', & rsqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) & *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, bzgradnum write(3,*) write(3,*) ' ********************** Average Rotation Values &***********************' write(3,*) 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 if ((nrall.gt.1).and.(nlall.gt.1)) then !rcm polarization product write(3,*) 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,*) write(3,726) real(rsnscPAavg*100), & dasin(rsrotnscavg/rsnscPAavg)*1d3, & rsqrt((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotavg/rsnscPAavg**2)**2)/ & rsqrt(1-(rsrotavg/rPAavg)**2)*1d3 write(3,725) real(lsnscPAavg*100), & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt((unclsnsc/lsnscPAavg)**2+ & (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/ & rsqrt(1-(lsrotavg/lPAavg)**2)*1d3 write(3,732)dasin(rsrotnscavg/rsnscPAavg)*1d3- & dasin(lsrotnscavg/lsnscPAavg)*1d3, & rsqrt( ((uncrsnsc/rsnscPAavg)**2+ & (uncrsnscPA*rsrotavg/rsnscPAavg**2)**2)/ & (1-(rsrotavg/rsnscPAavg)**2) + ((unclsnsc/lsnscPAavg)**2+ & (unclsnscPA*lsrotavg/lsnscPAavg**2)**2)/ & (1-(lsrotavg/lPAavg)**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 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) then !rcm 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 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)/nlts 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) then 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,*) write(3,*) ' ***************************************************', &'******************' write(3,*) 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 if ((nderw.gt.1).and.(ndelw.gt.1)) then write(3,802) rwde/nderw, & rsqrt((rwdesq-rwde**2/nderw)/nderw/(nderw-1)) write(3,801) lwde/ndelw, & rsqrt((lwdesq-lwde**2/ndelw)/ndelw/(ndelw-1)) end if 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,*) 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)) end if 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: ', & dncross write(3,*) ' Number crossed over after target area: ', & 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 ', & dble(dncross)/cot write(3,*) ' Frac. of det hits that crossed after targets ', & 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(6,*) ' argsqrt<0 ',B write(3,*) ' 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