Including St2 in the new geometry segmentation (Christian)
[u/mrichter/AliRoot.git] / Flugg / WrapLookMG.cxx
CommitLineData
26911512 1
2// Flugg tag
3
4///////////////////////////////////////////////////////////////////
5//
6// WrapLookMG.hh - Sara Vanini 26/X/99
7//
8// Wrapper for localisation of particle for magnetic field tracking
9//
10// modified 13/IV/00: check if point belongs to one of the lattice
11// histories stored in jrLtGeant
12//
13// modified 24.10.01: by I. Hrivnacova
14// functions declarations separated from implementation
15// (moved to Wrappers.hh);
16// modified 17/06/02: by I. Gonzalez. STL migration
17//
18///////////////////////////////////////////////////////////////////
19
20
21#include "Wrappers.hh"
22#include "FGeometryInit.hh"
23#include "NavHistWithCount.hh"
24#include "G4PhysicalVolumeStore.hh"
25#include "G4NormalNavigation.hh"
26#include "G4VoxelNavigation.hh"
27#include "G4ParameterisedNavigation.hh"
28#include "G4ReplicaNavigation.hh"
29#include "globals.hh"
30
31
32//auxiliary function declarations
33G4bool PointLocate(const G4NavigationHistory &,
34 const G4ThreeVector &,
35 const G4ThreeVector *,
36 G4int &);
37EVolume CharacteriseDaughters(const G4LogicalVolume *);
38
39
40void lkmgwr(G4double& pSx, G4double& pSy, G4double& pSz,
41 G4double* pV, const G4int& oldReg, const G4int& oldLttc,
42 G4int& flagErr, G4int& newReg, G4int& newLttc)
43{
44//flag
45#ifdef G4GEOMETRY_DEBUG
46 G4cout<<"============= LKMGWR =============="<<G4endl;
47#endif
48
49 //Geoinit, Navigator, etc. pointers
50 static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
51 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
52
53 //setting variables (and dimension: Fluka uses cm.!)
54 G4ThreeVector globalPointcm(pSx,pSy,pSz);
55 G4ThreeVector globalPoint = globalPointcm * 10.; //in mm.
56 G4ThreeVector globalDirection(pV[0],pV[1],pV[2]);
57
58 //get jrLtGeant array and initialize variables
59 G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
60 G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
61 G4bool belongToVolume = false;
62 G4int i = LttcFlagGeant;
63 G4int regionIn = 0;
64 newReg = -1;
65 newLttc = -1;
66
67 while(!belongToVolume && i>=0) {
68 //get history from jrLtGeant
69 NavHistWithCount * ptrNavHistCount =
70 reinterpret_cast<NavHistWithCount*>(jrLtGeant[i]);
71 const G4NavigationHistory * ptrNavHist =
72 ptrNavHistCount->GetNavHistPtr();
73
74 //check if globalPoint belongs to volume (can't call
75 //LocateGlobalPoint... because of flag settings)
76 belongToVolume = PointLocate(*ptrNavHist,globalPoint,
77 &globalDirection,regionIn);
78
79
80 //if point belongs to surface, check scalar product direction*normal
81 if(regionIn==-100) {
82#ifdef G4GEOMETRY_DEBUG
83 G4cout<<"On surface!"<<G4endl;
84#endif
85
86 //if entering, then point belongs to volume,
87 //if exiting, point doesn't belong to volume.
88 G4int oldReg,flagErr,reg;
89 G4double x,y,z,px,py,pz;
90 G4double * norml = new G4double[3];
91 x = globalPoint.x();
92 y = globalPoint.y();
93 z = globalPoint.z();
94 px = globalDirection.x();
95 py = globalDirection.y();
96 pz = globalDirection.z();
97
98 nrmlwr(x,y,z,px,py,pz,norml,oldReg,reg,flagErr);
99
100 G4ThreeVector normal(norml[0],norml[1],norml[2]);
101#ifdef G4GEOMETRY_DEBUG
102 // G4cout<<"Scalar product="<<globalDirection.dot(normal)<<G4endl;
103#endif
104 if(globalDirection.dot(normal)>0) {
105#ifdef G4GEOMETRY_DEBUG
106 G4cout<<"entering volume!"<<G4endl;
107#endif
108 G4int volIndex = ~0;
109 for (unsigned int i=0; i<pVolStore->size(); i++)
110 if ((*pVolStore)[i] == ptrNavHist->GetTopVolume())
111 volIndex = i;
112 if (volIndex==(~0)) {
113 G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
114 exit(-999);
115 }
116 //regionIn = G4int(pVolStore->index(ptrNavHist->GetTopVolume()))+1;
117 regionIn = volIndex+1;
118 belongToVolume = true;
119 }
120 }
121
122 i -= 1;
123 }
124
125 //output variables
126 if(belongToVolume) {
127 newReg = regionIn;
128 newLttc = jrLtGeant[i+1];
129 }
130
131 flagErr = pVolStore->size() + 2; //flagErr=fluka region number + 1
132
133
134#ifdef G4GEOMETRY_DEBUG
135 G4cout << "Global point (cm): " << globalPointcm << G4endl;
136 G4cout << "Direction: " << globalDirection << G4endl;
137 if(newReg!=-1)
138 G4cout << "Point belongs to region " << newReg
139 << " ptr history=" << newLttc << G4endl;
140 else
141 G4cout << "No containing volume found!" << G4endl;
142#endif
143}
144
145
146
147
148//returns false if the point doesn't belong to history top volume (can either
149//belong to one of its daughter volumes or none), otherwise returns true.
150
151G4bool PointLocate(const G4NavigationHistory & blockedHistConst,
152 const G4ThreeVector & globalPoint,
153 const G4ThreeVector * pGlobalDirection,
154 G4int & reg)
155{
156
157 //variables
158 G4NavigationHistory * fHistory = new G4NavigationHistory(blockedHistConst);
159 G4bool belongVolume;
160
161 //G4 flags resetted (see: ResetStackAndState)
162 G4bool fExiting = false;
163 G4VPhysicalVolume * fBlockedPhysicalVolume=0;
164 G4int fBlockedReplicaNo=-1;
165 G4bool fLocatedOnEdge=false;
166
167 G4bool notKnownContained=true,noResult;
168 G4VPhysicalVolume *targetPhysical;
169 G4LogicalVolume *targetLogical;
170 G4VSolid *targetSolid;
171 G4ThreeVector localPoint,localDirection;
172 EInside insideCode;
173
174
175 //Helpers/Utility classes
176 G4NormalNavigation fnormalNav;
177 G4VoxelNavigation fvoxelNav;
178 G4ParameterisedNavigation fparamNav;
179 G4ReplicaNavigation freplicaNav;
180 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
181
182 //local variables
183 localPoint=fHistory->GetTopTransform().TransformPoint(globalPoint);
184 localDirection=fHistory->GetTopTransform().TransformPoint(*pGlobalDirection);
185
186 //search if volume contains point
187 if (fHistory->GetTopVolumeType()!=kReplica) {
188 targetSolid=fHistory->GetTopVolume()->GetLogicalVolume()->GetSolid();
189 insideCode=targetSolid->Inside(localPoint);
190 }
191 else {
192 insideCode=freplicaNav.BackLocate(*fHistory,globalPoint,
193 localPoint,fExiting,notKnownContained);
194 // !CARE! if notKnownContained returns false then the point is within
195 // the containing placement volume of the replica(s). If insidecode
196 // will result in the history being backed up one level, then the
197 // local point returned is the point in the system of this new level
198#ifdef G4GEOMETRY_DEBUG
199 G4cout << "Replica: fExiting=" << fExiting << G4endl;
200 G4cout << " notKnownContained=" << notKnownContained << G4endl;
201#endif
202 }
203
204 if (insideCode==kOutside)
205 belongVolume = false;
206
207 else if (insideCode==kSurface) {
208 belongVolume = false;
209 reg = -100;
210 }
211 else {
212 // Search downwards in daughter volumes
213 // until deepest containing volume found
214 //
215 // 3 Cases:
216 //
217 // o Parameterised daughters
218 // =>Must be one G4PVParameterised daughter & voxels
219 // o Positioned daughters & voxels
220 // o Positioned daughters & no voxels
221
222 belongVolume = true;
223 noResult=true;
224
225 do {
226 // Determine `type' of current mother volume
227 targetPhysical=fHistory->GetTopVolume();
228 targetLogical=targetPhysical->GetLogicalVolume();
229 switch(CharacteriseDaughters(targetLogical)) {
230 case kNormal:
231 if (targetLogical->GetVoxelHeader()) {
232 noResult=fvoxelNav.LevelLocate(*fHistory,
233 fBlockedPhysicalVolume,
234 fBlockedReplicaNo,
235 globalPoint,
236 pGlobalDirection,
237 fLocatedOnEdge,
238 localPoint);
239 }
240 else {
241 noResult=fnormalNav.LevelLocate(*fHistory,
242 fBlockedPhysicalVolume,
243 fBlockedReplicaNo,
244 globalPoint,
245 pGlobalDirection,
246 fLocatedOnEdge,
247 localPoint);
248 }
249 break;
250 case kReplica:
251 noResult=freplicaNav.LevelLocate(*fHistory,
252 fBlockedPhysicalVolume,
253 fBlockedReplicaNo,
254 globalPoint,
255 pGlobalDirection,
256 fLocatedOnEdge,
257 localPoint);
258 break;
259 case kParameterised:
260 noResult=fparamNav.LevelLocate(*fHistory,
261 fBlockedPhysicalVolume,
262 fBlockedReplicaNo,
263 globalPoint,
264 pGlobalDirection,
265 fLocatedOnEdge,
266 localPoint);
267 break;
268 } //switch
269
270 // LevelLocate search in the first daughter level.
271 // LevelLocate returns noResult=true if it finds a daughter volume
272 // in which globalPoint is inside (or on the surface). So point
273 // doesn't belong only to mother volume ==> belongVolume=false.
274
275 if (noResult) {
276 // The blocked volume no longer valid - it was for another level
277 fBlockedPhysicalVolume= 0;
278 fBlockedReplicaNo= -1;
279 belongVolume=false;
280 }
281 } while (noResult);
282
283 //on exit targetPhysical is the volume globalPoint belongs to;
284 G4int volIndex = ~0;
285 for (unsigned int i=0; i<pVolStore->size(); i++)
286 if ((*pVolStore)[i] == targetPhysical)
287 volIndex = i;
288 if (volIndex==(~0)) {
289 G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
290 exit(-999);
291 }
292 //reg = G4int(pVolStore->index(targetPhysical))+1;
293 reg = volIndex+1;
294
295 }
296
297 delete fHistory;
298 return belongVolume;
299}
300
301EVolume CharacteriseDaughters(const G4LogicalVolume *pLog)
302{
303 EVolume type;
304 EAxis axis;
305 G4int nReplicas;
306 G4double width,offset;
307 G4bool consuming;
308 G4VPhysicalVolume *pVol;
309
310 if (pLog->GetNoDaughters()==1) {
311 pVol=pLog->GetDaughter(0);
312 if (pVol->IsReplicated()) {
313 pVol->GetReplicationData(axis,nReplicas,width,offset,consuming);
314 type=(consuming) ? kReplica : kParameterised;
315 }
316 else
317 type=kNormal;
318 }
319 else
320 type=kNormal;
321 return type;
322}
323
324
325