Static data member fgTop replaced with the static function GetTop();
[u/mrichter/AliRoot.git] / Flugg / WrapNorml.cxx
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
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)
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
66   G4bool valid;
67   normalLoc=ptrNavig->GetLocalExitNormal(&valid);
68   if(valid) {
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