Implementation of ITS tracking for HLT. The tracking is based on the off-line AliITSt...
[u/mrichter/AliRoot.git] / Flugg / WrapLookMG.cxx
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
33 G4bool PointLocate(const G4NavigationHistory &,
34                           const G4ThreeVector &, 
35                           const G4ThreeVector *,
36                           G4int &);
37 EVolume CharacteriseDaughters(const G4LogicalVolume *);
38
39
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)
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
151 G4bool 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
301 EVolume 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