Warning, /eic-opticks/u4/InstrumentedG4OpBoundaryProcess.rst is written in an unsupported language. File is not indexed.
0001 InstrumentedG4OpBoundaryProcess
0002 ===================================
0003
0004
0005 capture theGlobalExitNormal before any flips
0006 ------------------------------------------------
0007
0008 ::
0009
0010 434 G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
0011 435 std::vector<G4Navigator*>::iterator iNav =
0012 436 G4TransportationManager::GetTransportationManager()->
0013 437 GetActiveNavigatorsIterator();
0014 438 theGlobalNormal =
0015 439 (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
0016 440
0017 441 theGlobalExitNormal = theGlobalNormal ; // BEFORE ANY FLIPS
0018
0019
0020 Making sense of the normals
0021 ------------------------------
0022
0023
0024 CustomART::
0025
0026 1223 //dbg.nrm = oriented_normal ;
0027 1224 //dbg.nrm = surface_normal ; // verified that surface_normal always outwards
0028 1225 dbg.nrm = theGlobalExitNormal ; // inwards first, the rest outwards : oriented into direction of incident photon
0029 1226 //dbg.nrm = theGlobalNormal ; // this has been oriented : outwards first, the rest inwards
0030
0031
0032 Even theGlobalExitNormal has been flipped by the G4Navigator to point into the 2nd volume
0033
0034 This means::
0035
0036 theGlobalExitNormal*OldMomentum always +ve
0037
0038
0039 Is it possible to reliably recover the unflipped global normal somehow ?
0040 ----------------------------------------------------------------------------
0041
0042
0043
0044
0045
0046 theGlobalNormal and flips
0047 ----------------------------
0048
0049 ::
0050
0051 367 G4VParticleChange* InstrumentedG4OpBoundaryProcess::PostStepDoIt_(const G4Track& aTrack, const G4Step& aStep)
0052 368 {
0053 ...
0054 426 //[theGlobalNorml
0055 427 theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
0056 428
0057 429 G4bool valid;
0058 430 // Use the new method for Exit Normal in global coordinates,
0059 431 // which provides the normal more reliably.
0060 432 // ID of Navigator which limits step
0061 433
0062 434 G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
0063 435 std::vector<G4Navigator*>::iterator iNav =
0064 436 G4TransportationManager::GetTransportationManager()->
0065 437 GetActiveNavigatorsIterator();
0066 438 theGlobalNormal =
0067 439 (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
0068 440 #ifdef DEBUG_PIDX
0069 441 {
0070 442 quad2& prd = SEvt::Get()->current_prd ;
0071 443 prd.q0.f.x = theGlobalNormal.x() ;
0072 444 prd.q0.f.y = theGlobalNormal.y() ;
0073 445 prd.q0.f.z = theGlobalNormal.z() ;
0074 446
0075 447 // TRY USING PRE->POST POSITION CHANGE TO GET THE PRD DISTANCE ?
0076 448 G4ThreeVector theGlobalPoint_pre = pStep->GetPreStepPoint()->GetPosition();
0077 449 G4ThreeVector theGlobalPoint_delta = theGlobalPoint - theGlobalPoint_pre ;
0078 450 prd.q0.f.w = theGlobalPoint_delta.mag() ;
0079 451
0080 452 // HMM: PRD intersect identity ? how to mimic what Opticks does ?
0081 453 }
0082 454 #endif
0083 455
0084 456 if (valid)
0085 457 {
0086 458 theGlobalNormal = -theGlobalNormal;
0087 459 }
0088
0089
0090
0091 G4Navigator not providing the simple consistent Geometry Normal Needed for CustomBoundaryART
0092 ----------------------------------------------------------------------------------------------
0093
0094
0095 g4-cls G4Navigator::
0096
0097 227 virtual G4ThreeVector GetGlobalExitNormal(const G4ThreeVector& point,
0098 228 G4bool* valid);
0099 229 // Return Exit Surface Normal and validity too.
0100 230 // Can only be called if the Navigator's last Step has crossed a
0101 231 // volume geometrical boundary.
0102 232 // It returns the Normal to the surface pointing out of the volume that
0103 233 // was left behind and/or into the volume that was entered.
0104 234 // Convention:
0105 235 // The *local* normal is in the coordinate system of the *final* volume.
0106 236 // Restriction:
0107 237 // Normals are not available for replica volumes (returns valid= false)
0108 238 // These methods takes full care about how to calculate this normal,
0109 239 // but if the surfaces are not convex it will return valid=false.
0110 240
0111
0112
0113 1647 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,
0114 1648 &validNormal);
0115 1649 *pNormalCalculated = fCalculatedExitNormal;
0116 1650
0117 1651 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
0118 1652 globalNormal = localToGlobal.TransformAxis( localNormal );
0119 1653 }
0120
0121
0122 1329 // ********************************************************************
0123 1330 // GetLocalExitNormal
0124 1331 //
0125 1332 // Obtains the Normal vector to a surface (in local coordinates)
0126 1333 // pointing out of previous volume and into current volume
0127 1334 // ********************************************************************
0128 1335 //
0129 1336 G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
0130 1337 {
0131 1338 G4ThreeVector ExitNormal(0.,0.,0.);
0132 1339 G4VSolid *currentSolid=0;
0133 1340 G4LogicalVolume *candidateLogical;
0134 1341
0135 1342 if ( fLastTriedStepComputation )
0136 1343 {
0137 1344 // use fLastLocatedPointLocal and next candidate volume
0138 1345 //
0139 1346 G4ThreeVector nextSolidExitNormal(0.,0.,0.);
0140 1347
0141 1348 if( fEntering && (fBlockedPhysicalVolume!=0) )
0142 1349 {
0143 1350 candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume();
0144 1351 if( candidateLogical )
0145 1352 {
0146 1353 // fLastStepEndPointLocal is in the coordinates of the mother
0147 1354 // we need it in the daughter's coordinate system.
0148 1355
0149 1356 // The following code should also work in case of Replica
0150 1357 {
0151 1358 // First transform fLastLocatedPointLocal to the new daughter
0152 1359 // coordinates
0153 1360 //
0154 1361 G4AffineTransform MotherToDaughterTransform=
0155 1362 GetMotherToDaughterTransform( fBlockedPhysicalVolume,
0156 1363 fBlockedReplicaNo,
0157 1364 VolumeType(fBlockedPhysicalVolume) );
0158 1365 G4ThreeVector daughterPointOwnLocal=
0159 1366 MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
0160 1367
0161 1368 // OK if it is a parameterised volume
0162 1369 //
0163
0164
0165 1454 if ( EnteredDaughterVolume() )
0166 1455 {
0167 1456 G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume()
0168 1457 ->GetSolid();
0169 1458 ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
0170
0171
0172
0173 * daughterSolid->SurfaceNormal is outwards, but that gets flipped inwards
0174 when the incident photon is going from mother to daughter (Pyrex->Vacuum)
0175
0176 * so this is flipping based on volume hierarchy
0177
0178
0179
0180
0181
0182 Easiest to use G4VSolid::SurfaceNormal(localPoint) : to avoid all the flipping that hides the real normal
0183 -------------------------------------------------------------------------------------------------------------
0184
0185 But that normal is inherently local. Preferable to avoid having to transform into local
0186 and back.
0187
0188
0189
0190
0191 ::
0192
0193 127 virtual G4ThreeVector SurfaceNormal(const G4ThreeVector& p) const = 0;
0194 128 // Returns the outwards pointing unit normal of the shape for the
0195 129 // surface closest to the point at offset p.
0196
0197
0198 276 //////////////////////////////////////////////////////////////////////////
0199 277 //
0200 278 // Return unit normal of surface closest to p
0201 279
0202 280 G4ThreeVector G4Orb::SurfaceNormal( const G4ThreeVector& p ) const
0203 281 {
0204 282 return (1/p.mag())*p;
0205 283 }
0206
0207 Simply the gradient operator applied to the implicit eqn of the ellipsoid,
0208 gives the normal::
0209
0210 290 ///////////////////////////////////////////////////////////////////////////////
0211 291 //
0212 292 // Return unit normal of surface closest to p not protected against p=0
0213 293
0214 294 G4ThreeVector G4Ellipsoid::SurfaceNormal( const G4ThreeVector& p) const
0215 295 {
0216 296 G4double distR, distZBottom, distZTop;
0217 297
0218 298 // normal vector with special magnitude: parallel to normal, units 1/length
0219 299 // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside
0220 300 //
0221 301 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
0222 302 p.y()/(ySemiAxis*ySemiAxis),
0223 303 p.z()/(zSemiAxis*zSemiAxis));
0224 304 G4double radius = 1.0/norm.mag();
0225
0226 /// simple sphere case the axes are all 1 -> norm(x,y,z) radius=1
0227 /// so norm is the same as p, making distR zero for p on surface
0228
0229 305
0230 306 // approximate distance to curved surface
0231 307 //
0232 308 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
0233 309
0234 310 // Distance to z-cut plane
0235 311 //
0236 312 distZBottom = std::fabs( p.z() - zBottomCut );
0237 313 distZTop = std::fabs( p.z() - zTopCut );
0238
0239 /// extrema cuts are -r +r
0240
0241 314
0242 315 if ( (distZBottom < distR) || (distZTop < distR) )
0243 316 {
0244 317 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
0245 318 }
0246 319 return ( norm *= radius );
0247 320 }
0248
0249
