******************************************************************************* program nuscan1 ******************************************************************************* * * program to scan the SuperK contained event data for signs of neutrino * oscillations. The idea is to employ the contained events with mostly * visible energy and determined directions and to Fourier transform the * oscillation phases. * * first draft: jgl 6 May 1997 * * modified over 20-24 June 97 for L/E plots... JGL * * implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), * senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, + zmin,zmax,ringermin,ringermax * input data for special plots in loeplots and loescan data dm2_test/0.0025/ ! eV^2 best fit for vacuum osc in '04 c data dm2_test/0.0048/ ! eV^2 best fit for vacuum osc in '98 c data dm2_test/0.0135/ ! eV^2 best fit for HPS data s2th_test/1.0/ C* for decaying neutrino fits c data dm2_test/7.0e-5/ ! eV^2 decay gnus c data s2th_test/0.1313/ ! for decaying neutrinos c data dm2_test/2.24e-2/ ! eV^2 decay gnus, sterile c data s2th_test/0.725/ ! for decaying neutrinos, sterile c data dm2_test/2.0e-3/ ! eV^2 decay gnus, second min c data s2th_test/0.780/ ! for decaying neutrinos, second min call get_data call get_cuts call do_data_cuts call get_mcdata call do_mc_cuts call loeplots(dm2_test, s2th_test) call loescan(s2th_test) call excl_plot call asym_plot(dm2_test, s2th_test) stop end ************************************************************************** subroutine loeplots(dm2_test, s2th_test) ************************************************************************** implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), * senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts c common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, c + zmin,zmax,ringermin,ringermax dimension fdmu(10), fdel(10), fmmu(10), fmel(10) dimension ffdmu(10), ffdel(10), ffmmu(10), ffmel(10) dimension fdmuer(10), fdeler(10), fmmuer(10), fmeler(10) dimension ffdmuer(10),ffdeler(10),ffmmuer(10),ffmeler(10) dimension rmu(10), rel(10) dimension frmu(10), frel(10) dimension rmuer(10), reler(10) dimension frmuer(10), freler(10) do ib = 1,10 fdmu(ib) = 0.0 fdel(ib) = 0.0 fmmu(ib) = 0.0 fmel(ib) = 0.0 rmu(ib) = 0.0 rel(ib) = 0.0 ffdmu(ib) = 0.0 ffdel(ib) = 0.0 ffmmu(ib) = 0.0 ffmel(ib) = 0.0 frmu(ib) = 0.0 frel(ib) = 0.0 fdmuer(ib) = 0.0 fdeler(ib) = 0.0 fmmuer(ib) = 0.0 fmeler(ib) = 0.0 rmuer(ib) = 0.0 reler(ib) = 0.0 ffdmuer(ib) = 0.0 ffdeler(ib) = 0.0 ffmmuer(ib) = 0.0 ffmeler(ib) = 0.0 frmuer(ib) = 0.0 freler(ib) = 0.0 end do * make loe histogram for data nused = 0 ndels = 0 ndmus = 0 do ievt = 1, nevts phi = eloed(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. pass(ievt) .gt. 0.0) then nused = nused + 1 alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) if(adpid(ievt) .lt. 0.0) then ndmus = ndmus + 1 fdmu(ialphi) = fdmu(ialphi) + 1. else ndels = ndels + 1 fdel(ialphi) = fdel(ialphi) + 1. end if end if end do * make loe histogram for mc without osc mused = 0 nmels = 0 nmmus = 0 do ievt = 1, mcevts phi = eloem(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. spass(ievt) .gt. 0.0) then mused = mused + 1 alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) if(sadpid(ievt) .lt. 0.0) then nmmus = nmmus + 1 fmmu(ialphi) = fmmu(ialphi) + 1 else nmels = nmels + 1 fmel(ialphi) = fmel(ialphi) + 1 end if end if end do * make loe histogram for fake osc signal in MC feffels = 0.0 feffmus = 0.0 do ievt = 1, mcevts if(abs(sinter(ievt)) .lt. 2.1) then ntyp = 1 ! for electrons cc and nc else ntyp = 2 ! for muons, ditto end if phi = eloem(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. spass(ievt) .gt. 0.0) then alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) rmuel = rfmuel(ievt) phit = eloemt(ievt) call posc(phit,ntyp,dm2_test,s2th_test,rmuel,plep) if(sadpid(ievt) .gt. 0) then ! take it for an electron ffdel(ialphi) = ffdel(ialphi) + plep feffels = feffels + plep else ! an assumed muon ffdmu(ialphi) = ffdmu(ialphi) + plep feffmus = feffmus + plep end if end if end do close(6) open(unit=6, file='loe.prt', type='new') open(unit=7, file='loer.prt', type='new') c write(6,1900) c 1900 format(' output of L/E plots'/) c fnorm = float(nused)/float(mused) fnorm = float(ndels)/float(nmels) c fmorm = fnorm*feff/float(mused) fmorm = fnorm * use this to normalize the mc muons for L/E with oscillations gnorm = float(ndels)/feffels ! data electrons/mc osc els gmorm = gnorm ! and this for electrons feff = feffels + feffmus neff = feff neffels = feffels neffmus = feffmus type *,' data events used = ',nused type *,' mc events used = ',mused type *,' osc mc evts used = ',neff type *,' data muons = ',ndmus type *,' data elecs = ',ndels type *,' mc muons = ',nmmus type *,' mc elecs = ',nmels type *,' mc osc muons = ',neffmus type *,' mc osc elecs = ',neffels type *,' fnorm = ',fnorm,' for norm mc mus' type *,' fmorm = ',fmorm,' for norm mc els' type *,' gnorm = ',gnorm,' for norm osc mc mus' type *,' gmorm = ',gmorm,' for norm osc mc els' * do the ratios and print out elbin = -0.25 do ibin = 1, 10 if(fmel(ibin) .gt. 0.0 .and. fdel(ibin) .gt. 0.0) then rel(ibin) = fdel(ibin)/(fnorm*fmel(ibin)) reler(ibin) = rel(ibin)* * sqrt(1/fdel(ibin)+ 1/fmel(ibin)) else rel(ibin) = 0.0 reler(ibin) = 0.0 if(fdel(ibin) .eq. 0.0) reler(ibin) = 1.0 end if fdeler(ibin)= sqrt(fdel(ibin)) fmeler(ibin)= fnorm*sqrt(fmel(ibin)) ffmel(ibin) = fmel(ibin) fmel(ibin) = fnorm*fmel(ibin) if(fmmu(ibin) .gt. 0.0 .and. fdmu(ibin) .gt. 0.0) then rmu(ibin) = fdmu(ibin)/(fnorm*fmmu(ibin)) rmuer(ibin) = rmu(ibin)* * sqrt(1/fdmu(ibin)+ 1/fmmu(ibin)) else rmu(ibin) = 0.0 rmuer(ibin) = 0.0 if(fdmu(ibin) .eq. 0.0) rmuer(ibin) = 1.0 end if fdmuer(ibin) = sqrt(fdmu(ibin)) fmmuer(ibin) = fnorm*sqrt(fmmu(ibin)) ffmmu(ibin) = fmmu(ibin) fmmu(ibin) = fnorm*fmmu(ibin) relnorm = gmorm/fmorm if(ffmel(ibin) .gt. 0.0 .and. ffdel(ibin) .gt. 0.0) then frel(ibin) = relnorm*ffdel(ibin)/ffmel(ibin) freler(ibin) = frel(ibin)* * sqrt(1/ffdel(ibin)+1/ffmel(ibin)) else frel(ibin) = 0.0 freler(ibin) = 0.0 if(ffdel(ibin) .eq. 0.0) freler(ibin) = 1.0 end if ffdeler(ibin)= gmorm*sqrt(ffdel(ibin)) ffmeler(ibin)= fmorm*sqrt(ffmel(ibin)) ffdel(ibin) = gmorm*ffdel(ibin) ffmel(ibin) = fmorm*ffmel(ibin) rmunorm = gnorm/fnorm if(ffmmu(ibin) .gt. 0.0 .and. ffdmu(ibin) .gt. 0.0) then frmu(ibin) = rmunorm*ffdmu(ibin)/ffmmu(ibin) frmuer(ibin) = frmu(ibin)* * sqrt(1/ffdmu(ibin)+1/ffmmu(ibin)) else frmu(ibin) = 0.0 frmuer(ibin) = 0.0 if(ffdmu(ibin) .eq. 0.0) frmuer(ibin) = 1.0 end if ffdmuer(ibin) = gnorm*sqrt(ffdmu(ibin)) ffmmuer(ibin) = fmorm*sqrt(ffmmu(ibin)) ffdmu(ibin) = gnorm*ffdmu(ibin) ffmmu(ibin) = fmorm*ffmmu(ibin) ib = ibin elbin = elbin + 0.5 write(6,1000) elbin, * fdel(ib), fdmu(ib), fmel(ib), fmmu(ib), * ffdel(ib),ffdmu(ib),ffmel(ib),ffmmu(ib), * rel(ib), rmu(ib), frel(ib), frmu(ib) 1000 format(13(1pe10.2)) elbiner = 0.25 write(7,1000) elbiner, * fdeler(ib), fdmuer(ib), fmeler(ib), fmmuer(ib), * ffdeler(ib),ffdmuer(ib),ffmeler(ib),ffmmuer(ib), * reler(ib), rmuer(ib), freler(ib), frmuer(ib) end do close(6) close(7) return end ************************************************************************** subroutine loescan(s2th_test) ************************************************************************** * * routine to scan through test dm^2 values, testing for chisquare and * printing out list of chisquare values. Ripped off from loeplot routine * and modified appropriately. 6/25/97 jgl. * implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), * senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts c common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, c + zmin,zmax,ringermin,ringermax dimension fdmu(10), fdel(10), fmmu(10), fmel(10) dimension ffdmu(10), ffdel(10), ffmmu(10), ffmel(10) dimension fdmuer(10), fdeler(10), fmmuer(10), fmeler(10) dimension ffdmuer(10),ffdeler(10),ffmmuer(10),ffmeler(10) dimension rmu(10), rel(10) dimension frmu(10), frel(10) dimension rmuer(10), reler(10) dimension frmuer(10), freler(10) dimension Rosc(10), Roscer(10) dimension dm2(100), chisqel(100), chisqmu(100), chisqmup(100) dimension rlikel(100), alrlikel(100) dimension rlikmu(100), alrlikmu(100) dimension rlikmup(100), alrlikmup(100) * changed to agree with other scan ranges, 5/98 data dm2_min/0.00001/ data dm2_max/1.0/ data ndm2/100/ c data s2th_test/1.0/ dm2_step = (dm2_max/dm2_min)**(1./float(ndm2)) close(6) open(unit=6, file='dm2.prt', type='new') do 2000 idm = 1, ndm2 dm2_test = dm2_min*dm2_step**idm dm2(idm) = dm2_test if(idm .eq. 1) then dm2_test = 0.0 dm2(idm) = dm2_test end if chisqel(idm) = 0.0 chisqmu(idm) = 0.0 chisqmup(idm) = 0.0 rlikel(idm) = 1.0 rlikmu(idm) = 1.0 rlikmup(idm) = 1.0 alrlikel(idm) = 0.0 alrlikmu(idm) = 0.0 alrlikmup(idm) = 0.0 Rosc(idm) = 0.0 Roscer(idm) = 0.0 * rezero everything: do ib = 1,10 fdmu(ib) = 0.0 fdel(ib) = 0.0 fmmu(ib) = 0.0 fmel(ib) = 0.0 rmu(ib) = 0.0 rel(ib) = 0.0 ffdmu(ib) = 0.0 ffdel(ib) = 0.0 ffmmu(ib) = 0.0 ffmel(ib) = 0.0 frmu(ib) = 0.0 frel(ib) = 0.0 fdmuer(ib) = 0.0 fdeler(ib) = 0.0 fmmuer(ib) = 0.0 fmeler(ib) = 0.0 rmuer(ib) = 0.0 reler(ib) = 0.0 ffdmuer(ib) = 0.0 ffdeler(ib) = 0.0 ffmmuer(ib) = 0.0 ffmeler(ib) = 0.0 frmuer(ib) = 0.0 freler(ib) = 0.0 end do * make loe histogram for data nused = 0 ndels = 0 ndmus = 0 do ievt = 1, nevts phi = eloed(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. pass(ievt) .gt. 0.0) then nused = nused + 1 alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) if(adpid(ievt) .lt. 0.0) then ndmus = ndmus + 1 fdmu(ialphi) = fdmu(ialphi) + 1. else ndels = ndels + 1 fdel(ialphi) = fdel(ialphi) + 1. end if end if end do * make loe histogram for mc without osc mused = 0 nmels = 0 nmmus = 0 do ievt = 1, mcevts phi = eloem(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. spass(ievt) .gt. 0.0) then mused = mused + 1 alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) if(sadpid(ievt) .lt. 0.0) then nmmus = nmmus + 1 fmmu(ialphi) = fmmu(ialphi) + 1 else nmels = nmels + 1 fmel(ialphi) = fmel(ialphi) + 1 end if end if end do * make loe histogram for fake osc signal in MC feffels = 0.0 feffmus = 0.0 do ievt = 1, mcevts if(abs(sinter(ievt)) .lt. 2.1) then ntyp = 1 ! for electrons cc and nc else ntyp = 2 ! for muons, ditto end if phi = eloem(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. spass(ievt) .gt. 0.0) then alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) * oscillate both electrons and muons appropriately rmuel = rfmuel(ievt) phit = eloemt(ievt) call posc(phit,ntyp,dm2(idm),s2th_test,rmuel,plep) * but bin as observed if(sadpid(ievt) .gt. 0.0) then ! a declared electron ffdel(ialphi) = ffdel(ialphi) + plep feffels = feffels + plep else ffdmu(ialphi) = ffdmu(ialphi) + plep feffmus = feffmus + plep end if end if end do feff = feffels + feffmus * force the muon normalization at each dm2 so it will really test shape. * and electrons get normalized to the muons, so that tests R. fnorm = float(ndmus)/feffmus * use this later for electron normalization test gnorm = float(ndels)/feffels fmorm = fnorm neff = feff c type *,' data events used = ',nused c type *,' mc events used = ',mused c type *,' osc mc evts used = ',neff c type *,' data muons = ',ndmus c type *,' data elecs = ',ndels c type *,' mc muons = ',nmmus c type *,' mc elecs = ',nmels c type *,' fnorm = ',fnorm c type *,' fmorm = ',fmorm * little calculation for R value with oscillations, and its error if( ndels .gt. 0 .and. feffmus .gt. 0 .and. + ndmus .gt. 0 .and. feffels .gt. 0) then Rosc(idm) = (float(ndmus)/float(ndels))/(feffmus/feffels) Roscer(idm) = Rosc(idm)* * sqrt(1./ndmus + 1./ndels + 1./feffmus + 1/feffels) end if * do the ratios and calculate the chisquare between data and oscillated mc elbin = -0.25 do ibin = 1, 10 if(fmel(ibin) .gt. 0.0 .and. fdel(ibin) .gt. 0.0) then rel(ibin) = fdel(ibin)/(fnorm*fmel(ibin)) reler(ibin) = rel(ibin)* * sqrt(1/fdel(ibin)+ 1/fmel(ibin)) else rel(ibin) = 0.0 reler(ibin) = 0.0 if(fdel(ibin) .eq. 0.0) reler(ibin) = 1.0 end if fdeler(ibin)= sqrt(fdel(ibin)) fmeler(ibin)= fnorm*sqrt(fmel(ibin)) ffmel(ibin) = fmel(ibin) fmel(ibin) = fnorm*fmel(ibin) if(fmmu(ibin) .gt. 0.0 .and. fdmu(ibin) .gt. 0.0) then rmu(ibin) = fdmu(ibin)/(fnorm*fmmu(ibin)) rmuer(ibin) = rmu(ibin)* * sqrt(1/fdmu(ibin)+ 1/fmmu(ibin)) else rmu(ibin) = 0.0 rmuer(ibin) = 0.0 if(fdmu(ibin) .eq. 0.0) rmuer(ibin) = 1.0 end if fdmuer(ibin) = sqrt(fdmu(ibin)) fmmuer(ibin) = fnorm*sqrt(fmmu(ibin)) ffmmu(ibin) = fmmu(ibin) fmmu(ibin) = fnorm*fmmu(ibin) if(ffmel(ibin) .gt. 0.0 .and. ffdel(ibin) .gt. 0.0) then frel(ibin) = ffdel(ibin)/ffmel(ibin) freler(ibin) = frel(ibin)* * sqrt(1/ffdel(ibin)+1/ffmel(ibin)) else frel(ibin) = 0.0 freler(ibin) = 0.0 if(ffdel(ibin) .eq. 0.0) freler(ibin) = 1.0 end if ffdeler(ibin)= fmorm*sqrt(ffdel(ibin)) ffmeler(ibin)= fmorm*sqrt(ffmel(ibin)) ffdel(ibin) = fmorm*ffdel(ibin) ffmel(ibin) = fmorm*ffmel(ibin) if(ffmmu(ibin) .gt. 0.0 .and. ffdmu(ibin) .gt. 0.0) then frmu(ibin) = ffdmu(ibin)/ffmmu(ibin) frmuer(ibin) = frmu(ibin)* * sqrt(1/ffdmu(ibin)+1/ffmmu(ibin)) else frmu(ibin) = 0.0 frmuer(ibin) = 0.0 if(ffdmu(ibin) .eq. 0.0) frmuer(ibin) = 1.0 end if ffdmuer(ibin) = fmorm*sqrt(ffdmu(ibin)) ffmmuer(ibin) = fmorm*sqrt(ffmmu(ibin)) ffdmu(ibin) = fmorm*ffdmu(ibin) ffmmu(ibin) = fmorm*ffmmu(ibin) ib = ibin * make the chisquared of the differences between data and mc with osc c fnum = (fdel(ib)+fmorm*ffdel(ib)) fnum = 1.0 + ffdel(ib) if(fnum .gt. 0.0) then chisqel(idm) = chisqel(idm) + * (fdel(ib)-ffdel(ib))**2/fnum end if c fnum = (fdmu(ib)+fmorm*ffdmu(ib)) fnum = 1.0 + ffdmu(ib) if(fnum .gt. 0.0) then !first with mu norm, then el norm chisqmu(idm) = chisqmu(idm) + * (fdmu(ib)-ffdmu(ib))**2/fnum chisqmup(idm) = chisqmup(idm) + * (fdmu(ib)-gnorm*ffdmu(ib)/fnorm)**2/fnum end if data pi/3.1415927/ * do a simple relative liklihood calculation too * first for electrons with muon normalization call poisson(fdel(ib),ffdel(ib),prob) if(prob .le. 1.0e-10) prob = 1.0e-10 rlikel(idm) = rlikel(idm)*prob * this to take care of changing normalization from point to point rlikel(idm) = rlikel(idm)*sqrt(2.0*pi*(1.+ffdel(ib))) * ditto for muons with muon normalization call poisson(fdmu(ib),ffdmu(ib),prob) if(prob .le. 1.0e-10) prob = 1.0e-10 rlikmu(idm) = rlikmu(idm)*prob rlikmu(idm) = rlikmu(idm)*sqrt(2.0*pi*(1.+ffdmu(ib))) * and for muons with electron normalization ffdmup = ffdmu(ib)*gnorm/fnorm call poisson(fdmu(ib),ffdmup,prob) if(prob .le. 1.0e-10) prob = 1.0e-10 rlikmup(idm) = rlikmup(idm)*prob rlikmup(idm) = rlikmup(idm)* * sqrt(2.0*pi*(1.+ffdmu(ib)*gnorm/fnorm)) end do ! end of L/E bin loop c* make per degree of freedom c c chisqel(idm) = chisqel(idm)/10. c chisqmu(idm) = chisqmu(idm)/10. c chisqmup(idm) = chisqmup(idm)/10. if(rlikel(idm).gt.0.0) alrlikel(idm) = -dlog(rlikel(idm)) if(rlikmu(idm).gt.0.0) alrlikmu(idm) = -dlog(rlikmu(idm)) if(rlikmup(idm).gt.0.0) alrlikmup(idm) = -dlog(rlikmup(idm)) if(dm2(idm) .le. 0.0) then aldm2 = 0.0 else aldm2 = dlog10(dm2(idm)) end if c type *, idm,dm2(idm),aldm2,chisqel(idm),chisqmu(idm) write(6,1000) idm,dm2(idm),aldm2,chisqel(idm),chisqmu(idm), * chisqmup(idm),rlikel(idm),rlikmu(idm),rlikmup(idm), * alrlikel(idm),alrlikmu(idm),alrlikmup(idm), * Rosc(idm), Roscer(idm) 1000 format(i10,13(1pe11.2)) 2000 continue ! end of dm2 loop close(6) return end ************************************************************************** subroutine excl_plot ************************************************************************** * * routine to scan through test dm^2 values, testing for chisquare and * printing out list of chisquare values. Ripped off from loeplot routine * and modified appropriately. 6/25/97 jgl. * * further modified, 7/10/97 to add a sweep in s2th to make exclusion plot * implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), *senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts c common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, c + zmin,zmax,ringermin,ringermax dimension fdmu(10), fdel(10), fmmu(10), fmel(10) dimension ffdmu(10), ffdel(10), ffmmu(10), ffmel(10) dimension fdmuer(10), fdeler(10), fmmuer(10), fmeler(10) dimension ffdmuer(10),ffdeler(10),ffmmuer(10),ffmeler(10) dimension rmu(10), rel(10) dimension frmu(10), frel(10) dimension rmuer(10), reler(10) dimension frmuer(10), freler(10) dimension dm2(100), chisqel(100), chisqmu(100), chisqmup(100) dimension rlikel(100), alrlikel(100) dimension rlikmu(100), alrlikmu(100) dimension rlikmup(100), alrlikmup(100) dimension chi2(100,20) * changed to agree with other scans, 5/98 data dm2_min/0.00001/ data dm2_max/1.0/ data ndm2/100/ data ns2th/20/ data s2th_min/0.0/ data s2th_max/1.0/ data chisq_min/100.0/ ! initial min test for chisq min c data s2th_test/1.0/ dm2_step = (dm2_max/dm2_min)**(1./float(ndm2)) s2th_step = (s2th_max - s2th_min)/float(ns2th) close(6) open(unit=6, file='excl.prt', type='new') do 3000 is2 = 1, ns2th s2th_test = s2th_min + (is2 - 0.5)*s2th_step do 2000 idm = 1, ndm2 dm2_test = dm2_min*dm2_step**idm dm2(idm) = dm2_test chisqel(idm) = 0.0 chisqmu(idm) = 0.0 chisqmup(idm) = 0.0 rlikel(idm) = 1.0 rlikmu(idm) = 1.0 rlikmup(idm) = 1.0 alrlikel(idm) = 0.0 alrlikmu(idm) = 0.0 alrlikmup(idm) = 0.0 * rezero everything: do ib = 1,10 fdmu(ib) = 0.0 fdel(ib) = 0.0 fmmu(ib) = 0.0 fmel(ib) = 0.0 rmu(ib) = 0.0 rel(ib) = 0.0 ffdmu(ib) = 0.0 ffdel(ib) = 0.0 ffmmu(ib) = 0.0 ffmel(ib) = 0.0 frmu(ib) = 0.0 frel(ib) = 0.0 fdmuer(ib) = 0.0 fdeler(ib) = 0.0 fmmuer(ib) = 0.0 fmeler(ib) = 0.0 rmuer(ib) = 0.0 reler(ib) = 0.0 ffdmuer(ib) = 0.0 ffdeler(ib) = 0.0 ffmmuer(ib) = 0.0 ffmeler(ib) = 0.0 frmuer(ib) = 0.0 freler(ib) = 0.0 end do * make loe histogram for data nused = 0 ndels = 0 ndmus = 0 do ievt = 1, nevts phi = eloed(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. pass(ievt) .gt. 0.0) then nused = nused + 1 alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) if(adpid(ievt) .lt. 0.0) then ndmus = ndmus + 1 fdmu(ialphi) = fdmu(ialphi) + 1. else ndels = ndels + 1 fdel(ialphi) = fdel(ialphi) + 1. end if end if end do * make loe histogram for mc without osc mused = 0 nmels = 0 nmmus = 0 do ievt = 1, mcevts phi = eloem(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. spass(ievt) .gt. 0.0) then mused = mused + 1 alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) if(sadpid(ievt) .lt. 0.0) then nmmus = nmmus + 1 fmmu(ialphi) = fmmu(ialphi) + 1 else nmels = nmels + 1 fmel(ialphi) = fmel(ialphi) + 1 end if end if end do * make loe histogram for fake osc signal in MC feffels = 0.0 feffmus = 0.0 do ievt = 1, mcevts if(abs(sinter(ievt)) .lt. 2.1) then ntyp = 1 ! for electrons cc and nc else ntyp = 2 ! for muons, ditto end if phi = eloem(ievt) if(phi .gt. 1.0 .and. phi .lt. 1.0e5 + .and. spass(ievt) .gt. 0.0) then alphi = dlog10(phi) ialphi = 1 + int(alphi*2.0) ! in bins of 1/2 in log10(L/E) * oscillate both electrons and muons appropriately rmuel = rfmuel(ievt) phit = eloemt(ievt) call posc(phit,ntyp,dm2(idm),s2th_test,rmuel,plep) * but bin as observed if(sadpid(ievt) .gt. 0.0) then ffdel(ialphi) = ffdel(ialphi) + plep feffels = feffels + plep else ffdmu(ialphi) = ffdmu(ialphi) + plep feffmus = feffmus + plep end if end if end do feff = feffels + feffmus * force the muon normalization at each dm2 so it will really test shape. * and elctrons get normalized to the muons, so that tests R. fnorm = float(ndmus)/feffmus * use this later for electron normalization test gnorm = float(ndels)/feffels fmorm = fnorm neff = feff c type *,' data events used = ',nused c type *,' mc events used = ',mused c type *,' osc mc evts used = ',neff c type *,' data muons = ',ndmus c type *,' data elecs = ',ndels c type *,' mc muons = ',nmmus c type *,' mc elecs = ',nmels c type *,' fnorm = ',fnorm c type *,' fmorm = ',fmorm * do the ratios and calculate the chisquare between data and oscillated mc elbin = -0.25 do ibin = 1, 10 if(fmel(ibin) .gt. 0.0 .and. fdel(ibin) .gt. 0.0) then rel(ibin) = fdel(ibin)/(fnorm*fmel(ibin)) reler(ibin) = rel(ibin)* * sqrt(1/fdel(ibin)+ 1/fmel(ibin)) else rel(ibin) = 0.0 reler(ibin) = 0.0 if(fdel(ibin) .eq. 0.0) reler(ibin) = 1.0 end if fdeler(ibin)= sqrt(fdel(ibin)) fmeler(ibin)= fnorm*sqrt(fmel(ibin)) ffmel(ibin) = fmel(ibin) fmel(ibin) = fnorm*fmel(ibin) if(fmmu(ibin) .gt. 0.0 .and. fdmu(ibin) .gt. 0.0) then rmu(ibin) = fdmu(ibin)/(fnorm*fmmu(ibin)) rmuer(ibin) = rmu(ibin)* * sqrt(1/fdmu(ibin)+ 1/fmmu(ibin)) else rmu(ibin) = 0.0 rmuer(ibin) = 0.0 if(fdmu(ibin) .eq. 0.0) rmuer(ibin) = 1.0 end if fdmuer(ibin) = sqrt(fdmu(ibin)) fmmuer(ibin) = fnorm*sqrt(fmmu(ibin)) ffmmu(ibin) = fmmu(ibin) fmmu(ibin) = fnorm*fmmu(ibin) if(ffmel(ibin) .gt. 0.0 .and. ffdel(ibin) .gt. 0.0) then frel(ibin) = ffdel(ibin)/ffmel(ibin) freler(ibin) = frel(ibin)* * sqrt(1/ffdel(ibin)+1/ffmel(ibin)) else frel(ibin) = 0.0 freler(ibin) = 0.0 if(ffdel(ibin) .eq. 0.0) freler(ibin) = 1.0 end if ffdeler(ibin)= fmorm*sqrt(ffdel(ibin)) ffmeler(ibin)= fmorm*sqrt(ffmel(ibin)) ffdel(ibin) = fmorm*ffdel(ibin) ffmel(ibin) = fmorm*ffmel(ibin) if(ffmmu(ibin) .gt. 0.0 .and. ffdmu(ibin) .gt. 0.0) then frmu(ibin) = ffdmu(ibin)/ffmmu(ibin) frmuer(ibin) = frmu(ibin)* * sqrt(1/ffdmu(ibin)+1/ffmmu(ibin)) else frmu(ibin) = 0.0 frmuer(ibin) = 0.0 if(ffdmu(ibin) .eq. 0.0) frmuer(ibin) = 1.0 end if ffdmuer(ibin) = fmorm*sqrt(ffdmu(ibin)) ffmmuer(ibin) = fmorm*sqrt(ffmmu(ibin)) ffdmu(ibin) = fmorm*ffdmu(ibin) ffmmu(ibin) = fmorm*ffmmu(ibin) ib = ibin * make the chisquared of the differences between data and mc with osc c fnum = (fdel(ib)+fmorm*ffdel(ib)) fnum = 1.0 + ffdel(ib) if(fnum .gt. 0.0) then chisqel(idm) = chisqel(idm) + * (fdel(ib)-ffdel(ib))**2/fnum end if c fnum = (fdmu(ib)+fmorm*ffdmu(ib)) fnum = 1.0 + ffdmu(ib) if(fnum .gt. 0.0) then !first with mu norm, then el norm chisqmu(idm) = chisqmu(idm) + * (fdmu(ib)-ffdmu(ib))**2/fnum chisqmup(idm) = chisqmup(idm) + * (fdmu(ib)-gnorm*ffdmu(ib)/fnorm)**2/fnum end if data pi/3.1415927/ * do a simple relative liklihood calculation too * first for electrons with muon normalization call poisson(fdel(ib),ffdel(ib),prob) if(prob .le. 1.0e-10) prob = 1.0e-10 rlikel(idm) = rlikel(idm)*prob * this to take care of changing normalization from point to point rlikel(idm) = rlikel(idm)*sqrt(2.0*pi*(1.+ffdel(ib))) * ditto for muons with muon normalization call poisson(fdmu(ib),ffdmu(ib),prob) if(prob .le. 1.0e-10) prob = 1.0e-10 rlikmu(idm) = rlikmu(idm)*prob rlikmu(idm) = rlikmu(idm)*sqrt(2.0*pi*(1.+ffdmu(ib))) * and for muons with electron normalization ffdmup = ffdmu(ib)*gnorm/fnorm call poisson(fdmu(ib),ffdmup,prob) if(prob .le. 1.0e-10) prob = 1.0e-10 rlikmup(idm) = rlikmup(idm)*prob rlikmup(idm) = rlikmup(idm)* * sqrt(2.0*pi*(1.+ffdmu(ib)*gnorm/fnorm)) end do ! end of L/E bin loop c* make per degree of freedom c c chisqel(idm) = chisqel(idm)/10. c chisqmu(idm) = chisqmu(idm)/10. c chisqmup(idm) = chisqmup(idm)/10. * load the matrix of chi2s chi2(idm,is2) = chisqmup(idm) if(rlikel(idm).gt.0.0) alrlikel(idm) = -dlog(rlikel(idm)) if(rlikmu(idm).gt.0.0) alrlikmu(idm) = -dlog(rlikmu(idm)) if(rlikmup(idm).gt.0.0) alrlikmup(idm) = -dlog(rlikmup(idm)) aldm2 = dlog10(dm2(idm)) c type *, idm,dm2(idm),aldm2,chisqel(idm),chisqmu(idm) if(chisqmup(idm) .le. chisq_min) then ! save new value chisq_min = chisqmup(idm) s2th_minchisq = s2th_test dm2_minchisq = dm2_test end if 2000 continue ! end of dm2 loop 3000 continue ! end of sin^2(2 theta) loop type *,' exclusion sweep complete' type *,' min chisquare of ',chisq_min,' found at' type *,' sin^2(2th) = ',s2th_minchisq type *,' dm2 = ',dm2_minchisq,' eV^2' type *,' ' * now print out the matrix: do idm = 1,100 write(6,1000) (chi2(idm,is2), is2 = 1,20) 1000 format(20(1pe11.2)) end do close(6) return end ****************************************************************************** subroutine poisson(enn,emm,prob) ****************************************************************************** * * poisson prob of getting enn when given expectation emm * implicit real*8 (a-h,o-z) data pi/3.14159265/,nctmax/30/,ict/0/,elpmax/80.0/ * elpmax is the maximum exp value that the machine will handle. * updated for q format, 9/90 by jgl c type *,' poisson called with en = ',enn,' and emm = ',emm en = enn em = emm prob = 0.0 * if mean is zero return zero probability if( em .eq. 0.0 ) return * if negative values bomb out if( en .lt. 0 .or. em .lt. 0.0 ) go to 2000 * special case of n=0, do just exp(-m) if( en .eq. 0.0 ) go to 1500 * for n > 30 we can use sterling's approximation if( en .gt. 30 ) go to 1000 * thus in this branch 1 <= n <= 30 n = en fact = 1.0 do 100 i = 1, n fact = fact * i 100 continue if(em .gt. 80.0) go to 700 if(em .lt. 1.0e-6) go to 700 expm = exp(-em) prob = ( (em**en) /fact ) * expm return * for small values of poisson use exponential form 700 continue elp = en * dlog(em) -em -dlog(fact) if( elp .lt. -elpmax ) return prob = 1.0 if( elp .gt. 0.0 ) return prob = exp( elp ) return * here for values of n big enough to use stirling's formula 1000 continue elp = en * dlog(em/en) -em + en -(1.0/(12.0*en)) if( elp .lt. -elpmax ) return prob = 1.0 / sqrt( 2.0 * pi * en ) if( elp .gt. 0.0 ) return prob = prob * exp(elp) return * here for special case of n=0 1500 continue if( em .gt. elpmax ) return prob = exp(-em) return * bomb out to here, return zero with warning first nxtmax times 2000 continue nct = nct + 1 if( nct .gt. nctmax ) return write(6,2100) en, em 2100 format(' poisson called with n = ',1pe15.5, + ' and m = ',1pe15.5) return end ******************************************************************************* subroutine get_data ******************************************************************************* implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts data rhodet/1690.0/ ! nominal detector radius in x-y plane, cm data zdet/1810.0/ ! nominal detector height above and below ref plane open(unit=5, file='he1489_1.dat', type='old') ievt = 0 100 continue ievt = ievt + 1 c old c read(5,*,err=1000)irun,ievent,x(ievt),y(ievt),z(ievt), c * dx(ievt),dy(ievt),dz(ievt),evis(ievt),adpid(ievt) read(5,*,err=1000)irun,ievent,x(ievt),y(ievt),z(ievt), * dx(ievt),dy(ievt),dz(ievt),evis(ievt),ip c 200 format(2i,8e) c 200 format(2i,8e14.6) * ip(1) is on-site particle id for track one, 1=gam, 2=el, 3=mu * adpid should be .lt. 0 for mus adpid(ievt) = 2.5 - ip if(ievt .lt. 100) type *,irun,ievent,dz(ievt),evis(ievt),ip run(ievt) = irun event(ievt) = ievent * add energy corrections to evis for muons and electrons * per John Flanagan's plots of 6/28/97 * from Dave Casper 7/7/97: * = 1.00*Evis + 487 (for identified showers) and * = 1.10*Evis + 360 (for identified muons). if(adpid(ievt) .lt. 0.0) then c evis(ievt) = evis(ievt) + 530 ! JF's suggestion c evis(ievt) = evis(ievt) + 158. c evis(ievt) = evis(ievt) + 258. enuest(ievt) = 1.1*evis(ievt) + 360. ! Casper ntyp = 2 else c evis(ievt) = evis(ievt) + 259 ! JF c evis(ievt) = evis(ievt) + 0.0 c evis(ievt) = evis(ievt) + 100. enuest(ievt) = 1.0*evis(ievt) + 487. ! Casper ntyp = 1 end if * make estimate of L/E call elp(dz(ievt),ntyp,ell) ell = ell * 1000.0 ! km egev = enuest(ievt)/1000.0 ! GeV eloed(ievt) = 0.0 if(egev .gt. 0.0) eloed(ievt) = ell/egev ! km/GeV * calc rwall rhoevent = sqrt(x(ievt)**2 + y(ievt)**2) delrho = rhodet - rhoevent delz = zdet - dabs(z(ievt)) if( delrho .lt. delz ) then rwall(ievt) = delrho else rwall(ievt) = delz end if go to 100 1000 continue nevts = ievt - 1 type *,' got ',nevts,' from first data file' close(5) open(unit=5, file='he1489_2.dat', type='old') ievt = 0 250 continue ievt = ievt + 1 c read(5,*,err=2000) anisx(ievt),anisy(ievt),anisz(ievt), c * iyrmthday,ihrminsec,ringer(ievt),dinter(ievt) read(5,*,err=2000) anisx(ievt),anisy(ievt),anisz(ievt), * iyr,imth,iday,ihr,imin,isec,nring c 300 format(3e,2i,2e) c 300 format(3e14.6,2i,2e14.6) ringer(ievt) = nring yrmthday(ievt) = iyr*10000.0 +imth*100 +iday hrminsec(ievt) = ihr*10000.0 +imin*100 +isec if(ievt .lt. 100) type *,' ievt =',ievt,' ringer =',ringer(ievt) go to 250 2000 continue mevts = ievt -1 close(5) type *,' got ',nevts,' data events' if(nevts .ne. mevts) then type *,' different numbers of events from he1489_1.dat' type *,' and from he1489_2.dat',nevts,mevts stop end if return end ******************************************************************************* subroutine get_mcdata ******************************************************************************* implicit real*8 (a-h,o-z) parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), * senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts data rhodet/1690.0/ ! nominal detector radius in x-y plane, cm data zdet/1810.0/ ! nominal detector height above and below ref plane open(unit=5, file='hemc_1.dat', type='old') ievt = 0 100 continue ievt = ievt + 1 c read(5,*,err=1000)isrun,isevent, c * sx(ievt),sy(ievt),sz(ievt), c * sdx(ievt),sdy(ievt),sdz(ievt),sevis(ievt),sadpid(ievt) read(5,*,err=1000)isrun,isevent, * sx(ievt),sy(ievt),sz(ievt), * sdx(ievt),sdy(ievt),sdz(ievt),sevis(ievt),ip * ip(1) is on-site particle id for track one, 1=gam, 2=el, 3=mu * sadpid should be .lt. 0 for mus sadpid(ievt) = 2.5 - ip ***>>> kluge because of error in fcmc adir(3) ??? sdz(ievt) = -sdz(ievt) * senergy = neutrino energy * sdzn = neutrino zenith angle c 200 format(2i,8e) c 200 format(2i,8e14.6) srun(ievt) = isrun sevent(ievt) = isevent if(ievt .lt. 100) type *,ievt,sdz(ievt),sevis(ievt) * add unseen energy to evis for muons if(sadpid(ievt) .lt. 0.0) then c sevis(ievt) = sevis(ievt) + 158. c sevis(ievt) = sevis(ievt) + 530. ! JohnF's offsets c sevis(ievt) = sevis(ievt) + 258. senuest(ievt) = 1.1*sevis(ievt) + 360. ! Casper ntyp = 2 else c sevis(ievt) = sevis(ievt) + 0.0 c sevis(ievt) = sevis(ievt) + 259. ! JohnF's offsets c sevis(ievt) = sevis(ievt) + 100. senuest(ievt) = 1.0*sevis(ievt) + 487. ! Casper ntyp = 1 end if * make estimated L/E call elp(sdz(ievt),ntyp,ell) ell = ell * 1000.0 ! km egev = senuest(ievt)/1000.0 ! GeV eloem(ievt) = 0.0 if(egev .gt. 0.0) eloem(ievt) = ell/egev ! km/GeV * calc srwall rhoevent = sqrt(sx(ievt)**2 + sy(ievt)**2) delrho = rhodet - rhoevent delz = zdet - dabs(sz(ievt)) if( delrho .lt. delz ) then srwall(ievt) = delrho else srwall(ievt) = delz end if go to 100 1000 continue mcevts = ievt - 1 type *,' got first half of MC evts, total = ',mcevts ievt= 0 close(5) open(unit=5, file='hemc_2.dat', type='old') 250 continue ievt= ievt + 1 c read(5,*,err=2000) sanisx(ievt),sanisy(ievt),sanisz(ievt), c * senergy(ievt),sdzn(ievt), sringer(ievt), isinter read(5,*,err=2000) sanisx(ievt),sanisy(ievt),sanisz(ievt), * pnu, sdzn(ievt), nring, ipnu c type *,' for second MC event, ievt =',ievt c 300 format(7e) c 300 format(7e14.6) sringer(ievt) = nring * shift the interacting particle code if(ipnu .eq. 12) sinter(ievt) = 1 ! nu_e if(ipnu .eq. -12) sinter(ievt) = 2 ! nu_ebar if(ipnu .eq. 14) sinter(ievt) = 3 ! nu_mu if(ipnu .eq. -14) sinter(ievt) = 4 ! nu_mubar if(ipnu .eq. 16) sinter(ievt) = 5 ! nu_tau if(ipnu .eq. -16) sinter(ievt) = 6 ! nu_taubar * now do true quantities for mc data, for oscillating events if(abs(sinter(ievt)) .lt. 2.1) then ntyp = 1 ! electron, cc or nc, gnu or gnubar else ntyp = 2 ! muon end if call elp(sdzn(ievt),ntyp,ell) ell = ell * 1000.0 ! km senergy(ievt) = pnu*1000.0 ! MeV segev = senergy(ievt)/1000.0 ! GeV eloemt(ievt) = 0.0 if(segev .gt. 0.0) eloemt(ievt) = ell/segev ! km/GeV costheta = -sdzn(ievt) call Fmuel(senergy(ievt), costheta, rmuel) rfmuel(ievt) = rmuel go to 250 2000 continue ncevts = ievt - 1 type *,' got ',mcevts,' monte carlo events' if(ncevts .ne. mcevts) then type *,' numbers of events different in hemc_1.dat and' type *,' hemc_2.dat',ncevts, mcevts stop end if return end ******************************************************************************* subroutine get_cuts ******************************************************************************* implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, + zmin,zmax,ringermin,ringermax * default cuts * muons data adpidmin/-1.0e6/ c data adpidmax/0.0/ * electrons c data adpidmin/0.0/ data adpidmax/+1.0e6/ data evismin/200./ ! MeV c data evismin/400./ ! MeV c data evismin/600./ ! MeV c data evismin/1200./ ! MeV data evismax/10000./ ! MeV * rwall cut data rwallmin/200./ !cm distance of vertex from nearest wall c data rwallmin/500./ !cm distance of vertex from nearest wall * zcuts for making top and bottom scans data zmin/-1810.0/ ! cm, dist from central ref plane c data zmax/0.0/ ! cm, dist from central ref plane data zmax/+1810.0/ * ringer cuts here data ringermin/0.5/ ! min for single ring c data ringermax/14.0/ ! max for single ring data ringermax/1.5/ ! max for single ring ****>>> put in on-line controls here type *,' ' type *,' CUTS USED IN THIS RUN: ' type *,' ' type *,' adpid < 0 for non-showering tracks' type *,' adpid min = ',adpidmin type *,' adpid max = ',adpidmax type *,' the following is for corrected energy:' type *,' evis min = ',evismin,' MeV' type *,' evis max = ',evismax,' MeV' type *,' ' type *,' rwall min = ',rwallmin,' cm' type *,' ' type *,' z vertex min = ',zmin,' cm' type *,' z vertex max = ',zmax,' cm' type *,' ' type *,' ringer min = ',ringermin type *,' ringer max = ',ringermax type *,' ' return end ******************************************************************************* subroutine do_data_cuts ******************************************************************************* implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, + zmin,zmax,ringermin,ringermax * make cuts npass = 0 npid_fail = 0 ne_fail = 0 nrwall_fail = 0 nz_fail = 0 nring_fail= 0 do ievt = 1,nevts pass(ievt) = +1.0 ! pass * test for failure on any criterium if(adpid(ievt) .lt. adpidmin .or. + adpid(ievt) .gt. adpidmax) then pass(ievt) = -1.0 npid_fail = npid_fail + 1 end if if(evis(ievt) .lt. evismin .or. + evis(ievt) .gt. evismax) then pass(ievt) = -1.0 ne_fail = ne_fail + 1 end if c take this out as it is done in datadump anyway, and something is wrong too? c if(rwall(ievt) .lt. rwallmin) then c pass(ievt) = -1.0 c nrwall_fail = nrwall_fail + 1 c end if c might as well take out z cut too while at it. c if(z(ievt) .lt. zmin .or. c + z(ievt) .gt. zmax) then c pass(ievt) = -1.0 c nz_fail = nz_fail + 1 c end if if(ringer(ievt) .lt. ringermin .or. + ringer(ievt) .gt. ringermax) then pass(ievt) = -1.0 nring_fail = nring_fail + 1 end if * count passes if(pass(ievt) .gt. 0.0) npass = npass + 1 end do nfail = npid_fail + ne_fail + nrwall_fail + nz_fail + nring_fail type *,' data events passing cuts = ',npass type *,' failing cuts = ',nfail type *,' fail pid = ',npid_fail type *,' fail evis = ',ne_fail type *,' fail rwall = ',nrwall_fail type *,' fail z = ',nz_fail type *,' fail ring = ',nring_fail return end ******************************************************************************* subroutine do_mc_cuts ******************************************************************************* implicit real*8 (a-h,o-z) parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), * senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, + zmin,zmax,ringermin,ringermax * make cuts npass = 0 do ievt = 1, mcevts spass(ievt) = +1.0 ! pass * test for failure on any criterium if(sadpid(ievt) .lt. adpidmin) spass(ievt) = -1.0 if(sadpid(ievt) .gt. adpidmax) spass(ievt) = -1.0 if(sevis(ievt) .lt. evismin) spass(ievt) = -1.0 if(sevis(ievt) .gt. evismax) spass(ievt) = -1.0 if(srwall(ievt) .lt. rwallmin) spass(ievt) = -1.0 if(sz(ievt) .lt. zmin) spass(ievt) = -1.0 if(sz(ievt) .gt. zmax) spass(ievt) = -1.0 if(sringer(ievt) .lt. ringermin) spass(ievt) = -1.0 if(sringer(ievt) .gt. ringermax) spass(ievt) = -1.0 * count passes if(spass(ievt) .gt. 0.0) npass = npass + 1 end do type *,' mc events passing cuts = ',npass return end ************************************************************************** function el(dz) ************************************************************************** * * for given vertical direction cosine dz, returns source flight distance * in 1000km units * implicit real*8 (a-h,o-z) logical firstcall data firstcall/.true./ data Re/6.4/ ! earth radius in 1000's of km data D/0.02/ ! depth below production for neutrinos, in 10^3 km if(firstcall) then Rm = Re-D H = D*(2.0*Re-D) ! the horizontal distance to production firstcall = .false. end if * dz is zenith direction of motion cosine (- is downgoing, opposite of CR * convention) Rmdz = Rm*dz Rmdz2 = Rmdz**2 el = sqrt( Rmdz2 + H ) + Rmdz return end ************************************************************************** subroutine elp(dz,itype,ell) ************************************************************************** * * for given vertical direction cosine dz, returns source flight distance * in 1000km units * * call with dz = dir cos, -1 dn; itype = 1 for electrons, 2 for muons * returns distance ell * * this version has an approximation to the Stanev MC generated distances * hacked by jgl, 28 June 1997 * implicit real*8 (a-h,o-z) data Rearth/6378.0/ ! km * constant is thickness of atmosphere/some absorption length data const/10.0/ * roughly the scale height of atmosphere data H0/7.0/ ! km, if(itype .eq. 1) then Depth = 14.0 ! km for electrons else Depth = 16.0 ! km for muons end if Reff = Rearth + Depth ! km * get initial guess at local zenith angle, theta_star costheta = -dz sintheta = sqrt(1.0-costheta**2) sintheta_star = sintheta*Rearth/Reff if(sintheta_star .lt. 1.0) then costheta_star = sqrt(1.0-sintheta_star**2) else costheta_star = 0.0 end if * try simple approximation to changing height: if(itype .eq. 1) then H = H0 * 0.875 else H = H0 end if tDepth = H*dlog(const/costheta_star) Reff = Rearth + tDepth * iterate once on the theta* and tDepth sintheta_star = sintheta*Rearth/Reff if(sintheta_star .lt. 1.0) then costheta_star = sqrt(1.0-sintheta_star**2) else costheta_star = 0.0 end if * note that costheta_star is always positive (all particles heading into earth) tDepth = H*dlog(const/costheta_star) Re = Reff/1000.0 D = tDepth/1000.0 Rm = Rearth/1000.0 Hor2 = D*(Re + Rm) ! the horizontal distance to production, sqrd * dz is zenith direction of motion cosine (- is downgoing, opposite of CR * convention) Rmdz = Rm*dz Rmdz2 = Rmdz**2 ell = sqrt( Rmdz2 + Hor2 ) + Rmdz return end **************************************************************************** subroutine posc(eloe,ntyp,dm2,s2th,rmuel,prob) **************************************************************************** * * calculates oscillation probabilities for numu and nue for various * scenarios * call with eloe = L/E in km/GeV, * ntyp = 1 for electrons, 2 for muons * dm2 mass square diff in eV^2 * rmuel = nu_mu/nu_e flux ratio at that energy and angle (calculated) * and the routine returns the prob of a muon being a muon and same * for electron... can be greater than 1 if being fed.... * implicit real*8 (a-h,o-z) logical firstcall data firstcall/.true./ data cmptkm/1.0e5/ ! cm/km data eVpGeV/1.0e+9/ ! eV/GeV data hbarc/1.97e-5/ ! eV cm data isol/1/ ! 1 for solar min,2 for solar max fluxes c data scalefactor/1.0/ ! scale factor for using dm2 as alpha if(firstcall) then cc = cmptkm/(4.0*eVpGeV*hbarc) firstcall = .false. s1 = 8.0/9.0 s2 = 4.0/9.0 end if phi = cc*dm2*eloe ! radians ppp = (sin(phi))**2 c for decaying neutrinos use c alpha = dm2*scalefactor ! puts alpha in range so that dm2=1e-3 c ! <-> alpha = 1/1000 km/GeV c ppp = exp(-alpha*eloe) c ppp = exp(-alpha*eloe/2.0) ! for sterile version) c* Decaying nu_mu: c c if(ntyp .eq. 1) then ! nothing happens to electrons c prob = 1.0 c else ! but the muons may go away c s4 = s2th**2 c c4 = (1-s2th)**2 c prob = s4 + c4*ppp c prob = (s2th + (1.0-s2th)*ppp)**2 ! for sterile version c end if c* Harrison, Perkins, Scott oscillations: c c if(ntyp .eq. 1) then c prob = 1.0 -s1*ppp + s2*ppp*rmuel c else c prob = 1.0 -s1*ppp + s2*ppp/rmuel c end if * Vaccum nu_mu -> nutau oscillations: if(ntyp .eq. 1) then ! nothing happens to electrons prob = 1.0 else ! but the muons may go away prob = 1.0 - s2th*ppp end if return end *********************************************************************** subroutine Fmuel(enumev, costheta, rmuel) *********************************************************************** * * get rmuel = ratio of total muons neutrino flux to electron neutrinos * at an energy enumev in MeV, * and cosine zenith angle costheta (negative is upcoming... cr convention) * using Stanev's atmospheric neutrino interpolation routine * implicit real*8 (a-h,o-z) data isol/1/ rmuel = 2.0 if(enumev .le. 50.) return enu = enumev/1000.0 rmuel = 30.0 if(enu .gt. 1000.) return * get electron neutrino flux itype = 1 call kamnus(enu,costheta,itype,flnu,isol) fl1 = flnu * and electron anti-neutrino flux itype = 2 call kamnus(enu,costheta,itype,flnu,isol) fl2 = flnu * get muon neutrino flux itype = 3 call kamnus(enu,costheta,itype,flnu,isol) fl3 = flnu * and muon anti-neutrino flux itype = 4 call kamnus(enu,costheta,itype,flnu,isol) fl4 = flnu totmus = fl3 + fl4 totels = fl1 + fl2 if(totels .gt. 0.0) rmuel = totmus/totels return end ***************************************************************************** *Date: 15 Oct 1996 17:15:18 -0400 *From: stanev@bartol.udel.edu *Content-transfer-encoding: 7BIT * * 1 WRITE(*,*)' ENTER COSTH,ITYPE' * READ (*,*) COSTH,ITYPE * IF (ITYPE .EQ. 0) STOP * ENU = 0.055 * DELTAE = 10.**0.1 * DO I=1,60 * CALL KAMNUS(ENU,COSTH,ITYPE,FLNU) * WRITE(86,*) ENU,FLNU * ENU = DELTAE*ENU * WRITE(*,*) ENU,COSTH,ITYPE,FLNU * ENDDO * GO TO 1 * END C ********************************************************************* SUBROUTINE KAMNUS(ENE,COSTH,ITYPE,FLUX,ISOL) C ********************************************************************* C RETURNS THE FLUXES FOR ELECTRON NEUTRINOS AND ANTINEUTRINOS AND * C MUON NEUTRINOS AND ANTINEUTRINOS (ITYPE RESPECTIVELY FROM 1 TO 4 * C OF ATMOSPHERIC ORIGIN FOR ENERGIES BETWEEN 50 MeV AND 7x10^4 GeV * C FOR THE LOCATION OF KAMIOKA * C QUANTITY IS dN/d(lnE) per cm^2.s.sr. * C SEE KAMIOKA_NU FOR FURTHER INFORMATION. TSS,Oct '96 * C ********************************************************************* implicit real*8 (a-h,o-z) COMMON /FLXINT/ FLX_INT(12,4,63) COMMON /FLXDIF/ FLX_DIF(12,4,63) DIMENSION ANGLES(12),FLN(4),EN(3),FL(3) DATA ANGLES/1.00,0.75,0.50,0.25,0.15,0.05, * -0.05,-0.15,-0.25,-0.50,-0.75,-1.00/ DATA ICALL/0/ IF (ICALL .EQ. 0) THEN ISOLAR = ISOL CALL INIT_FLX(ISOLAR) ICALL = 177 ENDIF ALE = LOG10(ENE) IF (ALE .LT. -1.3) THEN WRITE(*,*)' ENERGY TOO LOW,BELOW 50 MeV' FLUX = 0. RETURN ENDIF IF (ALE .GT. 4.85) THEN WRITE(*,*)' ENERGY TOO HIGH,ABOVE 70000 GeV' FLUX = 0. RETURN ENDIF DO I=1,12 IF (COSTH .GT. ANGLES(I)) GO TO 101 ENDDO 101 IA = I - 1 IF (IA .LT. 2) IA = 2 IF (IA .GT. 11) IA = 11 IE = 10.*(ALE + 1.3) IF (IE .LT. 1) IE = 1 IF (IE .GT. 71) IE = 61 C INTERPOLATE INTEGRAL FLUXES IN ANGLE IRUN = 0 F1 = ANGLES(IA - 1) F2 = ANGLES(IA) F3 = ANGLES(IA + 1) DO IEE = IE,IE + 3 IRUN = IRUN + 1 FF1 = LOG(FLX_INT(IA - 1,ITYPE,IEE)) FF2 = LOG(FLX_INT(IA ,ITYPE,IEE)) FF3 = LOG(FLX_INT(IA + 1,ITYPE,IEE)) FF = QUAD_INT(COSTH,F1,F2,F3,FF1,FF2,FF3) FLN(IRUN) = EXP(FF) ENDDO IRUN = 0 DO IEE = IE,IE + 2 IRUN = IRUN + 1 FF1 = LOG(FLX_DIF(IA - 1,ITYPE,IEE)) FF2 = LOG(FLX_DIF(IA ,ITYPE,IEE)) FF3 = LOG(FLX_DIF(IA + 1,ITYPE,IEE)) FF = QUAD_INT(COSTH,F1,F2,F3,FF1,FF2,FF3) FL(IRUN) = EXP(FF) ENDDO C DONE. FLN(1:4) CONTAINS INTEGRAL NEUTRINO FLUXES INTERPOLATED IN ANGLE. C FL(1:3) CONTAINS THE INTERPOLATED DIFFERENTIAL FLUXES. C C NOW FIND THE LOCAL SPECTRAL INDEX AND ASSIGN THE CORRECT AVERAGE C ENERGIES FOR THE BINS DEL = 10.**0.1 DO I=1,3 RAT = FLN(I)/FLN(I+1) GAMMA = LOG(RAT)/0.230259 + 1. ccc debug 8/04: hanging up on following statement. Set to 1 to see what happens.. c corr = 1.0 CORR = (GAMMA - 1.)*(1.- DEL**(2.- GAMMA))/ * (GAMMA - 2.)/(1.- DEL**(1.- GAMMA)) C ANORM = (IE + I - 2)/10. ANORM = (IE + I - 15)/10. EN(I) = CORR*10.**ANORM ENDDO C NOW GET THE FINAL FLUX VALUES,INTERPOLATE IN ENERGY F1 = LOG10(EN(1)) F2 = LOG10(EN(2)) F3 = LOG10(EN(3)) FF1 = LOG(FL(1)) FF2 = LOG(FL(2)) FF3 = LOG(FL(3)) FF = QUAD_INT(ALE,F1,F2,F3,FF1,FF2,FF3) FLUX = EXP(FF) RETURN END C ********************************************************************* SUBROUTINE INIT_FLX(ISOL) C ********************************************************************* C ADDS THE DIFFERENTIAL NEUTRINO FLUXES TO OBTAIN THE INTEGRAL * C FLUX FOR INTERPOLATION IN ANGLE AND ENERGY. TSS, Nov '93 * C ********************************************************************* implicit real*8 (a-h,o-z) COMMON /FLXINT/ FLX_INT(12,4,63) COMMON /FLXDIF/ FLX_DIF(12,4,63) COMMON /FLXNUS/ FLX_NU(12,4,63), * FLX_NU_MAX(12,4,23),FLX_NU_MIN(12,4,23) DATA ICALL/0/ IF (ICALL .EQ. 0) THEN * WRITE(*,*)' ENTER 1 FOR SOLAR MINIMUM' * WRITE(*,*)' 2 FOR SOLAR MAXIMUM' * READ (*,*) ISOL ICALL = 171 C MAKE DOWNWARD GOING > 10 GeV NEUTRINOS SAME AS UPWARD C GOING FOR THE 6 ANGULAR BINS,NO GEOMAGNETIC EFFECTS. DO J=1,4 DO I=7,12 II = 13 - I DO K=24,63 FLX_NU(I,J,K) = FLX_NU(II,J,K) ENDDO ENDDO ENDDO DO KK = 1,12 DO KL = 1,4 FLX_INT(KK,KL,63) = FLX_NU(KK,KL,63) FLX_DIF(KK,KL,63) = FLX_NU(KK,KL,63) DO I=1,39 II = 63 - I FLX_DIF(KK,KL,II) = FLX_NU(KK,KL,II) FLX_INT(KK,KL,II) = FLX_INT(KK,KL,II+1) + + FLX_NU(KK,KL,II) ENDDO DO I=1,23 II = 24 - I IF (ISOL .EQ. 1) THEN FF = FLX_NU_MIN(KK,KL,II) ELSE FF = FLX_NU_MAX(KK,KL,II) ENDIF FLX_INT(KK,KL,II) = FLX_INT(KK,KL,II+1) + FF FLX_DIF(KK,KL,II) = FF ENDDO ENDDO ENDDO ELSE WRITE(*,*)' What the fuck...................?' WRITE(*,*)' I have done my bit,go away' STOP ENDIF RETURN END C ********************************************************************* BLOCK DATA KAMIOKA_NU C ********************************************************************* C CONTAINS ATMOSPHERIC NEUTRINO AND ANTINEUTRINO FLUXES CALCULATED * C WITH `TARGETNEWT' AND WITH ACCOUNT FOR THE GEOMAGNETIC CUTOFF AT * C THE LOCATION OF KAMIOKA. THE QUANTITY IS dN/d(lnE) per cm^2.s.sr. * C FIRST INDEX (1:12) REFERS TO ZENITH ANGLE [cos(theta) = 1.00,0.75, * C 0.50,0.25,0.15,0.05,-0.05 ,... -1.00]. * C SECOND INDEX IS 1 FOR ELECTRON NEUTRINOS, * C 2 FOR ELECTRON ANTINEUTRINOS, * C 3 FOR MUON NEUTRINOS, * C 4 FOR MUON ANTINEUTRINOS. * C THIRD IS THE ENERGY INDEX -- FIRST BIN FOR log10(E/GeV) FROM * C -1.3 TO -1.2,i.e. from 50 MeV to 63 MeV,etc. up to 100 TeV. * C NOTE THE INPUT TABLES CONTAIN SEPARATELY NEUTRINOS BELOW AND ABOVE * C 10 GeV. GEOMAGNETIC EFFECTS AND SOLAR MODULATION ARE ASSUMED * C NEGLIGIBLE FOR NEUTRINOS ABOVE 10 GeV. SOLAR MAXIMUM AND MINIMUM * C ARE GIVEN SEPARATELY FOR LOWER ENERGY NEUTRINOS. * C REFERENCES: * C Atmospheric neutrino flux above 1 GeV,V. Agraval,T.K. Gaisser, * C P. Lipari and T. Stanev,Phys. Rev. D {\bf 53},1314 (1996). * C An Improved Calculation of the Atmospheric Neutrino Flux, * C T.K.~Gaisser and T.~Stanev,in {Proc. 24th ICRC},(Rome 1995) * C TSS, Oct '96 * C ********************************************************************* implicit real*8 (a-h,o-z) COMMON /FLXNUS/ FLX_NU(12,4,63), * FLX_NU_MAX(12,4,23),FLX_NU_MIN(12,4,23) DATA (FLX_NU( 1,1,LW),LW=24,63)/ * 3.38E-05,1.84E-05,1.02E-05,6.16E-06,3.40E-06,1.89E-06,1.01E-06, * 5.79E-07,3.60E-07,1.98E-07,1.16E-07,6.58E-08,3.66E-08,1.90E-08, * 1.38E-08,7.67E-09,4.41E-09,2.55E-09,1.41E-09,8.42E-10,4.44E-10, * 2.50E-10,1.41E-10,7.72E-11,4.24E-11,2.36E-11,1.30E-11,6.96E-12, * 3.76E-12,2.04E-12,1.12E-12,6.06E-13,3.27E-13,1.76E-13,9.45E-14, * 4.92E-14,2.57E-14,1.34E-14,7.32E-15,3.89E-15/ DATA (FLX_NU( 2,1,LW),LW=24,63)/ * 4.27E-05,2.52E-05,1.32E-05,8.05E-06,4.56E-06,2.40E-06,1.38E-06, * 8.15E-07,4.27E-07,2.46E-07,1.38E-07,8.30E-08,4.80E-08,2.81E-08, * 1.37E-08,9.51E-09,5.00E-09,3.02E-09,1.68E-09,1.19E-09,5.42E-10, * 3.10E-10,1.73E-10,9.80E-11,5.47E-11,3.01E-11,1.66E-11,8.94E-12, * 4.86E-12,2.70E-12,1.44E-12,7.74E-13,4.18E-13,2.26E-13,1.20E-13, * 6.29E-14,3.33E-14,1.81E-14,9.12E-15,4.89E-15/ DATA (FLX_NU( 3,1,LW),LW=24,63)/ * 6.82E-05,3.84E-05,2.05E-05,1.26E-05,7.16E-06,3.75E-06,2.00E-06, * 1.18E-06,6.70E-07,3.48E-07,2.11E-07,1.14E-07,6.45E-08,4.07E-08, * 2.59E-08,1.26E-08,7.17E-09,4.10E-09,2.54E-09,1.44E-09,7.27E-10, * 4.20E-10,2.42E-10,1.36E-10,7.70E-11,4.36E-11,2.42E-11,1.32E-11, * 7.24E-12,3.95E-12,2.12E-12,1.15E-12,6.33E-13,3.42E-13,1.83E-13, * 9.71E-14,5.23E-14,2.73E-14,1.42E-14,7.87E-15/ DATA (FLX_NU( 4,1,LW),LW=24,63)/ * 1.16E-04,6.60E-05,3.91E-05,2.28E-05,1.37E-05,7.15E-06,4.14E-06, * 2.26E-06,1.26E-06,6.78E-07,3.78E-07,2.21E-07,1.15E-07,6.63E-08, * 3.51E-08,2.17E-08,1.19E-08,6.43E-09,3.85E-09,2.14E-09,1.16E-09, * 6.27E-10,3.66E-10,2.11E-10,1.21E-10,6.91E-11,3.92E-11,2.20E-11, * 1.22E-11,6.66E-12,3.73E-12,2.04E-12,1.09E-12,6.00E-13,3.24E-13, * 1.72E-13,8.98E-14,4.92E-14,2.57E-14,1.43E-14/ DATA (FLX_NU( 5,1,LW),LW=24,63)/ * 1.49E-04,8.91E-05,5.33E-05,3.15E-05,1.82E-05,1.04E-05,6.06E-06, * 3.22E-06,1.83E-06,1.03E-06,5.71E-07,3.23E-07,1.68E-07,8.86E-08, * 5.07E-08,2.81E-08,1.54E-08,8.46E-09,4.84E-09,2.41E-09,1.62E-09, * 9.16E-10,5.18E-10,3.02E-10,1.55E-10,8.61E-11,4.86E-11,2.73E-11, * 1.52E-11,8.44E-12,4.63E-12,2.51E-12,1.35E-12,7.51E-13,3.98E-13, * 2.16E-13,1.15E-13,6.13E-14,3.24E-14,1.70E-14/ DATA (FLX_NU( 6,1,LW),LW=24,63)/ * 1.98E-04,1.25E-04,7.72E-05,4.71E-05,3.02E-05,1.73E-05,9.79E-06, * 5.68E-06,3.42E-06,1.88E-06,1.09E-06,5.79E-07,3.25E-07,1.89E-07, * 9.15E-08,5.16E-08,2.77E-08,1.53E-08,8.57E-09,4.82E-09,2.49E-09, * 1.46E-09,7.26E-10,4.05E-10,2.37E-10,1.43E-10,9.30E-11,3.75E-11, * 2.10E-11,1.17E-11,6.84E-12,3.75E-12,2.03E-12,1.12E-12,6.08E-13, * 3.27E-13,1.75E-13,9.40E-14,5.01E-14,2.66E-14/ DATA (FLX_NU( 1,2,LW),LW=24,63)/ * 2.88E-05,1.58E-05,9.03E-06,5.04E-06,2.98E-06,1.51E-06,8.68E-07, * 4.78E-07,2.72E-07,1.73E-07,8.64E-08,4.62E-08,3.00E-08,1.89E-08, * 8.95E-09,5.14E-09,3.44E-09,1.56E-09,1.01E-09,5.09E-10,2.83E-10, * 1.60E-10,8.79E-11,4.78E-11,2.58E-11,1.45E-11,7.70E-12,4.10E-12, * 2.22E-12,1.18E-12,6.47E-13,3.53E-13,1.90E-13,1.01E-13,5.27E-14, * 2.85E-14,1.46E-14,7.56E-15,3.96E-15,2.06E-15/ DATA (FLX_NU( 2,2,LW),LW=24,63)/ * 3.62E-05,1.93E-05,1.15E-05,6.56E-06,3.73E-06,2.00E-06,1.09E-06, * 5.95E-07,3.34E-07,2.18E-07,1.16E-07,6.92E-08,3.44E-08,2.04E-08, * 1.31E-08,6.52E-09,3.60E-09,2.02E-09,1.16E-09,6.39E-10,3.52E-10, * 1.97E-10,1.09E-10,6.04E-11,3.32E-11,1.79E-11,9.79E-12,5.32E-12, * 2.89E-12,1.62E-12,8.48E-13,4.48E-13,2.41E-13,1.28E-13,6.82E-14, * 3.55E-14,1.88E-14,1.03E-14,5.19E-15,2.75E-15/ DATA (FLX_NU( 3,2,LW),LW=24,63)/ * 5.67E-05,3.13E-05,1.74E-05,1.02E-05,5.55E-06,3.12E-06,1.75E-06, * 9.13E-07,5.09E-07,2.65E-07,1.55E-07,9.54E-08,4.71E-08,3.14E-08, * 1.54E-08,8.12E-09,5.70E-09,2.87E-09,1.43E-09,8.83E-10,4.83E-10, * 2.78E-10,1.56E-10,8.93E-11,4.95E-11,2.70E-11,1.48E-11,8.10E-12, * 4.30E-12,2.33E-12,1.27E-12,6.87E-13,3.65E-13,1.95E-13,1.06E-13, * 5.63E-14,2.94E-14,1.59E-14,8.16E-15,4.23E-15/ DATA (FLX_NU( 4,2,LW),LW=24,63)/ * 9.57E-05,5.59E-05,3.28E-05,1.83E-05,1.07E-05,6.26E-06,3.52E-06, * 1.81E-06,9.87E-07,5.91E-07,3.19E-07,1.81E-07,9.60E-08,5.23E-08, * 2.88E-08,1.74E-08,8.81E-09,5.48E-09,2.91E-09,1.43E-09,7.38E-10, * 4.27E-10,2.46E-10,1.38E-10,7.70E-11,4.42E-11,2.47E-11,1.37E-11, * 7.38E-12,4.08E-12,2.25E-12,1.21E-12,6.45E-13,3.44E-13,1.86E-13, * 9.98E-14,5.14E-14,2.73E-14,1.47E-14,8.02E-15/ DATA (FLX_NU( 5,2,LW),LW=24,63)/ * 1.22E-04,7.50E-05,4.45E-05,2.52E-05,1.52E-05,8.68E-06,4.96E-06, * 2.81E-06,1.57E-06,8.55E-07,4.88E-07,2.66E-07,1.29E-07,7.10E-08, * 4.11E-08,2.22E-08,1.25E-08,6.40E-09,3.73E-09,2.06E-09,1.14E-09, * 7.23E-10,3.65E-10,1.81E-10,1.03E-10,5.48E-11,3.09E-11,1.72E-11, * 9.42E-12,5.16E-12,2.76E-12,1.49E-12,8.04E-13,4.40E-13,2.32E-13, * 1.22E-13,6.60E-14,3.37E-14,1.81E-14,9.85E-15/ DATA (FLX_NU( 6,2,LW),LW=24,63)/ * 1.66E-04,1.03E-04,6.28E-05,3.76E-05,2.35E-05,1.40E-05,8.32E-06, * 4.87E-06,2.79E-06,1.48E-06,8.82E-07,4.92E-07,2.73E-07,1.41E-07, * 7.53E-08,4.26E-08,2.34E-08,1.21E-08,6.76E-09,3.54E-09,2.01E-09, * 1.03E-09,5.18E-10,3.20E-10,1.61E-10,9.99E-11,5.42E-11,2.43E-11, * 1.33E-11,7.51E-12,4.12E-12,2.28E-12,1.23E-12,6.71E-13,3.60E-13, * 1.91E-13,1.02E-13,5.34E-14,2.79E-14,1.47E-14/ DATA (FLX_NU( 1,3,LW),LW=24,63)/ * 2.22E-04,1.46E-04,9.11E-05,5.90E-05,3.74E-05,2.35E-05,1.48E-05, * 9.34E-06,5.85E-06,3.62E-06,2.26E-06,1.38E-06,8.79E-07,5.37E-07, * 3.13E-07,1.94E-07,1.21E-07,6.97E-08,4.15E-08,2.44E-08,1.42E-08, * 8.03E-09,4.58E-09,2.59E-09,1.44E-09,8.28E-10,4.55E-10,2.55E-10, * 1.37E-10,7.41E-11,4.02E-11,2.19E-11,1.22E-11,6.58E-12,3.41E-12, * 1.85E-12,1.01E-12,5.21E-13,2.66E-13,1.44E-13/ DATA (FLX_NU( 2,3,LW),LW=24,63)/ * 2.42E-04,1.55E-04,9.66E-05,6.32E-05,4.13E-05,2.61E-05,1.63E-05, * 1.04E-05,6.38E-06,4.10E-06,2.48E-06,1.56E-06,9.28E-07,6.06E-07, * 3.80E-07,2.25E-07,1.32E-07,8.09E-08,4.83E-08,3.03E-08,1.69E-08, * 9.75E-09,5.58E-09,3.22E-09,1.80E-09,1.06E-09,5.81E-10,3.29E-10, * 1.80E-10,9.67E-11,5.26E-11,2.91E-11,1.55E-11,8.34E-12,4.48E-12, * 2.41E-12,1.31E-12,6.85E-13,3.64E-13,1.95E-13/ DATA (FLX_NU( 3,3,LW),LW=24,63)/ * 2.64E-04,1.75E-04,1.13E-04,7.10E-05,4.59E-05,2.95E-05,1.84E-05, * 1.18E-05,7.40E-06,4.72E-06,2.87E-06,1.82E-06,1.12E-06,6.99E-07, * 4.38E-07,2.63E-07,1.68E-07,9.84E-08,6.08E-08,3.75E-08,2.11E-08, * 1.26E-08,7.51E-09,4.34E-09,2.47E-09,1.44E-09,8.29E-10,4.55E-10, * 2.46E-10,1.39E-10,7.79E-11,4.28E-11,2.26E-11,1.22E-11,6.76E-12, * 3.68E-12,1.97E-12,1.03E-12,5.46E-13,2.85E-13/ DATA (FLX_NU( 4,3,LW),LW=24,63)/ * 3.18E-04,2.04E-04,1.35E-04,8.70E-05,5.39E-05,3.55E-05,2.29E-05, * 1.43E-05,9.20E-06,5.80E-06,3.64E-06,2.25E-06,1.45E-06,8.86E-07, * 5.60E-07,3.48E-07,2.19E-07,1.33E-07,8.30E-08,5.03E-08,2.92E-08, * 1.76E-08,1.05E-08,6.32E-09,3.71E-09,2.17E-09,1.26E-09,7.25E-10, * 4.07E-10,2.35E-10,1.32E-10,7.37E-11,4.01E-11,2.16E-11,1.18E-11, * 6.46E-12,3.44E-12,1.73E-12,1.00E-12,5.33E-13/ DATA (FLX_NU( 5,3,LW),LW=24,63)/ * 3.46E-04,2.23E-04,1.48E-04,9.39E-05,6.14E-05,3.92E-05,2.50E-05, * 1.60E-05,1.00E-05,6.46E-06,4.04E-06,2.54E-06,1.60E-06,9.94E-07, * 6.21E-07,3.85E-07,2.40E-07,1.45E-07,9.14E-08,5.51E-08,3.41E-08, * 2.04E-08,1.24E-08,7.52E-09,4.52E-09,2.60E-09,1.55E-09,8.79E-10, * 4.90E-10,2.91E-10,1.59E-10,8.88E-11,4.91E-11,2.65E-11,1.47E-11, * 8.12E-12,4.40E-12,2.28E-12,1.21E-12,6.45E-13/ DATA (FLX_NU( 6,3,LW),LW=24,63)/ * 3.92E-04,2.56E-04,1.68E-04,1.09E-04,7.01E-05,4.52E-05,2.98E-05, * 1.90E-05,1.21E-05,7.61E-06,4.89E-06,3.05E-06,1.94E-06,1.21E-06, * 7.60E-07,4.64E-07,2.86E-07,1.74E-07,1.10E-07,6.60E-08,4.06E-08, * 2.47E-08,1.53E-08,9.04E-09,5.64E-09,3.32E-09,1.99E-09,1.18E-09, * 6.74E-10,3.94E-10,2.27E-10,1.28E-10,7.04E-11,4.04E-11,2.12E-11, * 1.19E-11,6.56E-12,3.42E-12,1.83E-12,1.01E-12/ DATA (FLX_NU( 1,4,LW),LW=24,63)/ * 1.78E-04,1.13E-04,7.23E-05,4.38E-05,2.70E-05,1.70E-05,1.02E-05, * 6.38E-06,3.63E-06,2.32E-06,1.36E-06,8.26E-07,4.89E-07,2.91E-07, * 1.68E-07,1.01E-07,5.87E-08,3.43E-08,2.00E-08,1.15E-08,6.47E-09, * 3.65E-09,2.05E-09,1.17E-09,6.46E-10,3.58E-10,1.99E-10,1.06E-10, * 5.67E-11,2.96E-11,1.72E-11,9.05E-12,4.90E-12,2.53E-12,1.33E-12, * 7.25E-13,3.84E-13,2.04E-13,1.08E-13,5.78E-14/ DATA (FLX_NU( 2,4,LW),LW=24,63)/ * 1.95E-04,1.21E-04,7.90E-05,5.00E-05,2.97E-05,1.90E-05,1.19E-05, * 7.17E-06,4.31E-06,2.61E-06,1.59E-06,9.33E-07,5.73E-07,3.39E-07, * 2.10E-07,1.20E-07,6.99E-08,4.10E-08,2.45E-08,1.46E-08,8.04E-09, * 4.64E-09,2.62E-09,1.45E-09,8.10E-10,4.55E-10,2.51E-10,1.37E-10, * 7.47E-11,3.99E-11,2.19E-11,1.20E-11,6.36E-12,3.29E-12,1.76E-12, * 9.23E-13,5.24E-13,2.73E-13,1.36E-13,7.47E-14/ DATA (FLX_NU( 3,4,LW),LW=24,63)/ * 2.26E-04,1.42E-04,9.22E-05,5.84E-05,3.57E-05,2.22E-05,1.37E-05, * 8.67E-06,5.30E-06,3.10E-06,1.93E-06,1.19E-06,6.93E-07,4.24E-07, * 2.55E-07,1.49E-07,8.62E-08,5.28E-08,3.18E-08,1.82E-08,1.05E-08, * 6.08E-09,3.51E-09,2.03E-09,1.17E-09,6.42E-10,3.58E-10,1.99E-10, * 1.08E-10,6.09E-11,3.29E-11,1.75E-11,9.29E-12,5.03E-12,2.73E-12, * 1.42E-12,7.82E-13,4.32E-13,2.23E-13,1.13E-13/ DATA (FLX_NU( 4,4,LW),LW=24,63)/ * 2.87E-04,1.80E-04,1.17E-04,7.32E-05,4.56E-05,2.87E-05,1.74E-05, * 1.11E-05,6.74E-06,4.13E-06,2.57E-06,1.57E-06,9.52E-07,5.60E-07, * 3.46E-07,2.20E-07,1.30E-07,7.84E-08,4.66E-08,2.59E-08,1.49E-08, * 8.87E-09,5.17E-09,3.07E-09,1.71E-09,9.92E-10,5.76E-10,3.22E-10, * 1.81E-10,9.86E-11,5.42E-11,2.99E-11,1.62E-11,9.03E-12,4.84E-12, * 2.58E-12,1.39E-12,7.64E-13,3.76E-13,1.99E-13/ DATA (FLX_NU( 5,4,LW),LW=24,63)/ * 3.25E-04,2.08E-04,1.32E-04,8.19E-05,5.25E-05,3.29E-05,2.06E-05, * 1.28E-05,7.86E-06,4.83E-06,3.01E-06,1.83E-06,1.12E-06,6.58E-07, * 3.99E-07,2.41E-07,1.44E-07,8.59E-08,5.18E-08,3.04E-08,1.81E-08, * 1.06E-08,6.20E-09,3.75E-09,2.05E-09,1.19E-09,6.87E-10,3.88E-10, * 2.23E-10,1.20E-10,6.84E-11,3.81E-11,2.08E-11,1.09E-11,6.22E-12, * 3.30E-12,1.71E-12,9.22E-13,4.77E-13,2.38E-13/ DATA (FLX_NU( 6,4,LW),LW=24,63)/ * 3.77E-04,2.42E-04,1.57E-04,1.01E-04,6.42E-05,4.08E-05,2.55E-05, * 1.58E-05,1.00E-05,6.31E-06,3.83E-06,2.39E-06,1.45E-06,8.69E-07, * 5.28E-07,3.12E-07,1.91E-07,1.12E-07,6.70E-08,3.92E-08,2.39E-08, * 1.37E-08,8.08E-09,4.77E-09,2.84E-09,1.73E-09,9.39E-10,5.34E-10, * 3.06E-10,1.72E-10,9.66E-11,5.22E-11,2.96E-11,1.62E-11,8.85E-12, * 4.77E-12,2.57E-12,1.39E-12,7.49E-13,3.86E-13/ DATA (FLX_NU_MAX( 1,1,LW),LW=1,23)/ * 2.32E-02,2.34E-02,2.38E-02,2.45E-02,2.40E-02,2.30E-02,2.13E-02, * 1.93E-02,1.68E-02,1.39E-02,1.16E-02,8.78E-03,6.59E-03,4.74E-03, * 3.21E-03,2.17E-03,1.50E-03,8.78E-04,5.55E-04,3.20E-04,1.81E-04, * 1.04E-04,5.95E-05/ DATA (FLX_NU_MAX( 2,1,LW),LW=1,23)/ * 2.03E-02,2.05E-02,2.17E-02,2.24E-02,2.23E-02,2.13E-02,2.05E-02, * 1.86E-02,1.65E-02,1.42E-02,1.15E-02,9.40E-03,7.17E-03,5.35E-03, * 3.73E-03,2.57E-03,1.68E-03,1.11E-03,6.81E-04,4.18E-04,2.53E-04, * 1.34E-04,8.40E-05/ DATA (FLX_NU_MAX( 3,1,LW),LW=1,23)/ * 1.99E-02,2.09E-02,2.21E-02,2.29E-02,2.32E-02,2.27E-02,2.17E-02, * 2.03E-02,1.81E-02,1.58E-02,1.30E-02,1.05E-02,8.20E-03,6.24E-03, * 4.57E-03,3.17E-03,2.21E-03,1.42E-03,9.03E-04,5.80E-04,3.51E-04, * 2.03E-04,1.26E-04/ DATA (FLX_NU_MAX( 4,1,LW),LW=1,23)/ * 1.65E-02,1.76E-02,1.91E-02,1.99E-02,2.04E-02,2.04E-02,1.97E-02, * 1.85E-02,1.71E-02,1.50E-02,1.29E-02,1.08E-02,8.66E-03,6.89E-03, * 5.23E-03,3.81E-03,2.69E-03,1.78E-03,1.22E-03,7.67E-04,4.92E-04, * 3.17E-04,1.96E-04/ DATA (FLX_NU_MAX( 5,1,LW),LW=1,23)/ * 1.54E-02,1.69E-02,1.81E-02,1.92E-02,1.93E-02,1.94E-02,1.88E-02, * 1.76E-02,1.64E-02,1.46E-02,1.28E-02,1.08E-02,8.65E-03,6.81E-03, * 5.32E-03,3.97E-03,2.84E-03,1.99E-03,1.39E-03,9.16E-04,6.11E-04, * 3.82E-04,2.43E-04/ DATA (FLX_NU_MAX( 6,1,LW),LW=1,23)/ * 1.42E-02,1.57E-02,1.69E-02,1.78E-02,1.83E-02,1.81E-02,1.76E-02, * 1.66E-02,1.56E-02,1.39E-02,1.23E-02,1.01E-02,8.51E-03,6.63E-03, * 5.30E-03,4.02E-03,2.95E-03,2.10E-03,1.54E-03,1.05E-03,7.04E-04, * 4.95E-04,2.92E-04/ DATA (FLX_NU_MAX( 7,1,LW),LW=1,23)/ * 1.31E-02,1.45E-02,1.57E-02,1.65E-02,1.70E-02,1.69E-02,1.64E-02, * 1.56E-02,1.47E-02,1.31E-02,1.16E-02,9.62E-03,8.14E-03,6.39E-03, * 5.13E-03,3.92E-03,2.89E-03,2.07E-03,1.53E-03,1.05E-03,7.03E-04, * 4.94E-04,2.92E-04/ DATA (FLX_NU_MAX( 8,1,LW),LW=1,23)/ * 1.57E-02,1.71E-02,1.84E-02,1.93E-02,1.93E-02,1.92E-02,1.85E-02, * 1.74E-02,1.61E-02,1.42E-02,1.24E-02,1.05E-02,8.38E-03,6.62E-03, * 5.16E-03,3.88E-03,2.79E-03,1.96E-03,1.37E-03,9.10E-04,6.09E-04, * 3.81E-04,2.43E-04/ DATA (FLX_NU_MAX( 9,1,LW),LW=1,23)/ * 2.46E-02,2.55E-02,2.68E-02,2.69E-02,2.66E-02,2.54E-02,2.35E-02, * 2.14E-02,1.93E-02,1.64E-02,1.38E-02,1.13E-02,8.93E-03,6.95E-03, * 5.23E-03,3.78E-03,2.66E-03,1.75E-03,1.21E-03,7.59E-04,4.87E-04, * 3.15E-04,1.95E-04/ DATA (FLX_NU_MAX(10,1,LW),LW=1,23)/ * 3.32E-02,3.32E-02,3.40E-02,3.38E-02,3.29E-02,3.07E-02,2.82E-02, * 2.53E-02,2.19E-02,1.86E-02,1.49E-02,1.17E-02,8.92E-03,6.61E-03, * 4.75E-03,3.25E-03,2.23E-03,1.43E-03,9.01E-04,5.77E-04,3.49E-04, * 2.03E-04,1.25E-04/ DATA (FLX_NU_MAX(11,1,LW),LW=1,23)/ * 4.43E-02,4.25E-02,4.32E-02,4.26E-02,3.99E-02,3.65E-02,3.28E-02, * 2.84E-02,2.41E-02,1.97E-02,1.54E-02,1.20E-02,8.65E-03,6.26E-03, * 4.23E-03,2.78E-03,1.79E-03,1.15E-03,6.89E-04,4.18E-04,2.54E-04, * 1.34E-04,8.40E-05/ DATA (FLX_NU_MAX(12,1,LW),LW=1,23)/ * 3.69E-02,3.60E-02,3.58E-02,3.59E-02,3.43E-02,3.19E-02,2.85E-02, * 2.51E-02,2.13E-02,1.71E-02,1.38E-02,1.03E-02,7.41E-03,5.23E-03, * 3.43E-03,2.29E-03,1.54E-03,8.87E-04,5.56E-04,3.21E-04,1.81E-04, * 1.04E-04,5.95E-05/ DATA (FLX_NU_MAX( 1,2,LW),LW=1,23)/ * 2.24E-02,2.29E-02,2.34E-02,2.35E-02,2.32E-02,2.22E-02,2.05E-02, * 1.83E-02,1.58E-02,1.32E-02,1.05E-02,7.95E-03,5.88E-03,4.02E-03, * 2.92E-03,1.83E-03,1.17E-03,7.05E-04,4.33E-04,2.75E-04,1.51E-04, * 9.30E-05,4.53E-05/ DATA (FLX_NU_MAX( 2,2,LW),LW=1,23)/ * 1.96E-02,2.02E-02,2.11E-02,2.16E-02,2.16E-02,2.06E-02,1.97E-02, * 1.79E-02,1.57E-02,1.32E-02,1.09E-02,8.48E-03,6.31E-03,4.61E-03, * 3.16E-03,2.17E-03,1.42E-03,9.35E-04,5.70E-04,3.29E-04,2.02E-04, * 1.15E-04,6.39E-05/ DATA (FLX_NU_MAX( 3,2,LW),LW=1,23)/ * 1.93E-02,2.02E-02,2.16E-02,2.22E-02,2.23E-02,2.18E-02,2.07E-02, * 1.91E-02,1.71E-02,1.46E-02,1.20E-02,9.52E-03,7.42E-03,5.40E-03, * 3.92E-03,2.79E-03,1.89E-03,1.17E-03,7.47E-04,4.49E-04,2.80E-04, * 1.61E-04,9.18E-05/ DATA (FLX_NU_MAX( 4,2,LW),LW=1,23)/ * 1.60E-02,1.72E-02,1.84E-02,1.94E-02,1.98E-02,1.97E-02,1.87E-02, * 1.77E-02,1.61E-02,1.43E-02,1.18E-02,9.77E-03,7.85E-03,6.07E-03, * 4.47E-03,3.21E-03,2.26E-03,1.52E-03,9.85E-04,6.56E-04,4.01E-04, * 2.62E-04,1.61E-04/ DATA (FLX_NU_MAX( 5,2,LW),LW=1,23)/ * 1.50E-02,1.64E-02,1.75E-02,1.84E-02,1.86E-02,1.86E-02,1.80E-02, * 1.71E-02,1.54E-02,1.37E-02,1.19E-02,9.73E-03,7.79E-03,6.06E-03, * 4.65E-03,3.42E-03,2.36E-03,1.67E-03,1.17E-03,7.64E-04,4.90E-04, * 3.14E-04,1.98E-04/ DATA (FLX_NU_MAX( 6,2,LW),LW=1,23)/ * 1.40E-02,1.53E-02,1.64E-02,1.72E-02,1.76E-02,1.74E-02,1.72E-02, * 1.60E-02,1.48E-02,1.31E-02,1.11E-02,9.36E-03,7.59E-03,5.97E-03, * 4.56E-03,3.46E-03,2.51E-03,1.74E-03,1.32E-03,8.53E-04,5.75E-04, * 4.02E-04,2.39E-04/ DATA (FLX_NU_MAX( 7,2,LW),LW=1,23)/ * 1.30E-02,1.42E-02,1.53E-02,1.60E-02,1.64E-02,1.62E-02,1.61E-02, * 1.50E-02,1.40E-02,1.24E-02,1.06E-02,8.95E-03,7.31E-03,5.78E-03, * 4.45E-03,3.39E-03,2.48E-03,1.73E-03,1.31E-03,8.49E-04,5.74E-04, * 4.02E-04,2.39E-04/ DATA (FLX_NU_MAX( 8,2,LW),LW=1,23)/ * 1.53E-02,1.66E-02,1.77E-02,1.84E-02,1.86E-02,1.84E-02,1.77E-02, * 1.67E-02,1.50E-02,1.33E-02,1.15E-02,9.44E-03,7.55E-03,5.90E-03, * 4.54E-03,3.35E-03,2.32E-03,1.65E-03,1.16E-03,7.61E-04,4.89E-04, * 3.14E-04,1.98E-04/ DATA (FLX_NU_MAX( 9,2,LW),LW=1,23)/ * 2.32E-02,2.41E-02,2.50E-02,2.53E-02,2.49E-02,2.39E-02,2.18E-02, * 1.98E-02,1.75E-02,1.52E-02,1.23E-02,1.00E-02,7.94E-03,6.08E-03, * 4.45E-03,3.18E-03,2.23E-03,1.50E-03,9.76E-04,6.50E-04,3.99E-04, * 2.61E-04,1.61E-04/ DATA (FLX_NU_MAX(10,2,LW),LW=1,23)/ * 3.06E-02,3.07E-02,3.13E-02,3.12E-02,3.01E-02,2.83E-02,2.58E-02, * 2.30E-02,2.00E-02,1.66E-02,1.32E-02,1.03E-02,7.87E-03,5.64E-03, * 4.01E-03,2.82E-03,1.90E-03,1.17E-03,7.47E-04,4.48E-04,2.79E-04, * 1.60E-04,9.16E-05/ DATA (FLX_NU_MAX(11,2,LW),LW=1,23)/ * 4.07E-02,3.98E-02,3.97E-02,3.89E-02,3.70E-02,3.32E-02,3.01E-02, * 2.59E-02,2.17E-02,1.74E-02,1.38E-02,1.02E-02,7.32E-03,5.14E-03, * 3.46E-03,2.30E-03,1.47E-03,9.47E-04,5.78E-04,3.32E-04,2.01E-04, * 1.15E-04,6.39E-05/ DATA (FLX_NU_MAX(12,2,LW),LW=1,23)/ * 3.47E-02,3.43E-02,3.45E-02,3.37E-02,3.23E-02,2.97E-02,2.67E-02, * 2.32E-02,1.93E-02,1.57E-02,1.21E-02,8.89E-03,6.40E-03,4.31E-03, * 3.04E-03,1.89E-03,1.20E-03,7.10E-04,4.33E-04,2.76E-04,1.51E-04, * 9.30E-05,4.53E-05/ DATA (FLX_NU_MAX( 1,3,LW),LW=1,23)/ * 4.56E-02,4.62E-02,4.88E-02,5.07E-02,4.97E-02,4.86E-02,4.64E-02, * 4.25E-02,3.76E-02,3.24E-02,2.65E-02,2.10E-02,1.63E-02,1.20E-02, * 8.59E-03,6.10E-03,4.11E-03,2.77E-03,1.83E-03,1.20E-03,8.08E-04, * 5.41E-04,3.33E-04/ DATA (FLX_NU_MAX( 2,3,LW),LW=1,23)/ * 3.93E-02,4.07E-02,4.34E-02,4.48E-02,4.49E-02,4.43E-02,4.23E-02, * 3.95E-02,3.52E-02,3.03E-02,2.56E-02,2.04E-02,1.60E-02,1.22E-02, * 8.83E-03,6.33E-03,4.23E-03,2.94E-03,2.02E-03,1.33E-03,8.48E-04, * 5.65E-04,3.68E-04/ DATA (FLX_NU_MAX( 3,3,LW),LW=1,23)/ * 3.95E-02,4.15E-02,4.42E-02,4.61E-02,4.61E-02,4.53E-02,4.35E-02, * 4.11E-02,3.67E-02,3.14E-02,2.65E-02,2.15E-02,1.69E-02,1.28E-02, * 9.35E-03,6.77E-03,4.79E-03,3.26E-03,2.18E-03,1.44E-03,9.70E-04, * 6.08E-04,4.14E-04/ DATA (FLX_NU_MAX( 4,3,LW),LW=1,23)/ * 3.34E-02,3.58E-02,3.83E-02,4.01E-02,4.10E-02,4.05E-02,3.91E-02, * 3.67E-02,3.33E-02,2.93E-02,2.48E-02,2.06E-02,1.64E-02,1.27E-02, * 9.48E-03,7.01E-03,5.06E-03,3.52E-03,2.37E-03,1.63E-03,1.09E-03, * 7.25E-04,4.82E-04/ DATA (FLX_NU_MAX( 5,3,LW),LW=1,23)/ * 3.08E-02,3.38E-02,3.66E-02,3.82E-02,3.89E-02,3.83E-02,3.71E-02, * 3.49E-02,3.19E-02,2.81E-02,2.39E-02,2.00E-02,1.63E-02,1.25E-02, * 9.30E-03,7.06E-03,5.07E-03,3.57E-03,2.51E-03,1.73E-03,1.18E-03, * 7.93E-04,5.24E-04/ DATA (FLX_NU_MAX( 6,3,LW),LW=1,23)/ * 2.88E-02,3.17E-02,3.43E-02,3.59E-02,3.61E-02,3.60E-02,3.48E-02, * 3.32E-02,3.03E-02,2.65E-02,2.29E-02,1.91E-02,1.55E-02,1.20E-02, * 9.26E-03,6.85E-03,5.02E-03,3.63E-03,2.61E-03,1.83E-03,1.26E-03, * 8.54E-04,5.77E-04/ DATA (FLX_NU_MAX( 7,3,LW),LW=1,23)/ * 2.66E-02,2.93E-02,3.17E-02,3.33E-02,3.36E-02,3.35E-02,3.26E-02, * 3.12E-02,2.85E-02,2.51E-02,2.17E-02,1.83E-02,1.49E-02,1.16E-02, * 9.02E-03,6.72E-03,4.95E-03,3.60E-03,2.59E-03,1.83E-03,1.26E-03, * 8.54E-04,5.77E-04/ DATA (FLX_NU_MAX( 8,3,LW),LW=1,23)/ * 3.15E-02,3.43E-02,3.70E-02,3.85E-02,3.87E-02,3.80E-02,3.66E-02, * 3.42E-02,3.11E-02,2.73E-02,2.32E-02,1.94E-02,1.58E-02,1.21E-02, * 9.09E-03,6.92E-03,4.99E-03,3.53E-03,2.49E-03,1.73E-03,1.17E-03, * 7.93E-04,5.24E-04/ DATA (FLX_NU_MAX( 9,3,LW),LW=1,23)/ * 4.97E-02,5.16E-02,5.34E-02,5.37E-02,5.23E-02,4.95E-02,4.59E-02, * 4.15E-02,3.65E-02,3.13E-02,2.60E-02,2.11E-02,1.66E-02,1.27E-02, * 9.40E-03,6.92E-03,4.98E-03,3.47E-03,2.35E-03,1.61E-03,1.08E-03, * 7.21E-04,4.80E-04/ DATA (FLX_NU_MAX(10,3,LW),LW=1,23)/ * 6.58E-02,6.56E-02,6.74E-02,6.73E-02,6.40E-02,5.99E-02,5.53E-02, * 5.02E-02,4.32E-02,3.59E-02,2.95E-02,2.34E-02,1.79E-02,1.32E-02, * 9.55E-03,6.83E-03,4.81E-03,3.25E-03,2.17E-03,1.43E-03,9.68E-04, * 6.07E-04,4.14E-04/ DATA (FLX_NU_MAX(11,3,LW),LW=1,23)/ * 8.73E-02,8.55E-02,8.72E-02,8.52E-02,7.98E-02,7.36E-02,6.61E-02, * 5.84E-02,4.92E-02,4.03E-02,3.25E-02,2.48E-02,1.86E-02,1.36E-02, * 9.51E-03,6.66E-03,4.35E-03,2.98E-03,2.03E-03,1.33E-03,8.49E-04, * 5.65E-04,3.68E-04/ DATA (FLX_NU_MAX(12,3,LW),LW=1,23)/ * 7.36E-02,7.23E-02,7.46E-02,7.47E-02,7.07E-02,6.65E-02,6.11E-02, * 5.43E-02,4.63E-02,3.87E-02,3.06E-02,2.36E-02,1.77E-02,1.27E-02, * 8.97E-03,6.24E-03,4.14E-03,2.78E-03,1.83E-03,1.20E-03,8.07E-04, * 5.41E-04,3.33E-04/ DATA (FLX_NU_MAX( 1,4,LW),LW=1,23)/ * 4.60E-02,4.66E-02,4.92E-02,5.04E-02,5.05E-02,4.87E-02,4.61E-02, * 4.25E-02,3.71E-02,3.15E-02,2.57E-02,2.07E-02,1.59E-02,1.14E-02, * 8.24E-03,5.65E-03,3.92E-03,2.52E-03,1.67E-03,1.04E-03,6.86E-04, * 4.50E-04,2.75E-04/ DATA (FLX_NU_MAX( 2,4,LW),LW=1,23)/ * 3.93E-02,4.10E-02,4.38E-02,4.51E-02,4.49E-02,4.43E-02,4.25E-02, * 3.90E-02,3.48E-02,3.03E-02,2.53E-02,1.99E-02,1.56E-02,1.16E-02, * 8.55E-03,5.80E-03,4.10E-03,2.75E-03,1.84E-03,1.19E-03,7.49E-04, * 4.81E-04,3.15E-04/ DATA (FLX_NU_MAX( 3,4,LW),LW=1,23)/ * 3.94E-02,4.15E-02,4.48E-02,4.60E-02,4.65E-02,4.59E-02,4.38E-02, * 4.07E-02,3.65E-02,3.15E-02,2.64E-02,2.12E-02,1.64E-02,1.25E-02, * 9.28E-03,6.57E-03,4.42E-03,3.14E-03,2.07E-03,1.35E-03,8.93E-04, * 5.45E-04,3.43E-04/ DATA (FLX_NU_MAX( 4,4,LW),LW=1,23)/ * 3.32E-02,3.57E-02,3.85E-02,4.02E-02,4.10E-02,4.04E-02,3.94E-02, * 3.70E-02,3.31E-02,2.94E-02,2.48E-02,2.04E-02,1.62E-02,1.26E-02, * 9.47E-03,7.01E-03,4.99E-03,3.38E-03,2.35E-03,1.58E-03,1.04E-03, * 6.86E-04,4.28E-04/ DATA (FLX_NU_MAX( 5,4,LW),LW=1,23)/ * 3.10E-02,3.39E-02,3.63E-02,3.85E-02,3.87E-02,3.85E-02,3.72E-02, * 3.51E-02,3.21E-02,2.80E-02,2.41E-02,1.97E-02,1.58E-02,1.25E-02, * 9.45E-03,7.01E-03,5.06E-03,3.59E-03,2.52E-03,1.69E-03,1.14E-03, * 7.66E-04,4.97E-04/ DATA (FLX_NU_MAX( 6,4,LW),LW=1,23)/ * 2.89E-02,3.17E-02,3.44E-02,3.60E-02,3.64E-02,3.62E-02,3.49E-02, * 3.27E-02,2.99E-02,2.66E-02,2.28E-02,1.91E-02,1.51E-02,1.20E-02, * 9.30E-03,6.95E-03,5.04E-03,3.65E-03,2.56E-03,1.84E-03,1.25E-03, * 8.74E-04,5.60E-04/ DATA (FLX_NU_MAX( 7,4,LW),LW=1,23)/ * 2.67E-02,2.93E-02,3.18E-02,3.34E-02,3.39E-02,3.38E-02,3.27E-02, * 3.07E-02,2.81E-02,2.52E-02,2.17E-02,1.83E-02,1.45E-02,1.17E-02, * 9.05E-03,6.81E-03,4.96E-03,3.62E-03,2.54E-03,1.83E-03,1.25E-03, * 8.73E-04,5.60E-04/ DATA (FLX_NU_MAX( 8,4,LW),LW=1,23)/ * 3.17E-02,3.45E-02,3.67E-02,3.87E-02,3.86E-02,3.81E-02,3.66E-02, * 3.43E-02,3.12E-02,2.73E-02,2.34E-02,1.91E-02,1.54E-02,1.22E-02, * 9.23E-03,6.86E-03,4.98E-03,3.55E-03,2.51E-03,1.69E-03,1.14E-03, * 7.65E-04,4.97E-04/ DATA (FLX_NU_MAX( 9,4,LW),LW=1,23)/ * 4.95E-02,5.13E-02,5.34E-02,5.37E-02,5.20E-02,4.92E-02,4.61E-02, * 4.18E-02,3.64E-02,3.15E-02,2.60E-02,2.10E-02,1.64E-02,1.26E-02, * 9.42E-03,6.93E-03,4.91E-03,3.34E-03,2.32E-03,1.57E-03,1.04E-03, * 6.83E-04,4.28E-04/ DATA (FLX_NU_MAX(10,4,LW),LW=1,23)/ * 6.59E-02,6.56E-02,6.78E-02,6.72E-02,6.44E-02,6.05E-02,5.54E-02, * 4.95E-02,4.32E-02,3.60E-02,2.94E-02,2.30E-02,1.74E-02,1.31E-02, * 9.53E-03,6.65E-03,4.45E-03,3.14E-03,2.07E-03,1.35E-03,8.90E-04, * 5.45E-04,3.43E-04/ DATA (FLX_NU_MAX(11,4,LW),LW=1,23)/ * 8.78E-02,8.58E-02,8.71E-02,8.52E-02,7.93E-02,7.31E-02,6.61E-02, * 5.77E-02,4.88E-02,4.01E-02,3.20E-02,2.42E-02,1.81E-02,1.30E-02, * 9.31E-03,6.12E-03,4.23E-03,2.80E-03,1.85E-03,1.19E-03,7.49E-04, * 4.81E-04,3.15E-04/ DATA (FLX_NU_MAX(12,4,LW),LW=1,23)/ * 7.40E-02,7.26E-02,7.48E-02,7.46E-02,7.12E-02,6.65E-02,6.05E-02, * 5.41E-02,4.55E-02,3.75E-02,2.96E-02,2.32E-02,1.73E-02,1.21E-02, * 8.59E-03,5.80E-03,3.96E-03,2.54E-03,1.68E-03,1.04E-03,6.86E-04, * 4.50E-04,2.75E-04/ DATA (FLX_NU_MIN( 1,1,LW),LW=1,23)/ * 2.35E-02,2.37E-02,2.41E-02,2.48E-02,2.43E-02,2.32E-02,2.16E-02, * 1.95E-02,1.70E-02,1.41E-02,1.17E-02,8.87E-03,6.65E-03,4.78E-03, * 3.23E-03,2.18E-03,1.50E-03,8.81E-04,5.56E-04,3.21E-04,1.81E-04, * 1.04E-04,5.95E-05/ DATA (FLX_NU_MIN( 2,1,LW),LW=1,23)/ * 2.06E-02,2.08E-02,2.20E-02,2.27E-02,2.26E-02,2.16E-02,2.07E-02, * 1.88E-02,1.67E-02,1.43E-02,1.16E-02,9.49E-03,7.23E-03,5.38E-03, * 3.76E-03,2.58E-03,1.69E-03,1.12E-03,6.82E-04,4.18E-04,2.53E-04, * 1.34E-04,8.40E-05/ DATA (FLX_NU_MIN( 3,1,LW),LW=1,23)/ * 2.07E-02,2.17E-02,2.30E-02,2.38E-02,2.41E-02,2.35E-02,2.24E-02, * 2.10E-02,1.87E-02,1.62E-02,1.34E-02,1.08E-02,8.37E-03,6.35E-03, * 4.65E-03,3.21E-03,2.22E-03,1.43E-03,9.04E-04,5.80E-04,3.51E-04, * 2.03E-04,1.26E-04/ DATA (FLX_NU_MIN( 4,1,LW),LW=1,23)/ * 1.74E-02,1.85E-02,2.00E-02,2.07E-02,2.13E-02,2.12E-02,2.04E-02, * 1.92E-02,1.77E-02,1.55E-02,1.33E-02,1.10E-02,8.85E-03,7.02E-03, * 5.30E-03,3.85E-03,2.71E-03,1.79E-03,1.23E-03,7.68E-04,4.92E-04, * 3.17E-04,1.96E-04/ DATA (FLX_NU_MIN( 5,1,LW),LW=1,23)/ * 1.62E-02,1.77E-02,1.90E-02,2.01E-02,2.02E-02,2.02E-02,1.95E-02, * 1.83E-02,1.70E-02,1.50E-02,1.32E-02,1.11E-02,8.84E-03,6.93E-03, * 5.40E-03,4.01E-03,2.87E-03,2.00E-03,1.39E-03,9.17E-04,6.11E-04, * 3.82E-04,2.43E-04/ DATA (FLX_NU_MIN( 6,1,LW),LW=1,23)/ * 1.49E-02,1.64E-02,1.77E-02,1.86E-02,1.91E-02,1.89E-02,1.83E-02, * 1.72E-02,1.62E-02,1.44E-02,1.26E-02,1.04E-02,8.71E-03,6.75E-03, * 5.38E-03,4.07E-03,2.97E-03,2.10E-03,1.55E-03,1.05E-03,7.04E-04, * 4.95E-04,2.92E-04/ DATA (FLX_NU_MIN( 7,1,LW),LW=1,23)/ * 1.37E-02,1.51E-02,1.63E-02,1.72E-02,1.77E-02,1.76E-02,1.70E-02, * 1.61E-02,1.51E-02,1.35E-02,1.19E-02,9.86E-03,8.30E-03,6.49E-03, * 5.20E-03,3.95E-03,2.90E-03,2.07E-03,1.53E-03,1.05E-03,7.03E-04, * 4.94E-04,2.92E-04/ DATA (FLX_NU_MIN( 8,1,LW),LW=1,23)/ * 1.77E-02,1.92E-02,2.05E-02,2.14E-02,2.14E-02,2.11E-02,2.01E-02, * 1.87E-02,1.72E-02,1.51E-02,1.31E-02,1.09E-02,8.67E-03,6.80E-03, * 5.26E-03,3.93E-03,2.82E-03,1.97E-03,1.37E-03,9.11E-04,6.09E-04, * 3.81E-04,2.43E-04/ DATA (FLX_NU_MIN( 9,1,LW),LW=1,23)/ * 3.75E-02,3.82E-02,3.91E-02,3.81E-02,3.64E-02,3.33E-02,2.96E-02, * 2.59E-02,2.26E-02,1.87E-02,1.54E-02,1.23E-02,9.55E-03,7.28E-03, * 5.41E-03,3.86E-03,2.69E-03,1.77E-03,1.21E-03,7.59E-04,4.87E-04, * 3.15E-04,1.95E-04/ DATA (FLX_NU_MIN(10,1,LW),LW=1,23)/ * 5.61E-02,5.32E-02,5.21E-02,4.97E-02,4.63E-02,4.14E-02,3.61E-02, * 3.10E-02,2.59E-02,2.13E-02,1.67E-02,1.29E-02,9.64E-03,6.99E-03, * 4.95E-03,3.36E-03,2.27E-03,1.44E-03,9.03E-04,5.77E-04,3.49E-04, * 2.03E-04,1.25E-04/ DATA (FLX_NU_MIN(11,1,LW),LW=1,23)/ * 7.43E-02,6.88E-02,6.78E-02,6.46E-02,5.82E-02,5.11E-02,4.37E-02, * 3.62E-02,2.98E-02,2.36E-02,1.79E-02,1.35E-02,9.48E-03,6.75E-03, * 4.49E-03,2.87E-03,1.83E-03,1.17E-03,6.91E-04,4.18E-04,2.54E-04, * 1.34E-04,8.40E-05/ DATA (FLX_NU_MIN(12,1,LW),LW=1,23)/ * 4.94E-02,4.74E-02,4.67E-02,4.61E-02,4.32E-02,3.95E-02,3.45E-02, * 2.98E-02,2.49E-02,1.96E-02,1.56E-02,1.13E-02,8.01E-03,5.59E-03, * 3.59E-03,2.38E-03,1.58E-03,8.94E-04,5.58E-04,3.21E-04,1.81E-04, * 1.04E-04,5.95E-05/ DATA (FLX_NU_MIN( 1,2,LW),LW=1,23)/ * 2.26E-02,2.31E-02,2.36E-02,2.37E-02,2.34E-02,2.23E-02,2.07E-02, * 1.85E-02,1.59E-02,1.33E-02,1.05E-02,7.98E-03,5.90E-03,4.02E-03, * 2.92E-03,1.83E-03,1.17E-03,7.05E-04,4.33E-04,2.75E-04,1.51E-04, * 9.30E-05,4.53E-05/ DATA (FLX_NU_MIN( 2,2,LW),LW=1,23)/ * 1.98E-02,2.04E-02,2.13E-02,2.18E-02,2.18E-02,2.08E-02,1.99E-02, * 1.81E-02,1.58E-02,1.33E-02,1.10E-02,8.49E-03,6.32E-03,4.61E-03, * 3.16E-03,2.17E-03,1.42E-03,9.34E-04,5.70E-04,3.29E-04,2.02E-04, * 1.15E-04,6.39E-05/ DATA (FLX_NU_MIN( 3,2,LW),LW=1,23)/ * 2.00E-02,2.09E-02,2.23E-02,2.30E-02,2.30E-02,2.25E-02,2.13E-02, * 1.96E-02,1.76E-02,1.49E-02,1.22E-02,9.64E-03,7.52E-03,5.45E-03, * 3.94E-03,2.79E-03,1.90E-03,1.17E-03,7.49E-04,4.50E-04,2.80E-04, * 1.61E-04,9.18E-05/ DATA (FLX_NU_MIN( 4,2,LW),LW=1,23)/ * 1.67E-02,1.79E-02,1.92E-02,2.02E-02,2.06E-02,2.04E-02,1.93E-02, * 1.82E-02,1.65E-02,1.46E-02,1.20E-02,9.93E-03,7.96E-03,6.14E-03, * 4.50E-03,3.23E-03,2.27E-03,1.52E-03,9.85E-04,6.56E-04,4.01E-04, * 2.62E-04,1.61E-04/ DATA (FLX_NU_MIN( 5,2,LW),LW=1,23)/ * 1.56E-02,1.71E-02,1.82E-02,1.91E-02,1.93E-02,1.93E-02,1.87E-02, * 1.76E-02,1.58E-02,1.40E-02,1.21E-02,9.89E-03,7.90E-03,6.12E-03, * 4.69E-03,3.44E-03,2.36E-03,1.68E-03,1.17E-03,7.65E-04,4.90E-04, * 3.14E-04,1.98E-04/ DATA (FLX_NU_MIN( 6,2,LW),LW=1,23)/ * 1.47E-02,1.60E-02,1.71E-02,1.79E-02,1.83E-02,1.80E-02,1.78E-02, * 1.65E-02,1.53E-02,1.34E-02,1.14E-02,9.52E-03,7.69E-03,6.04E-03, * 4.59E-03,3.48E-03,2.52E-03,1.75E-03,1.32E-03,8.53E-04,5.75E-04, * 4.02E-04,2.39E-04/ DATA (FLX_NU_MIN( 7,2,LW),LW=1,23)/ * 1.35E-02,1.47E-02,1.58E-02,1.66E-02,1.70E-02,1.68E-02,1.66E-02, * 1.54E-02,1.43E-02,1.27E-02,1.08E-02,9.07E-03,7.39E-03,5.83E-03, * 4.47E-03,3.40E-03,2.48E-03,1.73E-03,1.31E-03,8.49E-04,5.74E-04, * 4.02E-04,2.39E-04/ DATA (FLX_NU_MIN( 8,2,LW),LW=1,23)/ * 1.70E-02,1.84E-02,1.96E-02,2.02E-02,2.03E-02,1.99E-02,1.90E-02, * 1.78E-02,1.59E-02,1.39E-02,1.19E-02,9.71E-03,7.72E-03,6.00E-03, * 4.59E-03,3.37E-03,2.33E-03,1.66E-03,1.16E-03,7.62E-04,4.89E-04, * 3.14E-04,1.98E-04/ DATA (FLX_NU_MIN( 9,2,LW),LW=1,23)/ * 3.25E-02,3.33E-02,3.37E-02,3.32E-02,3.17E-02,2.94E-02,2.60E-02, * 2.29E-02,1.98E-02,1.67E-02,1.33E-02,1.06E-02,8.30E-03,6.28E-03, * 4.54E-03,3.22E-03,2.25E-03,1.50E-03,9.77E-04,6.50E-04,3.99E-04, * 2.61E-04,1.61E-04/ DATA (FLX_NU_MIN(10,2,LW),LW=1,23)/ * 4.43E-02,4.31E-02,4.26E-02,4.13E-02,3.85E-02,3.50E-02,3.09E-02, * 2.67E-02,2.27E-02,1.83E-02,1.44E-02,1.10E-02,8.28E-03,5.86E-03, * 4.11E-03,2.85E-03,1.92E-03,1.18E-03,7.50E-04,4.48E-04,2.79E-04, * 1.60E-04,9.16E-05/ DATA (FLX_NU_MIN(11,2,LW),LW=1,23)/ * 6.09E-02,5.79E-02,5.62E-02,5.39E-02,4.95E-02,4.29E-02,3.75E-02, * 3.13E-02,2.56E-02,2.00E-02,1.53E-02,1.11E-02,7.81E-03,5.38E-03, * 3.60E-03,2.35E-03,1.48E-03,9.49E-04,5.81E-04,3.32E-04,2.01E-04, * 1.15E-04,6.39E-05/ DATA (FLX_NU_MIN(12,2,LW),LW=1,23)/ * 4.50E-02,4.36E-02,4.35E-02,4.19E-02,3.94E-02,3.54E-02,3.13E-02, * 2.67E-02,2.18E-02,1.74E-02,1.31E-02,9.51E-03,6.74E-03,4.51E-03, * 3.12E-03,1.93E-03,1.22E-03,7.13E-04,4.33E-04,2.77E-04,1.51E-04, * 9.30E-05,4.53E-05/ DATA (FLX_NU_MIN( 1,3,LW),LW=1,23)/ * 4.61E-02,4.67E-02,4.93E-02,5.12E-02,5.02E-02,4.91E-02,4.68E-02, * 4.29E-02,3.79E-02,3.27E-02,2.67E-02,2.11E-02,1.64E-02,1.20E-02, * 8.62E-03,6.12E-03,4.11E-03,2.78E-03,1.83E-03,1.20E-03,8.08E-04, * 5.41E-04,3.33E-04/ DATA (FLX_NU_MIN( 2,3,LW),LW=1,23)/ * 3.98E-02,4.12E-02,4.39E-02,4.53E-02,4.54E-02,4.47E-02,4.27E-02, * 3.99E-02,3.55E-02,3.05E-02,2.58E-02,2.05E-02,1.61E-02,1.22E-02, * 8.85E-03,6.34E-03,4.23E-03,2.94E-03,2.02E-03,1.33E-03,8.48E-04, * 5.65E-04,3.68E-04/ DATA (FLX_NU_MIN( 3,3,LW),LW=1,23)/ * 4.10E-02,4.32E-02,4.59E-02,4.78E-02,4.78E-02,4.69E-02,4.49E-02, * 4.23E-02,3.77E-02,3.22E-02,2.71E-02,2.19E-02,1.72E-02,1.29E-02, * 9.43E-03,6.81E-03,4.81E-03,3.27E-03,2.19E-03,1.44E-03,9.71E-04, * 6.08E-04,4.14E-04/ DATA (FLX_NU_MIN( 4,3,LW),LW=1,23)/ * 3.50E-02,3.74E-02,4.01E-02,4.19E-02,4.27E-02,4.21E-02,4.06E-02, * 3.80E-02,3.43E-02,3.01E-02,2.54E-02,2.11E-02,1.67E-02,1.29E-02, * 9.56E-03,7.06E-03,5.07E-03,3.53E-03,2.37E-03,1.63E-03,1.09E-03, * 7.25E-04,4.82E-04/ DATA (FLX_NU_MIN( 5,3,LW),LW=1,23)/ * 3.23E-02,3.53E-02,3.83E-02,3.99E-02,4.06E-02,3.98E-02,3.85E-02, * 3.61E-02,3.29E-02,2.89E-02,2.45E-02,2.05E-02,1.65E-02,1.27E-02, * 9.38E-03,7.10E-03,5.09E-03,3.58E-03,2.51E-03,1.73E-03,1.18E-03, * 7.93E-04,5.24E-04/ DATA (FLX_NU_MIN( 6,3,LW),LW=1,23)/ * 3.02E-02,3.32E-02,3.59E-02,3.75E-02,3.76E-02,3.74E-02,3.61E-02, * 3.44E-02,3.12E-02,2.73E-02,2.34E-02,1.95E-02,1.58E-02,1.22E-02, * 9.34E-03,6.89E-03,5.04E-03,3.64E-03,2.61E-03,1.83E-03,1.26E-03, * 8.55E-04,5.77E-04/ DATA (FLX_NU_MIN( 7,3,LW),LW=1,23)/ * 2.77E-02,3.05E-02,3.31E-02,3.46E-02,3.48E-02,3.47E-02,3.36E-02, * 3.21E-02,2.93E-02,2.57E-02,2.22E-02,1.86E-02,1.51E-02,1.18E-02, * 9.08E-03,6.75E-03,4.96E-03,3.60E-03,2.59E-03,1.83E-03,1.26E-03, * 8.54E-04,5.77E-04/ DATA (FLX_NU_MIN( 8,3,LW),LW=1,23)/ * 3.55E-02,3.84E-02,4.13E-02,4.26E-02,4.24E-02,4.13E-02,3.94E-02, * 3.66E-02,3.29E-02,2.87E-02,2.41E-02,2.01E-02,1.61E-02,1.24E-02, * 9.19E-03,6.97E-03,5.01E-03,3.54E-03,2.49E-03,1.73E-03,1.18E-03, * 7.93E-04,5.24E-04/ DATA (FLX_NU_MIN( 9,3,LW),LW=1,23)/ * 7.38E-02,7.48E-02,7.57E-02,7.36E-02,6.86E-02,6.22E-02,5.56E-02, * 4.86E-02,4.14E-02,3.46E-02,2.81E-02,2.25E-02,1.74E-02,1.31E-02, * 9.57E-03,7.00E-03,5.00E-03,3.48E-03,2.35E-03,1.61E-03,1.08E-03, * 7.21E-04,4.80E-04/ DATA (FLX_NU_MIN(10,3,LW),LW=1,23)/ * 1.05E-01,9.88E-02,9.84E-02,9.50E-02,8.60E-02,7.63E-02,6.76E-02, * 5.91E-02,4.94E-02,4.01E-02,3.23E-02,2.52E-02,1.89E-02,1.37E-02, * 9.79E-03,6.93E-03,4.83E-03,3.26E-03,2.18E-03,1.43E-03,9.69E-04, * 6.07E-04,4.14E-04/ DATA (FLX_NU_MIN(11,3,LW),LW=1,23)/ * 1.42E-01,1.34E-01,1.33E-01,1.25E-01,1.12E-01,9.83E-02,8.46E-02, * 7.19E-02,5.85E-02,4.65E-02,3.66E-02,2.71E-02,2.00E-02,1.42E-02, * 9.81E-03,6.77E-03,4.39E-03,3.00E-03,2.03E-03,1.33E-03,8.50E-04, * 5.66E-04,3.68E-04/ DATA (FLX_NU_MIN(12,3,LW),LW=1,23)/ * 9.84E-02,9.52E-02,9.71E-02,9.52E-02,8.81E-02,8.09E-02,7.28E-02, * 6.33E-02,5.28E-02,4.32E-02,3.36E-02,2.55E-02,1.87E-02,1.32E-02, * 9.24E-03,6.34E-03,4.16E-03,2.79E-03,1.83E-03,1.20E-03,8.07E-04, * 5.41E-04,3.33E-04/ DATA (FLX_NU_MIN( 1,4,LW),LW=1,23)/ * 4.66E-02,4.71E-02,4.97E-02,5.09E-02,5.10E-02,4.91E-02,4.65E-02, * 4.29E-02,3.74E-02,3.17E-02,2.58E-02,2.08E-02,1.59E-02,1.14E-02, * 8.26E-03,5.67E-03,3.92E-03,2.52E-03,1.67E-03,1.04E-03,6.86E-04, * 4.50E-04,2.75E-04/ DATA (FLX_NU_MIN( 2,4,LW),LW=1,23)/ * 3.98E-02,4.15E-02,4.43E-02,4.56E-02,4.54E-02,4.48E-02,4.29E-02, * 3.94E-02,3.50E-02,3.05E-02,2.54E-02,2.01E-02,1.56E-02,1.16E-02, * 8.58E-03,5.81E-03,4.11E-03,2.76E-03,1.84E-03,1.19E-03,7.49E-04, * 4.81E-04,3.15E-04/ DATA (FLX_NU_MIN( 3,4,LW),LW=1,23)/ * 4.10E-02,4.31E-02,4.65E-02,4.77E-02,4.82E-02,4.74E-02,4.52E-02, * 4.19E-02,3.75E-02,3.23E-02,2.70E-02,2.16E-02,1.67E-02,1.27E-02, * 9.37E-03,6.61E-03,4.44E-03,3.14E-03,2.07E-03,1.35E-03,8.93E-04, * 5.45E-04,3.43E-04/ DATA (FLX_NU_MIN( 4,4,LW),LW=1,23)/ * 3.48E-02,3.74E-02,4.02E-02,4.19E-02,4.27E-02,4.20E-02,4.09E-02, * 3.82E-02,3.41E-02,3.03E-02,2.54E-02,2.08E-02,1.64E-02,1.27E-02, * 9.56E-03,7.06E-03,5.01E-03,3.38E-03,2.35E-03,1.59E-03,1.04E-03, * 6.86E-04,4.28E-04/ DATA (FLX_NU_MIN( 5,4,LW),LW=1,23)/ * 3.25E-02,3.56E-02,3.79E-02,4.02E-02,4.03E-02,3.99E-02,3.86E-02, * 3.63E-02,3.30E-02,2.88E-02,2.46E-02,2.01E-02,1.61E-02,1.27E-02, * 9.56E-03,7.06E-03,5.08E-03,3.60E-03,2.53E-03,1.69E-03,1.14E-03, * 7.66E-04,4.97E-04/ DATA (FLX_NU_MIN( 6,4,LW),LW=1,23)/ * 3.03E-02,3.32E-02,3.60E-02,3.76E-02,3.80E-02,3.77E-02,3.62E-02, * 3.38E-02,3.08E-02,2.74E-02,2.33E-02,1.95E-02,1.54E-02,1.22E-02, * 9.40E-03,7.00E-03,5.06E-03,3.66E-03,2.56E-03,1.84E-03,1.25E-03, * 8.74E-04,5.60E-04/ DATA (FLX_NU_MIN( 7,4,LW),LW=1,23)/ * 2.79E-02,3.06E-02,3.32E-02,3.47E-02,3.52E-02,3.50E-02,3.38E-02, * 3.16E-02,2.89E-02,2.58E-02,2.21E-02,1.86E-02,1.47E-02,1.18E-02, * 9.13E-03,6.85E-03,4.98E-03,3.62E-03,2.54E-03,1.83E-03,1.25E-03, * 8.73E-04,5.60E-04/ DATA (FLX_NU_MIN( 8,4,LW),LW=1,23)/ * 3.57E-02,3.87E-02,4.08E-02,4.28E-02,4.23E-02,4.14E-02,3.94E-02, * 3.67E-02,3.31E-02,2.87E-02,2.43E-02,1.98E-02,1.58E-02,1.24E-02, * 9.36E-03,6.92E-03,5.00E-03,3.56E-03,2.51E-03,1.69E-03,1.14E-03, * 7.65E-04,4.97E-04/ DATA (FLX_NU_MIN( 9,4,LW),LW=1,23)/ * 7.36E-02,7.44E-02,7.53E-02,7.36E-02,6.83E-02,6.22E-02,5.61E-02, * 4.91E-02,4.16E-02,3.51E-02,2.84E-02,2.24E-02,1.72E-02,1.31E-02, * 9.64E-03,7.02E-03,4.95E-03,3.35E-03,2.32E-03,1.57E-03,1.04E-03, * 6.83E-04,4.28E-04/ DATA (FLX_NU_MIN(10,4,LW),LW=1,23)/ * 1.05E-01,9.94E-02,9.89E-02,9.45E-02,8.66E-02,7.74E-02,6.81E-02, * 5.87E-02,4.98E-02,4.03E-02,3.23E-02,2.48E-02,1.84E-02,1.36E-02, * 9.80E-03,6.75E-03,4.49E-03,3.15E-03,2.07E-03,1.35E-03,8.90E-04, * 5.45E-04,3.43E-04/ DATA (FLX_NU_MIN(11,4,LW),LW=1,23)/ * 1.43E-01,1.34E-01,1.32E-01,1.25E-01,1.11E-01,9.77E-02,8.46E-02, * 7.15E-02,5.84E-02,4.62E-02,3.59E-02,2.66E-02,1.94E-02,1.37E-02, * 9.65E-03,6.25E-03,4.28E-03,2.82E-03,1.85E-03,1.19E-03,7.49E-04, * 4.81E-04,3.15E-04/ DATA (FLX_NU_MIN(12,4,LW),LW=1,23)/ * 9.89E-02,9.53E-02,9.70E-02,9.51E-02,8.82E-02,8.08E-02,7.18E-02, * 6.30E-02,5.18E-02,4.19E-02,3.24E-02,2.49E-02,1.83E-02,1.26E-02, * 8.83E-03,5.90E-03,4.00E-03,2.55E-03,1.68E-03,1.04E-03,6.86E-04, * 4.50E-04,2.75E-04/ END ************************************************************************* function QUAD_INT (R,X0,X1,X2,V0,V1,V2) ************************************************************************* implicit real*8 (a-h,o-z) R0 = R -X0 R1 = R -X1 R2 = R -X2 S0 = X0 -X1 S1 = X0 -X2 S2 = X1 -X2 QUAD_INT = V0*R1*R2/S0/S1 - * V1*R0*R2/S0/S2 + * V2*R0*R1/S1/S2 RETURN END ************************************************************************** subroutine asym_plot(dm2_test, s2th_test) ************************************************************************** * * make up/dn asymmetry plot versus momentum for data, mc and ms w/osc * implicit real*8 (a-h,o-z) parameter nm = 10000 common /data/run(nm),event(nm),x(nm),y(nm),z(nm), * dx(nm),dy(nm),dz(nm),evis(nm),enuest(nm),adpid(nm),rwall(nm), * anisx(nm),anisy(nm),anisz(nm),yrmthday(nm),hrminsec(nm), * pass(nm), eloed(nm), ringer(nm), dinter(nm), nevts parameter mm = 250000 common /mcdata/srun(mm),sevent(mm),sx(mm),sy(mm),sz(mm), * sdx(mm),sdy(mm),sdz(mm),sevis(mm), * senuest(mm),sadpid(mm),srwall(mm), * sanisx(mm),sanisy(mm),sanisz(mm),senergy(mm),sdzn(mm), * spass(mm), rfmuel(mm), eloem(mm), eloemt(mm), * sringer(mm), sinter(mm), mcevts c common /cuts/adpidmin,adpidmax,evismin,evismax,rwallmin, c + zmin,zmax,ringermin,ringermax dimension upe_dat(10), dne_dat(10), ase_dat(10) dimension upe_mcn(10), dne_mcn(10), ase_mcn(10) dimension upe_mco(10), dne_mco(10), ase_mco(10) dimension upm_dat(10), dnm_dat(10), asm_dat(10) dimension upm_mcn(10), dnm_mcn(10), asm_mcn(10) dimension upm_mco(10), dnm_mco(10), asm_mco(10) dimension ase_date(10) dimension ase_mcne(10) dimension ase_mcoe(10) dimension asm_date(10) dimension asm_mcne(10) dimension asm_mcoe(10) dimension ebin(10), flebin(10) do ia = 1,10 upe_dat(ia) = 0.0 dne_dat(ia) = 0.0 ase_dat(ia) = 0.0 upe_mcn(ia) = 0.0 dne_mcn(ia) = 0.0 ase_mcn(ia) = 0.0 upe_mco(ia) = 0.0 dne_mco(ia) = 0.0 ase_mco(ia) = 0.0 upm_dat(ia) = 0.0 dnm_dat(ia) = 0.0 asm_dat(ia) = 0.0 upm_mcn(ia) = 0.0 dnm_mcn(ia) = 0.0 asm_mcn(ia) = 0.0 upm_mco(ia) = 0.0 dnm_mco(ia) = 0.0 asm_mco(ia) = 0.0 ase_date(ia) = 0.0 ase_mcne(ia) = 0.0 ase_mcoe(ia) = 0.0 asm_date(ia) = 0.0 asm_mcne(ia) = 0.0 asm_mcoe(ia) = 0.0 ebin(ia) = 0.0 flebin(ia) = 0.0 end do * make asym histogram for data nused = 0 ndels = 0 ndmus = 0 do ievt = 1, nevts if(pass(ievt) .gt. 0.0) then * binning energy in 6 log steps from 100 MeV to 10 GeV if(enuest(ievt) .gt. 0.0) then alpha = dlog10(enuest(ievt)) else alpha = -10.0 end if ia = 1 + int((alpha-2.0)*3.0) if(ia .ge. 1 .and. ia .le. 6) then nused = nused + 1 * now bin em by up-dn and mu-e if(dz(ievt) .lt. 0.0) then if(adpid(ievt) .lt. 0.0) then ndmus = ndmus + 1 dnm_dat(ia) = dnm_dat(ia) + 1. else ndels = ndels + 1 dne_dat(ia) = dne_dat(ia) + 1. end if else if(adpid(ievt) .lt. 0.0) then ndmus = ndmus + 1 upm_dat(ia) = upm_dat(ia) + 1. else ndels = ndels + 1 upe_dat(ia) = upe_dat(ia) + 1. end if end if end if end if end do * make asym histogram for mc without osc mused = 0 nmels = 0 nmmus = 0 do ievt = 1, mcevts if(spass(ievt) .gt. 0.0) then * binning energy in 6 log steps from 100 MeV to 10 GeV if(senuest(ievt) .gt. 0.0) then alpha = dlog10(senuest(ievt)) else alpha = -10.0 end if ia = 1 + int((alpha-2.0)*3.0) if(ia .ge. 1 .and. ia .le. 6) then mused = mused + 1 * now bin em by up-dn and mu-e if(sdz(ievt) .lt. 0.0) then if(sadpid(ievt) .lt. 0.0) then nmmus = nmmus + 1 dnm_mcn(ia) = dnm_mcn(ia) + 1. else nmels = nmels + 1 dne_mcn(ia) = dne_mcn(ia) + 1. end if else if(sadpid(ievt) .lt. 0.0) then nmmus = nmmus + 1 upm_mcn(ia) = upm_mcn(ia) + 1. else nmels = nmels + 1 upe_mcn(ia) = upe_mcn(ia) + 1. end if end if end if end if end do * make asym histogram for fake osc signal in MC feffels = 0.0 feffmus = 0.0 do ievt = 1, mcevts if(abs(sinter(ievt)) .lt. 2.1) then ntyp = 1 ! for electrons cc and nc else ntyp = 2 ! for muons, ditto end if phi = eloem(ievt) rmuel = rfmuel(ievt) phit = eloemt(ievt) call posc(phit,ntyp,dm2_test,s2th_test,rmuel,plep) if(spass(ievt) .gt. 0.0) then * binning energy in 6 log steps from 100 MeV to 10 GeV if(senuest(ievt) .gt. 0.0) then alpha = dlog10(senuest(ievt)) else alpha = -10.0 end if ia = 1 + int((alpha-2.0)*3.0) if(ia .ge. 1 .and. ia .le. 6) then mused = mused + 1 * now bin em by up-dn and mu-e if(sdz(ievt) .lt. 0.0) then if(sadpid(ievt) .lt. 0.0) then feffmus = feffmus + plep dnm_mco(ia) = dnm_mco(ia) + plep else feffels = feffels + plep dne_mco(ia) = dne_mco(ia) + plep end if else if(sadpid(ievt) .lt. 0.0) then feffmus = feffmus + plep upm_mco(ia) = upm_mco(ia) + plep else feffels = feffels + plep upe_mco(ia) = upe_mco(ia) + plep end if end if end if end if end do * now calculate asymmetries, errors and print out close(6) open(unit=6, file='asym.prt', type='new') * write(6,1900) * 1900 format(' output of asym plots'/) do ia = 1,6 * put the energy in the middle of the bin (the +0.5 below) flebin(ia) = 2.0 + (float(ia)-1.0+0.5)/3.0 ebin(ia) = 10**flebin(ia) !MeV en_up = upe_dat(ia) en_dn = dne_dat(ia) en_tot = en_up + en_dn en_dif = en_dn - en_up if(en_tot .gt. 0) then ase_dat(ia) = en_dif/en_tot ase_date(ia) = 2.0*sqrt(en_up*en_dn/en_tot**3) end if en_up = upm_dat(ia) en_dn = dnm_dat(ia) en_tot = en_up + en_dn en_dif = en_dn - en_up if(en_tot .gt. 0) then asm_dat(ia) = en_dif/en_tot asm_date(ia) = 2.0*sqrt(en_up*en_dn/en_tot**3) end if en_up = upe_mcn(ia) en_dn = dne_mcn(ia) en_tot = en_up + en_dn en_dif = en_dn - en_up if(en_tot .gt. 0) then ase_mcn(ia) = en_dif/en_tot ase_mcne(ia) = 2.0*sqrt(en_up*en_dn/en_tot**3) end if en_up = upm_mcn(ia) en_dn = dnm_mcn(ia) en_tot = en_up + en_dn en_dif = en_dn - en_up if(en_tot .gt. 0) then asm_mcn(ia) = en_dif/en_tot asm_mcne(ia) = 2.0*sqrt(en_up*en_dn/en_tot**3) end if en_up = upe_mco(ia) en_dn = dne_mco(ia) en_tot = en_up + en_dn en_dif = en_dn - en_up if(en_tot .gt. 0) then ase_mco(ia) = en_dif/en_tot ase_mcoe(ia) = 2.0*sqrt(en_up*en_dn/en_tot**3) end if en_up = upm_mco(ia) en_dn = dnm_mco(ia) en_tot = en_up + en_dn en_dif = en_dn - en_up if(en_tot .gt. 0) then asm_mco(ia) = en_dif/en_tot asm_mcoe(ia) = 2.0*sqrt(en_up*en_dn/en_tot**3) end if write(6,2000) ebin(ia), flebin(ia), + ase_dat(ia), ase_date(ia), + asm_dat(ia), asm_date(ia), + ase_mcn(ia), ase_mcne(ia), + asm_mcn(ia), asm_mcne(ia), + ase_mco(ia), ase_mcoe(ia), + asm_mco(ia), asm_mcoe(ia) 2000 format(14e11.4) end do * calc chisquare value chisq_fit = 0.0 do ia = 3,6 chisq_fit = chisq_fit + ((asm_dat(ia)-asm_mco(ia))/ + (asm_date(ia)+asm_mcoe(ia)))**2 end do type *, ' ' type *,' chisq of asym fit data and model', chisq_fit close(6) return end