dateBH=Import["C:\Users\Laurentiu\Google Drive\Proiecte\FractalSky\2micron\2RedshiftLimit\BHCatalog\DataBHFilteredFluxRedshiftNew.csv"];
Length[dateBH]

alldistances=Flatten[Reap[Do[Sow[dateBH[[j,15]]],{j,Length[dateBH]}]][[2]],1];
Max[alldistances]

nrintervale=11;
masedebh=Range[Length[dateBH]];Do[masedebh[[i]]=dateBH[[i,16]],{i,Length[dateBH]}];
numar=Table[0,{nrintervale}];distance=numar;
Min[masedebh];maxdemase=2.6*10^8;alfacoef=2.843431601453441`*^11;betacoef=1.9054888536775378`;
intervale=Sort[Drop[Table[10^i,{i,Log[10,Min[masedebh]],Log[10,Max[masedebh]],(Log[10,Max[masedebh]]-Log[10,Min[masedebh]])/(nrintervale)}],-1],Greater]
listaredshift=intervale;
Do[listaredshift[[i]]=zz/.Solve[intervale[[i]]==alfacoef*zz^betacoef,zz],{i,Length[intervale]}]
distanceofcut=Sort[Flatten[listaredshift]*2.99792458*10^5/73,Greater]
Do[If[intervale[[i]]>maxdemase,distanceofcut[[i]]=Max[alldistances]],{i,Length[distanceofcut]}]
Do[Do[If[dateBH[[j,16]]⩾intervale[[i]],If[dateBH[[j,15]]⩽distanceofcut[[i]],numar[[i]]++]],{j,Length[dateBH]}],{i,Length[intervale]}];
numar
newnumar=numar;
Do[newnumar[[i]]=numar[[i]]+Total[Drop[numar,i]],{i,Length[numar]}]
(*numar=newnumar*)
distance=distanceofcut

deprintat=Transpose[N[{intervale,numar/((4/3)Pi*(distance^3))}]]

deplotateroriup=deprintat;error=Range[Length[deprintat]];
Flatten[Reap[Do[error[[i]]=((((((Sqrt[numar[[i]]])))^(-1))/(numar[[i]]))-(((4/3)Pi*(0.2)*3)/((4/3)Pi)))*(numar[[i]]/((4/3)Pi*((distance[[i]])^3)));Sow[error[[i]]],{i,Length[error]}]][[2]],1];
newerror=error;
Do[newerror[[i]]=0.434*error[[i]]/deprintat[[i,2]],{i,Length[error]}]
Do[deplotateroriup[[i]]=ListPlot[{{Log[10,deprintat[[i,1]]],Log[10,deprintat[[i,2]]]-newerror[[i]]},{Log[10,deprintat[[i,1]]],Log[10,deprintat[[i,2]]]+newerror[[i]]}},PlotStyle->Blue,Joined->True],{i,Length[deprintat]}]
N[error]
newerror

alfagasit=alfa/.FindFit[Log[10,deprintat],alfa+beta*s,{alfa,beta},s]
betagasit=beta/.FindFit[Log[10,deprintat],alfa+beta*s,{alfa,beta},s]

betagasittwo=-2;
alfagasittwo=12.9;
betagasittrei=-1;
alfagasittrei=4.76;
alfagasifour=-8.9;
betagasitfour=15.19;

alfagasittrei+betagasittrei*Log[10,10^7]
alfagasittwo+betagasittwo*Log[10,10^8]

poza=ListPlot[Log[10,Drop[deprintat,0]],Frame->True,FrameTicksStyle->Thick,FrameLabel->{"Log (Subscript[M, BH]) [(Subscript[M, sol])]", "Log ϕ((Subscript[M, BH])) [(Mpc^(-3))]"},BaseStyle->{FontFamily->"Times",FontSize->30},PlotStyle->{PointSize[0.01]},PlotRange->{{1,10},{-8,5}},ImageSize->950,Axes->False] ;
poza2=Plot[alfagasit+betagasit*s,{s,Log[10,intervale[[nrintervale-1]]]-.5,.9+Log[10,intervale[[2]]]}]
poza2=Plot[alfagasit+betagasit*s,{s,2.5,9.5}]
poza3=Plot[alfagasittwo+betagasittwo*s,{s,Log[10,1*10^7],Log[10,2*10^9]},PlotStyle->Directive[Red,Thickness[.005]]]
poza4=Plot[alfagasittrei+betagasittrei*s,{s,Log[10,1*10^6],Log[10,6*10^8]},PlotStyle->Directive[Blue,Thickness[.005]]]
poza5=Plot[betagasitfour+alfagasifour*Log[s],{s,Log[10,1*10^6],Log[10,6*10^8]},PlotStyle->Directive[Black,Thick]];
(*alfagasit=1.17475*10^9;
betagasit=-258.533;*)
poza6=Plot[Log[10,-3*betagasit^2*(E^-(10^s/alfagasit)^(1/3) (-(1/(3 10^s))-1/(6*(10^s)^(1/3)*alfagasit^(2/3))+1/(6*(10^s)^(2/3)*alfagasit^(1/3)))-ExpIntegralEi[-(10^s/alfagasit)^(1/3)]/(6 alfagasit))],{s,5.5,Log[10,2*10^9]},PlotRange->All];
poza1=Show[poza,Drop[deplotateroriup,0],poza2,poza3,poza4]

Export["/home/lau/Lucru/Proiecte/FractalSky/2micron/FinalList/BHMassF/PlotsBHMFBHSelected/BHMassFunctionCorrectedWithErrorsAndFitJune29DeVocuCutCorrectedForHubbleTypeWithFits.eps",poza1]

nlm=NonlinearModelFit[Log[10,deprintat],alfa+beta*s,{alfa,beta},s];
nlm[{"BestFit","FitResiduals","ParameterTable"}]

Export["/home/lau/Lucru/Proiecte/Publicatii/TheMassRatioSMBH/BHMassCountWithFit.eps",poza1]