PPPC Machine: functions only (run this)

Flip Tanedo
24 Nov 2014: I retired the original “PPPCMachine.nb” because it was too long. I’m saving it as PPPCMachine_Historical because it contains lots of comments, but it’s basically un-readable. This file used to be PPPCMachineFuncitonsOnly and is now the main file. 
Instructions: just run the whole notebook, it shouldn’t take long. 

6 Nov 2014: Added some new helper functions for doing runs. Fixed error in p3ϕ2gam functions. Fixed error in p3ϕ2gamfaster: WeightedSpectrum[Egam_] was missing a factor of ΔEϕ. 
4 Nov 2014: this is a streamlined version of PPPC Machine; the idea is that you can just run the notebook to load all relevant functions. Based on the 3 November version. I’ve removed binning functions for now. Added DC’s SaveIt and ReadIt functions.
3 November 2014: Minor bug fixes.
2 Nov 2014: 2-3 spectra, faster 2-2 interpolations 
30 October 2014: kinematics and functions for on-shell mediators
21 October 2014

Machine for fitting FERMI diffuse γ-ray spectra to simple models of DM annihilation to (i) 2 quarks, (ii) 2 (axial-)vector mediators that decay to quarks, (iii) 3 pseudoscalar mediators that decay to quarks. We assume that the vectors have universal coupling and the pseudoscalar has mass-weighted coupling so that it is well approximated by a b-philic coupling.

Based on 2to2scanrevised.nb, HeavyHooperonPPPC.nb, 2to3Amplitude.nb

Goal: 
1. output PPPC binned spectra
2. take binned data and compare to PPPC spectra
`	a) compare to a histogram
	b) compare to a histogram with error bars
	c) compare to a histogram with degeneracies
	
Required: PPPC files. http://www.marcocirelli.net/PPPC4DMID.html. Download these somewhere and define the path below. See www.marcocirelli.net/PPPC4DMID.html for more details

ClearAll["Global`*"] (* This is just good form *)

User-specific input

Input your PPPC location and run

Go to Insert > File Path to insert the path to dldNdlxEW.m

Get["/Users/fliptomato/Documents/Work/PostLUX/pppc/dlNdlxEW.m"];

Define functions

Photon spectrum

Instructions from PPPC sample file

Fluxes are given as (Subscript[Log, 10])[((d N)/(d (Subscript[Log, 10])x))]  (x = K/(Subscript[m, DM]), with K the kinetic energy in GeV)  normalized per single DM annihilation. 
The format is

	dlNdlxIEW[primary -> secondary][mass, lx]

where

	primary = eL, eR, e, μL, μR, μ, τL, τR, τ, q, c, b, t, WL, WT, W, ZL, ZT, Z, g, γ, h, νe, νμ, ντ, V->e, V->μ, V->τ
	secondary = e, p, γ, d, νe, νμ, ντ	
	mass = (Subscript[m, DM])in GeV, on the range (Subscript[m, DM])= 5 GeV -> 100 TeV
	lx = Log[10, x ], on the range x = (10^(-9))-> 1 except for the channels V->e, V->μ, V->τ for which x = (10^(-5))-> 1

Warning:	Extrapolations are performed close to the kinematical threshold for massive primaries (e.g. in the case DM DM -> ZZ with (Subscript[m, DM])=93 GeV): the output fluxes should be reliable, but instabilities are possible. Also, instabilities very close to the edges of the domain in x are possible.

Generic primary to photon spectrum function

pE is the primary energy, i.e. the mass of the DM for pairwise annihilation to SM pairs.

p2gam[primary_,pE_,photonE_]:=If[photonE<pE,
((10^(dlNdlxIEW[primary->"γ"][pE,Log[10,(photonE/pE)]]))/(Log[10.] photonE )),
0 (* if photon energy is larger than injection energy *)
] (* usage: primary is of the form "τ", WITH the quotation marks *)

Discussion: Boosted Spectra

We give a brief discussion of boosted spectra.

The photon spectrum is given by

(((Subscript[dN, γ])^tot)/(Subscript[dE, γ]))=∫(Subscript[dE, SM]) (((Subscript[dN, SM])^tot)/(Subscript[dE, SM])) (((Subscript[dN, γ])((Subscript[E, SM])))/(Subscript[dE, γ]))

Where the spectrum of Standard Model primaries is a delta function for the case of annihilation of two DM particles into two SM particles. This is fixed purely from kinematics.

In the case of annihilation to on-shell mediators, however, we know that (((Subscript[dN, SM])^tot)/(Subscript[dE, SM])) has its own spectrum that is a box in the case of two on-shell mediators.

(((Subscript[dN, SM])^tot)/(Subscript[dE, SM]))=∫(Subscript[dE, ϕ]) ((Subscript[dN, ϕ])/(Subscript[dE, ϕ])) (Subscript[((((Subscript[dN, SM])((Subscript[E, ϕ])))/(Subscript[dE, SM]))), lab])

Note that (Subscript[((((Subscript[dN, SM])((Subscript[E, ϕ])))/(Subscript[dE, SM]))), lab]) is a box spectrum for a fixed (Subscript[E, ϕ]). This is because (Subscript[E, ϕ]) in the lab frame is

(Subscript[E, ϕ])=γ ((Subscript[m, ϕ])/2)+γβ(Sqrt[(((Subscript[m, ϕ])^2)/4)-((Subscript[m, SM])^2)])Cos[θ]

Where Cos[θ] is the angle with respect to the boost direction. We have to average over this angle to account for the isotropy of the DM annihilation. The boost parameter γ is the boost from the ϕ rest frame to the lab frame. In particular,

γ = ((Subscript[E, ϕ])/(Subscript[m, ϕ]))
γβ=(Sqrt[(((((Subscript[E, ϕ])/(Subscript[m, ϕ]))))^2)-1])

Then, in terms of the unit box function (box of unit length and height, centered around zero,

(Subscript[((((Subscript[dN, SM])((Subscript[E, ϕ])))/(Subscript[dE, SM]))), lab])= (Norm) UnitBox[(( ((Subscript[E, SM])-((Subscript[E, ϕ])/2)))/(2(Sqrt[(((((Subscript[E, ϕ])/(Subscript[m, ϕ]))))^2)-1])(Sqrt[(((Subscript[m, ϕ])^2)/4)-((Subscript[m, SM])^2)])))]
(Norm) =((Number of SM particles per ϕ)/(Width of the box))=(2/(2(Sqrt[(((((Subscript[E, ϕ])/(Subscript[m, ϕ]))))^2)-1])(Sqrt[(((Subscript[m, ϕ])^2)/4)-((Subscript[m, SM])^2)])))

So ultimately we have

(((Subscript[dN, γ])^tot)/(Subscript[dE, γ]))=∫(Subscript[dE, SM]) (∫(Subscript[dE, ϕ]) ((Subscript[dN, ϕ])/(Subscript[dE, ϕ])) (Norm) UnitBox[(( ((Subscript[E, SM])-((Subscript[E, ϕ])/2)))/(2(Sqrt[(((((Subscript[E, ϕ])/(Subscript[m, ϕ]))))^2)-1])(Sqrt[(((Subscript[m, ϕ])^2)/4)-((Subscript[m, SM])^2)])))]) (((Subscript[dN, γ])((Subscript[E, SM])))/(Subscript[dE, γ]))

Boosted Spectra: 2 to 2

Now we define a set of tools for annihilation modes to two on-shell [vector] mediators which subsequently decay into pairs of SM fermions. 

Notes: I’m not building in any automatic “zero if the kinematics don’t allow on-shell mediators.” Users should be careful.

Functions

Following my HeavyHooperonPPPC code, I’ll define a full funciton (p2pϕ2gam) which does a numerical integral over the SM final state energies (whose range comes from the boosted mediator). For scans it’s more useful to define an interpolating function for the spectrum

p2ϕ2gam[primary_,mprimary_,mχ_,mϕ_,photonE_]:=Module[
{finalstate=primary,mSM=mprimary,mchi=mχ, mphi=mϕ,eγ=photonE,γ,γβ, E0, norm,EfsHi,EfsLo},
γ=(mchi/mphi);
γβ=(Sqrt[((((mchi/mphi)))^2)-1]);
E0=mphi/2;
norm=(1/(γβ(Sqrt[((((mphi/2)))^2)-(mSM^2)])));
(*Just a sanity check, the (Subscript[dN, ϕ])/(Subscript[dE, ϕ]) distribution is*)
(*norm UnitBox[((eγ(*temp*) - E0)/(2 γβ(Sqrt[((((mphi/2)))^2)-(mSM^2)])))];*)
EfsHi=γ E0+γβ (Sqrt[((((mphi/2)))^2)-(mSM^2)]);
EfsLo=γ E0-γβ (Sqrt[((((mphi/2)))^2)-(mSM^2)]);
norm NIntegrate[p2gam[finalstate,Efs,eγ],{Efs,EfsLo,EfsHi}]
]

(* alternate version with "SymbolicProcessing->0", factor of 8 faster *)
p2ϕ2gamquick[primary_,mprimary_,mχ_,mϕ_,photonE_]:=Module[
{finalstate=primary,mSM=mprimary,mchi=mχ, mphi=mϕ,eγ=photonE,γ,γβ, E0, norm,EfsHi,EfsLo},
γ=(mchi/mphi);
γβ=(Sqrt[((((mchi/mphi)))^2)-1]);
E0=mphi/2;
norm=(1/(γβ(Sqrt[((((mphi/2)))^2)-(mSM^2)])));
(*Just a sanity check, the (Subscript[dN, ϕ])/(Subscript[dE, ϕ]) distribution is*)
(*norm UnitBox[((eγ(*temp*) - E0)/(2 γβ(Sqrt[((((mphi/2)))^2)-(mSM^2)])))];*)
EfsHi=γ E0+γβ (Sqrt[((((mphi/2)))^2)-(mSM^2)]);
EfsLo=γ E0-γβ (Sqrt[((((mphi/2)))^2)-(mSM^2)]);
norm NIntegrate[p2gam[finalstate,Efs,eγ],{Efs,EfsLo,EfsHi},Method->{Automatic,"SymbolicProcessing"->0}]//Quiet (* seems to give slow convergence warnings, only off at 4th decimal? *)
]

Fast Functions

p2ϕ2gam takes a long time to plot and, on top of that, can be tempramental when making log plots. So we’ll follow our approach in HeavyHooperonPPPC.nb and create an interpolating function so that the PPPC functions don’t have to be called over and over again. (I believe this saves a lot of time since it doesn’t have to re-interpolate?) Either way, this has been shown to speed things up nicely. Points are spaced evernly on a log scale. I found that 100 points work well, Anthony found that fewer points are also ok.

minE and maxE is the range of photon energies over which we’d like to interpolate
numpoints is the number of points to sample, I use 100

p2ϕ2gamfast[primary_,mprim_,mchi_,mphi_,minE_,maxE_,numpoints_]:=Module[
{prim=primary,mSM=mprim,mχ=mchi,mϕ=mphi,mine=minE,maxe=maxE, samplepoints},

samplepoints=(10^(#))&/@Range[Log[10,minE],Log[10,maxE],((Log[10,maxE]-Log[10,minE])/numpoints)];

(* Output *)
Interpolation[{#,p2ϕ2gam[prim,mSM,mχ,mϕ,#]}&/@samplepoints]
] (* usage: f=b2gamfast[...]. f[(Subscript[E, γ])] gives the number of photons with E=(Subscript[E, γ]). *)



(* faster version that uses p2ϕ2gamquick, defined with "SymbolicProcessing"->0 *)
p2ϕ2gamfastest[primary_,mprim_,mchi_,mphi_,minE_,maxE_,numpoints_]:=Module[
{prim=primary,mSM=mprim,mχ=mchi,mϕ=mphi,mine=minE,maxe=maxE, samplepoints},

samplepoints=(10^(#))&/@Range[Log[10,minE],Log[10,maxE],((Log[10,maxE]-Log[10,minE])/numpoints)];

(* Output *)
Interpolation[{#,p2ϕ2gamquick[prim,mSM,mχ,mϕ,#]}&/@samplepoints]
] (* usage: f=b2gamfast[...]. f[(Subscript[E, γ])] gives the number of photons with E=(Subscript[E, γ]). *)

Boosted Spectra: 2-to-3

This section is based on 2to3Amplitude.nb. Here we assume annihilation into three mediators is the dominant annihilation mode. We’ve worked out the amplitude and everything. See the other notebook for testing.

The Amplitude

The amplitude is scaled by (g^3). Note that there’s a factor of 2 to account for C12 = C21, etc.

p[Ei_,Mi_]:=(Sqrt[(Ei^2)-(Mi^2)])
C12[E1_,E2_,cos_,Mϕ_,m_]:=E1 E2 -p[E1,Mϕ]p[E2,Mϕ] cos
C13[E1_,E2_,cos_,Mϕ_,m_]:=2m E1 - ((Mϕ)^2)-C12[E1,E2,cos,Mϕ,m]
C23[E1_,E2_,cos_,Mϕ_,m_]:=2m E2- ((Mϕ)^2)-C12[E1,E2,cos,Mϕ,m]
D12[E1_,E2_,cos_,Mϕ_,m_]:=(((Mϕ)^2)-2 m E1)(((Mϕ)^2)-2 m E2)
D13[E1_,E2_,cos_,Mϕ_,m_]:=(((Mϕ)^2)-2 m E1)(((Mϕ)^2)-2 m (2 m - E1 - E2))
D23[E1_,E2_,cos_,Mϕ_,m_]:=(((Mϕ)^2)-2 m E2)(((Mϕ)^2)-2 m (2 m - E1 - E2))

M[E1_,E2_,cos_,Mϕ_,m_]:=-(2m) 2(((C12[E1,E2,cos,Mϕ,m])/(D12[E1,E2,cos,Mϕ,m]))+((C13[E1,E2,cos,Mϕ,m])/(D13[E1,E2,cos,Mϕ,m]))+((C23[E1,E2,cos,Mϕ,m])/(D23[E1,E2,cos,Mϕ,m]))) (* times (g^3) *)

Kinematics notes, not used

These correspond to acceptable of cosine theta. First, the range of (Subscript[E, 2]) depends on the value of (Subscript[E, 1]) and the CM energy 2m as follows:

((Subscript[E, 1])^max)=(((2m))/2)-((3((Mϕ)^2))/(2 (2m)))

let ξ(E1,E2,cθ)=E1+E2+(Sqrt[(p1^2)+(p2^2)+p1 p2 cθ + ((Mϕ)^2)]),
the total energy of the final states

if ξ(E1,0,-1) > (2m), then the limits are the solutions to 
ξ(E1,E2,-1)=2m

otherwise, the limits are solutions to 
ξ(E1,E2,+1)=2m
ξ(E1,E2,-1)=2m

ξ[E1_,E2_,cθ_]:=E1+E2+(Sqrt[((E1^2)-((Mϕ)^2))+((E2^2)-((Mϕ)^2))+2(Sqrt[((E1^2)-((Mϕ)^2))])(Sqrt[(E2^2)-((Mϕ)^2)])cθ +((Mϕ)^2)])

Define some useful functions

Etot[E1_,E2_,cθ_,Mϕ_]:=E1+E2+(Sqrt[((E1^2)-((Mϕ)^2))+((E2^2)-((Mϕ)^2))+2(Sqrt[((E1^2)-((Mϕ)^2))])(Sqrt[(E2^2)-((Mϕ)^2)])cθ +((Mϕ)^2)])
MaxE1[m_,Mϕ_]:=((2 m)/2)-((3 ((Mϕ)^2))/(2 2m))
IsCaseI[E1_,m_,Mϕ_]:=Etot[E1,Mϕ,-1,Mϕ]>2 m
MaxE2[E1_,m_,Mϕ_]:=E2x/.Solve[Etot[E1,E2x,-1,Mϕ]==2m,E2x]//Max
MinE2[E1_,m_,Mϕ_]:=If[IsCaseI[E1,m,Mϕ],
E2x/.Solve[Etot[E1,E2x,-1,Mϕ]==2m,E2x]//Min
,
E2x/.Solve[Etot[E1,E2x,+1,Mϕ]==2m,E2x]//Min]

(* Amplitude depends on Cos[θ] which is fixed by kinematics *)
costheta[E1_,E2_,m_,Mϕ_]:=((((2m-E1-E2))^2)-((E1^2)-((Mϕ)^2))-((E2^2)-((Mϕ)^2))-((Mϕ)^2))(1/(2 (Sqrt[(E1^2)-((Mϕ)^2)]) (Sqrt[(E2^2)-((Mϕ)^2)])))

Do Momentum Integral

Factor of 3 coming from the multiplicity of final states. Note that this is NOT dN/dE, this is d<sigma v> / dE.

Symbolic Processing ->0 trick: http://mathematica.stackexchange.com/questions/8768/methods-to-speed-up-numerical-ndsolve-nintegrate

dσvdE[E1_,m_,Mϕ_]:=3NIntegrate[((M[E1,E2,costheta[E1,E2,m,Mϕ],Mϕ,m])^2),{E2,MinE2[E1,m,Mϕ],MaxE2[E1,m,Mϕ]}]

(* slightly faster if I use SymbolicProcessing -> 0  *)
dσvdEfast[E1_,m_,Mϕ_]:=3NIntegrate[((M[E1,E2,costheta[E1,E2,m,Mϕ],Mϕ,m])^2),{E2,MinE2[E1,m,Mϕ],MaxE2[E1,m,Mϕ]},Method->{Automatic,"SymbolicProcessing"->0}]

Doing the cross section integral: not used in p3ϕ2gam fast

The above plots are not normalized---we ultimately want dN/dE. We’ll need to integrate them to normalize. Note that it’s hard to directly NIntegrate dNdE because Mathematica gets confused by the symolic limits of the E2 integral. We fix this by defining a new function which only takes numbers. See: http://stackoverflow.com/questions/18903010/mathematica-wont-nintegrate-a-conditional-function

dσvdEnum[E1_?NumericQ,m_?NumericQ,Mϕ_?NumericQ]:=dσvdE[E1,m, Mϕ]
(*NIntegrate[dσvdEnum[x,120,45],{x,45,MaxE1[120,45]}]//Timing*)

The total spectrum: nice, but not used... too slow

This might be a helpful tip: http://mathematica.stackexchange.com/questions/10533/nested-nintegrate

Nov 6 , the normalization is off: have to account for 3 final states. This info which is in dσvdEnum, but which is washed out when I normalized by Totalσv. Then I have to account for the fact that the boxnorm is normalized to 2 final state mediators.

p3ϕ2gam[primary_,mprimary_,mχ_,mϕ_,photonE_]:=Module[
{finalstate=primary,mSM=mprimary,mchi=mχ, mphi=mϕ,eγ=photonE,γ,γβ, E0, boxnorm, Totalσv,Ephi,EfsHi,EfsLo,InnerIntegrand,OuterIntegrand},
γ[Eϕ_]:=((Eϕ)/mphi);
γβ[Eϕ_]:=(Sqrt[(((((Eϕ)/mphi)))^2)-1]);
E0=mphi/2;
boxnorm[Eϕ_]:=(1/(2 (*normalize*)γβ[Eϕ](Sqrt[((((mphi/2)))^2)-(mSM^2)])));
Totalσv=NIntegrate[dσvdEnum[Ephi,mchi,mphi],{Ephi,mphi,MaxE1[mchi,mphi]}];
EfsHi[Eϕ_]:=γ[Eϕ] E0+γβ[Eϕ] (Sqrt[((((mphi/2)))^2)-(mSM^2)]);
EfsLo[Eϕ_]:=γ[Eϕ] E0-γβ[Eϕ] (Sqrt[((((mphi/2)))^2)-(mSM^2)]);
InnerIntegrand[Efs_?NumericQ]:=p2gam[finalstate,Efs,eγ];
OuterIntegrand[Ephi_?NumericQ]:=boxnorm[Ephi]NIntegrate[InnerIntegrand[Efs],{Efs,EfsLo[Ephi],EfsHi[Ephi]}];
((3 (* norm *))/(Totalσv))NIntegrate[OuterIntegrand[Ephi],{Ephi,mphi,MaxE1[mchi,mphi]}]

(* the versions below work, but they give a bunch of warnings... they end up being off at the 3rd decimal *)
(*(1/(Totalσv))NIntegrate[boxnorm[Ephi] p2gam[finalstate,Efs,eγ],{Ephi,mphi,MaxE1[mchi,mphi]},{Efs,EfsLo[Ephi],EfsHi[Ephi]}]*)

(*real version:*)
(*(1/(Totalσv))NIntegrate[boxnorm[Ephi]NIntegrate[p2gam[finalstate,Efs,eγ],{Efs,EfsLo[Ephi],EfsHi[Ephi]}],{Ephi,mphi,MaxE1[mchi,mphi]}]*)
]

2 to 3 Faster

First define a discrete version of p3ϕ2gam where we don’t need to calculate Totalσv or do a full numerical integra over Eϕ. Further, we’ll just appeal to p2ϕ2gamFast. The idea is that p3ϕ2gam can be approximated as the weighted average of many p2ϕ2gam functions where mχ is varied to match Eϕ.

Some testing: nϕbins = 5 is too small, nϕbins=20 agrees to the third decimal for a reasonable test value. From previous tests, Nγ should be at least 50.

p3ϕ2gamfaster[SMprimary_,mprimary_,mχ_,mϕ_,minEγ_,maxEγ_,nϕbins_,nγbins_]:=Module[
{primary=SMprimary,mSM=mprimary,mchi=mχ,mphi=mϕ,minEgam=minEγ,maxEgam=maxEγ,Nphi=nϕbins,Ngam=nγbins,
minEphi,maxEphi,ΔEϕ,Eϕmidpoints,dσvdEsamples,totalσv,WeightedSpectrum,samplepoints},
(* Identify points to sample in Eϕ *)
minEphi=mϕ;
maxEphi=MaxE1[mchi,mphi];
ΔEϕ=((maxEphi-minEphi)/Nphi);
Eϕmidpoints=((#[[1]]+#[[2]])/2)&/@Partition[Range[minEphi,maxEphi,ΔEϕ],2,1];

(* Sample dσv/dE, also do box approximation to intergal *)
dσvdEsamples={#,dσvdEfast[#,mchi,mphi]}&/@Eϕmidpoints;

totalσv=ΔEϕ Total[#[[2]]&/@dσvdEsamples]; 
(* note that totalσv includes the factor of 3 for the 3x2 final state SM particles vs. only 2 final state SM particles. I have to multiply by 3 again down the road. *)
WeightedSpectrum[Egam_]:=Total[(((3 (* 3x2 final particles *) #[[2]] ΔEϕ)/(totalσv))((p2ϕ2gamquick[primary,mSM,#[[1]],mphi,Egam])/(2 (* because p2ϕ2gam is normed to 2 mediators *))))&/@dσvdEsamples];
(*NB: one way to make this even faster is to define a p2ϕ2gamquickest that samples instead of NIntegrates *)


(* Identify log sampling in Eγ *)
samplepoints=(10^(#))&/@Range[Log[10,minEgam],Log[10,maxEgam],((Log[10,maxEgam]-Log[10,minEgam])/Ngam)];

(* Output an interpolated function *)
Interpolation[{#,WeightedSpectrum[#]}&/@samplepoints]
]

DC’s SaveIt and ReadIt functions

http://www.terpconnect.umd.edu/%7Edcurtin1/mathematica.html
See also DumpSave, http://mathematica.stackexchange.com/questions/1959/how-do-i-save-a-variable-or-function-definition-to-a-file
... but DC says this carries the risk of overwriting things because Mathematica saves the initial variable assignment, not just the object itself.

(* to avoid memory bloat *)
ClearMemory:=Module[{},Unprotect[In,Out];
Clear[In,Out];
Protect[In,Out];
ClearSystemCache[];];

SaveIt[filename_,expr_]:=Module[{output},output=Export[filename<>".dat",ToString[expr//InputForm],"String"];
ClearMemory;
output];
SaveIt[varnamestring_]:=Module[{output},output=Export[varnamestring<>".dat",ToString[ToExpression[varnamestring]//InputForm],"String"];
ClearMemory;
output];
ReadIt[filename_]:=Module[{output},output=ToExpression[Import[StringReplace[filename,".dat"->""]<>".dat","String"]];
ClearMemory;
output];

Output one data file per function

Makes it easier if I have to halt at some mid point.

Interpol3phi[SMprimary_,mprimary_,mχ_,mϕ_,minEγ_,maxEγ_,nϕbins_,nγbins_,filename_]:=Module[
{primary=SMprimary,mSM=mprimary,mchi=mχ,mphi=mϕ,minEgam=minEγ,maxEgam=maxEγ,Nphi=nϕbins,Ngam=nγbins,path=filename},
SaveIt[path<>"_"<>ToString[mchi]<>"_"<>ToString[mphi],{{primary,mSM,mchi,mphi,minEgam,maxEgam,Nphi,Ngam,path},p3ϕ2gamfaster[primary,mSM,mchi,mphi,minEgam,maxEgam,Nphi,Ngam]}]
]

Interpol2V[SMprimary_,mprimary_,mχ_,mϕ_,minEγ_,maxEγ_,nγbins_,filename_]:=Module[
{primary=SMprimary,mSM=mprimary,mchi=mχ,mphi=mϕ,minEgam=minEγ,maxEgam=maxEγ,Ngam=nγbins,path=filename},
SaveIt[path<>"_"<>ToString[mchi]<>"_"<>ToString[mphi],{{primary,mSM,mchi,mphi,minEgam,maxEgam,Ngam,path},p2ϕ2gamfastest[primary,mSM,mchi,mphi,minEgam,maxEgam,Ngam]}]
]