Import data and build interpolating functions

Needs["PlotLegends`"]
Off[NIntegrate::"ncvb",NIntegrate::"precw",Power::"infy",∞::"indet",NIntegrate::"inumr",NIntegrate::"inumri",NIntegrate::"slwcon",NIntegrate::"nlim",Table::"iterb"]
DataFolder =  NotebookDirectory[]<>"/data/"; (*insert the path to the folder where the FITS files are stored*)
energyPC = Import[DataFolder <>"energyPC1"  <>".fits","Data"];
lifetimePC = Import[DataFolder <>"lifetimePCA"  <>".fits","Data"];

(*read data for different energies and lifetimes*)
(*We use PCA for the Planck experiment*)
(*lifetime PCA is for ~30 MeV e+e- pairs, and is calibrated to match the MCMC result for decay of long-lifetime DM (see 1610.06933)*)
elecarray13=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,3,i]]},{i,1,41}];
photarray13=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,3,i]]},{i,42,82}];
elecarray135=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,4,i]]},{i,1,41}];
photarray135=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,4,i]]},{i,42,82}];
elecarray14=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,5,i]]},{i,1,41}];
photarray14=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,5,i]]},{i,42,82}];
elecarray145=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,6,i]]},{i,1,41}];
photarray145=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,6,i]]},{i,42,82}];
elecarray15=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,7,i]]},{i,1,41}];
photarray15=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,7,i]]},{i,42,82}];
elecarray27=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,8,i]]},{i,1,41}];
photarray27=Table[{energyPC [[1,1,1,i]],energyPC [[1,1,8,i]]},{i,42,82}];
lifetimearray=Table[{lifetimePC [[1,1,1,i]],lifetimePC [[1,1,4,i]]},{i,1,26}];
(*build interpolating functions*)
feffe13 = Interpolation[elecarray13];
feffp13 = Interpolation[photarray13];
feffe135 = Interpolation[elecarray135];
feffp135 = Interpolation[photarray135];
feffe14 = Interpolation[elecarray14];
feffp14 = Interpolation[photarray14];
feffe145 = Interpolation[elecarray145];
feffp145 = Interpolation[photarray145];
feffe15 = Interpolation[elecarray15];
feffp15 = Interpolation[photarray15];
feffe27 = Interpolation[elecarray27];
feffp27 = Interpolation[photarray27];
fefflife= Interpolation[lifetimearray];
feff2de=Interpolation[Join[Table[{{13,elecarray13[[i,1]]},elecarray13[[i,2]]},{i,1,Length[elecarray13]}],Table[{{13.5,elecarray135[[i,1]]},elecarray135[[i,2]]},{i,1,Length[elecarray135]}],Table[{{14,elecarray14[[i,1]]},elecarray14[[i,2]]},{i,1,Length[elecarray14]}],Table[{{14.5,elecarray145[[i,1]]},elecarray145[[i,2]]},{i,1,Length[elecarray145]}],
Table[{{15,elecarray15[[i,1]]},elecarray15[[i,2]]},{i,1,Length[elecarray15]}],
Table[{{27,elecarray27[[i,1]]},elecarray27[[i,2]]},{i,1,Length[elecarray27]}]  ],InterpolationOrder->1];
feff2dp=Interpolation[Join[Table[{{13,photarray13[[i,1]]},photarray13[[i,2]]},{i,1,Length[photarray13]}],Table[{{13.5,photarray135[[i,1]]},photarray135[[i,2]]},{i,1,Length[photarray135]}],Table[{{14,photarray14[[i,1]]},photarray14[[i,2]]},{i,1,Length[photarray14]}],Table[{{14.5,photarray145[[i,1]]},photarray145[[i,2]]},{i,1,Length[photarray145]}],
Table[{{15,photarray15[[i,1]]},photarray15[[i,2]]},{i,1,Length[photarray15]}],
Table[{{27,photarray27[[i,1]]},photarray27[[i,2]]},{i,1,Length[photarray27]}]  ],InterpolationOrder -> 1];
(*make plots of the PC for long-lifetime decay to e+e- pairs, long-lifetime decay to photon pairs, and decay with variable lifetime to ~30 MeV e+e- pairs, and the 2D contour plots of the PCs as a function of energy and lifetime (for both photons and e+e- pairs)*)
Plot[feffe27[x],{x,4,12},PlotStyle->{Thick,Blue},Frame->True,FrameLabel->{{"",""},{"(Subscript[Log, 10])(Energy(eV))","Electron"}},LabelStyle->{Bold,12}]
Plot[feffp27[x],{x,4,12},PlotStyle->{Thick,Blue},Frame->True,FrameLabel->{{"",""},{"(Subscript[Log, 10])(Energy(eV))","Photon"}},LabelStyle->{Bold,12}]
Plot[fefflife[x],{x,18,13},PlotStyle->{Thick,Blue},Frame->True,FrameLabel->{{"τ/(Subscript[f, x]) (sec)",""},{"(Subscript[Log, 10])(lifetime(sec))",""}},LabelStyle->{Bold,12}]
ContourPlot[feff2de[x,y],{x,13,17},{y,4,12},FrameLabel->{{"(Subscript[Log, 10])(Energy(eV))",""},{"(Subscript[Log, 10])(lifetime(sec))","2D plot for Electron"}},LabelStyle->{Bold,12}]
ContourPlot[feff2dp[x,y],{x,13,17},{y,4,12},FrameLabel->{{"(Subscript[Log, 10])(Energy(eV))",""},{"(Subscript[Log, 10])(lifetime(sec))","2D plot for Photon"}},LabelStyle->{Bold,12}]

Example: Planck constraints on decay to muons

(*Muon spectrum for DM >> muon mass*)

MuonElectronSpectrum[eng_,muoneng_] := (12/muoneng)*((1/4)(1 - (eng/muoneng)^2) - (1/9)(1 - (eng/muoneng)^3)) 
(*calculate effective f *)
FForMuons[lifetime_,logmass_] :=2 Log[10]NIntegrate[(10^logeng)^2 *MuonElectronSpectrum[10^logeng,10^logmass/2]*feff2de[lifetime,Log10[10^logeng - 511000]],{logeng,Log10[512000],Log10[10^logmass/2]}] / 10^logmass
(*bound on decaying mass fraction fX for different lifetimes*)
taubound13=  Table[{10^logmass/10^9, (10^13)/(FForMuons[13,logmass]* 2.6 *10^25/feffe13[7.5] * fefflife[13]/ fefflife[18] )},{logmass,9,13,0.2}];
taubound14=  Table[{10^logmass/10^9, (10^14)/(FForMuons[14,logmass]* 2.6 *10^25/feffe14[7.5] * fefflife[14]/ fefflife[18] )},{logmass,9,13,0.2}];
taubound15=  Table[{10^logmass/10^9, (10^15)/(FForMuons[15,logmass]* 2.6 *10^25/feffe15[7.5] * fefflife[15]/ fefflife[18] )},{logmass,9,13,0.2}];
taubound16=  Table[{10^logmass/10^9, (10^16)/(FForMuons[16,logmass]* 2.6 *10^25/feff2de[16,7.5] * fefflife[16]/ fefflife[18] )},{logmass,9,13,0.2}];
taubound17=  Table[{10^logmass/10^9, (10^17)/(FForMuons[17,logmass]* 2.6 *10^25/feff2de[17,7.5] * fefflife[17]/ fefflife[18] )},{logmass,9,13,0.2}];
taubound18=  Table[{10^logmass/10^9, (10^18)/(FForMuons[18,logmass]* 2.6 *10^25/feff2de[18,7.5] * fefflife[18]/ fefflife[18] )},{logmass,9,13,0.2}];

ListLogLogPlot[{taubound13,taubound14,taubound15,taubound16,taubound17,taubound18},Joined -> True,PlotRange->{{1,(10^4)},{(10^(-12)),(10^(-4))}},Frame->True,PlotStyle->{{Thick,Black},{Thick,Dashed,Red},{Thick,Yellow},{Thick,Green},{Thick,Blue},{Thick,Purple}},FrameLabel->{{"fraction (Subscript[f, x])",""},{"DM(GeV)",""}},LabelStyle->{Bold,18},PlotLegend->{Style["(10^13)sec"],Style["(10^14)sec"],Style["(10^15)sec"],Style["(10^16)sec"],Style["(10^17)sec"],Style["(10^18)sec"]},LegendShadow->None,LegendSize->0.6,LegendPosition->{0.8,0.1},ImageSize->600]