#!/sw/bin/gnuplot -persist # # # G N U P L O T # Version 4.0 patchlevel 0 # last modified Thu Apr 15 14:44:22 CEST 2004 # System: Darwin 8.10.1 # # Copyright (C) 1986 - 1993, 1998, 2004 # Thomas Williams, Colin Kelley and many others # # This is gnuplot version 4.0. Please refer to the documentation # for command syntax changes. The old syntax will be accepted # throughout the 4.0 series, but all save files use the new syntax. # # Type `help` to access the on-line reference manual. # The gnuplot FAQ is available from # http://www.gnuplot.info/faq/ # # Send comments and requests for help to # # Send bugs, suggestions and mods to # # # set terminal x11 # set output set grid set title "Kin4: Ratio N_{obs} / P_{obs} for Nitrogen, Helium-3, and Hydrogen" 0.000000,0.000000 font "" set xlabel "V1L Rate (Hz)" 0.000000,0.000000 font "" set ylabel "Observed Neutrons / Observed Protons" 0.000000,0.000000 font "" set fit errorvariables rdiff(a,b)=(a-b)/a # correct raw scaler rates for deadtime rate(obs,dt)=obs/(1-obs*dt) # slightly imprecise form -- probability = rate*time rold(np,k1,k2,T1)=np*(1.-(k1*dtn+k2*dtn-k1*k2*dtn*dtn*T1)*T1)/(np*(k1*dtp+k2*dtp-k1*k2*dtp*dtp*T1)*T1+(1.-k1*k2*T1*T1*edt*edt)) po_old(np,k1,k2,T1)=np*(k1*dtp+k2*dtp-k1*k2*dtp*dtp*T1)*T1+(1.-k1*k2*T1*T1*edt*edt) no_old(np,k1,k2,T1)=np*(1.-(k1*dtn+k2*dtn-k1*k2*dtn*dtn*T1)*T1) # using Poisson distribution -- probability for 1 or more hits = 1-exp(-rate*time) prob(mu)=1-exp(-mu) po(np,k1,k2,T1)=np*(prob((k1+k2)*dtp*T1))+(1.-prob(k1*edt*T1)*prob(k2*edt*T1)) no(np,k1,k2,T1)=np*(1.-prob((k1+k2)*dtn*T1)) r(np,k1,k2,T1)=no(np,k1,k2,T1)/po(np,k1,k2,T1) # Now, instead solve in terms of ratio at arbitrary rate, to ease # the error analysis, here np = ratio at the desired rate f. # # this is just the expanded verion of 'r' above, where instead of 'np' being # at zero random rate, it is at the desired reference freq. 'refrate' rref(np,k1,k2,T1,f)=(np*(1-prob(k1*edt*f)*prob(k2*edt*f))*(1-prob((k1+k2)*dtn*T1))) / ( (1-prob(k1*edt*f)*prob(k2*edt*f))*np*prob((k1+k2)*dtp*T1) + (1-prob((k1+k2)*dtn*f)-np*prob((k1+k2)*dtp*f))*(1-prob(k1*edt*T1)*prob(k2*edt*T1))) # set 'refrate' to the point you want the ratio's (in Hz) refrate=80.e3 # helper functions to handle fitting multiple datasets simultaneously. # normal 'np' is at 0 f(x,y)=(y>=0&&y<1)*(r(nph,vt,vt,x))+(y>=1&&y<2)*(r(nphe,vt,vt,x))+(y>=2)*(r(npn,vt,vt,x)) fr(x,y)=(y>=0&&y<1)*(rref(nph,vt,vt,x,refrate))+(y>=1&&y<2)*(rref(nphe,vt,vt,x,refrate))+(y>=2)*(rref(npn,vt,vt,x,refrate)) # possible values for constants.... dtn = 1.2e-07 dtp = 2e-08 edt = 0.0 nph = 0.0442458448311674 vt = 20.3989970190169 nphe = 0.110223385094474 npn = 0.186028710162187 dt = 1e-08 np = 0.15 k1_he = 28.584 dk1_he = 6.455 np_h = 0.104885195945144 dnp_h = 0.0057 np_he = 0.18081135099027 dnp_he = 0.01915 np_n = 0.303814526024512 dnp_n = 0.0194 k1 = -3.23076433728644 dk1 = 6.455 r0 = 80000.0 np_h0 = 0.0422421881184299 v0 = 0.170123744463467 v1 = 0.180742772447311 v2 = 0.149472178311587 dnp_h0 = 0.00594170264448289 np_n0 = 0.170123744463467 dnp_n0 = 0.0232217772758063 nph_err = 0.00237470524248336 nphe_err = 0.000443538062105081 npn_err = 0.00838503746602879 vt_err = 3.40964406390994 fit fr(x,y) 'ylds2_kin4.txt' using ($14*1.e6):-2:9:10 via nph,nphe,npn,vt dn = (2*nph*npn-nph*nphe - nphe*npn)/(nphe*(nph-npn)) edn = sqrt( (2*npn*(nphe-npn)/(nphe*(nph-npn)**2)*nph_err)**2 + (2*nph*(nph-nphe)/(nphe*(nph-npn)**2)*npn_err)**2 + (2*nph*npn/(nphe**2 * (nph-npn))*nphe_err)**2 ) dp = 2*(npn-nphe)/(npn-nph) edp = sqrt( (2*(npn-nphe)/(npn-nph)**2 * nph_err)**2 + (2*(nphe-nph)/(npn-nph)**2 * npn_err)**2 + (2/(npn-nph) * nphe_err)**2 ) print "Dn = ",dn," +/- ",edn print "Dp = ",dp," +/- ",edp set label 20 "Evaluated at V1L Rate = %3.0g Hz :",refrate at 40000, .41 set label 21 "nph = %3.5g",nph ," +- %3.5g",nph_err at 40000, .38 set label 22 "nphe = %3.5g",nphe," +- %3.5g",nphe_err at 40000, .35 set label 23 "npn = %3.5g",npn ," +- %3.5g",npn_err at 40000, .32 set label 24 "slp = %3.5g",vt ," +- %3.5g",vt_err at 40000, .29 set label 25 "Dn = %3.5g",dn ," +- %3.5g",edn at 60000, .26 set label 26 "Dp = %3.5g",dp ," +- %3.5g",edp at 60000, .23 plot [0.:130000] 'ylds2_kin4.txt' index 0 using ($14*1.e6):9:10 t 'H_2 data' with yerrorbars lw 2, 'ylds2_kin4.txt' index 1 using ($14*1.e6):9:10 t '3He data' with yerrorbars lw 2, 'ylds2_kin4.txt' index 2 using ($14*1.e6):9:10 t 'BeO/C data' with yerrorbars lw 2., 'ylds2_kin4.txt' index 3 using ($14*1.e6):9:10 t 'N_2 data' with yerrorbars lw 2., fr(x,1.1) lw 2, fr(x,0.1) lw 2, fr(x,2.1) lw 2 ## fit fr(x,y) 'ylds2_kin4.txt' using ($14*1.e6):-2:9:10 via nph,nphe,npn,vt # EOF