header = "/int/users/paulrom/" <> oheader

file = Import[header <> "freezeout.dat"];

xL = Length[file];

bigN=1;

rfunc[x_] = 
  Normal[Fit[file[[All, 1]], {Sqrt[xL^2-x^2],1, x, x^2, x^3, x^4,x^5, x^6,x^7,x^8}, x]];

tfunc[x_] = 
  Normal[Fit[file[[All, 2]], {1, x}, x]];


Print["xL=",xL];

dtdxfunc[x_] = D[tfunc[x], x];

drdxfunc[x_] = D[rfunc[x], x];

utfunc = Interpolation[file[[All, 3]]];

urfunc = Interpolation[file[[All, 4]]];

pirf=Interpolation[file[[All,5]]];

pief = Interpolation[file[[All, 6]]];

A:=Sqrt[p^2+m^2]/T utfunc[x] n

B:=p/T urfunc[x] n

mt:=Sqrt[p^2+m^2]

v:=urfunc[x]/utfunc[x]

L1:=2 p mt v pirf[
      x] (mt (BesselK[0,A]+BesselK[2,A]) BesselI[1,B] drdxfunc[x]-
        p dtdxfunc[x] BesselK[1,A] (BesselI[2,B]+BesselI[0,B]))

L2:=-p^2 pirf[
      x] (2 mt BesselK[1,A] BesselI[2,B] drdxfunc[x]-
        p dtdxfunc[x] BesselK[0,A] (BesselI[3,B]+BesselI[1,B]))

L3:=-v^2 mt^2 pirf[
      x] (mt/2 (BesselK[3,A]+3BesselK[1,A]) BesselI[0,B] drdxfunc[x]-
        p dtdxfunc[x] (BesselK[0,A]+BesselK[2,A]) BesselI[1,B])

L4:=-v^2 p^2 pirf[
      x] (mt BesselK[1,A](BesselI[0,B]-BesselI[2,B]) drdxfunc[x]-
        p dtdxfunc[x]/2 BesselK[0,A](BesselI[1,B]-BesselI[3,B]))

L5:=mt^2 pief[
      x](mt 2 BesselK[1,A] BesselI[0,B] drdxfunc[x]-
        p dtdxfunc[x] 2 BesselK[0,A] BesselI[1,B])

L6:=p^2 pief[
      x] (mt BesselK[1,A] (BesselI[0,B]-BesselI[2,B]) drdxfunc[x]-
        p dtdxfunc[x]/2 BesselK[0,A](BesselI[1,B]-BesselI[3,B]))

L7:=-mt^2 pief[
      x](mt/2 (BesselK[3,A]+3 BesselK[1,A]) BesselI[0,B] drdxfunc[x]-
        p dtdxfunc[x] (BesselK[0,A]+BesselK[2,A]) BesselI[1,B])

Ifunc[pp_,mm_,TT_,
      xx_,nn_]:=((L1+L2+L3+L4+L5+L6+L7))/.x\[Rule]xx/.T\[Rule]TT/.m\[Rule]
        mm/.p\[Rule]pp/.n->nn

ptarr={};

Do[ptarr=Append[ptarr,i/40.+0.01],{i,0,99}];


dNdp[p_, m_, T_,sig_] := -
  2/(2 Pi)^2 5.06842^3NIntegrate[Sum[(sig)^(n+1)
      rfunc[x] tfunc[x](Sqrt[p^2 + m^2] BesselI[0, p/T urfunc[x] n] BesselK[1, 
                Sqrt[p^2 + m^2]/T utfunc[x] n] drdxfunc[x]- 
            dtdxfunc[x] p BesselI[1, p/T urfunc[x] n] BesselK[0, 
							      Sqrt[p^2 + m^2]/T*utfunc[x] n]),{n,1,bigN}], {x, 1, xL}, MinRecursion->4, MaxRecursion -> 20];

ddNdp[p_, m_, T_,sig_] := -
  1/(2 Pi)^2 5.06842^3 1/2/T^2NIntegrate[
					 rfunc[x] tfunc[x]Sum[(sig)^(n+1) Ifunc[p, m, T, x,n],{n,1,bigN}], {x, 1, xL}, MinRecursion->4,MaxRecursion -> 20];

genidspec[m_,sig_,T_]:=Module[{res},
			      res=Table[{ptarr[[i]], dNdp[ptarr[[i]], m, T,sig]}, {i,1,Length[ptarr]}];
			      res];


genspec[m_,sig_,T_]:=Module[{res},
			      res=Table[{ptarr[[i]], dNdp[ptarr[[i]], m, T,sig]+ddNdp[ptarr[[i]],m,T,sig]}, {i, 1, Length[ptarr]}];
			    res];




model[T_, AA_, m0_, pt_] := 
Exp[(-Sqrt[pt^2 + m0^2] + m0)/Sqrt[T^2]]/T/(T + m0)*AA;

diff[T_, AA_, m0_, MIN_, MAX_] := If[AA>0,
  Sum[If[(dat[[i, 1]] > MIN) && (dat[[i, 1]] < MAX), (Log[dat[[i, 2]]] - 
						      Log[model[T, AA, m0, dat[[i, 1]]]])^2, 0], {i, 1, Length[dat]}],10^6];

bestfit[m0_,MIN_,MAX_]:=Module[{res,mT,mA},
			       bf=FindMinimum[diff[T,AA,m0,MIN,MAX],{T,0.1,0.2},{AA,1,2}];
			       mT=T/.bf[[2]];
			       mA=AA/.bf[[2]];
			       Print["Best fit found for T=",mT," and A=",mA];
			       res=Table[{p,model[mT,mA,m0,p]},{p,MIN,MAX,0.05}];
res]
