macro mkbg * * make background spectrum of T2K nue appearance using SK atmospheric * neutrino MC. * * ikam=1 : background at kamioka (for osc. calc) * 2 : background at korea * mess 'Make sure you change the kumac for the right parameter' exec init deg=25 exec makehist ikam=2 deg=25 return ********************************************************************** macro init deg=1 mess deg [deg] close 0 h/del * * open SK-I DST *exec /home/atmpd/paw/dstopen.kumac *uwfunc //fcmc/1 ana.inc exec ./dstopen.kumac uwfunc //fcmc/1 myana.inc close 0 * import T2K flux histogram * for neutrino run, numu and nue * for anti-neutrino run, numu-bar and nue-bar h/file 1 ../prob2event_4gev/fluxkorea[deg]_hist.hbk hrin 10,30 0 10000 | numu and nue close 1 h/file 1 ../prob2event_4gev/fluxkorea[deg]_hist.hbk hrin 10,30 0 11000 | numu-bar and nue-bar close 1 * import SK atmospheric neutrino flux histogram h/file 1 flux.hbk hrin 101,102,103,104 0 20000 | nue,nueb,numu,numub close 1 * make (T2K)/(ATMNU) flux ratio histogram div 10010 20103 30001 | numu div 11010 20104 30002 | numu-bar div 10030 20101 30003 | nue div 11030 20102 30004 | nue-bar h/file 3 ratio.hbk ! N hrout 30001 hrout 30002 hrout 30003 hrout 30004 close 3 return ********************************************************************** macro makehist ikam=1 deg=2 mess ikam [ikam] deg [deg] xmin=0.0 xmax=4000.0 nbin=400 h/cre/1d 11 'nu BG total (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 15 'numu BG CC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 16 'numu BG NC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 17 'nue BG CC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 18 'nue BG NC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 21 'nu-bar BG total (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 25 'numub BG CC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 26 'numub BG NC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 27 'nueb BG CC (Erecon)' [nbin] [xmin] [xmax] h/cre/1d 28 'nueb BG NC (Erecon)' [nbin] [xmin] [xmax] *Copy ratio histograms to vectors vec/cr v1(401) vec/cr v2(401) vec/cr v3(401) vec/cr v4(401) h/get/con 30001 v1 h/get/con 30002 v2 h/get/con 30003 v3 h/get/con 30004 v4 Appli comis quit include 'enue.f' function fillhist(deg,ikam) include 'ana.inc' logical lfcfv, lelike, llike logical lcc integer deg,ikam real sin22th23, dm integer idoff real like_weight real nue1,nue2,nue3,nue4 real numu1,numu2,numu3,numu4 real nc1,nc2,nc3,nc4 c vector v1 c vector v2 c vector v3 c vector v4 parameter ( sin22th23=1.0 ) parameter ( dm=2.5E-3 ) C FCFV cut c write(*,*)'deg=',deg,' ikam=',ikam, ' mode',mode lfcfv = (wall.gt.200.).and.(evis.gt.30.).and.(nhitac.le.9) C 1-ring e-like no-decaye cut lelike = (nring.eq.1).and.(ip(1).eq.2).and.(ndcy.eq.0) C I will use a weight for the likelihood: C update for later make different table from different OA. C 1.0OA: C Likelihood eff for nue: 84.4% 80.9% 77.3% 73.5% C Likelihood eff for numu:25.7% 36.1% 32.4% 13.3% C Likelihood eff for nc: 8.9% 22.1% 27.4% 38.5% C 1.5OA: C Likelihood eff for nue: 84.6% 81.2% 78.1% 73.2% C Likelihood eff for numu:22.9% 31.2% 30.6% 13.5% C Likelihood eff for nc: 10.5% 22.1% 28.5% 38.3% C 2.0OA: C Likelihood eff for nue: 85.1% 81.4% 78.4% 73.3% C Likelihood eff for numu:17.4% 30.5% 30.4% 14.3% C Likelihood eff for nc: 10.8% 21.9% 27.4% 38.0% C 2.5OA: C Likelihood eff for nue: 85.3% 81.7% 78.4% 73.3% C Likelihood eff for numu:15.2% 30.4% 31.2% 15.7% C Likelihood eff for nc: 10.5% 21.5% 27.4% 37.7% if (deg.eq.1) then nue1=0.844 nue2=0.809 nue3=0.773 nue4=0.735 numu1=0.257 numu2=0.361 numu3=0.324 numu4=0.133 nc1=0.089 nc2=0.221 nc3=0.274 nc4=0.385 else if (deg.eq.15) then nue1=0.846 nue2=0.812 nue3=0.781 nue4=0.732 numu1=0.229 numu2=0.312 numu3=0.306 numu4=0.135 nc1=0.105 nc2=0.221 nc3=0.285 nc4=0.385 else if (deg.eq.2) then nue1=0.851 nue2=0.817 nue3=0.784 nue4=0.733 numu1=0.174 numu2=0.305 numu3=0.304 numu4=0.143 nc1=0.108 nc2=0.219 nc3=0.274 nc4=0.380 else if (deg.eq.25) then nue1=0.853 nue2=0.817 nue3=0.784 nue4=0.733 numu1=0.152 numu2=0.304 numu3=0.312 numu4=0.157 nc1=0.105 nc2=0.225 nc3=0.274 nc4=0.377 endif if ((abs(ipnu(1)).eq.12).and.(abs(mode).lt.30)) then if (pnu(1).lt.0.35) then like_weight=nue1 else if (pnu(1).lt.0.85) then like_weight=nue2 else if (pnu(1).lt.1.5) then like_weight=nue3 else like_weight=nue4 endif else if ((abs(ipnu(1)).eq.14).and.(abs(mode).lt.30)) then if (pnu(1).lt.0.35) then like_weight=numu1 else if (pnu(1).lt.0.85) then like_weight=numu2 else if (pnu(1).lt.1.5) then like_weight=numu3 else like_weight=numu4 endif else if (abs(mode).gt.30) then if (pnu(1).lt.0.35) then like_weight=nc1 else if (pnu(1).lt.0.85) then like_weight=nc2 else if (pnu(1).lt.1.5) then like_weight=nc3 else like_weight=nc4 endif endif C likelihood cut llike = .true. C reject if ( .not. (lfcfv.and.lelike.and.llike )) goto 999 C check CC or NC if ( abs(mode).ge.30 ) then lcc=.false. else lcc=.true. endif C check neutrino or anti-neutrino if ( ipnu(1).gt.0 ) then lanti=.false. idoff=10 else lanti=.true. idoff=20 endif C numu oscillation if ( abs(ipnu(1)).eq.14 .and. lcc ) then if ( ikam.eq.1 ) then xlen=295. else xlen=1050. endif posc = 1.0-sin22th23 * sin(1.27*dm*xlen/pnu(1))**2. else posc = 1.0 endif C energy etrue = pnu(1)*1000. ereco = enue(1)*1000. C flux ratio ratio = flux_ratio(ipnu(1),etrue) C total histogram weight c Using the nomralization factor I computed myself: weight=ratio*posc*like_weight*0.00693 c print *,ratio C if ( abs(ipnu(1)).eq.14 .and. lcc ) then C print *,ratio,posc,weight,ereco C endif C numu if ( abs(ipnu(1)).eq.14 ) then if ( lcc ) then call hf1(idoff+5, ereco, weight) else call hf1(idoff+6, ereco, weight) endif C nue else if ( abs(ipnu(1)).eq.12 ) then if ( lcc ) then call hf1(idoff+7, ereco, weight) else call hf1(idoff+8, ereco, weight) endif endif 999 continue fillhist=1.0 return end C********************************************************************** C********************************************************************** real function flux_ratio(ip,etrue) integer hid, ip real etrue, egev,hx,test integer bin vector v1 vector v2 vector v3 vector v4 if ( etrue.gt.4000. ) then egev=4. else egev = etrue/1000. endif bin=int(egev*100)+1 if ( ip.eq.14 ) then hid=30001 flux_ratio=v1(bin) else if ( ip.eq.-14 ) then hid=30002 flux_ratio=v2(bin) else if ( ip.eq.12 ) then hid=30003 flux_ratio=v3(bin) else if ( ip.eq.-12 ) then hid=30004 flux_ratio=v4(bin) endif *Was used before wne running on solaris: c flux_ratio = hx(hid,egev) c print *,flux_ratio, egev, hid return end quit h/file 2 ratio.hbk nt/loop //fcmc/1 fillhist([deg],[ikam]) close 2 add 15 11 11 1.0 1.0 add 16 11 11 1.0 1.0 add 17 11 11 1.0 1.0 add 18 11 11 1.0 1.0 add 25 21 21 1.0 1.0 add 26 21 21 1.0 1.0 add 27 21 21 1.0 1.0 add 28 21 21 1.0 1.0 return