Quit[]

Import data, define energy-dependent effective f function

DataFolder = "data/"; (*insert the path to the folder where the FITS files are stored*)
Species = "elec"; (*adjust species to "phot" or "elec" here*)
AssumedIonizationHistory = "3keV"; (*set this to "3keV" or "SSCK", use "3keV" if you want to compare to Planck constraints*)
FEffFile = DataFolder <>"feff_summary_" <> Species <> "_" <> AssumedIonizationHistory <>".fits";
TempGrid=Import[FEffFile,"Data"];
{redshift, weightfn, log10energy, feffbest} = { TempGrid[[1,1,1]],TempGrid[[1,1,2]], TempGrid[[1,1,3]],TempGrid[[1,1,4]]};
fefftable = Table[{log10energy[[i]],feffbest[[i]]},{i,1,40}];
feff = Interpolation[fefftable];
Show[ListPlot[fefftable],Plot[feff[x],{x,3,13}]]

Example of computing effective f for DM annihilating to muons

MuonElectronSpectrum[eng_,muoneng_] := (12/muoneng)*((1/4)(1 - (eng/muoneng)^2) - (1/9)(1 - (eng/muoneng)^3)) 
(*Spectrum of electrons (dN/dE) produced by decay of muons from DM annihilation/decay. Approximate form okay for mDM >> muon mass. Normalized per DM particle.*)

(*First set Species to "elec" and run the "Import data" panel*)
FForMuons[mass_] :=NIntegrate[Exp[logeng]^2 MuonElectronSpectrum[Exp[logeng],mass]feff[Log10[Exp[logeng] - 511000]],{logeng,Log[512000],Log[mass]}] / mass(*DM mass must be in eV*)

Example of computing Planck constraints on a DM model

ftable=Table[{10^logmass,FForMuons[10^logmass]},{logmass,9,13,0.5}]

sigmavband = Table[{ftable[[i,1]]/10^9,4.1*10^(-28)*(ftable[[i,1]]/10^9)/ftable[[i,2]]},{i,1,Length[ftable]}]

ListLogLogPlot[sigmavband,Joined -> True,PlotRange -> {10^(-27),10^(-21)},Filling -> Top,AxesLabel -> {"DM mass (GeV)","<σv> ((cm^3)/s)"}]

Using the universal DM annihilation weighting function

(*Once an f-grid has been derived, we can determine the relevant f_eff value for any photon or electron energy (e.g. for comparison with the constraints presented by Planck in 1502.01589), by integrating over the weighting function derived by a principal component analysis. These results are summarized in the FITS files and their use is demonstrated in the "Simple effective f analysis" part of the notebook above.*)
WeightingFunction = Interpolation[Table[{redshift[[i]],weightfn[[i]]},{i,1,Length[weightfn]}]];
FFunction[z_] := Exp[-z/600]; (*Artificial placeholder f(z) function, for true f(z) functions read in results from supplemental materials for Paper I.*)
effectivef = NIntegrate[WeightingFunction[Exp[logz]] FFunction[Exp[logz]],{logz,Log[10],Log[1200]}]

(*Alternate version, which will only work if EnergyAbsorptionCalculator_General.nb has already been executed to obtain AggregateTotalF.*)
WeightingFunction = Interpolation[Table[{redshift[[i]],weightfn[[i]]},{i,1,Length[weightfn]}]];
FFunction = Interpolation[Flatten[AggregateTotalF,1]];
EffectiveF[energy_] := NIntegrate[WeightingFunction[Exp[logz]] FFunction[Log10[energy],Log10[Exp[logz]]],{logz,Log[10],Log[1200]}]
(*Note that IonizationPrescription (set in defining AggregateTotalF) should be set depending on the constraint to which one is comparing; for comparison to the Planck constraints, IonizationPrescription = "3keV" should be set.*)
ListPlot[Table[{log10energy,EffectiveF[10^log10energy]},{log10energy,3.5,13,0.1}],AxesLabel -> {"Log10(Energy (eV))","f_eff"}]