interection.py

import math

def intersec(vx,vy,vz,px,py,pz):
    lightSpeed = 29979245800;
    #Intersection with cylinder with following radius and length (in cm as vertex)
    radius =  130;
    length = 300;

    energy = math.sqrt(px*px+py*py+pz*pz);
    #First work out intersection with cylinder (barrel)
    a = (px*px + py*py)*lightSpeed*lightSpeed/(energy*energy);
    b = 2*(vx*px + vy*py)*lightSpeed/energy;
    c = (vx*vx + vy*vy) - radius*radius;

    sqrt_disc = math.sqrt(b*b - 4*a*c);
    tCircle1 = (-b + sqrt_disc)/(2*a);
    tCircle2 = (-b - sqrt_disc)/(2*a);
    #If intersection in the past it doesn't count
    if tCircle1 < 0: tCircle1 = 1E9;
    if tCircle2 < 0: tCircle2 = 1E9;
    #If the intsersection occurs outside the barrel length it doesn't count
    zPosCircle1 = tCircle1*(pz/energy)*lightSpeed + vz;
    zPosCircle2 = tCircle2*(pz/energy)*lightSpeed + vz;
    if zPosCircle1 > length:
        tCircle1 = 1E9;
    if zPosCircle2 > length:
        tCircle2 = 1E9;
   
    #Now work out if it intersects the endcap
    tPlane1 = (length-vz)*energy/(pz*lightSpeed);
    tPlane2 = (-length-vz)*energy/(pz*lightSpeed);
    
    #If intersection in the past it doesn't count
    if tPlane1 < 0: tPlane1 = 1E9;
    if tPlane2 < 0: tPlane2 = 1E9;

    xPosPlane1 = tPlane1*(px/energy)*lightSpeed + vx;
    yPosPlane1 = tPlane1*(py/energy)*lightSpeed + vy;
    xPosPlane2 = tPlane2*(px/energy)*lightSpeed + vx;
    yPosPlane2 = tPlane2*(py/energy)*lightSpeed + vy;
   
    #If the intsersection occurs outside the endcap radius it doesn't count
    if math.sqrt(xPosPlane1*xPosPlane1+yPosPlane1*yPosPlane1) > radius:
        tPlane1 = 1E9;
    if math.sqrt(xPosPlane2*xPosPlane2+yPosPlane2*yPosPlane2) > radius:
        tPlane2 = 1E9;
  
    #Find the first intersection
    tInter = min([tCircle1,tCircle2,tPlane1,tPlane2])
    
    #Return 1000,1000 if not intersection with barrel or endcap
    if tInter > 1E6:
        return [1000,1000];

    #Find position of intersection
    xPos = tInter*(px/energy)*lightSpeed + vx;
    yPos = tInter*(py/energy)*lightSpeed + vy;
    zPos = tInter*(pz/energy)*lightSpeed + vz;

    #Find eta/phi of intersection
    phi = math.atan2(yPos,xPos);
    theta = math.acos(zPos/math.sqrt(xPos*xPos+yPos*yPos+zPos*zPos));
    eta = -math.log(math.tan(theta/2.));
    return (eta,phi);

print (intersec(0,100,280,0.5,0.5,1))