#include #include #include #include #include using namespace RAT; SNOSVBellyPlate::SNOSVBellyPlate(const G4double rmin, const G4double rmax, const G4ThreeVector &xAxis, const G4ThreeVector &yAxis, const G4ThreeVector &zAxis, const G4double zInner, const G4double zOuter): fXaxis(xAxis), fYaxis(yAxis),fZaxis(zAxis){ //y axis goes through the center of the belly plate; x axis perpendicular and z vertical fThreshold = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); fRmin = rmin; fRmax = rmax; fZInner = zInner; fZOuter = zOuter; fRhomin = sqrt(fRmin*fRmin - fZInner*fZInner); fRhomax = sqrt(fRmax*fRmax - fZOuter*fZOuter); fRhoSlope = (fRhomax - fRhomin) / (fZOuter - fZInner); fZOffset = fZInner - fRhomin / fRhoSlope; } G4double SNOSVBellyPlate::Rho(G4double z)const { //returns cylindrical radius for a given z G4double rho; if (z<0)z = -z; if (zfZOuter)rho = sqrt(fRmax*fRmax - z*z); else rho = (z - fZOffset)*fRhoSlope; return rho; } G4double SNOSVBellyPlate::Radius(G4double z)const { //returns spherical radius for a given z G4double radius; double rho; if (z<0)z = -z; if (zfZOuter)radius = fRmax; else { rho = (z - fZOffset)*fRhoSlope; radius = sqrt(rho*rho + z*z); } return radius; } EInside SNOSVBellyPlate::Inside(const G4ThreeVector &p)const { EInside in; double x = p*fXaxis; double z = p*fZaxis; double r = Radius(x); double r2 = Radius(z); if (r2>r)r = r2; // always pick the maximum radius to deal with the curvature in both x and z directions. double delta = p.mag() - r; if (delta< -fThreshold)in = kInside; else if (delta>fThreshold)in = kOutside; else in = kSurface; return in; } G4ThreeVector SNOSVBellyPlate::SurfaceNormal(const G4ThreeVector &p, G4bool &isValid)const{ G4ThreeVector normal; isValid = true; // we set true , to allow us to return at any point. However, if we find that the point is not "within" the belly plate, we need to set it false. double x = p*fXaxis; double z = p*fZaxis; double r1 = Radius(x); double r2 = Radius(z); double r = r2>r1 ? r2 : r1; // always pick the maximum radius to deal with the curvature in both x and z directions. if (fabs(r - fRmin)r2){ //on one of the x-cones. For a cone oriented along the x-axis, with the point at the origin and we have rho=alpha*x=sqrt(y*y+z*z); // The normal to a point on the cone (x,y,z)=(x,rho*cos(phi),rho*sin(phi))=(x,alpha*x*cos(phi),alpha*x*sin(phi)) // is (x,-x/alpha*cos(phi),-x/alpha*sin(phi)); double phi = atan2(p.z(), p.y()); if (x>0){ double xx = p.x() - fZOffset; normal = G4ThreeVector(-xx, xx / fRhoSlope*cos(phi), xx / fRhoSlope*sin(phi)); } else{ double xx = p.x() + fZOffset; normal = G4ThreeVector(xx, -xx / fRhoSlope*cos(phi), -xx / fRhoSlope*sin(phi)); } } else{ //on one of the z-cones double phi = p.phi(); if (z>0){ double zz = p.z() - fZOffset; normal = G4ThreeVector(zz / fRhoSlope*cos(phi), zz / fRhoSlope*sin(phi), zz); } else{ double zz = p.z() + fZOffset; normal = G4ThreeVector(-zz / fRhoSlope*cos(phi), zz / fRhoSlope*sin(phi), -zz); } } return normal.unit(); } G4double SNOSVBellyPlate::DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &n, const G4double *, const G4double dOuter[2])const { // find the intersections with the surface of the belly plate. //Check the intersections with the inner radius // Note this routine ignores the proposed distance-- it finds the qualified intersection with all the surfaces of the bellyplate, and resets qualifiedDistance if smaller than before. double dist,dtest; G4ThreeVector test; int i; double x, z; // note that dInner for the SNOSV class refers to the id of the NCD's. Loosing irrelevant generality, // we make use of that fact by recalculating dinner for the Bellyplate. double dInner[2]; SphericalSubregionList::CalculateSphereIntersection(p,n,fRmin,dInner); if(dInner[0]>fThreshold)dtest=dInner[0]; else dtest=dInner[1]; dist=kInfinity; if (dtest>fThreshold){ test = p + dtest * n; x = fabs(test*fXaxis); z = fabs(test*fZaxis); double rr = test.mag(); if (fabs(rr - fRmin)fThreshold)dtest=dOuter[0];else dtest=dOuter[1]; if (dtest>fThreshold&&dtestfZOuter &&z >fZOuter)dist=dtest; } //calculate conical intersections : first Z-cones since this doesn't require any rotation double a, b, c, q, d; double a2 = fRhoSlope*fRhoSlope; z = p.z() - fZOffset; for (i = -1; i<2; i += 2){ z = p.z() + i*fZOffset; b = 2 * (p.x()*n.x() + p.y()*n.y() - a2*n.z()*z); c = p.x()*p.x() + p.y()*p.y() - a2*z*z; a = 1 - (1 + a2)*n.z()*n.z(); q = b*b - 4 * a*c; if (q >= 0){ q = sqrt(q); d = (-b + q) / (2 * a); if (d>fThreshold && dfThreshold && d= 0){ q = sqrt(q); d = (-b + q) / (2 * a); if (d>fThreshold && dfThreshold && d CLHEP::pi)phmax -= CLHEP::twopi; int i; bool retval=false; for(i=0;!retval && i<4;i++){ retval=(corners[i].cosTheta()>ctmin&&corners[i].cosTheta()phmin &&corners[i].phi()