inline std::tuple<TVector3,TVector3,double> find_PC_PathLength(const TVector3& x1, const TVector3& x2) {
    /*
        Function that finds the path length between anode and cathode surfaces, both one-sheet hyperboloids of form (x*x+y*y)/(a*a) - (z*z)/(c*c) = 1
        * path length found for a given particle moving along a certain direction from x1 to x2
        * Typical arguments here will be x1=r_rhoMin, x2=qqqevent.pos or sx3event.pos
    */
    TVector3 dx = x2-x1;// direction vector
    double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz 
    auto onesheet_hyperboloid_intersect = [&](double a, double c) {
		auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c);
		auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c);
		auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0;
		double disc = B*B - 4*A*C;
		if(disc<0) 
			return TVector3(0,0,54321);
		else {
			double tsol1 = (-B + TMath::Sqrt(disc))/(2*A);
			double tsol2 = (-B - TMath::Sqrt(disc))/(2*A);
			if(tsol1 >= 0  && tsol1 <= t2) 
				return x1+tsol1*dx;
			else if(tsol2>=0 && tsol2 <= t2)
				return x1+tsol2*dx;
			else
				return TVector3(0,0,54321);
		}
	};

	//TODO: Magic numbers here describing waist 'a', and flare 'c' will need better treatment. 
	//Currently, these are derived by fitting the crossover points to R^2/a^2 - z^2/c^2 = 1 for anodes
	// Cathode a, c values are found by scaling up the anode waist by 43/37, the ratio of the outermost radii
	TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895);
	TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895);
	TVector3 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895);
	if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321)
		return std::tuple(cathode_intersect,gw_intersect,(cathode_intersect-gw_intersect).Mag()*0.1);
	else
		return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321);
}