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 Nenscat1,Ntar1,Nenscat2,Ntar2,Ntar12,NEdet,Nthscat1,Nthscat2 integer encnt, Nsdet,Nthdet parameter(nevts=3e8) 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 tar2flg,enflg1,enflg2,thflg1,thflg2 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 double precision lsrot,lsrotsq, lsrotsc,lsrotscsq double precision lsrotscle,lsrotsclesq,lsde,lsdesq double precision rsrot,rsrotsq, rsrotsc,rsrotscsq double precision rsrotscle,rsrotsclesq,rsde,rsdesq integer nr,nrs,nrsl,nl,nls,nlsl,nder,nders,ndel,ndels 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 call hbook1(306,'w sim(1/psec)',100,0,4.0,0.) ! bec test q choices call hbook1(307,'Eo sim(meV)',3000,0.,40.0,0.) ! bec test q choices call hbook1(308,'Ep sim(meV)',3000,0.,40.0,0.) ! bec test q choices call hbook1(309,'ko sim(1/Ang)',100,0,4.0,0.) ! bec test q choices call hbook1(310,'kp sim(1/Ang)',100,0,4.0,0.) ! bec test q choices 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(20,'E (meV) start B',2000,0.,40.,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.) call hbook1(1005,'E (meV)/ non-scat both',2000,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.) call hbook1(1012,'E (meV)/ scat UpB',2000,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.) call hbook1(2005,'E (meV)/ non-scat DnB',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 DnB',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) DetB',2000,0.,40.,0.) call hbook1(4005,'Lambda (Angs) DetB',1000,0.,50.,0.) call hbook1(3006,'Path length (cm) DeB',200,477.18,477.34,0.) call hbook1(3007,'th_rotation (mrad) DetB',200,-10.,50.,0.) call hbook1(3008,'Path length (cm) DetR',200,477.18,477.34,0.) call hbook1(3009,'th_rotation (mrad) DetR',200,-10.,50.,0.) call hbook1(3010,'Path length (cm) DetL',200,477.18,477.34,0.) call hbook1(3011,'th_rotation (mrad) DetL',200,-10.,50.,0.) call hbook1(3012,'E (meV) DetR',2000,0.,40.,0.) call hbook1(3013,'E (meV) DetL',2000,0.,40.,0.) call hbook1(3014,'E (meV) scat DetB',2000,0.,40.,0.) call hbook1(3020,'E_loss (meV) scat DetB',2000,0.,10.,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,'sig tot (b/atom)',200,0.,1.,0.) c 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.95,210.15,0.) call hbook1(6002,'th_rotation (mrad) enter tar1 B',200,-10.,30.,0.) call hbook1(6003,'Path length (cm) Det scat B',200,477.18,477.34,0.) call hbook1(6004,'th_rotation (mrad) Det scat B',200,-10.,50.,0.) call hbook1(6005,'Path length (cm) Det scat L',200,477.18,477.34,0.) call hbook1(6006,'th_rotation (mrad) Det scat L',200,-10.,50.,0.) call hbook1(6007,'Path length (cm) Det scat R',200,477.18,477.34,0.) call hbook1(6008,'th_rotation (mrad) Det scat R',200,-10.,50.,0.) call hbook1(6009,'Path length (cm) Det no scat B',200,477.18,477.34,0.) call hbook1(6010,'th_rotation (mrad) Det no scat B',200,-10.,50.,0.) call hbook1(6011,'Path length (cm) Det no scat L',200,477.18,477.34,0.) call hbook1(6012,'th_rotation (mrad) Det no scat L',200,-10.,50.,0.) call hbook1(6013,'Path length (cm) Det no scat R',200,477.18,477.34,0.) call hbook1(6014,'th_rotation (mrad) Det no scat R',200,-10.,50.,0.) call hbook2(6020,'th_rot (mrad) vs. E (meV) detB',1000,0.,40., & 200,-10.,50.,0.) c call hbook2(6021,'th rot (mrad) vs. Path len (cm) detB',200,477.18,477.34, c & 200,0.,6000.,0.) c call hbook2(6022,'th rot (mrad) vs. E (meV) detR',1000,0.,40., c & 200,-10.,50.,0.) c call hbook2(6023,'th rot (mrad) vs. Path len (cm) detR',200,477.18,477.34, c & 200,-10.,50.,0.) c call hbook2(6024,'th rot (mrad) vs. E (meV) detL',1000,0.,40., c & 200,-10.,50.,0.) c call hbook2(6025,'th rot (mrad) vs. Path len (cm) detL',200,477.18,477.34, c & 200,-10.,50.,0.) c call hbook2(6026,'E (meV) vs. Path len (cm) detR',200,477.18,477.34, c & 1000,0.,40.,0.) c call hbook2(6027,'E (meV) vs. Path len (cm) detL',200,477.18,477.34, c & 1000,0.,40.,0.) call hbook2(8011,'Thx vs. E (meV) scat DetB',1000,0.,40.,200,0.,0.1,0.) call hbook2(8012,'Thx vs. E-loss (meV) scat DetB',1000,0.,40.,200, & 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 Nthdet=0 lsrot=0;lsrotsq=0;lsrotsc=0;lsrotscsq=0;lsrotscle=0;lsrotsclesq=0; lsde=0;lsdesq=0;rsrot=0;rsrotsq=0;rsrotsc=0;rsrotscsq=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 tar2flg = .false. ! bec c got = gin c do 101 i = 1, nevts c Event = 1 pltot = 0.0 thtot = 0.0 enflg1=.false. thflg1=.false. enflg2=.false. thflg2=.false. 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 c call ranlux(vect,4) ! temporary en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2 c en0 = 50.*vect(1) ! temporary call hf1(20,real(en0),1.) if (en0.le.1.0) encnt =encnt+1 c c Generates the initial trajectory of the neutron. c 321 call ranlux(vect,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 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.) 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 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). 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 * if( (x3.ge.tx1.and.x3.le.tx2) ) then sig = 0.0 ; sigsom = 0.0 ; sigsq =0.0 ; thscat = 0.0 ;thp = 0.0; php = 0.0 call target(q,T, en0, x3, y3, z3, thx3, thy3, tx1, tx2, & ty1, ty2, tz, en1, x4, y4, z4, thx4, thy4 & ,sig,thscat,sigsom,sigsq,enflg1,thflg1,thp,php,Ntar1) ! bec c make sure really went in target -- we know x-position is ok, check y-position if (y3.ge.ty1.and.y3.le.ty2) 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 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 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 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.) c using 1/2 gap length -- rotate to halfway point, reverse, rotate the other half pl_del = 0.5*sqrt((x5-x4)**2+(y5-y4)**2+(z5-z4)**2) c 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 c if (i.lt.100) write(17,*) ' half ', pltot,thtot c reverse thtot in center of gap thtot = - thtot 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 (enflg1.or.thflg1) then ! if scat in t1, now scat in both targets Ntar12=Ntar12+1 c call hf1(303,real(q),1.) endif sig = 0.0 ; sigsom = 0.0 ; sigsq =0.0 ; thscat = 0.0 ;thp = 0.0; php = 0.0 call target(q,T, en0, x5, y5, z5, thx5, thy5, tx3, tx4, & ty1, ty2, tz, en2, x6, y6, z6, thx6, thy6 & ,sig,thscat,sigsom,sigsq,enflg2,thflg2,thp,php,Ntar2) ! bec c make sure really went in target -- we know x-position is ok, check y-position if (y5.ge.ty1.and.y5.le.ty2) 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(4010,real(sig),1.) c 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.) end if c print *, sig, intlen if (en1.ne.en0) write(6,*) ' tar2 en1 NE en0 en1,2,3 x2,3,4' if (en1.ne.en0) write(6,*) en0, en1, en2 if (en1.ne.en0) write(6,*) x2,x3,x4 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 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 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(x10.ge.0.) then call hf1(3008,real(pltot),1.) call hf1(3009,real(thtot),1.) call hf1(3012,real(en2),1.) rsrot = rsrot + dsin(thtot/1000.) rsrotsq = rsrotsq + dsin(thtot/1000.)**2 nr = nr + 1 else call hf1(3010,real(pltot),1.) call hf1(3011,real(thtot),1.) call hf1(3013,real(en2),1.) lsrot = lsrot + dsin(thtot/1000.) lsrotsq = lsrotsq + dsin(thtot/1000.)**2 nl = nl + 1 endif 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.) c call hf2(6024,real(en2),real(thtot),1.) c call hf2(6025,real(pltot),real(thtot),1.) c call hf2(6027,real(pltot),real(en2),1.) lsrotsc = lsrotsc + dsin(thtot/1000.) lsrotscsq = lsrotscsq + dsin(thtot/1000.)**2 nls = nls + 1 if ((en1.le.0.35).or.(en2.le.0.35)) then lsrotscle = lsrotscle + dsin(thtot/1000.) lsrotsclesq = lsrotsclesq + dsin(thtot/1000.)**2 nlsl = nlsl + 1 end if if(en1.ne.en0) then call hf1(3020,real(en0-en1),1.) call hf2(8011,real(en2),real(thx10),1.) call hf2(8012,real(en0-en1),real(thx10),1.) if (thx10.gt.0.1) write(6,*) ' thx, e, eloss', thx10,en2,en0-en1 lsde = lsde + en0-en1 lsdesq = lsdesq + (en0-en1)**2 ndel = ndel + 1 endif if(en2.ne.en1) then call hf1(3020,real(en1-en2),1.) call hf2(8011,real(en2),real(thx10),1.) call hf2(8012,real(en1-en2),real(thx10),1.) if (thx10.gt.0.1) write(6,*) ' thx, e, eloss', thx10,en2,en1-en2 lsde = lsde + en1-en2 lsdesq = lsdesq + (en1-en2)**2 ndel = ndel + 1 endif else call hf1(6007, real(pltot),1.) call hf1(6008, real(thtot),1.) c call hf2(6022,real(en2),real(thtot),1.) c call hf2(6023,real(pltot),real(thtot),1.) c call hf2(6026,real(pltot),real(en2),1.) rsrotsc = rsrotsc + dsin(thtot/1000.) rsrotscsq = rsrotscsq + dsin(thtot/1000.)**2 nrs = nrs + 1 if ((en1.le.0.35).or.(en2.le.0.35)) then rsrotscle = rsrotscle + dsin(thtot/1000.) rsrotsclesq = rsrotsclesq + dsin(thtot/1000.)**2 nrsl = nrsl + 1 end if if(en1.ne.en0) then call hf1(3020,real(en0-en1),1.) call hf2(8011,real(en2),real(thx10),1.) call hf2(8012,real(en0-en1),real(thx10),1.) if (thx10.gt.0.1) write(6,*) ' thx, e, eloss', thx10,en2,en0-en1 rsde = rsde + en0-en1 rsdesq = rsdesq + (en0-en1)**2 nder = nder + 1 endif if(en2.ne.en1) then call hf1(3020,real(en1-en2),1.) call hf2(8011,real(en2),real(thx10),1.) call hf2(8012,real(en1-en2),real(thx10),1.) if (thx10.gt.0.1) write(6,*) ' thx, e, eloss', thx10,en2,en1-en2 rsde = rsde + en1-en2 rsdesq = rsdesq + (en1-en2)**2 nder = nder + 1 endif 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.) c call hf2(6024,real(en2),real(thtot),1.) c call hf2(6025,real(pltot),real(thtot),1.) c call hf2(6027,real(pltot),real(en2),1.) else call hf1(6013, real(pltot),1.) call hf1(6014, real(thtot),1.) c call hf2(6022,real(en2),real(thtot),1.) c call hf2(6023,real(pltot),real(thtot),1.) c call hf2(6026,real(pltot),real(en2),1.) endif endif c if (en2.lt.0.36) then c write(16,888) ' Edet<0.36 -- E0, E1, E2 ', en0, en1, en2 c write(16,889) ' E0-1, E1-2 ', en0-en1, en1-en2 c write(16,888) ' thx, thy', thx10, thy10 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 call hf1(3014,real(en2),1.) Nsdet=Nsdet+1 endif call hf2(6020,real(en2),real(thtot),1.) c call hf2(6021,real(pltot),real(thtot),1.) c if (i.lt.100) write(17,*) ' end ', pltot,thtot ******************************************************************** * 34 format(i12,2x,a6) 100 do j=1,10 lim = j*nevts/10. 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 with Ef.ne.Ei tar1 ", & float(Nenscat1)/float(Ntar1) c print *, " Frac # of scatters with thsc>0 tar1 ", c & float(Nthscat1)/float(Ntar1) print *, " Frac # of scatters with Ef.ne.Ei tar2 ", & float(Nenscat2)/float(Ntar2) c print *, " Frac # of scatters with thsc>0 tar2 ", c & float(Nthscat2)/float(Ntar2) c print *, " Frac num of tar1 entering tar2 ", c & float(Ntar12)/float(Ntar1) c print *, " Number of scattered reaching det ",Nsdet print *, " # of scat reaching det/Nevts ", & float(Nsdet)/float(Nevts) print *, " Frac of det hits that were scattered ", & float(Nsdet)/float(cot) c print *, " Number of det hits with Ef.ne.Ei ",NEdet c print *, " Frac of det hits that were scattered with thsc>0 ", c & float(Nthdet)/float(cot) print *, " Frac of det hits that were scattered with Ef.ne.Ei ", & float(NEdet)/float(cot) write(6,*) write(6,700) mag_B write(6,701) rsde/float(nder),sqrt((rsdesq-rsde**2/nder)/nder/(nder-1)) write(6,702) lsde/float(ndel),sqrt((lsdesq-lsde**2/ndel)/ndel/(ndel-1)) write(6,703) 1000*dasin(rsrot/float(nr)), & 1000*dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1))) write(6,704) 1000*dasin(lsrot/float(nl)), & 1000*dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1))) write(6,705) 1000*dasin(rsrotsc/float(nrs)), & 1000*dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1))) write(6,706) 1000*dasin(lsrotsc/float(nls)), & 1000*dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1))) if (nrsl.gt.1) then write(6,707) 1000*dasin(rsrotscle/float(nrsl)), & 1000*dasin(sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1))) else write(6,707) 0.,0. endif if(nlsl.gt.1) then write(6,708) 1000*dasin(lsrotscle/float(nlsl)), & 1000*dasin(sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1))) else write(6,708) 0.,0. end if write(6,709) rsrotsc/rsrot,rsrotscle/rsrot write(6,710) lsrotsc/lsrot,lsrotscle/lsrot write(6,711) 1000*(dasin(lsrot/float(nl))-dasin(rsrot/float(nr))), & 1000*sqrt(dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))**2 & + dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))**2) write(6,712) 1000*(dasin(lsrotsc/float(nls))-dasin(rsrotsc/float(nrs))), & 1000*sqrt(dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))**2 & + dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))**2) write(6,713) float(nl), float(nr) write(6,714) float(nl)/float(nr) write(6,715) float(nls), float(nrs) write(6,716) float(nls)/float(nrs) 700 format(3x,' Rotations for entire beamline ~477cm with B(mGauss)=',1f9.4) 701 format(3x,' Avg. Energy loss Right ',1pe11.4,' +- ',1pe11.4) 702 format(3x,' Avg. Energy loss Left ',1pe11.4,' +- ',1pe11.4) 703 format(3x,' Avg. Rot Right (mrad) ',1pe11.4,' +- ',1pe11.4) 704 format(3x,' Avg. Rot Left (mrad) ',1pe11.4,' +- ',1pe11.4) 705 format(3x,' Avg. Rot Scat Right (mrad) ',1pe11.4,' +- ',1pe11.4) 706 format(3x,' Avg. Rot Scat Left (mrad) ',1pe11.4,' +- ',1pe11.4) 707 format(3x,' Avg. Rot Scat Right E<0.35meV (mrad) ',1pe11.4,' +- ',1pe11.4) 708 format(3x,' Avg. Rot Scat Left E<0.35meV (mrad) ',1pe11.4,' +- ',1pe11.4) 709 format(3x,' Frac of rotation from scat Right ',1pe11.4, & ', oflow En ',1pe11.4) 710 format(3x,' Frac of rotation from scat Left ',1pe11.4, & ', oflow En ',1pe11.4) 711 format(3x,' Rotation UpL-DownR (mrad) ',1pe11.4,' +- ',1pe11.4) 712 format(3x,' Rotation UpL-DownR scat (mrad) ',1pe11.4,' +- ',1pe11.4) 713 format(3x,' # in DetL ',1pe11.4,2x,', # in Det R ',1pe11.4) 714 format(3x,' Ratio # DetL/DetR ,'1pe11.4) 715 format(3x,' # scat in DetL ',1pe11.4,2x,', # scat in Det R ',1pe11.4) 716 format(3x,' Ratio # scat DetL/DetR ,'1pe11.4) c write(3,*) ' End time: ', date, time write(3,*) ' Frac # of scatters with Ef.ne.Ei tar1 ', & float(Nenscat1)/float(Ntar1) write(3,*) ' Frac # of scatters with th>0 tar1 ', & float(Nthscat1)/float(Ntar2) write(3,*) ' Frac # of scatters with Ef.ne.Ei tar2 ', & float(Nenscat2)/float(Ntar2) write(3,*) ' Frac # of scatters with th>0 tar2 ', & float(Nthscat2)/float(Ntar2) write(3,*) ' Frac num of tar1 entering tar2 ', & float(Ntar12)/float(Ntar1) c 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) c write(3,*) ' Number of det hits with Ef.ne.Ei ',NEdet 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,*) write(3,700) mag_B write(3,701) rsde/float(nder),sqrt((rsdesq-rsde**2/nder)/nder/(nder-1)) write(3,702) lsde/float(nder),sqrt((lsdesq-lsde**2/nder)/nder/(nder-1)) write(3,703) 1000*dasin(rsrot/float(nr)), & 1000*dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1))) write(3,704) 1000*dasin(lsrot/float(nl)), & 1000*dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1))) write(3,705) 1000*dasin(rsrotsc/float(nrs)), & 1000*dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1))) write(3,706) 1000*dasin(lsrotsc/float(nls)), & 1000*dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1))) if (nrsl.gt.1) then write(6,707) 1000*dasin(rsrotscle/float(nrsl)), & 1000*dasin(sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1))) else write(6,707) 0.,0. endif if(nlsl.gt.1) then write(6,708) 1000*dasin(lsrotscle/float(nlsl)), & 1000*dasin(sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1))) else write(6,708) 0.,0. end if write(3,709) rsrotsc/rsrot,rsrotscle/rsrot write(3,710) lsrotsc/lsrot,lsrotscle/lsrot write(3,711) 1000*(dasin(lsrot/float(nl))-dasin(rsrot/float(nr))), & 1000*sqrt(dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))**2 & + dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))**2) write(3,712) 1000*(dasin(lsrotsc/float(nl))-dasin(rsrotsc/float(nr))), & 1000*sqrt(dasin(sqrt((lsrotscsq-lsrotsc**2/nl)/nl/(nl-1)))**2 & + dasin(sqrt((rsrotscsq-rsrotsc**2/nr)/nr/(nr-1)))**2) write(3,713) float(nl), float(nr) write(3,714) float(nl)/float(nr) write(3,715) float(nls), float(nrs) write(3,716) float(nls)/float(nrs) close(unit=3) c bec return end