program main_EUp C Targe 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, Nairin1 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 double precision Nemp1in,Nemp2in,Nemp1out,Nemp2out, NtarWEin 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 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 pl1_ns, th1_ns, pl2_ns, th2_ns, pl_ns, th_ns double precision pl1a_sc, pl1b_sc, th1a_sc, th1b_sc double precision pl2a_sc, pl2b_sc, th2a_sc, th2b_sc double precision pl1_sc, th1_sc, pl2_sc, th2_sc, pl3, th3, phi double precision pltot,thtot,pl_del double precision lsrot,lsrotsq, lsrottsc,lsrottscsq double precision lsrotgsc,lsrotgscsq,lsde,lsdesq double precision rsrot,rsrotsq, 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 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(699,'hbooks/sqw.hbook','') ! New S(q,w) cbec10 call hrget(700,'hbooks/sqw.hbook','') ! New S(q,w) c----------------------------------------------------------------------- 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(2,'x (cm) guide out',100,-5.5,5.5,0.) call hbook1(3,'x (cm) incoil out',100,-5.5,5.5,0.) call hbook1(4,'x (cm) target in',100,-5.5,5.5,0.) call hbook1(5,'x (cm) outcoil out',100,-5.5,5.5,0.) call hbook1(6,'x (cm) ASM out',100,-5.5,5.5,0.) call hbook1(20,'E (meV) start Both',2000,0.,40.,0.) call hbook1(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,480.50, & 480.90,0.) call hbook1(6007,'Path length (cm) Det scat East',300,480.50, & 480.90,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,480.50, & 480.90,0.) call hbook1(6013,'Path length (cm) Det no scat West',300,480.50, & 480.90,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.) call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., & 300,-800.,0.,0.) call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., & 300,-800.,0.,0.) call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., & 300,-800.,0.,0.) call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, & 480.50,480.90,300,-800.,0.,0.) call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., & 300,-800.,0.,0.) call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, & 480.50,480.90,300,-800.,0.,0.) call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', & 800,0.,40.,300,-800.,0.,0.) call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', & 300,480.50,480.90,300,-800.,0.,0.) call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', & 800,0.,40.,300,-800.,0.,0.) call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', & 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.) call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., & 300,-800.,0.,0.) call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., & 300,-800.,0.,0.) call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., & 300,0.,60.,0.) call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, & 480.50,480.90,300,0.,60.,0.) call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., & 300,0.,60.,0.) call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, & 480.50,480.90,300,0.,60.,0.) call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', & 800,0.,40.,300,0.,60.,0.) call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', & 300,480.50,480.90,300,0.,60.,0.) call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', & 800,0.,40.,300,0.,60.,0.) call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', & 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.) call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., & 300,-5.,50.,0.) call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., & 300,-5.,50.,0.) call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., & 300,-10.,10.,0.) call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, & 480.50,480.90,300,-10.,10.,0.) call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., & 300,-10.,10.,0.) call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, & 480.50,480.90,300,-10.,10.,0.) call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', & 800,0.,40.,300,-10.,10.,0.) call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', & 300,480.50,480.90,300,-10.,10.,0.) call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', & 800,0.,40.,300,-10.,10.,0.) call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', & 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.) call hbook2(6030,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., & 300,-5.,50.,0.) call hbook2(6031,'th_rotsc (mrad) vs. E(meV) det East',800,0.,40., & 300,-5.,50.,0.) call hbook2(6020,'th_rot (mrad) vs. E (meV) det West',800,0.,40., & 300,-5.,50.,0.) call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, & 480.50,480.90,300,-5.,50.,0.) call hbook2(6022,'th_rot (mrad) vs. E (meV) det East',800,0.,40., & 300,-5.,50.,0.) call hbook2(6023,'th rot (mrad) vs. Path len (cm) det East',300, & 480.50,480.90,300,-5.,50.,0.) call hbook2(6040,'Cr-ov. th_rot (mrad) vs. E (meV) det West', & 800,0.,40.,300,-5.,50.,0.) call hbook2(6041,'Cr-ov. th rot (mrad) vs. Path len (cm) det West', & 300,480.50,480.90,300,-5.,50.,0.) call hbook2(6042,'Cr-ov. th_rot (mrad) vs. E (meV) det East', & 800,0.,40.,300,-5.,50.,0.) call hbook2(6043,'Cr-ov. th rot (mrad) vs. Path len (cm) det East', & 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.) call hbook2(8011,'Thx vs. E (meV) scat DetB',100,0.,40., & 200,-0.1,0.1,0.) call hbook2(8012,'Thx vs. E-loss (meV) scat DetB',100, & -0.4,0.4,200,-0.1,0.1,0.) c call hbook2(9999,'PA vs.WaveLength (Ang)',160,0.,16.,102, c & 0.,1.02,0.) !-KG call hbook2(9998,'omega (1/psec) vs. q (1/Ang)',100,0.,0.2,100, & -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 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 nder=0;ndel=0 Nalscat1=0;Nalscat2=0;Nalout1=0;Nalout2=0; Nalin1=0 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 c do 101 i = 1, nevts c 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 c xlam = HRNDM1(99) !New initial spectrum after PSM -KG call hf1(101,real(xlam),1.) !-KG if(xlam.lt.1.) goto 100 cbec2010temp c if (i.eq.1) then c write(*,*) ' All Lambda = 5Ang ' c write(3,*) ' All Lambda = 5Ang ' c endif c xlam = 10.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' en0 = 0.0000001 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 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 xtest = xin*(2.*vect(3)-1.) !-KG if (i.eq.1) then write(*,*) ' ! Running with No Position (x) gradient !' write(3,*) ' ! Running with No Position (x) gradient !' 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 !! ' end if x0 = xtest xstart = x0 !-KG 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 square cross section y0 = yin*(2.*vect(4)-1.) z0= 0.0d0 if (i.eq.1) then 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 write(6,*) write(3,*) end if 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 = datan(dtan(th0)*dcos(ph0)) c thy0 = datan(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) be thx0 = datan(dtan(th0)*dcos(ph0)) thy0 = datan(dtan(th0)*dsin(ph0)) 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(1)-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**************************************************************** thxinput = thx0 thyinput = thy0 c NPSMout = NPSMout + 1 ! bec count exiting PSM call hf1(1,real(x0),1.) ******************************************************************** * Neutrons passing thruogh 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 = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot ! mrad if (pltot.gt.500) write(6,*) '0th air gap',pltot Nairout1 = Nairout1 + 1 if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'guide ',pltot Nalout1 = Nalout1 + 1 if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if ******************************************************************** c Neutrons passing thruogh 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 call guide(en0,x0,y0,z0,thx0,thy0,x1,y1,z1, & thx1,thy1,ngbx,ngby) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 call hf1(2,real(x1),1.) if (ngbx.eq.0) ngbx0=ngbx0+1 if (ngbx.eq.1) ngbx1=ngbx1+1 if (ngbx.eq.2) ngbx2=ngbx2+1 if (ngbx.eq.3) ngbx3=ngbx3+1 if (ngby.eq.0) ngby0=ngby0+1 if (ngby.eq.1) ngby1=ngby1+1 if (ngby.eq.2) ngby2=ngby2+1 if (ngby.eq.3) ngby3=ngby3+1 c This path length is now OK for reflected neutrons -KG pl_del = gl/dcos(datan(sqrt((dtan(thx0))**2+(dtan(thy0))**2))) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' guide thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if * if ((thx1.eq.-1.0*thx0).or.(thy1.eq.-1.0*thy0)) bounceflg = .true. ******************************************************************** ******************************************************************** * Neutrons passing thruogh 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 call incoil(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1,nibx,niby) c call collimator(x1,y1) !!1cm center septum -KG if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 Nincoilout = Nincoilout + 1 ! bec counttransmission through guide and Input coil call hf1(3,real(x1),1.) if (nibx.eq.0) nibx0=nibx0+1 if (nibx.eq.1) nibx1=nibx1+1 if (nibx.eq.2) nibx2=nibx2+1 if (nibx.eq.3) nibx3=nibx3+1 if (niby.eq.0) niby0=niby0+1 if (niby.eq.1) niby1=niby1+1 if (niby.eq.2) niby2=niby2+1 if (niby.eq.3) niby3=niby3+1 c This path length is now OK for reflected neutrons -KG pl_del = iz0/dcos(datan(sqrt((dtan(thx0))**2+(dtan(thy0))**2))) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'incoil ',pltot if (thtot.gt.1e9) then write(6,*) ' incoil thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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. ******************************************************************** ******************************************************************** * Neutrons passing thruogh the gap (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,tairz2, en1, x1, y1, z1, & thx1,thy1,airflg1) 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 = datan(dtan(thy1)/dtan(thx1)) elseif(thx1.lt.0..and.thy1.ge.0.) then phi = pi + datan(dtan(thy1)/dtan(thx1)) elseif(thx1.lt.0..and.thy1.lt.0.) then phi = pi + datan(dtan(thy1)/dtan(thx1)) elseif(thx1.ge.0..and.thy1.lt.0.) then phi = 2.*pi + datan(dtan(thy1)/dtan(thx1)) endif c call hf1(5021,real(180./3.14159*phi),1.) pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(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 (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if * ******************************************************************** * Neutrons passing through 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,talz1, 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 = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if ******************************************************************** NtarWEin = NtarWEin + 1 xtar1 = x1 ******************************************************************** c print *, " FIRST half target" * Neutrons passing thruogh 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 c Given current energy determine q from S(q) and delE from omega(q) k0 = sqrt(2*m_n*en0)/h_par 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 = sqrt(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, & pl1a_sc,pl1b_sc,scatflg1,nrfl1) ! 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.tary3.and.y0.le.tary4) Nemp1in=Nemp1in+wgt call empty(en0, x0, y0, z0, thx0, thy0, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1,empflg1,pl1a_sc,pl1b_sc) 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) if(en0.eq.en1) then pl1_ns = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) if (bzgradnum.eq.0) then th1_ns = gy*mag_b*pl1_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.1) then th1_ns = gy*4.866*pl1_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.2) then th1_ns = gy*8.471*pl1_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.3) then th1_ns = gy*(-0.008942)*pl1_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.4) then th1_ns = gy*(-5.722/1.8)*pl1_ns*0.1/sqrt(2.*en1/m_n) end if c write(16,*) th1_ns pl1_sc = 0. th1_sc = 0. pltot = pl1_ns + pltot thtot = th1_ns + thtot call hf1(1005,real(en1),wgt) c call hf1(1006,real(pl1_ns),wgt) c call hf1(1007,real(th1_ns),wgt) if(x0.gt.0..and.x1.gt.0.) then c call hf1(1008,real(pl1_ns),wgt) c call hf1(1009,real(th1_ns),wgt) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(1010,real(pl1_ns),wgt) c call hf1(1011,real(th1_ns),wgt) endif else ************Bz GRADIENT******************************* * if (bzgradnum.eq.0) then th1a_sc = gy*mag_b*pl1a_sc*0.1/sqrt(2.*en0/m_n) th1b_sc = gy*mag_b*pl1b_sc*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.1) then **************************************************************** c Field profile for gradient systematic runs. Elog Sim. post 13 -KG th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl1a_sc*dcos(th0)+23.7,IntBz) if(IntBz.lt.165.01.or.IntBz.gt.2189.39) then write(*,*) 'IntBz=', IntBz,' out of range 165.01~2189.39' endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416-10*pl1a_sc*dcos(th0))=0' goto 101 endif th1a_sc=gy*((IntBz-165.01)/(10.0*pl1a_sc*dcos(th0)))* & pl1a_sc*0.1/sqrt(2.*en0/m_n) th1b_sc=gy*((2189.39-IntBz)/(416.0-10.*pl1a_sc*dcos(th0)))* & pl1b_sc*0.1/sqrt(2.*en1/m_n) c IntBz=165.01(mG*mm)@23.7mm where 1st Tgt begins, 2189.39 @439.7mm at the end else if (bzgradnum.eq.2) then **************************************************************** c Field profile for large uniform field systematic runs. -KG c th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl1a_sc*dcos(th0)+23.7,IntBz) if(IntBz.lt.190.59.or.IntBz.gt.3714.64) then write(*,*) 'IntBz=', IntBz,' out of range 190.59~3714.64' endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416-10*pl1a_sc*dcos(th0))=0' goto 101 endif th1a_sc=gy*((IntBz-190.59)/(10.0*pl1a_sc*dcos(th0)))*pl1a_sc*0.1/ &sqrt(2.*en0/m_n) th1b_sc=gy*((3714.64-IntBz)/(416.0-10.*pl1a_sc*dcos(th0)))*pl1b_sc &*0.1/sqrt(2.*en1/m_n) c IntBz=190.59(mG*mm)@23.7mm where 1st Tgt begins, 3714.64 @439.7mm at the end else if (bzgradnum.eq.3) then *********************************************************************** c Gradient field profile for data run 1154 (biggest recorded gradient) th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl1a_sc*dcos(th0)+23.7,IntBz) cc if(IntBz.lt.165.01.or.IntBz.gt.2189.39) then cc write(*,*) 'IntBz=', IntBz,' out of range 165.01~2189.39' cc endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416-10*pl1a_sc*dcos(th0))=0' goto 101 endif th1a_sc=gy*((IntBz+35.84)/(10.0*pl1a_sc*dcos(th0)))*pl1a_sc*0.1/ &sqrt(2.*en0/m_n) th1b_sc=gy*((-39.54-IntBz)/(416.0-10.*pl1a_sc*dcos(th0)))*pl1b_sc &*0.1/sqrt(2.*en1/m_n) cc IntBz=-35.84(mG*mm)@23.7mm where 1st Tgt begins, -39.54 @439.7mm at the end ************************************************************************** else if (bzgradnum.eq.4) then c Gradient field profile for data run 1449 (flipping gradient) th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl1a_sc*dcos(th0)+23.7,IntBz) if(IntBz.lt.165.01.or.IntBz.gt.2189.39) then write(*,*) 'IntBz=', IntBz,' out of range 165.01~2189.39' endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416-10*pl1a_sc*dcos(th0))=0' goto 101 endif th1a_sc=gy*(((-IntBz+70.5959)/1.8)/(10.0*pl1a_sc*dcos(th0))) &*pl1a_sc*0.1/sqrt(2.*en0/m_n) th1b_sc=gy*(((-2450.759+IntBz)/1.8)/(416.0-10.*pl1a_sc*dcos(th0))) &*pl1b_sc*0.1/sqrt(2.*en1/m_n) cc IntBz=-70.5959(mG*mm)@23.7mm where 1st Tgt begins, -2450.75898 @439.7mm at the end *********************************************************************** end if pl1_ns = 0. th1_ns = 0. pl1_sc = pl1a_sc + pl1b_sc th1_sc = th1a_sc + th1b_sc pltot = pl1_sc + pltot thtot = th1_sc + thtot call hf1(1012,real(en1),wgt) c call hf1(1013,real(pl1_sc),wgt) c call hf1(1014,real(th1_sc),wgt) if(x0.gt.0..and.x1.gt.0.) then c call hf1(1015,real(pl1_sc),wgt) c call hf1(1016,real(th1_sc),wgt) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(1017,real(pl1_sc),wgt) c call hf1(1018,real(th1_sc),wgt) endif endif if (pltot.gt.500) write(6,*) 'target 1 ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl1_ns ',pl1_ns write(6,*) ' pl1_sc ',pl1_sc,', en ',en1 write(6,*) ' en0 ',en0, pltot end if c write(*,*) 'tar 1 ends at z = ',z1 !-KG * ******************************************************************** * Neutrons passing through 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,talz1, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if ******************************************************************** ctemp if (wgt.ne.1.0) write(17,*) ' tar1 out wgt ',wgt,tar1flg,tar2flg ******************************************************************** 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 = datan(dtan(thy1)/dtan(thx1)) elseif(thx1.lt.0..and.thy1.ge.0.) then phi = pi + datan(dtan(thy1)/dtan(thx1)) elseif(thx1.lt.0..and.thy1.lt.0.) then phi = pi + datan(dtan(thy1)/dtan(thx1)) elseif(thx1.ge.0..and.thy1.lt.0.) then phi = 2.*pi + datan(dtan(thy1)/dtan(thx1)) endif c call hf1(5022,real(180./3.14159*phi),wgt) 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*sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot if (bzgradnum.eq.0) then thtot = gy*mag_b*pl_del*0.1/sqrt(2.*en1/m_n) + thtot else if (bzgradnum.eq.1) then thtot = gy*0.371*pl_del*0.1/sqrt(2.*en1/m_n) + thtot else if (bzgradnum.eq.2) then thtot = gy*10.165*pl_del*0.1/sqrt(2.*en1/m_n) + thtot end if if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl1_ns ',pl1_ns write(6,*) ' pl1_sc ',pl1_sc,', en ',en1 write(6,*) ' en0 ',en0, pltot write(6,*) 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 if (bzgradnum.eq.0) then thtot = gy*mag_b*pl_del*0.1/sqrt(2.*en1/m_n) + thtot else if (bzgradnum.eq.1) then thtot = gy*(-0.803)*pl_del*0.1/sqrt(2.*en1/m_n) + thtot else if (bzgradnum.eq.2) then thtot = gy*10.167*pl_del*0.1/sqrt(2.*en1/m_n) + thtot end if if (pltot.gt.500) write(6,*) 'gap 2 ',pltot if (thtot.gt.1e9) then write(6,*) ' targap thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot write(6,*) end if * ******************************************************************** * Neutrons passing through 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,talz1, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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,pl2a_sc,pl2b_sc) 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, target 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 = sqrt(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 = sqrt(2*m_n*en0)/h_par kp2 = sqrt(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, & pl2a_sc, pl2b_sc,scatflg2, nrfl2) ! 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) c call hf1(2003,real(thx1),wgt) c call hf1(2004,real(thy1),wgt) if(en1.eq.en0) then pl2_ns = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) if (bzgradnum.eq.0) then th2_ns = gy*mag_b*pl2_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.1) then th2_ns = gy*(-4.864)*pl2_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.2) then th2_ns = gy*8.326*pl2_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.3) then th2_ns = gy*0.07834*pl2_ns*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.4) then th2_ns = gy*(5.671/1.8)*pl2_ns*0.1/sqrt(2.*en1/m_n) end if pl2_sc = 0. th2_sc = 0. pltot = pl2_ns + pltot thtot = th2_ns + thtot call hf1(2005,real(en1),wgt) c call hf1(2006,real(pl2_ns),wgt) c call hf1(2007,real(th2_ns),wgt) if(x0.gt.0..and.x1.gt.0.) then c call hf1(2008,real(pl2_ns),wgt) c call hf1(2009,real(th2_ns),wgt) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(2010,real(pl2_ns),wgt) c call hf1(2011,real(th2_ns),wgt) endif else enflg2=.true. if (bzgradnum.eq.0) then th2a_sc = gy*mag_b*pl2a_sc*0.1/sqrt(2.*en0/m_n) th2b_sc = gy*mag_b*pl2b_sc*0.1/sqrt(2.*en1/m_n) else if (bzgradnum.eq.1) then ******************************************************************** c Field profile for gradient systematic runs. Elog Sim. post 13 -KG th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl2a_sc*dcos(th0)+521.7,IntBz) if(IntBz.lt.145.89.or.IntBz.gt.2169.38) then write(*,*) 'IntBz=',IntBz,' out of range 145.89~2169.38' endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416.0-10*pl2a_sc*dcos(th0))=0' goto 101 endif th2a_sc =gy*((IntBz-2169.38)/(10.*pl2a_sc*dcos(th0)))* & pl2a_sc*0.1/sqrt(2.*en0/m_n) th2b_sc = gy*((145.89-IntBz)/(416.-10.*pl2a_sc*dcos(th0)))* & pl2b_sc*0.1/sqrt(2.*en1/m_n) c IntBz=2169.38(mG*mm)@521.7mm where 2nd Tgt begins, 145.89 @937.7mm at the end else if (bzgradnum.eq.2) then ******************************************************************** c Field profile for large uniform field systematic runs. -KG th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl2a_sc*dcos(th0)+521.7,IntBz) if(IntBz.lt.4548.25.or.IntBz.gt.8011.66) then write(*,*) 'IntBz=',IntBz,' out of range 4548.25~8011.66' endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416.0-10*pl2a_sc*dcos(th0))=0' goto 101 endif th2a_sc =gy*((IntBz-4548.25)/(10.*pl2a_sc*dcos(th0)))*pl2a_sc*0.1 &/sqrt(2.*en0/m_n) th2b_sc =gy*((8011.66-IntBz)/(416.-10.*pl2a_sc*dcos(th0)))*pl2b_sc &*0.1/sqrt(2.*en1/m_n) c IntBz=4548.25(mG*mm)@521.7mm where 2nd Tgt begins, 8011.66 @937.7mm at the end else if (bzgradnum.eq.3) then *********************************************************************** c Gradient field profile for data run 1154 (biggest recorded gradient) -KG th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl2a_sc*dcos(th0)+521.7,IntBz) c if(IntBz.lt.145.89.or.IntBz.gt.2169.38) then c write(*,*) 'IntBz=',IntBz,' out of range 145.89~2169.38' c endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416.0-10*pl2a_sc*dcos(th0))=0' goto 101 endif th2a_sc =gy*((IntBz+40.95)/(10.*pl2a_sc*dcos(th0)))*pl2a_sc*0.1 &/sqrt(2.*en0/m_n) th2b_sc = gy*((-8.36-IntBz)/(416.-10.*pl2a_sc*dcos(th0)))*pl2b_sc &*0.1/sqrt(2.*en1/m_n) cc IntBz=-40.95(mG*mm)@521.7mm where 2nd Tgt begins, -8.36 @937.7mm at the end else if (bzgradnum.eq.4) then ******************************************************************************* c Gradient field profile for data run 1449 (flipping gradient) -KG th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) call int_Bz(pl2a_sc*dcos(th0)+521.7,IntBz) if(IntBz.lt.145.89.or.IntBz.gt.2169.38) then write(*,*) 'IntBz=',IntBz,' out of range 145.89~2169.38' endif if(dcos(th0).eq.0.) then write(*,*) 'cos(th0)=0' goto 101 endif if((416.-10.*pl1a_sc*dcos(th0)).eq.0.)then write(*,*)'(416.0-10*pl2a_sc*dcos(th0))=0' goto 101 endif th2a_sc =gy*(((-IntBz+2418.265)/1.8)/(10.*pl2a_sc*dcos(th0))) &*pl2a_sc*0.1/sqrt(2.*en0/m_n) th2b_sc = gy*(((-59.181+IntBz)/1.8)/(416.-10.*pl2a_sc*dcos(th0))) &*pl2b_sc*0.1/sqrt(2.*en1/m_n) cc IntBz=-2418.265(mG*mm)@521.7mm where 2nd Tgt begins, -59.181 @937.7mm at the end ******************************************************************************* end if pl2_ns = 0. th2_ns = 0. pl2_sc = pl2a_sc + pl2b_sc th2_sc = th2a_sc + th2b_sc pltot = pl2_sc + pltot thtot = th2_sc + thtot call hf1(2012,real(en1),wgt) c call hf1(2013,real(pl2_sc),wgt) c call hf1(2014,real(th2_sc),wgt) if(x0.gt.0..and.x1.gt.0.) then c call hf1(2015,real(pl2_sc),wgt) c call hf1(2016,real(th2_sc),wgt) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(2017,real(pl2_sc),wgt) c call hf1(2018,real(th2_sc),wgt) endif endif if (pltot.gt.500) write(6,*) 'target 2 ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl2_ns ',pl2_ns write(6,*) ' pl2_sc ',pl2_sc,', en ',en1 write(6,*) ' en0 ',en0, pltot write(6,*) thx0,thy0, thx1, thy1 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 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 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,talz1, en1, x1, y1, z1, & thx1,thy1,alflg1) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 if(thx1.ge.0..and.thy1.ge.0.) then phi = datan(dtan(thy1)/dtan(thx1)) elseif(thx1.lt.0..and.thy1.ge.0.) then phi = pi + datan(dtan(thy1)/dtan(thx1)) elseif(thx1.lt.0..and.thy1.lt.0.) then phi = pi + datan(dtan(thy1)/dtan(thx1)) elseif(thx1.ge.0..and.thy1.lt.0.) then phi = 2.*pi + datan(dtan(thy1)/dtan(thx1)) endif c call hf1(5023,real(180./3.14159*phi),1.) pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'gap 3 ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if ******************************************************************** ******************************************************************** * Neutrons passing thruogh 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 outcoil(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1,nobx,noby) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 call hf1(5,real(x1),wgt) if (nobx.eq.0) nobx0=nobx0+1 if (nobx.eq.1) nobx1=nobx1+1 if (nobx.eq.2) nobx2=nobx2+1 if (nobx.eq.3) nobx3=nobx3+1 if (noby.eq.0) noby0=noby0+1 if (noby.eq.1) noby1=noby1+1 if (noby.eq.2) noby2=noby2+1 if (noby.eq.3) noby3=noby3+1 c This path length is not OK for reflected neutrons -KG c pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pl_del = oz/dcos(datan(sqrt((dtan(thx0))**2+(dtan(thy0))**2))) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'out coil ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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. * ******************************************************************** ******************************************************************** * 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) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'gap 5 ',pltot NOutcoilout = NOutcoilout + wgt if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if ******************************************************************** ******************************************************************** * Neutrons passing thruogh the supper mirror." c print *," Supper 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 = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) pltot = pl_del + pltot thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot if (pltot.gt.500) write(6,*) 'ASM ',pltot if (thtot.gt.1e9) then write(6,*) ' ASM thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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 Apply the polarization product -KG lambda=(2.*pi*h_par)/sqrt(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(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 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 if(vect(6).gt.PA) goto 100 call hf2(9999,real(lambda),real(PA),wgt) * ******************************************************************** *************** End of Beamline ********************************** * c write(6,*) gy*mag_B*480.6*0.1/sqrt(2.*en1/m_n) , c & gy*mag_B*pltot*0.1/sqrt(2.*en1/m_n) ! mrad c & , pltot,thtot if((xstart*x1).lt.0.) Nxo=Nxo+wgt !count cross-over neutrons -KG 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 hf1(3003,real(thx1),wgt) call hf1(3004,real(thy1),wgt) pl3 = pl1_ns + pl2_ns + pl1_sc + pl2_sc th3 = th1_ns + th2_ns + th1_sc + th2_sc call hf1(3005,real(en1),wgt) c lambda=(2.*pi*h_par)/sqrt(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 + thtot rallsq = rallsq + thtot**2 nrall = nrall + 1 end if rsrot = rsrot + wgt*thtot rsrotsq = rsrotsq + wgt*thtot**2 nr = nr +wgt nrsq = nrsq + wgt**2 if (en1.eq.enstart) then rsrotnsc = rsrotnsc + wgt*thtot rsrotnscsq = rsrotnscsq + wgt*thtot**2 nrnsc = nrnsc + wgt nrnscsq = nrnscsq + wgt**2 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 + thtot lallsq = lallsq + thtot**2 nlall = nlall + 1 end if lsrot = lsrot + wgt*thtot lsrotsq = lsrotsq + wgt*thtot**2 nl = nl + wgt nlsq = nlsq + wgt**2 if (en1.eq.enstart) then lsrotnsc = lsrotnsc + wgt*thtot lsrotnscsq = lsrotnscsq + wgt*thtot**2 nlnsc = nlnsc + wgt nlnscsq = nlnscsq + wgt**2 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 + thtot lscsq = lscsq + thtot**2 nlsc = nlsc + 1 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*thtot lsrottscsq = lsrottscsq + wgt*thtot**2 nltsraw = nltsraw + 1 nlts = nlts + wgt nltssq = nltssq + wgt**2 ctemp if (wgt.eq.1.0) write(6,*) ' tar2 wgt=1 ' else lsrotgsc = lsrotgsc + wgt*thtot lsrotgscsq = lsrotgscsq + wgt*thtot**2 nlgsraw = nlgsraw + 1 nlgs = nlgs + wgt nlgssq = nlgssq + wgt**2 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) call hf2(8011,real(en1),real(thx1),wgt) 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 + thtot rscsq = rscsq + thtot**2 nrsc = nrsc + 1 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) rsrottsc = rsrottsc + wgt*thtot rsrottscsq = rsrottscsq + wgt*thtot**2 nrtsraw = nrtsraw + 1 nrts = nrts + wgt nrtssq = nrtssq + wgt**2 ctemp if (wgt.eq.1.0) write(6,*) ' tar1 wgt=1 ' else rsrotgsc = rsrotgsc + wgt*thtot rsrotgscsq = rsrotgscsq + wgt*thtot**2 nrgsraw = nrgsraw + 1 nrgs = nrgs + wgt nrgssq = nrgssq + wgt**2 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) call hf2(8011,real(en1),real(thx1),wgt) 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/sqrt(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,' nr=',nr !-KG c write(*,*) 'lsrot=',lsrot,' lsrotsq=',lsrotsq,' nl=',nl !-KG 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 write(*,*) 'Air gap 1 out ',Nairout1 write(*,*) 'Al window 1 in ',Nalin1 write(*,*) 'Al window 1 out ',Nalout1 write(*,*) 'Input Coil output ',Nincoilout write(*,*) 'WE target input ',NtarWEin write(*,*) 'Up East target IN ',Ntar1 write(*,*) 'Up East target OUT ',Ntar1out write(*,*) 'Up West empty IN ',Nemp1in write(*,*) 'Up West empty Out ',Nemp1out write(*,*) 'Dn East empty IN ',Nemp2in write(*,*) 'DN East empty Out ',Nemp2out write(*,*) 'Dn West target IN ',Ntar2 write(*,*) 'Dn West target OUT ',Ntar2out write(*,*) 'E-W Target Outputs ',NTargetout write(*,*) 'Output Coil output ',NOutcoilout write(*,*) 'ASM output (geom. only) ', cot write(*,*) 'ASM out scattered ',NEdet write(*,*) write(3,*) 'input ',nevts write(3,*) 'PSM output ',NPSMout write(3,*) 'Air gap 1 out ',Nairout1 write(3,*) 'Al window 1 out ',Nalout1 write(3,*) 'Input Coil output ',Nincoilout write(3,*) 'WE target input ',NtarWEin write(3,*) 'Up East target IN ',Ntar1 write(3,*) 'Up East target OUT ',Ntar1out write(3,*) 'Up West empty IN ',Nemp1in write(3,*) 'Up West empty Out ',Nemp1out write(3,*) 'Dn East empty IN ',Nemp2in write(3,*) 'DN East empty Out ',Nemp2out write(3,*) 'Dn West target IN ',Ntar2 write(3,*) 'Dn West target OUT ',Ntar2out write(3,*) 'E-W Target Outputs ',NTargetout write(3,*) 'Output Coil output ',NOutcoilout write(3,*) 'ASM output (geom. only) ', cot write(3,*) 'ASM out scattered ',NEdet write(3,*) write(*,*) 'Number scattered in Al window 1 ',Nalscat1 write(3,*) 'Number scattered in Al window 1 ',Nalscat1 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 East ', & sqrt((rsrotsq-rsrot**2/nr)/nr) write(6,*) ' St Err All East ', & sqrt((rsrotsq-rsrot**2/nr)/nr)*Sqrt(nrsq)/nr c write(6,*) ' nr ',nr , ' nrsq ',nrsq c write(6,*) 'sqrt(nrsq)/nr ', Sqrt(nrsq)/nr if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then c write(6,*) write(6,*) ' StDev Tar Scat East ', & sqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) write(6,*) ' St Err All East ', & sqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*Sqrt(nrtssq)/nrts c write(6,*) ' nrts ',nrts , ' nrtssq ',nrtssq c write(6,*) 'sqrt(nrtssq)/nrts ', Sqrt(nrtssq)/nrts c write(6,*) end if if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then write(6,*) ' StDev Gas Scat East ', & sqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) write(6,*) ' St Err All East ', & sqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*Sqrt(nrgssq)/nrgs c write(6,*) ' nrgs ',nrgs , ' nrgssq ',nrgssq c write(6,*) 'sqrt(nrgssq)/nrgs ', Sqrt(nrgssq)/nrgs end if write(6,*) write(6,*) write(6,713) nr, nl write(6,714) nr/nl write(6,717) nr-nl,Sqrt(nlsq+nrsq) write(6,*) write(6,770) nrnsc , nlnsc write(6,*) write(6,700) mag_B, bzgradnum write(6,704) rsrot/nr, & sqrt((rsrotsq-rsrot**2/nr)/nr)*Sqrt(nrsq)/nr write(6,703) lsrot/nl, & sqrt((lsrotsq-lsrot**2/nl)/nl)*Sqrt(nlsq)/nl write(6,711) rsrot/nr-lsrot/nl, & sqrt( (lsrotsq-lsrot**2/nl)/nl*nlsq/nl**2 & + (rsrotsq-rsrot**2/nr)/nr*nrsq/nr**2 ) if ((nrall.gt.1).and.(nlall.gt.1)) then write(6,*) write(6,804) rall/float(nrall), ! alternate weighting method & sqrt((rallsq-rall**2/nrall)/nrall/(nrall-1)) write(6,803) lall/nlall, & sqrt((lallsq-lall**2/nlall)/nlall/(nlall-1)) write(6,811) rall/nrall-lall/nlall, & sqrt( (lallsq-lall**2/nlall)/nlall/(nlall-1) & + (rallsq-rall**2/nrall)/nrall/(nlall-1) ) else write(6,*) 'No alternate weighting data ' end if write(6,*) write(6,726) rsrotnsc/nrnsc, & sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*Sqrt(nrnscsq)/nrnsc write(6,725) lsrotnsc/nlnsc, & sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc)*Sqrt(nlnscsq)/nlnsc write(6,732) rsrotnsc/nrnsc - lsrotnsc/nlnsc, & sqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc*nlnscsq/nlnsc**2 & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc*nrnscsq/nrnsc**2) write(6,*) if ((nrsc.gt.1).and.(nlsc.gt.1)) then write(6,806) rsc/float(nrsc), ! alternate weighting method & sqrt((rscsq-rsc**2/nrsc)/nrsc/(nrsc-1)) write(6,805) lsc/nlsc, & sqrt((lscsq-lsc**2/nlsc)/nlsc/(nlsc-1)) write(6,812) rsc/nrsc-lsc/nlsc, & sqrt( (lscsq-lsc**2/nlsc)/nlsc/(nlsc-1) & + (rscsq-rsc**2/nrsc)/nrsc/(nlsc-1) ) else write(6,*) ' No alternate weighting scatters' end if write(6,*) if (nlts.gt.0) then write(6,706) rsrottsc/nrts, & sqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*Sqrt(nrtssq)/nrts write(6,705) lsrottsc/nlts, & sqrt((lsrottscsq-lsrottsc**2/nlts)/nlts)*Sqrt(nltssq)/nrts write(6,712) rsrottsc/nrts - lsrottsc/nlts , & sqrt( (lsrottscsq-lsrottsc**2/nlts)/nlts*nltssq/nlts**2 & + (rsrottscsq-rsrottsc**2/nrts)/nrts*nrtssq/nrts**2 ) write(6,*) if (nlgs.gt.0) then write(6,7006) rsrotgsc/nrgs, & sqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*Sqrt(nrgssq)/nrgs write(6,7005) lsrotgsc/nlgs, & sqrt((lsrotgscsq-lsrotgsc**2/nlgs)/nlgs)*Sqrt(nlgssq)/nrgs write(6,7012) rsrotgsc/nrgs - lsrotgsc/nlgs , & sqrt( (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs*nlgssq/nlgs**2 & + (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs*nrgssq/nrgs**2 ) end if write(6,*) write(6,702) rsde/nder, & sqrt((rsdesq-rsde**2/nder)/nder)*Sqrt(ndersq)/nder write(6,701) lsde/ndel, & sqrt((lsdesq-lsde**2/ndel)/ndel)*Sqrt(ndelsq)/ndel if ((nderw.gt.1).and.(ndelw.gt.1)) then write(6,802) rwde/nderw, & sqrt((rwdesq-rwde**2/nderw)/nderw/(nderw-1)) write(6,801) lwde/ndelw, & sqrt((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,Sqrt(nltssq+nrtssq) write(6,782) nrtsraw,nltsraw write(6,*) write(6,7015) nrgs, nlgs write(6,7016) nrgs/nlgs write(6,7018) nrgs-nlgs,Sqrt(nlgssq+nrgssq) write(6,7082) nrgsraw,nlgsraw write(6,*) if ((nrsc.gt.0).and.(nlsc.gt.0)) then write(6,815) float(nrsc), float(nlsc) write(6,816) float(nrsc)/float(nlsc) write(6,818) float(nrsc)-float(nlsc),Sqrt(float(nlsc+nrsc)) end if 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 480.6~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,' Avg. Rot West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 704 format(3x,' Avg. Rot East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 803 format(3x,' Avg. Rot West Dn (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 804 format(3x,' Avg. Rot East Up (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 705 format(3x,' Avg. Rot Target Scat West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 706 format(3x,' Avg. Rot Target Scat East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 805 format(3x,' Avg. Rot Scat West Dn (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 806 format(3x,' Avg. Rot Scat East Up (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 7005 format(3x,' Avg. Rot Gas Scat West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 7006 format(3x,' Avg. Rot Gas Scat East Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 725 format(3x,' Avg. Rot NOScat West Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 726 format(3x,' Avg. Rot 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,' Rotation UpE-DnW (mrad) ',1pe11.4, & ' +- ',1pe11.4) 811 format(3x,' Rotation UpE-DnW (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 712 format(3x,' Rotation UpE-DnW Tar scat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 812 format(3x,' Rotation UpE-DnW Scat (alt wgt) (mrad) ',1pe11.4, & ' +- ',1pe11.4) 7012 format(3x,' Rotation UpE-DnW Gas scat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 732 format(3x,' Rotation UpE-DnW NOScat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 713 format(3x,' # in DetE Up',1pe13.6,2x,', # in Det W 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 West: ',1pe11.4, & ' East: ',1pe11.4) 755 format(1x,' Frac # of scatters with Ef.ne.Ei tar1 ',1pe11.4) 756 format(1x,' Frac # of scatters with Ef.ne.Ei tar2 ',1pe11.4) 757 format(1x,' Frac num of tar1 entering tar2 ',1pe11.4) 758 format(1x,' Frac # of total reaching det ',1pe11.4) 759 format(1x,' Frac 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) 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 write(3,760) nlts/nl write(3,761) nrgs/nr write(3,762) nlgs/nl write(3,*) write(3,*) ' StDev All East ', & sqrt((rsrotsq-rsrot**2/nr)/nr) write(3,*) ' St Err All East ', & sqrt((rsrotsq-rsrot**2/nr)/nr)*Sqrt(nrsq)/nr c write(3,*) ' nr ',nr , ' nrsq ',nrsq c write(3,*) 'sqrt(nrsq)/nr ', Sqrt(nrsq)/nr if ((rsrottscsq-rsrottsc**2/nrts)/nrts.gt.0) then c write(3,*) write(3,*) ' StDev Tar Scat East ', & sqrt((rsrottscsq-rsrottsc**2/nrts)/nrts) write(3,*) ' St Err All East ', & sqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*Sqrt(nrtssq)/nrts c write(3,*) ' nrts ',nrts , ' nrtssq ',nrtssq c write(3,*) 'sqrt(nrtssq)/nrts ', Sqrt(nrtssq)/nrts end if if ((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs.gt.0) then c write(3,*) write(3,*) ' StDev Gas Scat East ', & sqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs) write(3,*) ' St Err All East ', & sqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*Sqrt(nrgssq)/nrgs c write(3,*) ' nrgs ',nrgs , ' nrgssq ',nrgssq c write(3,*) 'sqrt(nrgssq)/nrgs ', Sqrt(nrgssq)/nrgs end if write(3,*) write(3,*) write(3,713) nr, nl write(3,714) nr/nl write(3,717) nr-nl,Sqrt(nlsq+nrsq) write(3,*) write(3,770) nrnsc , nlnsc write(3,*) write(3,700) mag_B, bzgradnum write(3,704) rsrot/nr, & sqrt((rsrotsq-rsrot**2/nr)/nr)*Sqrt(nrsq)/nr write(3,703) lsrot/nl, & sqrt((lsrotsq-lsrot**2/nl)/nl)*Sqrt(nlsq)/nl write(3,711) rsrot/nr-lsrot/nl, & sqrt( (lsrotsq-lsrot**2/nl)/nl*nlsq/nl**2 & + (rsrotsq-rsrot**2/nr)/nr*nrsq/nr**2 ) if ((nrall.gt.1).and.(nlall.gt.1)) then write(3,*) write(3,804) rall/float(nrall), ! alternate weighting method & sqrt((rallsq-rall**2/nrall)/nrall/(nrall-1)) write(3,803) lall/nlall, & sqrt((lallsq-lall**2/nlall)/nlall/(nlall-1)) write(3,811) rall/nrall-lall/nlall, & sqrt( (lallsq-lall**2/nlall)/nlall/(nlall-1) & + (rallsq-rall**2/nrall)/nrall/(nlall-1) ) else write(3,*) 'No alternate weighting data ' end if write(3,*) write(3,726) rsrotnsc/nrnsc, & sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc)*Sqrt(nrnscsq)/nrnsc write(3,725) lsrotnsc/nlnsc, & sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc)*Sqrt(nlnscsq)/nlnsc write(3,732) rsrotnsc/nrnsc - lsrotnsc/nlnsc, & sqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc*nlnscsq/nlnsc**2 & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc*nrnscsq/nrnsc**2) write(3,*) if ((nrsc.gt.1).and.(nlsc.gt.1)) then write(3,806) rsc/float(nrsc), ! alternate weighting method & sqrt((rscsq-rsc**2/nrsc)/nrsc/(nrsc-1)) write(3,805) lsc/nlsc, & sqrt((lscsq-lsc**2/nlsc)/nlsc/(nlsc-1)) write(3,812) rsc/nrsc-lsc/nlsc, & sqrt( (lscsq-lsc**2/nlsc)/nlsc/(nlsc-1) & + (rscsq-rsc**2/nrsc)/nrsc/(nlsc-1) ) else write(3,*) ' No alternate weighting scatters' end if write(3,*) if (nlts.gt.0) then write(3,706) rsrottsc/nrts, & sqrt((rsrottscsq-rsrottsc**2/nrts)/nrts)*Sqrt(nrtssq)/nrts write(3,705) lsrottsc/nlts, & sqrt((lsrottscsq-lsrottsc**2/nlts)/nlts)*Sqrt(nltssq)/nrts write(3,712) rsrottsc/nrts - lsrottsc/nlts , & sqrt( (lsrottscsq-lsrottsc**2/nlts)/nlts*nltssq/nlts**2 & + (rsrottscsq-rsrottsc**2/nrts)/nrts*nrtssq/nrts**2 ) write(3,*) if (nlgs.gt.0) then write(3,7006) rsrotgsc/nrgs, & sqrt((rsrotgscsq-rsrotgsc**2/nrgs)/nrgs)*Sqrt(nrgssq)/nrgs write(3,7005) lsrotgsc/nlgs, & sqrt((lsrotgscsq-lsrotgsc**2/nlgs)/nlgs)*Sqrt(nlgssq)/nrgs write(3,7012) rsrotgsc/nrgs - lsrotgsc/nlgs , & sqrt( (lsrotgscsq-lsrotgsc**2/nlgs)/nlgs*nlgssq/nlgs**2 & + (rsrotgscsq-rsrotgsc**2/nrgs)/nrgs*nrgssq/nrgs**2 ) end if write(3,*) write(3,702) rsde/nder, & sqrt((rsdesq-rsde**2/nder)/nder)*Sqrt(ndersq)/nder write(3,701) lsde/ndel, & sqrt((lsdesq-lsde**2/ndel)/ndel)*Sqrt(ndelsq)/ndel if ((nderw.gt.1).and.(ndelw.gt.1)) then write(3,802) rwde/nderw, & sqrt((rwdesq-rwde**2/nderw)/nderw/(nderw-1)) write(3,801) lwde/ndelw, & sqrt((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,Sqrt(nltssq+nrtssq) write(3,782) nrtsraw,nltsraw write(3,*) write(3,7015) nrgs, nlgs write(3,7016) nrgs/nlgs write(3,7018) nrgs-nlgs,Sqrt(nlgssq+nrgssq) write(3,7082) nrgsraw,nlgsraw write(3,*) if ((nrsc.gt.0).and.(nlsc.gt.0)) then write(3,815) float(nrsc), float(nlsc) write(3,816) float(nrsc)/float(nlsc) write(3,818) float(nrsc)-float(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(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(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(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(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(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(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