First implementation. Needs cleanup.
[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   if(ptrNavig->IsExitNormalValid()) {
67     normalLoc=ptrNavig->GetLocalExitNormal();
68     normalLoc *= -1;        
69     
70     //global cooordinates normal
71     normalGlob = ptrNavig->GetLocalToGlobalTransform().
72       TransformAxis(normalLoc);
73   }
74   else {
75     G4VPhysicalVolume *touchvolume;
76     G4ThreeVector theLocalPoint;
77     
78     if(ptrNavig->EnteredDaughterVolume()) {
79       //volume from touchable history
80       G4TouchableHistory *ptrTouchableHistory=ptrGeoInit->
81         GetTouchableHistory();
82       touchvolume = ptrTouchableHistory->GetVolume();
83       
84       //local point from navigator and normal
85       theLocalPoint = ptrNavig->GetCurrentLocalCoordinate(); 
86       normalLoc = touchvolume->GetLogicalVolume()->GetSolid()->
87         SurfaceNormal(theLocalPoint);
88       
89       //global cooordinates normal
90       normalGlob = ptrNavig->GetLocalToGlobalTransform().
91         TransformAxis(normalLoc);
92     }
93     else {
94       //volume from old history
95       const G4NavigationHistory * ptrOldNavHist = 
96         ptrGeoInit->GetOldNavHist()->GetHistory();
97       touchvolume = ptrOldNavHist->GetTopVolume();
98       
99       if(!touchvolume) {
100         // Old history has been reseted by LOOKZ relocation,
101         // so is necessary to track back and forward to find
102         // the right histories.
103         
104         
105         ////////////  COMPUTE STEP BACKWARD //////////////////
106         G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
107           TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
108         
109         //compute step and new location
110         G4ThreeVector pVec(pVx,pVy,pVz);
111         pVec=-pVec;
112         G4double appStep = 0;
113         G4double safety = 0;
114         G4bool onBoundary = false;
115         G4double physStep = 1000000000;
116         G4int newRegStep;
117         G4ThreeVector partLoc =  theGlobalPoint;
118         G4bool fErr=false;
119         
120 #ifdef G4GEOMETRY_DEBUG
121         G4cout << "Old history not found" << G4endl;
122         G4cout << "* NRML needs boundary-crossing: computing step backward..."
123                << G4endl;
124 #endif
125         
126         //compute step and location 
127         newRegStep=StepAndLocation(partLoc,pVec,physStep,
128                                    appStep,safety,onBoundary,
129                                    fErr,oldReg);
130         
131         G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::
132           GetInstance();
133         
134         if(appStep<physStep && newRegStep!=G4int(pVolStore->size())+1) {
135           //end-step point is on boundary between volumes;
136           //==> update step-histories
137           
138 #ifdef G4GEOMETRY_DEBUG
139           G4cout << "* updating step-histories" << G4endl;
140 #endif
141           
142           G4TouchableHistory * ptrTouchHist = 
143             ptrGeoInit->GetTouchableHistory();
144           ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
145         }
146         else {
147 #ifdef G4GEOMETRY_DEBUG
148           G4cout << "ERROR! Boundary not found" << G4endl;
149 #endif
150         }
151         
152 #ifdef G4GEOMETRY_DEBUG
153         G4cout << "* computing step forward..." << G4endl;
154 #endif
155         pVec=-pVec;
156         safety = 0;
157         onBoundary = false;
158         
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
165           
166 #ifdef G4GEOMETRY_DEBUG
167           G4cout << "* updating step-histories" << G4endl;
168 #endif
169           
170           G4TouchableHistory * ptrTouchHist = 
171             ptrGeoInit->GetTouchableHistory();
172           ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
173         }
174       }
175       
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.
182       
183       touchvolume = ptrOldNavHist->GetTopVolume();
184       
185       G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
186         TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
187       theLocalPoint = ptrOldNavHist->GetTopTransform().
188         TransformPoint(theGlobalPoint);
189       normalLoc = (touchvolume->GetLogicalVolume()->GetSolid()->
190                    SurfaceNormal(theLocalPoint));
191       normalLoc *= -1; 
192       
193       //global cooordinates normal
194       normalGlob = ptrOldNavHist->GetTopTransform().
195         Inverse().TransformAxis(normalLoc);         
196     }
197   }
198   
199   //return normal:
200   norml[0]=G4double(normalGlob.x());
201   norml[1]=G4double(normalGlob.y());
202   norml[2]=G4double(normalGlob.z());
203   
204   //for debugging
205 #ifdef G4GEOMETRY_DEBUG
206   G4cout << "Normal: " << normalGlob << G4endl;
207 #endif
208 }
209
210