4 ////////////////////////////////////////////////////////////////////
6 // WrapNorml.hh - Sara Vanini
8 // Wrapper for computing normal unit-vector in global coordinates.
10 // Fluka requires normal vector exiting from final position (of
11 // particle) volume, that is: entering in volume of initial position.
12 // Geant4 always computes the normal vector exiting from the volume.
13 // In GetLocalExitNormal() call the volume is the pre-step volume
14 // (so G4 normal vector sign is opposite of fluka-required normal).
15 // If IsLocalExitNormalValid=false, normal value is computed from
16 // init-step volume (in this case sign must change), or from
17 // end-step volume (the sign is the same). Normal vector is computed
18 // always on boundary between volumes, in global coordinates (to take
19 // rotation of parameterised volumes in hierarchy in consideration).
20 // So: nrmlwr returns inwards pointing unit normal of the shape for
21 // surface closest to the point returned by the navigator (last step
26 // modified 7/VI/00 for boundary-crossing in case of relocation
27 // modified 5/VII/00 geometry error on boundary fixed
28 // modified 24.10.01: by I. Hrivnacova
29 // functions declarations separated from implementation
30 // (moved to Wrappers.hh);
32 ////////////////////////////////////////////////////////////////////
34 #include "Wrappers.hh"
35 #include "WrapUtils.hh"
36 #include "FGeometryInit.hh"
37 #include "G4VPhysicalVolume.hh"
38 #include "FluggNavigator.hh"
39 #include "G4ThreeVector.hh"
43 void nrmlwr(G4double& pSx, G4double& pSy, G4double& pSz,
44 G4double& pVx, G4double& pVy, G4double& pVz,
45 G4double* norml, const G4int& oldReg,
46 const G4int& newReg, G4int& flagErr)
49 #ifdef G4GEOMETRY_DEBUG
50 G4cout << "============ NRMLWR-DBG =============" << G4endl;
57 static FGeometryInit * ptrGeoInit;
58 ptrGeoInit = FGeometryInit::GetInstance();
59 FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
62 G4ThreeVector normalLoc;
63 G4ThreeVector normalGlob;
66 if(ptrNavig->IsExitNormalValid()) {
67 normalLoc=ptrNavig->GetLocalExitNormal();
70 //global cooordinates normal
71 normalGlob = ptrNavig->GetLocalToGlobalTransform().
72 TransformAxis(normalLoc);
75 G4VPhysicalVolume *touchvolume;
76 G4ThreeVector theLocalPoint;
78 if(ptrNavig->EnteredDaughterVolume()) {
79 //volume from touchable history
80 G4TouchableHistory *ptrTouchableHistory=ptrGeoInit->
81 GetTouchableHistory();
82 touchvolume = ptrTouchableHistory->GetVolume();
84 //local point from navigator and normal
85 theLocalPoint = ptrNavig->GetCurrentLocalCoordinate();
86 normalLoc = touchvolume->GetLogicalVolume()->GetSolid()->
87 SurfaceNormal(theLocalPoint);
89 //global cooordinates normal
90 normalGlob = ptrNavig->GetLocalToGlobalTransform().
91 TransformAxis(normalLoc);
94 //volume from old history
95 const G4NavigationHistory * ptrOldNavHist =
96 ptrGeoInit->GetOldNavHist()->GetHistory();
97 touchvolume = ptrOldNavHist->GetTopVolume();
100 // Old history has been reseted by LOOKZ relocation,
101 // so is necessary to track back and forward to find
102 // the right histories.
105 //////////// COMPUTE STEP BACKWARD //////////////////
106 G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
107 TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
109 //compute step and new location
110 G4ThreeVector pVec(pVx,pVy,pVz);
112 G4double appStep = 0;
114 G4bool onBoundary = false;
115 G4double physStep = 1000000000;
117 G4ThreeVector partLoc = theGlobalPoint;
120 #ifdef G4GEOMETRY_DEBUG
121 G4cout << "Old history not found" << G4endl;
122 G4cout << "* NRML needs boundary-crossing: computing step backward..."
126 //compute step and location
127 newRegStep=StepAndLocation(partLoc,pVec,physStep,
128 appStep,safety,onBoundary,
131 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::
134 if(appStep<physStep && newRegStep!=G4int(pVolStore->size())+1) {
135 //end-step point is on boundary between volumes;
136 //==> update step-histories
138 #ifdef G4GEOMETRY_DEBUG
139 G4cout << "* updating step-histories" << G4endl;
142 G4TouchableHistory * ptrTouchHist =
143 ptrGeoInit->GetTouchableHistory();
144 ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
147 #ifdef G4GEOMETRY_DEBUG
148 G4cout << "ERROR! Boundary not found" << G4endl;
152 #ifdef G4GEOMETRY_DEBUG
153 G4cout << "* computing step forward..." << G4endl;
159 //compute step and location for boundary crossing
160 newRegStep=StepAndLocation(partLoc,pVec,physStep,
161 appStep,safety,onBoundary,fErr,oldReg);
162 if(appStep<physStep) {
163 //end-step point is on boundary between volumes;
164 //==> update step-histories
166 #ifdef G4GEOMETRY_DEBUG
167 G4cout << "* updating step-histories" << G4endl;
170 G4TouchableHistory * ptrTouchHist =
171 ptrGeoInit->GetTouchableHistory();
172 ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
176 // now touchvolume exist.
177 // local point from navigator and global point
178 // N.B. if particle has exited world volume,
179 // FluggNavigator doesn't update lastLocatedPoint.
180 // So be carefull in building geometry always to have a
181 // big world volume that fluka won't exit.
183 touchvolume = ptrOldNavHist->GetTopVolume();
185 G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
186 TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
187 theLocalPoint = ptrOldNavHist->GetTopTransform().
188 TransformPoint(theGlobalPoint);
189 normalLoc = (touchvolume->GetLogicalVolume()->GetSolid()->
190 SurfaceNormal(theLocalPoint));
193 //global cooordinates normal
194 normalGlob = ptrOldNavHist->GetTopTransform().
195 Inverse().TransformAxis(normalLoc);
200 norml[0]=G4double(normalGlob.x());
201 norml[1]=G4double(normalGlob.y());
202 norml[2]=G4double(normalGlob.z());
205 #ifdef G4GEOMETRY_DEBUG
206 G4cout << "Normal: " << normalGlob << G4endl;