fitting - Confidence Regions for NonLinearFit with “complicated” integral model

I have a set of data $y_i (z_i)$ with errors $\Delta y_i$, that can be found here

Here

Thus first, I call the data

dataSN = Import["Cat3_0mod.txt", "Table"];
sigdata = dataSN[[All, 4]];
lendata = Length[zdata];

fitSN = Table[{dataSN[[i, 2]], dataSN[[i, 3]]}, {i, 1, lendata}];

I want to fit these data to a complicated model

$$ Y(z) = 5\log \Bigl((1+z)\int_0^z \frac{dz'}{a+bz'+cz'^2}\Bigr)+25; $$

With help of a previous post NonLinearFit with "complicated" integral model

I used 'NonLinerFit'

ModelSN1[z_] = 
  Simplify[5 Log[10, (1 + z) Integrate[1/(a + b x + c x^2), {x, 0, z}]] + 25, 
    Assumptions -> {z > 0, -b^2 + 4 a c > 0}];

fit = 
  NonlinearModelFit[fitSN, {ModelSN1[z], -b^2 + 4 a c > 0}, 
   {{a, 0.000254}, {b, .000016}, {c, .000016}}, z, Weights -> 1/(sigdata)^2];

Even the fit looks fine

The reality is that the fit is complex

I Have some problems making the fit

1. Even I put conditions over the parameters $a,b,c$ the fitted model is complex. I tried to put as condition that the whole term inside $\log$ should be positive, but it does not work.


2. 
   

   I want to make confidence regions on this fit using

   
   

    band1SN[z_] = fit["MeanPredictionBands", ConfidenceLevel -> 0.68]; 
   
    band2SN[z_] = fit["MeanPredictionBands", ConfidenceLevel -> 0.95];
   

but Mathematica does not make it, an error related with the fact of assuming a constrained model appears. Of course, if the fit is complex, the confidence regions could not be calculated. I appreciate your help.