[u/mrichter/AliRoot.git] / Flugg / WrapLookMG.cxx


// Flugg tag 

///////////////////////////////////////////////////////////////////
//
// WrapLookMG.hh - Sara Vanini 26/X/99
//
// Wrapper for localisation of particle for magnetic field tracking 
//
// modified 13/IV/00: check if point belongs to one of the lattice 
// histories stored in jrLtGeant 
//
// modified 24.10.01: by I. Hrivnacova
//   functions declarations separated from implementation
//   (moved to Wrappers.hh);
// modified 17/06/02: by I. Gonzalez. STL migration
//
///////////////////////////////////////////////////////////////////


#include "Wrappers.hh"
#include "FGeometryInit.hh"
#include "NavHistWithCount.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4NormalNavigation.hh"
#include "G4VoxelNavigation.hh"
#include "G4ParameterisedNavigation.hh"
#include "G4ReplicaNavigation.hh"
#include "globals.hh"


//auxiliary function declarations
G4bool PointLocate(const G4NavigationHistory &,
			  const G4ThreeVector &, 
			  const G4ThreeVector *,
			  G4int &);
EVolume CharacteriseDaughters(const G4LogicalVolume *);


void lkmgwr(G4double& pSx, G4double& pSy, G4double& pSz,
            G4double* pV, const G4int& oldReg, const G4int& oldLttc,
	    G4int& flagErr, G4int& newReg, G4int& newLttc)
{
//flag
#ifdef G4GEOMETRY_DEBUG
    G4cout<<"============= LKMGWR =============="<<G4endl;
#endif

    //Geoinit, Navigator, etc. pointers
    static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
    G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
    
    //setting variables (and dimension: Fluka uses cm.!)
    G4ThreeVector globalPointcm(pSx,pSy,pSz);
    G4ThreeVector globalPoint =  globalPointcm * 10.; //in mm.
    G4ThreeVector globalDirection(pV[0],pV[1],pV[2]);
    
    //get jrLtGeant array and initialize variables
    G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
    G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
    G4bool belongToVolume = false;
    G4int i = LttcFlagGeant;
    G4int regionIn = 0;
    newReg = -1;
    newLttc = -1;
    
    while(!belongToVolume && i>=0) {
      //get history from jrLtGeant
      NavHistWithCount * ptrNavHistCount = 
	reinterpret_cast<NavHistWithCount*>(jrLtGeant[i]);
      const G4NavigationHistory * ptrNavHist =
	ptrNavHistCount->GetNavHistPtr();
      
      //check if globalPoint belongs to volume (can't call 
      //LocateGlobalPoint... because of flag settings)
      belongToVolume = PointLocate(*ptrNavHist,globalPoint,
				   &globalDirection,regionIn);
      
      
      //if point belongs to surface, check scalar product direction*normal
      if(regionIn==-100) {
#ifdef G4GEOMETRY_DEBUG
	G4cout<<"On surface!"<<G4endl;
#endif
	
	//if entering, then point belongs to volume, 
	//if exiting, point doesn't belong to volume.
	G4int oldReg,flagErr,reg;
	G4double x,y,z,px,py,pz;
	G4double * norml = new G4double[3];
	x = globalPoint.x();
	y = globalPoint.y();
	z = globalPoint.z();
	px = globalDirection.x();
	py = globalDirection.y();
	pz = globalDirection.z();
	
	nrmlwr(x,y,z,px,py,pz,norml,oldReg,reg,flagErr);
	
	G4ThreeVector normal(norml[0],norml[1],norml[2]);
#ifdef G4GEOMETRY_DEBUG
	// G4cout<<"Scalar product="<<globalDirection.dot(normal)<<G4endl;
#endif
	if(globalDirection.dot(normal)>0) {
#ifdef G4GEOMETRY_DEBUG
	  G4cout<<"entering volume!"<<G4endl;
#endif
	  G4int volIndex = ~0;
	  for (unsigned int i=0; i<pVolStore->size(); i++)
	    if ((*pVolStore)[i] == ptrNavHist->GetTopVolume()) 
	      volIndex = i;
	  if (volIndex==(~0)) {
	    G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
	    exit(-999);
	  }
	  //regionIn = G4int(pVolStore->index(ptrNavHist->GetTopVolume()))+1;
	  regionIn = volIndex+1; 
	  belongToVolume = true;
	}
      }
      
      i -= 1;
    }
    
    //output variables
    if(belongToVolume) {
      newReg = regionIn;
      newLttc = jrLtGeant[i+1];
    }

    flagErr = pVolStore->size() + 2;      //flagErr=fluka region number + 1
    
    
#ifdef G4GEOMETRY_DEBUG
    G4cout << "Global point (cm): " << globalPointcm << G4endl;
    G4cout << "Direction: " << globalDirection << G4endl;
    if(newReg!=-1) 
      G4cout << "Point belongs to region " << newReg
	     << " ptr history=" << newLttc << G4endl;
    else 
      G4cout << "No containing volume found!" << G4endl;
#endif
}


//returns false if the point doesn't belong to history top volume (can either 
//belong to one of its daughter volumes or none), otherwise returns true.

G4bool PointLocate(const G4NavigationHistory & blockedHistConst,
		   const G4ThreeVector & globalPoint, 
		   const G4ThreeVector * pGlobalDirection,
		   G4int & reg)
{
  
  //variables
  G4NavigationHistory * fHistory = new G4NavigationHistory(blockedHistConst);
  G4bool belongVolume;
  
  //G4 flags resetted (see: ResetStackAndState)
  G4bool fExiting = false; 
  G4VPhysicalVolume * fBlockedPhysicalVolume=0;
  G4int fBlockedReplicaNo=-1;
  G4bool fLocatedOnEdge=false;  
  
  G4bool notKnownContained=true,noResult;
  G4VPhysicalVolume *targetPhysical;
  G4LogicalVolume *targetLogical;
  G4VSolid *targetSolid;
  G4ThreeVector localPoint,localDirection;
  EInside insideCode;
  
  
  //Helpers/Utility classes
  G4NormalNavigation  fnormalNav;
  G4VoxelNavigation fvoxelNav;
  G4ParameterisedNavigation fparamNav;
  G4ReplicaNavigation freplicaNav;	  
  G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
  
  //local variables
  localPoint=fHistory->GetTopTransform().TransformPoint(globalPoint);
  localDirection=fHistory->GetTopTransform().TransformPoint(*pGlobalDirection);
  
  //search if volume contains point
  if (fHistory->GetTopVolumeType()!=kReplica) {
    targetSolid=fHistory->GetTopVolume()->GetLogicalVolume()->GetSolid();
    insideCode=targetSolid->Inside(localPoint);
  }
  else {
    insideCode=freplicaNav.BackLocate(*fHistory,globalPoint,
				      localPoint,fExiting,notKnownContained);
    // !CARE! if notKnownContained returns false then the point is within
    // the containing placement volume of the replica(s). If insidecode
    // will result in the history being backed up one level, then the
    // local point returned is the point in the system of this new level
#ifdef G4GEOMETRY_DEBUG
    G4cout << "Replica: fExiting=" << fExiting << G4endl;
    G4cout << "         notKnownContained=" << notKnownContained << G4endl;
#endif
  }
  
  if (insideCode==kOutside) 
    belongVolume = false;

  else if (insideCode==kSurface) {
    belongVolume = false;
    reg = -100;
  }
  else {
    // Search downwards in daughter volumes
    // until deepest containing volume found
    //
    // 3 Cases:
    //
    // o Parameterised daughters
    //   =>Must be one G4PVParameterised daughter & voxels
    // o Positioned daughters & voxels
    // o Positioned daughters & no voxels
    
    belongVolume = true;    
    noResult=true;  
    
    do {
      // Determine `type' of current mother volume
      targetPhysical=fHistory->GetTopVolume();
      targetLogical=targetPhysical->GetLogicalVolume();
      switch(CharacteriseDaughters(targetLogical)) {
      case kNormal:
	if (targetLogical->GetVoxelHeader()) {
	  noResult=fvoxelNav.LevelLocate(*fHistory,
					 fBlockedPhysicalVolume,
					 fBlockedReplicaNo,
					 globalPoint,
					 pGlobalDirection,
					 fLocatedOnEdge,
					 localPoint);
	}
	else {
	  noResult=fnormalNav.LevelLocate(*fHistory,
					  fBlockedPhysicalVolume,
					  fBlockedReplicaNo,
					  globalPoint,
					  pGlobalDirection,
					  fLocatedOnEdge,
					  localPoint);
	}
	break;
      case kReplica:
	noResult=freplicaNav.LevelLocate(*fHistory,
					 fBlockedPhysicalVolume,
					 fBlockedReplicaNo,
					 globalPoint,
					 pGlobalDirection,
					 fLocatedOnEdge,
					 localPoint);
	break;
      case kParameterised:
	noResult=fparamNav.LevelLocate(*fHistory,
				       fBlockedPhysicalVolume,
				       fBlockedReplicaNo,
				       globalPoint,
				       pGlobalDirection,
				       fLocatedOnEdge,
				       localPoint);
	break;
      } //switch
      
      // LevelLocate search in the first daughter level. 
      // LevelLocate returns noResult=true if it finds a daughter volume
      // in which globalPoint is inside (or on the surface). So point
      // doesn't belong only to mother volume ==> belongVolume=false.
      
      if (noResult) {
	// The blocked volume no longer valid - it was for another level
	fBlockedPhysicalVolume= 0;
	fBlockedReplicaNo= -1;
	belongVolume=false;
      }
    } while (noResult);
    
    //on exit targetPhysical is the volume globalPoint belongs to;
    G4int volIndex = ~0;
    for (unsigned int i=0; i<pVolStore->size(); i++)
      if ((*pVolStore)[i] == targetPhysical) 
	volIndex = i;     
    if (volIndex==(~0)) {
      G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
      exit(-999);
    }
    //reg = G4int(pVolStore->index(targetPhysical))+1;
    reg = volIndex+1;
    
  }
  
  delete fHistory;
  return belongVolume;
}

EVolume CharacteriseDaughters(const G4LogicalVolume *pLog)
{
  EVolume type;
  EAxis axis;
  G4int nReplicas;
  G4double width,offset;
  G4bool consuming;
  G4VPhysicalVolume *pVol;
  
  if (pLog->GetNoDaughters()==1) {
    pVol=pLog->GetDaughter(0);
    if (pVol->IsReplicated()) {
      pVol->GetReplicationData(axis,nReplicas,width,offset,consuming);
      type=(consuming) ? kReplica : kParameterised;
    }
    else
      type=kNormal;
  }
  else
    type=kNormal;
  return type;
}
Commit	Line	Data
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
	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