Read in the LBNF neutrino beam spectrum.  And calculate the event rate approximately.  

Let’s learn a little about units of energy. 

1 electron charge is  1.60*10^19 Coulomb  

Recall:     Current * Volt  =  Power ==>     Charge * Voltage = Energy  

1 electroncharge * 1 volt =  1 eV of energy  =>    1.60*10^-19 Coulomb * 1 Volt  = 1.60*10^-19 Joule.  

In Eletrical engineering, physics,  astonoimy  etc.  eV  is used for energy.   


1 eV  
1000 eV = 1 keV  
10^6 eV  = 1 MeV  
10^9  eV = 1 GeV    

etc.

electronmass = UnitConvert[Quantity[1,"ElectronMass?"],"kg"] //N

electroncharge = UnitConvert[Quantity[1,"ElectronCharge"],"Coulomb"] //N

What is the speed of an electron that goes through 200 volts of potential difference. 
Recall volts  = Joules/Coulomb  And so the electron gets a kinetic energy of 200 Joules/Coulomb*electron charge. 
1/2

UnitConvert[Quantity[200,"Volts"],"Joules/Coulomb"]

electronkinenergy = UnitConvert[Quantity[200,"Volts"],"Joules/Coulomb"]* UnitConvert[Quantity[1,"ElectronCharge"],"Coulomb"] ; 
electronvelocity = UnitConvert[Sqrt[2*electronkinenergy/electronmass],"meters/sec"]

electronvelocity/Quantity[1,"SpeedOfLight"]

The electron is so light that even after 200 volts it is quite fast.

SetDirectory["/Users/diwan/Documents/talks/renew-lectures/units/unit6-neutrinos"]

<<nuosc`

Read the neutrino flux in. The first line is text.

neutrinoflux = OpenRead["histos_g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_neutrino_LBNEFD_globes_flux.txt"] ; 
firstline= Read[neutrinoflux, String]; 
fnudata= ReadList[neutrinoflux, {Real, Real, Real, Real, Real, Real, Real}] ; 
Close[neutrinoflux]; 
antineutrinoflux = OpenRead["histos_g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_antineutrino_LBNEFD_globes_flux.txt"]; 
firstaline= Read[antineutrinoflux, String]; 
fanudata= ReadList[antineutrinoflux, {Real, Real, Real, Real, Real, Real, Real}] ; 
Close[antineutrinoflux];

Let’s reproduce the oscillation curve from the slides.

LogLinearPlot[{nuosc[1, 1, 0, len, 1,2, 2], 1-nuosc[1,1,0,len, 1,2,3]},{len, 10., 4000}, Frame-> True, FrameStyle -> Bold, LabelStyle->{Bold, Large, Italic}, PlotRange -> {{10., 4000.},{0, 1.}}, PlotStyle -> {Thick, Bold}, FrameLabel-> {"Baseline Length (km)", "Probabilty for 1 GeV muon neutrino"}   ]

The units for the flux are neutrinos per GeV per meter^2 per proton on target.

TableForm[fnudata,TableHeadings->{None, {"E", "nue","numu","nutau","nuebar", "numubar", "nutaubar"}}]

TableForm[fanudata,TableHeadings->{None, {"E", "nue","numu","nutau","nuebar", "numubar", "nutaubar"}}]

epl  = ListLogPlot[Transpose[{Transpose[fnudata][[1]], Transpose[fnudata][[2]]}],PlotRange -> {{0, 40},All} ,Joined->  True, PlotStyle -> {Blue, Dashed} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Flux /Gev/m^2/POT"} ]

mupl  = ListLogPlot[Transpose[{Transpose[fnudata][[1]], Transpose[fnudata][[3]]}],PlotRange -> {{0, 40},All} ,Joined->  True, PlotStyle -> {Red, Bold} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Flux /Gev/m^2/POT"} ]

Show[mupl,epl]

Calculate the numbers of protons per year assuming 120 GeV and 1 MW operation for 2*10^7 seconds per year.  
10^6 watts/120 *10^9 volts / 1.6*10^-19 coulomb/proton     * 2*10^7 second

potperyear = 10^6/(120*10^9)/(1.6*10^-19)  * 2*10^7

What is the cross section for neutrinos and antineutrinos per nucleon in the target. Units are cm^2/GeV and enu is in GeV.  
this is a first order approximation.   In practice, the cross section turns on and has a big rise at low energies before it is linear in energy.  This is good enough above 1 GeV.

signu[enu_]:= 0.8*enu*10^-38 ;  siganu[enu_]:= 0.35*enu*10^-38 ; 
signuqe[enu_]:= Piecewise[{{0.8*enu*10^-38, enu< 1}, {0.8*10^-38, enu≥ 1}}];   (* this is not really correct but close enough *)  
siganuqe[enu_]:= Piecewise[{{0.35*enu*10^-38, enu< 1}, {0.35*10^-38, enu≥ 1}}];

Plot[{signu[enu],siganu[enu], signuqe[enu],siganuqe[enu]},{enu, 0, 5}]

How many nucleons are there in a kiloton of liquid argon or any other detector.

targetsperkt = 6.022*10^23 * 1000*1000*1000 ;    (* 1000 tons * 1000 kg * 1000 gm *)  
totalmass = 40;  (* 40 ktons *)  
totalyrs = 5;

ebins = Transpose[fnudata][[1]]; 
binsize = ebins[[2]]-ebins[[1]]; 
muflux = Transpose[fnudata][[3]];  (* flux is per meter =^2 *) 
eflux =  Transpose[fnudata][[2]];
amuflux = Transpose[fanudata][[6]];  (* flux is per meter =^2 *) 
aeflux =  Transpose[fanudata][[5]];

mueventspec = Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
eeventspec =  Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*eflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
amueventspec = Table[totalmass*totalyrs*targetsperkt*potperyear*siganu[ebins[[i]]]*amuflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
aeeventspec =  Table[totalmass*totalyrs*targetsperkt*potperyear*siganu[ebins[[i]]]*aeflux[[i]]*10^-4,  {i,Length[fnudata]}]; 

mueventspecqe = Table[totalmass*totalyrs*targetsperkt*potperyear*signuqe[ebins[[i]]]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
eeventspecqe =  Table[totalmass*totalyrs*targetsperkt*potperyear*signuqe[ebins[[i]]]*eflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
amueventspecqe = Table[totalmass*totalyrs*targetsperkt*potperyear*siganuqe[ebins[[i]]]*amuflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
aeeventspecqe =  Table[totalmass*totalyrs*targetsperkt*potperyear*siganuqe[ebins[[i]]]*aeflux[[i]]*10^-4,  {i,Length[fnudata]}]; 

totalmus = Total[mueventspec]*binsize; 
totales = Total[eeventspec]*binsize; 
totalamus = Total[amueventspec]*binsize; 
totalaes = Total[aeeventspec]*binsize; 

totalmusqe = Total[mueventspecqe]*binsize; 
totalesqe = Total[eeventspecqe]*binsize; 
totalamusqe = Total[amueventspecqe]*binsize; 
totalaesqe = Total[aeeventspecqe]*binsize;

{totalmus,totales,totalamus,totalaes}

{totalmusqe,totalesqe,totalamusqe,totalaesqe}

Plot the number of expected events for muon neutrinos.

pl1  =ListPlot[Transpose[{ebins,mueventspec}],PlotRange -> {{0, 20},All} ,Joined->  True, PlotStyle -> {Blue, Bold} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Events/Gev/m^2/POT"} ]

Apply disappearance oscillation probability to the muon neutrino expectation.

oscmueventspec = Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*nuosc[ebins[[i]],1, 0, 1300, 1,2,2]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
pl2 = ListPlot[Transpose[{ebins,oscmueventspec}],PlotRange -> {{0, 10},All} ,Joined->  True, PlotStyle -> {Red, Bold} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Events/Gev/m^2/POT"} ]

oscmueventspecimh = Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*nuosc[ebins[[i]],-1, 0, 1300, 1,2,2]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
pl2imh = ListPlot[Transpose[{ebins,oscmueventspecimh}],PlotRange -> {{0, 10},All} ,Joined->  True, PlotStyle -> {Black, Bold} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Events/Gev/m^2/POT"} ]

Show[pl2,pl2imh]

?nuosc

Show[pl1,pl2]

osceeventspecm90 = Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*nuosc[ebins[[i]],1,-90,1300, 1,2,1]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
pl3= ListPlot[Transpose[{ebins,osceeventspecm90}],PlotRange -> {{0, 10},All} ,Joined->  True, PlotStyle -> {Green, Bold} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Events/Gev/m^2/POT"} ]

osceeventspecp90 = Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*nuosc[ebins[[i]],1,90,1300, 1,2,1]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]; 
pl4= ListPlot[Transpose[{ebins,osceeventspecp90}],PlotRange -> {{0, 10},All} ,Joined->  True, PlotStyle -> {Magenta, Bold} ,Frame -> True, LabelStyle -> {Bold, Large,Italic},  FrameStyle -> Bold,  FrameLabel-> {"Energy/GeV", "Events/Gev/m^2/POT"} ]

Re[{Total[mueventspec]*0.25, Total[oscmueventspec]*0.25 ,  Total[osceeventspecm90]*0.25, Total[osceeventspecp90]*0.25} ]

{44821.75055663915`,17171.973777125288`,1867.8294449386408`,2798.1134211836647`}

Show[pl4,pl3]

oscecprangenmh = Table[Total[Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*Re[nuosc[ebins[[i]],1,cp,1300, 1,2,1]]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]]*binsize, {cp, 0, 360, 30} ]

oscecprangeimh = Table[Total[Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*Re[nuosc[ebins[[i]],-1,cp,1300, 1,2,1]]*muflux[[i]]*10^-4,  {i,Length[fnudata]}]]*binsize, {cp, 0, 360, 30} ]

oscaecprangenmh = Table[Total[Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*Re[nuosc[ebins[[i]],1,cp,1300, -1,2,1]]*amuflux[[i]]*10^-4,  {i,Length[fnudata]}]]*binsize, {cp, 0, 360, 30} ]

oscaecprangeimh = Table[Total[Table[totalmass*totalyrs*targetsperkt*potperyear*signu[ebins[[i]]]*Re[nuosc[ebins[[i]],-1,cp,1300, -1,2,1]]*amuflux[[i]]*10^-4,  {i,Length[fnudata]}]]*binsize, {cp, 0, 360, 30} ]

ListPlot[{Transpose[{oscecprangenmh,oscaecprangenmh}],Transpose[{oscecprangeimh,oscaecprangeimh}]}, Frame -> True, FrameStyle -> Bold, LabelStyle -> {Bold, Large, Italic}, Joined -> True, PlotStyle-> Thick, FrameLabel->{"muon neutrino -> electron neutrino  events", "muon antineutrino -> electron antineutrino events"},PlotRange->{{0, 3000},{0, 3000}}]