In[1]:= data = {{3501, 19}, {3502, 48}, {3503, 145}, {3504,
318}, {3505, 554}, {3506, 1017}, {3507, 1405}, {3508,
1829}, {3509, 2132}, {3510, 2255}, {3511, 1824}, {3512,
1438}, {3513, 971}, {3514, 548}, {3515, 244}, {3516, 136}, {3517,
52}, {3518, 12}};

In[2]:= G[x_, a_, b_, c_] := a*Exp[-(((x) - b)^2/(2 c^2))];

In[3]:= NLM = NonlinearModelFit[data, G[q, a, b, c]), {a, b, c}, q, Method -> NMinimize] // Normal


Out[3]= 2186.87 E^(-0.067269 (-3509.5 + q)^2)

In[4]:= NLM["ParameterTable"]

Out[4]= (2186.87 E^(-0.067269 (-3509.5 + q)^2))["ParameterTable"]