MapQ = Compile[{{r0, _Real, 2}, {r, _Real, 2}, xc, yc, {w, _Real, 1}},
   Block[{ M = ({
         {0., 0., 0., 0.},
         {0., 0., 0., 0.},
          {0., 0., 0., 0.},
          {0., 0., 0., 0.}
        }), x0, y0, z0, x, y, z, w1, w2, r2, c1, c2, c3, c4, x0x, y0y,
      z0z, sys, [Theta], temp},
    Do[
     x0 = r0[[i, 1]] - xc; x = r[[i, 1]] - xc;
          y0 = r0[[i, 2]] - yc; y = r[[i, 2]] - yc;
        z0 = r0[[i, 3]]; z = r[[i, 3]];

     (*// matrix elements//*)
                 w1 = w[[i]];
     w2 = 2*w1;
     r2 = x^2 + y^2 + z^2 + x0^2 + y0^2 + z0^2;
     x0x = x0 x;
     y0y = y0 y;
     z0z = z0 z;
     c1 = (r2 - 2 (x0x + y0y + z0z));
     c2 = (r2 + 2 (-x0x + y0y + z0z));
     c3 = (r2 + 2 (x0x - y0y + z0z));
     c4 = (r2 + 2 (x0x + y0y - z0z));
     M[[1, 1]] += w1 c1;
     M[[1, 2]] += w2 (y z0 - z y0);
     M[[1, 3]] += w2 (-x z0 + z x0);
     M[[1, 4]] += w2 (x y0 - y x0);
     M[[2, 2]] += w1 c2;
     M[[2, 3]] -= w2 (x y0 + y x0);
     M[[2, 4]] -= w2 (x z0 + z x0);
     M[[3, 3]] += w1 c3;
     M[[3, 4]] -= w2 (y z0 + z y0);
     M[[4, 4]] += w1 c4;

     , {i, Length[r0]}];
    M = ({
       {M[[1, 1]], M[[1, 2]], M[[1, 3]], M[[1, 4]]},
       {M[[1, 2]], M[[2, 2]], M[[2, 3]], M[[2, 4]]},
       {M[[1, 3]], M[[2, 3]], M[[3, 3]], M[[3, 4]]},
       {M[[1, 4]], M[[2, 4]], M[[3, 4]], M[[4, 4]]}
      })

    ], RuntimeAttributes -> {Listable}, Parallelization -> True,
   CompilationTarget -> "C"];



AngleQ[M_] := Block[{sys, [Theta]},
       sys = Eigensystem[M];
  If[sys[[2, 4, 1]] < 0, sys[[2, 4]] = -sys[[2, 4]], None];
       sys[[2, 4]] = Normalize[sys[[2, 4]]];
  [Theta] = 2 ArcCos[sys[[2, 4, 1]]];
  If[sys[[2, 4, 4]]/Sin[[Theta]/2] < 0, [Theta] = -[Theta], None];
  [Theta]
  ]