Using TGeo to retrieve the mean material budget between two points (M.Ivanov)
[u/mrichter/AliRoot.git] / Flugg / WrapNorml.cxx
CommitLineData
26911512 1
2// Flugg tag
3
4////////////////////////////////////////////////////////////////////
5//
6// WrapNorml.hh - Sara Vanini
7//
8// Wrapper for computing normal unit-vector in global coordinates.
9//
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
22// end-point).
23//
24// modified 10/III/99
25// modified 25/V/00
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);
31//
32////////////////////////////////////////////////////////////////////
33
34#include "Wrappers.hh"
35#include "WrapUtils.hh"
36#include "FGeometryInit.hh"
37#include "G4VPhysicalVolume.hh"
38#include "FluggNavigator.hh"
39#include "G4ThreeVector.hh"
40#include "globals.hh"
41
42
c5a7ed1d 43void nrmlwr(G4double& /*pSx*/, G4double& /*pSy*/, G4double& /*pSz*/,
26911512 44 G4double& pVx, G4double& pVy, G4double& pVz,
45 G4double* norml, const G4int& oldReg,
c5a7ed1d 46 const G4int& /*newReg*/, G4int& flagErr)
26911512 47{
48 //flag
49#ifdef G4GEOMETRY_DEBUG
50 G4cout << "============ NRMLWR-DBG =============" << G4endl;
51#endif
52
53 //dummy variables
54 flagErr=0;
55
56 //navigator pointer
57 static FGeometryInit * ptrGeoInit;
58 ptrGeoInit = FGeometryInit::GetInstance();
59 FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
60
61 //variables
62 G4ThreeVector normalLoc;
63 G4ThreeVector normalGlob;
64
65 //normal computing
c5a7ed1d 66 G4bool valid;
67 normalLoc=ptrNavig->GetLocalExitNormal(&valid);
0edf14e7 68 if(valid) {
26911512 69 normalLoc *= -1;
70
71 //global cooordinates normal
72 normalGlob = ptrNavig->GetLocalToGlobalTransform().
73 TransformAxis(normalLoc);
74 }
75 else {
76 G4VPhysicalVolume *touchvolume;
77 G4ThreeVector theLocalPoint;
78
79 if(ptrNavig->EnteredDaughterVolume()) {
80 //volume from touchable history
81 G4TouchableHistory *ptrTouchableHistory=ptrGeoInit->
82 GetTouchableHistory();
83 touchvolume = ptrTouchableHistory->GetVolume();
84
85 //local point from navigator and normal
86 theLocalPoint = ptrNavig->GetCurrentLocalCoordinate();
87 normalLoc = touchvolume->GetLogicalVolume()->GetSolid()->
88 SurfaceNormal(theLocalPoint);
89
90 //global cooordinates normal
91 normalGlob = ptrNavig->GetLocalToGlobalTransform().
92 TransformAxis(normalLoc);
93 }
94 else {
95 //volume from old history
96 const G4NavigationHistory * ptrOldNavHist =
97 ptrGeoInit->GetOldNavHist()->GetHistory();
98 touchvolume = ptrOldNavHist->GetTopVolume();
99
100 if(!touchvolume) {
101 // Old history has been reseted by LOOKZ relocation,
102 // so is necessary to track back and forward to find
103 // the right histories.
104
105
106 //////////// COMPUTE STEP BACKWARD //////////////////
107 G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
108 TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
109
110 //compute step and new location
111 G4ThreeVector pVec(pVx,pVy,pVz);
112 pVec=-pVec;
113 G4double appStep = 0;
114 G4double safety = 0;
115 G4bool onBoundary = false;
116 G4double physStep = 1000000000;
117 G4int newRegStep;
118 G4ThreeVector partLoc = theGlobalPoint;
119 G4bool fErr=false;
120
121#ifdef G4GEOMETRY_DEBUG
122 G4cout << "Old history not found" << G4endl;
123 G4cout << "* NRML needs boundary-crossing: computing step backward..."
124 << G4endl;
125#endif
126
127 //compute step and location
128 newRegStep=StepAndLocation(partLoc,pVec,physStep,
129 appStep,safety,onBoundary,
130 fErr,oldReg);
131
132 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::
133 GetInstance();
134
135 if(appStep<physStep && newRegStep!=G4int(pVolStore->size())+1) {
136 //end-step point is on boundary between volumes;
137 //==> update step-histories
138
139#ifdef G4GEOMETRY_DEBUG
140 G4cout << "* updating step-histories" << G4endl;
141#endif
142
143 G4TouchableHistory * ptrTouchHist =
144 ptrGeoInit->GetTouchableHistory();
145 ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
146 }
147 else {
148#ifdef G4GEOMETRY_DEBUG
149 G4cout << "ERROR! Boundary not found" << G4endl;
150#endif
151 }
152
153#ifdef G4GEOMETRY_DEBUG
154 G4cout << "* computing step forward..." << G4endl;
155#endif
156 pVec=-pVec;
157 safety = 0;
158 onBoundary = false;
159
160 //compute step and location for boundary crossing
161 newRegStep=StepAndLocation(partLoc,pVec,physStep,
162 appStep,safety,onBoundary,fErr,oldReg);
163 if(appStep<physStep) {
164 //end-step point is on boundary between volumes;
165 //==> update step-histories
166
167#ifdef G4GEOMETRY_DEBUG
168 G4cout << "* updating step-histories" << G4endl;
169#endif
170
171 G4TouchableHistory * ptrTouchHist =
172 ptrGeoInit->GetTouchableHistory();
173 ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
174 }
175 }
176
177 // now touchvolume exist.
178 // local point from navigator and global point
179 // N.B. if particle has exited world volume,
180 // FluggNavigator doesn't update lastLocatedPoint.
181 // So be carefull in building geometry always to have a
182 // big world volume that fluka won't exit.
183
184 touchvolume = ptrOldNavHist->GetTopVolume();
185
186 G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
187 TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
188 theLocalPoint = ptrOldNavHist->GetTopTransform().
189 TransformPoint(theGlobalPoint);
190 normalLoc = (touchvolume->GetLogicalVolume()->GetSolid()->
191 SurfaceNormal(theLocalPoint));
192 normalLoc *= -1;
193
194 //global cooordinates normal
195 normalGlob = ptrOldNavHist->GetTopTransform().
196 Inverse().TransformAxis(normalLoc);
197 }
198 }
199
200 //return normal:
201 norml[0]=G4double(normalGlob.x());
202 norml[1]=G4double(normalGlob.y());
203 norml[2]=G4double(normalGlob.z());
204
205 //for debugging
206#ifdef G4GEOMETRY_DEBUG
207 G4cout << "Normal: " << normalGlob << G4endl;
208#endif
209}
210
211