Quit[];

(* Set Directory path as the same where the Notebook is saved *)

SetDirectory[NotebookDirectory[]];

(* Load FeynCalc package *)

<<FeynCalc`

Useful definitions

Plot Style Options

MyFont="Times";
MyLabelSize=20;
MyTickSize=16;
MyLegendSize=14;

MyOptions=Sequence[Axes->None,Frame->True,GridLines->Automatic,GridLinesStyle->Lighter[Lighter[Gray]],LabelStyle->Directive[FontFamily->MyFont,FontSize->MyLabelSize],FrameStyle->Directive[FontFamily->MyFont,FontSize->MyLabelSize],FrameTicksStyle->Directive[FontFamily->MyFont,FontSize->MyTickSize],ImageSize->450];

Constants

αem = 1/137; (* QED coupling constant *)

(* Write everything in GeV *)
TeV=(10^3);MeV=(10^(-3));keV = (10^(-6));eV=(10^(-9));

(* barn converted to (GeV^(-2))*)
barn=2.56819 (10^(-3)) (MeV^(-2)); pico=(10^(-12));pb = pico barn;

hbarc = 197.327 (10^(-18)); (* in GeV*m *)

(* Belle II *)

γβ B2 = 0.28;

sB2 = (10.6^2); (* COM energy *)
LumB2 = 445; (* in (pb^(-1)) *)
LumB2 final = 50 (10^3); (* 50 (ab^(-1)) *)

Lmin = 0.1; Lmax = 3;(* Lengths are in meters *)

(* Angular acceptance (in the lab frame) *)
θmin = 37.3 Degree;θmax = 123 Degree;

(* Minimum angular separation of two photons *)
Δθ γ = 0.014; (* rad*)

(* Boost angle along propagation direction*)

LHE Analysis

These functions are useful to read .lhe (typical output of MadGraph) files in Mathematica

(*  Based on the  Chameleon  package  by  Philip  Shuster , Jesse  Thaler , \Natalia  Toro *)(*  Modified  by  Maxim  Perelstein  and  Andi  Weiler , \2008 -09 *)

(*<<Histograms `*)


(* format  of the  compulsory  event  vector *)eventInfo = {_, _, _, _, _, _};ClearAll[particle , finout ];(*PDG  numbers  see e.g. \http :// home.fnal.gov/~ maeshima/alignment/ORCA/PYTHIA_particle_codes .\ps*)particle [2] = "u";particle [-2] = "uBar";particle [1] =" d";
particle [-1] = "dBar";particle [3] = "s";particle [-3] = "sBar";particle [4] = "c";particle [-4] = "cBar";particle [5] = "b";particle [-5] = "bBar";particle [6] = "t";particle [-6] = "tBar";particle [11] = "e^-";particle [-11] = "e^+";particle [12] = "\!\(\* SubscriptBox [\"ν\", \"e\"]\) ";particle [-12] ="\!\(\* SuperscriptBox[SubscriptBox [\"\[ Vee]\", \"e\"], \"c\"]\)";particle [13] = "\!\(\* SuperscriptBox [\"μ\", \" -\"]\)";particle [-13] = "\!\(\* SuperscriptBox [\"μ\",  \"+\"]\) ";particle [14] = "\!\(\* SubscriptBox [\"ν\", \"μ\"]\) ";particle [-14] ="\!\(\* SuperscriptBox[SubscriptBox [\"ν\", \"μ\"], \"c\"]\)";particle [16] = "\!\(\* SubscriptBox [\"ν\", \"τ\"]\) ";particle [-16] ="\!\(\* SuperscriptBox[SubscriptBox [\"ν\", \"τ\"], \"c\"]\)";particle [21] = "gluon";
particle [22] = "a";particle [24] = "W^+";particle [1000006] = Subscript[OverTilde[t], 1];particle [ -1000006] =\!\(\* SuperscriptBox[SubscriptBox[RowBox [{" OverTilde", "[", "t", "]"}], "1"],  "*"]\);particle [9000005] = "dm";finout [-1] = "in";finout [1] = "out";finout [2]   = "decayed";

Totxsec[rawInput_]:=(rawInput[[Position[rawInput,{"<init>"}][[1]]+2,1]])[[1]];

ReadME[rawInput_]:=(
(*Search for beginning of events*)pos=Position[rawInput,{"</init>"}][[1,1]];(*Throw away crap at start of file,combine events in {} array and remove the XML commands and\compulsory eventinfo*)DeleteCases[Split[Drop[rawInput,pos],#=!= {"</event>"}&],{"<event>"}| {"</LesHouchesEvents>"}|{"</event>"}|eventInfo,2]);

EventPrint[event_]:=(Print[Join[{{"pid","in/out","mother1","mother2","color1","color2","px","py","pz","p0","mass","cτ","hel"}},event/.{x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_,x9_,x10_,x11_,x12_,x13_}->{particle[x1],finout[x2],x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13}]//MatrixForm]);

(*those are input or decayed particles*)decayedOrIn={_,-1,_,_,_,_,_,_,_,_,_,_,_}|{_,2_,_,_,_,_,_,_,_,_,_,_};EventOut[event_]:=(DeleteCases[event,decayedOrIn]);SetOut[event_]:=(DeleteCases[event,decayedOrIn,2]);

(*EffMassAll[objList_]:=Plus@@Map[ptOf,objList//Transpose,{2}];*)EffMass[event_]:=Plus@@Map[ptOf,event,{1}];pT[event_]:=Map[ptOf,event,{1}]//Flatten;eta[event_]:=Map[etaOf,event,{1}]//Flatten;theta[event_]:=Map[thetaOf,event,{1}]//Flatten;threeVector[event_]:=Map[ThreeVectorFrom,event,{1}];MissingEnergy[event_]:=Plus@@Map[enOf,event,{1}];ptOf[{_,_,_,_,_,_,px_,py_,___}]:=Sqrt[px^2+py^2];enOf[{_,_,_,_,_,_,_,_,_,En_,___}]:=En;thetaOf[{_,_,_,_,_,_,px_,py_,pz_,___}]:=ArcCos[pz/Sqrt[px^2+py^2+pz^2]];etaOf[{_,_,_,_,_,_,px_,py_,pz_,___}]:=-Log[Abs[Tan[ArcCos[pz/Sqrt[px^2+py^2+pz^2]]/2]]];FourVectorFrom[{_,_,_,_,_,_,px_,py_,pz_,En_,___}]:={En,px,py,pz};FourLength[{pe_,pz_,px_,py_}]:=Sqrt[Max[pe^2-pz^2-px^2-py^2,0.0]];

ThreeVectorFrom[{_,_,_,_,_,_,px_,py_,pz_,___}]:={px,py,pz};

CosThetabetween[event_,i_,j_]:=((ThreeVectorFrom[event[[i]]].ThreeVectorFrom[event[[j]]])/(Norm[ThreeVectorFrom[event[[i]]]] Norm[ThreeVectorFrom[event[[j]]]]));

CosthetaTwoJet[{jet1_,jet2_}]:=(a=ThreeVectorFrom[jet1];b=ThreeVectorFrom[jet2];a.b/Sqrt[a.a b.b]);GetAll[evt_]:=evt;all[evt_]:=True;Freq[evtList_,crit_,objsel_,plotfunc_,{min_,max_,nbins_}]:=BinCounts[Flatten[plotfunc/@Flatten4[objsel/@Select[evtList,crit]]],{min,max,(max-min)/nbins}];Hist[freqcombo__,{min_,max_,nbins_},HistogramOptions___]:=Histogram[freqcombo/.Bin[args__]->Freq[args,{min,max,nbins}],FrequencyData->True,HistogramCategories->Table[i,{i,min,max,(max-min)/nbins}],HistogramOptions];

Flatten4[List_]:=Flatten[List,Max[0,Depth[List]-4]];HGet[patt_,All][evt_]:=Cases[evt,patt];oJetlhe={2,1,_,_,_,_,_,_,_,_,_,_,_}|{-2,1,_,_,_,_,_,_,_,_,_,_,_}|{21,1,_,_,_,_,_,_,_,_,_,_,_};oMissingEnergy={8888,1,_,_,_,_,_,_,_,_,_,_,_};PTSort[event_]:=Sort[event,Greater[ptOf[#1],ptOf[#2]]&];and[f_][evt_]:=f[evt];and[f_,g__][evt_]:=f[evt]&&and[g][evt];or[f_][evt_]:=f[evt];or[f_,g__][evt_]:=f[evt]||or[g][evt];join[f_][evt_]:=f[evt];join[f_,g__][evt_]:=Join[f[evt],join[g[evt]]];

Belle II

ALP decay

Width and decay length

The Lagrangian is ℒ = ((Subscript[g, aγγ])/4) a F*(OverTilde[F])

(* a -> γγ decay width *)
Γ aγγ[gaγγ_,ma_]:=((((gaγγ)^2) (ma^3))/(128 π));

(* Boost factor γ*β = |p|/m *)
γβ a[Ea_,ma_]:=((Sqrt[(Ea^2) - (ma^2)])/ma) HeavisideTheta[Ea-ma];

(* Boosted decay length, converted to meters *)
γβcτ a[gaγγ_,Ea_,ma_]:= γβ a[Ea,ma] hbarc ((Γ aγγ[gaγγ,ma])^(-1));

Prompt vs Displaced plot

(* ALP produced via 2 -> 2 scattering will have energy Ea = ((s + (Subscript[m, a]^2))/(2 (Sqrt[s]))) *)Ea22[ma_]:=((sB2 + (ma^2))/(2 (Sqrt[sB2])));

Roughly, the ALP is prompt if γβcτ < (Subscript[L, min]), displaced if (Subscript[L, min])<γβcτ < (Subscript[L, max]) and invisible if γβcτ > (Subscript[L, max])

RegionPlot[{γβcτ a[gaγγ,Ea22[ma],ma]<Lmin, (* Prompt *)
Lmin<γβcτ a[gaγγ,Ea22[ma],ma]<Lmax, (* Displaced *)
γβcτ a[gaγγ,Ea22[ma],ma]>Lmax}, (* Invisible *)
{ma,(10^(-3)),10},{gaγγ,(10^(-10)),1},ScalingFunctions->{"Log","Log"},PlotLegends->{"Prompt","Displaced","Invisible"},Evaluate[MyOptions],FrameLabel->{"(Subscript[m, a]) (GeV)","(Subscript[g, aγγ]) ((GeV^(-1)))"}]

Prompt vs Displaced plot (complete)

(* Boosted separation angle vs ALP mass*)

Δθ[ma_]:=θ/.FindRoot[θ ==ArcTan[((ma/((Sqrt[sB2])/2)))((Sin@(π/2))/(Cos@(π/2)+1))],{θ,.1}];

γResolved limit = m/.FindRoot[Interpolation[Table[{m,Δθ[m]},{m,0.01,1,.01}] ][m] == 2 Δθ γ,{m,.1}]

RegionPlot[{
Piecewise[{{γβcτ a[gaγγ,Ea22[ma],ma],ma>=γResolved limit},{Null,ma<γResolved limit}}]<Lmin, (* Prompt (Resolved) *)
Piecewise[{{γβcτ a[gaγγ,Ea22[ma],ma],ma<γResolved limit},{Null,ma>=γResolved limit}}]<Lmin,(* Prompt (Merged) *)
Lmin<γβcτ a[gaγγ,Ea22[ma],ma]<Lmax, (* Displaced *)
γβcτ a[gaγγ,Ea22[ma],ma]>Lmax}, (* Invisible *)
{ma,(10^(-3)),10},{gaγγ,(10^(-10)),1},ScalingFunctions->{"Log","Log"},PlotLegends->{"Prompt (Resolved)","Prompt (Merged)","Displaced","Invisible"},Evaluate[MyOptions],FrameLabel->{"(Subscript[m, a]) (GeV)","(Subscript[g, aγγ]) ((GeV^(-1)))"}]

Associated ALP production

Photon mediated 2 -> 2 scattering: (SuperPlus[e])((Subscript[p, 1]))(SuperMinus[e])((Subscript[p, 2])) -> γ(k) a((Subscript[p, a])), assuming massless electrons

Amplitude and cross section calculation with FeynCalc

Define spinor term and amplitude

ClearAll[electron,ALP γ, iM];

electron[μ_] := SpinorUBar[p1].GA[μ].SpinorU[p2];

ALP γ[μ_,ν_] :=4 LeviCivita[μ,ν,ρ,τ] FV[q,ρ]FV[k,τ] PolarizationVector[k,ν] (*4 (FV[q,ν]FV[k,μ]-SP[q,k] MT[μ,ν])PolarizationVector[k,ν];*);

iM[μ_,ν_] := electron[μ] ALP γ[μ,ν]

Get amplitude squared

(* To compute the amplitude squared, we need to sum over final spin states (FermionSpinSum function) and the vector polarization states (DoPolarizationSums function). For the latter we need to specify the mass of the boson, thus we put it on shell by defining k*k = 0  *)

FCClearScalarProducts[];
SP[k,k]=0;
SP[p2,p2]=0;
SP[p1,p1]=0;

((((1/(SP[q,q]))))^2)DoPolarizationSums[FermionSpinSum[iM[μ,ν] ComplexConjugate[iM[μ1,ν1]],ExtraFactor->1/4],k,0]/.DiracTrace->Tr

Add COM kinematics and couplings to |(M^2)|

FCClearScalarProducts[];
SP[k,k]=0;
SP[p2,p2]=0;
SP[p1,p1]=0;

SP[p1,p2]=s/2;
SP[k,p1]=((Sqrt[s])/2)(((s-(ma^2))/(2 (Sqrt[s])))) (1 - Cos[θ]);
SP[k,p2]=((Sqrt[s])/2)(((s-(ma^2))/(2 (Sqrt[s])))) (1 + Cos[θ]);

M2AssALP = ((((((Sqrt[α 4 π])gaγγ)/4)))^2)FullSimplify[((((1/(SP[q,q]))))^2)DoPolarizationSums[FermionSpinSum[iM[μ,ν] ComplexConjugate[iM[μ1,ν1]],ExtraFactor->1/4],k,0]/.DiracTrace->Tr/.q->p1+p2//ScalarProductExpand//DiracSimplify,Trig->False]

(* ALP momentum |(Subscript[p, a])| *)
FullSimplify[(Sqrt[(((((s + (ma^2))/(2 (Sqrt[s])))))^2) - (ma^2)]),Assumptions->{ma>0,s>ma^2}]

Get differential cross section

((dσ)/(dcos(θ))) = (1/(8 π s (Sqrt[s]))) |(Subscript[p, a])| |M((|)^2)

FullSimplify[((1/(8π s (Sqrt[s]))))(((s-(ma^2))/(2 (Sqrt[s]))))M2AssALP/.Ea->((s + (ma^2))/(2 (Sqrt[s]))),Assumptions->{ma>0,s>(ma^2)},Trig->False]

dσProd dcosθ[gaγγ_,ma_,x_]:= (( αem ((gaγγ)^2))/32)((x^2) + 1)(((1 - ((ma^2)/sB2)))^3);

Differential cross section

ClearAll[dσProd dcosθ];

dσProd dcosθ[gaγγ_,ma_,x_]:= (( αem ((gaγγ)^2))/32)((x^2) + 1)(((1 - ((ma^2)/sB2)))^3); (* Changed variable to x = Cos[θ] *)

Plot[dσProd dcosθ[(10^(-3)),1,x],{x,-1,1},Evaluate[MyOptions],FrameLabel->{"x","(Subscript[σ, prod]) ((GeV^(-2)))"},PlotLabel->"(Subscript[m, a]) = 1 GeV, (Subscript[g, aγγ]) = (10^(-3)) (GeV^(-2))"]

Integrate[(( α ((gaγγ)^2))/32)((x^2) + 1)(((1 - ((ma^2)/s)))^3),x]

(* Integrate the differential cross section over the angular coverage of BelleII detector *)

Integrate[(( αem ((gaγγ)^2))/32)((x^2) + 1)(((1 - ((ma^2)/sB2)))^3),{x,Cos[θmax],Cos[θmin]}]//StandardForm

(* Total cross section *)
σProd[gaγγ_,ma_]:=0.00035623881567210114` ((gaγγ)^2) (((1-0.0088999644001424` (ma^2)))^3);

ContourPlot[σProd[g,m],{m,(10^(-4)),10},{g,(10^(-7)),1},ContourShading->None,Contours->Table[(10^i),{i,-20,-4,2}],ScalingFunctions->{"Log","Log"},Evaluate[MyOptions],FrameLabel->{"(Subscript[m, a]) (GeV)","(Subscript[g, aγγ]) ((GeV^(-1)))"}]

LogLogPlot[{σProd[g,.01],σProd[g,.1],σProd[g,1],σProd[g,6],σProd[g,9.5]},{g,(10^(-7)),1},Evaluate[MyOptions],FrameLabel->{"(Subscript[g, aγγ]) ((GeV^(-1)))","(Subscript[σ, prod]) ((GeV^(-2)))"},PlotLegends->{"(Subscript[m, a]) = 0.01 GeV","(Subscript[m, a]) = 0.1 GeV","(Subscript[m, a]) = 1 GeV","(Subscript[m, a]) = 6 GeV","(Subscript[m, a]) = 9.5 GeV"}]

Invisible region: bounds (assuming no backgrounds)

(* Here we convert σ to pb units *)

N invisible[Lum_,gaγγ_,ma_]:= Lum (σProd[gaγγ,ma]/pb) Exp[-(Lmax/(γβcτ a[gaγγ,Ea22[ma],ma]))];

Bounds ALP invisible =RegionPlot[{N invisible[LumB2,gaγγ,ma]>3,N invisible[LumB2 final,gaγγ,ma]>3},{ma,(10^(-3)),10},{gaγγ,(10^(-8)),1},ScalingFunctions->{"Log","Log"},Evaluate[MyOptions],FrameLabel->{"(Subscript[m, a]) (GeV)","(Subscript[g, aγγ]) ((GeV^(-1)))"},PlotLegends->{"ℒ = 445 (pb^(-1))","ℒ = 50 (ab^(-1))"}]//Quiet

Visible region

(SuperPlus[e])(SuperMinus[e])->γγγ background

This is a 2 -> 3 scattering process, not impossible to compute “by hand” but very long and messy, so let’s be smart (or lazy) and simulate it with MadGraph.

Boost angles to COM

η[θ_]:=-Log[Tan[θ/2]];

θminCM=θminCM/.FindRoot[θmin ==ArcTan[((1/(.28)))((Sin@θminCM)/(Cos@θminCM+1))],{θminCM,θmin}]//Quiet
%/Degree

θmaxCM=θmaxCM/.FindRoot[(π-θmax) ==ArcTan[((1/(.28)))((Sin@θmaxCM)/(-Cos@θmaxCM+1))],{θmaxCM,θmax}]//Quiet
%/Degree

η[θminCM]
η[θmaxCM]

Upload MG simulation

ee3γ MG = Import["ee_aaa.lhe","Table"];

xsec =Totxsec[ee3γ MG];

eventList=ReadME[ee3γ MG];

eventListOut=SetOut[eventList];

Nev = Length@eventList - 1;

Create lists of photon momenta and angles

(* The function FourVectorForm gives the 4-momentum of an event, which for convenience is split in first component (energy) and spatial components. The function thetaOf gives the angle with respect to the z-axis *)

γ1List = Table[{FourVectorFrom[eventListOut[[i,1]]][[1]],FourVectorFrom[eventListOut[[i,1]]][[2;;4]],thetaOf[eventListOut[[i,1]]]},{i,Nev}];

γ2List = Table[{FourVectorFrom[eventListOut[[i,2]]][[1]],FourVectorFrom[eventListOut[[i,2]]][[2;;4]],thetaOf[eventListOut[[i,2]]]},{i,Nev}];

γ3List = Table[{FourVectorFrom[eventListOut[[i,3]]][[1]],FourVectorFrom[eventListOut[[i,3]]][[2;;4]],thetaOf[eventListOut[[i,3]]]},{i,Nev}];

Belle II 3γ acceptance (from Fig. 6 of 1709.00009)

acceptance HighMass = Import["B2_acceptance_HighMass.csv","Data"];
acceptance HighMass = {((acceptance HighMass[[All,1]])^2),acceptance HighMass[[All,2]]}//Transpose;
acceptance HighMass Interp = Interpolation[acceptance HighMass,InterpolationOrder->1][x];

acceptance LowMass = Import["B2_acceptance_LowMass.csv","Data"];
acceptance LowMass = {((acceptance LowMass[[All,1]])^2),acceptance LowMass[[All,2]]}//Transpose;
acceptance LowMass Interp = Interpolation[acceptance LowMass,InterpolationOrder->1][x];

acceptance = Piecewise[{{acceptance HighMass Interp,4.5<=x<100},{acceptance LowMass Interp,3 ((γResolved limit)^2)<=x<4.5}}];

(* "fake" acceptance, it is similar to the one from the Belle II projection, but it grows linearly for low masses and then it cut sharply at high masses *)

l = 90;l2=100;
maxϵ = 0.55;

acceptance fake[x_] := Piecewise[{{0,x<((γResolved limit)^2)},{((maxϵ)/((Log[ (3.4/((γResolved limit)^2)) ])^2))((Log[ (x/((γResolved limit)^2)) ])^2),((γResolved limit)^2)<=x<3.4},{maxϵ, 3.4<=x<=l},{maxϵ + .06(l-x),l<x<=l2}}]

Show[LogLinearPlot[{acceptance,acceptance fake[x]},{x,((γResolved limit)^2),110},Exclusions->None,PlotRange->All,FrameLabel->{"(Subscript[m, a]^2) ((GeV^2))","ϵ"},Evaluate[MyOptions]],ListLogLinearPlot[{aa,bb},PlotStyle->LightBlue]]

Invariant mass

(* (Subscript[M, γ,recoil]^2) list *)

m1 recoil  = Table[sB2 - 2 (Sqrt[sB2])γ1List[[i,1]],{i,Length@γ1List}];

(* Define PDF function of recoil invariant mass *)

D1 recoil = PDF[HistogramDistribution[m1 recoil,300],x];

(* Define background as 3γ cross section x luminosity x acceptance *)
ClearAll[background,background fake];

background[Lum_] := Lum xsec D1 recoil acc(*acceptance Interp*);
background fake[Lum_]:= Lum xsec D1 recoil acceptance fake[x];

background Lower[Lum_] := (((.8))^2)Lum xsec D1 recoil acc(*acceptance Interp*);
background Lower fake[Lum_]:=(((.8))^2) Lum xsec D1 recoil acceptance fake[x];

Upload Belle II data

B2 data Upper =Select[ Import["B2_data.csv","Data"],#[[1]]>16 &];

B2 data stat Upper =Table[{B2 data Upper[[i,1]],Around[B2 data Upper[[i,2]],(Sqrt[B2 data Upper[[i,2]]])]},{i,Length@B2 data Upper}];

B2 data Lower = Delete[Import["B2_data_Lower.csv","Data"],1];
B2 data stat Lower =Table[{B2 data Lower[[i,1]],Around[B2 data Lower[[i,2]],(Sqrt[B2 data Lower[[i,2]]])]},{i,Length@B2 data Lower}];

B2 data =B2 data Lower~Join~ B2 data Upper;
B2 data stat = B2 data stat Upper[[All,2,1]]~Join~B2 data stat Lower[[All,2,1]];

Show[Plot[{background[LumB2],background fake[LumB2]},{x,16,110},Exclusions->None,Evaluate[MyOptions],PlotRange->All,FrameLabel->{"(Subscript[M, γ]^2) ((GeV^2))","Events ((GeV^(-2)))"}],ListPlot[B2 data stat Upper,PlotStyle->Black]]

Show[Plot[{background Lower[LumB2], background Lower fake[LumB2]},{x,0,20},Exclusions->None,Evaluate[MyOptions],PlotRange->All,FrameLabel->{"(Subscript[M, γ]^2) ((GeV^2))","Events ((GeV^(-2)))"}],ListPlot[B2 data stat Lower,PlotStyle->Black]]

Plot from 2007.13071

Add signal to background and get bounds

Add ALP signal to SM background

ClearAll[signal,signal bg];

signal[Lum_,g_,m_] := Lum (σProd[g,m]/pb) PDF[NormalDistribution[(m^2),.03 (m^2)],x] acc fake[x];

signal bg[Lum_,g_,m_] :=background fake[Lum] + signal[Lum,g,m];

signal bg Lower[Lum_,g_,m_] :=background Lower fake[Lum] + signal[Lum,g,m];

Show[
Plot[{
background fake[LumB2], (* BG only *)
signal bg[LumB2,3 (10^(-3)),5],(* BG + signal at (Subscript[m, a]) = 5 GeV *)
signal bg[LumB2,2 (10^(-2)),9](* BG + signal at (Subscript[m, a]) = 9 GeV *)},{x,16.1,110},Exclusions->None,Evaluate[MyOptions],PlotRange->All,FrameLabel->{"(Subscript[M, γ]^2) ((GeV^2))","Events ((GeV^(-2)))"}],ListPlot[B2 data stat Upper,PlotStyle->Black]]//Quiet

Show[Plot[{background Lower fake[LumB2],signal bg Lower[LumB2,5 (10^(-4)),1],signal bg Lower[LumB2, (10^(-3)),3]},{x,0,16},Exclusions->None,Evaluate[MyOptions],PlotRange->All,FrameLabel->{"(Subscript[M, γ]^2) ((GeV^2))","Events ((GeV^(-2)))"}],ListPlot[B2 data stat Lower,PlotStyle->Black]]//Quiet

(χ^2)  function

ClearAll[χ2 Upper,χ2 Lower,χ2 bgOnly];

(* Define (χ^2) for the Upper part of the data *)χ2 Upper[Lum_,g_,m_] :=Total@DeleteCases[Table[((((B2 data stat Upper[[i,2,1]] -(signal bg[LumB2,g,m]/.x->B2 data stat Upper[[i,1]]) )/((Sqrt[(LumB2/Lum)])(Sqrt[B2 data stat Upper[[i,2,1]]]))))^2),{i,Length@B2 data stat Upper}],{_,_}]//Quiet;

(* Define (χ^2) for the Lower part of the data *)χ2 Lower[Lum_,g_,m_] :=Total@DeleteCases[Table[((((B2 data stat Lower[[i,2,1]] -(signal bg Lower[LumB2,g,m]/.x->B2 data stat Lower[[i,1]]) )/((Sqrt[(LumB2/Lum)])(Sqrt[B2 data stat Lower[[i,2,1]]]))))^2),{i,Length@B2 data stat Lower}],{_,_}]//Quiet;

(* Define (χ^2) assuming Data == Background *)χ2 bgOnly[Lum_,g_,m_] :=Total@DeleteCases[Table[((((signal[LumB2,g,m]/.x->B2 data[[i,1]])/((Sqrt[(LumB2/Lum)])(Sqrt[B2 data stat[[i]]]))))^2),{i,Length@B2 data stat}],{_,_}]//Quiet;

χ2 Upper[LumB2,0,0]
χ2 Lower[LumB2,0,0]
χ2 bgOnly[LumB2,0,0]

Get limits

(* Create list of mass vs limit for the Upper part of the data *)Limit Upper LowLum =  Table[{m,g/.FindRoot[Abs[χ2 Upper[LumB2,g,m]-χ2 Upper[LumB2,0,0]]==3.8,{g,(10^(-3))}]},{m,4,10,.3}];//Quiet

(* Create list of mass vs limit for the Lower part of the data *)Limit Lower LowLum =  Table[{(Sqrt[B2 data Lower[[i,1]]]),g/.FindRoot[Abs[χ2 Lower[LumB2,g,(Sqrt[B2 data Lower[[i,1]]])]-χ2 Lower[LumB2,0,0]]==3.8,{g,(10^(-2))}]},{i,Length@B2 data Lower}];//Quiet

(* Join the two lists*)
Limit LowLum = Limit Lower LowLum~Join~Limit Upper LowLum;

(* Create list of mass vs limit assuming Data == Background *)Limit All LowLum =  Table[{(Sqrt[B2 data[[i,1]]]),g/.FindRoot[χ2 bgOnly[LumB2,g,(Sqrt[B2 data[[i,1]]])]-χ2 bgOnly[LumB2,0,0]==3.8,{g,(10^(-3))}]},{i,Length@B2 data}];//Quiet

(* Do the same as before for the 50 (ab^(-1)) luminosity case *)
Limit Upper HiLum =  Table[{m,g/.FindRoot[Abs[χ2 Upper[LumB2 final,g,m]-χ2 Upper[LumB2 final,0,0]]==3.8,{g,(10^(-3))}]},{m,4,10,.3}];//Quiet

Limit Lower HiLum =  Table[{(Sqrt[B2 data Lower[[i,1]]]),g/.FindRoot[Abs[χ2 Lower[LumB2 final,g,(Sqrt[B2 data Lower[[i,1]]])]-χ2 Lower[LumB2 final,0,0]]==3.8,{g,(10^(-3))}]},{i,Length@B2 data Lower}];//Quiet

(* Join the two lists*)
Limit HiLum = Limit Lower HiLum~Join~Limit Upper HiLum;

Limit All HiLum =  Table[{(Sqrt[B2 data[[i,1]]]),g/.FindRoot[χ2 bgOnly[LumB2 final,g,(Sqrt[B2 data[[i,1]]])]-χ2 bgOnly[LumB2 final,0,0]==3.8,{g,(10^(-3))}]},{i,Length@B2 data}];//Quiet

ListLogLogPlot[{Limit LowLum,Limit All LowLum,Limit HiLum,Limit All HiLum},PlotRange->{{(10^(-3)),10},{(10^(-7)),(10^(-1))}},Evaluate[MyOptions],Joined->True,Filling->Top,FrameLabel->{"(Subscript[m, a]) (GeV)","(Subscript[g, aγγ]) ((GeV^(-1)))"},PlotLegends->{"LowLum bound","LowLum bound (Data = BG)","HiLum bound","HiLum bound (Data = BG)"}]

Bounds ALP visible =ListLogLogPlot[{Limit LowLum,Limit HiLum},PlotRange->{{(10^(-3)),10},{(10^(-7)),(10^(-1))}},Evaluate[MyOptions],Joined->True,Filling->Top,FrameLabel->{"(Subscript[m, a]) (GeV)","(Subscript[g, aγγ]) ((GeV^(-1)))"}]

Combine plots

Show[Bounds ALP visible,Bounds ALP invisible]

Electron Beam Dump

Primakoff production cross section

Γ aγγ[gaγγ_,ma_]:=((((gaγγ)^2) (ma^3))/(128 π));

me = 511 keV;

(* Photon distribution, x = Eγ/(Subscript[E, beam]) ~ (Subscript[E, a])/(Subscript[E, beam]) *)
γDist[x_]:=((αem)/(2π)) (((1 +(((1-x))^2))/(4 (x^3) me)));

ap[Z_]:=112 (Z^(-1/3))/me;

tmin[ma_,Eγ_]:=((ma^4)/(4 ((Eγ)^2)));

d[A_]:=0.164 (A^(-2/3)) (* in (GeV^2)*); 

σPrim[A_,Z_,ma_,gaγγ_,Eγ_]:=((16 π αem)/(ma^3)) Γ aγγ[gaγγ,ma] ((Z^2)/2) (Log[((d[A])/(((ap[Z])^(-2)) + tmin[ma,Eγ]))] - 2);

E137 experiment

Ne OnTarget = 2 (10^20); E beam = 20 (* GeV*);L shield = 179;L decay = 204;Δx = 1.5;Δy = 1; (* All lengths are in meters *)

(* Target is made of Aluminium *)
Z Al = 13; A Al = 27;
X Al = 24.3 (* Radiation length, units of g/(cm^2)*);
ρ Al =2.699 (* Density of Al in g/(cm^3) *) ;
ρ Fe =7.874 (* Density of Al in g/(cm^3) *) ;
X Fe = 13.84 (* Radiation length, units of g/(cm^2)*);
m Al = 26.98 (* in g/mol*);
NA = 6.022 (10^23) (* Avogadro number, (mol^(-1)) *);

(* Detector length along z is "8 radiation lengths" of iron or aluminium *)8 X Al/ρ Al (* in cm *)
8 X Fe/ρ Fe (* in cm *)