program main_WUp c same as main_slim (main_cwn with no 1000s or 2000s histograms) except targets on opposite side c West up (x<0) and East dn (x>0) c calculating path length and rotation from entire flight path, not just while in target implicit none include 'const.inc' integer loop, numbounce double precision ep, wp, kp, k1, thp,php,qlast integer i, j, counter,lim,Event,cin,cot,trig !bec integer Nenscat1,Ntar1,Nenscat2,Ntar2,Ntar12,NEdet integer encnt, Nsdet,Nthdet, Ntar1out, Ntar2out integer Nalscat1, Nalscat2, Nthscat1, Nthscat2 integer Nalout1, Nalout2, Nalin1, Nairin1, Nairout1, Nairscat1 integer NPSMout, Nincoilout, NOutcoilout, NTargetout ! bec integer numtar1flg, numtar2flg, numendif1, numendif2,numendifdet c real xlam, HRNDM1, vect(6), lm(3), const, xinit, yinit, plr,xtest character*28 outfname double precision th0, ph0, thindex, enstart, entar1, entar2 double precision en0, x0, y0, z0, thx0, thy0 double precision en1, x1, y1, z1, thx1, thy1 double precision time1 logical tar2flg,enflg1,enflg2,thflg1,thflg2,alflg1,alflg2 logical airflg1, histflg, histflg2, empflg1, empflg2,bounceflg integer ranno integer seedc integer n1, n2, n3, Nempty1, Nempty2 C parameter(ranno=1) double precision r2d, lambda, sig, intlen, thscat, sigsom, sigsq 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, lsrotsc,lsrotscsq double precision lsrotscle,lsrotsclesq,lsde,lsdesq double precision rsrot,rsrotsq, rsrotsc,rsrotscsq double precision rsrotnsc,rsrotnscsq,lsrotnsc,lsrotnscsq double precision rsrotscle,rsrotsclesq,rsde,rsdesq double precision thxinput, thyinput integer nr,nrs,nrsl,nl,nls,nlsl,nder,nders,ndel,ndels,nb integer nrnsc,nlnsc,numE2W, numW2E,NtarEout,NtarWout integer numcross double precision IntBz,temp,norm !for gradient runs -KG CHARACTER day*2,month*2,hour*2,minute*2, filename*20 c----------------------------------- integer qmin, qmax double precision k0, qdummy, q c----------------------------------- c c ------------------------- c fill the ntuple 'spin.rz' //// ntuple variables .... integer max_paw, iquest integer istat, icycle, datetime(8) parameter(max_paw = 10000000) ! 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) open(unit=37,file='Bz.dat',status='old') 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 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 left, x<0. Dn on right ' write(3,*) ' Upstream target on left, x<0. Dn on right ' 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(100,'lam.hbook','') c call hrget(102,'xyspot.hbook','') call hrget(151,'gradx.hbook','') call hrget(152,'gradx.hbook','') call hrget(301,'sq.hbook','') call hcopyr(301,201,"",1,100,0.,0.,"") c----------------------------------------------------------------------- c read(*,*) ranno call rluxgo(3,ranno,0,0) ! RANLUX luxury level 2 cc call rluxgo(3,710,0,0) ! RANLUX luxury level 2 c----------------------------------------------------------------------- c Print *, "setting up Hbooks" call hbook1(305,'West q sim(1/Ang)',100,0,4.0,0.) ! bec test q choices call hbook1(306,'West w sim(1/psec)',100,-4.0,4.0,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,4.0,0.) ! bec test q choices call hbook1(206,'East w sim(1/psec)',100,-4.0,4.0,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) ',100,-3.5,3.5,0.) call hbook1(1,'x (cm) start',100,-3.5,3.5,0.) call hbook1(2,'x (cm) guide out',100,-3.5,3.5,0.) call hbook1(3,'x (cm) incoil in',100,-3.5,3.5,0.) call hbook1(4,'x (cm) incoil out',100,-3.5,3.5,0.) call hbook1(5,'x (cm) outcoil out',100,-3.5,3.5,0.) call hbook1(6,'x (cm) ASM out',100,-3.5,3.5,0.) call hbook1(20,'E (meV) start Both',2000,0.,40.,0.) call hbook1(1002,'y (cm) UpWtar',100,-3.,3.,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 WE',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,-3.5,3.5,0.) call hbook1(2002,'y (cm) Dn',100,-3.,3.,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,-3.,3.,0.) call hbook1(3003,'thx (mrad) Det',200,0.,0.1,0.) call hbook1(3004,'thy (mrad) Det',200,0.,0.1,0.) call hbook1(3005,'E (meV) DetWE',2000,0.,40.,0.) call hbook1(4005,'Lambda (Angs) DetWE',1000,0.,50.,0.) c call hbook1(3006,'Path length (cm) Det WE',200,480.58,480.76,0.) !KG c call hbook1(3007,'th_rotation (mrad) DetWE',200,-400.,0.0,0.) !KG call hbook1(3008,'Path length (cm) DetEast',300,480.68, & 481.0,0.) call hbook1(3010,'Path length (cm) DetWest',300,480.68, & 481.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.) call hbook1(3015,'E (meV) scat DetE',2000,0.,40.,0.) call hbook1(3016,'E (meV) scat DetW',2000,0.,40.,0.) call hbook1(3017,'E (meV) Non-scat DetE',2000,0.,40.,0.) call hbook1(3018,'E (meV) Non-scat DetW',2000,0.,40.,0.) c call hbook1(3020,'E_loss (meV) scat DetWE',2000,-3.,3.,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)',100,0.,4., & 200,0.,1.,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)',200,0.,180.,0.) call hbook1(4013,'West sig Sommers (b/atom)',200,0.,1.,0.) call hbook1(4014,'West sig Sq (b/atom)',200,0.,1.,0.) call hbook1(4023,'East sig Sommers (b/atom)',200,0.,1.,0.) call hbook1(4024,'East sig Sq (b/atom)',200,0.,1.,0.) call hbook1(4015,'theta-p (deg) U',200,0.,180.,0.) call hbook1(4016,'phi-p (deg) U',200,0.,360.,0.) call hbook1(4017,'th-0 (deg)',200,0.,10.,0.) 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.) 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.) call hbook1(5021,'phi (deg) af-guide',200,0.,360.,0.) call hbook1(5022,'phi (deg) af-target 1',200,0.,360.,0.) call hbook1(5023,'phi (deg) af-target 2',200,0.,360.,0.) call hbook1(6001,'Path length (cm) enter Up WE',200,213.0, & 215.0,0.) c call hbook1(6003,'Path length (cm) Det scat WE',300,480.68, c & 481.0,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.68, & 481.0,0.) call hbook1(6007,'Path length (cm) Det scat East',300,480.68, & 481.0,0.) c call hbook1(6009,'Path length (cm) Det no scat WE',300,480.68, c & 481.0,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.68, & 481.0,0.) call hbook1(6013,'Path length (cm) Det no scat West',300,480.68, & 481.0,0.) call hbook1(6015,'Crossover Path length (cm) Det scat West', & 300,480.68,481.0,0.) call hbook1(6017,'Crossover Path length (cm) Det scat East', & 300,480.68,481.0,0.) call hbook1(6035,'Crossover Path length (cm) DetEast', & 300,480.68,481.0,0.) call hbook1(6037,'Crossover Path length (cm) DetWest', & 300,480.68,481.0,0.) c c Bz gradient and Pi-coil adjustments c if (bzgradflg.and.picoilflg) then ! Bz gradient, Picoil ON call hbook1(3009,'th_rotation (mrad) DetEast',300,-800.,10.0,0.) 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,-5.,50.,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.68,481.0,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.68,481.0,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.68,481.0,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.68,481.0,300,-800.,0,0.) else if (bzgradflg.and.(.not.picoilflg)) then ! Bz gradient, Picoil OFF call hbook1(3009,'th_rotation (mrad) DetEast',300,-10.,80.0,0.) 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(6003,'th_rotsc (mrad) vs. E(meV) det West',800,0.,40., & 300,0.,60.,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,0.,60.,0.) call hbook2(6021,'th rot (mrad) vs. Path len (cm) det West',300, & 480.68,481.0,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.68,481.0,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.68,481.0,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.68,481.0,300,0.,60,0.) else if ((.not.bzgradflg).and.(picoilflg)) then ! NO Bz gradient, Picoil ON 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,-10.,10.,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.68,481.0,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.68,481.0,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.68,481.0,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.68,481.0,300,-10.,10,0.) else if ((.not.bzgradflg).and.(.not.picoilflg)) then ! 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.68,481.0,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.68,481.0,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.68,481.0,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.68,481.0,300,-5.,50,0.) end if c call hbook2(6022,'th rot (mrad) vs. E (meV) detW',1000,0.,40., c & 200,-600.,0.,0.) c call hbook2(6023,'th rot (mrad) vs. Path len (cm) detW',200,480.68,481.00, c & 200,-600.,0.,0.) c call hbook2(6024,'th rot (mrad) vs. E (meV) detE',1000,0.,40., c & 200,-600.,0.,0.) c call hbook2(6025,'th rot (mrad) vs. Path len (cm) detE',200,480.68,481.00, c & 200,-600.,0.,0.) c call hbook2(6026,'E (meV) vs. Path len (cm) detW',200,480.68,481.00, c & 1000,0.,40.,0.) c call hbook2(6027,'E (meV) vs. Path len (cm) detE',200,480.68,481.00, c & 1000,0.,40.,0.) c call hbook2(8011,'Thx vs. E (meV) scat DetWE',1000,0.,40., c & 200,0.,0.1,0.) c call hbook2(8012,'Thx vs. E-loss (meV) scat DetWE',1000,0.,40.,200, c & 0.,0.1,0.) 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 NTargetout=0; Ntar1out=0; Ntar2out=0 !bec Nthdet=0 NPSMout=0 ; Nincoilout=0 ; NOutcoilout = 0 !bec lsrot=0;lsrotsq=0;lsrotsc=0;lsrotscsq=0;lsrotscle=0;lsrotsclesq=0 lsde=0;lsdesq=0;rsrot=0;rsrotsq=0;rsrotsc=0;rsrotscsq=0 nrnsc=0;nlnsc=0;rsrotnsc=0;lsrotnsc=0;rsrotnscsq=0;lsrotnscsq=0 rsrotscle=0 rsrotsclesq=0;rsde=0;rsdesq=0;nr=0;nrs=0;nrsl=0;nl=0;nls=0;nlsl=0 nder=0;nders=0;ndel=0;ndels=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 ; numendifdet = 0 numE2W=0 ; numW2E=0 ; NtarEout = 0 ;NtarWout = 0 numcross=0 ; numbounce = 0 c c c write(21,*) '0.41 ' goto 321 end if th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec if (th0.lt.1.E-8) th0=1.E-8 ph0 = 2.0d0*pi*vect(2) ! azimuthal angle (radians) call hf1(4017,real(th0*180./3.14159),1.) call hf1(4018,real(ph0*180./3.14159),1.) c call HRNDM2(102,xinit,yinit) c x0=xinit/10. ! cm c y0=yinit/10. ! cm c write(18,*) x0, y0; goto 100 xtest = xin*(2.*vect(3)-1.) c if (i.eq.1) write(6,*) '!!!!!fake x gradient!!!!!!!!!!!!!!' c if (i.eq.1) write(3,*) '!!!!!fake x gradient!!!!!!!!!!!!!!' c if (vect(6).gt.0.3) then c xtest = xin*vect(3) ! positive x values c else c xtest = -xin*vect(3) ! negative x values c endif c x0=xtest c if (i.eq.1) write(*,*) ' ! Running with No x gradient !' if (i.eq.1) write(3,*) ' ! Running with No x gradient !' x0 = xtest c c if (i.eq.1) write(*,*) ' ! Running with R to L x gradient !' c if (i.eq.1) write(3,*) ' ! Running with R to L x gradient !' c x0 = 3.0d0 - 2.d0*HRNDM1(151) !bec R to L gradient c c if (i.eq.1) write(*,*) ' ! Running with L to R x gradient !' c if (i.eq.1) write(3,*) ' ! Running with L to R x gradient !' c x0 = 2.d0*HRNDM1(151) - 3.0d0 !bec L to R gradient c c this gives double gradient -- gradient across each half... c if(xtest.gt.0.) then c x0 = HRNDM1(151) ! gives gradient only acrros half cc x0 = 3.d0 - HRNDM1(151) c elseif(xtest.lt.0) then c x0 = HRNDM1(151) - 3. cc x0 = -HRNDM1(151) c endif c c y0 = 2.25*(2.*vect(4)-1.) c square cross section y0 = yin*(2.*vect(4)-1.) z0= 0.0d0 if (i.eq.1) then if (bzgradflg) then write(6,*) ' !!!!Running with Bz gradient in target region!!!!' write(3,*) ' !!!!Running with Bz gradient in target region!!!!' else write(6,*) ' ! NO Bz gradient in target region !' write(3,*) ' ! NO Bz gradient in target region !' 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 c Calculate the incident angles in the xz-plane and yz-plane c note: these are not "direction cosines"= angle between momentum and x (or y) c thx0 = datan(dtan(th0)*dcos(ph0)) thy0 = datan(dtan(th0)*dsin(ph0)) c thxinput = thx0 thyinput = thy0 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 * c x1 = x0 + tzair1*dtan(thx0) c y1 = y0 + tairz1*dtan(thy0) c z1 = z0 + tairz1 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 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(18,*) '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 c write(6,*) c write(6,*) ' PSM out thtot ', thtot, ', pl_del ',pl_del ******************************************************************** * 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 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(18,*) '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) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 call hf1(2,real(x1),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 = 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(18,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' guide thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if c write(6,*) ' Guideout thtot ', thtot, ', pl_del ',pl_del * ******************************************************************** ******************************************************************** * 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) 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.) 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 = 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(18,*) 'incoil ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if * c write(6,*) ' IC out thtot ', thtot, ', pl_del ',pl_del ******************************************************************** ******************************************************************** * 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 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(18,*) '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 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(18,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if c if ( ((x1.gt.0.0).and.(xtest.lt.0.0)).or. c & ((x1.lt.0.0).and.(xtest.gt.0.0)) ) then c numcross = numcross + 1 c goto 100 c endif ******************************************************************** ******************************************************************** 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. Once it enters target c target.f checks the y-direction too, so might have neutrons that get into target.f c that don't run target.f really, since the y-values are out of range, hence sig or thscat values c get 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 write(*,*) 'tar1 starts at z = ',z0 !-KG c Given current energy determine q from S(q) and delE from omega(q) k0 = sqrt(2*m_n*en0)/h_par qmin = 0 c qmax = int(k0*25.) ! for sq_1.hbook qmax = int(2.0*k0*25.) ! bec qmax = 2*k from p-consv. if((2.0*k0).gt.3.6) qmax=int(3.6*25.) c Check to see if we've been here before and thus created HIST 302 c histflg set false before loop, then always true if(histflg) call hdelet(302) histflg = .true. call hcopyr(301,302,"",qmin,qmax,0.,0.,"") c c qdummy=HRNDM1(301) ! bec c HRNDM1 replaces spec with integrated spec!, so next pass through the loop c hcopyr is copying the integrated spec into 302, not the original spec!! c without the qdummy=HRNDM1(301) line, the simulated q-spec (305) looks c much more reasonable! bec loop=0 123 q=HRNDM1(302) if (loop.gt.1000) then write(6,*) ' 1 q-loop ',loop,en0,q write(6,*) k0,kp,ep end if loop=loop+1 call wq(q,wp,en0,k0,vect(5)) ! wq modified to prevent E,p consv. probs ep = en0- (h_par*dabs(wp)) !-KG if (ep.le.0.) goto 123 kp=sqrt(2*m_n*ep)/h_par if (dabs((k0**2+kp**2-q**2)/(2.*k0*kp)).le.1.d0) goto 124 goto 123 !!remove the check right above here -KG 124 continue call hf1(305,real(q),1.) call hf1(306,real(wp),1.) call hf1(307,real(en0),1.) call hf1(308,real(ep),1.) call hf1(309,real(k0),1.) call hf1(310,real(kp),1.) c enflg1 = .false. thflg1 = .false. empflg1 = .false. if( (x0.ge.tarx3.and.x0.le.tarx4) ) then ! West=pos. x, target UPstream sig = 0.0 ; sigsom = 0.0 ; sigsq =0.0 ; thscat = 0.0 thp = 0.0; php = 0.0 call target(q, en0, x0, y0, z0, thx0, thy0, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1, sig, & thscat,sigsom,sigsq,enflg1,thflg1,thp,php,Ntar1, & pl1a_sc,pl1b_sc,vect(5)) ! bec 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+1 if (enflg1) Nenscat1 = Nenscat1+1 ! count scatters with energy change bec if (thflg1) Nthscat1 = Nthscat1+1 ! count scatters with theta change bec c capture diagnostics for target interaction stuff bec call hf1(4010,real(sig),1.) call hf2(8010,real(q),real(sig),1.) !temporary c call hf1(4011,real(intlen),1.) call hf1(4012,real(180./3.14159*thscat),1.) call hf1(4013,real(sigsom),1.) call hf1(4014,real(sigsq),1.) call hf1(4015,real(180./3.14159*thp),1.) call hf1(4016,real(180./3.14159*php),1.) end if if ((x1.ge.tarx3.and.x1.le.tarx4).and. & (y1.ge.tary1.and.y1.le.tary2)) & Ntar1out = Ntar1out+1 ! out of target west side only elseif( (x0.ge.tarx1.and.x0.le.tarx2) ) then ! East=neg. x, empty UP call empty(en0, x0, y0, z0, thx0, thy0, tarx1, tarx2, & en1, x1, y1, z1, thx1, thy1,empflg1,pl1a_sc,pl1b_sc) if (empflg1) Nempty1 = Nempty1+1 ! count scatters with energy change bec 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),1.) call hf1(1002,real(y1),1.) c call hf1(1003,real(thx1),1.) c call hf1(1004,real(thy1),1.) if(en0.eq.en1) then pl1_ns = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) if (bzgradflg) then ! Bz gradient c Replace mag_b with 4.866mG for gradient run -KG th1_ns = gy*4.866*pl1_ns*0.1/sqrt(2.*en1/m_n) else ! NO Bz gradient th1_ns = gy*mag_B*pl1_ns*0.1/sqrt(2.*en1/m_n) end if pl1_sc = 0. th1_sc = 0. pltot = pl1_ns + pltot thtot = th1_ns + thtot call hf1(1005,real(en1),1.) c call hf1(1006,real(pl1_ns),1.) c call hf1(1007,real(th1_ns),1.) if(x0.gt.0..and.x1.gt.0.) then c call hf1(1008,real(pl1_ns),1.) c call hf1(1009,real(th1_ns),1.) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(1010,real(pl1_ns),1.) c call hf1(1011,real(th1_ns),1.) endif else cc Replace mag_b with 4.862mG for gradient run -KG cc next two lines are old, before int_Bz.f cc th1a_sc = gy*4.862*pl1a_sc*0.1/sqrt(2.*en0/m_n) cc th1b_sc = gy*4.862*pl1b_sc*0.1/sqrt(2.*en1/m_n) c th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) c if (bzgradflg) then ! Bz gradient c 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) cccc IntBz=165.01(mG*mm)@23.7mm where 1st Tgt begins, 2189.39 @439.7mm at the end else ! NO Bz gradient 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) end if c 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),1.) c call hf1(1013,real(pl1_sc),1.) c call hf1(1014,real(th1_sc),1.) if(x0.gt.0..and.x1.gt.0.) then c call hf1(1015,real(pl1_sc),1.) c call hf1(1016,real(th1_sc),1.) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(1017,real(pl1_sc),1.) c call hf1(1018,real(th1_sc),1.) endif endif if (pltot.gt.500) write(18,*) '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 c write(6,*) ' T1 out thtot ', thtot, ', pl_del ',pl_del * ******************************************************************** * 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(18,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot 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),1.) call hf1(5020,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 call hf1(5022,real(180./3.14159*phi),1.) 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 c if (bzgradflg) then ! Bz gradient cc Replace mag_b with 0.371mG for gradient run -KG thtot = gy*0.371*pl_del*0.1/sqrt(2.*en1/m_n) + thtot else ! NO Bz gradient thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot end if if (thtot.gt.1e4) 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 (bzgradflg) then ! Bz gradient cc Replace mag_b with -0.803mG for gradient run -KG thtot = gy*(-0.803)*pl_del*0.1/sqrt(2.*en1/m_n) + thtot else ! NO Bz gradient thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot end if if (pltot.gt.500) write(18,*) 'gap 2 ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if c write(6,*) ' T1gap out thtot ', thtot, ', pl_del ',pl_del * ******************************************************************** * 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(18,*) 'guide ',pltot if (thtot.gt.1e9) then write(6,*) ' Al thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if ******************************************************************** ******************************************************************** 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.tarx3.and.x0.lt.tarx4) ) then ! West = pos. x, empty down call empty(en1, x0, y0, z0, thx0, thy0, tarx3, tarx4, & en1, x1, y1, z1, thx1, thy1, empflg2,pl2a_sc,pl2b_sc) if (empflg2) Nempty2 = Nempty2+1 ! count scatters with energy change bec elseif( (x0.gt.tarx1.and.x0.lt.tarx2) ) then ! East = neg. x, target DOWN if (enflg1.or.thflg1) then ! if scat in t1, now scat in both targets Ntar12=Ntar12+1 endif sig = 0.0 ; sigsom = 0.0 ; sigsq =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 qmin = 0 c qmax = int(k0*25.) ! for sq_1.hbook qmax = int(2.0*k0*25.) ! bec qmax = 2*k from p-consv. if((2.0*k0).gt.3.6) qmax=int(3.6*25.) if(histflg2) call hdelet(202) histflg2 = .true. call hcopyr(201,202,"",qmin,qmax,0.,0.,"") c loop=0 223 q=HRNDM1(202) if (loop.gt.1000) then write(6,*) ' 2 q-loop ',loop,en0,q write(6,*) k0,kp,ep end if loop=loop+1 call wq(q,wp,en0,k0,vect(5)) ! wq modified to prevent E,p consv. probs ep = en0- (h_par*dabs(wp)) !-KG if (ep.le.0.) goto 223 kp=sqrt(2*m_n*ep)/h_par if (dabs((k0**2+kp**2-q**2)/(2.*k0*kp)).le.1.d0) goto 224 goto 223 !!remove the check right above here -KG 224 continue call hf1(205,real(q),1.) call hf1(206,real(wp),1.) call hf1(207,real(en0),1.) call hf1(208,real(ep),1.) call hf1(209,real(k0),1.) call hf1(210,real(kp),1.) call target(q, en0, x0, y0, z0, thx0, thy0, tarx1, tarx2, & en1, x1, y1, z1, thx1, thy1, sig, & thscat,sigsom,sigsq,enflg2,thflg2,thp,php,Ntar2, & pl2a_sc,pl2b_sc,vect(5)) ! bec 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+1 if (enflg2) Nenscat2 = Nenscat2+1 ! count scatters with ef.ne.ei bec if (thflg2) Nthscat2 = Nthscat2+1 ! count scatters with thsc>0 bec c capture diagnostics for target interaction stuff bec call hf1(4020,real(sig),1.) c call hf1(4021,real(intlen),1.) call hf1(5012,real(180./3.14159*thscat),1.) call hf1(4023,real(sigsom),1.) call hf1(4024,real(sigsq),1.) call hf1(5015,real(180./3.14159*thp),1.) call hf1(5016,real(180./3.14159*php),1.) end if if ((x1.ge.tarx1.and.x1.le.tarx2).and. & (y1.ge.tary1.and.y1.le.tary2)) Ntar2out = Ntar2out+1 ! out of target east 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),1.) call hf1(2002,real(y1),1.) c call hf1(2003,real(thx1),1.) c call hf1(2004,real(thy1),1.) if(.not.enflg2) then pl2_ns = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2) if (bzgradflg) then ! Bz gradient c Replace mag_b with -4.864mG for gradient run -KG th2_ns = gy*(-4.864)*pl2_ns*0.1/sqrt(2.*en1/m_n) else ! NO Bz gradient th2_ns = gy*mag_B*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),1.) c call hf1(2006,real(pl2_ns),1.) c call hf1(2007,real(th2_ns),1.) if(x0.gt.0..and.x1.gt.0.) then c call hf1(2008,real(pl2_ns),1.) c call hf1(2009,real(th2_ns),1.) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(2010,real(pl2_ns),1.) c call hf1(2011,real(th2_ns),1.) endif else cc Replace mag_b with -4.857mG for gradient run -KG cc next two lines are old, before int_Bz.f cc th2a_sc = gy*(-4.857)*pl2a_sc*0.1/sqrt(2.*en0/m_n) cc th2b_sc = gy*(-4.857)*pl2b_sc*0.1/sqrt(2.*en1/m_n) th0 = datan(sqrt(dtan(thx0)**2 + dtan(thy0)**2)) if (bzgradflg) then ! Bz gradient 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) cccc IntBz=2169.38(mG*mm)@521.7mm where 2nd Tgt begins, 145.89 @937.7mm at the end else ! NO Bz gradient 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) 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),1.) c call hf1(2013,real(pl2_sc),1.) c call hf1(2014,real(th2_sc),1.) if(x0.gt.0..and.x1.gt.0.) then c call hf1(2015,real(pl2_sc),1.) c call hf1(2016,real(th2_sc),1.) elseif(x0.lt.0..and.x1.lt.0.) then c call hf1(2017,real(pl2_sc),1.) c call hf1(2018,real(th2_sc),1.) endif endif if (pltot.gt.500) write(18,*) 'target 2 ',pltot if (thtot.gt.1e4) 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(6,*) ' T2out thtot ', thtot, ', pl_del ',pl_del c write(*,*) 'tar 2 ends at z = ', z1 * ******************************************************************** NTargetout = NTargetout + 1 if ((x1.ge.tarx1.and.x1.le.tarx2).and. & (y1.ge.tary1.and.y1.le.tary2)) NtarEout = NtarEout+1 ! out of target east side only if ((x1.ge.tarx3.and.x1.le.tarx4).and. & (y1.ge.tary1.and.y1.le.tary2)) NtarWout = NtarWout+1 ! out of target east side only ******************************************************************** ******************************************************************** * 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(18,*) '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 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(18,*) '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) if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100 call hf1(5,real(x1),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(18,*) 'out coil ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if c write(6,*) ' OCout thtot ', thtot, ', pl_del ',pl_del * ******************************************************************** ******************************************************************** * 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(18,*) 'gap 5 ',pltot NOutcoilout = NOutcoilout + 1 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+1 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(18,*) 'ASM ',pltot if (thtot.gt.1e9) then write(6,*) ' thtot ', thtot, ', pl_del ',pl_del write(6,*) ' en ',en1, ', pltot ',pltot end if if ( ((thx1.ne.thxinput).or.(thy1.ne.thyinput)).and. &(en1.eq.enstart) ) then numbounce = numbounce + 1 bounceflg = .true. end if c write(6,*) ' ASMout thtot ', thtot, ', En ',en1 c write(6,*) * ******************************************************************** ******************************************************************** ******************************************************************** call hf1(3001,real(x1),1.) call hf1(3002,real(y1),1.) call hf1(3003,real(thx1),1.) call hf1(3004,real(thy1),1.) pl3 = pl1_ns + pl2_ns + pl1_sc + pl2_sc th3 = th1_ns + th2_ns + th1_sc + th2_sc call hf1(3005,real(en1),1.) lambda=(2.*pi*h_par)/sqrt(2*m_n*en1) call hf1(4005,real(lambda),1.) c call hf1(3006,real(pltot),1.) c call hf1(3007,real(thtot),1.) ************Rotations and Histograms for ALL neutrons**************** if(x1.ge.0.) then ! x>0 = West (=old right) call hf1(3010,real(pltot),1.) call hf1(3011,real(thtot),1.) call hf1(3013,real(en1),1.) call hf2(6020,real(enstart),real(thtot),1.) call hf2(6021,real(pltot),real(thtot),1.) if ((xtest.lt.0).and.(x1.gt.0)) then call hf1(6035, real(pltot),1.) call hf1(6036, real(thtot),1.) call hf2(6040,real(enstart),real(thtot),1.) call hf2(6041,real(pltot),real(thtot),1.) end if rsrot = rsrot + thtot rsrotsq = rsrotsq + thtot**2 nr = nr + 1 if (en1.eq.enstart) then rsrotnsc = rsrotnsc + thtot rsrotnscsq = rsrotnscsq + thtot**2 nrnsc = nrnsc + 1 endif else ! x<0 = East (=old left) call hf1(3008,real(pltot),1.) call hf1(3009,real(thtot),1.) call hf1(3012,real(en1),1.) call hf2(6022,real(enstart),real(thtot),1.) call hf2(6023,real(pltot),real(thtot),1.) if ((xtest.gt.0).and.(x1.lt.0)) then call hf1(6037, real(pltot),1.) call hf1(6038, real(thtot),1.) call hf2(6042,real(enstart),real(thtot),1.) call hf2(6043,real(pltot),real(thtot),1.) end if lsrot = lsrot + thtot lsrotsq = lsrotsq + thtot**2 nl = nl + 1 if (en1.eq.enstart) then lsrotnsc = lsrotnsc + thtot lsrotnscsq = lsrotnscsq + thtot**2 nlnsc = nlnsc + 1 endif endif if (entar1.ne.enstart) numendif1 = numendif1 + 1 if (enflg1) numtar1flg = numtar1flg + 1 if (entar2.ne.entar1) numendif2 = numendif2 + 1 if (enflg2) numtar2flg = numtar2flg + 1 if (en1.ne.enstart) numendifdet = numendifdet +1 ****************Rotations and Histograms for Scattered neutrons**************** if(en1.ne.enstart) then ! scattered and changed energy anywhere NEdet=NEdet+1 c call hf1(6003, real(pltot),1.) c call hf1(6004, real(thtot),1.) if (x1.lt.0.) then ! x<0 = East (=old left) call hf2(6031,real(enstart),real(thtot),1.) call hf1(6007, real(pltot),1.) call hf1(6008, real(thtot),1.) if ((xtest.gt.0).and.(x1.lt.0)) then call hf1(6017, real(pltot),1.) call hf1(6018, real(thtot),1.) end if call hf1(3015,real(en1),1.) c call hf2(6024,real(enstart),real(thtot),1.) c call hf2(6025,real(pltot),real(thtot),1.) c call hf2(6027,real(pltot),real(en1),1.) lsrotsc = lsrotsc + thtot lsrotscsq = lsrotscsq + thtot**2 nls = nls + 1 if (en1.le.0.35) then lsrotscle = lsrotscle + thtot lsrotsclesq=lsrotsclesq+thtot**2 nlsl = nlsl + 1 end if c call hf1(3020,real(enstart-en1),1.) call hf1(3021,real(enstart-en1),1.) call hf2(8011,real(en1),real(thx1),1.) call hf2(8012,real(enstart-en1),real(thx1),1.) if (thx1.gt.0.1) write(6,*) ' thx, e, eloss', & thx1,en1,enstart-en1 lsde = lsde + enstart-en1 lsdesq = lsdesq + (enstart-en1)**2 ndel = ndel + 1 c write(6,*) ' 1W eloss tar1 tar2 ', entar1-entar2 else ! x>0 = West (=old right) call hf2(6030,real(enstart),real(thtot),1.) call hf1(6005, real(pltot),1.) call hf1(6006, real(thtot),1.) if ((xtest.lt.0).and.(x1.gt.0)) then call hf1(6015, real(pltot),1.) call hf1(6016, real(thtot),1.) end if call hf1(3016,real(en1),1.) c call hf2(6022,real(enstart),real(thtot),1.) c call hf2(6023,real(pltot),real(thtot),1.) c call hf2(6026,real(pltot),real(en1),1.) rsrotsc = rsrotsc + thtot rsrotscsq = rsrotscsq + thtot**2 nrs = nrs + 1 if (en1.le.0.35) then rsrotscle = rsrotscle + thtot rsrotsclesq=rsrotsclesq+thtot**2 nrsl = nrsl + 1 c write(6,*) ' 1E eloss tar1 tar2 ', entar1-entar2 end if c call hf1(3020,real(enstart-en1),1.) call hf1(3022,real(enstart-en1),1.) call hf2(8011,real(en1),real(thx1),1.) call hf2(8012,real(enstart-en1),real(thx1),1.) if (thx1.gt.0.1) write(6,*) ' thx, e, eloss', & thx1,en1,enstart-en1 rsde = rsde + enstart-en1 rsdesq = rsdesq + (enstart-en1)**2 nder = nder + 1 c write(6,*) ' 2E eloss tar1 tar2 ', entar1-entar2 endif else c call hf1(6009, real(pltot),1.) c call hf1(6010, real(thtot),1.) if (x1.lt.0.) then call hf1(6011, real(pltot),1.) call hf1(6012, real(thtot),1.) call hf1(3018,real(en1),1.) c call hf2(6024,real(entar2),real(thtot),1.) c call hf2(6025,real(pltot),real(thtot),1.) c call hf2(6027,real(pltot),real(en1),1.) else call hf1(6013, real(pltot),1.) call hf1(6014, real(thtot),1.) call hf1(3017,real(en1),1.) c call hf2(6022,real(entar2),real(thtot),1.) c call hf2(6023,real(pltot),real(thtot),1.) c call hf2(6026,real(pltot),real(en1),1.) endif endif ******************************************************************** c if (en1.lt.0.36) then c write(16,888) ' Edet<0.36 -- E0, E1, E2 ', enstart, entar1, entar2 c write(16,889) ' E0-1, E1-2 ', enstart-entar1, entar1-entar2 c write(16,888) ' thx, thy', thx1, thy1 c write(16,889) ' q, thscat ', q, thscat c write(16,889) ' Th-rot, pltot ',thtot, pltot c write(16,*) c endif c 888 format(2x,a30,1x,1pe11.4,1x,1pe11.4,1x,1pe11.4) c 889 format(2x,a30,1x,1pe11.4,1x,1pe11.4) if(thflg1.or.thflg2) Nthdet=Nthdet+1 if(enflg1.or.enflg2.or.thflg1.or.thflg2) then ! bec c call hf1(3014,real(en1),1.) Nsdet=Nsdet+1 endif c if (i.lt.100) write(17,*) ' end ', pltot,thtot ******************************************************************** * if (abs(thtot).gt.300.0) then write(6,941) xtest,x1,enstart,en1,pltot,thtot write(19,941) xtest,x1,enstart,en1,pltot,thtot 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 + 1 if ((xtest.gt.0).and.(x1.lt.0)) numW2E = numW2E + 1 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 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 scat West ',nrs write(6,*) ' num scat East ',nls write(3,*) ' num scat West ',nrs write(3,*) ' num scat East ',nls 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 write(6,*) ' num crosstalk elim. ',numcross write(3,*) ' num crosstalk elim. ',numcross write(6,*) ; write(3,*) write(6,*) ' Fraction of det hits that bounced ', & float(numbounce)/float(cot) write(3,*) ' Fraction of det hits that bounced ', & float(numbounce)/float(cot) print * if (hegastar) then write(6,*) ' 4He gas target ' write(3,*) ' 4He gas target ' elseif (airtar) then write(6,*) ' Air target ' write(3,*) ' Air target ' 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(*,*) 'West target IN ',Ntar1 write(*,*) 'West target OUT ',Ntar1out write(*,*) 'East target IN ',Ntar2 write(*,*) 'East target OUT ',Ntar2out write(*,*) 'E-W Target Outputs ',NTargetout write(3,*) 'West target end out ',NtarWout write(3,*) 'East target end out ',NtarEout write(*,*) 'Output Coil output ',NOutcoilout write(*,*) 'ASM output (geom. only) ', cot 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,*) 'West target IN ',Ntar1 write(3,*) 'West target OUT ',Ntar1out write(3,*) 'East target IN ',Ntar2 write(3,*) 'East target OUT ',Ntar2out write(3,*) 'E-W Target Outputs ',NTargetout write(3,*) 'West target end out ',NtarWout write(3,*) 'East target end out ',NtarEout write(3,*) 'Output Coil output ',NOutcoilout write(3,*) 'ASM output (geom. only) ', cot 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 (En0.ne.En1)',Nenscat1 write(3,*) 'Number scattered in UpTar (En0.ne.En1)',Nenscat1 write(*,*) 'Number scattered in DnTar (En0.ne.En1)',Nenscat2 write(3,*) 'Number scattered in DnTar (En0.ne.En1)',Nenscat2 write(*,*) 'Number scattered in EmptyTar East ',Nempty1 write(*,*) 'Number scattered in EmptyTar West ',Nempty2 write(3,*) 'Number scattered in EmptyTar East ',Nempty1 write(3,*) 'Number scattered in EmptyTar West ',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,*) write(6,*) ' Frac # of TARGET scatters with Ef.ne.Ei tar1 ', & float(Nenscat1)/float(Ntar1) write(6,*) ' Frac # of TARGET scatters with Ef.ne.Ei tar2 ', & float(Nenscat2)/float(Ntar2) write(6,*) ' Frac num of tar1 entering tar2 ', & float(Ntar12)/float(Ntar1) write(6,*) ' Frac of scat in TARGET of total reaching det ', & float(Nsdet)/float(Nevts) write(6,*) ' Frac of det hits that were scattered in TARGET ', & float(Nsdet)/float(cot) write(6,*) ' Frac of det hits that were scattered with thsc>0 ', &'in TARGET ', & float(Nthdet)/float(cot) write(6,*) ' Frac of det hits that were scat in TAR with Ef.ne.Ei ', & float(NEdet)/float(cot) write(6,*) ' Frac of det hits that were scattered somewhere ', & float(numendifdet)/float(cot) write(6,713) float(nr), float(nl) write(6,714) float(nr)/float(nl) write(6,717) float(nr)-float(nl),Sqrt(float(nl+nr)) write(6,*) write(6,700) mag_B cbec write(6,704) 1000*dasin(rsrot/float(nr)), cbec & 1000*dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1))) cbec write(6,703) 1000*dasin(lsrot/float(nl)), cbec & 1000*dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1))) cbec write(6,711) 1000*(dasin(rsrot/float(nr))-dasin(lsrot/float(nl))), cbec & 1000*sqrt(dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))**2 cbec & + dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))**2) write(6,704) rsrot/float(nr), & sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)) write(6,703) lsrot/float(nl), & sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)) write(6,711) rsrot/float(nr)-lsrot/float(nl), & sqrt( (lsrotsq-lsrot**2/nl)/nl/(nl-1) & + (rsrotsq-rsrot**2/nr)/nr/(nr-1) ) write(6,*) cbec write(6,726) 1000*dasin(rsrotnsc/float(nrnsc)), cbec & 1000*dasin(sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1))) cbec write(6,725) 1000*dasin(lsrotnsc/float(nlnsc)), cbec & 1000*dasin(sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/(nlnsc-1))) cbec write(6,732) 1000*(dasin(rsrotnsc/float(nrnsc)) cbec & - dasin(lsrotnsc/float(nlnsc))), cbec & 1000*sqrt(dasin(sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/ cbec &(nlnsc-1)))**2 cbec & + dasin(sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1)))**2) write(6,726) rsrotnsc/float(nrnsc), & sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1)) write(6,725) lsrotnsc/float(nlnsc), & sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/(nlnsc-1)) write(6,732) rsrotnsc/float(nrnsc) - lsrotnsc/float(nlnsc), & sqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/(nlnsc-1) & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1) ) write(6,*) ' num NO scat det W ', nrnsc , & ' num NO scat det E ', nlnsc write(6,*) if (nls.gt.0) then c write(6,*) rsde, nder, lsde, ndel cbec write(6,706) 1000*dasin(rsrotsc/float(nrs)), cbec & 1000*dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1))) cbec write(6,705) 1000*dasin(lsrotsc/float(nls)), cbec & 1000*dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1))) cbec write(6,712) 1000*(dasin(rsrotsc/float(nrs)) cbec & - dasin(lsrotsc/float(nls))), cbec & 1000*sqrt(dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))**2 cbec & + dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))**2) write(6,706) rsrotsc/float(nrs), & sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)) write(6,705) lsrotsc/float(nls), & sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)) write(6,712) rsrotsc/float(nrs) - lsrotsc/float(nls) , & sqrt( (lsrotscsq-lsrotsc**2/nls)/nls/(nls-1) & + (rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1) ) write(6,*) write(6,702) rsde/float(nder), & sqrt((rsdesq-rsde**2/nder)/nder/(nder-1)) write(6,701) lsde/float(ndel), & sqrt((lsdesq-lsde**2/ndel)/ndel/(ndel-1)) write(6,*) if (nrsl.gt.1) then cbec write(6,708) 1000*dasin(rsrotscle/float(nrsl)), cbec & 1000*dasin(sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1))) write(6,708) rsrotscle/float(nrsl), & sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1)) else write(6,708) 0.,0. endif if(nlsl.gt.1) then cbec write(6,707) 1000*dasin(lsrotscle/float(nlsl)), cbec & 1000*dasin(sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1))) write(6,707) lsrotscle/float(nlsl), & sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1)) else write(6,707) 0.,0. end if write(6,710) rsrotsc/rsrot,rsrotscle/rsrot write(6,709) lsrotsc/lsrot,lsrotscle/lsrot write(6,*) write(6,715) float(nrs), float(nls) write(6,716) float(nrs)/float(nls) write(6,718) float(nrs)-float(nls),Sqrt(float(nls+nrs)) write(6,719) float(nrsl), float(nlsl) else write(6,*) ' NO scattering so limited output written ' endif 700 format(3x,' Rotations for Incoil Ouput to Outcoil input 114.6~cm with & B(mGauss)=',1f9.4) 701 format(3x,' Avg. Energy loss at Det East Dn (meV)',1pe11.4,' +- ' & ,1pe11.4) 702 format(3x,' Avg. Energy loss at Det West Up (meV)',1pe11.4,' +- ' & ,1pe11.4) 703 format(3x,' Avg. Rot East Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 704 format(3x,' Avg. Rot West Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 705 format(3x,' Avg. Rot Scat East Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 706 format(3x,' Avg. Rot Scat West Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 725 format(3x,' Avg. Rot NOScat East Dn (mrad) ',1pe11.4, & ' +- ',1pe11.4) 726 format(3x,' Avg. Rot NOScat West Up (mrad) ',1pe11.4, & ' +- ',1pe11.4) 707 format(3x,' Avg. Rot Scat East Dn E<0.35meV (mrad) ',1pe11.4, & ' +- ',1pe11.4) 708 format(3x,' Avg. Rot Scat West Up E<0.35meV (mrad) ',1pe11.4, & ' +- ',1pe11.4) 709 format(3x,' Frac of rotation from scat East Dn ',1pe11.4, & ', oflow En ',1pe11.4) 710 format(3x,' Frac of rotation from scat West Up ',1pe11.4, & ', oflow En ',1pe11.4) 711 format(3x,' Rotation UpW-DownE (mrad) ',1pe11.4, & ' +- ',1pe11.4) 712 format(3x,' Rotation UpqW-DownE scat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 732 format(3x,' Rotation UpqW-DownE NOScat (mrad) ',1pe11.4, & ' +- ',1pe11.4) 713 format(3x,' # in DetW Up',1pe13.6,2x,', # in Det E Dn ',1pe13.6) 714 format(3x,' Ratio # DetWUp/DetEDn ,'1pe13.6) 717 format(3x,' Det Counts UpW-DownE ',1pe11.4,' +- ',1pe11.4) 715 format(3x,' # scat in DetW Up ',1pe11.4,2x, & ', # scat in Det E Dn ',1pe11.4) 716 format(3x,' Ratio # scat DetWUp/DetEDn ,'1pe11.4) 718 format(3x,' Scat Det Counts UpW-DownE ',1pe11.4,' +- ',1pe11.4) 719 format(3x,'# scat w/ E<0.35meV at det West: ',1pe11.4, & ' East: ',1pe11.4) c write(3,*) ' End time: ', date, time write(3,*) ' HBOOK filename: ',outfname write(3,*) write(3,*) ' Frac # of scatters with Ef.ne.Ei tar1 ', & float(Nenscat1)/float(Ntar1) write(3,*) ' Frac # of scatters with Ef.ne.Ei tar2 ', & float(Nenscat2)/float(Ntar2) write(3,*) ' Frac num of tar1 entering tar2 ', & float(Ntar12)/float(Ntar1) write(3,*) ' Frac of scat of total reaching det ', & float(Nsdet)/float(Nevts) write(3,*) ' Frac of det hits that were scattered ', & float(Nsdet)/float(cot) write(3,*) ' Frac of det hits that were scattered with thsc>0 ', & float(Nthdet)/float(cot) write(3,*) ' Frac of det hits that were scattered with Ef.ne.Ei ', & float(NEdet)/float(cot) write(3,*) ' Frac of det hits that were scattered somewhere ', & float(numendifdet)/float(cot) write(3,713) float(nr), float(nl) write(3,714) float(nr)/float(nl) write(3,717) float(nr)-float(nl),Sqrt(float(nl+nr)) write(3,*) write(3,700) mag_B cbec write(3,704) 1000*dasin(rsrot/float(nr)), cbec & 1000*dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1))) cbec write(3,703) 1000*dasin(lsrot/float(nl)), cbec & 1000*dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1))) cbec write(3,711) 1000*(dasin(rsrot/float(nr))-dasin(lsrot/float(nl))), cbec & 1000*sqrt(dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))**2 cbec & + dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))**2) write(3,704) rsrot/float(nr), & sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)) write(3,703) lsrot/float(nl), & sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)) write(3,711) rsrot/float(nr)-lsrot/float(nl), & sqrt( (lsrotsq-lsrot**2/nl)/nl/(nl-1) & + (rsrotsq-rsrot**2/nr)/nr/(nr-1) ) write(3,*) cbec write(3,726) 1000*dasin(rsrotnsc/float(nrnsc)), cbec & 1000*dasin(sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1))) cbec write(3,725) 1000*dasin(lsrotnsc/float(nlnsc)), cbec & 1000*dasin(sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/(nlnsc-1))) cbec write(3,732) 1000*(dasin(rsrotnsc/float(nrnsc)) cbec & - dasin(lsrotnsc/float(nlnsc))), cbec & 1000*sqrt(dasin(sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/ cbec &(nlnsc-1)))**2 cbec & + dasin(sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1)))**2) write(3,726) rsrotnsc/float(nrnsc), & sqrt((rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1)) write(3,725) lsrotnsc/float(nlnsc), & sqrt((lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/(nlnsc-1)) write(3,732) rsrotnsc/float(nrnsc) - lsrotnsc/float(nlnsc), & sqrt( (lsrotnscsq-lsrotnsc**2/nlnsc)/nlnsc/(nlnsc-1) & + (rsrotnscsq-rsrotnsc**2/nrnsc)/nrnsc/(nrnsc-1) ) write(3,*) write(3,*) ' num NO scat det W ', nrnsc , & ' num NO scat det E ', nlnsc write(3,*) if (nls.gt.0) then c write(3,*) rsde, nder, lsde, ndel cbec write(3,706) 1000*dasin(rsrotsc/float(nrs)), cbec & 1000*dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1))) cbec write(3,705) 1000*dasin(lsrotsc/float(nls)), cbec & 1000*dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1))) cbec write(3,712) 1000*(dasin(rsrotsc/float(nrs)) cbec & - dasin(lsrotsc/float(nls))), cbec & 1000*sqrt(dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))**2 cbec & + dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))**2) write(3,706) rsrotsc/float(nrs), & sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)) write(3,705) lsrotsc/float(nls), & sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)) write(3,712) rsrotsc/float(nrs) - lsrotsc/float(nls) , & sqrt( (lsrotscsq-lsrotsc**2/nls)/nls/(nls-1) & + (rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1) ) write(3,*) write(3,702) rsde/float(nder), & sqrt((rsdesq-rsde**2/nder)/nder/(nder-1)) write(3,701) lsde/float(ndel), & sqrt((lsdesq-lsde**2/ndel)/ndel/(ndel-1)) write(3,*) if (nrsl.gt.1) then cbec write(3,708) 1000*dasin(rsrotscle/float(nrsl)), cbec & 1000*dasin(sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1))) write(3,708) rsrotscle/float(nrsl), & sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1)) else write(3,708) 0.,0. endif if(nlsl.gt.1) then cbec write(3,707) 1000*dasin(lsrotscle/float(nlsl)), cbec & 1000*dasin(sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1))) write(3,707) lsrotscle/float(nlsl), & sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1)) else write(3,707) 0.,0. end if write(3,*) write(3,710) rsrotsc/rsrot,rsrotscle/rsrot write(3,709) lsrotsc/lsrot,lsrotscle/lsrot write(3,715) float(nrs), float(nls) write(3,716) float(nrs)/float(nls) write(3,718) float(nrs)-float(nls),Sqrt(float(nls+nrs)) write(3,719) float(nrsl), float(nlsl) else write(3,*) ' NO scattering so limited output written ' endif c close(unit=3) c bec return end