NotebookEvaluate["/home/max/projects/mathematica/style.nb"]

ϵπ=20*10^9;
ϵ=87*10^6;
fracEnergyInHadrComp[lgE0_, β_,ecp_]:=(((10^lgE0)/ecp)^β )*(ecp/10^lgE0)

ContourPlot[fracEnergyInHadrComp[e, b,ϵπ],{e,15,21},{b,0.8,0.99},FrameLabel->{T["lg((Subscript[E, 0])/eV)"],T["β"]}, Evaluate@customStyle,PlotLegends->Automatic]

fracEnergyInHadrComp[18, 0.95]

Nμ[lgE0_,β_,A_,ecp_]:=((10^lgE0/ecp/A)^β A)
Neγ[lgE0_,β_,A_,ecp_,ec_]:=(10^lgE0-ecp(10^lgE0/ecp)^β*A^(1-β))/ec
DNeγ[lge_,b_,a_,ecp_,ec_]:=(((a^(-b)) (-1+b) (((((10^lge)/ecp)))^b) ecp)/ec)
FullSimplify[Neγ[Log10[e0],b,a,ecp,ec]]
FullSimplify[D[Neγ[Log10[e0],b,a,ecp,ec],a]/Neγ[Log10[e0],b,1,ecp,ec]]

LogPlot[
{Neγ[16, 0.96,a,ϵπ,ϵ],
Neγ[17, 0.96,a,ϵπ,ϵ],
Neγ[18, 0.96,a,ϵπ,ϵ],
Neγ[19, 0.96,a,ϵπ,ϵ]
},{a,1,56},
 PlotRange->All,Evaluate@customStyle,PlotLabels->Table[T["lg((Subscript[E, 0])/eV)= "<>ToString[e]],{e,{16,17,18,19}}],
FrameLabel->{T["A"],T["((Subscript[N, eγ])^((A)))"]}
]

Plot[
{DNeγ[16, 0.96,a,ϵπ,ϵ],
DNeγ[17, 0.96,a,ϵπ,ϵ],
DNeγ[18, 0.96,a,ϵπ,ϵ],
DNeγ[19, 0.96,a,ϵπ,ϵ]
}
,{a,1,56},PlotRange->All,Evaluate@customStyle,PlotLabels->Table[T["lg((Subscript[E, 0])/eV)= "<>ToString[e]],{e,{16,17,18,19}}],
FrameLabel->{T["A"],T["(∂((Subscript[N, eγ])^((A)))/∂A)/ %"]}
]

Plot[
{DNeγ[16, 0.96,a,ϵπ,ϵ]/Neγ[16, 0.96,1,ϵπ,ϵ]100,
DNeγ[17, 0.96,a,ϵπ,ϵ]/Neγ[17, 0.96,1,ϵπ,ϵ]100,
DNeγ[18, 0.96,a,ϵπ,ϵ]/Neγ[18, 0.96,1,ϵπ,ϵ]100,
DNeγ[19, 0.96,a,ϵπ,ϵ]/Neγ[19, 0.96,1,ϵπ,ϵ]100,
DNeγ[20, 0.96,a,ϵπ,ϵ]/Neγ[20, 0.96,1,ϵπ,ϵ]100
}
,{a,1,56},PlotRange->All,Evaluate@customStyle,PlotLabels->Table[T["lg((Subscript[E, 0])/eV)= "<>ToString[e]],{e,{16,17,18,19,20}}],
FrameLabel->{T["A"],T["1/((Subscript[N, eγ])^((1)))(∂((Subscript[N, eγ])^((A)))/∂A)/ %"]}
]

s[X_,Xmax_]:=3X/(X+2Xmax)
GreissenProfile[X_,Xmax_,X0_]:=0.31/Sqrt[Xmax/X0] Exp[X/X0(1-3/2Log[s[X,Xmax]])]
GaisserHillasProfile[X_,Xmax_,λ_,Nmax_,X1_]:=Nmax((X-X1)/(Xmax-X1))^((Xmax-X1)/λ) Exp[-(X-Xmax)/(λ)]

Plot[{GreissenProfile[x,700,37],
GaisserHillasProfile[x,700,1.5*37,0.31Exp[700/37]/Sqrt[700/37],1]},
{x,0,1200},PlotStyle->{Directive[Orange,Thickness[0.005]], Directive[Black,Thickness[0.005]]},Evaluate@customStyle,
PlotLegends->Placed[{T["Greisen"],T["Gaisser-Hillas"]},{0.25,0.8}],
FrameLabel->{T["X / (g (cm^(-2)))"],T["N"]}
]

(* why components? *)

Neγhad[lgE0_, beta_, ecp_,ec_,m_]:=(fracEnergyInHadrComp[lgE0, beta,ecp]*m*10^(lgE0))/ec

Plot[
{Neγhad[lge, 0.96,ϵπ,ϵ,0.15]/(Neγ[lge, 0.96,1,ϵπ,ϵ]+Neγhad[lge, 0.96,ϵπ,ϵ,0.15]) 100,
Neγhad[lge, 0.96,ϵπ,ϵ,0.1]/(Neγ[lge, 0.96,1,ϵπ,ϵ]+Neγhad[lge, 0.96,ϵπ,ϵ,0.1]) 100,
Neγhad[lge, 0.96,ϵπ,ϵ,0.05]/(Neγ[lge, 0.96,1,ϵπ,ϵ]+Neγhad[lge, 0.96,ϵπ,ϵ,0.05]) 100
},{lge,15,20.5},
 PlotRange->{0,25},Evaluate@customStyle,
PlotLabels->Table[T["m = "<>ToString[e]],{e,{0.15,0.10,0.05}}],
FrameLabel->{T["lg((Subscript[E, 0])/eV)"],T["(Subscript[N, eγhad])/(Subscript[N, eγtot]) in %"]}
]

Neγhad[Log10[e],be,ecp, ec,m]