I don't know how you constructed that figure so I'm ignoring it.

Mathematica has a command (DistributionFitTest) that should test for any distribution but it either takes too long or doesn't finish for a Poisson when the counts are very high as in your example. So here is a brute force approach.

data = {812, 827, 782, 841, 839, 803, 800, 818, 790, 838, 783, 768, 
   799, 854, 805, 832, 786, 838, 802, 805, 825, 784, 797, 808, 846, 
   809, 814, 857, 884, 855, 800, 855, 818, 865, 801, 823, 812, 861, 
   807, 830, 813, 827, 819, 796, 870, 825, 827, 872, 770, 849, 813, 
   828, 839, 804, 805, 833, 833, 871, 856, 819, 787, 819, 852, 786, 
   864, 841, 799, 804, 812, 831, 830, 810, 863, 856, 794, 808, 794, 
   838, 748, 757, 866, 778, 823, 818, 830, 794, 831, 838, 853, 836, 
   821, 836, 797, 863, 811, 867, 803, 821, 864, 819};

(* Determine bins which would hold approximately equal numbers of 
   observations if the distribution was Poisson *)
dist = PoissonDistribution[Mean[data] // N];
nBins = 12;
bins = Join[{-1}, Table[InverseCDF[dist, i/nBins], {i, nBins - 1}], {\[Infinity]}] // Quiet
(* {-1,783,794,803,810,816,822,828,834,841,850,862,\[Infinity]} *)

(* Observed counts *)
observed = Table[Length[Select[data, bins[[i]] < # <= bins[[i + 1]] &]], {i, nBins}]
(* {7,8,11,10,7,9,8,8,10,2,9,11} *)

(* Expected counts *)
n = Length[data];
expected = Table[n (CDF[dist, bins[[i + 1]]] - CDF[dist, bins[[i]]]) // N, {i, nBins}]
(* {8.86592, 7.98339, 9.1535, 8.54294, 8.01289, 8.31321, 8.25561, 
7.84997, 8.25645, 8.72493, 8.04893, 7.99226} *)

(* Chisquare statistic *)
chisq = Total[(observed - expected)^2/expected]
(* 8.00515 *)

(* P-value *)
pValue = 1 - CDF[ChiSquareDistribution[nBins - 2], chisq]
(* 0.628334 *)