4 ///////////////////////////////////////////////////////////////////
6 // WrapLookMG.hh - Sara Vanini 26/X/99
8 // Wrapper for localisation of particle for magnetic field tracking
10 // modified 13/IV/00: check if point belongs to one of the lattice
11 // histories stored in jrLtGeant
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
18 ///////////////////////////////////////////////////////////////////
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"
32 //auxiliary function declarations
33 G4bool PointLocate(const G4NavigationHistory &,
34 const G4ThreeVector &,
35 const G4ThreeVector *,
37 EVolume CharacteriseDaughters(const G4LogicalVolume *);
40 void lkmgwr(G4double& pSx, G4double& pSy, G4double& pSz,
41 G4double* pV, const G4int& oldReg, const G4int& oldLttc,
42 G4int& flagErr, G4int& newReg, G4int& newLttc)
45 #ifdef G4GEOMETRY_DEBUG
46 G4cout<<"============= LKMGWR =============="<<G4endl;
49 //Geoinit, Navigator, etc. pointers
50 static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
51 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
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]);
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;
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();
74 //check if globalPoint belongs to volume (can't call
75 //LocateGlobalPoint... because of flag settings)
76 belongToVolume = PointLocate(*ptrNavHist,globalPoint,
77 &globalDirection,regionIn);
80 //if point belongs to surface, check scalar product direction*normal
82 #ifdef G4GEOMETRY_DEBUG
83 G4cout<<"On surface!"<<G4endl;
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];
94 px = globalDirection.x();
95 py = globalDirection.y();
96 pz = globalDirection.z();
98 nrmlwr(x,y,z,px,py,pz,norml,oldReg,reg,flagErr);
100 G4ThreeVector normal(norml[0],norml[1],norml[2]);
101 #ifdef G4GEOMETRY_DEBUG
102 // G4cout<<"Scalar product="<<globalDirection.dot(normal)<<G4endl;
104 if(globalDirection.dot(normal)>0) {
105 #ifdef G4GEOMETRY_DEBUG
106 G4cout<<"entering volume!"<<G4endl;
109 for (unsigned int i=0; i<pVolStore->size(); i++)
110 if ((*pVolStore)[i] == ptrNavHist->GetTopVolume())
112 if (volIndex==(~0)) {
113 G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
116 //regionIn = G4int(pVolStore->index(ptrNavHist->GetTopVolume()))+1;
117 regionIn = volIndex+1;
118 belongToVolume = true;
128 newLttc = jrLtGeant[i+1];
131 flagErr = pVolStore->size() + 2; //flagErr=fluka region number + 1
134 #ifdef G4GEOMETRY_DEBUG
135 G4cout << "Global point (cm): " << globalPointcm << G4endl;
136 G4cout << "Direction: " << globalDirection << G4endl;
138 G4cout << "Point belongs to region " << newReg
139 << " ptr history=" << newLttc << G4endl;
141 G4cout << "No containing volume found!" << G4endl;
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.
151 G4bool PointLocate(const G4NavigationHistory & blockedHistConst,
152 const G4ThreeVector & globalPoint,
153 const G4ThreeVector * pGlobalDirection,
158 G4NavigationHistory * fHistory = new G4NavigationHistory(blockedHistConst);
161 //G4 flags resetted (see: ResetStackAndState)
162 G4bool fExiting = false;
163 G4VPhysicalVolume * fBlockedPhysicalVolume=0;
164 G4int fBlockedReplicaNo=-1;
165 G4bool fLocatedOnEdge=false;
167 G4bool notKnownContained=true,noResult;
168 G4VPhysicalVolume *targetPhysical;
169 G4LogicalVolume *targetLogical;
170 G4VSolid *targetSolid;
171 G4ThreeVector localPoint,localDirection;
175 //Helpers/Utility classes
176 G4NormalNavigation fnormalNav;
177 G4VoxelNavigation fvoxelNav;
178 G4ParameterisedNavigation fparamNav;
179 G4ReplicaNavigation freplicaNav;
180 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
183 localPoint=fHistory->GetTopTransform().TransformPoint(globalPoint);
184 localDirection=fHistory->GetTopTransform().TransformPoint(*pGlobalDirection);
186 //search if volume contains point
187 if (fHistory->GetTopVolumeType()!=kReplica) {
188 targetSolid=fHistory->GetTopVolume()->GetLogicalVolume()->GetSolid();
189 insideCode=targetSolid->Inside(localPoint);
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;
204 if (insideCode==kOutside)
205 belongVolume = false;
207 else if (insideCode==kSurface) {
208 belongVolume = false;
212 // Search downwards in daughter volumes
213 // until deepest containing volume found
217 // o Parameterised daughters
218 // =>Must be one G4PVParameterised daughter & voxels
219 // o Positioned daughters & voxels
220 // o Positioned daughters & no voxels
226 // Determine `type' of current mother volume
227 targetPhysical=fHistory->GetTopVolume();
228 targetLogical=targetPhysical->GetLogicalVolume();
229 switch(CharacteriseDaughters(targetLogical)) {
231 if (targetLogical->GetVoxelHeader()) {
232 noResult=fvoxelNav.LevelLocate(*fHistory,
233 fBlockedPhysicalVolume,
241 noResult=fnormalNav.LevelLocate(*fHistory,
242 fBlockedPhysicalVolume,
251 noResult=freplicaNav.LevelLocate(*fHistory,
252 fBlockedPhysicalVolume,
260 noResult=fparamNav.LevelLocate(*fHistory,
261 fBlockedPhysicalVolume,
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.
276 // The blocked volume no longer valid - it was for another level
277 fBlockedPhysicalVolume= 0;
278 fBlockedReplicaNo= -1;
283 //on exit targetPhysical is the volume globalPoint belongs to;
285 for (unsigned int i=0; i<pVolStore->size(); i++)
286 if ((*pVolStore)[i] == targetPhysical)
288 if (volIndex==(~0)) {
289 G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
292 //reg = G4int(pVolStore->index(targetPhysical))+1;
301 EVolume CharacteriseDaughters(const G4LogicalVolume *pLog)
306 G4double width,offset;
308 G4VPhysicalVolume *pVol;
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;