1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4VTwistSurface.cc,v 1.9 2007/05/31 13:52:48 gcosmo Exp $ |
---|
28 | // GEANT4 tag $Name: HEAD $ |
---|
29 | // |
---|
30 | // |
---|
31 | // -------------------------------------------------------------------- |
---|
32 | // GEANT 4 class source file |
---|
33 | // |
---|
34 | // |
---|
35 | // G4VTwistSurface.cc |
---|
36 | // |
---|
37 | // Author: |
---|
38 | // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp) |
---|
39 | // |
---|
40 | // History: |
---|
41 | // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4 |
---|
42 | // from original version in Jupiter-2.5.02 application. |
---|
43 | // -------------------------------------------------------------------- |
---|
44 | |
---|
45 | #include <iomanip> |
---|
46 | |
---|
47 | #include "G4VTwistSurface.hh" |
---|
48 | #include "G4GeometryTolerance.hh" |
---|
49 | |
---|
50 | const G4int G4VTwistSurface::sOutside = 0x00000000; |
---|
51 | const G4int G4VTwistSurface::sInside = 0x10000000; |
---|
52 | const G4int G4VTwistSurface::sBoundary = 0x20000000; |
---|
53 | const G4int G4VTwistSurface::sCorner = 0x40000000; |
---|
54 | const G4int G4VTwistSurface::sC0Min1Min = 0x40000101; |
---|
55 | const G4int G4VTwistSurface::sC0Max1Min = 0x40000201; |
---|
56 | const G4int G4VTwistSurface::sC0Max1Max = 0x40000202; |
---|
57 | const G4int G4VTwistSurface::sC0Min1Max = 0x40000102; |
---|
58 | const G4int G4VTwistSurface::sAxisMin = 0x00000101; |
---|
59 | const G4int G4VTwistSurface::sAxisMax = 0x00000202; |
---|
60 | const G4int G4VTwistSurface::sAxisX = 0x00000404; |
---|
61 | const G4int G4VTwistSurface::sAxisY = 0x00000808; |
---|
62 | const G4int G4VTwistSurface::sAxisZ = 0x00000C0C; |
---|
63 | const G4int G4VTwistSurface::sAxisRho = 0x00001010; |
---|
64 | const G4int G4VTwistSurface::sAxisPhi = 0x00001414; |
---|
65 | |
---|
66 | // mask |
---|
67 | const G4int G4VTwistSurface::sAxis0 = 0x0000FF00; |
---|
68 | const G4int G4VTwistSurface::sAxis1 = 0x000000FF; |
---|
69 | const G4int G4VTwistSurface::sSizeMask = 0x00000303; |
---|
70 | const G4int G4VTwistSurface::sAxisMask = 0x0000FCFC; |
---|
71 | const G4int G4VTwistSurface::sAreaMask = 0XF0000000; |
---|
72 | |
---|
73 | //===================================================================== |
---|
74 | //* constructors ------------------------------------------------------ |
---|
75 | |
---|
76 | G4VTwistSurface::G4VTwistSurface(const G4String &name) |
---|
77 | : fName(name) |
---|
78 | { |
---|
79 | |
---|
80 | fAxis[0] = kUndefined; |
---|
81 | fAxis[1] = kUndefined; |
---|
82 | fAxisMin[0] = kInfinity; |
---|
83 | fAxisMin[1] = kInfinity; |
---|
84 | fAxisMax[0] = kInfinity; |
---|
85 | fAxisMax[1] = kInfinity; |
---|
86 | fHandedness = 1; |
---|
87 | |
---|
88 | G4int i; |
---|
89 | for (i=0; i<4; i++) { |
---|
90 | fCorners[i].set(kInfinity, kInfinity, kInfinity); |
---|
91 | fNeighbours[i] = 0; |
---|
92 | } |
---|
93 | |
---|
94 | fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity); |
---|
95 | |
---|
96 | fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity); |
---|
97 | fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity); |
---|
98 | kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
99 | } |
---|
100 | |
---|
101 | G4VTwistSurface::G4VTwistSurface(const G4String &name, |
---|
102 | const G4RotationMatrix &rot, |
---|
103 | const G4ThreeVector &tlate, |
---|
104 | G4int handedness, |
---|
105 | const EAxis axis0 , |
---|
106 | const EAxis axis1 , |
---|
107 | G4double axis0min, |
---|
108 | G4double axis1min, |
---|
109 | G4double axis0max, |
---|
110 | G4double axis1max ) |
---|
111 | : fName(name) |
---|
112 | { |
---|
113 | fAxis[0] = axis0; |
---|
114 | fAxis[1] = axis1; |
---|
115 | fAxisMin[0] = axis0min; |
---|
116 | fAxisMin[1] = axis1min; |
---|
117 | fAxisMax[0] = axis0max; |
---|
118 | fAxisMax[1] = axis1max; |
---|
119 | fHandedness = handedness; |
---|
120 | fRot = rot; |
---|
121 | fTrans = tlate; |
---|
122 | |
---|
123 | G4int i; |
---|
124 | for (i=0; i<4; i++) { |
---|
125 | fCorners[i].set(kInfinity, kInfinity, kInfinity); |
---|
126 | fNeighbours[i] = 0; |
---|
127 | } |
---|
128 | |
---|
129 | fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity); |
---|
130 | |
---|
131 | fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity); |
---|
132 | fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity); |
---|
133 | kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
134 | } |
---|
135 | |
---|
136 | //===================================================================== |
---|
137 | //* Fake default constructor ------------------------------------------ |
---|
138 | |
---|
139 | G4VTwistSurface::G4VTwistSurface( __void__& ) |
---|
140 | : fName("") |
---|
141 | { |
---|
142 | } |
---|
143 | |
---|
144 | //===================================================================== |
---|
145 | //* destructor -------------------------------------------------------- |
---|
146 | |
---|
147 | G4VTwistSurface::~G4VTwistSurface() |
---|
148 | { |
---|
149 | } |
---|
150 | |
---|
151 | //===================================================================== |
---|
152 | //* AmIOnLeftSide ----------------------------------------------------- |
---|
153 | |
---|
154 | G4int G4VTwistSurface::AmIOnLeftSide(const G4ThreeVector &me, |
---|
155 | const G4ThreeVector &vec, |
---|
156 | G4bool withtol) |
---|
157 | { |
---|
158 | // AmIOnLeftSide returns phi-location of "me" |
---|
159 | // (phi relation between me and vec projected on z=0 plane). |
---|
160 | // If "me" is on -ve-phi-side of "vec", it returns 1. |
---|
161 | // On the other hand, if "me" is on +ve-phi-side of "vec", |
---|
162 | // it returns -1. |
---|
163 | // (The return value represents z-coordinate of normal vector |
---|
164 | // of me.cross(vec).) |
---|
165 | // If me is on boundary of vec, return 0. |
---|
166 | |
---|
167 | static const G4double kAngTolerance |
---|
168 | = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); |
---|
169 | |
---|
170 | static G4RotationMatrix unitrot; // unit matrix |
---|
171 | static const G4RotationMatrix rottol = unitrot.rotateZ(0.5*kAngTolerance); |
---|
172 | static const G4RotationMatrix invrottol = unitrot.rotateZ(-1.*kAngTolerance); |
---|
173 | |
---|
174 | if (fAmIOnLeftSide.me == me |
---|
175 | && fAmIOnLeftSide.vec == vec |
---|
176 | && fAmIOnLeftSide.withTol == withtol) { |
---|
177 | return fAmIOnLeftSide.amIOnLeftSide; |
---|
178 | } |
---|
179 | |
---|
180 | fAmIOnLeftSide.me = me; |
---|
181 | fAmIOnLeftSide.vec = vec; |
---|
182 | fAmIOnLeftSide.withTol = withtol; |
---|
183 | |
---|
184 | G4ThreeVector met = (G4ThreeVector(me.x(), me.y(), 0.)).unit(); |
---|
185 | G4ThreeVector vect = (G4ThreeVector(vec.x(), vec.y(), 0.)).unit(); |
---|
186 | |
---|
187 | G4ThreeVector ivect = invrottol * vect; |
---|
188 | G4ThreeVector rvect = rottol * vect; |
---|
189 | |
---|
190 | G4double metcrossvect = met.x() * vect.y() - met.y() * vect.x(); |
---|
191 | |
---|
192 | if (withtol) { |
---|
193 | if (met.x() * ivect.y() - met.y() * ivect.x() > 0 && |
---|
194 | metcrossvect >= 0) { |
---|
195 | fAmIOnLeftSide.amIOnLeftSide = 1; |
---|
196 | } else if (met.x() * rvect.y() - met.y() * rvect.x() < 0 && |
---|
197 | metcrossvect <= 0) { |
---|
198 | fAmIOnLeftSide.amIOnLeftSide = -1; |
---|
199 | } else { |
---|
200 | fAmIOnLeftSide.amIOnLeftSide = 0; |
---|
201 | } |
---|
202 | } else { |
---|
203 | if (metcrossvect > 0) { |
---|
204 | fAmIOnLeftSide.amIOnLeftSide = 1; |
---|
205 | } else if (metcrossvect < 0 ) { |
---|
206 | fAmIOnLeftSide.amIOnLeftSide = -1; |
---|
207 | } else { |
---|
208 | fAmIOnLeftSide.amIOnLeftSide = 0; |
---|
209 | } |
---|
210 | } |
---|
211 | |
---|
212 | #ifdef G4TWISTDEBUG |
---|
213 | G4cout << " === G4VTwistSurface::AmIOnLeftSide() ==============" |
---|
214 | << G4endl; |
---|
215 | G4cout << " Name , returncode : " << fName << " " |
---|
216 | << fAmIOnLeftSide.amIOnLeftSide << G4endl; |
---|
217 | G4cout << " me, vec : " << std::setprecision(14) << me |
---|
218 | << " " << vec << G4endl; |
---|
219 | G4cout << " met, vect : " << met << " " << vect << G4endl; |
---|
220 | G4cout << " ivec, rvec : " << ivect << " " << rvect << G4endl; |
---|
221 | G4cout << " met x vect : " << metcrossvect << G4endl; |
---|
222 | G4cout << " met x ivec : " << met.cross(ivect) << G4endl; |
---|
223 | G4cout << " met x rvec : " << met.cross(rvect) << G4endl; |
---|
224 | G4cout << " ==============================================" |
---|
225 | << G4endl; |
---|
226 | #endif |
---|
227 | |
---|
228 | return fAmIOnLeftSide.amIOnLeftSide; |
---|
229 | } |
---|
230 | |
---|
231 | //===================================================================== |
---|
232 | //* DistanceToBoundary ------------------------------------------------ |
---|
233 | |
---|
234 | G4double G4VTwistSurface::DistanceToBoundary(G4int areacode, |
---|
235 | G4ThreeVector &xx, |
---|
236 | const G4ThreeVector &p) |
---|
237 | { |
---|
238 | // DistanceToBoundary |
---|
239 | // |
---|
240 | // return distance to nearest boundary from arbitrary point p |
---|
241 | // in local coodinate. |
---|
242 | // Argument areacode must be one of them: |
---|
243 | // sAxis0 & sAxisMin, sAxis0 & sAxisMax, |
---|
244 | // sAxis1 & sAxisMin, sAxis1 & sAxisMax. |
---|
245 | // |
---|
246 | |
---|
247 | G4ThreeVector d; // direction vector of the boundary |
---|
248 | G4ThreeVector x0; // reference point of the boundary |
---|
249 | G4double dist = kInfinity; |
---|
250 | G4int boundarytype; |
---|
251 | |
---|
252 | if (IsAxis0(areacode) && IsAxis1(areacode)) { |
---|
253 | G4cerr << "ERROR - G4VTwistSurface::DistanceToBoundary()" << G4endl |
---|
254 | << " Point is in the corner area. This function returns" |
---|
255 | << G4endl |
---|
256 | << " a direction vector of a boundary line." << G4endl |
---|
257 | << " areacode = " << areacode << G4endl; |
---|
258 | G4Exception("G4VTwistSurface::DistanceToBoundary()", "InvalidSetup", |
---|
259 | FatalException, "Point is in the corner area."); |
---|
260 | } else if (IsAxis0(areacode) || IsAxis1(areacode)) { |
---|
261 | GetBoundaryParameters(areacode, d, x0, boundarytype); |
---|
262 | if (boundarytype == sAxisPhi) { |
---|
263 | G4double t = x0.getRho() / p.getRho(); |
---|
264 | xx.set(t*p.x(), t*p.y(), x0.z()); |
---|
265 | dist = (xx - p).mag(); |
---|
266 | } else { |
---|
267 | // linear boundary |
---|
268 | // sAxisX, sAxisY, sAxisZ, sAxisRho |
---|
269 | dist = DistanceToLine(p, x0, d, xx); |
---|
270 | } |
---|
271 | } else { |
---|
272 | G4cerr << "ERROR - G4VTwistSurface::DistanceToBoundary()" << G4endl |
---|
273 | << " areacode = " << areacode << G4endl; |
---|
274 | G4Exception("G4VTwistSurface::DistanceToBoundary()", "InvalidSetup", |
---|
275 | FatalException, "Bad areacode of boundary."); |
---|
276 | } |
---|
277 | return dist; |
---|
278 | } |
---|
279 | |
---|
280 | //===================================================================== |
---|
281 | //* DistanceToIn ------------------------------------------------------ |
---|
282 | |
---|
283 | G4double G4VTwistSurface::DistanceToIn(const G4ThreeVector &gp, |
---|
284 | const G4ThreeVector &gv, |
---|
285 | G4ThreeVector &gxxbest) |
---|
286 | { |
---|
287 | #ifdef G4TWISTDEBUG |
---|
288 | G4cout << " ~~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~~" << G4endl; |
---|
289 | G4cout << " Name : " << fName << G4endl; |
---|
290 | G4cout << " gp : " << gp << G4endl; |
---|
291 | G4cout << " gv : " << gv << G4endl; |
---|
292 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
293 | #endif |
---|
294 | |
---|
295 | G4ThreeVector gxx[G4VSURFACENXX]; |
---|
296 | G4double distance[G4VSURFACENXX] ; |
---|
297 | G4int areacode[G4VSURFACENXX] ; |
---|
298 | G4bool isvalid[G4VSURFACENXX] ; |
---|
299 | |
---|
300 | for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) { |
---|
301 | distance[i] = kInfinity ; |
---|
302 | areacode[i] = sOutside ; |
---|
303 | isvalid[i] = false ; |
---|
304 | } |
---|
305 | |
---|
306 | G4double bestdistance = kInfinity; |
---|
307 | G4int besti = -1; |
---|
308 | G4ThreeVector bestgxx(kInfinity, kInfinity, kInfinity); |
---|
309 | |
---|
310 | G4int nxx = DistanceToSurface(gp, gv, gxx, distance, areacode, |
---|
311 | isvalid, kValidateWithTol); |
---|
312 | |
---|
313 | for (G4int i=0; i< nxx; i++) { |
---|
314 | |
---|
315 | // skip this intersection if: |
---|
316 | // - invalid intersection |
---|
317 | // - particle goes outword the surface |
---|
318 | |
---|
319 | if (!isvalid[i]) { |
---|
320 | // xx[i] is sOutside or distance[i] < 0 |
---|
321 | continue; |
---|
322 | } |
---|
323 | |
---|
324 | G4ThreeVector normal = GetNormal(gxx[i], true); |
---|
325 | |
---|
326 | if ((normal * gv) >= 0) { |
---|
327 | |
---|
328 | #ifdef G4TWISTDEBUG |
---|
329 | G4cout << " G4VTwistSurface::DistanceToIn(p,v): " |
---|
330 | << "particle goes outword the surface." << G4endl; |
---|
331 | #endif |
---|
332 | continue; |
---|
333 | } |
---|
334 | |
---|
335 | // |
---|
336 | // accept this intersection if the intersection is inside. |
---|
337 | // |
---|
338 | |
---|
339 | if (IsInside(areacode[i])) { |
---|
340 | if (distance[i] < bestdistance) { |
---|
341 | bestdistance = distance[i]; |
---|
342 | bestgxx = gxx[i]; |
---|
343 | besti = i; |
---|
344 | |
---|
345 | #ifdef G4TWISTDEBUG |
---|
346 | G4cout << " G4VTwistSurface::DistanceToIn(p,v): " |
---|
347 | << " areacode sInside name, distance = " |
---|
348 | << fName << " "<< bestdistance << G4endl; |
---|
349 | #endif |
---|
350 | } |
---|
351 | |
---|
352 | // |
---|
353 | // else, the intersection is on boundary or corner. |
---|
354 | // |
---|
355 | |
---|
356 | } else { |
---|
357 | |
---|
358 | G4VTwistSurface *neighbours[2]; |
---|
359 | G4bool isaccepted[2] = {false, false}; |
---|
360 | G4int nneighbours = GetNeighbours(areacode[i], neighbours); |
---|
361 | |
---|
362 | for (G4int j=0; j< nneighbours; j++) { |
---|
363 | // if on corner, nneighbours = 2. |
---|
364 | // if on boundary, nneighbours = 1. |
---|
365 | |
---|
366 | G4ThreeVector tmpgxx[G4VSURFACENXX]; |
---|
367 | G4double tmpdist[G4VSURFACENXX] ; |
---|
368 | G4int tmpareacode[G4VSURFACENXX] ; |
---|
369 | G4bool tmpisvalid[G4VSURFACENXX] ; |
---|
370 | |
---|
371 | for (G4int l = 0 ; l<G4VSURFACENXX ; l++ ) { |
---|
372 | tmpdist[l] = kInfinity ; |
---|
373 | tmpareacode[l] = sOutside ; |
---|
374 | tmpisvalid[l] = false ; |
---|
375 | } |
---|
376 | |
---|
377 | G4int tmpnxx = neighbours[j]->DistanceToSurface( |
---|
378 | gp, gv, tmpgxx, tmpdist, |
---|
379 | tmpareacode, tmpisvalid, |
---|
380 | kValidateWithTol); |
---|
381 | G4ThreeVector neighbournormal; |
---|
382 | |
---|
383 | for (G4int k=0; k< tmpnxx; k++) { |
---|
384 | |
---|
385 | // |
---|
386 | // if tmpxx[k] is valid && sInside, the final winner must |
---|
387 | // be neighbour surface. return kInfinity. |
---|
388 | // else , choose tmpxx on same boundary of xx, then check normal |
---|
389 | // |
---|
390 | |
---|
391 | if (IsInside(tmpareacode[k])) { |
---|
392 | |
---|
393 | #ifdef G4TWISTDEBUG |
---|
394 | G4cout << " G4VTwistSurface:DistanceToIn(p,v): " |
---|
395 | << " intersection "<< tmpgxx[k] << G4endl |
---|
396 | << " is inside of neighbour surface of " << fName |
---|
397 | << " . returning kInfinity." << G4endl; |
---|
398 | G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" |
---|
399 | << G4endl; |
---|
400 | G4cout << " No intersections " << G4endl; |
---|
401 | G4cout << " Name : " << fName << G4endl; |
---|
402 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" |
---|
403 | << G4endl; |
---|
404 | #endif |
---|
405 | if (tmpisvalid[k]) return kInfinity; |
---|
406 | continue; |
---|
407 | |
---|
408 | // |
---|
409 | // if tmpxx[k] is valid && sInside, the final winner must |
---|
410 | // be neighbour surface. return . |
---|
411 | // |
---|
412 | |
---|
413 | } else if (IsSameBoundary(this,areacode[i], |
---|
414 | neighbours[j], tmpareacode[k])) { |
---|
415 | // tmpxx[k] is same boundary (or corner) of xx. |
---|
416 | |
---|
417 | neighbournormal = neighbours[j]->GetNormal(tmpgxx[k], true); |
---|
418 | if (neighbournormal * gv < 0) isaccepted[j] = true; |
---|
419 | } |
---|
420 | } |
---|
421 | |
---|
422 | // if nneighbours = 1, chabge isaccepted[1] before |
---|
423 | // exiting neighboursurface loop. |
---|
424 | |
---|
425 | if (nneighbours == 1) isaccepted[1] = true; |
---|
426 | |
---|
427 | } // neighboursurface loop end |
---|
428 | |
---|
429 | // now, we can accept xx intersection |
---|
430 | |
---|
431 | if (isaccepted[0] == true && isaccepted[1] == true) { |
---|
432 | if (distance[i] < bestdistance) { |
---|
433 | bestdistance = distance[i]; |
---|
434 | gxxbest = gxx[i]; |
---|
435 | besti = i; |
---|
436 | #ifdef G4TWISTDEBUG |
---|
437 | G4cout << " G4VTwistSurface::DistanceToIn(p,v): " |
---|
438 | << " areacode sBoundary & sBoundary distance = " |
---|
439 | << fName << " " << distance[i] << G4endl; |
---|
440 | #endif |
---|
441 | } |
---|
442 | } |
---|
443 | |
---|
444 | } // else end |
---|
445 | } // intersection loop end |
---|
446 | |
---|
447 | gxxbest = bestgxx; |
---|
448 | |
---|
449 | #ifdef G4TWISTDEBUG |
---|
450 | if (besti < 0) { |
---|
451 | G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" << G4endl; |
---|
452 | G4cout << " No intersections " << G4endl; |
---|
453 | G4cout << " Name : " << fName << G4endl; |
---|
454 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
455 | } else { |
---|
456 | G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" << G4endl; |
---|
457 | G4cout << " Name, i : " << fName << " , " << besti << G4endl; |
---|
458 | G4cout << " gxx[i] : " << gxxbest << G4endl; |
---|
459 | G4cout << " bestdist : " << bestdistance << G4endl; |
---|
460 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
461 | } |
---|
462 | |
---|
463 | #endif |
---|
464 | |
---|
465 | return bestdistance; |
---|
466 | } |
---|
467 | |
---|
468 | //===================================================================== |
---|
469 | //* DistanceToOut(p, v) ----------------------------------------------- |
---|
470 | |
---|
471 | G4double G4VTwistSurface::DistanceToOut(const G4ThreeVector &gp, |
---|
472 | const G4ThreeVector &gv, |
---|
473 | G4ThreeVector &gxxbest) |
---|
474 | { |
---|
475 | #ifdef G4TWISTDEBUG |
---|
476 | G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" << G4endl; |
---|
477 | G4cout << " Name : " << fName << G4endl; |
---|
478 | G4cout << " gp : " << gp << G4endl; |
---|
479 | G4cout << " gv : " << gv << G4endl; |
---|
480 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
481 | #endif |
---|
482 | |
---|
483 | G4ThreeVector gxx[G4VSURFACENXX]; |
---|
484 | G4double distance[G4VSURFACENXX] ; |
---|
485 | G4int areacode[G4VSURFACENXX] ; |
---|
486 | G4bool isvalid[G4VSURFACENXX] ; |
---|
487 | G4int i; |
---|
488 | |
---|
489 | for ( i = 0 ; i<G4VSURFACENXX ; i++ ) |
---|
490 | { |
---|
491 | distance[i] = kInfinity ; |
---|
492 | areacode[i] = sOutside ; |
---|
493 | isvalid[i] = false ; |
---|
494 | } |
---|
495 | |
---|
496 | G4int nxx; |
---|
497 | G4double bestdistance = kInfinity; |
---|
498 | G4int besti = -1; |
---|
499 | |
---|
500 | nxx = DistanceToSurface(gp, gv, gxx, distance, areacode, |
---|
501 | isvalid, kValidateWithTol); |
---|
502 | |
---|
503 | for (i=0; i<nxx; i++) { |
---|
504 | if (!(isvalid[i])) { |
---|
505 | continue; |
---|
506 | } |
---|
507 | |
---|
508 | G4ThreeVector normal = GetNormal(gxx[i], true); |
---|
509 | if (normal * gv <= 0) { |
---|
510 | // particle goes toword inside of solid, return kInfinity |
---|
511 | #ifdef G4TWISTDEBUG |
---|
512 | G4cout << " G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0, normal " |
---|
513 | << fName << " " << normal |
---|
514 | << G4endl; |
---|
515 | #endif |
---|
516 | } else { |
---|
517 | // gxx[i] is accepted. |
---|
518 | if (distance[i] < bestdistance) { |
---|
519 | bestdistance = distance[i]; |
---|
520 | gxxbest = gxx[i]; |
---|
521 | besti = i; |
---|
522 | } |
---|
523 | } |
---|
524 | } |
---|
525 | |
---|
526 | #ifdef G4TWISTDEBUG |
---|
527 | if (besti < 0) { |
---|
528 | G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~~" << G4endl; |
---|
529 | G4cout << " No intersections " << G4endl; |
---|
530 | G4cout << " Name : " << fName << G4endl; |
---|
531 | G4cout << " bestdist : " << bestdistance << G4endl; |
---|
532 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
533 | } else { |
---|
534 | G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~~" << G4endl; |
---|
535 | G4cout << " Name, i : " << fName << " , " << i << G4endl; |
---|
536 | G4cout << " gxx[i] : " << gxxbest << G4endl; |
---|
537 | G4cout << " bestdist : " << bestdistance << G4endl; |
---|
538 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
539 | } |
---|
540 | #endif |
---|
541 | |
---|
542 | return bestdistance; |
---|
543 | } |
---|
544 | |
---|
545 | //===================================================================== |
---|
546 | //* DistanceTo(p) ----------------------------------------------------- |
---|
547 | |
---|
548 | G4double G4VTwistSurface::DistanceTo(const G4ThreeVector &gp, |
---|
549 | G4ThreeVector &gxxbest) |
---|
550 | { |
---|
551 | #ifdef G4TWISTDEBUG |
---|
552 | G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" << G4endl; |
---|
553 | G4cout << " Name : " << fName << G4endl; |
---|
554 | G4cout << " gp : " << gp << G4endl; |
---|
555 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
556 | #endif |
---|
557 | |
---|
558 | |
---|
559 | G4ThreeVector gxx[G4VSURFACENXX]; |
---|
560 | G4double distance[G4VSURFACENXX] ; |
---|
561 | G4int areacode[G4VSURFACENXX] ; |
---|
562 | |
---|
563 | for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) { |
---|
564 | distance[i] = kInfinity ; |
---|
565 | areacode[i] = sOutside ; |
---|
566 | } |
---|
567 | |
---|
568 | G4int nxx; |
---|
569 | |
---|
570 | nxx = DistanceToSurface(gp, gxx, distance, areacode); |
---|
571 | gxxbest = gxx[0]; |
---|
572 | |
---|
573 | #ifdef G4TWISTDEBUG |
---|
574 | G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" << G4endl; |
---|
575 | G4cout << " Name : " << fName << G4endl; |
---|
576 | G4cout << " gxx : " << gxxbest << G4endl; |
---|
577 | G4cout << " bestdist : " << distance[0] << G4endl; |
---|
578 | G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; |
---|
579 | #endif |
---|
580 | |
---|
581 | return distance[0]; |
---|
582 | } |
---|
583 | |
---|
584 | //===================================================================== |
---|
585 | //* IsSameBoundary ---------------------------------------------------- |
---|
586 | |
---|
587 | G4bool G4VTwistSurface::IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1, |
---|
588 | G4VTwistSurface *surface2, G4int areacode2 ) const |
---|
589 | { |
---|
590 | // |
---|
591 | // IsSameBoundary |
---|
592 | // |
---|
593 | // checking tool whether two boundaries on different surfaces are same or not. |
---|
594 | // |
---|
595 | |
---|
596 | G4bool testbitmode = true; |
---|
597 | G4bool iscorner[2] = {IsCorner(areacode1, testbitmode), |
---|
598 | IsCorner(areacode2, testbitmode)}; |
---|
599 | |
---|
600 | if (iscorner[0] && iscorner[1]) { |
---|
601 | // on corner |
---|
602 | G4ThreeVector corner1 = |
---|
603 | surface1->ComputeGlobalPoint(surface1->GetCorner(areacode1)); |
---|
604 | G4ThreeVector corner2 = |
---|
605 | surface2->ComputeGlobalPoint(surface2->GetCorner(areacode2)); |
---|
606 | |
---|
607 | if ((corner1 - corner2).mag() < kCarTolerance) { |
---|
608 | return true; |
---|
609 | } else { |
---|
610 | return false; |
---|
611 | } |
---|
612 | |
---|
613 | } else if ((IsBoundary(areacode1, testbitmode) && (!iscorner[0])) && |
---|
614 | (IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) { |
---|
615 | // on boundary |
---|
616 | G4ThreeVector d1, d2, ld1, ld2; |
---|
617 | G4ThreeVector x01, x02, lx01, lx02; |
---|
618 | G4int type1, type2; |
---|
619 | surface1->GetBoundaryParameters(areacode1, ld1, lx01, type1); |
---|
620 | surface2->GetBoundaryParameters(areacode2, ld2, lx02, type2); |
---|
621 | |
---|
622 | x01 = surface1->ComputeGlobalPoint(lx01); |
---|
623 | x02 = surface2->ComputeGlobalPoint(lx02); |
---|
624 | d1 = surface1->ComputeGlobalDirection(ld1); |
---|
625 | d2 = surface2->ComputeGlobalDirection(ld2); |
---|
626 | |
---|
627 | if ((x01 - x02).mag() < kCarTolerance && |
---|
628 | (d1 - d2).mag() < kCarTolerance) { |
---|
629 | return true; |
---|
630 | } else { |
---|
631 | return false; |
---|
632 | } |
---|
633 | |
---|
634 | } else { |
---|
635 | return false; |
---|
636 | } |
---|
637 | } |
---|
638 | |
---|
639 | //===================================================================== |
---|
640 | //* GetBoundaryParameters --------------------------------------------- |
---|
641 | |
---|
642 | void G4VTwistSurface::GetBoundaryParameters(const G4int &areacode, |
---|
643 | G4ThreeVector &d, |
---|
644 | G4ThreeVector &x0, |
---|
645 | G4int &boundarytype) const |
---|
646 | { |
---|
647 | // areacode must be one of them: |
---|
648 | // sAxis0 & sAxisMin, sAxis0 & sAxisMax, |
---|
649 | // sAxis1 & sAxisMin, sAxis1 & sAxisMax. |
---|
650 | |
---|
651 | G4int i; |
---|
652 | for (i=0; i<4; i++) { |
---|
653 | if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0, |
---|
654 | boundarytype)) { |
---|
655 | return; |
---|
656 | } |
---|
657 | } |
---|
658 | |
---|
659 | G4cerr << "ERROR - G4VTwistSurface::GetBoundaryParameters()" << G4endl |
---|
660 | << " Boundary at areacode " << std::hex << areacode |
---|
661 | << std::dec << G4endl |
---|
662 | << " is not be registered." << G4endl; |
---|
663 | G4Exception("G4VTwistSurface::GetBoundaryParameters()", "InvalidSetup", |
---|
664 | FatalException, "Not registered boundary."); |
---|
665 | } |
---|
666 | |
---|
667 | //===================================================================== |
---|
668 | //* GetBoundaryAtPZ --------------------------------------------------- |
---|
669 | |
---|
670 | G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ(G4int areacode, |
---|
671 | const G4ThreeVector &p) const |
---|
672 | { |
---|
673 | // areacode must be one of them: |
---|
674 | // sAxis0 & sAxisMin, sAxis0 & sAxisMax, |
---|
675 | // sAxis1 & sAxisMin, sAxis1 & sAxisMax. |
---|
676 | |
---|
677 | if (areacode & sAxis0 && areacode & sAxis1) { |
---|
678 | G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl |
---|
679 | << " Point is in the corner area. This function returns" |
---|
680 | << G4endl |
---|
681 | << " a direction vector of a boundary line." << G4endl |
---|
682 | << " areacode = " << areacode << G4endl; |
---|
683 | G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidCondition", |
---|
684 | FatalException, "Point is in the corner area."); |
---|
685 | } |
---|
686 | |
---|
687 | G4ThreeVector d; |
---|
688 | G4ThreeVector x0; |
---|
689 | G4int boundarytype; |
---|
690 | G4bool found = false; |
---|
691 | |
---|
692 | for (G4int i=0; i<4; i++) { |
---|
693 | if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0, |
---|
694 | boundarytype)){ |
---|
695 | found = true; |
---|
696 | continue; |
---|
697 | } |
---|
698 | } |
---|
699 | |
---|
700 | if (!found) { |
---|
701 | G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl |
---|
702 | << " Boundary at areacode " << areacode << G4endl |
---|
703 | << " is not be registered." << G4endl; |
---|
704 | G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidSetup", |
---|
705 | FatalException, "Not registered boundary."); |
---|
706 | } |
---|
707 | |
---|
708 | if (((boundarytype & sAxisPhi) == sAxisPhi) || |
---|
709 | ((boundarytype & sAxisRho) == sAxisRho)) { |
---|
710 | G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAtPZ()" << G4endl |
---|
711 | << " Boundary at areacode " << areacode << G4endl |
---|
712 | << " is not a z-depended line." << G4endl; |
---|
713 | G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "InvalidSetup", |
---|
714 | FatalException, "Not a z-depended line boundary."); |
---|
715 | } |
---|
716 | return ((p.z() - x0.z()) / d.z()) * d + x0; |
---|
717 | } |
---|
718 | |
---|
719 | //===================================================================== |
---|
720 | //* SetCorner --------------------------------------------------------- |
---|
721 | |
---|
722 | void G4VTwistSurface::SetCorner(G4int areacode, G4double x, G4double y, G4double z) |
---|
723 | { |
---|
724 | if ((areacode & sCorner) != sCorner){ |
---|
725 | G4cerr << "ERROR - G4VTwistSurface::SetCorner()" << G4endl |
---|
726 | << " areacode " << areacode << G4endl; |
---|
727 | G4Exception("G4VTwistSurface::SetCorner()", "InvalidSetup", |
---|
728 | FatalException, "Area code must represents corner."); |
---|
729 | } |
---|
730 | |
---|
731 | if ((areacode & sC0Min1Min) == sC0Min1Min) { |
---|
732 | fCorners[0].set(x, y, z); |
---|
733 | } else if ((areacode & sC0Max1Min) == sC0Max1Min) { |
---|
734 | fCorners[1].set(x, y, z); |
---|
735 | } else if ((areacode & sC0Max1Max) == sC0Max1Max) { |
---|
736 | fCorners[2].set(x, y, z); |
---|
737 | } else if ((areacode & sC0Min1Max) == sC0Min1Max) { |
---|
738 | fCorners[3].set(x, y, z); |
---|
739 | } |
---|
740 | } |
---|
741 | |
---|
742 | //===================================================================== |
---|
743 | //* SetBoundaryAxis --------------------------------------------------- |
---|
744 | |
---|
745 | void G4VTwistSurface::GetBoundaryAxis(G4int areacode, EAxis axis[]) const |
---|
746 | { |
---|
747 | if ((areacode & sBoundary) != sBoundary) { |
---|
748 | G4Exception("G4VTwistSurface::GetBoundaryAxis()", "InvalidCondition", |
---|
749 | FatalException, "Not located on a boundary!"); |
---|
750 | } |
---|
751 | G4int i; |
---|
752 | for (i=0; i<2; i++) { |
---|
753 | |
---|
754 | G4int whichaxis = 0 ; |
---|
755 | if (i == 0) { |
---|
756 | whichaxis = sAxis0; |
---|
757 | } else if (i == 1) { |
---|
758 | whichaxis = sAxis1; |
---|
759 | } |
---|
760 | |
---|
761 | // extracted axiscode of whichaxis |
---|
762 | G4int axiscode = whichaxis & sAxisMask & areacode ; |
---|
763 | if (axiscode) { |
---|
764 | if (axiscode == (whichaxis & sAxisX)) { |
---|
765 | axis[i] = kXAxis; |
---|
766 | } else if (axiscode == (whichaxis & sAxisY)) { |
---|
767 | axis[i] = kYAxis; |
---|
768 | } else if (axiscode == (whichaxis & sAxisZ)) { |
---|
769 | axis[i] = kZAxis; |
---|
770 | } else if (axiscode == (whichaxis & sAxisRho)) { |
---|
771 | axis[i] = kRho; |
---|
772 | } else if (axiscode == (whichaxis & sAxisPhi)) { |
---|
773 | axis[i] = kPhi; |
---|
774 | } else { |
---|
775 | G4cerr << "ERROR - G4VTwistSurface::GetBoundaryAxis()" << G4endl |
---|
776 | << " areacode " << areacode << G4endl; |
---|
777 | G4Exception("G4VTwistSurface::GetBoundaryAxis()", "InvalidSetup", |
---|
778 | FatalException, "Not supported areacode."); |
---|
779 | } |
---|
780 | } |
---|
781 | } |
---|
782 | } |
---|
783 | |
---|
784 | //===================================================================== |
---|
785 | //* SetBoundaryLimit -------------------------------------------------- |
---|
786 | |
---|
787 | void G4VTwistSurface::GetBoundaryLimit(G4int areacode, G4double limit[]) const |
---|
788 | { |
---|
789 | if (areacode & sCorner) { |
---|
790 | if (areacode & sC0Min1Max) { |
---|
791 | limit[0] = fAxisMin[0]; |
---|
792 | limit[1] = fAxisMin[1]; |
---|
793 | } else if (areacode & sC0Max1Min) { |
---|
794 | limit[0] = fAxisMax[0]; |
---|
795 | limit[1] = fAxisMin[1]; |
---|
796 | } else if (areacode & sC0Max1Max) { |
---|
797 | limit[0] = fAxisMax[0]; |
---|
798 | limit[1] = fAxisMax[1]; |
---|
799 | } else if (areacode & sC0Min1Max) { |
---|
800 | limit[0] = fAxisMin[0]; |
---|
801 | limit[1] = fAxisMax[1]; |
---|
802 | } |
---|
803 | } else if (areacode & sBoundary) { |
---|
804 | if (areacode & (sAxis0 | sAxisMin)) { |
---|
805 | limit[0] = fAxisMin[0]; |
---|
806 | } else if (areacode & (sAxis1 | sAxisMin)) { |
---|
807 | limit[0] = fAxisMin[1]; |
---|
808 | } else if (areacode & (sAxis0 | sAxisMax)) { |
---|
809 | limit[0] = fAxisMax[0]; |
---|
810 | } else if (areacode & (sAxis1 | sAxisMax)) { |
---|
811 | limit[0] = fAxisMax[1]; |
---|
812 | } |
---|
813 | } else { |
---|
814 | G4cerr << "WARNING - G4VTwistSurface::GetBoundaryAxis()" << G4endl |
---|
815 | << " areacode " << areacode << G4endl; |
---|
816 | G4Exception("G4VTwistSurface::GetBoundaryLimit()", "InvalidCondition", |
---|
817 | JustWarning, "Not located on a boundary!"); |
---|
818 | } |
---|
819 | } |
---|
820 | |
---|
821 | //===================================================================== |
---|
822 | //* SetBoundary ------------------------------------------------------- |
---|
823 | |
---|
824 | void G4VTwistSurface::SetBoundary(const G4int &axiscode, |
---|
825 | const G4ThreeVector &direction, |
---|
826 | const G4ThreeVector &x0, |
---|
827 | const G4int &boundarytype) |
---|
828 | { |
---|
829 | G4int code = (~sAxisMask) & axiscode; |
---|
830 | if ((code == (sAxis0 & sAxisMin)) || |
---|
831 | (code == (sAxis0 & sAxisMax)) || |
---|
832 | (code == (sAxis1 & sAxisMin)) || |
---|
833 | (code == (sAxis1 & sAxisMax))) { |
---|
834 | |
---|
835 | G4int i; |
---|
836 | G4bool done = false; |
---|
837 | for (i=0; i<4; i++) { |
---|
838 | if (fBoundaries[i].IsEmpty()) { |
---|
839 | fBoundaries[i].SetFields(axiscode, direction, |
---|
840 | x0, boundarytype); |
---|
841 | done = true; |
---|
842 | break; |
---|
843 | } |
---|
844 | } |
---|
845 | |
---|
846 | if (!done) { |
---|
847 | G4Exception("G4VTwistSurface::SetBoundary()", "InvalidCondition", |
---|
848 | FatalException, "Number of boundary exceeding 4!"); |
---|
849 | } |
---|
850 | } else { |
---|
851 | G4cerr << "ERROR - G4VTwistSurface::SetBoundary()" << G4endl |
---|
852 | << " invalid axiscode. axiscode = " |
---|
853 | << std::hex << axiscode << std::dec << G4endl; |
---|
854 | G4Exception("G4VTwistSurface::SetBoundary()", "InvalidCondition", |
---|
855 | FatalException, "Invalid axis-code."); |
---|
856 | } |
---|
857 | } |
---|
858 | |
---|
859 | //===================================================================== |
---|
860 | //* GetFace ----------------------------------------------------------- |
---|
861 | |
---|
862 | G4int G4VTwistSurface::GetFace( G4int i, G4int j, G4int m, |
---|
863 | G4int n, G4int iside ) |
---|
864 | { |
---|
865 | // this is the face mapping function |
---|
866 | // (i,j) -> face number |
---|
867 | |
---|
868 | if ( iside == 0 ) { |
---|
869 | return i * ( m - 1 ) + j ; |
---|
870 | } |
---|
871 | |
---|
872 | else if ( iside == 1 ) { |
---|
873 | return (m-1)*(m-1) + i*(m-1) + j ; |
---|
874 | } |
---|
875 | |
---|
876 | else if ( iside == 2 ) { |
---|
877 | return 2*(m-1)*(m-1) + i*(m-1) + j ; |
---|
878 | } |
---|
879 | |
---|
880 | else if ( iside == 3 ) { |
---|
881 | return 2*(m-1)*(m-1) + (n-1)*(m-1) + i*(m-1) + j ; |
---|
882 | } |
---|
883 | |
---|
884 | else if ( iside == 4 ) { |
---|
885 | return 2*(m-1)*(m-1) + 2*(n-1)*(m-1) + i*(m-1) + j ; |
---|
886 | } |
---|
887 | |
---|
888 | else if ( iside == 5 ) { |
---|
889 | return 2*(m-1)*(m-1) + 3*(n-1)*(m-1) + i*(m-1) + j ; |
---|
890 | } |
---|
891 | |
---|
892 | else { |
---|
893 | |
---|
894 | G4cerr << "ERROR - G4VTwistSurface::GetFace(): " |
---|
895 | << GetName() << G4endl |
---|
896 | << " Not correct side number! - " << G4endl |
---|
897 | << "iside is " << iside << " but should be " |
---|
898 | << "0,1,2,3,4 or 5" << "." << G4endl ; |
---|
899 | G4Exception("G4TwistSurface::G4GetFace()", "InvalidSetup", |
---|
900 | FatalException, "Not correct side number."); |
---|
901 | |
---|
902 | |
---|
903 | } |
---|
904 | |
---|
905 | return -1 ; // wrong face |
---|
906 | } |
---|
907 | |
---|
908 | //===================================================================== |
---|
909 | //* GetNode ----------------------------------------------------------- |
---|
910 | |
---|
911 | G4int G4VTwistSurface::GetNode( G4int i, G4int j, G4int m, |
---|
912 | G4int n, G4int iside ) |
---|
913 | { |
---|
914 | // this is the node mapping function |
---|
915 | // (i,j) -> node number |
---|
916 | // Depends on the side iside and the used meshing of the surface |
---|
917 | |
---|
918 | if ( iside == 0 ) { |
---|
919 | // lower endcap is mxm squared. |
---|
920 | // n = m |
---|
921 | return i * m + j ; |
---|
922 | } |
---|
923 | |
---|
924 | if ( iside == 1 ) { |
---|
925 | // upper endcap is mxm squared. Shift by m*m |
---|
926 | // n = m |
---|
927 | return m*m + i*m + j ; |
---|
928 | } |
---|
929 | |
---|
930 | else if ( iside == 2 ) { |
---|
931 | // front side. |
---|
932 | if ( i == 0 ) { return j ; } |
---|
933 | else if ( i == n-1 ) { return m*m + j ; } |
---|
934 | else { return 2*m*m + 4*(i-1)*(m-1) + j ; } |
---|
935 | } |
---|
936 | |
---|
937 | else if ( iside == 3 ) { |
---|
938 | // right side |
---|
939 | if ( i == 0 ) { return (j+1)*m - 1 ; } |
---|
940 | else if ( i == n-1 ) { return m*m + (j+1)*m - 1 ; } |
---|
941 | else { return 2*m*m + 4*(i-1)*(m-1) + (m-1) + j ; } |
---|
942 | } |
---|
943 | else if ( iside == 4 ) { |
---|
944 | // back side. |
---|
945 | if ( i == 0 ) { return m*m - 1 - j ; } // reversed order |
---|
946 | else if ( i == n-1 ) { return 2*m*m - 1 - j ; } // reversed order |
---|
947 | else { return 2*m*m + 4*(i-1)*(m-1) + 2*(m-1) + j ; // normal order |
---|
948 | } |
---|
949 | } |
---|
950 | else if ( iside == 5 ) { |
---|
951 | // left side |
---|
952 | if ( i == 0 ) { return m*m - (j+1)*m ; } // reversed order |
---|
953 | else if ( i == n-1) { return 2*m*m - (j+1)*m ; } // reverded order |
---|
954 | else { |
---|
955 | if ( j == m-1 ) { return 2*m*m + 4*(i-1)*(m-1) ; } // special case |
---|
956 | else { return 2*m*m + 4*(i-1)*(m-1) + 3*(m-1) + j ; } // normal order |
---|
957 | } |
---|
958 | } |
---|
959 | |
---|
960 | else { |
---|
961 | |
---|
962 | G4cerr << "ERROR - G4VTwistSurface::GetNode(): " |
---|
963 | << GetName() << G4endl |
---|
964 | << " Not correct side number! - " << G4endl |
---|
965 | << "iside is " << iside << " but should be " |
---|
966 | << "0,1,2,3,4 or 5" << "." << G4endl ; |
---|
967 | G4Exception("G4TwistSurface::G4GetNode()", "InvalidSetup", |
---|
968 | FatalException, "Not correct side number."); |
---|
969 | } |
---|
970 | return -1 ; // wrong node |
---|
971 | } |
---|
972 | |
---|
973 | //===================================================================== |
---|
974 | //* GetEdgeVisiblility ------------------------------------------------ |
---|
975 | |
---|
976 | G4int G4VTwistSurface::GetEdgeVisibility( G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation) |
---|
977 | { |
---|
978 | |
---|
979 | // clockwise filling -> positive orientation |
---|
980 | // counter clockwise filling -> negative orientation |
---|
981 | |
---|
982 | // |
---|
983 | // d C c |
---|
984 | // +------+ |
---|
985 | // | | |
---|
986 | // | | |
---|
987 | // | | |
---|
988 | // D | |B |
---|
989 | // | | |
---|
990 | // | | |
---|
991 | // | | |
---|
992 | // +------+ |
---|
993 | // a A b |
---|
994 | // |
---|
995 | // a = +--+ A = ---+ |
---|
996 | // b = --++ B = --+- |
---|
997 | // c = -++- C = -+-- |
---|
998 | // d = ++-- D = +--- |
---|
999 | |
---|
1000 | |
---|
1001 | // check first invisible faces |
---|
1002 | |
---|
1003 | if ( ( i>0 && i<n-2 ) && ( j>0 && j<m-2 ) ) { |
---|
1004 | return -1 ; // always invisible, signs: ---- |
---|
1005 | } |
---|
1006 | |
---|
1007 | // change first the vertex number (depends on the orientation) |
---|
1008 | // 0,1,2,3 -> 3,2,1,0 |
---|
1009 | if ( orientation < 0 ) { number = ( 3 - number ) ; } |
---|
1010 | |
---|
1011 | // check true edges |
---|
1012 | if ( ( j>=1 && j<=m-3 ) ) { |
---|
1013 | |
---|
1014 | if ( i == 0 ) { // signs (A): ---+ |
---|
1015 | return ( number == 3 ) ? 1 : -1 ; |
---|
1016 | } |
---|
1017 | |
---|
1018 | else if ( i == n-2 ) { // signs (C): -+-- |
---|
1019 | return ( number == 1 ) ? 1 : -1 ; |
---|
1020 | } |
---|
1021 | |
---|
1022 | else { |
---|
1023 | G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): " |
---|
1024 | << GetName() << G4endl |
---|
1025 | << " Not correct face number! - " << G4endl ; |
---|
1026 | G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup", |
---|
1027 | FatalException, "Not correct face number."); |
---|
1028 | } |
---|
1029 | } |
---|
1030 | |
---|
1031 | if ( ( i>=1 && i<=n-3 ) ) { |
---|
1032 | |
---|
1033 | if ( j == 0 ) { // signs (D): +--- |
---|
1034 | return ( number == 0 ) ? 1 : -1 ; |
---|
1035 | } |
---|
1036 | |
---|
1037 | else if ( j == m-2 ) { // signs (B): --+- |
---|
1038 | return ( number == 2 ) ? 1 : -1 ; |
---|
1039 | } |
---|
1040 | |
---|
1041 | else { |
---|
1042 | G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): " |
---|
1043 | << GetName() << G4endl |
---|
1044 | << " Not correct face number! - " << G4endl ; |
---|
1045 | G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup", |
---|
1046 | FatalException, "Not correct face number."); |
---|
1047 | } |
---|
1048 | } |
---|
1049 | |
---|
1050 | // now the corners |
---|
1051 | if ( i == 0 && j == 0 ) { // signs (a) : +--+ |
---|
1052 | return ( number == 0 || number == 3 ) ? 1 : -1 ; |
---|
1053 | } |
---|
1054 | else if ( i == 0 && j == m-2 ) { // signs (b) : --++ |
---|
1055 | return ( number == 2 || number == 3 ) ? 1 : -1 ; |
---|
1056 | } |
---|
1057 | else if ( i == n-2 && j == m-2 ) { // signs (c) : -++- |
---|
1058 | return ( number == 1 || number == 2 ) ? 1 : -1 ; |
---|
1059 | } |
---|
1060 | else if ( i == n-2 && j == j ) { // signs (d) : ++-- |
---|
1061 | return ( number == 0 || number == 1 ) ? 1 : -1 ; |
---|
1062 | } |
---|
1063 | else { |
---|
1064 | G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): " |
---|
1065 | << GetName() << G4endl |
---|
1066 | << " Not correct face number! - " << G4endl ; |
---|
1067 | G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup", |
---|
1068 | FatalException, "Not correct face number."); |
---|
1069 | } |
---|
1070 | |
---|
1071 | G4cerr << "ERROR - G4VTwistSurface::GetEdgeVisibliity(): " |
---|
1072 | << GetName() << G4endl |
---|
1073 | << " Not correct face number! - " << G4endl ; |
---|
1074 | G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "InvalidSetup", |
---|
1075 | FatalException, "Not correct face number."); |
---|
1076 | |
---|
1077 | return 0 ; |
---|
1078 | |
---|
1079 | } |
---|
1080 | |
---|
1081 | |
---|
1082 | //===================================================================== |
---|
1083 | //* DebugPrint -------------------------------------------------------- |
---|
1084 | |
---|
1085 | void G4VTwistSurface::DebugPrint() const |
---|
1086 | { |
---|
1087 | G4ThreeVector A = fRot * GetCorner(sC0Min1Min) + fTrans; |
---|
1088 | G4ThreeVector B = fRot * GetCorner(sC0Max1Min) + fTrans; |
---|
1089 | G4ThreeVector C = fRot * GetCorner(sC0Max1Max) + fTrans; |
---|
1090 | G4ThreeVector D = fRot * GetCorner(sC0Min1Max) + fTrans; |
---|
1091 | |
---|
1092 | G4cout << "/* G4VTwistSurface::DebugPrint():-------------------------------" |
---|
1093 | << G4endl; |
---|
1094 | G4cout << "/* Name = " << fName << G4endl; |
---|
1095 | G4cout << "/* Axis = " << std::hex << fAxis[0] << " " |
---|
1096 | << std::hex << fAxis[1] |
---|
1097 | << " (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)" |
---|
1098 | << std::dec << G4endl; |
---|
1099 | G4cout << "/* BoundaryLimit(in local) fAxis0(min, max) = ("<<fAxisMin[0] |
---|
1100 | << ", " << fAxisMax[0] << ")" << G4endl; |
---|
1101 | G4cout << "/* BoundaryLimit(in local) fAxis1(min, max) = ("<<fAxisMin[1] |
---|
1102 | << ", " << fAxisMax[1] << ")" << G4endl; |
---|
1103 | G4cout << "/* Cornar point sC0Min1Min = " << A << G4endl; |
---|
1104 | G4cout << "/* Cornar point sC0Max1Min = " << B << G4endl; |
---|
1105 | G4cout << "/* Cornar point sC0Max1Max = " << C << G4endl; |
---|
1106 | G4cout << "/* Cornar point sC0Min1Max = " << D << G4endl; |
---|
1107 | G4cout << "/*---------------------------------------------------------" |
---|
1108 | << G4endl; |
---|
1109 | } |
---|
1110 | |
---|
1111 | //===================================================================== |
---|
1112 | // G4VTwistSurface::CurrentStatus class |
---|
1113 | //===================================================================== |
---|
1114 | |
---|
1115 | //===================================================================== |
---|
1116 | //* CurrentStatus::CurrentStatus -------------------------------------- |
---|
1117 | |
---|
1118 | G4VTwistSurface::CurrentStatus::CurrentStatus() |
---|
1119 | { |
---|
1120 | for (size_t i=0; i<G4VSURFACENXX; i++) |
---|
1121 | { |
---|
1122 | fDistance[i] = kInfinity; |
---|
1123 | fAreacode[i] = sOutside; |
---|
1124 | fIsValid[i] = false; |
---|
1125 | fXX[i].set(kInfinity, kInfinity, kInfinity); |
---|
1126 | } |
---|
1127 | fNXX = 0; |
---|
1128 | fLastp.set(kInfinity, kInfinity, kInfinity); |
---|
1129 | fLastv.set(kInfinity, kInfinity, kInfinity); |
---|
1130 | fLastValidate = kUninitialized; |
---|
1131 | fDone = false; |
---|
1132 | } |
---|
1133 | |
---|
1134 | //===================================================================== |
---|
1135 | //* CurrentStatus::~CurrentStatus ------------------------------------- |
---|
1136 | |
---|
1137 | G4VTwistSurface::CurrentStatus::~CurrentStatus() |
---|
1138 | { |
---|
1139 | } |
---|
1140 | |
---|
1141 | //===================================================================== |
---|
1142 | //* CurrentStatus::SetCurrentStatus ----------------------------------- |
---|
1143 | |
---|
1144 | void |
---|
1145 | G4VTwistSurface::CurrentStatus::SetCurrentStatus(G4int i, |
---|
1146 | G4ThreeVector &xx, |
---|
1147 | G4double &dist, |
---|
1148 | G4int &areacode, |
---|
1149 | G4bool &isvalid, |
---|
1150 | G4int nxx, |
---|
1151 | EValidate validate, |
---|
1152 | const G4ThreeVector *p, |
---|
1153 | const G4ThreeVector *v) |
---|
1154 | { |
---|
1155 | fDistance[i] = dist; |
---|
1156 | fAreacode[i] = areacode; |
---|
1157 | fIsValid[i] = isvalid; |
---|
1158 | fXX[i] = xx; |
---|
1159 | fNXX = nxx; |
---|
1160 | fLastValidate = validate; |
---|
1161 | if (p) |
---|
1162 | { |
---|
1163 | fLastp = *p; |
---|
1164 | } |
---|
1165 | else |
---|
1166 | { |
---|
1167 | G4Exception("G4VTwistSurface::CurrentStatus::CurrentStatus()", |
---|
1168 | "InvalidCondition", FatalException, |
---|
1169 | "SetCurrentStatus: p = 0!"); |
---|
1170 | } |
---|
1171 | if (v) |
---|
1172 | { |
---|
1173 | fLastv = *v; |
---|
1174 | } |
---|
1175 | else |
---|
1176 | { |
---|
1177 | fLastv.set(kInfinity, kInfinity, kInfinity); |
---|
1178 | } |
---|
1179 | fDone = true; |
---|
1180 | } |
---|
1181 | |
---|
1182 | //===================================================================== |
---|
1183 | //* CurrentStatus::ResetfDone ----------------------------------------- |
---|
1184 | |
---|
1185 | void |
---|
1186 | G4VTwistSurface::CurrentStatus::ResetfDone(EValidate validate, |
---|
1187 | const G4ThreeVector *p, |
---|
1188 | const G4ThreeVector *v) |
---|
1189 | |
---|
1190 | { |
---|
1191 | if (validate == fLastValidate && p && *p == fLastp) |
---|
1192 | { |
---|
1193 | if (!v || (v && *v == fLastv)) return; |
---|
1194 | } |
---|
1195 | G4ThreeVector xx(kInfinity, kInfinity, kInfinity); |
---|
1196 | for (size_t i=0; i<G4VSURFACENXX; i++) |
---|
1197 | { |
---|
1198 | fDistance[i] = kInfinity; |
---|
1199 | fAreacode[i] = sOutside; |
---|
1200 | fIsValid[i] = false; |
---|
1201 | fXX[i] = xx; // bug in old code ( was fXX[i] = xx[i] ) |
---|
1202 | } |
---|
1203 | fNXX = 0; |
---|
1204 | fLastp.set(kInfinity, kInfinity, kInfinity); |
---|
1205 | fLastv.set(kInfinity, kInfinity, kInfinity); |
---|
1206 | fLastValidate = kUninitialized; |
---|
1207 | fDone = false; |
---|
1208 | } |
---|
1209 | |
---|
1210 | //===================================================================== |
---|
1211 | //* CurrentStatus::DebugPrint ----------------------------------------- |
---|
1212 | |
---|
1213 | void |
---|
1214 | G4VTwistSurface::CurrentStatus::DebugPrint() const |
---|
1215 | { |
---|
1216 | G4cout << "CurrentStatus::Dist0,1= " << fDistance[0] |
---|
1217 | << " " << fDistance[1] << " areacode = " << fAreacode[0] |
---|
1218 | << " " << fAreacode[1] << G4endl; |
---|
1219 | } |
---|
1220 | |
---|
1221 | //===================================================================== |
---|
1222 | // G4VTwistSurface::Boundary class |
---|
1223 | //===================================================================== |
---|
1224 | |
---|
1225 | //===================================================================== |
---|
1226 | //* Boundary::Boundary ------------------------------------------------ |
---|
1227 | |
---|
1228 | G4VTwistSurface::Boundary::Boundary() |
---|
1229 | : fBoundaryAcode(-1), fBoundaryType(0) |
---|
1230 | { |
---|
1231 | } |
---|
1232 | |
---|
1233 | //===================================================================== |
---|
1234 | //* Boundary::~Boundary ----------------------------------------------- |
---|
1235 | |
---|
1236 | G4VTwistSurface::Boundary::~Boundary() |
---|
1237 | { |
---|
1238 | } |
---|
1239 | |
---|
1240 | //===================================================================== |
---|
1241 | //* Boundary::SetFields ----------------------------------------------- |
---|
1242 | |
---|
1243 | void |
---|
1244 | G4VTwistSurface::Boundary::SetFields(const G4int &areacode, |
---|
1245 | const G4ThreeVector &d, |
---|
1246 | const G4ThreeVector &x0, |
---|
1247 | const G4int &boundarytype) |
---|
1248 | { |
---|
1249 | fBoundaryAcode = areacode; |
---|
1250 | fBoundaryDirection = d; |
---|
1251 | fBoundaryX0 = x0; |
---|
1252 | fBoundaryType = boundarytype; |
---|
1253 | } |
---|
1254 | |
---|
1255 | //===================================================================== |
---|
1256 | //* Boundary::IsEmpty ------------------------------------------------- |
---|
1257 | |
---|
1258 | G4bool G4VTwistSurface::Boundary::IsEmpty() const |
---|
1259 | { |
---|
1260 | if (fBoundaryAcode == -1) return true; |
---|
1261 | return false; |
---|
1262 | } |
---|
1263 | |
---|
1264 | //===================================================================== |
---|
1265 | //* Boundary::GetBoundaryParameters ----------------------------------- |
---|
1266 | |
---|
1267 | G4bool |
---|
1268 | G4VTwistSurface::Boundary::GetBoundaryParameters(const G4int &areacode, |
---|
1269 | G4ThreeVector &d, |
---|
1270 | G4ThreeVector &x0, |
---|
1271 | G4int &boundarytype) const |
---|
1272 | { |
---|
1273 | // areacode must be one of them: |
---|
1274 | // sAxis0 & sAxisMin, sAxis0 & sAxisMax, |
---|
1275 | // sAxis1 & sAxisMin, sAxis1 & sAxisMax |
---|
1276 | if ((areacode & sAxis0) && (areacode & sAxis1)) |
---|
1277 | { |
---|
1278 | G4cerr << "ERROR - G4VTwistSurface::Boundary::GetBoundaryParameters()" |
---|
1279 | << G4endl |
---|
1280 | << " Located in the corner area. This function" |
---|
1281 | << " returns a direction vector of a boundary line." |
---|
1282 | << " areacode = " << areacode << G4endl; |
---|
1283 | G4Exception("G4VTwistSurface::Boundary::GetBoundaryParameters()", |
---|
1284 | "InvalidCondition", FatalException, |
---|
1285 | "Located in the corner area."); |
---|
1286 | } |
---|
1287 | if ((areacode & sSizeMask) != (fBoundaryAcode & sSizeMask)) |
---|
1288 | { |
---|
1289 | return false; |
---|
1290 | } |
---|
1291 | d = fBoundaryDirection; |
---|
1292 | x0 = fBoundaryX0; |
---|
1293 | boundarytype = fBoundaryType; |
---|
1294 | return true; |
---|
1295 | } |
---|