source: trunk/source/geometry/navigation/src/G4GeomTestVolume.cc@ 1349

Last change on this file since 1349 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 18.4 KB
Line 
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: G4GeomTestVolume.cc,v 1.6 2007/11/16 09:39:14 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// --------------------------------------------------------------------
31// GEANT 4 class source file
32//
33// G4GeomTestVolume
34//
35// Author: D.C.Williams, UCSC (davidw@scipp.ucsc.edu)
36// --------------------------------------------------------------------
37
38#include "G4GeomTestVolume.hh"
39
40#include "G4GeomTestLogger.hh"
41#include "G4GeomTestVolPoint.hh"
42#include "G4GeomTestSegment.hh"
43
44#include "G4VPhysicalVolume.hh"
45#include "G4LogicalVolume.hh"
46#include "G4VSolid.hh"
47
48#include <vector>
49#include <set>
50#include <algorithm>
51#include <iomanip>
52
53//
54// Constructor
55//
56G4GeomTestVolume::G4GeomTestVolume( const G4VPhysicalVolume *theTarget,
57 G4GeomTestLogger *theLogger,
58 G4double theTolerance )
59 : target(theTarget),
60 logger(theLogger),
61 tolerance(theTolerance),
62 extent(theTarget->GetLogicalVolume()->GetSolid()->GetExtent())
63{;}
64
65//
66// Destructor
67//
68G4GeomTestVolume::~G4GeomTestVolume() {;}
69
70//
71// Get error tolerance
72//
73G4double G4GeomTestVolume::GetTolerance() const
74{
75 return tolerance;
76}
77
78//
79// Set error tolerance
80//
81void G4GeomTestVolume::SetTolerance(G4double val)
82{
83 tolerance = val;
84}
85
86//
87// TestCartGridXYZ
88//
89void G4GeomTestVolume::TestCartGridXYZ( G4int nx, G4int ny, G4int nz )
90{
91 TestCartGridX( ny, nz );
92 TestCartGridY( nz, nx );
93 TestCartGridZ( nx, ny );
94}
95
96//
97// TestCartGridX
98//
99void G4GeomTestVolume::TestCartGridX( G4int ny, G4int nz )
100{
101 TestCartGrid( G4ThreeVector(0,1,0), G4ThreeVector(0,0,1),
102 G4ThreeVector(1,0,0), ny, nz );
103}
104
105//
106// TestCartGridY
107//
108void G4GeomTestVolume::TestCartGridY( G4int nz, G4int nx )
109{
110 TestCartGrid( G4ThreeVector(0,0,1), G4ThreeVector(1,0,0),
111 G4ThreeVector(0,1,0), nz, nx );
112}
113
114//
115// TestCartGridZ
116//
117void G4GeomTestVolume::TestCartGridZ( G4int nx, G4int ny )
118{
119 TestCartGrid( G4ThreeVector(1,0,0), G4ThreeVector(0,1,0),
120 G4ThreeVector(0,0,1), nx, ny );
121}
122
123//
124// TestRecursiveCartGrid
125//
126void G4GeomTestVolume::TestRecursiveCartGrid( G4int nx, G4int ny, G4int nz,
127 G4int slevel, G4int depth )
128{
129 // If reached requested level of depth (i.e. set to 0), exit.
130 // If not depth specified (i.e. set to -1), visit the whole tree.
131 // If requested initial level of depth is not zero, visit from beginning
132 //
133 if (depth == 0) return;
134 if (depth != -1) depth--;
135 if (slevel != 0) slevel--;
136
137 //
138 // As long as we aren't a replica and we reached the requested
139 // initial level of depth, test ourselves
140 //
141 if ( (!target->IsReplicated()) && (slevel==0) )
142 {
143 TestCartGridXYZ( nx, ny, nz );
144 ReportErrors();
145 }
146
147 //
148 // Loop over unique daughters
149 //
150 std::set<const G4LogicalVolume *> tested;
151
152 const G4LogicalVolume *logical = target->GetLogicalVolume();
153 G4int nDaughter = logical->GetNoDaughters();
154 G4int iDaughter;
155 for( iDaughter=0; iDaughter<nDaughter; ++iDaughter )
156 {
157 const G4VPhysicalVolume *daughter =
158 logical->GetDaughter(iDaughter);
159 const G4LogicalVolume *daughterLogical =
160 daughter->GetLogicalVolume();
161
162 //
163 // Skip empty daughters
164 //
165 if (daughterLogical->GetNoDaughters() == 0) continue;
166
167 //
168 // Tested already?
169 //
170 std::pair<std::set<const G4LogicalVolume *>::iterator,G4bool>
171 there = tested.insert(daughterLogical);
172 if (!there.second) continue;
173
174 //
175 // Recurse
176 //
177 G4GeomTestVolume vTest( daughter, logger, tolerance );
178 vTest.TestRecursiveCartGrid( nx,ny,nz,slevel,depth );
179 }
180}
181
182//
183// TestRecursiveCylinder
184//
185void
186G4GeomTestVolume::TestRecursiveCylinder( G4int nPhi, G4int nZ, G4int nRho,
187 G4double fracZ, G4double fracRho,
188 G4bool usePhi,
189 G4int slevel, G4int depth )
190{
191 // If reached requested level of depth (i.e. set to 0), exit.
192 // If not depth specified (i.e. set to -1), visit the whole tree.
193 // If requested initial level of depth is not zero, visit from beginning
194 //
195 if (depth == 0) return;
196 if (depth != -1) depth--;
197 if (slevel != 0) slevel--;
198
199 //
200 // As long as we aren't a replica and we reached the requested
201 // initial level of depth, test ourselves
202 //
203 if ( (!target->IsReplicated()) && (slevel==0) )
204 {
205 TestCylinder( nPhi, nZ, nRho, fracZ, fracRho, usePhi );
206 ReportErrors();
207 }
208
209 //
210 // Loop over unique daughters
211 //
212 std::set<const G4LogicalVolume *> tested;
213
214 const G4LogicalVolume *logical = target->GetLogicalVolume();
215 G4int nDaughter = logical->GetNoDaughters();
216 G4int iDaughter;
217 for( iDaughter=0; iDaughter<nDaughter; ++iDaughter )
218 {
219 const G4VPhysicalVolume *daughter =
220 logical->GetDaughter(iDaughter);
221 const G4LogicalVolume *daughterLogical =
222 daughter->GetLogicalVolume();
223
224 //
225 // Skip empty daughters
226 //
227 if (daughterLogical->GetNoDaughters() == 0) continue;
228
229 //
230 // Tested already?
231 //
232 std::pair<std::set<const G4LogicalVolume *>::iterator,G4bool>
233 there = tested.insert(daughterLogical);
234 if (!there.second) continue;
235
236 //
237 // Recurse
238 //
239 G4GeomTestVolume vTest( daughter, logger, tolerance );
240 vTest.TestRecursiveCylinder( nPhi, nZ, nRho,
241 fracZ, fracRho, usePhi,
242 slevel, depth );
243 }
244}
245
246//
247// TestCylinder
248//
249void G4GeomTestVolume::TestCylinder( G4int nPhi, G4int nZ, G4int nRho,
250 G4double fracZ, G4double fracRho,
251 G4bool usePhi )
252{
253 //
254 // Get size of our volume
255 //
256 G4double xMax = std::max(extent.GetXmax(),-extent.GetXmin());
257 G4double yMax = std::max(extent.GetYmax(),-extent.GetYmin());
258 G4double rhoMax = std::sqrt(xMax*xMax + yMax*yMax);
259
260 G4double zMax = extent.GetZmax();
261 G4double zMin = extent.GetZmin();
262 G4double zHalfLength = 0.5*(zMax-zMin);
263 G4double z0 = 0.5*(zMax+zMin);
264
265 //
266 // Loop over phi
267 //
268 G4double deltaPhi = 2*pi/G4double(nPhi);
269
270 G4double phi = 0;
271 G4int iPhi = nPhi;
272 if ((iPhi&1) == 0) iPhi++; // Also use odd number phi slices
273 do
274 {
275 G4double cosPhi = std::cos(phi);
276 G4double sinPhi = std::sin(phi);
277 G4ThreeVector vPhi1(sinPhi,-cosPhi,0);
278
279 //
280 // Loop over rho
281 //
282 G4double rho = rhoMax;
283 G4int iRho = nRho;
284 do
285 {
286 G4ThreeVector p(rho*cosPhi,rho*sinPhi,0);
287 static G4ThreeVector v(0,0,1);
288
289 TestOneLine( p, v );
290
291 if (usePhi)
292 {
293 //
294 // Loop over z
295 //
296 G4double zScale = 1.0;
297 G4int iZ=nZ;
298 do
299 {
300 p.setZ( z0 + zScale*zHalfLength );
301 TestOneLine(p,vPhi1);
302 p.setZ( z0 - zScale*zHalfLength );
303 TestOneLine(p,vPhi1);
304 } while( zScale *= fracZ, --iZ );
305 }
306 } while( rho *= fracRho, --iRho );
307
308 //
309 // Loop over z
310 //
311 G4ThreeVector p(0,0,0);
312 G4ThreeVector vPhi2(cosPhi,sinPhi,0);
313
314 G4double zScale = 1.0;
315 G4int iZ=nZ;
316 do
317 {
318 p.setZ( z0 + zScale*zHalfLength );
319
320 TestOneLine(p,vPhi2);
321
322 p.setZ( z0 - zScale*zHalfLength );
323
324 TestOneLine(p,vPhi2);
325 } while( zScale *= fracZ, --iZ );
326
327 } while( phi += deltaPhi, --iPhi );
328}
329
330//
331// TestCartGrid
332//
333void G4GeomTestVolume::TestCartGrid( const G4ThreeVector &theG1,
334 const G4ThreeVector &theG2,
335 const G4ThreeVector &theV,
336 G4int n1, G4int n2 )
337{
338 if (n1 <= 0 || n2 <= 0)
339 G4Exception( "G4GeomTestVolume::TestCartGrid()", "WrongArgumentValue",
340 FatalException, "Arguments n1 and n2 must be >= 1" );
341
342 G4ThreeVector xMin( extent.GetXmin(), extent.GetYmin(),
343 extent.GetZmin() );
344 G4ThreeVector xMax( extent.GetXmax(), extent.GetYmax(),
345 extent.GetZmax() );
346
347 G4ThreeVector g1(theG1.unit());
348 G4ThreeVector g2(theG2.unit());
349 G4ThreeVector v(theV.unit());
350
351 G4double gMin1 = g1.dot(xMin);
352 G4double gMax1 = g1.dot(xMax);
353
354 G4double gMin2 = g2.dot(xMin);
355 G4double gMax2 = g2.dot(xMax);
356
357 G4double delta1 = (gMax1-gMin1)/G4double(n1);
358 G4double delta2 = (gMax2-gMin2)/G4double(n2);
359
360 G4int i1, i2;
361
362 for(i1=0;i1<=n1;++i1)
363 {
364 G4ThreeVector p1 = (gMin1 + G4double(i1)*delta1)*g1;
365
366 for(i2=0;i2<=n2;++i2)
367 {
368 G4ThreeVector p2 = (gMin2 + G4double(i2)*delta2)*g2;
369
370 TestOneLine( p1+p2, v );
371 }
372 }
373}
374
375//
376// TestRecursiveLine
377//
378void
379G4GeomTestVolume::TestRecursiveLine( const G4ThreeVector& p,
380 const G4ThreeVector& v,
381 G4int slevel, G4int depth )
382{
383 // If reached requested level of depth (i.e. set to 0), exit.
384 // If not depth specified (i.e. set to -1), visit the whole tree.
385 // If requested initial level of depth is not zero, visit from beginning
386 //
387 if (depth == 0) return;
388 if (depth != -1) depth--;
389 if (slevel != 0) slevel--;
390
391 //
392 // As long as we aren't a replica and we reached the requested
393 // initial level of depth, test ourselves
394 //
395 if ( (!target->IsReplicated()) && (slevel==0) )
396 {
397 TestOneLine( p, v );
398 ReportErrors();
399 }
400
401 //
402 // Loop over unique daughters
403 //
404 std::set<const G4LogicalVolume *> tested;
405
406 const G4LogicalVolume *logical = target->GetLogicalVolume();
407 G4int nDaughter = logical->GetNoDaughters();
408 G4int iDaughter;
409 for( iDaughter=0; iDaughter<nDaughter; ++iDaughter )
410 {
411 const G4VPhysicalVolume *daughter =
412 logical->GetDaughter(iDaughter);
413 const G4LogicalVolume *daughterLogical =
414 daughter->GetLogicalVolume();
415
416 //
417 // Skip empty daughters
418 //
419 if (daughterLogical->GetNoDaughters() == 0) continue;
420
421 //
422 // Tested already?
423 //
424 std::pair<std::set<const G4LogicalVolume *>::iterator,G4bool>
425 there = tested.insert(daughterLogical);
426 if (!there.second) continue;
427
428 //
429 // Recurse
430 //
431 G4GeomTestVolume vTest( daughter, logger, tolerance );
432 vTest.TestRecursiveLine( p, v, slevel, depth );
433 }
434}
435
436//
437// TestOneLine
438//
439// Test geometry consistency along a single line
440//
441void G4GeomTestVolume::TestOneLine( const G4ThreeVector &p,
442 const G4ThreeVector &v )
443{
444 //
445 // Keep track of intersection points
446 //
447 std::vector<G4GeomTestVolPoint> points;
448
449 //
450 // Calculate intersections with the mother volume
451 //
452 G4GeomTestSegment targetSegment( target->GetLogicalVolume()->GetSolid(),
453 p, v, logger );
454
455 //
456 // Add them to our list
457 //
458 G4int n = targetSegment.GetNumberPoints();
459 G4int i;
460 for(i=0;i<n;++i)
461 {
462 points.push_back( G4GeomTestVolPoint( targetSegment.GetPoint(i), -1 ) );
463 }
464
465 //
466 // Loop over daughter volumes
467 //
468 const G4LogicalVolume *logical = target->GetLogicalVolume();
469 G4int nDaughter = logical->GetNoDaughters();
470 G4int iDaughter;
471 for( iDaughter=0; iDaughter<nDaughter; ++iDaughter)
472 {
473 const G4VPhysicalVolume *daughter =
474 logical->GetDaughter(iDaughter);
475
476 //
477 // Convert coordinates to daughter local coordinates
478 //
479 const G4RotationMatrix *rotation = daughter->GetFrameRotation();
480 const G4ThreeVector &translation = daughter->GetFrameTranslation();
481
482 G4ThreeVector pLocal = translation + p;
483 G4ThreeVector vLocal = v;
484
485 if (rotation)
486 {
487 pLocal = (*rotation)*pLocal;
488 vLocal = (*rotation)*vLocal;
489 }
490
491 //
492 // Find intersections
493 //
494 G4GeomTestSegment
495 daughterSegment( daughter->GetLogicalVolume()->GetSolid(),
496 pLocal, vLocal, logger );
497
498 //
499 // Add them to the list
500 //
501 n = daughterSegment.GetNumberPoints();
502 for(i=0;i<n;++i)
503 {
504 points.push_back( G4GeomTestVolPoint( daughterSegment.GetPoint(i),
505 iDaughter, translation, rotation ) );
506 }
507 }
508
509 //
510 // Now sort the list of intersection points
511 //
512 std::sort( points.begin(), points.end() );
513
514 //
515 // Search for problems:
516 // 1. Overlap daughters will be indicated by intersecting
517 // points that lie within another daughter
518 // 2. Daughter extending outside the mother will be
519 // indicated by intersecting points outside the mother
520 //
521 // The search method is always looking forward when
522 // resolving discrepencies due to roundoff error. Once
523 // one instance of a daughter intersection is analyzed,
524 // it is removed from further consideration
525 //
526 n = points.size();
527
528 //
529 // Set true if this point has been analyzed
530 //
531 std::vector<G4bool> checked( n, false );
532
533 for(i=0;i<n;++i)
534 {
535 if (checked[i]) continue;
536
537 G4int iDaug = points[i].GetDaughterIndex();
538 if (iDaug < 0) continue;
539
540 //
541 // Intersection found. Find matching pair.
542 //
543 G4double iS = points[i].GetDistance();
544 G4int j = i;
545 while(++j<n)
546 {
547 if (iDaug == points[j].GetDaughterIndex()) break;
548 }
549 if (j>=n)
550 {
551 //
552 // Unmatched? This shouldn't happen
553 //
554 logger->SolidProblem( logical->GetDaughter(iDaug)->
555 GetLogicalVolume()->GetSolid(),
556 "Unmatched intersection point",
557 points[i].GetPosition() );
558 continue;
559 }
560
561 checked[j] = true;
562
563 G4double jS = points[j].GetDistance();
564
565 //
566 // Otherwise, we could have a problem
567 //
568 G4int k = i;
569 while(++k<j)
570 {
571 if (checked[k]) continue;
572
573 G4bool kEntering = points[k].Entering();
574 G4double kS = points[k].GetDistance();
575
576
577 //
578 // Problem found: catagorize
579 //
580 G4int kDaug = points[k].GetDaughterIndex();
581 if (kDaug < 0)
582 {
583 //
584 // Ignore small overshoots if they are within tolerance
585 //
586 if (kS-iS < tolerance && kEntering ) continue;
587 if (jS-kS < tolerance && (!kEntering)) continue;
588
589 //
590 // We appear to extend outside the mother volume
591 //
592 std::map<G4long,G4GeomTestOvershootList>::iterator overshoot =
593 overshoots.find(iDaug);
594 if (overshoot == overshoots.end())
595 {
596 std::pair<std::map<G4long,G4GeomTestOvershootList>::iterator,G4bool>
597 result =
598 overshoots.insert( std::pair<const G4long,G4GeomTestOvershootList>
599 (iDaug,G4GeomTestOvershootList(target,iDaug)) );
600 assert(result.second);
601 overshoot = result.first;
602 }
603
604 if (kEntering)
605 (*overshoot).second.AddError( points[i].GetPosition(),
606 points[k].GetPosition() );
607 else
608 (*overshoot).second.AddError( points[k].GetPosition(),
609 points[j].GetPosition() );
610 }
611 else
612 {
613 //
614 // Ignore small overlaps if they are within tolerance
615 //
616 if (kS-iS < tolerance && (!kEntering)) continue;
617 if (jS-kS < tolerance && kEntering ) continue;
618
619 //
620 // We appear to overlap with another daughter
621 //
622 G4long key = iDaug < kDaug ?
623 (iDaug*nDaughter + kDaug) : (kDaug*nDaughter + iDaug);
624
625 std::map<G4long,G4GeomTestOverlapList>::iterator overlap =
626 overlaps.find(key);
627 if (overlap == overlaps.end())
628 {
629 std::pair<std::map<G4long,G4GeomTestOverlapList>::iterator,G4bool>
630 result =
631 overlaps.insert( std::pair<const G4long,G4GeomTestOverlapList>
632 (key,G4GeomTestOverlapList(target,iDaug,kDaug)) );
633 assert(result.second);
634 overlap = result.first;
635 }
636
637 if (points[k].Entering())
638 (*overlap).second.AddError( points[k].GetPosition(),
639 points[j].GetPosition() );
640 else
641 (*overlap).second.AddError( points[i].GetPosition(),
642 points[k].GetPosition() );
643 }
644 }
645 }
646}
647
648//
649// ReportErrors
650//
651void G4GeomTestVolume::ReportErrors()
652{
653 //
654 // Report overshoots
655 //
656 if (overshoots.empty())
657 logger->NoProblem("GeomTest: no daughter volume extending outside mother detected.");
658 else
659 {
660 std::map<G4long,G4GeomTestOvershootList>::iterator overshoot =
661 overshoots.begin();
662 while( overshoot != overshoots.end() )
663 {
664 logger->OvershootingDaughter( &(*overshoot).second );
665 ++overshoot;
666 }
667 }
668
669 //
670 // Report overlaps
671 //
672 if (overlaps.empty())
673 logger->NoProblem("GeomTest: no overlapping daughters detected.");
674 else
675 {
676 std::map<G4long,G4GeomTestOverlapList>::iterator overlap =
677 overlaps.begin();
678 while( overlap != overlaps.end() )
679 {
680 logger->OverlappingDaughters( &(*overlap).second );
681 ++overlap;
682 }
683 }
684}
685
686//
687// ClearErrors
688//
689void G4GeomTestVolume::ClearErrors()
690{
691 overshoots.clear();
692 overlaps.clear();
693}
Note: See TracBrowser for help on using the repository browser.