#include #include #include #include using namespace RAT; SNONCDAnchor::SNONCDAnchor(G4double rmin, G4double rmax, G4double radius, const G4ThreeVector &position):fPosition(position){ fDirection = position.unit(); fRmin = rmin; fRmax = rmax; fRadius = radius; fThreshold = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); } EInside SNONCDAnchor::Inside(const G4ThreeVector &p)const { // JP - Adding a surface possibility (NCD top) and handling the non-NCD surface EInside in; double radius = Radius(p); double z=p.mag(); z=p.mag(); //JP addition if(zfRmax+fThreshold) return kOutside; if (radius0){ if(z>fRmin + fThreshold) in=kOutside; else in=kSurface; } else if (radius>fRadius + fThreshold){ // JP mods - if not really on the NCD, still check for the SV std volume if(zfRmin - fThreshold && z < fRmax + fThreshold){ isValid = true; if (zunit(); } G4double SNONCDAnchor::DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &n, const G4double dInner[2], const G4double dOuter[2])const{ double dist=kInfinity; G4ThreeVector test = n.cross(fDirection); test = test.unit(); if (fabs(test*(p - fPosition))>fRadius)return dist; // no intersection possible double dtest, z; // JPmod instead of checking only one dInner, the check should loop over both // when both are positive the particle could be traveling from an NCD region to another one, missing the first NCD // and checking only one means missing the second intersection int nref; for(nref=0; nref<2; nref++){ dtest = dInner[nref]; if(dtest>fThreshold && dtest0?test.mag():-test.mag(); // JP modified the skew calculation to include the fact that test is at the top of the NCD and fposition at bottom G4ThreeVector skew = (test - (fRmin*fDirection)); double radius=skew.mag(); // JP changed this condition. It was requiring being on the surface of the cylinder instead of the top if(radiusfThreshold && d1fRmin&&zfThreshold && d2fRmin&&z 0){ //if dr<0, just return dist double dz = z-fRmin; if (dist<0)dist = dz < dr ? dz : dr; else dist = sqrt(dr*dr + dist*dist); } if(dist<0) return 0; else return dist; } G4double SNONCDAnchor::DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &n, const G4double dInner[2], const G4double dOuter[2])const{ double dist= DistanceToIn(p, n, dInner, dOuter); return dist; } G4double SNONCDAnchor::DistanceToOut(const G4ThreeVector &p)const { double z = p*fDirection; G4ThreeVector skew = (p - fPosition); skew = skew - (skew*fDirection)*fDirection; //vector perpendicular to cylinder axis, through p. double radius = skew.mag(); double dr = radius - fRadius; double dist=fRmin - z;; if (dr > 0){ //if dr<0, just return dist double dz = fRmax - z; if (dist<0)dist = dz < dr ? dz : dr; else { double d = sqrt(dr*dr + dz*dz); dist = dz < d ? dz : d; } } return dist; } G4bool SNONCDAnchor::OverlapsRegion(const G4ThreeVector corners[4])const { int i; double minang=-2; for(i=0;i<4;i++){ minang=std::max(minang,corners[i].unit()*fDirection); } if(minang>cos((acos(fRmin/fRmax)+fRadius/fRmax)*1.05)){ return true; } return false; }