program main_LUp c same as main_slim (main_cwn with no 1000s or 2000s histograms) except targets on opposite side c left up and right 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=1e8) 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 ='xportLUp-'//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','') 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.,3500.,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.,3500.,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.,3500.,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.,3500.,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.,3500.,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.,3500.,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 target(q,T, en0, x3, y3, z3, thx3, thy3, tx1, tx2, & 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.) elseif( (x3.ge.tx3.and.x3.le.tx4) ) then call empty(en0, x3, y3, z3, thx3, thy3, tx3, tx4, & ty1, ty2, tz, en1, x4, y4, z4, thx4, thy4) 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 call empty(en1, x5, y5, z5, thx5, thy5, tx1, tx2, & ty1, ty2, tz, en2, x6, y6, z6, thx6, thy6) elseif( (x5.gt.tx3.and.x5.lt.tx4) ) 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, tx3, tx4, & 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 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