program main_LUp
c same as main_slim (main_cwn with no 1000s or 2000s histograms) except targets on opposite side
c right up and left dn
c calculating path length and rotation from entire flight path, not just while in target
implicit none
include 'const.inc'
integer loop
double precision ep, wp, kp, k1, thp,php,qlast
integer i, j, nevts,counter,lim,Event,cin,cot,trig !bec
integer Nscat1,Ntar1,Nscat2,Ntar2,Ntar12,NEdet
integer encnt, Nsdet
parameter(nevts=1e9)
real xlam, HRNDM1, vect(5), lm(3), const, xinit, yinit, plr,xtest
character*28 outfname
double precision th0, ph0, thindex
double precision en0, x0, y0, z0, thx0, thy0
double precision en1, x1, y1, z1, thx1, thy1
double precision en2, x2, y2, z2, thx2, thy2
double precision x3, y3, z3, thx3, thy3
double precision x4, y4, z4, thx4, thy4
double precision x5, y5, z5, thx5, thy5
double precision x6, y6, z6, thx6, thy6
double precision x7, y7, z7, thx7, thy7
double precision x8, y8, z8, thx8, thy8
double precision x9, y9, z9, thx9, thy9
double precision x10, y10, z10, thx10, thy10
double precision time1, gl, gin, got
logical scatflg1, scatflg2, tar2flg
integer ranno
integer seedc
integer n1, n2, n3
C parameter(ranno=1)
double precision T
parameter (gl = 111.0, ! Guide length (cm)
+ gin = 0.001095, ! Input critical index (rad/Ang)
+ T = 2.10) ! helium temp. (K)
c + T = 4.25) ! helium temp. (K)
c + got = 0.001) ! Output critical index (rad/Ang)
double precision tx1, tx2, tx3, tx4, ty1, ty2, tz
parameter (tx1 = -3.0, ! cm
+ tx2 = -0.35, ! cm
+ tx3 = 0.35, ! cm
+ tx4 = 3.0, ! cm
+ ty1 = -2.5, ! cm
+ ty2 = 2.5, ! cm
+ tz = 41.6) ! cm
double precision gab1, gab2, gab3, gab4
parameter (gab1 = 10.0, ! gab between incoil and target (cm)
+ gab2 = 8.0, ! cm
+ gab3 = 13.0, gab4 = 10.0) ! cm
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 scd1, scd2, 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
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 = 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)
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 ='xportRUp-'//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 right, x<0. Dn on left '
write(3,*) ' Upstream target on right, x<0. Dn on left '
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','')
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,'q sim(1/Ang)',100,0,4.0,0.) ! bec test q choices
c call hbook1(304,'k beam (1/Ang)',100,0,4.0,0.) ! bec
c call hbook1(404,'k1 (1/Ang)',100,0,4.0,0.) ! bec
call hbook1(1001,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(1,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(2,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(3,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(4,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(5,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(6,'x (cm) Up',100,-3.5,3.5,0.)
call hbook1(1002,'y (cm) Up',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.)
c call hbook1(1005,'E (meV)/ non-scat both',1000,0.,40.,0.)
c call hbook1(1006,'Path len (cm)/ non-scat UpB',100,41.5,41.7,0.)
c call hbook1(1007,'th_rot (mrad)/ non-scat UpB',100,0.,1000.,0.)
c call hbook1(1008,'Path len (cm)/ non-scat UpR',100,41.5,41.7,0.)
c call hbook1(1009,'th_rot (mrad)/ non-scat UpR',100,0.,1000.,0.)
c call hbook1(1010,'Path len (cm)/ non-scat UL',100,41.5,41.7,0.)
c call hbook1(1011,'th_rot (mrad)/ non-scat UL',100,0.,1000.,0.)
c call hbook1(1012,'E (meV)/ scat UpB',100,0.,40.,0.)
c call hbook1(1013,'Path len (cm)/ scat UpB',100,41.,45.,0.)
c call hbook1(1014,'th_rot (mrad)/ scat UpB',100,0.,1000.,0.)
c call hbook1(1015,'Path len (cm)/ scat UpR',100,41.,45.,0.)
c call hbook1(1016,'th_rot (mrad)/ scat UpR',100,0.,1000.,0.)
c call hbook1(1017,'Path len (cm)/ scat UpL',100,41.,45.,0.)
c call hbook1(1018,'th_rot (mrad)/ scat UpL',100,0.,1000.,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.)
c call hbook1(2005,'E (meV)/ non-scat DnB',1000,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.)
c call hbook1(2012,'E (meV)/ scat DnB',100,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) DetB',1000,0.,40.,0.)
call hbook1(4005,'Lambda (Angs) DetB',1000,0.,20.,0.)
call hbook1(3006,'Path length (cm) DeB',200,0.,500.,0.)
call hbook1(3007,'th_rotation (mrad) DetB',200,0.,6000.,0.)
call hbook1(3008,'Path length (cm) DetR',200,0.,500.,0.)
call hbook1(3009,'th_rotation (mrad) DetR',200,0.,6000.,0.)
call hbook1(3010,'Path length (cm) DetL',200,0.,500.,0.)
call hbook1(3011,'th_rotation (mrad) DetL',200,0.,6000.,0.)
call hbook1(3012,'E (meV) DetR',200,0.,40.,0.)
call hbook1(3013,'E (meV) DetL',200,0.,40.,0.)
call hbook1(3014,'E (meV)/ scat DetB',200,0.,40.,0.)
c
c adding more spectra below for diagnostics bec
call hbook1(4010,'sig tot (b/atom)',200,0.,1.,0.)
call hbook1(4011,'interaction len (cm)',200,0.,2000.,0.)
call hbook1(4012,'theta scattered (deg)',200,0.,180.,0.)
call hbook1(4013,'sig Sommers (b/atom)',200,0.,1.,0.)
call hbook1(4014,'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 tar1 B',200,209.,211.,0.)
call hbook1(6002,'th_rotation (mrad) enter tar1 B',200,0.,1500.,0.)
call hbook1(6003,'Path length (cm) Det scat B',200,477.1,477.3,0.)
call hbook1(6004,'th_rotation (mrad) Det scat B',200,0.,6000.,0.)
call hbook1(6005,'Path length (cm) Det scat L',200,477.1,477.3,0.)
call hbook1(6006,'th_rotation (mrad) Det scat L',200,0.,6000.,0.)
call hbook1(6007,'Path length (cm) Det scat R',200,477.1,477.3,0.)
call hbook1(6008,'th_rotation (mrad) Det scat R',200,0.,6000.,0.)
call hbook1(6009,'Path length (cm) Det no scat B',200,477.1,477.3,0.)
call hbook1(6010,'th_rotation (mrad) Det no scat B',200,0.,6000.,0.)
call hbook1(6011,'Path length (cm) Det no scat L',200,477.1,477.3,0.)
call hbook1(6012,'th_rotation (mrad) Det no scat L',200,0.,6000.,0.)
call hbook1(6013,'Path length (cm) Det no scat R',200,477.1,477.3,0.)
call hbook1(6014,'th_rotation (mrad) Det no scat R',200,0.,6000.,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
Nscat1=0 !bec
Ntar1=0 !bec
Nscat2=0 !bec
Ntar2=0 !bec
Nsdet=0 !bec
NEdet=0 !bec
tar2flg = .false. ! bec
c
got = gin
c
do 101 i = 1, nevts
c
Event = 1
pltot = 0.0
thtot = 0.0
c print *, " Envent number ", i
c
c Generate initial neutron wavelength
c
xlam = HRNDM1(100)
c
c The energy calculated from the neutron's wavelength
c
if(xlam.lt.1.) goto 100
en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2
if (en0.le.1.0) encnt =encnt+1
c
c Generates the initial trajectory of the neutron.
c
321 call ranlux(vect,4)
thindex = 0.003*xlam ! cautious number thindex from NG6 closer to 1.7mrad/ang
if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) goto 321
th0 = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec
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 = 2.75d0*(2.*vect(3)-1.)
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 = 2.75*(2.*vect(4)-1.)
z0= 0.0d0
c print *, x0, y0, z0
c
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
k0 = sqrt(2*m_n*en0)/h_par
call hf1(304,real(k0),1.)
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(i.gt.1) call hdelet(302)
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)
loop=loop+1
call wq(q,wp,en0,k0) ! wq modified to prevent E,p consv. probs
ep = en0- (h_par*wp)
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
124 continue
call hf1(1,real(x0),1.)
call hf1(305,real(q),1.)
********************************************************************
********************************************************************
c Neutrons passing thruogh the guide"
*
call guide(gl,gin,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.)
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.*en0/m_n) + thtot
if (pltot.gt.500) write(18,*) 'guide ',pltot
*
********************************************************************
********************************************************************
* Neutrons passing thruogh the input coil
*
call incoil(gin,en0,x1,y1,z1,thx1,thy1,x2,y2,z2,thx2,thy2)
if(x2.lt.-1d5.or.y2.lt.-1d5) goto 100
call hf1(3,real(x2),1.)
pl_del = sqrt((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en0/m_n) + thtot
if (pltot.gt.500) write(18,*) 'incoil ',pltot
*
********************************************************************
********************************************************************
* Neutrons passing thruogh the gab (10.0 cm) between the input
* coil and the target.
*
x3 = x2 + gab1*dtan(thx2)
y3 = y2 + gab1*dtan(thy2)
z3 = z2 + gab1
thx3 = thx2
thy3 = thy2
call hf1(4,real(x3),1.)
call hf1(5017,real(180./3.14159*thx3),1.)
call hf1(5018,real(180./3.14159*thy3),1.)
if(thx3.ge.0..and.thy3.ge.0.) then
phi = datan(dtan(thy3)/dtan(thx3))
elseif(thx3.lt.0..and.thy3.ge.0.) then
phi = pi + datan(dtan(thy3)/dtan(thx3))
elseif(thx3.lt.0..and.thy3.lt.0.) then
phi = pi + datan(dtan(thy3)/dtan(thx3))
elseif(thx3.ge.0..and.thy3.lt.0.) then
phi = 2.*pi + datan(dtan(thy3)/dtan(thx3))
endif
call hf1(5021,real(180./3.14159*phi),1.)
pl_del = sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en0/m_n) + thtot
call hf1(6001,real(pltot),1.)
call hf1(6002,real(thtot),1.)
if (pltot.gt.500) write(18,*) 'gap 1 ',pltot
*
********************************************************************
********************************************************************
c print *, " first half target"
* Neutrons passing thruogh the 1st half of the target (41.6 cm).
*
scatflg1=.false.
if( (x3.ge.tx1.and.x3.le.tx2) ) then
call empty(en0, x3, y3, z3, thx3, thy3, tx1, tx2,
& ty1, ty2, tz, en1, x4, y4, z4, thx4, thy4)
elseif( (x3.ge.tx3.and.x3.le.tx4) ) then
call target(q,T, en0, x3, y3, z3, thx3, thy3, tx3, tx4,
& ty1, ty2, tz, en1, x4, y4, z4, thx4, thy4
& ,sig,intlen,thscat,sigsom,sigsq,scatflg1,thp,php,Ntar1) ! bec
Ntar1=Ntar1+1
if (scatflg1) Nscat1=Nscat1+1 ! count scatters bec
c capture diagnostics for target interaction stuff bec
call hf1(4010,real(sig),1.)
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.)
c write(11,*) 'en0:',en0,' en1:',en1
else
x4 = -1.d6
y4 = -1.d6
z4 = -1.d6
thx4 = -1.d6
thy4 = -1.d6
en1 = -1.d6
endif
*
if(x4.lt.-1d5.or.y4.lt.-1d5) goto 100
*
call hf1(1001,real(x4),1.)
call hf1(1002,real(y4),1.)
c call hf1(1003,real(thx4),1.)
c call hf1(1004,real(thy4),1.)
if(en0.eq.en1) then
pl1_ns = sqrt((x4-x3)**2+(y4-y3)**2+(z4-z3)**2)
th1_ns = gy*mag_B*pl1_ns*0.1/sqrt(2.*en1/m_n)
c write(16,*) th1_ns
pl1_sc = 0.
th1_sc = 0.
pltot = pl1_ns + pltot
thtot = th1_ns + thtot
c call hf1(1005,real(en1),1.)
c call hf1(1006,real(pl1_ns),1.)
c call hf1(1007,real(th1_ns),1.)
if(x3.gt.0..and.x4.gt.0.) then
c call hf1(1008,real(pl1_ns),1.)
c call hf1(1009,real(th1_ns),1.)
elseif(x3.lt.0..and.x4.lt.0.) then
c call hf1(1010,real(pl1_ns),1.)
c call hf1(1011,real(th1_ns),1.)
endif
else
scd1 = (x4-x3-(tz*dtan(thx4)))/(dtan(thx3)-dtan(thx4))
pl1a_sc = sqrt((scd1*dtan(thx3))**2+(scd1*dtan(thy3))**2+
& scd1**2)
pl1b_sc = sqrt( ((tz-scd1)*dtan(thx4))**2
& + ((tz-scd1)*dtan(thy4))**2+(tz-scd1)**2)
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)
c write(17,*) th1a_sc,th1b_sc
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
c call hf1(1012,real(en1),1.)
c call hf1(1013,real(pl1_sc),1.)
c call hf1(1014,real(th1_sc),1.)
if(x3.gt.0..and.x4.gt.0.) then
c call hf1(1015,real(pl1_sc),1.)
c call hf1(1016,real(th1_sc),1.)
elseif(x3.lt.0..and.x4.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
*
********************************************************************
********************************************************************
c print *," Gap btw targets"
* Neutrons passing thruogh the gab (pi-coil = 7.3 cm) between the
* two tagets = 8.0cm.
*
x5 = x4 + gab2*dtan(thx4)
y5 = y4 + gab2*dtan(thy4)
z5 = z4 + gab2
thx5 = thx4
thy5 = thy4
c do NOT allow neutrons from one target cell into opposite side cell
if( (x4.ge.tx1.and.x4.le.tx2).and.(x5.lt.tx1.or.x5.gt.tx2) )
& goto 100
if( (x4.ge.tx3.and.x4.le.tx4).and.(x5.lt.tx3.or.x5.gt.tx4) )
& goto 100
call hf1(5019,real(180./3.14159*thx5),1.)
call hf1(5020,real(180./3.14159*thy5),1.)
if(thx5.ge.0..and.thy5.ge.0.) then
phi = datan(dtan(thy5)/dtan(thx5))
elseif(thx5.lt.0..and.thy5.ge.0.) then
phi = pi + datan(dtan(thy5)/dtan(thx5))
elseif(thx5.lt.0..and.thy5.lt.0.) then
phi = pi + datan(dtan(thy5)/dtan(thx5))
elseif(thx5.ge.0..and.thy5.lt.0.) then
phi = 2.*pi + datan(dtan(thy5)/dtan(thx5))
endif
call hf1(5022,real(180./3.14159*phi),1.)
pl_del = sqrt((x5-x4)**2+(y5-y4)**2+(z5-z4)**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 2 ',pltot
*
********************************************************************
********************************************************************
c print *, " second half target"
* Neutrons passing thruogh the 2nd part of the target (41.6 cm).
*
* since no possibility for scattered neutron from 1st target to get into second target
* use original q value
*
c call hf1(203,real(x5),1.)
if( (x5.gt.tx1.and.x5.lt.tx2) ) then
if (scatflg1) then ! if scat in t1, now scat in both targets
Ntar12=Ntar12+1
call hf1(303,real(q),1.)
endif
scatflg2=.false.
call target(q,T, en1, x5, y5, z5, thx5, thy5, tx1, tx2,
& ty1, ty2, tz, en2, x6, y6, z6, thx6, thy6
& ,sig,intlen,thscat,sigsom,sigsq,scatflg2,thp,php,Ntar2) ! bec
Ntar2=Ntar2+1
if (scatflg2) Nscat2=Nscat2+1
c capture diagnostics for target interaction stuff bec
call hf1(4010,real(sig),1.)
call hf1(4011,real(intlen),1.)
call hf1(5012,real(180./3.14159*thscat),1.)
call hf1(4013,real(sigsom),1.)
call hf1(4014,real(sigsq),1.)
call hf1(5015,real(180./3.14159*thp),1.)
call hf1(5016,real(180./3.14159*php),1.)
c print *, sig, intlen
elseif( (x5.gt.tx3.and.x5.lt.tx4) ) then
call empty(en1, x5, y5, z5, thx5, thy5, tx3, tx4,
& ty1, ty2, tz, en2, x6, y6, z6, thx6, thy6)
else
x6 = -1.d6
y6 = -1.d6
z6 = -1.d6
endif
*
if(x6.lt.-1d5.or.y6.lt.-1d5) goto 100
*
call hf1(2001,real(x6),1.)
call hf1(2002,real(y6),1.)
c call hf1(2003,real(thx6),1.)
c call hf1(2004,real(thy6),1.)
if(en2.eq.en1) then
pl2_ns = sqrt((x6-x5)**2+(y6-y5)**2+(z6-z5)**2)
th2_ns = gy*mag_B*pl2_ns*0.1/sqrt(2.*en2/m_n)
c write(18,*) th2_ns
pl2_sc = 0.
th2_sc = 0.
pltot = pl2_ns + pltot
thtot = th2_ns + thtot
c call hf1(2005,real(en2),1.)
c call hf1(2006,real(pl2_ns),1.)
c call hf1(2007,real(th2_ns),1.)
if(x5.gt.0..and.x6.gt.0.) then
c call hf1(2008,real(pl2_ns),1.)
c call hf1(2009,real(th2_ns),1.)
elseif(x5.lt.0..and.x6.lt.0.) then
c call hf1(2010,real(pl2_ns),1.)
c call hf1(2011,real(th2_ns),1.)
endif
else
scd2 = (x6-x5-(tz*dtan(thx6)))/(dtan(thx5)-dtan(thx6))
pl2a_sc = sqrt((scd2*dtan(thx5))**2+(scd2*dtan(thy5))**2+
& scd2**2)
pl2b_sc = sqrt( ((tz-scd2)*dtan(thx6))**2
& + ((tz-scd2)*dtan(thy6))**2+(tz-scd2)**2)
th2a_sc = gy*mag_B*pl2a_sc*0.1/sqrt(2.*en1/m_n)
th2b_sc = gy*mag_B*pl2b_sc*0.1/sqrt(2.*en2/m_n)
c write(19,*) th2a_sc,th2b_sc
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
c call hf1(2012,real(en2),1.)
c call hf1(2013,real(pl2_sc),1.)
c call hf1(2014,real(th2_sc),1.)
if(x5.gt.0..and.x6.gt.0.) then
c call hf1(2015,real(pl2_sc),1.)
c call hf1(2016,real(th2_sc),1.)
elseif(x5.lt.0..and.x6.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
*
********************************************************************
********************************************************************
* Neutrons passing thruogh the gab (13.0 cm) between the taget and
* the output coil.
*
x7 = x6 + gab3*dtan(thx6)
y7 = y6 + gab3*dtan(thy6)
z7 = z6 + gab3
thx7 = thx6
thy7 = thy6
if(thx7.ge.0..and.thy7.ge.0.) then
phi = datan(dtan(thy7)/dtan(thx7))
elseif(thx7.lt.0..and.thy7.ge.0.) then
phi = pi + datan(dtan(thy7)/dtan(thx7))
elseif(thx7.lt.0..and.thy7.lt.0.) then
phi = pi + datan(dtan(thy7)/dtan(thx7))
elseif(thx7.ge.0..and.thy7.lt.0.) then
phi = 2.*pi + datan(dtan(thy7)/dtan(thx7))
endif
call hf1(5023,real(180./3.14159*phi),1.)
pl_del = sqrt((x7-x6)**2+(y7-y6)**2+(z7-z6)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en2/m_n) + thtot
if (pltot.gt.500) write(18,*) 'gap 3 ',pltot
*
********************************************************************
********************************************************************
* Neutrons passing thruogh the output coil.
*
call outcoil(got,en2,x7,y7,z7,thx7,thy7,x8,y8,z8,thx8,thy8)
if(x8.lt.-1d5.or.y8.lt.-1d5) goto 100
call hf1(5,real(x8),1.)
pl_del = sqrt((x8-x7)**2+(y8-y7)**2+(z8-z7)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en2/m_n) + thtot
if (pltot.gt.500) write(18,*) 'out coil ',pltot
*
********************************************************************
********************************************************************
* Neutrons passing through the gab (10.0 cm) between the output coil
* and the analyzing super mirror (ASM).
*
x9 = x8 + gab4*dtan(thx8)
y9 = y8 + gab4*dtan(thy8)
z9 = z8 + gab4
thx9 = thx8
thy9 = thy8
pl_del = sqrt((x9-x8)**2+(y9-y8)**2+(z9-z8)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en2/m_n) + thtot
if (pltot.gt.500) write(18,*) 'gap 4 ',pltot
********************************************************************
********************************************************************
* Neutrons passing thruogh the supper mirror."
c print *," Supper Mirror"
*
call sm(en2,x9,y9,z9,thx9,thy9,x10,y10,z10,thx10,thy10)
if(x10.lt.-1d5.or.y10.lt.-1d5) goto 100
cot=cot+1
call hf1(6,real(x10),1.)
pl_del = sqrt((x10-x9)**2+(y10-y9)**2+(z10-z9)**2)
pltot = pl_del + pltot
thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en2/m_n) + thtot
if (pltot.gt.500) write(18,*) 'ASM ',pltot
*
********************************************************************
********************************************************************
*
call hf1(3001,real(x10),1.)
call hf1(3002,real(y10),1.)
call hf1(3003,real(thx10),1.)
call hf1(3004,real(thy10),1.)
pl3 = pl1_ns + pl2_ns + pl1_sc + pl2_sc
th3 = th1_ns + th2_ns + th1_sc + th2_sc
c write(20,*) th3
c if(en0.ne.en1 .or. en1.ne.en2)
c & write(*,*) pl1_ns, pl2_ns, pl1_sc, pl2_sc
call hf1(3005,real(en2),1.)
lambda=(2.*pi*h_par)/sqrt(2*m_n*en2)
call hf1(4005,real(lambda),1.)
call hf1(3006,real(pltot),1.)
call hf1(3007,real(thtot),1.)
if(x3.gt.0..and.x6.gt.0.) then
call hf1(3008,real(pltot),1.)
call hf1(3009,real(thtot),1.)
call hf1(3012,real(en2),1.)
elseif(x3.lt.0..and.x6.lt.0.) then
call hf1(3010,real(pltot),1.)
call hf1(3011,real(thtot),1.)
call hf1(3013,real(en2),1.)
endif
c ! misses some since Tsipenyuk cross section is non-zero at q=0 (Escat=E)-- bec
if(en0.ne.en1 .or. en1.ne.en2) then
NEdet=NEdet+1
call hf1(6003, real(pltot),1.)
call hf1(6004, real(thtot),1.)
if (x10.lt.0.) then
call hf1(6005, real(pltot),1.)
call hf1(6006, real(thtot),1.)
else
call hf1(6007, real(pltot),1.)
call hf1(6008, real(thtot),1.)
endif
else
call hf1(6009, real(pltot),1.)
call hf1(6010, real(thtot),1.)
if (x10.lt.0.) then
call hf1(6011, real(pltot),1.)
call hf1(6012, real(thtot),1.)
else
call hf1(6013, real(pltot),1.)
call hf1(6014, real(thtot),1.)
endif
endif
if(scatflg1.or.scatflg2) then ! bec
call hf1(3014,real(en2),1.)
Nsdet=Nsdet+1
endif
if (i.lt.1000) write(17,*) pltot,thtot
********************************************************************
*
34 format(i12,2x,a6)
100 do j=1,10
lim = j*1e7
if(i.eq.lim) write(*,34) i, 'events'
enddo
101 continue
write(*,*) 'input ',nevts
write(*,*) 'tar1 in',Ntar1
write(*,*) 'tar2 in',Ntar2
write(*,*),'det in', cot
write(3,*) 'input ',nevts
write(3,*) 'tar1 in',Ntar1
write(3,*) 'tar2 in',Ntar2
write(3,*),'det in', cot
c print *, " through loop, writing to outfname"
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')
print *, outfname
write(3,*) outfname
c -------------------------
c Adding some printouts at end of run -- bec
Call Date_and_time(date,time,zone,datetime)
print *, " End time: ", date, time
print *, " Frac # of scatters tar1 ", float(Nscat1)/float(Ntar1)
print *, " Frac # of scatters tar2 ", float(Nscat2)/float(Ntar2)
print *, " Fraction of cases where En < 1 meV",
& float(encnt)/float(Nevts)
print *, " Frac num of tar1 entering tar2 ",
& float(Ntar12)/float(Ntar1)
print *, " Number of scattered reaching det ",Nsdet
print *, " Frac of scat of total reaching det ",
& float(Nsdet)/float(Nevts)
print *, " Frac of det hits that were scattered ",
& float(Nsdet)/float(cot)
print *, " Number of det hits with Ef.ne.Ei ",NEdet
print *, " Frac of det hits that were scattered with Ef.ne.Ei ",
& float(NEdet)/float(cot)
c
write(3,*) ' End time: ', date, time
write(3,*) ' Frac # of scatters tar1 ', float(Nscat1)/float(Ntar1)
write(3,*) ' Frac # of scatters tar2 ', float(Nscat2)/float(Ntar2)
write(3,*) ' Fraction of cases where En < 1 meV',
& float(encnt)/float(Nevts)
write(3,*) ' Frac num of tar1 entering tar2 ',
& float(Ntar12)/float(Ntar1)
write(3,*) ' Number of scattered reaching det ',Nsdet
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,*) ' Number of det hits with Ef.ne.Ei ',NEdet
write(3,*) ' Frac of det hits that were scattered with Ef.ne.Ei ',
& float(NEdet)/float(cot)
close(unit=3)
c bec return
end