Probability Sphere Origin in Tetrahedron

//make outside sphere for vis
float scale = ch('scale');
int origpt = addpoint(0, 0);
setpointattrib(0, "Alpha", origpt, 0.3);
int sphere = addprim(0, "sphere", origpt);
matrix scalem = maketransform(0, 0, {0,0,0}, {0,0,0}, scale, 0);
setprimintrinsic(0, "transform", sphere, matrix3(scalem));

//make a little sphere at the origin
int vispt = addpoint(0, 0);
int vissphere = addprim(0, "sphere", vispt);
matrix viss = maketransform(0, 0, {0,0,0}, {0,0,0}, scale * 0.02, 0);
setprimintrinsic(0, "transform", vissphere, matrix3(viss));

//loop generate a tetrahedron and calc if origin is inside it
int iter = chi('iterations');
iter = 15000;
int count = 0;
int make_geo = chi('make_geo');
int pt1, pt2, pt3, pt4;
for(int i = 0; i < iter; i += 1){
    float seed = ch('seed') + i;
    seed = 0.283 + i;


    vector pos1 = normalize(fit01(vector(rand(2.4 + seed)), -1, 1)) * scale;
    vector pos2 = normalize(fit01(vector(rand(3.4 + seed)), -1, 1)) * scale;
    vector pos3 = normalize(fit01(vector(rand(4.4 + seed)), -1, 1)) * scale;
    vector pos4 = normalize(fit01(vector(rand(5.4 + seed)), -1, 1)) * scale;
    if(make_geo == 1){
        pt1 = addpoint(0, pos1);
        pt2 = addpoint(0, pos2);
        pt3 = addpoint(0, pos3);
        pt4 = addpoint(0, pos4);

        int prim1 = addprim(0, "poly", pt1, pt2, pt3);
        int prim2 = addprim(0, "poly", pt2, pt3, pt4);
        int prim3 = addprim(0, "poly", pt3, pt4, pt1);
        int prim4 = addprim(0, "poly", pt4, pt1, pt2);
    }
    int intersect3D_RayTriangle( vector R[]; vector T[]){
        vector    u, v, n;              // triangle vectors
        vector    dir, w0, w;           // ray vectors
        float     r, a, b;              // params to calc ray-plane intersect

        // get triangle edge vectors and plane normal
        u = T[1] - T[0];
        v = T[2] - T[0];
        n = cross(u, v);              // cross product
        if (n == {0,0,0})             // triangle is degenerate
            return -1;                  // do not deal with this case

        dir = R[1] - R[0];              // ray direction vector
        w0 = R[0] - T[0];
        a = -dot(n,w0);
        b = dot(n,dir);
        if (abs(b) < 0.000001) {     // ray is  parallel to triangle plane
            if (a == 0)                 // ray lies in triangle plane
                return 2;
            else return 0;              // ray disjoint from plane
        }

        // get intersect point of ray with triangle plane
        r = a / b;
        if (r < 0.0)                    // ray goes away from triangle
            return 0;                   // => no intersect
        // for a segment, also test if (r > 1.0) => no intersect

        vector I = R[0] + r * dir;            // intersect point of ray and plane
        //addpoint(0, I);

        // is I inside T?
        float    uu, uv, vv, wu, wv, D;
        uu = dot(u,u);
        uv = dot(u,v);
        vv = dot(v,v);
        w = I - T[0];
        wu = dot(w,u);
        wv = dot(w,v);
        D = uv * uv - uu * vv;

        // get and test parametric coords
        float s, t;
        s = (uv * wv - vv * wu) / D;
        if (s < 0.0 || s > 1.0)         // I is outside T
            return 0;
        t = (uv * wu - uu * wv) / D;
        if (t < 0.0 || (s + t) > 1.0)  // I is outside T
            return 0;

        return 1;                       // I is in T
    }


    vector R[] = {{0,0,0}, {0,- 1,0}};
    vector T[] = array(pos1, pos2, pos3);

    int hit = intersect3D_RayTriangle(R, T);

    T = array(pos2, pos3, pos4);
    hit += intersect3D_RayTriangle(R, T);

    T = array(pos3, pos4, pos1);
    hit += intersect3D_RayTriangle(R, T);

    T = array(pos4, pos1, pos2);
    hit += intersect3D_RayTriangle(R, T);
    if(hit % 2 != 0){
        if(make_geo >0){
            setpointattrib(0, "Cd", pt1, {1, 0, 0});
            setpointattrib(0, "Cd", pt2, {1, 0, 0});
            setpointattrib(0, "Cd", pt3, {1, 0, 0});
            setpointattrib(0, "Cd", pt4, {1, 0, 0});
        }
        setpointattrib(0, "Cd", vispt, {1, 0, 0});
        count += 1;
    }

}
@result = float(count) / float(iter);
@ref = 1.0/8.0;