Integration of Interpolated function take an unacceptable amount of time

The code in the question takes about 40 seconds on my computer. Turning off SymbolicProcessing, as suggested by Tim Laska in a comment, reduces the run time to about 16 seconds. Although reducing MaxRecursion further reduces the run time, it also generates error messages that the integration is not converging to sufficient accuracy. (Whether this matters from a practical point of view depends on how the results will be used.)

Instead, let us take advantage of the fact that NIntegrate in Plot repeatedly performs the same integral over the range {2.66*^17, 3.78*^18} to break the integral into two pieces.

int = NIntegrate[ElectronNumberDensity[eta], {eta, 2.66*^17, 3.78*^18}]
Plot[int + NIntegrate[ElectronNumberDensity[eta], {eta, η, 2.66*^17}, 
    Method -> {Automatic, "SymbolicProcessing" -> 0}], {η, 1.47*^17, 2.66*^17}, 
    ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

which produces the desired plot in under 5 seconds.

enter image description here

The plot agrees well with that produced by the code in the question.

Addendum

Even faster is to use NDSolveValue instead of NIntegrate.

NDSolveValue[{g'[eta] == -ElectronNumberDensity[eta], g[2.66*^17] == int}, 
    g, {eta, 1.47*^17, 2.66*^17}]
Plot[%[η], {η, 1.47*^17, 2.66*^17}, ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

which produces the same curve in less than a second.


This

c = 2.99792*10^5;
A = 3.87624*10^-14;
FreeElectronFractionData = {...YourData...}; 
Clear[\[Eta]]; (* just in case there was a prior assignment to Eta *)
redShift = 6.64*^18^2/((c - Sqrt[c]*Sqrt[c - 2.*A*\[Eta]])/A)^2 - 1.;
FreeElectronFraction = Interpolation[FreeElectronFractionData, InterpolationOrder -> 1];
ElectronNumberDensity[\[Eta]_] := FreeElectronFraction[redShift ]*1.42*^-7*(1. +redShift^3);
Plot[NIntegrate[ElectronNumberDensity[eta], {eta, \[Eta], 3.78*^18}, MaxRecursion -> 15], {\[Eta], 1.47*^17, 2.66*^17}]
ListPlot[FreeElectronFractionData]

completes in a less than ten seconds, mostly I believe because I changed

FreeElectronFraction := Interpolation[...]

to

FreeElectronFraction = Interpolation[...]

which means the calculation of the interpolating function will be done only once and the result will then be repeatedly used instead of calculating the interpolation function thousands of times.

I also moved the location of your calculation of redShift which may or may not have been necessary, but I was getting recursion errors and fixed those by moving the redShift calculation.

Please check all this very carefully to make certain the result is correct.